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