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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], PetscReal x[])
192: {
193: const PetscReal q = 2 * p[0] - 1.;
194: const PetscInt maxIter = 100;
195: PetscReal ck[100], r = 0.;
197: PetscFunctionBeginHot;
198: /* Transform input to [-1, 1] since the code below computes the inverse error function */
199: for (PetscInt k = 0; k < maxIter; ++k) ck[k] = 0.;
200: ck[0] = 1;
201: r = ck[0] * (PetscSqrtReal(PETSC_PI) / 2.) * q;
202: for (PetscInt k = 1; k < maxIter; ++k) {
203: const PetscReal temp = 2. * k + 1.;
205: for (PetscInt m = 0; m <= k - 1; ++m) {
206: PetscReal denom = (m + 1.) * (2. * m + 1.);
208: ck[k] += (ck[m] * ck[k - 1 - m]) / denom;
209: }
210: r += (ck[k] / temp) * PetscPowReal((PetscSqrtReal(PETSC_PI) / 2.) * q, 2. * k + 1.);
211: }
212: /* Scale erfinv() by \sqrt{\pi/2} */
213: x[0] = PetscSqrtReal(PETSC_PI * 0.5) * r;
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /*@
218: PetscPDFGaussian2D - PDF for the Gaussian distribution in 2D
220: Not Collective
222: Input Parameters:
223: + x - Coordinate in $[-\infty, \infty]^2$
224: - unused - ignored
226: Output Parameter:
227: . p - The probability density at `x`
229: Level: beginner
231: .seealso: `PetscPDFSampleGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()`
232: @*/
233: PetscErrorCode PetscPDFGaussian2D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
234: {
235: p[0] = (1. / PETSC_PI) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1])));
236: return PETSC_SUCCESS;
237: }
239: /*@
240: PetscPDFSampleGaussian2D - Sample uniformly from a Gaussian distribution in 2D
241: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
243: Not Collective
245: Input Parameters:
246: + p - A uniform variable on $[0, 1]^2$
247: - unused - ignored
249: Output Parameter:
250: . x - Coordinate in $[-\infty, \infty]^2 $
252: Level: beginner
254: .seealso: `PetscPDFGaussian2D()`, `PetscPDFMaxwellBoltzmann3D()`
255: @*/
256: PetscErrorCode PetscPDFSampleGaussian2D(const PetscReal p[], const PetscReal unused[], PetscReal x[])
257: {
258: const PetscReal mag = PetscSqrtReal(-2.0 * PetscLogReal(p[0]));
259: x[0] = mag * PetscCosReal(2.0 * PETSC_PI * p[1]);
260: x[1] = mag * PetscSinReal(2.0 * PETSC_PI * p[1]);
261: return PETSC_SUCCESS;
262: }
264: /*@
265: PetscPDFGaussian3D - PDF for the Gaussian distribution in 3D
267: Not Collective
269: Input Parameters:
270: + x - Coordinate in $[-\infty, \infty]^3$
271: - unused - ignored
273: Output Parameter:
274: . p - The probability density at `x`
276: Level: beginner
278: .seealso: `PetscPDFSampleGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
279: @*/
280: PetscErrorCode PetscPDFGaussian3D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
281: {
282: p[0] = (1. / PETSC_PI * PetscSqrtReal(PETSC_PI)) * PetscExpReal(-0.5 * (PetscSqr(x[0]) + PetscSqr(x[1]) + PetscSqr(x[2])));
283: return PETSC_SUCCESS;
284: }
286: /*@
287: PetscPDFSampleGaussian3D - Sample uniformly from a Gaussian distribution in 3D
288: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
290: Not Collective
292: Input Parameters:
293: + p - A uniform variable on $[0, 1]^3$
294: - unused - ignored
296: Output Parameter:
297: . x - Coordinate in $[-\infty, \infty]^3$
299: Level: beginner
301: .seealso: `PetscPDFGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
302: @*/
303: PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal unused[], PetscReal x[])
304: {
305: PetscCall(PetscPDFSampleGaussian1D(p, unused, x));
306: PetscCall(PetscPDFSampleGaussian2D(&p[1], unused, &x[1]));
307: return PETSC_SUCCESS;
308: }
310: /*@
311: PetscPDFConstant1D - PDF for the uniform distribution in 1D
313: Not Collective
315: Input Parameters:
316: + x - Coordinate in $[-1, 1]$
317: - unused - Unused
319: Output Parameter:
320: . p - The probability density at `x`
322: Level: beginner
324: .seealso: `PetscCDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscPDFConstant2D()`, `PetscPDFConstant3D()`
325: @*/
326: PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
327: {
328: p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.;
329: return PETSC_SUCCESS;
330: }
332: /*@
333: PetscCDFConstant1D - CDF for the uniform distribution in 1D
335: Not Collective
337: Input Parameters:
338: + x - Coordinate in $[-1, 1]$
339: - unused - Unused
341: Output Parameter:
342: . p - The cumulative probability at `x`
344: Level: beginner
346: .seealso: `PetscPDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscCDFConstant2D()`, `PetscCDFConstant3D()`
347: @*/
348: PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
349: {
350: p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.));
351: return PETSC_SUCCESS;
352: }
354: /*@
355: PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D
357: Not Collective
359: Input Parameters:
360: + p - A uniform variable on $[0, 1]$
361: - unused - Unused
363: Output Parameter:
364: . x - Coordinate in $[-1, 1]$
366: Level: beginner
368: .seealso: `PetscPDFConstant1D()`, `PetscCDFConstant1D()`, `PetscPDFSampleConstant2D()`, `PetscPDFSampleConstant3D()`
369: @*/
370: PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal unused[], PetscReal x[])
371: {
372: x[0] = 2. * p[0] - 1.;
373: return PETSC_SUCCESS;
374: }
376: /*@
377: PetscPDFConstant2D - PDF for the uniform distribution in 2D
379: Not Collective
381: Input Parameters:
382: + x - Coordinate in $[-1, 1]^2$
383: - unused - Unused
385: Output Parameter:
386: . p - The probability density at `x`
388: Level: beginner
390: .seealso: `PetscCDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscPDFConstant1D()`, `PetscPDFConstant3D()`
391: @*/
392: PetscErrorCode PetscPDFConstant2D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
393: {
394: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. ? 0.25 : 0.;
395: return PETSC_SUCCESS;
396: }
398: /*@
399: PetscCDFConstant2D - CDF for the uniform distribution in 2D
401: Not Collective
403: Input Parameters:
404: + x - Coordinate in $[-1, 1]^2$
405: - unused - Unused
407: Output Parameter:
408: . p - The cumulative probability at `x`
410: Level: beginner
412: .seealso: `PetscPDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscCDFConstant1D()`, `PetscCDFConstant3D()`
413: @*/
414: PetscErrorCode PetscCDFConstant2D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
415: {
416: 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.));
417: return PETSC_SUCCESS;
418: }
420: /*@
421: PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on $[-1, 1]^2$ in 2D
423: Not Collective
425: Input Parameters:
426: + p - Two uniform variables on $[0, 1]$
427: - unused - Unused
429: Output Parameter:
430: . x - Coordinate in $[-1, 1]^2$
432: Level: beginner
434: .seealso: `PetscPDFConstant2D()`, `PetscCDFConstant2D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant3D()`
435: @*/
436: PetscErrorCode PetscPDFSampleConstant2D(const PetscReal p[], const PetscReal unused[], PetscReal x[])
437: {
438: x[0] = 2. * p[0] - 1.;
439: x[1] = 2. * p[1] - 1.;
440: return PETSC_SUCCESS;
441: }
443: /*@
444: PetscPDFConstant3D - PDF for the uniform distribution in 3D
446: Not Collective
448: Input Parameters:
449: + x - Coordinate in $[-1, 1]^3$
450: - unused - Unused
452: Output Parameter:
453: . p - The probability density at `x`
455: Level: beginner
457: .seealso: `PetscCDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
458: @*/
459: PetscErrorCode PetscPDFConstant3D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
460: {
461: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. && x[2] > -1. && x[2] <= 1. ? 0.125 : 0.;
462: return PETSC_SUCCESS;
463: }
465: /*@
466: PetscCDFConstant3D - CDF for the uniform distribution in 3D
468: Not Collective
470: Input Parameters:
471: + x - Coordinate in $[-1, 1]^3$
472: - unused - Unused
474: Output Parameter:
475: . p - The cumulative probability at `x`
477: Level: beginner
479: .seealso: `PetscPDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscCDFConstant1D()`, `PetscCDFConstant2D()`
480: @*/
481: PetscErrorCode PetscCDFConstant3D(const PetscReal x[], const PetscReal unused[], PetscReal p[])
482: {
483: 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.));
484: return PETSC_SUCCESS;
485: }
487: /*@
488: PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on $[-1, 1]^3$ in 3D
490: Not Collective
492: Input Parameters:
493: + p - Three uniform variables on $[0, 1]$
494: - unused - Unused
496: Output Parameter:
497: . x - Coordinate in $[-1, 1]^3$
499: Level: beginner
501: .seealso: `PetscPDFConstant3D()`, `PetscCDFConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
502: @*/
503: PetscErrorCode PetscPDFSampleConstant3D(const PetscReal p[], const PetscReal unused[], PetscReal x[])
504: {
505: x[0] = 2. * p[0] - 1.;
506: x[1] = 2. * p[1] - 1.;
507: x[2] = 2. * p[2] - 1.;
508: return PETSC_SUCCESS;
509: }
511: /*@C
512: PetscProbCreateFromOptions - Return the probability distribution specified by the arguments and options
514: Not Collective
516: Input Parameters:
517: + dim - The dimension of sample points
518: . prefix - The options prefix, or `NULL`
519: - name - The options database name for the probability distribution type
521: Output Parameters:
522: + pdf - The PDF of this type, or `NULL`
523: . cdf - The CDF of this type, or `NULL`
524: - sampler - The PDF sampler of this type, or `NULL`
526: Level: intermediate
528: .seealso: `PetscProbFn`, `PetscPDFMaxwellBoltzmann1D()`, `PetscPDFGaussian1D()`, `PetscPDFConstant1D()`
529: @*/
530: PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFn **pdf, PetscProbFn **cdf, PetscProbFn **sampler)
531: {
532: DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN;
534: PetscFunctionBegin;
535: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");
536: PetscCall(PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum)den, (PetscEnum *)&den, NULL));
537: PetscOptionsEnd();
539: if (pdf) {
540: PetscAssertPointer(pdf, 4);
541: *pdf = NULL;
542: }
543: if (cdf) {
544: PetscAssertPointer(cdf, 5);
545: *cdf = NULL;
546: }
547: if (sampler) {
548: PetscAssertPointer(sampler, 6);
549: *sampler = NULL;
550: }
551: switch (den) {
552: case DTPROB_DENSITY_CONSTANT:
553: switch (dim) {
554: case 1:
555: if (pdf) *pdf = PetscPDFConstant1D;
556: if (cdf) *cdf = PetscCDFConstant1D;
557: if (sampler) *sampler = PetscPDFSampleConstant1D;
558: break;
559: case 2:
560: if (pdf) *pdf = PetscPDFConstant2D;
561: if (cdf) *cdf = PetscCDFConstant2D;
562: if (sampler) *sampler = PetscPDFSampleConstant2D;
563: break;
564: case 3:
565: if (pdf) *pdf = PetscPDFConstant3D;
566: if (cdf) *cdf = PetscCDFConstant3D;
567: if (sampler) *sampler = PetscPDFSampleConstant3D;
568: break;
569: default:
570: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
571: }
572: break;
573: case DTPROB_DENSITY_GAUSSIAN:
574: switch (dim) {
575: case 1:
576: if (pdf) *pdf = PetscPDFGaussian1D;
577: if (cdf) *cdf = PetscCDFGaussian1D;
578: if (sampler) *sampler = PetscPDFSampleGaussian1D;
579: break;
580: case 2:
581: if (pdf) *pdf = PetscPDFGaussian2D;
582: if (sampler) *sampler = PetscPDFSampleGaussian2D;
583: break;
584: case 3:
585: if (pdf) *pdf = PetscPDFGaussian3D;
586: if (sampler) *sampler = PetscPDFSampleGaussian3D;
587: break;
588: default:
589: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
590: }
591: break;
592: case DTPROB_DENSITY_MAXWELL_BOLTZMANN:
593: switch (dim) {
594: case 1:
595: if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D;
596: if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D;
597: break;
598: case 2:
599: if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D;
600: if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D;
601: break;
602: case 3:
603: if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D;
604: if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D;
605: break;
606: default:
607: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
608: }
609: break;
610: default:
611: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]);
612: }
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: #ifdef PETSC_HAVE_KS
617: EXTERN_C_BEGIN
618: #include <KolmogorovSmirnovDist.h>
619: EXTERN_C_END
620: #endif
622: typedef enum {
623: NONE,
624: ASCII,
625: DRAW
626: } OutputType;
628: static PetscErrorCode KSViewerCreate(PetscObject obj, OutputType *outputType, PetscViewer *viewer)
629: {
630: PetscViewerFormat format;
631: PetscOptions options;
632: const char *prefix;
633: PetscBool flg;
634: MPI_Comm comm;
636: PetscFunctionBegin;
637: *outputType = NONE;
638: PetscCall(PetscObjectGetComm(obj, &comm));
639: PetscCall(PetscObjectGetOptionsPrefix(obj, &prefix));
640: PetscCall(PetscObjectGetOptions(obj, &options));
641: PetscCall(PetscOptionsCreateViewer(comm, options, prefix, "-ks_monitor", viewer, &format, &flg));
642: if (flg) {
643: PetscCall(PetscObjectTypeCompare((PetscObject)*viewer, PETSCVIEWERASCII, &flg));
644: if (flg) *outputType = ASCII;
645: PetscCall(PetscObjectTypeCompare((PetscObject)*viewer, PETSCVIEWERDRAW, &flg));
646: if (flg) *outputType = DRAW;
647: PetscCall(PetscViewerPushFormat(*viewer, format));
648: }
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: static PetscErrorCode KSViewerDestroy(PetscViewer *viewer)
653: {
654: PetscFunctionBegin;
655: if (*viewer) {
656: PetscCall(PetscViewerFlush(*viewer));
657: PetscCall(PetscViewerPopFormat(*viewer));
658: PetscCall(PetscViewerDestroy(viewer));
659: }
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: static PetscErrorCode PetscProbComputeKSStatistic_Internal(MPI_Comm comm, PetscInt n, PetscReal val[], PetscReal wgt[], PetscProbFn *cdf, PetscReal *alpha, OutputType outputType, PetscViewer viewer)
664: {
665: #if !defined(PETSC_HAVE_KS)
666: SETERRQ(comm, PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
667: #else
668: PetscDraw draw;
669: PetscDrawLG lg;
670: PetscDrawAxis axis;
671: const char *names[2] = {"Analytic", "Empirical"};
672: char title[PETSC_MAX_PATH_LEN];
673: PetscReal Fn = 0., Dn = PETSC_MIN_REAL;
674: PetscMPIInt size;
676: PetscFunctionBegin;
677: PetscCallMPI(MPI_Comm_size(comm, &size));
678: PetscCheck(size == 1, comm, PETSC_ERR_SUP, "Parallel K-S test not yet supported");
679: if (outputType == DRAW) {
680: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
681: PetscCall(PetscDrawLGCreate(draw, 2, &lg));
682: PetscCall(PetscDrawLGSetLegend(lg, names));
683: }
684: if (wgt) {
685: PetscReal *tmpv, *tmpw;
686: PetscInt *perm;
688: PetscCall(PetscMalloc3(n, &perm, n, &tmpv, n, &tmpw));
689: for (PetscInt i = 0; i < n; ++i) perm[i] = i;
690: PetscCall(PetscSortRealWithPermutation(n, val, perm));
691: for (PetscInt i = 0; i < n; ++i) {
692: tmpv[i] = val[perm[i]];
693: tmpw[i] = wgt[perm[i]];
694: }
695: for (PetscInt i = 0; i < n; ++i) {
696: val[i] = tmpv[i];
697: wgt[i] = tmpw[i];
698: }
699: PetscCall(PetscFree3(perm, tmpv, tmpw));
700: } else PetscCall(PetscSortReal(n, val));
701: // Compute empirical cumulative distribution F_n and discrepancy D_n
702: for (PetscInt p = 0; p < n; ++p) {
703: const PetscReal x = val[p];
704: const PetscReal w = wgt ? wgt[p] : 1. / n;
705: PetscReal F, vals[2];
707: Fn += w;
708: PetscCall(cdf(&x, NULL, &F));
709: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
710: switch (outputType) {
711: case ASCII:
712: PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
713: break;
714: case DRAW:
715: vals[0] = F;
716: vals[1] = Fn;
717: PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
718: break;
719: case NONE:
720: break;
721: }
722: }
723: if (outputType == DRAW) {
724: PetscCall(PetscDrawLGGetAxis(lg, &axis));
725: PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn));
726: PetscCall(PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)"));
727: PetscCall(PetscDrawLGDraw(lg));
728: PetscCall(PetscDrawLGDestroy(&lg));
729: }
730: *alpha = KSfbar((int)n, (double)Dn);
731: if (outputType == ASCII) PetscCall(PetscViewerASCIIPrintf(viewer, "KSfbar(%" PetscInt_FMT ", %.2g): %g\n", n, Dn, *alpha));
732: PetscFunctionReturn(PETSC_SUCCESS);
733: #endif
734: }
736: /*@C
737: PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.
739: Collective
741: Input Parameters:
742: + v - The data vector, blocksize is the sample dimension
743: - cdf - The analytic CDF
745: Output Parameter:
746: . alpha - The KS statistic
748: Level: advanced
750: Notes:
751: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
753: $$
754: D_n = \sup_x \left| F_n(x) - F(x) \right|
755: $$
757: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by
759: $$
760: F_n = # of samples <= x / n
761: $$
763: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
764: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
765: the largest absolute difference between the two distribution functions across all $x$ values.
767: The goodness-of-fit test, or Kolmogorov-Smirnov test, is constructed using the Kolmogorov
768: distribution. It rejects the null hypothesis at level $\alpha$ if
770: $$
771: \sqrt{n} D_{n} > K_{\alpha},
772: $$
774: where $K_\alpha$ is found from
776: $$
777: \operatorname{Pr}(K \leq K_{\alpha}) = 1 - \alpha.
778: $$
780: This means that getting a small alpha says that we have high confidence that the data did not come
781: from the input distribution, so we say that it rejects the null hypothesis.
783: .seealso: `PetscProbComputeKSStatisticWeighted()`, `PetscProbComputeKSStatisticMagnitude()`, `PetscProbFn`
784: @*/
785: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFn *cdf, PetscReal *alpha)
786: {
787: PetscViewer viewer = NULL;
788: OutputType outputType = NONE;
789: const PetscScalar *val;
790: PetscInt n;
792: PetscFunctionBegin;
793: PetscCall(KSViewerCreate((PetscObject)v, &outputType, &viewer));
794: PetscCall(VecGetLocalSize(v, &n));
795: PetscCall(VecGetArrayRead(v, &val));
796: PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "K-S test does not support complex");
797: PetscCall(PetscProbComputeKSStatistic_Internal(PetscObjectComm((PetscObject)v), n, (PetscReal *)val, NULL, cdf, alpha, outputType, viewer));
798: PetscCall(VecRestoreArrayRead(v, &val));
799: PetscCall(KSViewerDestroy(&viewer));
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: /*@C
804: PetscProbComputeKSStatisticWeighted - Compute the Kolmogorov-Smirnov statistic for the weighted empirical distribution for an input vector, compared to an analytic CDF.
806: Collective
808: Input Parameters:
809: + v - The data vector, blocksize is the sample dimension
810: . w - The vector of weights for each sample, instead of the default 1/n
811: - cdf - The analytic CDF
813: Output Parameter:
814: . alpha - The KS statistic
816: Level: advanced
818: Notes:
819: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
821: $$
822: D_n = \sup_x \left| F_n(x) - F(x) \right|
823: $$
825: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by
827: $$
828: F_n = # of samples <= x / n
829: $$
831: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
832: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
833: the largest absolute difference between the two distribution functions across all $x$ values.
835: The goodness-of-fit test, or Kolmogorov-Smirnov test, is constructed using the Kolmogorov
836: distribution. It rejects the null hypothesis at level $\alpha$ if
838: $$
839: \sqrt{n} D_{n} > K_{\alpha},
840: $$
842: where $K_\alpha$ is found from
844: $$
845: \operatorname{Pr}(K \leq K_{\alpha}) = 1 - \alpha.
846: $$
848: This means that getting a small alpha says that we have high confidence that the data did not come
849: from the input distribution, so we say that it rejects the null hypothesis.
851: .seealso: `PetscProbComputeKSStatistic()`, `PetscProbComputeKSStatisticMagnitude()`, `PetscProbFn`
852: @*/
853: PetscErrorCode PetscProbComputeKSStatisticWeighted(Vec v, Vec w, PetscProbFn *cdf, PetscReal *alpha)
854: {
855: PetscViewer viewer = NULL;
856: OutputType outputType = NONE;
857: const PetscScalar *val, *wgt;
858: PetscInt n;
860: PetscFunctionBegin;
861: PetscCall(KSViewerCreate((PetscObject)v, &outputType, &viewer));
862: PetscCall(VecGetLocalSize(v, &n));
863: PetscCall(VecGetArrayRead(v, &val));
864: PetscCall(VecGetArrayRead(w, &wgt));
865: PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "K-S test does not support complex");
866: PetscCall(PetscProbComputeKSStatistic_Internal(PetscObjectComm((PetscObject)v), n, (PetscReal *)val, (PetscReal *)wgt, cdf, alpha, outputType, viewer));
867: PetscCall(VecRestoreArrayRead(v, &val));
868: PetscCall(VecRestoreArrayRead(w, &wgt));
869: PetscCall(KSViewerDestroy(&viewer));
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@C
874: PetscProbComputeKSStatisticMagnitude - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for the magnitude over each block of an input vector, compared to an analytic CDF.
876: Collective
878: Input Parameters:
879: + v - The data vector, blocksize is the sample dimension
880: - cdf - The analytic CDF
882: Output Parameter:
883: . alpha - The KS statistic
885: Level: advanced
887: Notes:
888: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
890: $$
891: D_n = \sup_x \left| F_n(x) - F(x) \right|
892: $$
894: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by
896: $$
897: F_n = # of samples <= x / n
898: $$
900: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
901: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
902: the largest absolute difference between the two distribution functions across all $x$ values.
904: The goodness-of-fit test, or Kolmogorov-Smirnov test, is constructed using the Kolmogorov
905: distribution. It rejects the null hypothesis at level $\alpha$ if
907: $$
908: \sqrt{n} D_{n} > K_{\alpha},
909: $$
911: where $K_\alpha$ is found from
913: $$
914: \operatorname{Pr}(K \leq K_{\alpha}) = 1 - \alpha.
915: $$
917: This means that getting a small alpha says that we have high confidence that the data did not come
918: from the input distribution, so we say that it rejects the null hypothesis.
920: .seealso: `PetscProbComputeKSStatistic()`, `PetscProbComputeKSStatisticWeighted()`, `PetscProbFn`
921: @*/
922: PetscErrorCode PetscProbComputeKSStatisticMagnitude(Vec v, PetscProbFn *cdf, PetscReal *alpha)
923: {
924: PetscViewer viewer = NULL;
925: OutputType outputType = NONE;
926: const PetscScalar *a;
927: PetscReal *speed;
928: PetscInt n, dim;
930: PetscFunctionBegin;
931: PetscCall(KSViewerCreate((PetscObject)v, &outputType, &viewer));
932: // Convert to a scalar value
933: PetscCall(VecGetLocalSize(v, &n));
934: PetscCall(VecGetBlockSize(v, &dim));
935: n /= dim;
936: PetscCall(PetscMalloc1(n, &speed));
937: PetscCall(VecGetArrayRead(v, &a));
938: for (PetscInt p = 0; p < n; ++p) {
939: PetscReal mag = 0.;
941: for (PetscInt d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d]));
942: speed[p] = PetscSqrtReal(mag);
943: }
944: PetscCall(VecRestoreArrayRead(v, &a));
945: // Compute statistic
946: PetscCall(PetscProbComputeKSStatistic_Internal(PetscObjectComm((PetscObject)v), n, speed, NULL, cdf, alpha, outputType, viewer));
947: PetscCall(PetscFree(speed));
948: PetscCall(KSViewerDestroy(&viewer));
949: PetscFunctionReturn(PETSC_SUCCESS);
950: }