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