pynkowski.theory.utils_th

  1'''This submodule contains utilities for the theory module.'''
  2import warnings
  3import numpy as np
  4from scipy.special import eval_hermitenorm, gamma, comb, factorial, gammainc
  5from scipy.stats import norm, multivariate_normal
  6from .base_th import _prepare_lkc
  7
  8
  9def get_σ(cls):
 10    """Compute the variance of a field with an angular power spectrum 'cls'.
 11
 12    Parameters
 13    ----------
 14    cls : np.array
 15        The angular power spectrum of the Gaussian field.
 16    
 17    Returns
 18    -------
 19    σ² : float
 20        The variance of the field.
 21    
 22    """
 23    cls = np.array(cls, dtype=float)
 24    ℓ = np.arange(cls.shape[0])
 25    return np.sum(cls * (2.*ℓ+1.) / (4.*np.pi))
 26
 27
 28def get_μ(cls):
 29    """Compute the first derivative of the covariance function at the origin for a Gaussian field 
 30    defined on the sphere with angular power spectrum 'cls', which are normalised to unit variance.
 31
 32    Parameters
 33    ----------
 34    cls : np.array
 35        The angular power spectrum of the Gaussian field.
 36    
 37    Returns
 38    -------
 39    μ : float
 40        The first derivative of the covariance function of a field at the origin.
 41    
 42    """
 43    cls = np.array(cls, dtype=float)
 44    ℓ = np.arange(cls.shape[0])
 45    cls /= np.sum(cls * (2.*ℓ+1.) / (4.*np.pi))
 46    μ = np.sum(cls * (2.*ℓ+1.) * ℓ*(ℓ+1.) / (8.*np.pi))
 47    return μ
 48
 49def get_C2(cls):
 50    """Compute the second derivative of the covariance function at the origin for a Gaussian field 
 51    defined on the sphere with angular power spectrum 'cls', which are normalised to unit variance.
 52
 53    Parameters
 54    ----------
 55    cls : np.array
 56        The angular power spectrum of the Gaussian field.
 57    
 58    Returns
 59    -------
 60    C2 : float
 61        The second derivative of the covariance function of a field at the origin.
 62    
 63    """
 64    cls = np.array(cls, dtype=float)
 65    ℓ = np.arange(cls.shape[0])
 66    cls /= np.sum(cls * (2.*ℓ+1.) / (4.*np.pi))
 67    C2 = np.sum(cls * (2.*ℓ+1.) * (ℓ-1.)*ℓ*(ℓ+1.)*(ℓ+2.) / (32.*np.pi))
 68    return C2
 69
 70def flag(kj, k):
 71    """Compute flag coefficients.
 72
 73    Parameters
 74    ----------
 75    ki : int
 76        The upper index of the flag coefficient.
 77        
 78    k : int
 79        The lower index of the flag coefficient.
 80    
 81    Returns
 82    -------
 83    flag : float
 84        The value of the float coefficient.
 85    
 86    """
 87    assert np.all(kj>=k), 'The first argument must be larger than the second.'
 88    return comb(kj, k) * omega(kj) / (omega(k) * omega(kj-k))
 89    
 90def omega(j):
 91    """Compute the volume of the $j$—dimensional ball, $\omega_j$.
 92
 93    Parameters
 94    ----------
 95    j : int
 96        The dimension of the ball.
 97    
 98    Returns
 99    -------
100    omega : float
101        The volume of the $j$—dimensional ball.
102    
103    """
104    return np.pi**(j/2.) / gamma(j/2. +1.)
105
106
107def rho(k, us):
108    """Compute the density functions $\rho_k(us)$ for Gaussian fields needed for the Gaussian Kinematic Formula.
109
110    Parameters
111    ----------
112    k : int
113        The order of the density function.
114        
115    us : np.array
116        The thresholds where the density function is evaluated.
117    
118    Returns
119    -------
120    rho : np.array
121        The density function evaluated at the thresholds.
122    
123    """
124    if k==0:
125        return 1. - norm.cdf(us)
126    else:
127        return 1. / (2. * np.pi)**(k/2.) * norm.pdf(us) * eval_hermitenorm(k-1, us)
128    
129def rho_Chi2(k, dof, us):
130    """Compute the density functions $\rho_{k, \Chi^2}(us)$ for $\Chi^2_k$ fields (with $k$ degrees of freedom), needed for the Gaussian Kinematic Formula.
131
132    Parameters
133    ----------
134    k : int
135        The order of the density function.
136        
137    dof : int
138        The degrees of freedom of the $\Chi^2$.
139        
140    us : np.array
141        The thresholds where the density function is evaluated.
142    
143    Returns
144    -------
145    rho : np.array
146        The density function evaluated at the thresholds.
147    
148    """
149    if k==0:
150        return 1.- (us >= 0.) * gammainc(dof/2., us/2.)
151    else:
152        factor = (us >= 0.) * us**((dof-k)/2.) * np.exp(-us/2.) / ((2.*np.pi)**(k/2.) * gamma(dof/2.) * 2.**((dof-2)/2.))
153        summ = np.zeros_like(us)
154        for l in np.arange(0, ((k-1)//2)+1):
155            for m in np.arange(0,k-1-2*l+1):
156                summ += (dof>= k - m - 2*l) * comb(dof-1, k-1-m-2*l) * (-1)**(k-1+m+l) * factorial(k-1) / (factorial(m) * factorial(l) * 2**l ) * us**(m+l)
157        return factor*summ
158    
159def LKC(j, us, mu, dim=None, lkc_ambient=None):
160    """Compute the expected value of the Lipschitz–Killing Curvatures (LKC) of the excursion set for Gaussian Isotropic fields.
161
162    Parameters
163    ----------
164    j : int
165        The index of the LKC.
166        
167    us : np.array
168        The thresholds where the LKC is evaluated.
169
170    mu : float
171        The value of the derivative of the covariance function at the origin for the field.
172        
173    dim : int, optional
174        The dimension of the ambient manifold.
175        
176    lkc_ambient : np.array, optional
177        An array of the Lipschitz–Killing Curvatures of the ambient manifold. Its lenght must be `dim+1` if `dim` is also given.
178        
179    Returns
180    ----------
181    LKC : np.array
182        The expected value of the Lipschitz–Killing Curvatures at the thresholds.
183    """
184    dim, lkc_ambient = _prepare_lkc(dim, lkc_ambient)
185    result = np.zeros_like(us)
186    for k in np.arange(0,dim-j+1):
187        result += flag(k+j, k) * rho(k, us) * lkc_ambient[k+j] * mu**(k/2.)
188    return result
189    
190def LKC_Chi2(j, us, mu, dof=2, dim=None, lkc_ambient=None):
191    """Compute the expected value of the Lipschitz–Killing Curvatures (LKC) of the excursion set  for $\Chi^2_2$ fields.
192
193    Parameters
194    ----------
195    j : int
196        The index of the LKC.
197        
198    us : np.array
199        The thresholds where the LKC is evaluated.
200
201    mu : float
202        The value of the derivative of the covariance function at the origin for the original Gaussian fields.
203
204    dof : int, optional
205        Number of degrees of freedom of the $\Chi^2$ field (i.e., number of squared Gaussian distributions that are summed). 
206        Default: 2.
207        
208    dim : int, optional
209        The dimension of the ambient manifold.
210        
211    lkc_ambient : np.array, optional
212        An array of the Lipschitz–Killing Curvatures of the ambient manifold. Its lenght must be `dim+1` if `dim` is also given.
213        
214    Returns
215    ----------
216    LKC : np.array
217        The expected value of the Lipschitz–Killing Curvatures at the thresholds.
218    """
219    dim, lkc_ambient = _prepare_lkc(dim, lkc_ambient)
220    result = np.zeros_like(us)
221    for k in np.arange(0,dim-j+1):
222        result += flag(k+j, k) * rho_Chi2(k, dof, us) * lkc_ambient[k+j] * mu**(k/2.)
223    return result
224
225
226def egoe(dim, a, b):
227    """Compute the expected value of the Gaussian Orthogonal Ensamble, needed for the maxima distribution of Gaussian Isotropic fields.
228    Details can be found in Cheng and Schwartzman, 2015. This formula is only proven for `a>0`, although it has been hypothetised to work for `a<0` under some conditions.
229    
230    Parameters
231    ----------
232    sim : int
233        The dimension of the space (1 ≤ dim ≤ 3)
234        
235    a : float
236        Parameter `a` of the formula. 
237
238    b : float
239        Parameter `b` of the formula.
240
241    Returns
242    ----------
243    egoe : float
244        The expected value needed for the computation of the maxima distribution.
245    """
246    assert dim in [1, 2, 3], '`dim` must be between 1 and 3'
247    if a < 0:
248        a = a + 0.j
249        b = b + 0.j
250        warnings.warn('This formula is only proven for `a>=0`, although it has been hypothetised to work for `a<0` under some conditions. If the maxima distribution looks wrong, it could be because of this.')
251    # assert a > 0, '`a` must be positive'
252    if dim == 1:
253        return np.sqrt(4.*a+2.) * np.exp(-a*b**2./(2.*a +1.)) / (4.*a) + b*np.sqrt(np.pi/(2.*a)) * norm.cdf(b * np.sqrt(2.*a / (2.*a + 1)) )
254    if dim == 2:
255        return ( (1./a + 2.*b**2. - 1.) / np.sqrt(2.*a) * norm.cdf(b * np.sqrt(2.*a / (a + 1.))) +
256            np.sqrt((a + 1.) / (2.*np.pi)) * np.exp(-a*b**2./(a +1.)) * b / a +
257            np.sqrt(2./(2.*a+1.)) * np.exp(-a*b**2./(2.*a +1.)) * norm.cdf(a * b * np.sqrt(2. / ((2.*a + 1.)*(a + 1.))) ) )
258    if dim == 3:
259        Sigma1 = np.array([[1.5, -0.5], [-0.5, (1.+a)/(2.*a)]])
260        Sigma2 = np.array([[1.5, -1.], [-1., (1.+2.*a)/(2.*a)]])
261        term1 = (((24.*a**3. + 12.*a**2. + 6.*a + 1.) * b**2. / (2.*a * (2.*a+1)**2. ) +
262                (6.*a**2 + 3.*a +2.) / (4.*a**2. * (2.*a+1)) + 3./2.) / np.sqrt(2. * (2.*a +1.)) * 
263                np.exp(-a*b**2./(2.*a +1.)) * norm.cdf(2. * a * b * np.sqrt(2. / ((2.*a + 1.) * (2.*a +3.)) ) ) )
264        term2 = (((a+1.) * b**2. / (2.*a) + (1.-a) / (2.*a**2.) - 1.) / np.sqrt(2. * (a +1.)) *
265                np.exp(-a*b**2./(a +1.)) * norm.cdf(a * b * np.sqrt(2. / ((a + 1.) * (2.*a +3.)) ))  )
266        term3 = (6.*a + 1. + (28.*a**2. + 12.*a + 3.) / (2.*a * (2.*a+1)) * 
267                b / (2. * np.sqrt(2.*np.pi) * (2.*a+1.) * np.sqrt(2.*a +3.)) * np.exp(-3.*a*b**2./(2.*a +3.))  )
268        term4 = (b**2. + 3. * (1.-a) / (2.*a)) * np.sqrt(np.pi/2.) * b/a * (multivariate_normal.cdf([0,b], cov=Sigma1) + multivariate_normal.cdf([0,b], cov=Sigma2))
269        total = term1 + term2 + term3 + term4
270        return total
271    
272
273lkc_ambient_dict = {"2D":np.array([0., 0., 1.]),
274                    "3D":np.array([0., 0., 0., 1.]),
275                    "sphere":np.array([0., 0., 4.*np.pi]),  # / (4.*np.pi),
276                    "SO3":np.array([0., 6.*np.pi, 0., 4.*np.pi**2])}  # / (4.*np.pi**2)}
277"""Dictionary with the characterization of different spaces through their Lipschitz–Killing Curvatures"""
278
279
280__all__ = ["get_μ", "flag", "omega", "rho", "rho_Chi2", "LKC", "LKC_Chi2", "lkc_ambient_dict"]
281
282__docformat__ = "numpy"
283
284