source: issm/trunk-jpl/src/c/shared/Random/random.cpp@ 26656

Last change on this file since 26656 was 26621, checked in by vverjans, 3 years ago

CHG: making random number generation in ISSM dependent on classes of randomgenerator.h (should allow cross-platform tests of stochasticity)
  • Property svn:executable set to *
File size: 3.2 KB
Line 
1/*!\file: random
2 * \brief random number generating functions
3 */
4
5/*Headers*/
6/*{{{*/
7#include <stdio.h>
8#include <sys/types.h>
9#include <math.h>
10#include <float.h> /* DBL_EPSILON */
11#include <cstdarg>
12#include <iostream>
13
14#include "../Matrix/matrix.h"
15#include "../Exceptions/exceptions.h"
16#include "../MemOps/MemOps.h"
17#include "../io/io.h"
18#include "./randomgenerator.h"
19/*}}}*/
20
21void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev, int seed=-1) { /*{{{*/
22
23 /*Seed the pseudo-random number generator*/
24 rnd::linear_congruential_engine randomengine;
25 randomengine.seed(seed);
26 /*Normal distribution*/
27 rnd::normal_distribution distriNormal(mean,sdev);
28 *prand = distriNormal.generator(randomengine);
29} /*}}}*/
30void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
31
32 IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
33 IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
34 IssmDouble* Lchol = xNewZeroInit<IssmDouble>(dim*dim);
35
36 /*True randomness if seed<0, otherwise random seed is fixed at seed*/
37 /*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
38 rnd::linear_congruential_engine randomengine;
39 randomengine.seed(seed);
40 /*Normal distribution*/
41 rnd::normal_distribution distriNormal(0.0,1.0);
42 for(int i=0;i<dim;i++){
43 sampleStandardNormal[i] = distriNormal.generator(randomengine);
44 }
45
46 /*Cholsesky decomposition of the covariance matrix*/
47 CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
48
49 /*Matrix by vector multiplication*/
50 for(int i=0;i<dim;i++){
51 /*Entry-by-entry multiplication along matrix row*/
52 IssmDouble sum=0.;
53 for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j];
54 sampleMultivariateNormal[i] = mean+sum;
55 }
56
57 /*Assign output pointer and cleanup*/
58 *prand = sampleMultivariateNormal;
59 xDelete<IssmPDouble>(sampleStandardNormal);
60 xDelete<IssmDouble>(Lchol);
61} /*}}}*/
62void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
63
64 IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
65 IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
66 IssmDouble* Lchol = xNewZeroInit<IssmDouble>(dim*dim);
67
68 /*True randomness if seed<0, otherwise random seed is fixed at seed*/
69 /*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
70 rnd::linear_congruential_engine randomengine;
71 randomengine.seed(seed);
72 /*Normal distribution*/
73 rnd::normal_distribution distriNormal(0.0,1.0);
74 for(int i=0;i<dim;i++){
75 sampleStandardNormal[i] = distriNormal.generator(randomengine);
76 }
77
78 /*Cholsesky decomposition of the covariance matrix*/
79 CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
80
81 /*Matrix by vector multiplication*/
82 for(int i=0;i<dim;i++){
83 IssmDouble sum = 0.;
84 for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j];
85 sampleMultivariateNormal[i] = mean[i]+sum;
86 }
87
88 /*Assign output pointer and cleanup*/
89 *prand = sampleMultivariateNormal;
90 xDelete<IssmPDouble>(sampleStandardNormal);
91 xDelete<IssmDouble>(Lchol);
92} /*}}}*/
93
94
95
Note: See TracBrowser for help on using the repository browser.