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: Note:
186: See <http://www.mimirgames.com/articles/programming/approximations-of-the-inverse-error-function> and
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
242: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
244: Not Collective
246: Input Parameters:
247: + p - A uniform variable on $[0, 1]^2$
248: - dummy - ignored
250: Output Parameter:
251: . x - Coordinate in $[-\infty, \infty]^2 $
253: Level: beginner
255: .seealso: `PetscPDFGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()`
256: @*/
257: PetscErrorCode PetscPDFSampleGaussian2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
258: {
259: const PetscReal mag = PetscSqrtReal(-2.0 * PetscLogReal(p[0]));
260: x[0] = mag * PetscCosReal(2.0 * PETSC_PI * p[1]);
261: x[1] = mag * PetscSinReal(2.0 * PETSC_PI * p[1]);
262: return PETSC_SUCCESS;
263: }
265: /*@
266: PetscPDFGaussian3D - PDF for the Gaussian distribution in 3D
268: Not Collective
270: Input Parameters:
271: + x - Coordinate in $[-\infty, \infty]^3$
272: - dummy - ignored
274: Output Parameter:
275: . p - The probability density at `x`
277: Level: beginner
279: .seealso: `PetscPDFSampleGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
280: @*/
281: PetscErrorCode PetscPDFGaussian3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
282: {
283: p[0] = (1. / PETSC_PI * PetscSqrtReal(PETSC_PI)) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1]) + PetscSqr(x[2])));
284: return PETSC_SUCCESS;
285: }
287: /*@
288: PetscPDFSampleGaussian3D - Sample uniformly from a Gaussian distribution in 3D
289: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
291: Not Collective
293: Input Parameters:
294: + p - A uniform variable on $[0, 1]^3$
295: - dummy - ignored
297: Output Parameter:
298: . x - Coordinate in $[-\infty, \infty]^3$
300: Level: beginner
302: .seealso: `PetscPDFGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
303: @*/
304: PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
305: {
306: PetscCall(PetscPDFSampleGaussian1D(p, dummy, x));
307: PetscCall(PetscPDFSampleGaussian2D(&p[1], dummy, &x[1]));
308: return PETSC_SUCCESS;
309: }
311: /*@
312: PetscPDFConstant1D - PDF for the uniform distribution in 1D
314: Not Collective
316: Input Parameters:
317: + x - Coordinate in $[-1, 1]$
318: - dummy - Unused
320: Output Parameter:
321: . p - The probability density at `x`
323: Level: beginner
325: .seealso: `PetscCDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscPDFConstant2D()`, `PetscPDFConstant3D()`
326: @*/
327: PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
328: {
329: p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.;
330: return PETSC_SUCCESS;
331: }
333: /*@
334: PetscCDFConstant1D - CDF for the uniform distribution in 1D
336: Not Collective
338: Input Parameters:
339: + x - Coordinate in $[-1, 1]$
340: - dummy - Unused
342: Output Parameter:
343: . p - The cumulative probability at `x`
345: Level: beginner
347: .seealso: `PetscPDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscCDFConstant2D()`, `PetscCDFConstant3D()`
348: @*/
349: PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
350: {
351: p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.));
352: return PETSC_SUCCESS;
353: }
355: /*@
356: PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D
358: Not Collective
360: Input Parameters:
361: + p - A uniform variable on $[0, 1]$
362: - dummy - Unused
364: Output Parameter:
365: . x - Coordinate in $[-1, 1]$
367: Level: beginner
369: .seealso: `PetscPDFConstant1D()`, `PetscCDFConstant1D()`, `PetscPDFSampleConstant2D()`, `PetscPDFSampleConstant3D()`
370: @*/
371: PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
372: {
373: x[0] = 2. * p[0] - 1.;
374: return PETSC_SUCCESS;
375: }
377: /*@
378: PetscPDFConstant2D - PDF for the uniform distribution in 2D
380: Not Collective
382: Input Parameters:
383: + x - Coordinate in $[-1, 1]^2$
384: - dummy - Unused
386: Output Parameter:
387: . p - The probability density at `x`
389: Level: beginner
391: .seealso: `PetscCDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscPDFConstant1D()`, `PetscPDFConstant3D()`
392: @*/
393: PetscErrorCode PetscPDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
394: {
395: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. ? 0.25 : 0.;
396: return PETSC_SUCCESS;
397: }
399: /*@
400: PetscCDFConstant2D - CDF for the uniform distribution in 2D
402: Not Collective
404: Input Parameters:
405: + x - Coordinate in $[-1, 1]^2$
406: - dummy - Unused
408: Output Parameter:
409: . p - The cumulative probability at `x`
411: Level: beginner
413: .seealso: `PetscPDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscCDFConstant1D()`, `PetscCDFConstant3D()`
414: @*/
415: PetscErrorCode PetscCDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
416: {
417: 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.));
418: return PETSC_SUCCESS;
419: }
421: /*@
422: PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on $[-1, 1]^2$ in 2D
424: Not Collective
426: Input Parameters:
427: + p - Two uniform variables on $[0, 1]$
428: - dummy - Unused
430: Output Parameter:
431: . x - Coordinate in $[-1, 1]^2$
433: Level: beginner
435: .seealso: `PetscPDFConstant2D()`, `PetscCDFConstant2D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant3D()`
436: @*/
437: PetscErrorCode PetscPDFSampleConstant2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
438: {
439: x[0] = 2. * p[0] - 1.;
440: x[1] = 2. * p[1] - 1.;
441: return PETSC_SUCCESS;
442: }
444: /*@
445: PetscPDFConstant3D - PDF for the uniform distribution in 3D
447: Not Collective
449: Input Parameters:
450: + x - Coordinate in $[-1, 1]^3$
451: - dummy - Unused
453: Output Parameter:
454: . p - The probability density at `x`
456: Level: beginner
458: .seealso: `PetscCDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
459: @*/
460: PetscErrorCode PetscPDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
461: {
462: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. && x[2] > -1. && x[2] <= 1. ? 0.125 : 0.;
463: return PETSC_SUCCESS;
464: }
466: /*@
467: PetscCDFConstant3D - CDF for the uniform distribution in 3D
469: Not Collective
471: Input Parameters:
472: + x - Coordinate in $[-1, 1]^3$
473: - dummy - Unused
475: Output Parameter:
476: . p - The cumulative probability at `x`
478: Level: beginner
480: .seealso: `PetscPDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscCDFConstant1D()`, `PetscCDFConstant2D()`
481: @*/
482: PetscErrorCode PetscCDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
483: {
484: 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.));
485: return PETSC_SUCCESS;
486: }
488: /*@
489: PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on $[-1, 1]^3$ in 3D
491: Not Collective
493: Input Parameters:
494: + p - Three uniform variables on $[0, 1]$
495: - dummy - Unused
497: Output Parameter:
498: . x - Coordinate in $[-1, 1]^3$
500: Level: beginner
502: .seealso: `PetscPDFConstant3D()`, `PetscCDFConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
503: @*/
504: PetscErrorCode PetscPDFSampleConstant3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
505: {
506: x[0] = 2. * p[0] - 1.;
507: x[1] = 2. * p[1] - 1.;
508: x[2] = 2. * p[2] - 1.;
509: return PETSC_SUCCESS;
510: }
512: /*@C
513: PetscProbCreateFromOptions - Return the probability distribution specified by the arguments and options
515: Not Collective
517: Input Parameters:
518: + dim - The dimension of sample points
519: . prefix - The options prefix, or `NULL`
520: - name - The option name for the probability distribution type
522: Output Parameters:
523: + pdf - The PDF of this type
524: . cdf - The CDF of this type
525: - sampler - The PDF sampler of this type
527: Level: intermediate
529: .seealso: `PetscProbFunc`, `PetscPDFMaxwellBoltzmann1D()`, `PetscPDFGaussian1D()`, `PetscPDFConstant1D()`
530: @*/
531: PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler)
532: {
533: DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN;
535: PetscFunctionBegin;
536: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");
537: PetscCall(PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum)den, (PetscEnum *)&den, NULL));
538: PetscOptionsEnd();
540: if (pdf) {
541: PetscAssertPointer(pdf, 4);
542: *pdf = NULL;
543: }
544: if (cdf) {
545: PetscAssertPointer(cdf, 5);
546: *cdf = NULL;
547: }
548: if (sampler) {
549: PetscAssertPointer(sampler, 6);
550: *sampler = NULL;
551: }
552: switch (den) {
553: case DTPROB_DENSITY_CONSTANT:
554: switch (dim) {
555: case 1:
556: if (pdf) *pdf = PetscPDFConstant1D;
557: if (cdf) *cdf = PetscCDFConstant1D;
558: if (sampler) *sampler = PetscPDFSampleConstant1D;
559: break;
560: case 2:
561: if (pdf) *pdf = PetscPDFConstant2D;
562: if (cdf) *cdf = PetscCDFConstant2D;
563: if (sampler) *sampler = PetscPDFSampleConstant2D;
564: break;
565: case 3:
566: if (pdf) *pdf = PetscPDFConstant3D;
567: if (cdf) *cdf = PetscCDFConstant3D;
568: if (sampler) *sampler = PetscPDFSampleConstant3D;
569: break;
570: default:
571: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
572: }
573: break;
574: case DTPROB_DENSITY_GAUSSIAN:
575: switch (dim) {
576: case 1:
577: if (pdf) *pdf = PetscPDFGaussian1D;
578: if (cdf) *cdf = PetscCDFGaussian1D;
579: if (sampler) *sampler = PetscPDFSampleGaussian1D;
580: break;
581: case 2:
582: if (pdf) *pdf = PetscPDFGaussian2D;
583: if (sampler) *sampler = PetscPDFSampleGaussian2D;
584: break;
585: case 3:
586: if (pdf) *pdf = PetscPDFGaussian3D;
587: if (sampler) *sampler = PetscPDFSampleGaussian3D;
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: case DTPROB_DENSITY_MAXWELL_BOLTZMANN:
594: switch (dim) {
595: case 1:
596: if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D;
597: if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D;
598: break;
599: case 2:
600: if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D;
601: if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D;
602: break;
603: case 3:
604: if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D;
605: if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D;
606: break;
607: default:
608: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
609: }
610: break;
611: default:
612: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]);
613: }
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: #ifdef PETSC_HAVE_KS
618: EXTERN_C_BEGIN
619: #include <KolmogorovSmirnovDist.h>
620: EXTERN_C_END
621: #endif
623: /*@C
624: PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.
626: Collective
628: Input Parameters:
629: + v - The data vector, blocksize is the sample dimension
630: - cdf - The analytic CDF
632: Output Parameter:
633: . alpha - The KS statisic
635: Level: advanced
637: Notes:
638: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
640: $$
641: D_n = \sup_x \left| F_n(x) - F(x) \right|
642: $$
644: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by
646: $$
647: F_n = # of samples <= x / n
648: $$
650: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
651: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
652: the largest absolute difference between the two distribution functions across all $x$ values.
654: .seealso: `PetscProbFunc`
655: @*/
656: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha)
657: {
658: #if !defined(PETSC_HAVE_KS)
659: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
660: #else
661: PetscViewer viewer = NULL;
662: PetscViewerFormat format;
663: PetscDraw draw;
664: PetscDrawLG lg;
665: PetscDrawAxis axis;
666: const PetscScalar *a;
667: PetscReal *speed, Dn = PETSC_MIN_REAL;
668: PetscBool isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg;
669: PetscInt n, p, dim, d;
670: PetscMPIInt size;
671: const char *names[2] = {"Analytic", "Empirical"};
672: char title[PETSC_MAX_PATH_LEN];
673: PetscOptions options;
674: const char *prefix;
675: MPI_Comm comm;
677: PetscFunctionBegin;
678: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
679: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, &prefix));
680: PetscCall(PetscObjectGetOptions((PetscObject)v, &options));
681: PetscCall(PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg));
682: if (flg) {
683: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
684: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
685: }
686: if (isascii) {
687: PetscCall(PetscViewerPushFormat(viewer, format));
688: } else if (isdraw) {
689: PetscCall(PetscViewerPushFormat(viewer, format));
690: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
691: PetscCall(PetscDrawLGCreate(draw, 2, &lg));
692: PetscCall(PetscDrawLGSetLegend(lg, names));
693: }
695: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
696: PetscCall(VecGetLocalSize(v, &n));
697: PetscCall(VecGetBlockSize(v, &dim));
698: n /= dim;
699: /* TODO Parallel algorithm is harder */
700: if (size == 1) {
701: PetscCall(PetscMalloc1(n, &speed));
702: PetscCall(VecGetArrayRead(v, &a));
703: for (p = 0; p < n; ++p) {
704: PetscReal mag = 0.;
706: for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d]));
707: speed[p] = PetscSqrtReal(mag);
708: }
709: PetscCall(PetscSortReal(n, speed));
710: PetscCall(VecRestoreArrayRead(v, &a));
711: for (p = 0; p < n; ++p) {
712: const PetscReal x = speed[p], Fn = ((PetscReal)p) / n;
713: PetscReal F, vals[2];
715: PetscCall(cdf(&x, NULL, &F));
716: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
717: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
718: if (isdraw) {
719: vals[0] = F;
720: vals[1] = Fn;
721: PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
722: }
723: }
724: if (speed[n - 1] < 6.) {
725: const PetscReal k = (PetscInt)(6. - speed[n - 1]) + 1, dx = (6. - speed[n - 1]) / k;
726: for (p = 0; p <= k; ++p) {
727: const PetscReal x = speed[n - 1] + p * dx, Fn = 1.0;
728: PetscReal F, vals[2];
730: PetscCall(cdf(&x, NULL, &F));
731: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
732: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
733: if (isdraw) {
734: vals[0] = F;
735: vals[1] = Fn;
736: PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
737: }
738: }
739: }
740: PetscCall(PetscFree(speed));
741: }
742: if (isdraw) {
743: PetscCall(PetscDrawLGGetAxis(lg, &axis));
744: PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn));
745: PetscCall(PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)"));
746: PetscCall(PetscDrawLGDraw(lg));
747: PetscCall(PetscDrawLGDestroy(&lg));
748: }
749: if (viewer) {
750: PetscCall(PetscViewerFlush(viewer));
751: PetscCall(PetscViewerPopFormat(viewer));
752: PetscCall(PetscViewerDestroy(&viewer));
753: }
754: *alpha = KSfbar((int)n, (double)Dn);
755: PetscFunctionReturn(PETSC_SUCCESS);
756: #endif
757: }