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.;
    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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], PetscReal x[])
    305: {
    306:   PetscCall(PetscPDFSampleGaussian1D(p, unused, x));
    307:   PetscCall(PetscPDFSampleGaussian2D(&p[1], unused, &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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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: - unused - 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 unused[], 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 options database name for the probability distribution type

    522:   Output Parameters:
    523: + pdf     - The PDF of this type, or `NULL`
    524: . cdf     - The CDF of this type, or `NULL`
    525: - sampler - The PDF sampler of this type, or `NULL`

    527:   Level: intermediate

    529: .seealso: `PetscProbFn`, `PetscPDFMaxwellBoltzmann1D()`, `PetscPDFGaussian1D()`, `PetscPDFConstant1D()`
    530: @*/
    531: PetscErrorCode PetscProbCreateFromOptions(PetscInt dim, const char prefix[], const char name[], PetscProbFn **pdf, PetscProbFn **cdf, PetscProbFn **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: typedef enum {
    624:   NONE,
    625:   ASCII,
    626:   DRAW
    627: } OutputType;

    629: static PetscErrorCode KSViewerCreate(PetscObject obj, OutputType *outputType, PetscViewer *viewer)
    630: {
    631:   PetscViewerFormat format;
    632:   PetscOptions      options;
    633:   const char       *prefix;
    634:   PetscBool         flg;
    635:   MPI_Comm          comm;

    637:   PetscFunctionBegin;
    638:   *outputType = NONE;
    639:   PetscCall(PetscObjectGetComm(obj, &comm));
    640:   PetscCall(PetscObjectGetOptionsPrefix(obj, &prefix));
    641:   PetscCall(PetscObjectGetOptions(obj, &options));
    642:   PetscCall(PetscOptionsCreateViewer(comm, options, prefix, "-ks_monitor", viewer, &format, &flg));
    643:   if (flg) {
    644:     PetscCall(PetscObjectTypeCompare((PetscObject)*viewer, PETSCVIEWERASCII, &flg));
    645:     if (flg) *outputType = ASCII;
    646:     PetscCall(PetscObjectTypeCompare((PetscObject)*viewer, PETSCVIEWERDRAW, &flg));
    647:     if (flg) *outputType = DRAW;
    648:     PetscCall(PetscViewerPushFormat(*viewer, format));
    649:   }
    650:   PetscFunctionReturn(PETSC_SUCCESS);
    651: }

    653: static PetscErrorCode KSViewerDestroy(PetscViewer *viewer)
    654: {
    655:   PetscFunctionBegin;
    656:   if (*viewer) {
    657:     PetscCall(PetscViewerFlush(*viewer));
    658:     PetscCall(PetscViewerPopFormat(*viewer));
    659:     PetscCall(PetscViewerDestroy(viewer));
    660:   }
    661:   PetscFunctionReturn(PETSC_SUCCESS);
    662: }

    664: static PetscErrorCode PetscProbComputeKSStatistic_Internal(MPI_Comm comm, PetscInt n, PetscReal val[], PetscReal wgt[], PetscProbFn *cdf, PetscReal *alpha, OutputType outputType, PetscViewer viewer)
    665: {
    666: #if !defined(PETSC_HAVE_KS)
    667:   SETERRQ(comm, PETSC_ERR_SUP, "No support for Kolmogorov-Smirnov test.\nReconfigure using --download-ks");
    668: #else
    669:   PetscDraw     draw;
    670:   PetscDrawLG   lg;
    671:   PetscDrawAxis axis;
    672:   const char   *names[2] = {"Analytic", "Empirical"};
    673:   char          title[PETSC_MAX_PATH_LEN];
    674:   PetscReal     Fn = 0., Dn = PETSC_MIN_REAL;
    675:   PetscMPIInt   size;

    677:   PetscFunctionBegin;
    678:   PetscCallMPI(MPI_Comm_size(comm, &size));
    679:   PetscCheck(size == 1, comm, PETSC_ERR_SUP, "Parallel K-S test not yet supported");
    680:   if (outputType == DRAW) {
    681:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    682:     PetscCall(PetscDrawLGCreate(draw, 2, &lg));
    683:     PetscCall(PetscDrawLGSetLegend(lg, names));
    684:   }
    685:   if (wgt) {
    686:     PetscReal *tmpv, *tmpw;
    687:     PetscInt  *perm;

    689:     PetscCall(PetscMalloc3(n, &perm, n, &tmpv, n, &tmpw));
    690:     for (PetscInt i = 0; i < n; ++i) perm[i] = i;
    691:     PetscCall(PetscSortRealWithPermutation(n, val, perm));
    692:     for (PetscInt i = 0; i < n; ++i) {
    693:       tmpv[i] = val[perm[i]];
    694:       tmpw[i] = wgt[perm[i]];
    695:     }
    696:     for (PetscInt i = 0; i < n; ++i) {
    697:       val[i] = tmpv[i];
    698:       wgt[i] = tmpw[i];
    699:     }
    700:     PetscCall(PetscFree3(perm, tmpv, tmpw));
    701:   } else PetscCall(PetscSortReal(n, val));
    702:   // Compute empirical cumulative distribution F_n and discrepancy D_n
    703:   for (PetscInt p = 0; p < n; ++p) {
    704:     const PetscReal x = val[p];
    705:     const PetscReal w = wgt ? wgt[p] : 1. / n;
    706:     PetscReal       F, vals[2];

    708:     Fn += w;
    709:     PetscCall(cdf(&x, NULL, &F));
    710:     Dn = PetscMax(PetscAbsReal(Fn - F), Dn);
    711:     switch (outputType) {
    712:     case ASCII:
    713:       PetscCall(PetscViewerASCIIPrintf(viewer, "x: %g F: %g Fn: %g Dn: %.2g\n", x, F, Fn, Dn));
    714:       break;
    715:     case DRAW:
    716:       vals[0] = F;
    717:       vals[1] = Fn;
    718:       PetscCall(PetscDrawLGAddCommonPoint(lg, x, vals));
    719:       break;
    720:     case NONE:
    721:       break;
    722:     }
    723:   }
    724:   if (outputType == DRAW) {
    725:     PetscCall(PetscDrawLGGetAxis(lg, &axis));
    726:     PetscCall(PetscSNPrintf(title, PETSC_MAX_PATH_LEN, "Kolmogorov-Smirnov Test (Dn %.2g)", Dn));
    727:     PetscCall(PetscDrawAxisSetLabels(axis, title, "x", "CDF(x)"));
    728:     PetscCall(PetscDrawLGDraw(lg));
    729:     PetscCall(PetscDrawLGDestroy(&lg));
    730:   }
    731:   *alpha = KSfbar((int)n, (double)Dn);
    732:   if (outputType == ASCII) PetscCall(PetscViewerASCIIPrintf(viewer, "KSfbar(%" PetscInt_FMT ", %.2g): %g\n", n, Dn, *alpha));
    733:   PetscFunctionReturn(PETSC_SUCCESS);
    734: #endif
    735: }

    737: /*@C
    738:   PetscProbComputeKSStatistic - Compute the Kolmogorov-Smirnov statistic for the empirical distribution for an input vector, compared to an analytic CDF.

    740:   Collective

    742:   Input Parameters:
    743: + v   - The data vector, blocksize is the sample dimension
    744: - cdf - The analytic CDF

    746:   Output Parameter:
    747: . alpha - The KS statistic

    749:   Level: advanced

    751:   Notes:
    752:   The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is

    754:   $$
    755:   D_n = \sup_x \left| F_n(x) - F(x) \right|
    756:   $$

    758:   where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by

    760:   $$
    761:   F_n =  # of samples <= x / n
    762:   $$

    764:   The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
    765:   cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
    766:   the largest absolute difference between the two distribution functions across all $x$ values.

    768:   The goodness-of-fit test, or Kolmogorov-Smirnov test, is constructed using the Kolmogorov
    769:   distribution. It rejects the null hypothesis at level $\alpha$ if

    771:   $$
    772:   \sqrt{n} D_{n} > K_{\alpha},
    773:   $$

    775:   where $K_\alpha$ is found from

    777:   $$
    778:   \operatorname{Pr}(K \leq K_{\alpha}) = 1 - \alpha.
    779:   $$

    781:   This means that getting a small alpha says that we have high confidence that the data did not come
    782:   from the input distribution, so we say that it rejects the null hypothesis.

    784: .seealso: `PetscProbComputeKSStatisticWeighted()`, `PetscProbComputeKSStatisticMagnitude()`, `PetscProbFn`
    785: @*/
    786: PetscErrorCode PetscProbComputeKSStatistic(Vec v, PetscProbFn *cdf, PetscReal *alpha)
    787: {
    788:   PetscViewer        viewer     = NULL;
    789:   OutputType         outputType = NONE;
    790:   const PetscScalar *val;
    791:   PetscInt           n;

    793:   PetscFunctionBegin;
    794:   PetscCall(KSViewerCreate((PetscObject)v, &outputType, &viewer));
    795:   PetscCall(VecGetLocalSize(v, &n));
    796:   PetscCall(VecGetArrayRead(v, &val));
    797:   PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "K-S test does not support complex");
    798:   PetscCall(PetscProbComputeKSStatistic_Internal(PetscObjectComm((PetscObject)v), n, (PetscReal *)val, NULL, cdf, alpha, outputType, viewer));
    799:   PetscCall(VecRestoreArrayRead(v, &val));
    800:   PetscCall(KSViewerDestroy(&viewer));
    801:   PetscFunctionReturn(PETSC_SUCCESS);
    802: }

    804: /*@C
    805:   PetscProbComputeKSStatisticWeighted - Compute the Kolmogorov-Smirnov statistic for the weighted empirical distribution for an input vector, compared to an analytic CDF.

    807:   Collective

    809:   Input Parameters:
    810: + v   - The data vector, blocksize is the sample dimension
    811: . w   - The vector of weights for each sample, instead of the default 1/n
    812: - cdf - The analytic CDF

    814:   Output Parameter:
    815: . alpha - The KS statistic

    817:   Level: advanced

    819:   Notes:
    820:   The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is

    822:   $$
    823:   D_n = \sup_x \left| F_n(x) - F(x) \right|
    824:   $$

    826:   where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by

    828:   $$
    829:   F_n =  # of samples <= x / n
    830:   $$

    832:   The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
    833:   cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
    834:   the largest absolute difference between the two distribution functions across all $x$ values.

    836:   The goodness-of-fit test, or Kolmogorov-Smirnov test, is constructed using the Kolmogorov
    837:   distribution. It rejects the null hypothesis at level $\alpha$ if

    839:   $$
    840:   \sqrt{n} D_{n} > K_{\alpha},
    841:   $$

    843:   where $K_\alpha$ is found from

    845:   $$
    846:   \operatorname{Pr}(K \leq K_{\alpha}) = 1 - \alpha.
    847:   $$

    849:   This means that getting a small alpha says that we have high confidence that the data did not come
    850:   from the input distribution, so we say that it rejects the null hypothesis.

    852: .seealso: `PetscProbComputeKSStatistic()`, `PetscProbComputeKSStatisticMagnitude()`, `PetscProbFn`
    853: @*/
    854: PetscErrorCode PetscProbComputeKSStatisticWeighted(Vec v, Vec w, PetscProbFn *cdf, PetscReal *alpha)
    855: {
    856:   PetscViewer        viewer     = NULL;
    857:   OutputType         outputType = NONE;
    858:   const PetscScalar *val, *wgt;
    859:   PetscInt           n;

    861:   PetscFunctionBegin;
    862:   PetscCall(KSViewerCreate((PetscObject)v, &outputType, &viewer));
    863:   PetscCall(VecGetLocalSize(v, &n));
    864:   PetscCall(VecGetArrayRead(v, &val));
    865:   PetscCall(VecGetArrayRead(w, &wgt));
    866:   PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "K-S test does not support complex");
    867:   PetscCall(PetscProbComputeKSStatistic_Internal(PetscObjectComm((PetscObject)v), n, (PetscReal *)val, (PetscReal *)wgt, cdf, alpha, outputType, viewer));
    868:   PetscCall(VecRestoreArrayRead(v, &val));
    869:   PetscCall(VecRestoreArrayRead(w, &wgt));
    870:   PetscCall(KSViewerDestroy(&viewer));
    871:   PetscFunctionReturn(PETSC_SUCCESS);
    872: }

    874: /*@C
    875:   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.

    877:   Collective

    879:   Input Parameters:
    880: + v   - The data vector, blocksize is the sample dimension
    881: - cdf - The analytic CDF

    883:   Output Parameter:
    884: . alpha - The KS statistic

    886:   Level: advanced

    888:   Notes:
    889:   The Kolmogorov-Smirnov statistic for a given cumulative distribution function $F(x)$ is

    891:   $$
    892:   D_n = \sup_x \left| F_n(x) - F(x) \right|
    893:   $$

    895:   where $\sup_x$ is the supremum of the set of distances, and the empirical distribution function $F_n(x)$ is discrete, and given by

    897:   $$
    898:   F_n =  # of samples <= x / n
    899:   $$

    901:   The empirical distribution function $F_n(x)$ is discrete, and thus had a ``stairstep''
    902:   cumulative distribution, making $n$ the number of stairs. Intuitively, the statistic takes
    903:   the largest absolute difference between the two distribution functions across all $x$ values.

    905:   The goodness-of-fit test, or Kolmogorov-Smirnov test, is constructed using the Kolmogorov
    906:   distribution. It rejects the null hypothesis at level $\alpha$ if

    908:   $$
    909:   \sqrt{n} D_{n} > K_{\alpha},
    910:   $$

    912:   where $K_\alpha$ is found from

    914:   $$
    915:   \operatorname{Pr}(K \leq K_{\alpha}) = 1 - \alpha.
    916:   $$

    918:   This means that getting a small alpha says that we have high confidence that the data did not come
    919:   from the input distribution, so we say that it rejects the null hypothesis.

    921: .seealso: `PetscProbComputeKSStatistic()`, `PetscProbComputeKSStatisticWeighted()`, `PetscProbFn`
    922: @*/
    923: PetscErrorCode PetscProbComputeKSStatisticMagnitude(Vec v, PetscProbFn *cdf, PetscReal *alpha)
    924: {
    925:   PetscViewer        viewer     = NULL;
    926:   OutputType         outputType = NONE;
    927:   const PetscScalar *a;
    928:   PetscReal         *speed;
    929:   PetscInt           n, dim;

    931:   PetscFunctionBegin;
    932:   PetscCall(KSViewerCreate((PetscObject)v, &outputType, &viewer));
    933:   // Convert to a scalar value
    934:   PetscCall(VecGetLocalSize(v, &n));
    935:   PetscCall(VecGetBlockSize(v, &dim));
    936:   n /= dim;
    937:   PetscCall(PetscMalloc1(n, &speed));
    938:   PetscCall(VecGetArrayRead(v, &a));
    939:   for (PetscInt p = 0; p < n; ++p) {
    940:     PetscReal mag = 0.;

    942:     for (PetscInt d = 0; d < dim; ++d) mag += PetscSqr(PetscRealPart(a[p * dim + d]));
    943:     speed[p] = PetscSqrtReal(mag);
    944:   }
    945:   PetscCall(VecRestoreArrayRead(v, &a));
    946:   // Compute statistic
    947:   PetscCall(PetscProbComputeKSStatistic_Internal(PetscObjectComm((PetscObject)v), n, speed, NULL, cdf, alpha, outputType, viewer));
    948:   PetscCall(PetscFree(speed));
    949:   PetscCall(KSViewerDestroy(&viewer));
    950:   PetscFunctionReturn(PETSC_SUCCESS);
    951: }