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

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

CHG: deleting old autoregression scripts
  • Property svn:executable set to *
File size: 3.5 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 /*Assign output pointer and cleanup*/
29 *prand = distriNormal.generator(randomengine);
30 randomengine.free_resources();
31} /*}}}*/
32void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
33
34 IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
35 IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
36 IssmDouble* Lchol = xNewZeroInit<IssmDouble>(dim*dim);
37
38 /*True randomness if seed<0, otherwise random seed is fixed at seed*/
39 /*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
40 rnd::linear_congruential_engine randomengine;
41 randomengine.seed(seed);
42 /*Normal distribution*/
43 rnd::normal_distribution distriNormal(0.0,1.0);
44 for(int i=0;i<dim;i++){
45 sampleStandardNormal[i] = distriNormal.generator(randomengine);
46 //_printf_("VV i sampleStandardNormal[i]: "<<i<<" "<<sampleStandardNormal[i]<<'\n');
47 }
48
49 /*Cholsesky decomposition of the covariance matrix*/
50 CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
51
52 /*Matrix by vector multiplication*/
53 for(int i=0;i<dim;i++){
54 /*Entry-by-entry multiplication along matrix row*/
55 IssmDouble sum=0.;
56 for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j];
57 sampleMultivariateNormal[i] = mean+sum;
58 }
59
60 /*Assign output pointer and cleanup*/
61 *prand = sampleMultivariateNormal;
62 xDelete<IssmPDouble>(sampleStandardNormal);
63 xDelete<IssmDouble>(Lchol);
64 randomengine.free_resources();
65} /*}}}*/
66void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
67
68 IssmPDouble* sampleStandardNormal = xNew<IssmPDouble>(dim);
69 IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
70 IssmDouble* Lchol = xNewZeroInit<IssmDouble>(dim*dim);
71
72 /*True randomness if seed<0, otherwise random seed is fixed at seed*/
73 /*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
74 rnd::linear_congruential_engine randomengine;
75 randomengine.seed(seed);
76 /*Normal distribution*/
77 rnd::normal_distribution distriNormal(0.0,1.0);
78 for(int i=0;i<dim;i++){
79 sampleStandardNormal[i] = distriNormal.generator(randomengine);
80 //_printf_("VV i sampleStandardNormal[i]: "<<i<<" "<<sampleStandardNormal[i]<<'\n');
81 }
82
83 /*Cholsesky decomposition of the covariance matrix*/
84 CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
85
86 /*Matrix by vector multiplication*/
87 for(int i=0;i<dim;i++){
88 IssmDouble sum = 0.;
89 for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j];
90 sampleMultivariateNormal[i] = mean[i]+sum;
91 }
92
93 /*Assign output pointer and cleanup*/
94 *prand = sampleMultivariateNormal;
95 xDelete<IssmPDouble>(sampleStandardNormal);
96 xDelete<IssmDouble>(Lchol);
97 randomengine.free_resources();
98} /*}}}*/
99
100
101
Note: See TracBrowser for help on using the repository browser.