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 Parameters:
 15: + x     - Speed in [0, \infty]
 16: - dummy - Unused

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

 21:   Level: beginner

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

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

 34:   Not Collective

 36:   Input Parameters:
 37: + x     - Speed in [0, \infty]
 38: - dummy - Unused

 40:   Output Parameter:
 41: . p - The probability density at x

 43:   Level: beginner

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

 53: /*@
 54:   PetscPDFMaxwellBoltzmann2D - PDF for the Maxwell-Boltzmann distribution in 2D

 56:   Not Collective

 58:   Input Parameters:
 59: + x     - Speed in [0, \infty]
 60: - dummy - Unused

 62:   Output Parameter:
 63: . p - The probability density at x

 65:   Level: beginner

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

 75: /*@
 76:   PetscCDFMaxwellBoltzmann2D - CDF for the Maxwell-Boltzmann distribution in 2D

 78:   Not Collective

 80:   Input Parameters:
 81: + x     - Speed in [0, \infty]
 82: - dummy - Unused

 84:   Output Parameter:
 85: . p - The probability density at x

 87:   Level: beginner

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

 97: /*@
 98:   PetscPDFMaxwellBoltzmann3D - PDF for the Maxwell-Boltzmann distribution in 3D

100:   Not Collective

102:   Input Parameters:
103: + x     - Speed in [0, \infty]
104: - dummy - Unused

106:   Output Parameter:
107: . p - The probability density at x

109:   Level: beginner

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

119: /*@
120:   PetscCDFMaxwellBoltzmann3D - CDF for the Maxwell-Boltzmann distribution in 3D

122:   Not Collective

124:   Input Parameters:
125: + x     - Speed in [0, \infty]
126: - dummy - Unused

128:   Output Parameter:
129: . p - The probability density at x

131:   Level: beginner

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

141: /*@
142:   PetscPDFGaussian1D - PDF for the Gaussian distribution in 1D

144:   Not Collective

146:   Input Parameters:
147: + x     - Coordinate in [-\infty, \infty]
148: - scale - Scaling value

150:   Output Parameter:
151: . p - The probability density at x

153:   Level: beginner

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

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

171: /*@
172:   PetscPDFSampleGaussian1D - Sample uniformly from a Gaussian distribution in 1D

174:   Not Collective

176:   Input Parameters:
177: + p     - A uniform variable on [0, 1]
178: - dummy - ignored

180:   Output Parameter:
181: . x - Coordinate in [-\infty, \infty]

183:   Level: beginner

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

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;

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

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

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

218: /*@
219:   PetscPDFGaussian2D - PDF for the Gaussian distribution in 2D

221:   Not Collective

223:   Input Parameters:
224: + x     - Coordinate in [-\infty, \infty]^2
225: - dummy - ignored

227:   Output Parameter:
228: . p - The probability density at x

230:   Level: beginner

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

240: /*@
241:   PetscPDFSampleGaussian2D - Sample uniformly from a Gaussian distribution in 2D

243:   Not Collective

245:   Input Parameters:
246: + p     - A uniform variable on [0, 1]^2
247: - dummy - ignored

249:   Output Parameter:
250: . x - Coordinate in [-\infty, \infty]^2

252:   Level: beginner

254:   References:
255: . * - https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

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

267: /*@
268:   PetscPDFGaussian3D - PDF for the Gaussian distribution in 3D

270:   Not Collective

272:   Input Parameters:
273: + x     - Coordinate in [-\infty, \infty]^3
274: - dummy - ignored

276:   Output Parameter:
277: . p - The probability density at x

279:   Level: beginner

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

289: /*@
290:   PetscPDFSampleGaussian3D - Sample uniformly from a Gaussian distribution in 3D

292:   Not Collective

294:   Input Parameters:
295: + p     - A uniform variable on [0, 1]^3
296: - dummy - ignored

298:   Output Parameter:
299: . x - Coordinate in [-\infty, \infty]^3

301:   Level: beginner

303:   References:
304:   https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

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

315: /*@
316:   PetscPDFConstant1D - PDF for the uniform distribution in 1D

318:   Not Collective

320:   Input Parameters:
321: + x     - Coordinate in [-1, 1]
322: - dummy - Unused

324:   Output Parameter:
325: . p - The probability density at x

327:   Level: beginner

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

337: /*@
338:   PetscCDFConstant1D - CDF for the uniform distribution in 1D

340:   Not Collective

342:   Input Parameters:
343: + x     - Coordinate in [-1, 1]
344: - dummy - Unused

346:   Output Parameter:
347: . p - The cumulative probability at x

349:   Level: beginner

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

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

362:   Not Collective

364:   Input Parameters:
365: + p     - A uniform variable on [0, 1]
366: - dummy - Unused

368:   Output Parameter:
369: . x - Coordinate in [-1, 1]

371:   Level: beginner

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

381: /*@
382:   PetscPDFConstant2D - PDF for the uniform distribution in 2D

384:   Not Collective

386:   Input Parameters:
387: + x     - Coordinate in [-1, 1] x [-1, 1]
388: - dummy - Unused

390:   Output Parameter:
391: . p - The probability density at x

393:   Level: beginner

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

403: /*@
404:   PetscCDFConstant2D - CDF for the uniform distribution in 2D

406:   Not Collective

408:   Input Parameters:
409: + x     - Coordinate in [-1, 1] x [-1, 1]
410: - dummy - Unused

412:   Output Parameter:
413: . p - The cumulative probability at x

415:   Level: beginner

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

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

428:   Not Collective

430:   Input Parameters:
431: + p     - Two uniform variables on [0, 1]
432: - dummy - Unused

434:   Output Parameter:
435: . x - Coordinate in [-1, 1] x [-1, 1]

437:   Level: beginner

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

448: /*@
449:   PetscPDFConstant3D - PDF for the uniform distribution in 3D

451:   Not Collective

453:   Input Parameters:
454: + x     - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
455: - dummy - Unused

457:   Output Parameter:
458: . p - The probability density at x

460:   Level: beginner

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

470: /*@
471:   PetscCDFConstant3D - CDF for the uniform distribution in 3D

473:   Not Collective

475:   Input Parameters:
476: + x     - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
477: - dummy - Unused

479:   Output Parameter:
480: . p - The cumulative probability at x

482:   Level: beginner

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

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

495:   Not Collective

497:   Input Parameters:
498: + p     - Three uniform variables on [0, 1]
499: - dummy - Unused

501:   Output Parameter:
502: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]

504:   Level: beginner

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

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

519:   Not Collective

521:   Input Parameters:
522: + dim    - The dimension of sample points
523: . prefix - The options prefix, or NULL
524: - name   - The option name for the probility distribution type

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

531:   Level: intermediate

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;

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

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

621: #ifdef PETSC_HAVE_KS
622: EXTERN_C_BEGIN
623:   #include <KolmogorovSmirnovDist.h>
624: EXTERN_C_END
625: #endif

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

630:   Collective

632:   Input Parameters:
633: + v   - The data vector, blocksize is the sample dimension
634: - cdf - The analytic CDF

636:   Output Parameter:
637: . alpha - The KS statisic

639:   Level: advanced

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|

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

648:     F_n =  # of samples <= x / n
649: .ve

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.

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;

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

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

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

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

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