pynkowski.theory.gaussian

  1"""This submodule contains the definition for general Isotropic Gaussian fields, as well as several classes of 
  2Isotropic Gaussian fields defined on particular spaces.
  3"""
  4import warnings
  5import numpy as np
  6from scipy.special import comb, gamma
  7from scipy.stats import norm
  8from scipy.interpolate import interp1d
  9from scipy.integrate import quad
 10from .base_th import TheoryField
 11from .utils_th import get_μ, get_σ, LKC, lkc_ambient_dict, egoe, get_C2
 12
 13
 14class Gaussian(TheoryField):
 15    """General class for Isotropic Gaussian fields, to be used directly or as the base for specific Gaussian fields.
 16
 17    Parameters
 18    ----------
 19    dim : int
 20        Dimension of the space where the field is defined.
 21    
 22    sigma : float, optional
 23        The standard deviation of the field.
 24        Default : 1.
 25        
 26    mu : float, optional
 27        The derivative of the covariance function at the origin, times $-2$. Equal to the variance of the first derivatives of the field.
 28        Default : 1.
 29
 30    nu : float, optional
 31        The second derivative of the covariance function at the origin.
 32        Default : 1.
 33        
 34    lkc_ambient : np.array or None, optional
 35        The values for the Lipschitz–Killing Curvatures of the ambient space. If `None`, it is assumed that the volume is 1 and the rest is 0. This is exact for many spaces like Euclidean spaces or the sphere, and exact to leading order in μ for the rest.
 36        Default : None
 37        
 38    Attributes
 39    ----------
 40    dim : int
 41        Dimension of the space where the field is defined.
 42    
 43    name : str
 44        Name of the field, `"Isotropic Gaussian"` by default.
 45    
 46    sigma : float
 47        The standard deviation of the field.
 48        
 49    mu : float
 50        The derivative of the covariance function at the origin, times $-2$. Equal to the variance of the first derivatives of the field.
 51
 52    nu : float
 53        The second derivative of the covariance function at the origin. Equal to $12$ times the variance of the second derivatives of the field ($4$ times for the cross derivative).
 54        
 55    lkc_ambient : np.array or None
 56        The values for the Lipschitz–Killing Curvatures of the ambient space.
 57        
 58    """   
 59    def __init__(self, dim, sigma=1., mu=1., nu=1., lkc_ambient=None):
 60        super().__init__(dim, name='Isotropic Gaussian', sigma=sigma, mu=mu, nu=nu, lkc_ambient=lkc_ambient)
 61        self._warn_non_euclidean = True
 62        
 63    def LKC(self, j, us):
 64        """Compute the expected values of the Lipschitz–Killing Curvatures of the excursion sets at thresholds `us`, $\mathbb{L}_j(A_u(f))$.
 65
 66        Parameters
 67        ----------
 68        j : int
 69            Index of the LKC to compute, `0 < j < dim`.
 70            
 71        us : np.array
 72            The thresholds considered for the computation.
 73
 74        Returns
 75        -------
 76        lkc : np.array
 77            Expected value of LKC evaluated at the thresholds.
 78
 79        """
 80        us /= self.sigma
 81        return LKC(j, us, self.mu, dim=self.dim, lkc_ambient=self.lkc_ambient)/self.lkc_ambient[-1]
 82    
 83    def MF(self, j, us):
 84        """Compute the expected values of the Minkowski Functionals of the excursion sets at thresholds `us`, $V_j(A_u(f))$.
 85
 86        Parameters
 87        ----------
 88        j : int
 89            Index of the MF to compute, `0 < j < dim`.
 90            
 91        us : np.array
 92            The thresholds considered for the computation.
 93
 94        Returns
 95        -------
 96        mf : np.array
 97            Expected value of MFs evaluated at the thresholds.
 98
 99        """
100        return 1./comb(self.dim,j)*self.LKC(self.dim-j, us)
101    
102    def V0(self, us):
103        """Compute the expected values of the first Minkowski Functionals of the excursion sets at thresholds `us`, $V_0(A_u(f))$.
104
105        Parameters
106        ----------
107        us : np.array
108            The thresholds considered for the computation.
109
110        Returns
111        -------
112        v0 : np.array
113            Expected value of $V_0$ evaluated at the thresholds.
114
115        """
116        us /= self.sigma
117        return 1. - norm.cdf(us)
118    
119    def V1(self, us):
120        """Compute the expected values of the first Minkowski Functionals of the excursion sets at thresholds `us`, $V_1(A_u(f))$.
121
122        Parameters
123        ----------
124        us : np.array
125            The thresholds considered for the computation.
126
127        Returns
128        -------
129        v1 : np.array
130            Expected value of $V_1$ evaluated at the thresholds.
131
132        """
133        us /= self.sigma
134        return self.MF(1, us)
135    
136    def V2(self, us):
137        """Compute the expected values of the first Minkowski Functionals of the excursion sets at thresholds `us`, $V_2(A_u(f))$. Only for fields with `dim ≥ 2`.
138
139        Parameters
140        ----------
141        us : np.array
142            The thresholds considered for the computation.
143
144        Returns
145        -------
146        v2 : np.array
147            Expected value of $V_2$ evaluated at the thresholds.
148
149        """
150        assert self.dim>=2, 'V2 is defined only for fields with dim≥2'
151        us /= self.sigma
152        return self.MF(2, us)
153    
154    def V3(self, us):
155        """Compute the expected values of the first Minkowski Functionals of the excursion sets at thresholds `us`, $V_3(A_u(f))$. Only for fields with `dim ≥ 3`.
156
157        Parameters
158        ----------
159        us : np.array
160            The thresholds considered for the computation.
161
162        Returns
163        -------
164        v3 : np.array
165            Expected value of $V_3$ evaluated at the thresholds.
166
167        """
168        assert self.dim>=3, 'V3 is defined only for fields with dim≥3'
169        us /= self.sigma
170        return self.MF(3, us)
171    
172    # def V4(self, us):
173    #     assert self.dim>=4, 'V4 is defined only for fields with dim≥4'
174    #     us /= self.sigma
175    #     return self.MF(4, us)
176
177    def maxima_total(self):
178        """Compute the expected values of local maxima of the field.
179
180        Returns
181        -------
182        number_total : float
183            Expected value of the number of local maxima.
184
185        """
186        if self._warn_non_euclidean:
187            warnings.warn('If the ambient space is not euclidean, this function is an approximation. You can use `SphericalGaussian` or `EuclideanGaussian` instead.')
188        rho1 = -0.5*self.mu
189        return (2./np.pi)**((self.dim+1.)/2.) * gamma((self.dim+1.)/2.) * (-self.nu/rho1)**(self.dim/2.) * egoe(self.dim, 1., 0.) * self.lkc_ambient[-1]
190
191    def maxima_dist(self, us):
192        """Compute the expected distribution of local maxima of the field, as a function of threshold.
193
194        Parameters
195        ----------
196        us : np.array
197            The thresholds considered for the computation.
198
199        Returns
200        -------
201        density : np.array
202            Density distribution of the local maxima.
203
204        """
205        rho1 = -0.5*self.mu
206        κ = -rho1 / np.sqrt(self.nu)
207        assert κ <= (self.dim + 2.)/self.dim, '`0.5*mu/sqrt(nu)` must be $≤ (dim+2)/dim$'
208        if self._warn_non_euclidean:
209            warnings.warn('If the ambient space is not euclidean, this function is an approximation. You can use `SphericalGaussian` or `EuclideanGaussian` instead.')
210        return np.real(np.sqrt(1./(1.-κ**2.+ 0.j)) * norm.pdf(us) * egoe(self.dim, 1./(1.-κ**2.), κ*us/np.sqrt(2.)) / egoe(self.dim, 1., 0.))
211
212    def minima_total(self):
213        """Compute the expected values of local minima of the field.
214
215        Returns
216        -------
217        number_total : float
218            Expected value of the number of local minima.
219
220        """
221        return self.maxima_total()
222
223    def minima_dist(self, us):
224        """Compute the expected distribution of local minima of the field, as a function of threshold.
225
226        Parameters
227        ----------
228        us : np.array
229            The thresholds considered for the computation.
230
231        Returns
232        -------
233        density : np.array
234            Density distribution of the local minima.
235
236        """
237        return self.maxima_dist(-us)
238    
239    
240class SphericalGaussian(Gaussian):
241    """Class for Spherical Isotropic Gaussian fields.
242
243    Parameters
244    ----------
245    cls : np.array
246        Angular power spectrum of the field.
247    
248    normalise : bool, optional
249        If `True`, normalise the field to unit variance.
250        Default : True
251
252    fsky : float, optional
253        Fraction of the sky covered by the field, `0<fsky<=1`.
254        Default : 1.
255    
256    Attributes
257    ----------
258    cls : np.array
259        Angular Power Spectrum of the field.
260    
261    fsky : float
262        Fraction of the sky covered by the field.
263
264    dim : int
265        Dimension of the space where the field is defined, in this case this is 2.
266    
267    name : str
268        Name of the field, `"Spherical Isotropic Gaussian"` by default.
269    
270    sigma : float
271        The standard deviation of the field.
272
273    mu : float
274        The derivative of the covariance function at the origin, times $-2$. Equal to the variance of the first derivatives of the field.
275
276    nu : float
277        The second derivative of the covariance function at the origin. 
278        
279    C2 : float
280        The second derivative of the angular covariance function at 1. 
281        
282    lkc_ambient : np.array or None
283        The values for the Lipschitz–Killing Curvatures of the ambient space.
284        
285    """   
286    def __init__(self, cls, normalise=True, fsky=1.):
287        if normalise:
288            cls /= get_σ(cls)
289            self.sigma = 1.
290        else:
291            self.sigma = np.sqrt(get_σ(cls))
292        self.cls = cls
293        self.fsky = fsky
294        self.mu = get_μ(cls)
295        self.C2 = get_C2(cls)
296        self.nu = self.C2/4. - self.mu/24.
297        super().__init__(2, sigma=self.sigma, mu=self.mu, nu=self.nu, lkc_ambient=lkc_ambient_dict["sphere"]*self.fsky)
298        self.name = 'Spherical Isotropic Gaussian'
299
300    def maxima_total(self):
301        """Compute the expected values of local maxima of the field.
302
303        Returns
304        -------
305        number_total : float
306            Expected value of the number of local maxima.
307
308        """
309        k1 = self.mu/self.C2
310        return np.sqrt(2./(1.+k1)) / np.pi**((self.dim+1.)/2.) * k1**(-self.dim/2) * gamma((self.dim+1.)/2.) * egoe(self.dim,1./(1.+k1),0.) * self.lkc_ambient[-1]
311
312    def maxima_dist(self, us):
313        """Compute the expected distribution of local maxima of the field, as a function of threshold.
314
315        Parameters
316        ----------
317        us : np.array
318            The thresholds considered for the computation.
319
320        Returns
321        -------
322        density : np.array
323            Density distribution of the local maxima.
324
325        """
326        k1 = self.mu/self.C2
327        k2 = self.mu**2./self.C2
328        return np.real(np.sqrt((1.+k1)/(1.+k1-k2) +0.j) * norm.pdf(us) * egoe(self.dim, 1./(1.+k1-k2), np.sqrt(k2/2.)*us) / egoe(self.dim, 1./(1.+k1), 0.))
329
330# class EuclideanGaussian(Gaussian):
331#     """Class for Isotropic Gaussian fields defined in an Euclidean space (such as a 2D plane or a 3D cube). 
332
333#     Parameters
334#     ----------
335#     dim : int
336#         Dimension of the space where the field is defined.
337
338#     Pk : np.array or function
339#         Power spectrum of the field. It can be an array with the values of the power spectrum at wavenumbers `k`, or a function Pk(k) which returns the power spectrum at wavenumbers `k`. The latter is more accurate.
340
341#     k : np.array
342#         If `Pk` is an array, wavenumbers corresponding to the power spectrum. If `Pk` is a function, `k` is the array `[kmin, kmax]`.
343    
344#     normalise : bool, optional
345#         If `True`, normalise the field to unit variance.
346#         Default : True
347
348#     volume : float, optional
349#         Hausdorff measure of the space where the field is defined (area if `dim=2`, volume if `dim=3`, and so on).
350#         Default : 1.
351    
352#     Attributes
353#     ----------
354#     dim : int
355#         Dimension of the space where the field is defined.
356
357#     Pk : function
358#         Power spectrum of the field, to be evaluated at wavenumbers `k`.
359
360#     klim : np.array
361#         Minimum and maximum wavenumbers to be considered for the power spectrum `Pk`.
362    
363#     name : str
364#         Name of the field, `"Euclidean Isotropic Gaussian"` by default.
365    
366#     sigma : float
367#         The standard deviation of the field.
368
369#     mu : float
370#         The derivative of the covariance function at the origin, times $-2$. Equal to the variance of the first derivatives of the field.
371
372#     nu : float
373#         The second derivative of the covariance function at the origin. Equal to $12$ times the variance of the second derivatives of the field ($4$ times for the cross derivative).
374        
375#     lkc_ambient : np.array
376#         The values for the Lipschitz–Killing Curvatures of the ambient space.
377        
378#     """   
379#     def __init__(self, dim, Pk, k, normalise=True, volume=1.):
380#         if type(Pk) is np.ndarray:
381#             self.Pk = interp1d(k, Pk)
382#             self.klim = np.array([k.min(), k.max()])
383#         else:
384#             assert callable(Pk), 'The power spectrum must be a function or an array'
385#             assert k.shape == (2,), 'The power spectrum must be an array with two elements: `kmin, kmax`'
386#             self.klim = k
387#             self.Pk = Pk
388        
389#         if dim != 3:
390#             warnings.warn('The formulae for the computation of sigma, mu, and nu are valid for `dim=3`. Please verify and possibly redefine these quantities manually.')
391#             # Are these formulae valid also for other dimensions?
392#         k2Pk = lambda k: self.Pk(k)*k**2.
393#         sigma = np.sqrt(quad(k2Pk, self.klim[0], self.klim[1])[0] / (2.*np.pi**2.)) #* volume
394
395#         k4Pk = lambda k: self.Pk(k)*k**4.
396#         mu = quad(k4Pk, self.klim[0], self.klim[1])[0] / (6.*np.pi**2.) #* volume
397
398#         k6Pk = lambda k: self.Pk(k)*k**6.
399#         nu = quad(k6Pk, self.klim[0], self.klim[1])[0] / (120.*np.pi**2.) #* volume
400
401#         if normalise:
402#             mu = mu/sigma**2.
403#             nu = nu/sigma**2.
404#             self.Pk = lambda k: self.Pk(k)/sigma**2.
405#             sigma = 1.
406        
407#         lkc_ambient = np.zeros(dim+1)
408#         lkc_ambient[-1] = volume
409#         super().__init__(dim=dim, sigma=sigma, mu=mu, nu=nu, lkc_ambient=lkc_ambient)
410#         self.name = 'Euclidean Isotropic Gaussian'
411#         self._warn_non_euclidean = False
412
413#     def maxima_total(self):
414#         """Compute the expected values of local maxima of the field.
415
416#         Returns
417#         -------
418#         number_total : float
419#             Expected value of the number of local maxima.
420
421#         """
422#         rho1 = -0.5*self.mu
423#         return (2./np.pi)**((self.dim+1.)/2.) * gamma((self.dim+1.)/2.) * (-self.nu/rho1)**(self.dim/2.) * egoe(self.dim, 1., 0.) * self.lkc_ambient[-1]
424
425#     def maxima_dist(self, us):
426#         """Compute the expected distribution of local maxima of the field, as a function of threshold.
427
428#         Parameters
429#         ----------
430#         us : np.array
431#             The thresholds considered for the computation.
432
433#         Returns
434#         -------
435#         density : np.array
436#             Density distribution of the local maxima.
437
438#         """
439#         rho1 = -0.5*self.mu
440#         κ = -rho1 / np.sqrt(self.nu)
441#         assert κ <= (self.dim + 2.)/self.dim, '`0.5*mu/sqrt(nu)` must be $≤ (dim+2)/dim$'
442#         return np.real(np.sqrt(1./(1.-κ**2.+ 0.j))) * norm.pdf(us) * egoe(self.dim, 1./(1.-κ**2.), κ*us/np.sqrt(2.)) / egoe(self.dim, 1., 0.)
443
444
445    
446
447__all__ = ["SphericalGaussian", "Gaussian"] #, "EuclideanGaussian", 
448
449__docformat__ = "numpy"
450
451