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

1"""

2Library of SPAM random field generation functions based on R.

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

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

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

8any later version.

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

11ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

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

13more details.

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

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

17"""

19import numpy

20import tifffile

21import spam.helpers

24def simulateRandomField(lengths=(1, 1, 1), nNodes=(10, 10, 10), covariance={'type': 'stable', 'alpha': 2},

25 seed='NA', method=None, maxGB=1, dim=2, nRea=1, data=None, given=None,

26 tifFile=None, fieldFile=None, vtkFile=None,

27 verbose=False, RprintLevel=0):

28 """Generate a Correlated Random Field based on the package RandomFields (cran-R). Wrapping is handled by the module rpy2.

30 * cran-R: https://cran.r-project.org/

31 * Package RandomFields (Martin Schlather): https://cran.r-project.org/web/packages/RandomFields/index.html

32 * rpy2: https://pypi.python.org/pypi/rpy2

34 Parameters

35 ----------

36 lengths : array of float, default=(1,1,1)

37 Length of the definition domain in every directions (should be size of dim).

38 nNodes : array of int, default=(10,10,10)

39 Number of nodes in every directions (should be size of dim).

40 covariance : dictionary (default {'type':'stable', 'alpha':2})

41 Sets the type of covariance and its parameters. It defines the 'RMmodel'. Note: stable with alpha = 2 is the Gaussian covariance function.

42 seed : string, default='NA'

43 Defines the seed for the generator of random numbers. If seed=0, all realisations are the same.

44 method : string, default=None

45 Describe the RP method for R. If `None` given the default behavior is used (R looks for the best one). Works with method=`circulant` (FFT needs a lot of memory) method=`tbm` (turning band method)

46 dim : int

47 spatial dimention

48 nRea: int, default=1

49 number of realisations

50 given: string

51 string input for given in RFsimulation for conditional simulations. Ex. given='c(0,1)' sets node at x=0 and x=1.

52 data: string

53 string input for data in RFsimulation for conditional simulations. Ex. data='c(0,0)' sets values for node 1 and 2 at 0.

54 maxGB: float,default=1

55 set the maximum memory that R can use.

56 RprintLevel: int, default=0

57 The print level of R. The higher it is the more information is printed.

58 verbose: bool, default=False

59 Verbose mode if True.

60 tifFile: string, default=None

61 If not None, the base name of the tif files saved.

62 fieldFile: string, default=None

63 If not None, the base name of the field files saved. One file is saved for each realisation.

64 These file can be used as input for ``spam.mesh.projection``.

65 vtkFile: string, default=None

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

68 Returns

69 -------

70 array:

71 The realisations of the random fields. ``shape = (nx, ny, nz, nRea)``

73 """

75 # rpy2

76 import rpy2.robjects.packages as rpackages

77 from rpy2.robjects import r

79 if verbose:

80 print('*** Welcome to spam.excursions.randomfields.simulateRandomField()')

81 # step 0: defines shortcuts

82 if dim > 1:

83 le = (lengths, lengths, lengths) if isinstance(lengths, float) else lengths

84 n = (nNodes, nNodes, nNodes) if isinstance(nNodes, int) else nNodes

85 nCells = [_ - 1 for _ in n]

86 else:

87 le = lengths if isinstance(lengths, (int, float)) else lengths[0]

88 n = nNodes if isinstance(nNodes, (int, float)) else nNodes[0]

89 nCells = n - 1

91 if verbose:

92 print('*\tLengths : {}'.format(le))

93 print('*\tNumber of nodes: {}'.format(n))

94 print('*\tNumber of cells: {}'.format(nCells))

95 print('*\tCovariance:')

96 for key, value in covariance.items():

97 print('*\t\t {}: {}'.format(key, value))

98 print('*')

100 # Step 1: R and RandomFields preparation

101 # Step 1.1: import package

102 if verbose:

103 print('* Import R package RandomFields')

104 try:

105 rpackages.importr("RandomFields")

106 except BaseException: # no cover

107 if verbose:

108 print('*\tFail import.')

109 if verbose:

110 print('*\tTrying to install package RandomFields.')

111 # from rpy2.robjects.vectors import StrVector

112 utils = rpackages.importr('utils')

113 utils.chooseCRANmirror(ind=1)

114 utils.install_packages('RandomFields')

115 rpackages.importr('RandomFields')

116 if verbose:

117 print('*\tPackage successfully installed and loaded.')

118

119 # Step 1.2: RMmodel

120 try:

121 if covariance['type'] == 'stable':

122 if 'variance' not in covariance:

123 covariance['variance'] = 1.0

124 if 'correlation_length' not in covariance:

125 covariance['correlation_length'] = 1.0

126 if 'alpha' not in covariance:

127 covariance['alpha'] = 2

128 _r_RMmodel = 'RMstable(var={variance}, scale={correlation_length}, alpha={alpha})'.format(**covariance)

129 elif covariance['type'] == 'matern':

130 if 'variance' not in covariance:

131 covariance['variance'] = 1

132 if 'correlation_length' not in covariance:

133 covariance['correlation_length'] = 1

134 if 'nu' not in covariance:

135 covariance['nu'] = 2

136 _r_RMmodel = 'RMmatern(nu={nu}, var={variance}, scale={correlation_length})'.format(**covariance)

137 elif covariance['type'] == 'bessel':

138 if 'variance' not in covariance:

139 covariance['variance'] = 1

140 if 'correlation_length' not in covariance:

141 covariance['correlation_length'] = 1

142 if 'nu' not in covariance:

143 covariance['nu'] = 2

144 _r_RMmodel = 'RMbessel(nu={nu}, var={variance}, scale={correlation_length})'.format(**covariance)

145 elif covariance['type'] == 'dampedcos':

146 if 'variance' not in covariance:

147 covariance['variance'] = 1

148 if 'correlation_length' not in covariance:

149 covariance['correlation_length'] = 1

150 if 'lambda' not in covariance:

151 covariance['lambda'] = 0

152 _r_RMmodel = 'RMdampedcos(lambda={lambda},var={variance}, scale={correlation_length})'.format(**covariance)

153 elif covariance['type'] == 'sinepower':

154 if 'variance' not in covariance:

155 covariance['variance'] = 1

156 if 'correlation_length' not in covariance:

157 covariance['correlation_length'] = 1

158 if 'alpha' not in covariance:

159 covariance['alpha'] = 2

160 _r_RMmodel = 'RMstable(var={variance}, scale={correlation_length}, alpha={alpha})'.format(**covariance)

161 except BaseException:

162 _r_RMmodel = covariance

163 if verbose:

164 print('*\tRMmodel: {}'.format(_r_RMmodel))

165

166 # Step 1.3: RPmodel (add circulant if necessary)

167 if method is not None:

168 _r_RPmodel = "RP{method}({rm})".format(method=method, rm=_r_RMmodel)

169 else:

170 _r_RPmodel = _r_RMmodel

171 if verbose:

172 print('*\tRPmodel: {}'.format(_r_RPmodel))

173

174 # Step 1.4: Options

175 _r_printlevel = RprintLevel

176 _r_options = 'RFoptions(seed={},spConform=FALSE,cPrintlevel={})'.format(

177 seed, _r_printlevel)

178 if verbose:

179 print('*\tOptions: {}'.format(_r_options))

180

181 # Step 1.5: Grid, RFsimulate and sequence

182 _r_grid = 'NULL'

183 _r_RFsimulate = 'RFsimulate'

184 if dim == 1:

185 _r_x = 'seq(0.0, {}, {})'.format(le, float(le) / (n - 1))

186 else:

187 _r_x = 'seq(0.0, {}, {})'.format(le[0], float(le[0]) / (n[0] - 1))

188 if dim > 1:

189 _r_y = 'seq(0.0, {}, {})'.format(le[1], float(le[1]) / (n[1] - 1))

190 if dim > 2:

191 _r_z = 'seq(0.0, {}, {})'.format(le[2], float(le[2]) / (n[2] - 1))

192

193 # step 2: transform _r_string into r objects

194 if verbose:

195 print('* Convert R object')

196 rpmo = r(_r_RPmodel)

197 rfsi = r(_r_RFsimulate)

198 opti = r(_r_options)

199 grid = r(_r_grid)

200 if dim > 0:

201 x = r(_r_x)

202 if dim > 1:

203 y = r(_r_y)

204 if dim > 2:

205 z = r(_r_z)

206

207 if data is not None:

208 data = r(data)

209 else:

210 data = r('NULL')

211 if given is not None:

212 given = r(given)

213 else:

214 given = r('NULL')

215

216 # Step 3: launch r command

217 if verbose:

218 print('* Simulate Random Field')

219 if verbose:

220 print('<R>')

221 if dim == 1:

222 rf = numpy.array(rfsi(model=rpmo, x=x, grid=grid, n=nRea, given=given, data=data)).astype('<f4')

223 if dim == 2:

224 rf = numpy.array(rfsi(model=rpmo, x=x, y=y, grid=grid, n=nRea)).astype('<f4')

225 if dim == 3:

226 rf = numpy.array(rfsi(model=rpmo, x=x, y=y, z=z, grid=grid, n=nRea, maxGB=maxGB)).astype('<f4')

227

228 if verbose:

229 print('</R>')

230

231 # save one tif for each realisation

232 if tifFile:

233 if verbose:

234 print('* Save image(s):')

235 for i in range(nRea):

236 currentName = '{tifFile}_{rea:05d}.tif'.format(tifFile=tifFile, rea=i)

237 if verbose:

238 print('*\t{}'.format(currentName))

239 if dim == 1:

240 if verbose:

241 print('# /!\ Impossible do make tif with 1D Fields')

242 elif dim == 2:

243 if nRea == 1:

244 tifffile.imwrite(currentName, rf.astype('<f4'))

245 else:

246 tifffile.imwrite(currentName, rf[:, :, i].astype('<f4'))

247 elif dim == 3:

248 if nRea == 1:

249 tifffile.imwrite(currentName, rf.astype('<f4'))

250 else:

251 tifffile.imwrite(currentName, rf[:, :, :, i].astype('<f4'))

252

253 # save dat files with good format for projmorpho ## only 3D

254 if fieldFile and dim == 3:

255 if verbose:

256 print('* Save realisations in field files:')

257 for i in range(nRea):

258 currentName = '{o}_{i:05d}'.format(o=fieldFile, i=i)

259 if verbose:

260 print('*\t{}.dat'.format(currentName))

261 with open('{}.dat'.format(currentName), 'w') as f:

262 f.write('{}, {}, {}\n'.format(*le))

263 f.write('0,0,0\n')

264 f.write('{}, {}, {}\n'.format(*nCells))

265 if nRea == 1:

266 for val in rf.flatten():

267 f.write('{}\n'.format(str(val)))

268 else:

269 for val in rf[:, :, :, i].flatten():

270 f.write('{}\n'.format(str(val)))

271

272 # save dat files with good format for projmorpho ## only 3D

273 if vtkFile and dim == 3:

274 if verbose:

275 print('* Save realisations in several vtk files:')

276 for i in range(nRea):

277 currentName = '{o}_{i:05d}'.format(o=vtkFile, i=i)

278 if verbose:

279 print('*\t{}.vtk'.format(currentName))

280

281 if nRea == 1:

282 # tmp = rf.ravel()

283 tmp = rf

284 else:

285 # tmp = rf[:, :, :, i].ravel()

286 tmp = rf[:, :, :, i]

287 spam.helpers.writeStructuredVTK(aspectRatio=[float(a) / float(b) for a, b in zip(le, nCells)], pointData={'RandomField': tmp}, fileName="{}.vtk".format(currentName))

288

289 if verbose:

290 print('*')

291 print('*** See you soon.')

292

293 return rf

294

295

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

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

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

299

300 Parameters

301 ----------

302 muLog: float

303 Mean value of the log normal distribution.

304 stdLog: float

305 Standard deviation of the log normal distribution.

306 lcLog: float, default=1

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

308

309 Returns

310 -------

311 muGauss: float

312 Mean value of the Gaussian distribution.

313 stdGauss: float

314 Standard deviation of the Gaussian distribution.

315 lcGauss: float

316 Correlation length of the underlying Gaussian correlated Random Field.

317

318

319 """

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

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

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

323

324 return muGauss, stdGauss, lcGauss

325

326

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

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

329

330 Parameters

331 ----------

332 g: array

333 The values of the Gaussian Random Fields.

334 a: float, default=0.0

335 The minimum borne of the uniform distribution.

336 b: float, default=1.0

337 The maximum borne of the uniform distribution.

338

339 Return

340 ------

341 array:

342 The uniform Random Field.

343 """

344 from scipy.special import erf

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

346

347

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

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

350

351 Parameters

352 ----------

353 g: array

354 The values of the uniform Random Fields.

355 a: float, default=0.0

356 The minimum borne of the uniform distribution.

357 b: float, default=1.0

358 The maximum borne of the uniform distribution.

359

360 Return

361 ------

362 array:

363 The Gaussian Random Field.

364 """

365 from scipy.special import erfinv

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