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

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:   return 0;
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 0;
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 0;
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 0;
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:   Note: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

294: .seealso: `PetscPDFGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
295: @*/
296: PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
297: {
298:   PetscPDFSampleGaussian1D(p, dummy, x);
299:   PetscPDFSampleGaussian2D(&p[1], dummy, &x[1]);
300:   return 0;
301: }

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

306:   Not collective

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

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

314:   Level: beginner

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

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

327:   Not collective

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

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

335:   Level: beginner

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

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

348:   Not collective

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

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

356:   Level: beginner

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

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

369:   Not collective

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

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

377:   Level: beginner

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

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

390:   Not collective

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

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

398:   Level: beginner

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

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

411:   Not collective

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

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

419:   Level: beginner

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

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

433:   Not collective

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

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

441:   Level: beginner

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

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

454:   Not collective

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

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

462:   Level: beginner

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

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

475:   Not collective

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

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

483:   Level: beginner

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

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

498:   Not collective

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

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

510:   Level: intermediate

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;

518:   PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");
519:   PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum)den, (PetscEnum *)&den, NULL);
520:   PetscOptionsEnd();

522:   if (pdf) {
524:     *pdf = NULL;
525:   }
526:   if (cdf) {
528:     *cdf = NULL;
529:   }
530:   if (sampler) {
532:     *sampler = NULL;
533:   }
534:   switch (den) {
535:   case DTPROB_DENSITY_CONSTANT:
536:     switch (dim) {
537:     case 1:
538:       if (pdf) *pdf = PetscPDFConstant1D;
539:       if (cdf) *cdf = PetscCDFConstant1D;
540:       if (sampler) *sampler = PetscPDFSampleConstant1D;
541:       break;
542:     case 2:
543:       if (pdf) *pdf = PetscPDFConstant2D;
544:       if (cdf) *cdf = PetscCDFConstant2D;
545:       if (sampler) *sampler = PetscPDFSampleConstant2D;
546:       break;
547:     case 3:
548:       if (pdf) *pdf = PetscPDFConstant3D;
549:       if (cdf) *cdf = PetscCDFConstant3D;
550:       if (sampler) *sampler = PetscPDFSampleConstant3D;
551:       break;
552:     default:
553:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
554:     }
555:     break;
556:   case DTPROB_DENSITY_GAUSSIAN:
557:     switch (dim) {
558:     case 1:
559:       if (pdf) *pdf = PetscPDFGaussian1D;
560:       if (cdf) *cdf = PetscCDFGaussian1D;
561:       if (sampler) *sampler = PetscPDFSampleGaussian1D;
562:       break;
563:     case 2:
564:       if (pdf) *pdf = PetscPDFGaussian2D;
565:       if (sampler) *sampler = PetscPDFSampleGaussian2D;
566:       break;
567:     case 3:
568:       if (pdf) *pdf = PetscPDFGaussian3D;
569:       if (sampler) *sampler = PetscPDFSampleGaussian3D;
570:       break;
571:     default:
572:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
573:     }
574:     break;
575:   case DTPROB_DENSITY_MAXWELL_BOLTZMANN:
576:     switch (dim) {
577:     case 1:
578:       if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D;
579:       if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D;
580:       break;
581:     case 2:
582:       if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D;
583:       if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D;
584:       break;
585:     case 3:
586:       if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D;
587:       if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D;
588:       break;
589:     default:
590:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
591:     }
592:     break;
593:   default:
594:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]);
595:   }
596:   return 0;
597: }

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

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

608:   Collective on v

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

614:   Output Parameter:
615: . alpha - The KS statisic

617:   Level: advanced

619:   Note: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
620: $
621: $    D_n = \sup_x \left| F_n(x) - F(x) \right|
622: $
623: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution
624: function $F_n(x)$ is discrete, and given by
625: $
626: $    F_n =  # of samples <= x / n
627: $
628: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
629: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
630: the largest absolute difference between the two distribution functions across all $x$ values.

632: .seealso: `PetscProbFunc`
633: @*/
634: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha)
635: {
636: #if !defined(PETSC_HAVE_KS)
637:   SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
638: #else
639:   PetscViewer        viewer = NULL;
640:   PetscViewerFormat  format;
641:   PetscDraw          draw;
642:   PetscDrawLG        lg;
643:   PetscDrawAxis      axis;
644:   const PetscScalar *a;
645:   PetscReal         *speed, Dn = PETSC_MIN_REAL;
646:   PetscBool          isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg;
647:   PetscInt           n, p, dim, d;
648:   PetscMPIInt        size;
649:   const char        *names[2] = {"Analytic", "Empirical"};
650:   char               title[PETSC_MAX_PATH_LEN];
651:   PetscOptions       options;
652:   const char        *prefix;
653:   MPI_Comm           comm;

655:   PetscObjectGetComm((PetscObject)v, &comm);
656:   PetscObjectGetOptionsPrefix((PetscObject)v, &prefix);
657:   PetscObjectGetOptions((PetscObject)v, &options);
658:   PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg);
659:   if (flg) {
660:     PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
661:     PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
662:   }
663:   if (isascii) {
664:     PetscViewerPushFormat(viewer, format);
665:   } else if (isdraw) {
666:     PetscViewerPushFormat(viewer, format);
667:     PetscViewerDrawGetDraw(viewer, 0, &draw);
668:     PetscDrawLGCreate(draw, 2, &lg);
669:     PetscDrawLGSetLegend(lg, names);
670:   }

672:   MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
673:   VecGetLocalSize(v, &n);
674:   VecGetBlockSize(v, &dim);
675:   n /= dim;
676:   /* TODO Parallel algorithm is harder */
677:   if (size == 1) {
678:     PetscMalloc1(n, &speed);
679:     VecGetArrayRead(v, &a);
680:     for (p = 0; p < n; ++p) {
681:       PetscReal mag = 0.;

683:       for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d]));
684:       speed[p] = PetscSqrtReal(mag);
685:     }
686:     PetscSortReal(n, speed);
687:     VecRestoreArrayRead(v, &a);
688:     for (p = 0; p < n; ++p) {
689:       const PetscReal x = speed[p], Fn = ((PetscReal)p) / n;
690:       PetscReal       F, vals[2];

692:       cdf(&x, NULL, &F);
693:       Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
694:       if (isascii) PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);
695:       if (isdraw) {
696:         vals[0] = F;
697:         vals[1] = Fn;
698:         PetscDrawLGAddCommonPoint(lg, x, vals);
699:       }
700:     }
701:     if (speed[n - 1] < 6.) {
702:       const PetscReal k = (PetscInt)(6. - speed[n - 1]) + 1, dx = (6. - speed[n - 1]) / k;
703:       for (p = 0; p <= k; ++p) {
704:         const PetscReal x = speed[n - 1] + p * dx, Fn = 1.0;
705:         PetscReal       F, vals[2];

707:         cdf(&x, NULL, &F);
708:         Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
709:         if (isascii) PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);
710:         if (isdraw) {
711:           vals[0] = F;
712:           vals[1] = Fn;
713:           PetscDrawLGAddCommonPoint(lg, x, vals);
714:         }
715:       }
716:     }
717:     PetscFree(speed);
718:   }
719:   if (isdraw) {
720:     PetscDrawLGGetAxis(lg, &axis);
721:     PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn);
722:     PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)");
723:     PetscDrawLGDraw(lg);
724:     PetscDrawLGDestroy(&lg);
725:   }
726:   if (viewer) {
727:     PetscViewerFlush(viewer);
728:     PetscViewerPopFormat(viewer);
729:     PetscViewerDestroy(&viewer);
730:   }
731:   *alpha = KSfbar((int)n, (double)Dn);
732:   return 0;
733: #endif
734: }