Actual source code: dtprob.c

  1: #include <petscdt.h>

  3: #include <petscvec.h>
  4: #include <petscdraw.h>
  5: #include <petsc/private/petscimpl.h>

  7: const char *const DTProbDensityTypes[] = {"constant", "gaussian", "maxwell_boltzmann", "DTProb Density", "DTPROB_DENSITY_", NULL};

  9: /*@
 10:   PetscPDFMaxwellBoltzmann1D - PDF for the Maxwell-Boltzmann distribution in 1D

 12:   Not Collective

 14:   Input Parameter:
 15: . x - Speed in [0, \infty]

 17:   Output Parameter:
 18: . p - The probability density at x

 20:   Level: beginner

 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: }

 30: /*@
 31:   PetscCDFMaxwellBoltzmann1D - CDF for the Maxwell-Boltzmann distribution in 1D

 33:   Not Collective

 35:   Input Parameter:
 36: . x - Speed in [0, \infty]

 38:   Output Parameter:
 39: . p - The probability density at x

 41:   Level: beginner

 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: }

 51: /*@
 52:   PetscPDFMaxwellBoltzmann2D - PDF for the Maxwell-Boltzmann distribution in 2D

 54:   Not Collective

 56:   Input Parameter:
 57: . x - Speed in [0, \infty]

 59:   Output Parameter:
 60: . p - The probability density at x

 62:   Level: beginner

 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: }

 72: /*@
 73:   PetscCDFMaxwellBoltzmann2D - CDF for the Maxwell-Boltzmann distribution in 2D

 75:   Not Collective

 77:   Input Parameter:
 78: . x - Speed in [0, \infty]

 80:   Output Parameter:
 81: . p - The probability density at x

 83:   Level: beginner

 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: }

 93: /*@
 94:   PetscPDFMaxwellBoltzmann3D - PDF for the Maxwell-Boltzmann distribution in 3D

 96:   Not Collective

 98:   Input Parameter:
 99: . x - Speed in [0, \infty]

101:   Output Parameter:
102: . p - The probability density at x

104:   Level: beginner

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: }

114: /*@
115:   PetscCDFMaxwellBoltzmann3D - CDF for the Maxwell-Boltzmann distribution in 3D

117:   Not Collective

119:   Input Parameter:
120: . x - Speed in [0, \infty]

122:   Output Parameter:
123: . p - The probability density at x

125:   Level: beginner

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: }

135: /*@
136:   PetscPDFGaussian1D - PDF for the Gaussian distribution in 1D

138:   Not Collective

140:   Input Parameter:
141: . x - Coordinate in [-\infty, \infty]

143:   Output Parameter:
144: . p - The probability density at x

146:   Level: beginner

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: }

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: }

164: /*@
165:   PetscPDFSampleGaussian1D - Sample uniformly from a Gaussian distribution in 1D

167:   Not Collective

169:   Input Parameters:
170: + p - A uniform variable on [0, 1]
171: - dummy - ignored

173:   Output Parameter:
174: . x - Coordinate in [-\infty, \infty]

176:   Level: beginner

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;

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.;

196:     for (m = 0; m <= k - 1; ++m) {
197:       PetscReal denom = (m + 1.) * (2. * m + 1.);

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: }

208: /*@
209:   PetscPDFGaussian2D - PDF for the Gaussian distribution in 2D

211:   Not Collective

213:   Input Parameters:
214: + x - Coordinate in [-\infty, \infty]^2
215: - dummy - ignored

217:   Output Parameter:
218: . p - The probability density at x

220:   Level: beginner

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: }

230: /*@
231:   PetscPDFSampleGaussian2D - Sample uniformly from a Gaussian distribution in 2D

233:   Not Collective

235:   Input Parameters:
236: + p - A uniform variable on [0, 1]^2
237: - dummy - ignored

239:   Output Parameter:
240: . x - Coordinate in [-\infty, \infty]^2

242:   Level: beginner

244:   Note: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

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: }

256: /*@
257:   PetscPDFGaussian3D - PDF for the Gaussian distribution in 3D

259:   Not Collective

261:   Input Parameters:
262: + x - Coordinate in [-\infty, \infty]^3
263: - dummy - ignored

265:   Output Parameter:
266: . p - The probability density at x

268:   Level: beginner

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: }

278: /*@
279:   PetscPDFSampleGaussian3D - Sample uniformly from a Gaussian distribution in 3D

281:   Not Collective

283:   Input Parameters:
284: + p - A uniform variable on [0, 1]^3
285: - dummy - ignored

287:   Output Parameter:
288: . x - Coordinate in [-\infty, \infty]^3

290:   Level: beginner

292:   Reference:
293:   https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

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: }

304: /*@
305:   PetscPDFConstant1D - PDF for the uniform distribution in 1D

307:   Not Collective

309:   Input Parameter:
310: . x - Coordinate in [-1, 1]

312:   Output Parameter:
313: . p - The probability density at x

315:   Level: beginner

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: }

325: /*@
326:   PetscCDFConstant1D - CDF for the uniform distribution in 1D

328:   Not Collective

330:   Input Parameter:
331: . x - Coordinate in [-1, 1]

333:   Output Parameter:
334: . p - The cumulative probability at x

336:   Level: beginner

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: }

346: /*@
347:   PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D

349:   Not Collective

351:   Input Parameter:
352: . p - A uniform variable on [0, 1]

354:   Output Parameter:
355: . x - Coordinate in [-1, 1]

357:   Level: beginner

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: }

367: /*@
368:   PetscPDFConstant2D - PDF for the uniform distribution in 2D

370:   Not Collective

372:   Input Parameter:
373: . x - Coordinate in [-1, 1] x [-1, 1]

375:   Output Parameter:
376: . p - The probability density at x

378:   Level: beginner

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: }

388: /*@
389:   PetscCDFConstant2D - CDF for the uniform distribution in 2D

391:   Not Collective

393:   Input Parameter:
394: . x - Coordinate in [-1, 1] x [-1, 1]

396:   Output Parameter:
397: . p - The cumulative probability at x

399:   Level: beginner

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: }

409: /*@
410:   PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 2D

412:   Not Collective

414:   Input Parameter:
415: . p - Two uniform variables on [0, 1]

417:   Output Parameter:
418: . x - Coordinate in [-1, 1] x [-1, 1]

420:   Level: beginner

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: }

431: /*@
432:   PetscPDFConstant3D - PDF for the uniform distribution in 3D

434:   Not Collective

436:   Input Parameter:
437: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]

439:   Output Parameter:
440: . p - The probability density at x

442:   Level: beginner

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: }

452: /*@
453:   PetscCDFConstant3D - CDF for the uniform distribution in 3D

455:   Not Collective

457:   Input Parameter:
458: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]

460:   Output Parameter:
461: . p - The cumulative probability at x

463:   Level: beginner

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: }

473: /*@
474:   PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 3D

476:   Not Collective

478:   Input Parameter:
479: . p - Three uniform variables on [0, 1]

481:   Output Parameter:
482: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]

484:   Level: beginner

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: }

496: /*@C
497:   PetscProbCreateFromOptions - Return the probability distribution specified by the arguments and options

499:   Not Collective

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

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

511:   Level: intermediate

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;

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();

524:   if (pdf) {
526:     *pdf = NULL;
527:   }
528:   if (cdf) {
530:     *cdf = NULL;
531:   }
532:   if (sampler) {
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: }

601: #ifdef PETSC_HAVE_KS
602: EXTERN_C_BEGIN
603:   #include <KolmogorovSmirnovDist.h>
604: EXTERN_C_END
605: #endif

607: /*@C
608:   PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.

610:   Collective

612:   Input Parameters:
613: + v   - The data vector, blocksize is the sample dimension
614: - cdf - The analytic CDF

616:   Output Parameter:
617: . alpha - The KS statisic

619:   Level: advanced

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|

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

628:     F_n =  # of samples <= x / n
629: .ve

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.

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;

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:   }

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.;

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];

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];

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: }