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: PetscPDFConstant1D - PDF for the uniform distribution in 1D
259: Not collective
261: Input Parameter:
262: . x - Coordinate in [-1, 1]
264: Output Parameter:
265: . p - The probability density at x
267: Level: beginner
269: .seealso: PetscCDFConstant1D()
270: @*/
271: PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
272: {
273: p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.;
274: return 0;
275: }
277: PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
278: {
279: p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5*(x[0] + 1.));
280: return 0;
281: }
283: /*@
284: PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D
286: Not collective
288: Input Parameter:
289: . p - A uniform variable on [0, 1]
291: Output Parameter:
292: . x - Coordinate in [-1, 1]
294: Level: beginner
296: .seealso: PetscPDFConstant1D(), PetscCDFConstant1D()
297: @*/
298: PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
299: {
300: x[0] = 2.*p[1] - 1.;
301: return 0;
302: }
304: /*@C
305: PetscProbCreateFromOptions - Return the probability distribution specified by the argumetns and options
307: Not collective
309: Input Parameters:
310: + dim - The dimension of sample points
311: . prefix - The options prefix, or NULL
312: - name - The option name for the probility distribution type
314: Output Parameters:
315: + pdf - The PDF of this type
316: . cdf - The CDF of this type
317: - sampler - The PDF sampler of this type
319: Level: intermediate
321: .seealso: PetscPDFMaxwellBoltzmann1D(), PetscPDFGaussian1D(), PetscPDFConstant1D()
322: @*/
323: PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler)
324: {
325: DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN;
326: PetscErrorCode ierr;
328: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");
329: PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum) den, (PetscEnum *) &den, NULL);
330: PetscOptionsEnd();
335: switch (den) {
336: case DTPROB_DENSITY_CONSTANT:
337: switch (dim) {
338: case 1:
339: if (pdf) *pdf = PetscPDFConstant1D;
340: if (cdf) *cdf = PetscCDFConstant1D;
341: if (sampler) *sampler = PetscPDFSampleConstant1D;
342: break;
343: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
344: }
345: break;
346: case DTPROB_DENSITY_GAUSSIAN:
347: switch (dim) {
348: case 1:
349: if (pdf) *pdf = PetscPDFGaussian1D;
350: if (cdf) *cdf = PetscCDFGaussian1D;
351: if (sampler) *sampler = PetscPDFSampleGaussian1D;
352: break;
353: case 2:
354: if (pdf) *pdf = PetscPDFGaussian2D;
355: if (sampler) *sampler = PetscPDFSampleGaussian2D;
356: break;
357: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
358: }
359: break;
360: case DTPROB_DENSITY_MAXWELL_BOLTZMANN:
361: switch (dim) {
362: case 1:
363: if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D;
364: if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D;
365: break;
366: case 2:
367: if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D;
368: if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D;
369: break;
370: case 3:
371: if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D;
372: if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D;
373: break;
374: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
375: }
376: break;
377: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]);
378: }
379: return 0;
380: }
382: #ifdef PETSC_HAVE_KS
383: EXTERN_C_BEGIN
384: #include <KolmogorovSmirnovDist.h>
385: EXTERN_C_END
386: #endif
388: /*@C
389: PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.
391: Collective on v
393: Input Parameters:
394: + v - The data vector, blocksize is the sample dimension
395: - cdf - The analytic CDF
397: Output Parameter:
398: . alpha - The KS statisic
400: Level: advanced
402: Note: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
403: $
404: $ D_n = \sup_x \left| F_n(x) - F(x) \right|
405: $
406: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution
407: function $F_n(x)$ is discrete, and given by
408: $
409: $ F_n = # of samples <= x / n
410: $
411: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
412: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
413: the largest absolute difference between the two distribution functions across all $x$ values.
415: .seealso: PetscProbFunc
416: @*/
417: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha)
418: {
419: #if !defined(PETSC_HAVE_KS)
420: SETERRQ(PetscObjectComm((PetscObject) v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
421: #else
422: PetscViewer viewer = NULL;
423: PetscViewerFormat format;
424: PetscDraw draw;
425: PetscDrawLG lg;
426: PetscDrawAxis axis;
427: const PetscScalar *a;
428: PetscReal *speed, Dn = PETSC_MIN_REAL;
429: PetscBool isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg;
430: PetscInt n, p, dim, d;
431: PetscMPIInt size;
432: const char *names[2] = {"Analytic", "Empirical"};
433: char title[PETSC_MAX_PATH_LEN];
434: PetscOptions options;
435: const char *prefix;
436: MPI_Comm comm;
438: PetscObjectGetComm((PetscObject) v, &comm);
439: PetscObjectGetOptionsPrefix((PetscObject) v, &prefix);
440: PetscObjectGetOptions((PetscObject) v, &options);
441: PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg);
442: if (flg) {
443: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
444: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);
445: }
446: if (isascii) {
447: PetscViewerPushFormat(viewer, format);
448: } else if (isdraw) {
449: PetscViewerPushFormat(viewer, format);
450: PetscViewerDrawGetDraw(viewer, 0, &draw);
451: PetscDrawLGCreate(draw, 2, &lg);
452: PetscDrawLGSetLegend(lg, names);
453: }
455: MPI_Comm_size(PetscObjectComm((PetscObject) v), &size);
456: VecGetLocalSize(v, &n);
457: VecGetBlockSize(v, &dim);
458: n /= dim;
459: /* TODO Parallel algorithm is harder */
460: if (size == 1) {
461: PetscMalloc1(n, &speed);
462: VecGetArrayRead(v, &a);
463: for (p = 0; p < n; ++p) {
464: PetscReal mag = 0.;
466: for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p*dim+d]));
467: speed[p] = PetscSqrtReal(mag);
468: }
469: PetscSortReal(n, speed);
470: VecRestoreArrayRead(v, &a);
471: for (p = 0; p < n; ++p) {
472: const PetscReal x = speed[p], Fn = ((PetscReal) p) / n;
473: PetscReal F, vals[2];
475: cdf(&x, NULL, &F);
476: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
477: if (isascii) PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);
478: if (isdraw) {vals[0] = F; vals[1] = Fn; PetscDrawLGAddCommonPoint(lg, x, vals);}
479: }
480: if (speed[n-1] < 6.) {
481: const PetscReal k = (PetscInt) (6. - speed[n-1]) + 1, dx = (6. - speed[n-1])/k;
482: for (p = 0; p <= k; ++p) {
483: const PetscReal x = speed[n-1] + p*dx, Fn = 1.0;
484: PetscReal F, vals[2];
486: cdf(&x, NULL, &F);
487: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
488: if (isascii) PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn);
489: if (isdraw) {vals[0] = F; vals[1] = Fn; PetscDrawLGAddCommonPoint(lg, x, vals);}
490: }
491: }
492: PetscFree(speed);
493: }
494: if (isdraw) {
495: PetscDrawLGGetAxis(lg, &axis);
496: PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn);
497: PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)");
498: PetscDrawLGDraw(lg);
499: PetscDrawLGDestroy(&lg);
500: }
501: if (viewer) {
502: PetscViewerFlush(viewer);
503: PetscViewerPopFormat(viewer);
504: PetscViewerDestroy(&viewer);
505: }
506: *alpha = KSfbar((int) n, (double) Dn);
507: return 0;
508: #endif
509: }