1 #include <petscdt.h> /*I "petscdt.h" I*/ 2 3 #include <petscvec.h> 4 #include <petscdraw.h> 5 #include <petsc/private/petscimpl.h> 6 7 const char * const DTProbDensityTypes[] = {"constant", "gaussian", "maxwell_boltzmann", "DTProb Density", "DTPROB_DENSITY_", NULL}; 8 9 /*@ 10 PetscPDFMaxwellBoltzmann1D - PDF for the Maxwell-Boltzmann distribution in 1D 11 12 Not collective 13 14 Input Parameter: 15 . x - Speed in [0, \infty] 16 17 Output Parameter: 18 . p - The probability density at x 19 20 Level: beginner 21 22 .seealso: PetscCDFMaxwellBoltzmann1D 23 @*/ 24 PetscErrorCode PetscPDFMaxwellBoltzmann1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 25 { 26 p[0] = PetscSqrtReal(2. / PETSC_PI) * PetscExpReal(-0.5 * PetscSqr(x[0])); 27 return 0; 28 } 29 30 /*@ 31 PetscCDFMaxwellBoltzmann1D - CDF for the Maxwell-Boltzmann distribution in 1D 32 33 Not collective 34 35 Input Parameter: 36 . x - Speed in [0, \infty] 37 38 Output Parameter: 39 . p - The probability density at x 40 41 Level: beginner 42 43 .seealso: PetscPDFMaxwellBoltzmann1D 44 @*/ 45 PetscErrorCode PetscCDFMaxwellBoltzmann1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 46 { 47 p[0] = PetscErfReal(x[0] / PETSC_SQRT2); 48 return 0; 49 } 50 51 /*@ 52 PetscPDFMaxwellBoltzmann2D - PDF for the Maxwell-Boltzmann distribution in 2D 53 54 Not collective 55 56 Input Parameter: 57 . x - Speed in [0, \infty] 58 59 Output Parameter: 60 . p - The probability density at x 61 62 Level: beginner 63 64 .seealso: PetscCDFMaxwellBoltzmann2D 65 @*/ 66 PetscErrorCode PetscPDFMaxwellBoltzmann2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 67 { 68 p[0] = x[0] * PetscExpReal(-0.5 * PetscSqr(x[0])); 69 return 0; 70 } 71 72 /*@ 73 PetscCDFMaxwellBoltzmann2D - CDF for the Maxwell-Boltzmann distribution in 2D 74 75 Not collective 76 77 Input Parameter: 78 . x - Speed in [0, \infty] 79 80 Output Parameter: 81 . p - The probability density at x 82 83 Level: beginner 84 85 .seealso: PetscPDFMaxwellBoltzmann2D 86 @*/ 87 PetscErrorCode PetscCDFMaxwellBoltzmann2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 88 { 89 p[0] = 1. - PetscExpReal(-0.5 * PetscSqr(x[0])); 90 return 0; 91 } 92 93 /*@ 94 PetscPDFMaxwellBoltzmann3D - PDF for the Maxwell-Boltzmann distribution in 3D 95 96 Not collective 97 98 Input Parameter: 99 . x - Speed in [0, \infty] 100 101 Output Parameter: 102 . p - The probability density at x 103 104 Level: beginner 105 106 .seealso: PetscCDFMaxwellBoltzmann3D 107 @*/ 108 PetscErrorCode PetscPDFMaxwellBoltzmann3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 109 { 110 p[0] = PetscSqrtReal(2. / PETSC_PI) * PetscSqr(x[0]) * PetscExpReal(-0.5 * PetscSqr(x[0])); 111 return 0; 112 } 113 114 /*@ 115 PetscCDFMaxwellBoltzmann3D - CDF for the Maxwell-Boltzmann distribution in 3D 116 117 Not collective 118 119 Input Parameter: 120 . x - Speed in [0, \infty] 121 122 Output Parameter: 123 . p - The probability density at x 124 125 Level: beginner 126 127 .seealso: PetscPDFMaxwellBoltzmann3D 128 @*/ 129 PetscErrorCode PetscCDFMaxwellBoltzmann3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 130 { 131 p[0] = PetscErfReal(x[0] / PETSC_SQRT2) - PetscSqrtReal(2. / PETSC_PI) * x[0] * PetscExpReal(-0.5 * PetscSqr(x[0])); 132 return 0; 133 } 134 135 /*@ 136 PetscPDFGaussian1D - PDF for the Gaussian distribution in 1D 137 138 Not collective 139 140 Input Parameter: 141 . x - Coordinate in [-\infty, \infty] 142 143 Output Parameter: 144 . p - The probability density at x 145 146 Level: beginner 147 148 .seealso: PetscPDFMaxwellBoltzmann3D 149 @*/ 150 PetscErrorCode PetscPDFGaussian1D(const PetscReal x[], const PetscReal scale[], PetscReal p[]) 151 { 152 const PetscReal sigma = scale ? scale[0] : 1.; 153 p[0] = PetscSqrtReal(1. / (2.*PETSC_PI)) * PetscExpReal(-0.5 * PetscSqr(x[0] / sigma)) / sigma; 154 return 0; 155 } 156 157 PetscErrorCode PetscCDFGaussian1D(const PetscReal x[], const PetscReal scale[], PetscReal p[]) 158 { 159 const PetscReal sigma = scale ? scale[0] : 1.; 160 p[0] = 0.5 * (1. + PetscErfReal(x[0] / PETSC_SQRT2 / sigma)); 161 return 0; 162 } 163 164 /*@ 165 PetscPDFSampleGaussian1D - Sample uniformly from a Gaussian distribution in 1D 166 167 Not collective 168 169 Input Parameter: 170 . p - A uniform variable on [0, 1] 171 172 Output Parameter: 173 . x - Coordinate in [-\infty, \infty] 174 175 Level: beginner 176 177 Note: http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/ 178 https://stackoverflow.com/questions/27229371/inverse-error-function-in-c 179 @*/ 180 PetscErrorCode PetscPDFSampleGaussian1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 181 { 182 const PetscReal q = 2*p[0] - 1.; 183 const PetscInt maxIter = 100; 184 PetscReal ck[100], r = 0.; 185 PetscInt k, m; 186 187 PetscFunctionBeginHot; 188 /* Transform input to [-1, 1] since the code below computes the inverse error function */ 189 for (k = 0; k < maxIter; ++k) ck[k] = 0.; 190 ck[0] = 1; 191 r = ck[0]* (PetscSqrtReal(PETSC_PI) / 2.) * q; 192 for (k = 1; k < maxIter; ++k) { 193 const PetscReal temp = 2.*k + 1.; 194 195 for (m = 0; m <= k-1; ++m) { 196 PetscReal denom = (m + 1.)*(2.*m + 1.); 197 198 ck[k] += (ck[m]*ck[k-1-m])/denom; 199 } 200 r += (ck[k]/temp) * PetscPowReal((PetscSqrtReal(PETSC_PI)/2.) * q, 2.*k+1.); 201 } 202 /* Scale erfinv() by \sqrt{\pi/2} */ 203 x[0] = PetscSqrtReal(PETSC_PI * 0.5) * r; 204 PetscFunctionReturn(0); 205 } 206 207 /*@ 208 PetscPDFGaussian2D - PDF for the Gaussian distribution in 2D 209 210 Not collective 211 212 Input Parameter: 213 . x - Coordinate in [-\infty, \infty]^2 214 215 Output Parameter: 216 . p - The probability density at x 217 218 Level: beginner 219 220 .seealso: PetscPDFSampleGaussian2D(), PetscPDFMaxwellBoltzmann3D() 221 @*/ 222 PetscErrorCode PetscPDFGaussian2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 223 { 224 p[0] = (1. / PETSC_PI) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1]))); 225 return 0; 226 } 227 228 /*@ 229 PetscPDFSampleGaussian2D - Sample uniformly from a Gaussian distribution in 2D 230 231 Not collective 232 233 Input Parameter: 234 . p - A uniform variable on [0, 1]^2 235 236 Output Parameter: 237 . x - Coordinate in [-\infty, \infty]^2 238 239 Level: beginner 240 241 Note: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform 242 243 .seealso: PetscPDFGaussian2D(), PetscPDFMaxwellBoltzmann3D() 244 @*/ 245 PetscErrorCode PetscPDFSampleGaussian2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 246 { 247 const PetscReal mag = PetscSqrtReal(-2.0 * PetscLogReal(p[0])); 248 x[0] = mag * PetscCosReal(2.0 * PETSC_PI * p[1]); 249 x[1] = mag * PetscSinReal(2.0 * PETSC_PI * p[1]); 250 return 0; 251 } 252 253 /*@ 254 PetscPDFConstant1D - PDF for the uniform distribution in 1D 255 256 Not collective 257 258 Input Parameter: 259 . x - Coordinate in [-1, 1] 260 261 Output Parameter: 262 . p - The probability density at x 263 264 Level: beginner 265 266 .seealso: PetscCDFConstant1D() 267 @*/ 268 PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 269 { 270 p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.; 271 return 0; 272 } 273 274 PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 275 { 276 p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5*(x[0] + 1.)); 277 return 0; 278 } 279 280 /*@ 281 PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D 282 283 Not collective 284 285 Input Parameter: 286 . p - A uniform variable on [0, 1] 287 288 Output Parameter: 289 . x - Coordinate in [-1, 1] 290 291 Level: beginner 292 293 .seealso: PetscPDFConstant1D(), PetscCDFConstant1D() 294 @*/ 295 PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 296 { 297 x[0] = 2.*p[1] - 1.; 298 return 0; 299 } 300 301 /*@C 302 PetscProbCreateFromOptions - Return the probability distribution specified by the argumetns and options 303 304 Not collective 305 306 Input Parameters: 307 + dim - The dimension of sample points 308 . prefix - The options prefix, or NULL 309 - name - The option name for the probility distribution type 310 311 Output Parameters: 312 + pdf - The PDF of this type 313 . cdf - The CDF of this type 314 - sampler - The PDF sampler of this type 315 316 Level: intermediate 317 318 .seealso: PetscPDFMaxwellBoltzmann1D(), PetscPDFGaussian1D(), PetscPDFConstant1D() 319 @*/ 320 PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler) 321 { 322 DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 ierr = PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");CHKERRQ(ierr); 327 ierr = PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum) den, (PetscEnum *) &den, NULL);CHKERRQ(ierr); 328 ierr = PetscOptionsEnd();CHKERRQ(ierr); 329 330 if (pdf) {PetscValidPointer(pdf, 4); *pdf = NULL;} 331 if (cdf) {PetscValidPointer(cdf, 5); *cdf = NULL;} 332 if (sampler) {PetscValidPointer(sampler, 6); *sampler = NULL;} 333 switch (den) { 334 case DTPROB_DENSITY_CONSTANT: 335 switch (dim) { 336 case 1: 337 if (pdf) *pdf = PetscPDFConstant1D; 338 if (cdf) *cdf = PetscCDFConstant1D; 339 if (sampler) *sampler = PetscPDFSampleConstant1D; 340 break; 341 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 342 } 343 break; 344 case DTPROB_DENSITY_GAUSSIAN: 345 switch (dim) { 346 case 1: 347 if (pdf) *pdf = PetscPDFGaussian1D; 348 if (cdf) *cdf = PetscCDFGaussian1D; 349 if (sampler) *sampler = PetscPDFSampleGaussian1D; 350 break; 351 case 2: 352 if (pdf) *pdf = PetscPDFGaussian2D; 353 if (sampler) *sampler = PetscPDFSampleGaussian2D; 354 break; 355 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 356 } 357 break; 358 case DTPROB_DENSITY_MAXWELL_BOLTZMANN: 359 switch (dim) { 360 case 1: 361 if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D; 362 if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D; 363 break; 364 case 2: 365 if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D; 366 if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D; 367 break; 368 case 3: 369 if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D; 370 if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D; 371 break; 372 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 373 } 374 break; 375 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]); 376 } 377 PetscFunctionReturn(0); 378 } 379 380 #ifdef PETSC_HAVE_KS 381 EXTERN_C_BEGIN 382 #include <KolmogorovSmirnovDist.h> 383 EXTERN_C_END 384 #endif 385 386 /*@C 387 PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF. 388 389 Collective on v 390 391 Input Parameters: 392 + v - The data vector, blocksize is the sample dimension 393 - cdf - The analytic CDF 394 395 Output Parameter: 396 . alpha - The KS statisic 397 398 Level: advanced 399 400 Note: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is 401 $ 402 $ D_n = \sup_x \left| F_n(x) - F(x) \right| 403 $ 404 where $\sup_x$ is the supremum of the set of distances, and the empirical distribution 405 function $F_n(x)$ is discrete, and given by 406 $ 407 $ F_n = # of samples <= x / n 408 $ 409 The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep'' 410 cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes 411 the largest absolute difference between the two distribution functions across all $x$ values. 412 413 .seealso: PetscProbFunc 414 @*/ 415 PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha) 416 { 417 #if !defined(PETSC_HAVE_KS) 418 SETERRQ(PetscObjectComm((PetscObject) v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks"); 419 #else 420 PetscViewer viewer = NULL; 421 PetscViewerFormat format; 422 PetscDraw draw; 423 PetscDrawLG lg; 424 PetscDrawAxis axis; 425 const PetscScalar *a; 426 PetscReal *speed, Dn = PETSC_MIN_REAL; 427 PetscBool isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg; 428 PetscInt n, p, dim, d; 429 PetscMPIInt size; 430 const char *names[2] = {"Analytic", "Empirical"}; 431 char title[PETSC_MAX_PATH_LEN]; 432 PetscOptions options; 433 const char *prefix; 434 MPI_Comm comm; 435 PetscErrorCode ierr; 436 437 PetscFunctionBegin; 438 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr); 439 ierr = PetscObjectGetOptionsPrefix((PetscObject) v, &prefix);CHKERRQ(ierr); 440 ierr = PetscObjectGetOptions((PetscObject) v, &options);CHKERRQ(ierr); 441 ierr = PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg);CHKERRQ(ierr); 442 if (flg) { 443 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 444 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);CHKERRQ(ierr); 445 } 446 if (isascii) { 447 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 448 } else if (isdraw) { 449 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 450 ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);CHKERRQ(ierr); 451 ierr = PetscDrawLGCreate(draw, 2, &lg);CHKERRQ(ierr); 452 ierr = PetscDrawLGSetLegend(lg, names);CHKERRQ(ierr); 453 } 454 455 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) v), &size);CHKERRMPI(ierr); 456 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 457 ierr = VecGetBlockSize(v, &dim);CHKERRQ(ierr); 458 n /= dim; 459 /* TODO Parallel algorithm is harder */ 460 if (size == 1) { 461 ierr = PetscMalloc1(n, &speed);CHKERRQ(ierr); 462 ierr = VecGetArrayRead(v, &a);CHKERRQ(ierr); 463 for (p = 0; p < n; ++p) { 464 PetscReal mag = 0.; 465 466 for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p*dim+d])); 467 speed[p] = PetscSqrtReal(mag); 468 } 469 ierr = PetscSortReal(n, speed);CHKERRQ(ierr); 470 ierr = VecRestoreArrayRead(v, &a);CHKERRQ(ierr); 471 for (p = 0; p < n; ++p) { 472 const PetscReal x = speed[p], Fn = ((PetscReal) p) / n; 473 PetscReal F, vals[2]; 474 475 ierr = cdf(&x, NULL, &F);CHKERRQ(ierr); 476 Dn = PetscMax(PetscAbsReal(Fn - F), Dn); 477 if (isascii) {ierr = PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);CHKERRQ(ierr);} 478 if (isdraw) {vals[0] = F; vals[1] = Fn; ierr = PetscDrawLGAddCommonPoint(lg, x, vals);CHKERRQ(ierr);} 479 } 480 if (speed[n-1] < 6.) { 481 const PetscReal k = (PetscInt) (6. - speed[n-1]) + 1, dx = (6. - speed[n-1])/k; 482 for (p = 0; p <= k; ++p) { 483 const PetscReal x = speed[n-1] + p*dx, Fn = 1.0; 484 PetscReal F, vals[2]; 485 486 ierr = cdf(&x, NULL, &F);CHKERRQ(ierr); 487 Dn = PetscMax(PetscAbsReal(Fn - F), Dn); 488 if (isascii) {ierr = PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);CHKERRQ(ierr);} 489 if (isdraw) {vals[0] = F; vals[1] = Fn; ierr = PetscDrawLGAddCommonPoint(lg, x, vals);CHKERRQ(ierr);} 490 } 491 } 492 ierr = PetscFree(speed);CHKERRQ(ierr); 493 } 494 if (isdraw) { 495 ierr = PetscDrawLGGetAxis(lg, &axis);CHKERRQ(ierr); 496 ierr = PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn);CHKERRQ(ierr); 497 ierr = PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)");CHKERRQ(ierr); 498 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 499 ierr = PetscDrawLGDestroy(&lg);CHKERRQ(ierr); 500 } 501 if (viewer) { 502 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 503 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 504 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 505 } 506 *alpha = KSfbar((int) n, (double) Dn); 507 PetscFunctionReturn(0); 508 #endif 509 } 510