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 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(0); 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 0; 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 0; 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 0; 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 Note: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform 293 294 .seealso: PetscPDFGaussian3D(), PetscPDFMaxwellBoltzmann3D() 295 @*/ 296 PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 297 { 298 PetscCall(PetscPDFSampleGaussian1D(p, dummy, x)); 299 PetscCall(PetscPDFSampleGaussian2D(&p[1], dummy, &x[1])); 300 return 0; 301 } 302 303 /*@ 304 PetscPDFConstant1D - PDF for the uniform distribution in 1D 305 306 Not collective 307 308 Input Parameter: 309 . x - Coordinate in [-1, 1] 310 311 Output Parameter: 312 . p - The probability density at x 313 314 Level: beginner 315 316 .seealso: PetscCDFConstant1D(), PetscPDFSampleConstant1D(), PetscPDFConstant2D(), PetscPDFConstant3D() 317 @*/ 318 PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 319 { 320 p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.; 321 return 0; 322 } 323 324 /*@ 325 PetscCDFConstant1D - CDF for the uniform distribution in 1D 326 327 Not collective 328 329 Input Parameter: 330 . x - Coordinate in [-1, 1] 331 332 Output Parameter: 333 . p - The cumulative probability at x 334 335 Level: beginner 336 337 .seealso: PetscPDFConstant1D(), PetscPDFSampleConstant1D(), PetscCDFConstant2D(), PetscCDFConstant3D() 338 @*/ 339 PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 340 { 341 p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5*(x[0] + 1.)); 342 return 0; 343 } 344 345 /*@ 346 PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D 347 348 Not collective 349 350 Input Parameter: 351 . p - A uniform variable on [0, 1] 352 353 Output Parameter: 354 . x - Coordinate in [-1, 1] 355 356 Level: beginner 357 358 .seealso: PetscPDFConstant1D(), PetscCDFConstant1D(), PetscPDFSampleConstant2D(), PetscPDFSampleConstant3D() 359 @*/ 360 PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 361 { 362 x[0] = 2.*p[0] - 1.; 363 return 0; 364 } 365 366 /*@ 367 PetscPDFConstant2D - PDF for the uniform distribution in 2D 368 369 Not collective 370 371 Input Parameter: 372 . x - Coordinate in [-1, 1] x [-1, 1] 373 374 Output Parameter: 375 . p - The probability density at x 376 377 Level: beginner 378 379 .seealso: PetscCDFConstant2D(), PetscPDFSampleConstant2D(), PetscPDFConstant1D(), PetscPDFConstant3D() 380 @*/ 381 PetscErrorCode PetscPDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 382 { 383 p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. ? 0.25 : 0.; 384 return 0; 385 } 386 387 /*@ 388 PetscCDFConstant2D - CDF for the uniform distribution in 2D 389 390 Not collective 391 392 Input Parameter: 393 . x - Coordinate in [-1, 1] x [-1, 1] 394 395 Output Parameter: 396 . p - The cumulative probability at x 397 398 Level: beginner 399 400 .seealso: PetscPDFConstant2D(), PetscPDFSampleConstant2D(), PetscCDFConstant1D(), PetscCDFConstant3D() 401 @*/ 402 PetscErrorCode PetscCDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 403 { 404 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.)); 405 return 0; 406 } 407 408 /*@ 409 PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 2D 410 411 Not collective 412 413 Input Parameter: 414 . p - Two uniform variables on [0, 1] 415 416 Output Parameter: 417 . x - Coordinate in [-1, 1] x [-1, 1] 418 419 Level: beginner 420 421 .seealso: PetscPDFConstant2D(), PetscCDFConstant2D(), PetscPDFSampleConstant1D(), PetscPDFSampleConstant3D() 422 @*/ 423 PetscErrorCode PetscPDFSampleConstant2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 424 { 425 x[0] = 2.*p[0] - 1.; 426 x[1] = 2.*p[1] - 1.; 427 return 0; 428 } 429 430 /*@ 431 PetscPDFConstant3D - PDF for the uniform distribution in 3D 432 433 Not collective 434 435 Input Parameter: 436 . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1] 437 438 Output Parameter: 439 . p - The probability density at x 440 441 Level: beginner 442 443 .seealso: PetscCDFConstant3D(), PetscPDFSampleConstant3D(), PetscPDFSampleConstant1D(), PetscPDFSampleConstant2D() 444 @*/ 445 PetscErrorCode PetscPDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 446 { 447 p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. && x[2] > -1. && x[2] <= 1. ? 0.125 : 0.; 448 return 0; 449 } 450 451 /*@ 452 PetscCDFConstant3D - CDF for the uniform distribution in 3D 453 454 Not collective 455 456 Input Parameter: 457 . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1] 458 459 Output Parameter: 460 . p - The cumulative probability at x 461 462 Level: beginner 463 464 .seealso: PetscPDFConstant3D(), PetscPDFSampleConstant3D(), PetscCDFConstant1D(), PetscCDFConstant2D() 465 @*/ 466 PetscErrorCode PetscCDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[]) 467 { 468 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.)); 469 return 0; 470 } 471 472 /*@ 473 PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 3D 474 475 Not collective 476 477 Input Parameter: 478 . p - Three uniform variables on [0, 1] 479 480 Output Parameter: 481 . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1] 482 483 Level: beginner 484 485 .seealso: PetscPDFConstant3D(), PetscCDFConstant3D(), PetscPDFSampleConstant1D(), PetscPDFSampleConstant2D() 486 @*/ 487 PetscErrorCode PetscPDFSampleConstant3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[]) 488 { 489 x[0] = 2.*p[0] - 1.; 490 x[1] = 2.*p[1] - 1.; 491 x[2] = 2.*p[2] - 1.; 492 return 0; 493 } 494 495 /*@C 496 PetscProbCreateFromOptions - Return the probability distribution specified by the argumetns and options 497 498 Not collective 499 500 Input Parameters: 501 + dim - The dimension of sample points 502 . prefix - The options prefix, or NULL 503 - name - The option name for the probility distribution type 504 505 Output Parameters: 506 + pdf - The PDF of this type 507 . cdf - The CDF of this type 508 - sampler - The PDF sampler of this type 509 510 Level: intermediate 511 512 .seealso: PetscPDFMaxwellBoltzmann1D(), PetscPDFGaussian1D(), PetscPDFConstant1D() 513 @*/ 514 PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler) 515 { 516 DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN; 517 PetscErrorCode ierr; 518 519 PetscFunctionBegin; 520 ierr = PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");PetscCall(ierr); 521 PetscCall(PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum) den, (PetscEnum *) &den, NULL)); 522 ierr = PetscOptionsEnd();PetscCall(ierr); 523 524 if (pdf) {PetscValidPointer(pdf, 4); *pdf = NULL;} 525 if (cdf) {PetscValidPointer(cdf, 5); *cdf = NULL;} 526 if (sampler) {PetscValidPointer(sampler, 6); *sampler = NULL;} 527 switch (den) { 528 case DTPROB_DENSITY_CONSTANT: 529 switch (dim) { 530 case 1: 531 if (pdf) *pdf = PetscPDFConstant1D; 532 if (cdf) *cdf = PetscCDFConstant1D; 533 if (sampler) *sampler = PetscPDFSampleConstant1D; 534 break; 535 case 2: 536 if (pdf) *pdf = PetscPDFConstant2D; 537 if (cdf) *cdf = PetscCDFConstant2D; 538 if (sampler) *sampler = PetscPDFSampleConstant2D; 539 break; 540 case 3: 541 if (pdf) *pdf = PetscPDFConstant3D; 542 if (cdf) *cdf = PetscCDFConstant3D; 543 if (sampler) *sampler = PetscPDFSampleConstant3D; 544 break; 545 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 546 } 547 break; 548 case DTPROB_DENSITY_GAUSSIAN: 549 switch (dim) { 550 case 1: 551 if (pdf) *pdf = PetscPDFGaussian1D; 552 if (cdf) *cdf = PetscCDFGaussian1D; 553 if (sampler) *sampler = PetscPDFSampleGaussian1D; 554 break; 555 case 2: 556 if (pdf) *pdf = PetscPDFGaussian2D; 557 if (sampler) *sampler = PetscPDFSampleGaussian2D; 558 break; 559 case 3: 560 if (pdf) *pdf = PetscPDFGaussian3D; 561 if (sampler) *sampler = PetscPDFSampleGaussian3D; 562 break; 563 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 564 } 565 break; 566 case DTPROB_DENSITY_MAXWELL_BOLTZMANN: 567 switch (dim) { 568 case 1: 569 if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D; 570 if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D; 571 break; 572 case 2: 573 if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D; 574 if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D; 575 break; 576 case 3: 577 if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D; 578 if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D; 579 break; 580 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]); 581 } 582 break; 583 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]); 584 } 585 PetscFunctionReturn(0); 586 } 587 588 #ifdef PETSC_HAVE_KS 589 EXTERN_C_BEGIN 590 #include <KolmogorovSmirnovDist.h> 591 EXTERN_C_END 592 #endif 593 594 /*@C 595 PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF. 596 597 Collective on v 598 599 Input Parameters: 600 + v - The data vector, blocksize is the sample dimension 601 - cdf - The analytic CDF 602 603 Output Parameter: 604 . alpha - The KS statisic 605 606 Level: advanced 607 608 Note: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is 609 $ 610 $ D_n = \sup_x \left| F_n(x) - F(x) \right| 611 $ 612 where $\sup_x$ is the supremum of the set of distances, and the empirical distribution 613 function $F_n(x)$ is discrete, and given by 614 $ 615 $ F_n = # of samples <= x / n 616 $ 617 The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep'' 618 cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes 619 the largest absolute difference between the two distribution functions across all $x$ values. 620 621 .seealso: PetscProbFunc 622 @*/ 623 PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha) 624 { 625 #if !defined(PETSC_HAVE_KS) 626 SETERRQ(PetscObjectComm((PetscObject) v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks"); 627 #else 628 PetscViewer viewer = NULL; 629 PetscViewerFormat format; 630 PetscDraw draw; 631 PetscDrawLG lg; 632 PetscDrawAxis axis; 633 const PetscScalar *a; 634 PetscReal *speed, Dn = PETSC_MIN_REAL; 635 PetscBool isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg; 636 PetscInt n, p, dim, d; 637 PetscMPIInt size; 638 const char *names[2] = {"Analytic", "Empirical"}; 639 char title[PETSC_MAX_PATH_LEN]; 640 PetscOptions options; 641 const char *prefix; 642 MPI_Comm comm; 643 644 PetscFunctionBegin; 645 PetscCall(PetscObjectGetComm((PetscObject) v, &comm)); 646 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) v, &prefix)); 647 PetscCall(PetscObjectGetOptions((PetscObject) v, &options)); 648 PetscCall(PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg)); 649 if (flg) { 650 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii)); 651 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw)); 652 } 653 if (isascii) { 654 PetscCall(PetscViewerPushFormat(viewer, format)); 655 } else if (isdraw) { 656 PetscCall(PetscViewerPushFormat(viewer, format)); 657 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 658 PetscCall(PetscDrawLGCreate(draw, 2, &lg)); 659 PetscCall(PetscDrawLGSetLegend(lg, names)); 660 } 661 662 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) v), &size)); 663 PetscCall(VecGetLocalSize(v, &n)); 664 PetscCall(VecGetBlockSize(v, &dim)); 665 n /= dim; 666 /* TODO Parallel algorithm is harder */ 667 if (size == 1) { 668 PetscCall(PetscMalloc1(n, &speed)); 669 PetscCall(VecGetArrayRead(v, &a)); 670 for (p = 0; p < n; ++p) { 671 PetscReal mag = 0.; 672 673 for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p*dim+d])); 674 speed[p] = PetscSqrtReal(mag); 675 } 676 PetscCall(PetscSortReal(n, speed)); 677 PetscCall(VecRestoreArrayRead(v, &a)); 678 for (p = 0; p < n; ++p) { 679 const PetscReal x = speed[p], Fn = ((PetscReal) p) / n; 680 PetscReal F, vals[2]; 681 682 PetscCall(cdf(&x, NULL, &F)); 683 Dn = PetscMax(PetscAbsReal(Fn - F), Dn); 684 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn)); 685 if (isdraw) {vals[0] = F; vals[1] = Fn; PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));} 686 } 687 if (speed[n-1] < 6.) { 688 const PetscReal k = (PetscInt) (6. - speed[n-1]) + 1, dx = (6. - speed[n-1])/k; 689 for (p = 0; p <= k; ++p) { 690 const PetscReal x = speed[n-1] + p*dx, Fn = 1.0; 691 PetscReal F, vals[2]; 692 693 PetscCall(cdf(&x, NULL, &F)); 694 Dn = PetscMax(PetscAbsReal(Fn - F), Dn); 695 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn)); 696 if (isdraw) {vals[0] = F; vals[1] = Fn; PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));} 697 } 698 } 699 PetscCall(PetscFree(speed)); 700 } 701 if (isdraw) { 702 PetscCall(PetscDrawLGGetAxis(lg, &axis)); 703 PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn)); 704 PetscCall(PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)")); 705 PetscCall(PetscDrawLGDraw(lg)); 706 PetscCall(PetscDrawLGDestroy(&lg)); 707 } 708 if (viewer) { 709 PetscCall(PetscViewerFlush(viewer)); 710 PetscCall(PetscViewerPopFormat(viewer)); 711 PetscCall(PetscViewerDestroy(&viewer)); 712 } 713 *alpha = KSfbar((int) n, (double) Dn); 714 PetscFunctionReturn(0); 715 #endif 716 } 717