Coverage for /usr/local/lib/python3.8/site-packages/spam/excursions/elkc.py: 87%
93 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1"""
2Library of SPAM functions for computation Lipschitz-Killing Curvatures.
3Copyright (C) 2020 SPAM Contributors
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
20from scipy.constants import pi
22def expectedMesures(kappa, j, n, hittingSet='tail', distributionType='gauss', mu=0.0, std=1.0, correlationType='gauss', lc=1.0, nu=2.0, a=1.0):
23 """
24 Compute the Lipschitz-Killing Curvatures E(LKC)
26 Parameters
27 -----------
28 kappa : float or list
29 value of the threshold
30 j : int
31 number of the functionnal
32 n : int
33 spatial dimension
34 hittingSet : string
35 the hitting set
36 distributionType : string
37 type of distribution (see gaussianMinkowskiFunctionals)
38 mu : float
39 mean value of the distribution
40 std : float
41 standard deviation of the distribution
42 correlationType : string
43 type of correlation function
44 lc : float
45 correlation length
46 nu : float
47 parameter used for correlationType='matern'
48 a : float or list of floats
49 size(s) of the object
50 if float, the object is considered as a cube
52 Returns
53 --------
54 Same type as kappa
55 the Gaussian Minkowski functionnal
56 """
58 e = 0.0 * kappa
59 # compute second spetral moment
60 lam2 = secondSpectralMoment(std, lc, correlationType=correlationType, nu=nu)
61 # loop over i
62 for i in range(n - j + 1):
63 # comute LKC of the cube
64 lkc = lkcCuboid(i + j, n, a)
65 # compute GMF
66 gmf = gaussianMinkowskiFunctionals(kappa, i, hittingSet=hittingSet, distributionType=distributionType, mu=mu, std=std)
67 # compute
68 fla = flag(i + j, i)
69 e = e + fla * lkc * gmf * (lam2 / (2.0 * pi))**(i / 2.0)
71 return e
74def gaussianMinkowskiFunctionals(kappa, j, hittingSet='tail', distributionType='gauss', mu=0.0, std=1.0):
75 """
76 Evaluate the Gaussian Minkowski functionals j
78 Parameters
79 -----------
80 kappa : float or list of floats
81 value(s) of the threshold(s)
82 j : int
83 number of the functionnal
84 hittingSet : string
85 hitting set
86 distributionType : string
87 type of the underlying distribution. Implemented: ``gauss`` and ``log``
88 mu : float
89 mean value of the distribution
90 std : float
91 standard deviation of the distribution
93 Returns
94 --------
95 mink : same type as kappa
96 the Gaussian Minkowski functionnal
97 """
98 from scipy.special import erf
100 if hittingSet == 'tail':
101 sign = 1.0
102 elif hittingSet == 'head':
103 sign = (-1.0)**(j + 1)
104 else:
105 print('error in spam.excursions.elkc.gaussianMinkowskiFunctionals(): hitting set {} not implemented.'.format(hittingSet))
106 return 0
108 # STEP 2: change variable
109 if distributionType == 'gauss':
110 k = (kappa - mu) / std
111 elif distributionType == 'log':
112 k = (numpy.log(kappa) - mu) / std
113 else:
114 print('error in spam.excursions.elkc.gaussianMinkowskiFunctionals(): distribution type {} not implemented.'.format(distributionType))
115 return 0
117 if j == 0:
118 mink = 0.5 * (1.0 - sign * erf(k / 2**0.5))
119 elif j > 0:
120 mink = sign * numpy.exp(-k**2.0 / 2.0) * hermitePolynomial(k, j - 1) / (std**j * (2.0 * pi)**0.5)
121 # mink = - (k**3 - 3*k)*numpy.exp(-k**2.0/2.0) / (std**4*(2.0*pi)**0.5)
122 else:
123 print('error in spam.excursions.elkc.gaussianMinkowskiFunctionals(): j should be > 0. {} given.'.format(j))
124 return 0
126 return mink
129def secondSpectralMoment(std, lc, correlationType='gauss', nu=2.0):
130 """
131 Compute the second spectral moment for a given covariance function:
133 .. math::
134 {\lambda}_2 = C''(h)|_{h=0}
136 Parameters
137 -----------
138 std : float
139 standard deviation
140 lc : float
141 correlation length
142 correlationType : string
143 type of correlation function
144 nu : float
145 parameter used for correlationType='matern'
147 Returns
148 --------
149 float
150 the second spectral moment
151 """
152 v = std**2.0
153 if correlationType == 'gauss':
154 return 2.0 * v / lc**2.0
155 elif correlationType == 'matern':
156 if nu <= 1:
157 print('error in spam.excursions.elkc.secondSpectralMoment(): nu should be > 1. {} given.'.format(nu))
158 return 0
159 return (v * nu) / (lc**2.0 * (nu - 1.0))
160 else:
161 print('error in spam.excursions.elkc.secondSpectralMoment(): correlation type {} not implemented.'.format(correlationType))
162 return 0
165def flag(n, j):
166 """
167 Compute the flag coefficients ``[ n, j ] = ( n, j ) w_n / w_{n-j}w_j``
169 Parameters
170 -----------
171 n : int
172 first parameter of the binom
173 j : int
174 second parameter of the binom
176 Returns
177 --------
178 float
179 the flag coefficient [n , j]
180 """
181 from scipy.special import binom
183 if j > n:
184 print('error in spam.excursions.elkc.flag(): n should be >= j. Here n={}, j={}'.format(n, j))
185 return 0
186 vn = ballVolume(n)
187 vj = ballVolume(j)
188 vnj = ballVolume(n - j)
190 return binom(n, j) * vn / (vnj * vj)
193def hermitePolynomial(x, n):
194 """
195 Evaluate a x the probabilitic Hermite polynomial n
197 Parameters
198 -----------
199 x : float
200 point of evaluation
201 n : int
202 number of Hermite polynomia
204 Returns
205 --------
206 float
207 ``He_n(x)`` the value of the polynomial
209 """
210 from numpy.polynomial.hermite_e import hermeval
212 coefs = [0] * (n + 1)
213 coefs[n] = 1
215 return hermeval(x, coefs)
218def ballVolume(n, r=1.0):
219 """
220 Compute the volume of a nth dimensional ball
222 Parameters
223 -----------
224 n : int
225 spatial dimension
226 r : float
227 radius of the ball
229 Returns
230 --------
231 float
232 volume of the ball
233 """
234 from scipy.special import gamma
236 return float(r)**float(n) * pi**(0.5 * n) / gamma(1.0 + 0.5 * n)
239def lkcCuboid(i, n, a):
240 """
241 Compute the Lipschitz-Killing Curvatures of a parallelepiped
243 Parameters
244 -----------
245 i : int
246 number of the LKC. ``0 <= i <= n``
247 n : int
248 spatial dimension
249 a : float or list of float
250 size(s) of the cuboid. If float, the object is considered to be a cube.
252 Returns
253 --------
254 float
255 The Lipschitz-Killing Curvature
256 """
257 from scipy.special import binom
259 if i > n or i < 0:
260 print('error in spam.excursions.elkc.lkcCuboid(): i should be between 0 and n. Here n={}, i={}'.format(n, i))
261 return 0
263 # cube
264 if isinstance(a, (int, float)):
265 # general formula for every i and n
266 lkc = binom(n, i) * a**i
267 return lkc
269 # cuboid
270 else:
271 # check if a is well defined
272 if len(a) != n:
273 print('error in spam.excursions.elkc.lkcCuboid(): Spatial dimension doesn\'t match length definition. len(a) = {} and n = {}'.format(len(a), n))
274 return 0
276 # 2D
277 if n == 2:
278 switcher = {
279 0: 1,
280 1: a[0] + a[1],
281 2: a[0] * a[1],
282 }
283 return switcher.get(i, -1)
285 # 3D
286 elif n == 3:
287 switcher = {
288 0: 1,
289 1: a[0] + a[1] + a[2],
290 2: a[0] * a[1] + a[0] * a[2] + a[2] * a[1],
291 3: a[0] * a[1] * a[2],
292 }
293 return switcher.get(i, -1)
295 else:
296 print('error in spam.excursions.elkc.lkcCuboid(): Spatial dimension n = {} not implemented.'.format(n))
297 return 0
299#
300# def phi(x, s):
301# return numpy.exp(-x**2/2.0)/(s*numpy.sqrt(2.0*pi))
302#
303#
304# def bigPhi(x, s):
305# from scipy.special import erf
306# return phi(x, s) - 0.5*x*(1.0+erf(-x/numpy.sqrt(2.0)))
309# def expectedMesuresLinearThreshold(alpha, beta, j, n, hittingSet='tail', distributionType='gauss', mu=0.0, std=1.0, correlationType='gauss', lc=1.0, nu=2.0, a=1.0):
310# """
311# Compute the Lipschitz-Killing Curvatures E(LKC)
312#
313# Parameters
314# -----------
315# kappa : float or list
316# value of the threshold
317# j : int
318# number of the functionnal
319# n : int
320# spatial dimension
321# hittingSet : string
322# the hitting set
323# distributionType : string
324# type of distribution (see gaussianMinkowskiFunctionals)
325# mu : float
326# mean value of the distribution
327# std : float
328# standard deviation of the distribution
329# correlationType : string
330# type of correlation function
331# lc : float
332# correlation length
333# nu : float
334# parameter used for correlationType='matern'
335# a : float or list of floats
336# size(s) of the object
337# if float, the object is considered as a cube
338#
339# Returns
340# --------
341# mink : same type as kappa
342# the Gaussian Minkowski functionnal
343# """
344# # HISTORY:
345# # 2016-04-08: First version of the function
346# #
347#
348# e = 0.0*alpha
349# # compute second spetral moment
350# lam2 = secondSpectralMoment(std, lc, correlationType=correlationType, nu=2.0)
351# # loop over i
352# # if prlv>5: print 'IN ELKC: n = {}, j = {}'.format( n, j )
353# # for i in range( n-j+1 ):
354# # if prlv>5: print '\t\ti = {}'.format( i )
355# # # comute LKC of the cube
356# # lkc = lkcCuboid( i+j, n, a )
357# # # compute GMF
358# # gmf = gaussianMinkowskiFunctionalsLinearThreshold( alpha, beta, a, i, hittingSet=hittingSet, lam2=lam2 )
359#
360# # # compute
361# # fla = flag( i+j , i )
362# # e = e + fla*gmf*lkc*( lam2/(2.0*pi) )**(i/2.0)
363#
364# e = 0.0
365# if n == 1:
366# from scipy.special import erf, erfc
367# cstA = alpha*a/(numpy.sqrt(2.0)*std)
368# cstB = beta/(numpy.sqrt(2.0)*std)
369# cstAB = cstA+cstB
370#
371# # intOfMo = (( numpy.sqrt(2.0/pi) * std * ( numpy.exp(-cstAB**2.0) - numpy.exp(-cstB**2) ) + (alpha*a+beta) * erf(cstAB) - beta*erf(cstB) )/(2.0*alpha) + 0.5*a) # tail
372# intOfMo = ((numpy.sqrt(2.0/pi) * std * (-numpy.exp(-cstAB**2.0) + numpy.exp(-cstB**2)) +
373# (alpha*a) * erfc(cstAB) - beta*erf(cstAB) + + beta*erf(cstB))/(2.0*alpha)) # head
374#
375# if j == 1:
376# e = intOfMo
377# elif j == 0:
378# c = numpy.sqrt(lam2) * bigPhi(alpha/(numpy.sqrt(lam2)*std), std)
379# intOfUpCrossed = c * (erf(cstAB) - erf(cstB))/(2.0*alpha)
380#
381# e = intOfMo/a + intOfUpCrossed
382# else:
383# print("ELKC not implemented for n={} and j={}-> exiting".format(n, j))
384#
385# else:
386# print("ELKC not implemented for n={} -> exiting".format(n))
387#
388# return e
389#
390#
391# def gaussianMinkowskiFunctionalsLinearThreshold(alpha, beta, a, j, hittingSet='tail', distributionType='gauss', mu=0.0, std=1.0, lam2=1.0):
392# """
393# Compute the Gaussian Minkowski functionals j
394#
395# Parameters
396# -----------
397# kappa : float or list
398# value of the threshold
399# j : int
400# number of the functionnal
401# hittingSet : string
402# the hitting set
403# distributionType : string
404# type of distribution (see gaussianMinkowskiFunctionals)
405# implemented: distributionType='gauss' and distributionType='log'
406# mu : float
407# mean value of the distribution
408# std : float
409# standard deviation of the distribution
410#
411# Returns:
412# mink : same type as kappa
413# the Gaussian Minkowski functionnal
414# """
415# # HISTORY:
416# # 2016-04-08: First version of the function
417# #
418#
419# from scipy.special import erf, erfc
420# cstA = alpha*a/(numpy.sqrt(2.0)*std)
421# cstB = beta/(numpy.sqrt(2.0)*std)
422# cstAB = cstA+cstB
423#
424# if j == 0:
425# if hittingSet == 'tail':
426# mink = ((numpy.sqrt(2.0/pi) * std * (numpy.exp(-cstAB**2.0) - numpy.exp(-cstB**2)) +
427# (alpha*a+beta) * erf(cstAB) - beta*erf(cstB))/(2.0*alpha) + 0.5*a)/a
428# elif hittingSet == 'head':
429# mink = ((numpy.sqrt(2.0/pi) * std * (-numpy.exp(-cstAB**2.0) + numpy.exp(-cstB**2)) +
430# (alpha*a) * erfc(cstAB) - beta*erf(cstAB) + + beta*erf(cstB))/(2.0*alpha))/a
431#
432# elif j == 1:
433# if hittingSet == 'tail':
434# mink = (erf(cstAB)-erf(cstB))/(2*alpha*a)
435# elif hittingSet == 'head':
436# x = alpha/(numpy.sqrt(lam2)*std)
437# mink = (erf(cstAB)-erf(cstB))/(2*alpha*a) * bigPhi(x, std)
438# else:
439# print("Error GMF not implemented", j)
440# exit()
441# return mink
444def getThresholdFromVolume(targetVolume, initialThreshold, hittingSet='tail', distributionType='gauss', mu=0.0, std=1.0, correlationType='gauss', nu=2.0, a=1.0):
445 import scipy.optimize
447 def minFunc(t):
448 vol = expectedMesures(t, 3, 3, hittingSet=hittingSet, distributionType=distributionType, mu=mu, std=std, correlationType=correlationType, nu=nu, a=a)
449 # print("t = {}, vol = {}, diff = {}".format(t, vol, vol-targetVolume))
450 return numpy.square(vol - targetVolume)
452 res = scipy.optimize.minimize(minFunc, initialThreshold, method='BFGS')
453 return res.x[0]
456def getCorrelationLengthFormArea(targetArea, initialLength, threshold, hittingSet='tail', distributionType='gauss', mu=0.0, std=1.0, correlationType='gauss', nu=2.0, a=1.0):
457 import scipy.optimize
459 def minFunc(l):
460 ar = 2.0 * expectedMesures(threshold, 2, 3, hittingSet=hittingSet, distributionType=distributionType, mu=mu, std=std, lc=l, correlationType=correlationType, nu=nu, a=a)
461 # print("l = {}, ar = {}, diff = {}".format(l, ar, ar-targetArea))
462 return numpy.square(ar - targetArea)
464 res = scipy.optimize.minimize(minFunc, initialLength, method='BFGS')
465 return res.x[0]