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