Coverage for /usr/local/lib/python3.10/site-packages/spam/excursions/randomFields.py: 100%

1# Library of SPAM random field generation functions based on R. 

2# Copyright (C) 2020 SPAM Contributors 

3#  

4# This program is free software: you can redistribute it and/or modify it 

5# under the terms of the GNU General Public License as published by the Free 

6# Software Foundation, either version 3 of the License, or (at your option) 

7# any later version. 

8#  

9# This program is distributed in the hope that it will be useful, but WITHOUT 

10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

11# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

12# more details. 

13#  

14# You should have received a copy of the GNU General Public License along with 

15# this program. If not, see <http://www.gnu.org/licenses/>. 

16""" 

17This module is a wrapper around the `python library gstools <https://github.com/GeoStat-Framework/GSTools>`_ made to intuitively generate realisations of correlated Random Fields on a regular grid with a given size, discretisation and characteristic length. 

18 

19""" 

20import numpy 

21import gstools 

22import spam.helpers 

23import time 

24 

25def simulateRandomField( 

26 lengths=1.0, 

27 nNodes=100, 

28 covarianceModel='Gaussian', 

29 covarianceParameters={}, 

30 dim=3, 

31 nRea=1, 

32 seed=None, 

33 vtkFile=None, 

34): 

35 """Generates realisations of Correlated Random Field based on the `python library gstools <https://github.com/GeoStat-Framework/GSTools>`_. 

36  

37 Parameters 

38 ---------- 

39 lengths: float | list of floats, default=1.0 

40 Length of the definition domain in every directions (should be size of dim). If a single float is given it will be used in each directions. 

41 nNodes: int | list of int, default=100 

42 Number of nodes in every directions. If a single float is given it will be used in each directions. 

43 covarianceModel: string, default="Gaussian" 

44 The name of the covariance model used by gstools. 

45 To be chosen among `this list of models <https://geostat-framework.readthedocs.io/projects/gstools/en/stable/api/gstools.covmodel.html>`_. 

46 covarianceParameters: dict, default={} 

47 A dictionnary of the keyword arguments of the covariance model used (see gstools documentation). 

48 dim: int default=3 

49 Spatial dimention 

50 nRea: int, default=1 

51 number of realisations 

52 seed: int, default=None 

53 Defines the seed for the random numbers generator. Setting a seed to a certain integer yields the same realisations. 

54 vtkFile: string, default=None 

55 If not None, the base name of the vtk files saved. 

56 

57 Returns 

58 ------- 

59 List of `nRea` numpy arrays if shape nNodes: 

60 The realisations of the random fields. 

61 

62 Example 

63 ------- 

64 >>> from spam.excursions import simulateRandomField 

65 >>> 

66 >>> # generation of two realisations of a Gaussian Correlated Random field 

67 >>> realisations = simulateRandomField( 

68 >>> lengths=25, 

69 >>> nNodes=100, 

70 >>> dim=3, 

71 >>> nRea=2, 

72 >>> covarianceModel="Gaussian", 

73 >>> covarianceParameters={"len_scale": 10}, 

74 >>> vtkFile="gaussian_randomfield" 

75 >>> ) 

76 >>> 

77 >>> for realisation in realisations: 

78 >>> print(realisation.shape) 

79 >>> 

80 >>> # generation of one realisation of a Matern Correlated Random field 

81 >>> realisations = simulateRandomField( 

82 >>> lengths=25, 

83 >>> nNodes=100, 

84 >>> dim=3, 

85 >>> covarianceModel="Matern", 

86 >>> covarianceParameters={"len_scale": 10, "nu": 0.2}, 

87 >>> vtkFile="matern_randomfield" 

88 >>> ) 

89 

90 Warning 

91 ------- 

92 If no rescaling factors are given in the covariance parameters it is set to 1 which deviates from gstools default behavior. 

93  

94 """ 

95 

96 def _castToList(v, d, types): 

97 return [v] * dim if isinstance(v, types) else v 

98 

99 # cast lengths into a list 

100 lengths = _castToList(lengths, dim, (int, float)) 

101 nNodes = _castToList(nNodes, dim, int) 

102 

103 if len(lengths) != dim: 

104 raise ValueError(f"lengths doesn't have the good size with regards to dim {len(lengths)} != {dim}.") 

105 

106 # Parameters 

107 # print("Random fields parameters") 

108 # print(f'\tLengths {lengths}') 

109 # print(f'\tNumber of nodes {nNodes}') 

110 # print(f'\tNumber of realisations {nRea}') 

111 

112 # model 

113 if "rescale" not in covarianceParameters: 

114 covarianceParameters["rescale"] = 1.0 

115 model = getattr(gstools.covmodel, covarianceModel)(dim=dim, **covarianceParameters) 

116 # print(f'Covariance model ({covarianceModel})') 

117 # print(f'\t{model}') 

118 

119 # generator 

120 srf = gstools.SRF(model) 

121 pos = [numpy.linspace(0, l, n) for n, l in zip(nNodes, lengths)] 

122 srf.set_pos(pos, "structured") 

123 seed = gstools.random.MasterRNG(seed) 

124 

125 for i in range(nRea): 

126 tic = time.perf_counter() 

127 name = f'{vtkFile if vtkFile else "Field"}_{i:04d}' 

128 print(f"Generating {name}...", end=" ") 

129 srf(seed=seed(), store=name) 

130 print(f"{time.perf_counter() - tic:.2f} seconds") 

131 if vtkFile is not None: 

132 srf.vtk_export(name, name) 

133 

134 return srf.all_fields 

135 

136 

137def parametersLogToGauss(muLog, stdLog, lcLog=1): 

138 """ Gives the underlying Gaussian standard deviation and mean value 

139 of the log normal distribution of standard deviation and mean value. 

140 

141 Parameters 

142 ---------- 

143 muLog: float 

144 Mean value of the log normal distribution. 

145 stdLog: float 

146 Standard deviation of the log normal distribution. 

147 lcLog: float, default=1 

148 Correlation length of the underlying log normal correlated Random Field. 

149 

150 Returns 

151 ------- 

152 muGauss: float 

153 Mean value of the Gaussian distribution. 

154 stdGauss: float 

155 Standard deviation of the Gaussian distribution. 

156 lcGauss: float 

157 Correlation length of the underlying Gaussian correlated Random Field. 

158 

159 

160 """ 

161 stdGauss = numpy.sqrt(numpy.log(1.0 + stdLog**2 / muLog**2)) 

162 muGauss = numpy.log(muLog) - 0.5 * stdGauss**2 

163 lcGauss = (-1.0 / numpy.log(numpy.log(1 + stdLog**2 * numpy.exp(-1.0 / lcLog) / muLog**2) / stdLog**2))**(0.5) 

164 

165 return muGauss, stdGauss, lcGauss 

166 

167 

168def fieldGaussToUniform(g, a=0.0, b=1.0): 

169 """Transforms a Gaussian Random Field into a uniform Random Field. 

170 

171 Parameters 

172 ---------- 

173 g: array 

174 The values of the Gaussian Random Fields. 

175 a: float, default=0.0 

176 The minimum borne of the uniform distribution. 

177 b: float, default=1.0 

178 The maximum borne of the uniform distribution. 

179 

180 Return 

181 ------ 

182 array: 

183 The uniform Random Field. 

184 """ 

185 from scipy.special import erf 

186 return float(a) + 0.5 * float(b) * (1.0 + erf(g / 2.0**0.5)) 

187 

188 

189def fieldUniformToGauss(g, a=0.0, b=1.0): 

190 """Transforms a uniform Random Field into a Gaussian Random Field. 

191 

192 Parameters 

193 ---------- 

194 g: array 

195 The values of the uniform Random Fields. 

196 a: float, default=0.0 

197 The minimum borne of the uniform distribution. 

198 b: float, default=1.0 

199 The maximum borne of the uniform distribution. 

200 

201 Return 

202 ------ 

203 array: 

204 The Gaussian Random Field. 

205 """ 

206 from scipy.special import erfinv 

207 return 2.0**0.5 * erfinv(2.0 * (g - float(a)) / float(b) - 1.0)