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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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;
188: PetscFunctionBeginHot;
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: PetscFunctionReturn(PETSC_SUCCESS);
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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: Reference:
293: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
295: .seealso: `PetscPDFGaussian3D()`, `PetscPDFMaxwellBoltzmann3D()`
296: @*/
297: PetscErrorCode PetscPDFSampleGaussian3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
298: {
299: PetscCall(PetscPDFSampleGaussian1D(p, dummy, x));
300: PetscCall(PetscPDFSampleGaussian2D(&p[1], dummy, &x[1]));
301: return PETSC_SUCCESS;
302: }
304: /*@
305: PetscPDFConstant1D - PDF for the uniform distribution in 1D
307: Not Collective
309: Input Parameter:
310: . x - Coordinate in [-1, 1]
312: Output Parameter:
313: . p - The probability density at x
315: Level: beginner
317: .seealso: `PetscCDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscPDFConstant2D()`, `PetscPDFConstant3D()`
318: @*/
319: PetscErrorCode PetscPDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
320: {
321: p[0] = x[0] > -1. && x[0] <= 1. ? 0.5 : 0.;
322: return PETSC_SUCCESS;
323: }
325: /*@
326: PetscCDFConstant1D - CDF for the uniform distribution in 1D
328: Not Collective
330: Input Parameter:
331: . x - Coordinate in [-1, 1]
333: Output Parameter:
334: . p - The cumulative probability at x
336: Level: beginner
338: .seealso: `PetscPDFConstant1D()`, `PetscPDFSampleConstant1D()`, `PetscCDFConstant2D()`, `PetscCDFConstant3D()`
339: @*/
340: PetscErrorCode PetscCDFConstant1D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
341: {
342: p[0] = x[0] <= -1. ? 0. : (x[0] > 1. ? 1. : 0.5 * (x[0] + 1.));
343: return PETSC_SUCCESS;
344: }
346: /*@
347: PetscPDFSampleConstant1D - Sample uniformly from a uniform distribution on [-1, 1] in 1D
349: Not Collective
351: Input Parameter:
352: . p - A uniform variable on [0, 1]
354: Output Parameter:
355: . x - Coordinate in [-1, 1]
357: Level: beginner
359: .seealso: `PetscPDFConstant1D()`, `PetscCDFConstant1D()`, `PetscPDFSampleConstant2D()`, `PetscPDFSampleConstant3D()`
360: @*/
361: PetscErrorCode PetscPDFSampleConstant1D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
362: {
363: x[0] = 2. * p[0] - 1.;
364: return PETSC_SUCCESS;
365: }
367: /*@
368: PetscPDFConstant2D - PDF for the uniform distribution in 2D
370: Not Collective
372: Input Parameter:
373: . x - Coordinate in [-1, 1] x [-1, 1]
375: Output Parameter:
376: . p - The probability density at x
378: Level: beginner
380: .seealso: `PetscCDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscPDFConstant1D()`, `PetscPDFConstant3D()`
381: @*/
382: PetscErrorCode PetscPDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
383: {
384: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. ? 0.25 : 0.;
385: return PETSC_SUCCESS;
386: }
388: /*@
389: PetscCDFConstant2D - CDF for the uniform distribution in 2D
391: Not Collective
393: Input Parameter:
394: . x - Coordinate in [-1, 1] x [-1, 1]
396: Output Parameter:
397: . p - The cumulative probability at x
399: Level: beginner
401: .seealso: `PetscPDFConstant2D()`, `PetscPDFSampleConstant2D()`, `PetscCDFConstant1D()`, `PetscCDFConstant3D()`
402: @*/
403: PetscErrorCode PetscCDFConstant2D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
404: {
405: 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.));
406: return PETSC_SUCCESS;
407: }
409: /*@
410: PetscPDFSampleConstant2D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 2D
412: Not Collective
414: Input Parameter:
415: . p - Two uniform variables on [0, 1]
417: Output Parameter:
418: . x - Coordinate in [-1, 1] x [-1, 1]
420: Level: beginner
422: .seealso: `PetscPDFConstant2D()`, `PetscCDFConstant2D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant3D()`
423: @*/
424: PetscErrorCode PetscPDFSampleConstant2D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
425: {
426: x[0] = 2. * p[0] - 1.;
427: x[1] = 2. * p[1] - 1.;
428: return PETSC_SUCCESS;
429: }
431: /*@
432: PetscPDFConstant3D - PDF for the uniform distribution in 3D
434: Not Collective
436: Input Parameter:
437: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
439: Output Parameter:
440: . p - The probability density at x
442: Level: beginner
444: .seealso: `PetscCDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
445: @*/
446: PetscErrorCode PetscPDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
447: {
448: p[0] = x[0] > -1. && x[0] <= 1. && x[1] > -1. && x[1] <= 1. && x[2] > -1. && x[2] <= 1. ? 0.125 : 0.;
449: return PETSC_SUCCESS;
450: }
452: /*@
453: PetscCDFConstant3D - CDF for the uniform distribution in 3D
455: Not Collective
457: Input Parameter:
458: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
460: Output Parameter:
461: . p - The cumulative probability at x
463: Level: beginner
465: .seealso: `PetscPDFConstant3D()`, `PetscPDFSampleConstant3D()`, `PetscCDFConstant1D()`, `PetscCDFConstant2D()`
466: @*/
467: PetscErrorCode PetscCDFConstant3D(const PetscReal x[], const PetscReal dummy[], PetscReal p[])
468: {
469: 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.));
470: return PETSC_SUCCESS;
471: }
473: /*@
474: PetscPDFSampleConstant3D - Sample uniformly from a uniform distribution on [-1, 1] x [-1, 1] in 3D
476: Not Collective
478: Input Parameter:
479: . p - Three uniform variables on [0, 1]
481: Output Parameter:
482: . x - Coordinate in [-1, 1] x [-1, 1] x [-1, 1]
484: Level: beginner
486: .seealso: `PetscPDFConstant3D()`, `PetscCDFConstant3D()`, `PetscPDFSampleConstant1D()`, `PetscPDFSampleConstant2D()`
487: @*/
488: PetscErrorCode PetscPDFSampleConstant3D(const PetscReal p[], const PetscReal dummy[], PetscReal x[])
489: {
490: x[0] = 2. * p[0] - 1.;
491: x[1] = 2. * p[1] - 1.;
492: x[2] = 2. * p[2] - 1.;
493: return PETSC_SUCCESS;
494: }
496: /*@C
497: PetscProbCreateFromOptions - Return the probability distribution specified by the arguments and options
499: Not Collective
501: Input Parameters:
502: + dim - The dimension of sample points
503: . prefix - The options prefix, or NULL
504: - name - The option name for the probility distribution type
506: Output Parameters:
507: + pdf - The PDF of this type
508: . cdf - The CDF of this type
509: - sampler - The PDF sampler of this type
511: Level: intermediate
513: .seealso: `PetscProbFunc`, `PetscPDFMaxwellBoltzmann1D()`, `PetscPDFGaussian1D()`, `PetscPDFConstant1D()`
514: @*/
515: PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFunc *pdf, PetscProbFunc *cdf, PetscProbFunc *sampler)
516: {
517: DTProbDensityType den = DTPROB_DENSITY_GAUSSIAN;
519: PetscFunctionBegin;
520: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "PetscProb Options", "DT");
521: PetscCall(PetscOptionsEnum(name, "Method to compute PDF <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum)den, (PetscEnum *)&den, NULL));
522: PetscOptionsEnd();
524: if (pdf) {
526: *pdf = NULL;
527: }
528: if (cdf) {
530: *cdf = NULL;
531: }
532: if (sampler) {
534: *sampler = NULL;
535: }
536: switch (den) {
537: case DTPROB_DENSITY_CONSTANT:
538: switch (dim) {
539: case 1:
540: if (pdf) *pdf = PetscPDFConstant1D;
541: if (cdf) *cdf = PetscCDFConstant1D;
542: if (sampler) *sampler = PetscPDFSampleConstant1D;
543: break;
544: case 2:
545: if (pdf) *pdf = PetscPDFConstant2D;
546: if (cdf) *cdf = PetscCDFConstant2D;
547: if (sampler) *sampler = PetscPDFSampleConstant2D;
548: break;
549: case 3:
550: if (pdf) *pdf = PetscPDFConstant3D;
551: if (cdf) *cdf = PetscCDFConstant3D;
552: if (sampler) *sampler = PetscPDFSampleConstant3D;
553: break;
554: default:
555: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
556: }
557: break;
558: case DTPROB_DENSITY_GAUSSIAN:
559: switch (dim) {
560: case 1:
561: if (pdf) *pdf = PetscPDFGaussian1D;
562: if (cdf) *cdf = PetscCDFGaussian1D;
563: if (sampler) *sampler = PetscPDFSampleGaussian1D;
564: break;
565: case 2:
566: if (pdf) *pdf = PetscPDFGaussian2D;
567: if (sampler) *sampler = PetscPDFSampleGaussian2D;
568: break;
569: case 3:
570: if (pdf) *pdf = PetscPDFGaussian3D;
571: if (sampler) *sampler = PetscPDFSampleGaussian3D;
572: break;
573: default:
574: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
575: }
576: break;
577: case DTPROB_DENSITY_MAXWELL_BOLTZMANN:
578: switch (dim) {
579: case 1:
580: if (pdf) *pdf = PetscPDFMaxwellBoltzmann1D;
581: if (cdf) *cdf = PetscCDFMaxwellBoltzmann1D;
582: break;
583: case 2:
584: if (pdf) *pdf = PetscPDFMaxwellBoltzmann2D;
585: if (cdf) *cdf = PetscCDFMaxwellBoltzmann2D;
586: break;
587: case 3:
588: if (pdf) *pdf = PetscPDFMaxwellBoltzmann3D;
589: if (cdf) *cdf = PetscCDFMaxwellBoltzmann3D;
590: break;
591: default:
592: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %" PetscInt_FMT " not supported for density type %s", dim, DTProbDensityTypes[den]);
593: }
594: break;
595: default:
596: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Density type %s is not supported", DTProbDensityTypes[PetscMax(0, PetscMin(den, DTPROB_NUM_DENSITY))]);
597: }
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: #ifdef PETSC_HAVE_KS
602: EXTERN_C_BEGIN
603: #include <KolmogorovSmirnovDist.h>
604: EXTERN_C_END
605: #endif
607: /*@C
608: PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.
610: Collective
612: Input Parameters:
613: + v - The data vector, blocksize is the sample dimension
614: - cdf - The analytic CDF
616: Output Parameter:
617: . alpha - The KS statisic
619: Level: advanced
621: Notes:
622: The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is
623: .vb
624: D_n = \sup_x \left| F_n(x) - F(x) \right|
626: where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by
628: F_n = # of samples <= x / n
629: .ve
631: The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
632: cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
633: the largest absolute difference between the two distribution functions across all $x$ values.
635: .seealso: `PetscProbFunc`
636: @*/
637: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFunc cdf, PetscReal *alpha)
638: {
639: #if !defined(PETSC_HAVE_KS)
640: SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
641: #else
642: PetscViewer viewer = NULL;
643: PetscViewerFormat format;
644: PetscDraw draw;
645: PetscDrawLG lg;
646: PetscDrawAxis axis;
647: const PetscScalar *a;
648: PetscReal *speed, Dn = PETSC_MIN_REAL;
649: PetscBool isascii = PETSC_FALSE, isdraw = PETSC_FALSE, flg;
650: PetscInt n, p, dim, d;
651: PetscMPIInt size;
652: const char *names[2] = {"Analytic", "Empirical"};
653: char title[PETSC_MAX_PATH_LEN];
654: PetscOptions options;
655: const char *prefix;
656: MPI_Comm comm;
658: PetscFunctionBegin;
659: PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
660: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, &prefix));
661: PetscCall(PetscObjectGetOptions((PetscObject)v, &options));
662: PetscCall(PetscOptionsGetViewer(comm, options, prefix, "-ks_monitor", &viewer, &format, &flg));
663: if (flg) {
664: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
665: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
666: }
667: if (isascii) {
668: PetscCall(PetscViewerPushFormat(viewer, format));
669: } else if (isdraw) {
670: PetscCall(PetscViewerPushFormat(viewer, format));
671: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
672: PetscCall(PetscDrawLGCreate(draw, 2, &lg));
673: PetscCall(PetscDrawLGSetLegend(lg, names));
674: }
676: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
677: PetscCall(VecGetLocalSize(v, &n));
678: PetscCall(VecGetBlockSize(v, &dim));
679: n /= dim;
680: /* TODO Parallel algorithm is harder */
681: if (size == 1) {
682: PetscCall(PetscMalloc1(n, &speed));
683: PetscCall(VecGetArrayRead(v, &a));
684: for (p = 0; p < n; ++p) {
685: PetscReal mag = 0.;
687: for (d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d]));
688: speed[p] = PetscSqrtReal(mag);
689: }
690: PetscCall(PetscSortReal(n, speed));
691: PetscCall(VecRestoreArrayRead(v, &a));
692: for (p = 0; p < n; ++p) {
693: const PetscReal x = speed[p], Fn = ((PetscReal)p) / n;
694: PetscReal F, vals[2];
696: PetscCall(cdf(&x, NULL, &F));
697: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
698: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
699: if (isdraw) {
700: vals[0] = F;
701: vals[1] = Fn;
702: PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
703: }
704: }
705: if (speed[n - 1] < 6.) {
706: const PetscReal k = (PetscInt)(6. - speed[n - 1]) + 1, dx = (6. - speed[n - 1]) / k;
707: for (p = 0; p <= k; ++p) {
708: const PetscReal x = speed[n - 1] + p * dx, Fn = 1.0;
709: PetscReal F, vals[2];
711: PetscCall(cdf(&x, NULL, &F));
712: Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
713: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
714: if (isdraw) {
715: vals[0] = F;
716: vals[1] = Fn;
717: PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
718: }
719: }
720: }
721: PetscCall(PetscFree(speed));
722: }
723: if (isdraw) {
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: if (viewer) {
731: PetscCall(PetscViewerFlush(viewer));
732: PetscCall(PetscViewerPopFormat(viewer));
733: PetscCall(PetscViewerDestroy(&viewer));
734: }
735: *alpha = KSfbar((int)n, (double)Dn);
736: PetscFunctionReturn(PETSC_SUCCESS);
737: #endif
738: }