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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 162 } 163 164 /*@ 165 PetscPDFSampleGaussian1D - Sample uniformly from a Gaussian distribution in 1D 166 167 Not collective 168 169 Input Parameters: 170 + p - A uniform variable on [0, 1] 171 - dummy - ignored 172 173 Output Parameter: 174 . x - Coordinate in [-\infty, \infty] 175 176 Level: beginner 177 178 Note: http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function/ 179 https://stackoverflow.com/questions/27229371/inverse-error-function-in-c 180 @*/ 181 PetscErrorCode PetscPDFSampleGaussian1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 182 { 183 const PetscReal q = 2 * p[0] - 1.; 184 const PetscInt maxIter = 100; 185 PetscReal ck[100], r = 0.; 186 PetscInt k, m; 187 188 PetscFunctionBeginHot; 189 /* Transform input to [-1, 1] since the code below computes the inverse error function */ 190 for (k = 0; k < maxIter; ++k) ck[k] = 0.; 191 ck[0] = 1; 192 r = ck[0] * (PetscSqrtReal(PETSC_PI) / 2.) * q; 193 for (k = 1; k < maxIter; ++k) { 194 const PetscReal temp = 2. * k + 1.; 195 196 for (m = 0; m <= k - 1; ++m) { 197 PetscReal denom = (m + 1.) * (2. * m + 1.); 198 199 ck[k] += (ck[m] * ck[k - 1 - m]) / denom; 200 } 201 r += (ck[k] / temp) * PetscPowReal((PetscSqrtReal(PETSC_PI) / 2.) * q, 2. * k + 1.); 202 } 203 /* Scale erfinv() by \sqrt{\pi/2} */ 204 x[0] = PetscSqrtReal(PETSC_PI * 0.5) * r; 205 PetscFunctionReturn(PETSC_SUCCESS); 206 } 207 208 /*@ 209 PetscPDFGaussian2D - PDF for the Gaussian distribution in 2D 210 211 Not collective 212 213 Input Parameters: 214 + x - Coordinate in [-\infty, \infty]^2 215 - dummy - ignored 216 217 Output Parameter: 218 . p - The probability density at x 219 220 Level: beginner 221 222 .seealso: `PetscPDFSampleGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()` 223 @*/ 224 PetscErrorCode PetscPDFGaussian2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 225 { 226 p[0] = (1. / PETSC_PI) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1]))); 227 return PETSC_SUCCESS; 228 } 229 230 /*@ 231 PetscPDFSampleGaussian2D - Sample uniformly from a Gaussian distribution in 2D 232 233 Not collective 234 235 Input Parameters: 236 + p - A uniform variable on [0, 1]^2 237 - dummy - ignored 238 239 Output Parameter: 240 . x - Coordinate in [-\infty, \infty]^2 241 242 Level: beginner 243 244 Note: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform 245 246 .seealso: `PetscPDFGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()` 247 @*/ 248 PetscErrorCode PetscPDFSampleGaussian2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 249 { 250 const PetscReal mag = PetscSqrtReal(-2.0 * PetscLogReal(p[0])); 251 x[0] = mag * PetscCosReal(2.0 * PETSC_PI * p[1]); 252 x[1] = mag * PetscSinReal(2.0 * PETSC_PI * p[1]); 253 return PETSC_SUCCESS; 254 } 255 256 /*@ 257 PetscPDFGaussian3D - PDF for the Gaussian distribution in 3D 258 259 Not collective 260 261 Input Parameters: 262 + x - Coordinate in [-\infty, \infty]^3 263 - dummy - ignored 264 265 Output Parameter: 266 . p - The probability density at x 267 268 Level: beginner 269 270 .seealso: `PetscPDFSampleGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()` 271 @*/ 272 PetscErrorCode PetscPDFGaussian3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 273 { 274 p[0] = (1. / PETSC_PI * PetscSqrtReal(PETSC_PI)) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1]) + PetscSqr(x[2]))); 275 return PETSC_SUCCESS; 276 } 277 278 /*@ 279 PetscPDFSampleGaussian3D - Sample uniformly from a Gaussian distribution in 3D 280 281 Not collective 282 283 Input Parameters: 284 + p - A uniform variable on [0, 1]^3 285 - dummy - ignored 286 287 Output Parameter: 288 . x - Coordinate in [-\infty, \infty]^3 289 290 Level: beginner 291 292 Reference: 293 https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform 294 295 .seealso: `PetscPDFGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()` 296 @*/ 297 PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 298 { 299 PetscCall(PetscPDFSampleGaussian1D(p, dummy, x)); 300 PetscCall(PetscPDFSampleGaussian2D(&p[1], dummy, &x[1])); 301 return PETSC_SUCCESS; 302 } 303 304 /*@ 305 PetscPDFConstant1D - PDF for the uniform distribution in 1D 306 307 Not collective 308 309 Input Parameter: 310 . x - Coordinate in [-1, 1] 311 312 Output Parameter: 313 . p - The probability density at x 314 315 Level: beginner 316 317 .seealso: `PetscCDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscPDFConstant2D()`, `PetscPDFConstant3D()` 318 @*/ 319 PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 320 { 321 p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.; 322 return PETSC_SUCCESS; 323 } 324 325 /*@ 326 PetscCDFConstant1D - CDF for the uniform distribution in 1D 327 328 Not collective 329 330 Input Parameter: 331 . x - Coordinate in [-1, 1] 332 333 Output Parameter: 334 . p - The cumulative probability at x 335 336 Level: beginner 337 338 .seealso: `PetscPDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscCDFConstant2D()`, `PetscCDFConstant3D()` 339 @*/ 340 PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 341 { 342 p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.)); 343 return PETSC_SUCCESS; 344 } 345 346 /*@ 347 PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D 348 349 Not collective 350 351 Input Parameter: 352 . p - A uniform variable on [0, 1] 353 354 Output Parameter: 355 . x - Coordinate in [-1, 1] 356 357 Level: beginner 358 359 .seealso: `PetscPDFConstant1D()`, `PetscCDFConstant1D()`, `PetscPDFSampleConstant2D()`, `PetscPDFSampleConstant3D()` 360 @*/ 361 PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 362 { 363 x[0] = 2. * p[0] - 1.; 364 return PETSC_SUCCESS; 365 } 366 367 /*@ 368 PetscPDFConstant2D - PDF for the uniform distribution in 2D 369 370 Not collective 371 372 Input Parameter: 373 . x - Coordinate in [-1, 1] x [-1, 1] 374 375 Output Parameter: 376 . p - The probability density at x 377 378 Level: beginner 379 380 .seealso: `PetscCDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscPDFConstant1D()`, `PetscPDFConstant3D()` 381 @*/ 382 PetscErrorCode PetscPDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 383 { 384 p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. ? 0.25 : 0.; 385 return PETSC_SUCCESS; 386 } 387 388 /*@ 389 PetscCDFConstant2D - CDF for the uniform distribution in 2D 390 391 Not collective 392 393 Input Parameter: 394 . x - Coordinate in [-1, 1] x [-1, 1] 395 396 Output Parameter: 397 . p - The cumulative probability at x 398 399 Level: beginner 400 401 .seealso: `PetscPDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscCDFConstant1D()`, `PetscCDFConstant3D()` 402 @*/ 403 PetscErrorCode PetscCDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 404 { 405 p[0] = x[0] <= -1. || x[1] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.)) * (x[1] > 1. ? 1. : 0.5 * (x[1] + 1.)); 406 return PETSC_SUCCESS; 407 } 408 409 /*@ 410 PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 2D 411 412 Not collective 413 414 Input Parameter: 415 . p - Two uniform variables on [0, 1] 416 417 Output Parameter: 418 . x - Coordinate in [-1, 1] x [-1, 1] 419 420 Level: beginner 421 422 .seealso: `PetscPDFConstant2D()`, `PetscCDFConstant2D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant3D()` 423 @*/ 424 PetscErrorCode PetscPDFSampleConstant2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 425 { 426 x[0] = 2. * p[0] - 1.; 427 x[1] = 2. * p[1] - 1.; 428 return PETSC_SUCCESS; 429 } 430 431 /*@ 432 PetscPDFConstant3D - PDF for the uniform distribution in 3D 433 434 Not collective 435 436 Input Parameter: 437 . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1] 438 439 Output Parameter: 440 . p - The probability density at x 441 442 Level: beginner 443 444 .seealso: `PetscCDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()` 445 @*/ 446 PetscErrorCode PetscPDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 447 { 448 p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. && x[2] > -1. && x[2] <= 1. ? 0.125 : 0.; 449 return PETSC_SUCCESS; 450 } 451 452 /*@ 453 PetscCDFConstant3D - CDF for the uniform distribution in 3D 454 455 Not collective 456 457 Input Parameter: 458 . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1] 459 460 Output Parameter: 461 . p - The cumulative probability at x 462 463 Level: beginner 464 465 .seealso: `PetscPDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscCDFConstant1D()`, `PetscCDFConstant2D()` 466 @*/ 467 PetscErrorCode PetscCDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 468 { 469 p[0] = x[0] <= -1. || x[1] <= -1. || x[2] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.)) * (x[1] > 1. ? 1. : 0.5 * (x[1] + 1.)) * (x[2] > 1. ? 1. : 0.5 * (x[2] + 1.)); 470 return PETSC_SUCCESS; 471 } 472 473 /*@ 474 PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 3D 475 476 Not collective 477 478 Input Parameter: 479 . p - Three uniform variables on [0, 1] 480 481 Output Parameter: 482 . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1] 483 484 Level: beginner 485 486 .seealso: `PetscPDFConstant3D()`, `PetscCDFConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()` 487 @*/ 488 PetscErrorCode PetscPDFSampleConstant3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 489 { 490 x[0] = 2. * p[0] - 1.; 491 x[1] = 2. * p[1] - 1.; 492 x[2] = 2. * p[2] - 1.; 493 return PETSC_SUCCESS; 494 } 495 496 /*@C 497 PetscProbCreateFromOptions - Return the probability distribution specified by the arguments and options 498 499 Not collective 500 501 Input Parameters: 502 + dim - The dimension of sample points 503 . prefix - The options prefix, or NULL 504 - name - The option name for the probility distribution type 505 506 Output Parameters: 507 + pdf - The PDF of this type 508 . cdf - The CDF of this type 509 - sampler - The PDF sampler of this type 510 511 Level: intermediate 512 513 .seealso: `PetscProbFunc`, `PetscPDFMaxwellBoltzmann1D()`, `PetscPDFGaussian1D()`, `PetscPDFConstant1D()` 514 @*/ 515 PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler) 516 { 517 DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN; 518 519 PetscFunctionBegin; 520 PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT"); 521 PetscCall(PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum)den, (PetscEnum *)&den, NULL)); 522 PetscOptionsEnd(); 523 524 if (pdf) { 525 PetscValidPointer(pdf, 4); 526 *pdf = NULL; 527 } 528 if (cdf) { 529 PetscValidPointer(cdf, 5); 530 *cdf = NULL; 531 } 532 if (sampler) { 533 PetscValidPointer(sampler, 6); 534 *sampler = NULL; 535 } 536 switch (den) { 537 case DTPROB_DENSITY_CONSTANT: 538 switch (dim) { 539 case 1: 540 if (pdf) *pdf = PetscPDFConstant1D; 541 if (cdf) *cdf = PetscCDFConstant1D; 542 if (sampler) *sampler = PetscPDFSampleConstant1D; 543 break; 544 case 2: 545 if (pdf) *pdf = PetscPDFConstant2D; 546 if (cdf) *cdf = PetscCDFConstant2D; 547 if (sampler) *sampler = PetscPDFSampleConstant2D; 548 break; 549 case 3: 550 if (pdf) *pdf = PetscPDFConstant3D; 551 if (cdf) *cdf = PetscCDFConstant3D; 552 if (sampler) *sampler = PetscPDFSampleConstant3D; 553 break; 554 default: 555 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 556 } 557 break; 558 case DTPROB_DENSITY_GAUSSIAN: 559 switch (dim) { 560 case 1: 561 if (pdf) *pdf = PetscPDFGaussian1D; 562 if (cdf) *cdf = PetscCDFGaussian1D; 563 if (sampler) *sampler = PetscPDFSampleGaussian1D; 564 break; 565 case 2: 566 if (pdf) *pdf = PetscPDFGaussian2D; 567 if (sampler) *sampler = PetscPDFSampleGaussian2D; 568 break; 569 case 3: 570 if (pdf) *pdf = PetscPDFGaussian3D; 571 if (sampler) *sampler = PetscPDFSampleGaussian3D; 572 break; 573 default: 574 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 575 } 576 break; 577 case DTPROB_DENSITY_MAXWELL_BOLTZMANN: 578 switch (dim) { 579 case 1: 580 if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D; 581 if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D; 582 break; 583 case 2: 584 if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D; 585 if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D; 586 break; 587 case 3: 588 if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D; 589 if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D; 590 break; 591 default: 592 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 593 } 594 break; 595 default: 596 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]); 597 } 598 PetscFunctionReturn(PETSC_SUCCESS); 599 } 600 601 #ifdef PETSC_HAVE_KS 602 EXTERN_C_BEGIN 603 #include <KolmogorovSmirnovDist.h> 604 EXTERN_C_END 605 #endif 606 607 /*@C 608 PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF. 609 610 Collective on v 611 612 Input Parameters: 613 + v - The data vector, blocksize is the sample dimension 614 - cdf - The analytic CDF 615 616 Output Parameter: 617 . alpha - The KS statisic 618 619 Level: advanced 620 621 Notes: 622 The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is 623 .vb 624 D_n = \sup_x \left| F_n(x) - F(x) \right| 625 626 where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by 627 628 F_n = # of samples <= x / n 629 .ve 630 631 The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep'' 632 cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes 633 the largest absolute difference between the two distribution functions across all $x$ values. 634 635 .seealso: `PetscProbFunc` 636 @*/ 637 PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha) 638 { 639 #if !defined(PETSC_HAVE_KS) 640 SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks"); 641 #else 642 PetscViewer viewer = NULL; 643 PetscViewerFormat format; 644 PetscDraw draw; 645 PetscDrawLG lg; 646 PetscDrawAxis axis; 647 const PetscScalar *a; 648 PetscReal *speed, Dn = PETSC_MIN_REAL; 649 PetscBool isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg; 650 PetscInt n, p, dim, d; 651 PetscMPIInt size; 652 const char *names[2] = {"Analytic", "Empirical"}; 653 char title[PETSC_MAX_PATH_LEN]; 654 PetscOptions options; 655 const char *prefix; 656 MPI_Comm comm; 657 658 PetscFunctionBegin; 659 PetscCall(PetscObjectGetComm((PetscObject)v, &comm)); 660 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, &prefix)); 661 PetscCall(PetscObjectGetOptions((PetscObject)v, &options)); 662 PetscCall(PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg)); 663 if (flg) { 664 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 665 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 666 } 667 if (isascii) { 668 PetscCall(PetscViewerPushFormat(viewer, format)); 669 } else if (isdraw) { 670 PetscCall(PetscViewerPushFormat(viewer, format)); 671 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 672 PetscCall(PetscDrawLGCreate(draw, 2, &lg)); 673 PetscCall(PetscDrawLGSetLegend(lg, names)); 674 } 675 676 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size)); 677 PetscCall(VecGetLocalSize(v, &n)); 678 PetscCall(VecGetBlockSize(v, &dim)); 679 n /= dim; 680 /* TODO Parallel algorithm is harder */ 681 if (size == 1) { 682 PetscCall(PetscMalloc1(n, &speed)); 683 PetscCall(VecGetArrayRead(v, &a)); 684 for (p = 0; p < n; ++p) { 685 PetscReal mag = 0.; 686 687 for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d])); 688 speed[p] = PetscSqrtReal(mag); 689 } 690 PetscCall(PetscSortReal(n, speed)); 691 PetscCall(VecRestoreArrayRead(v, &a)); 692 for (p = 0; p < n; ++p) { 693 const PetscReal x = speed[p], Fn = ((PetscReal)p) / n; 694 PetscReal F, vals[2]; 695 696 PetscCall(cdf(&x, NULL, &F)); 697 Dn = PetscMax(PetscAbsReal(Fn - F), Dn); 698 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn)); 699 if (isdraw) { 700 vals[0] = F; 701 vals[1] = Fn; 702 PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals)); 703 } 704 } 705 if (speed[n - 1] < 6.) { 706 const PetscReal k = (PetscInt)(6. - speed[n - 1]) + 1, dx = (6. - speed[n - 1]) / k; 707 for (p = 0; p <= k; ++p) { 708 const PetscReal x = speed[n - 1] + p * dx, Fn = 1.0; 709 PetscReal F, vals[2]; 710 711 PetscCall(cdf(&x, NULL, &F)); 712 Dn = PetscMax(PetscAbsReal(Fn - F), Dn); 713 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn)); 714 if (isdraw) { 715 vals[0] = F; 716 vals[1] = Fn; 717 PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals)); 718 } 719 } 720 } 721 PetscCall(PetscFree(speed)); 722 } 723 if (isdraw) { 724 PetscCall(PetscDrawLGGetAxis(lg, &axis)); 725 PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn)); 726 PetscCall(PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)")); 727 PetscCall(PetscDrawLGDraw(lg)); 728 PetscCall(PetscDrawLGDestroy(&lg)); 729 } 730 if (viewer) { 731 PetscCall(PetscViewerFlush(viewer)); 732 PetscCall(PetscViewerPopFormat(viewer)); 733 PetscCall(PetscViewerDestroy(&viewer)); 734 } 735 *alpha = KSfbar((int)n, (double)Dn); 736 PetscFunctionReturn(PETSC_SUCCESS); 737 #endif 738 } 739