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

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

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