pynkowski.stats.minkowski

  1'''This submodule contains the functions to compute the Minkowski Functionals of a field.'''
  2import numpy as np
  3from scipy.special import comb, gamma
  4from .utils_st import subsample_us, define_ubins
  5from ..data.base_da import DataField
  6from ..theory.base_th import TheoryField
  7
  8try:
  9    from tqdm.auto import tqdm
 10except:
 11    tqdm = lambda x: x
 12    print('tqdm not loaded')
 13
 14def _MF_prefactor(d,j):
 15    """Compute the prefactor in the definition of Minkowski Functionals. This factor multiplies the integral of the curvatures.
 16
 17    Parameters
 18    ----------
 19    d : int
 20        The dimension of the ambient space (e.g., 2 for the sphere).
 21        
 22    j : int
 23        The index of the Minkowski functional (`j>=1`).
 24    
 25    Returns
 26    -------
 27    prefactor : float
 28        The prefactor in the definition of MFs.
 29    
 30    """
 31    return 1./ (comb(d,j) * j * np.pi**(j/2.) / gamma(j/2. +1.))
 32    
 33
 34def V0(field, us, edges=False, verbose=True):
 35    """Compute the first Minkowski Functional, $V_0$, normalized by the volume of the space.
 36
 37    Parameters
 38    ----------
 39    field : DataField or TheoryField
 40        Field on which to compute $V_0$, which can be a theoretical field or a data field.
 41        
 42    us : np.array
 43        The thresholds where $V_0$ is computed.
 44        
 45    edges : bool, optional
 46        If False (default), the given `us` is assumed to be an array of uniformly distributed thresholds, which are taken as the central values of the bins.
 47        If True, input `us` is assumed to be a monotonically increasing array of bin edges, including the rightmost edge, allowing for non-uniform distributions of thresholds. 
 48        In the latter case, the effective thresholds are the central value of the given bins.
 49        
 50    verbose : bool, optional
 51        If True (default), progress bars are shown for the computations on data.
 52    
 53    Returns
 54    -------
 55    V0 : np.array
 56        The values of the first Minkowski Functional at the given thresholds.
 57    
 58    """
 59    us = np.atleast_1d(us)
 60    us, dus = define_ubins(us, edges)
 61    
 62    if isinstance(field, TheoryField):
 63        us = subsample_us(us, dus)
 64        try:
 65            return np.nanmean(field.V0(us), axis=1)
 66        except AttributeError:
 67            raise NotImplementedError(f"The theoretical expectation of V0 for {field.name} fields is not implemented. If you know an expression, please get in touch.")
 68        
 69    if isinstance(field, DataField):
 70        try:
 71            return field._V0(us, dus, verbose=verbose)
 72        except AttributeError:                
 73            stat = np.zeros_like(us)
 74            for ii in tqdm(np.arange(us.shape[0]), disable=not verbose):
 75                stat[ii] = np.mean((field.field > us[ii])[field.mask])
 76            return stat
 77        
 78    else:
 79        raise TypeError(f"The field must be either TheoryField or DataField (or a subclass).")
 80    
 81
 82def V1(field, us, edges=False, verbose=True):
 83    """Compute the second Minkowski Functional, $V_1$, normalized by the volume of the space.
 84
 85    Parameters
 86    ----------
 87    field : DataField or TheoryField
 88        Field on which to compute $V_1$, which can be a theoretical field or a data field.
 89        
 90    us : np.array
 91        The thresholds where $V_1$ is computed.
 92        
 93    edges : bool, optional
 94        If False (default), the given `us` is assumed to be an array of uniformly distributed thresholds, which are taken as the central values of the bins.
 95        If True, input `us` is assumed to be a monotonically increasing array of bin edges, including the rightmost edge, allowing for non-uniform distributions of thresholds. 
 96        In the latter case, the effective thresholds are the central value of the given bins.
 97        
 98    verbose : bool, optional
 99        If True (default), progress bars are shown for the computations on data.
100    
101    Returns
102    -------
103    V1 : np.array
104        The values of the second Minkowski Functional at the given thresholds.
105    
106    """
107    us = np.atleast_1d(us)
108    us, dus = define_ubins(us, edges)
109    
110    if isinstance(field, TheoryField):
111        us = subsample_us(us, dus)
112        try:
113            return np.nanmean(field.V1(us), axis=1)
114        except AttributeError:
115            raise NotImplementedError(f"The theoretical expectation of V1 for {field.name} fields is not implemented. If you know an expression, please get in touch.") from None
116        
117    elif isinstance(field, DataField):
118        try:
119            return field._V1(us, dus, verbose=verbose)
120        except AttributeError:
121            if field.first_der is None:
122                field.get_first_der()
123                
124            stat = np.zeros_like(us)
125            grad_modulus = np.sqrt(np.sum((field.first_der)**2., axis=0))
126            for ii in tqdm(np.arange(us.shape[0]), disable=not verbose):
127                this_mask = (us[ii] + dus[ii]/2. > field.field) & (us[ii] - dus[ii]/2. <= field.field)
128                stat[ii] = np.mean((grad_modulus*this_mask/dus[ii])[field.mask])
129            return _MF_prefactor(field.dim, 1) * stat
130        
131    else:
132        raise TypeError(f"The field must be either TheoryField or DataField (or a subclass).")
133        
134        
135  
136def general_curvature(field, order):
137    """Computes the symmetrised polynomial of the principal curvatures of order `order`, multiplied by the modulus of the gradient, at every pixel.
138
139    Parameters
140    ----------
141    field : DataField
142        Field on which to compute the general curvatures.
143        
144    order : int
145        Order of the polynomial of the curvatures.
146    
147    Returns
148    -------
149    curvature : np.array
150        The value of the polynomial at each pixel. 
151        
152        For example:
153        `field.dim = 2`, `order=1` : geodesic curvature $\kappa$
154        `field.dim = 3`, `order=1` : mean curvature (H = ($\kappa_1$ + $\kappa_2$)/2 )
155        `field.dim = 3`, `order=2` : Gaussian curvature (K = $\kappa_1 \cdot \kappa_2$)
156    
157    """
158    assert isinstance(field, DataField), "The field must be a DataField (or subclass)"
159    assert field.dim > 1, "The field must be at least 2D"
160    if order >= field.dim:
161        raise ValueError("The order of the polynomial has to be less than the dimension of the field.")
162        
163    if order == 0:
164        return np.sqrt(np.sum(field.first_der**2., axis=0))
165    
166    if field.dim == 2:
167        if order == 1: # Geodesic curvature κ  (times |∇f|)
168            num = 2.*field.first_der[0]*field.first_der[1]*field.second_der[2] - field.first_der[0]**2.*field.second_der[1] - field.first_der[1]**2.*field.second_der[0]
169            den = field.first_der[0]**2. + field.first_der[1]**2.
170            return num / den
171            
172    if field.dim == 3:
173        if order == 1: # Mean curvature (times 2), 2H = κ_1 + κ_2   (times |∇f|)
174            mod_grad = np.sqrt(np.sum(field.first_der**2., axis=0))
175            hessian = np.array([[field.second_der[0], field.second_der[3], field.second_der[4]],
176                                [field.second_der[3], field.second_der[1], field.second_der[5]],
177                                [field.second_der[4], field.second_der[5], field.second_der[2]]])
178            norm_grad = (field.first_der / mod_grad)
179            mean_curvature = np.einsum('j...,jk...,k...->...', norm_grad, hessian, norm_grad) - np.trace(hessian, axis1=0, axis2=1)
180            return mean_curvature
181            
182        if order == 2: # Gaussian curvature, K = κ_1 · κ_2   (times |∇f|)
183            mod_grad = np.sqrt(np.sum(field.first_der**2., axis=0))
184            norm_grad = (field.first_der / mod_grad)
185            norm_second_der = (field.second_der / mod_grad)
186            extended_hessian_det = np.linalg.det(np.array([[norm_second_der[0], norm_second_der[3], norm_second_der[4], norm_grad[0]],
187                                                           [norm_second_der[3], norm_second_der[1], norm_second_der[5], norm_grad[1]],
188                                                           [norm_second_der[4], norm_second_der[5], norm_second_der[2], norm_grad[2]],
189                                                           [norm_grad[0], norm_grad[1], norm_grad[2], np.zeros_like(norm_grad[0])]]).T).T
190            return -extended_hessian_det*mod_grad
191
192        
193    if field.dim > 3:
194        raise NotImplementedError("The computation of principal curvatures in not implemented for spaces with 4 or more dimensions. If you need this, please get in touch.")
195        # I am not sure if there is a easy expression for the principal curvatures of implicit surfaces in higher dimensions. Typically, one would need to create
196        # the second fundamental form of the surface with the orthonormal bais of the tangent space; its eigenvalues are the principal curvatures at that point.
197        
198
199    
200def V2(field, us, edges=False, verbose=True):
201    """Compute the third Minkowski Functional, $V_2$, normalized by the volume of the space.
202
203    Parameters
204    ----------
205    field : DataField or TheoryField
206        Field on which to compute $V_2$, which can be a theoretical field or a data field. Its dimension must be at least 2.
207        
208    us : np.array
209        The thresholds where $V_2$ is computed.
210        
211    edges : bool, optional
212        If False (default), the given `us` is assumed to be an array of uniformly distributed thresholds, which are taken as the central values of the bins.
213        If True, input `us` is assumed to be a monotonically increasing array of bin edges, including the rightmost edge, allowing for non-uniform distributions of thresholds. 
214        In the latter case, the effective thresholds are the central value of the given bins.
215        
216    verbose : bool, optional
217        If True (default), progress bars are shown for the computations on data.
218    
219    Returns
220    -------
221    V2 : np.array
222        The values of the second Minkowski Functional at the given thresholds.
223    
224    """
225    if field.dim<2:
226        raise ValueError(f"V2 is defined for fields with at least 2 dimensions, but this field is {field.dim}D.")
227        
228    us = np.atleast_1d(us)
229    us, dus = define_ubins(us, edges)
230    
231    if isinstance(field, TheoryField):
232        us = subsample_us(us, dus)
233        try:
234            return np.nanmean(field.V2(us), axis=1)
235        except AttributeError:
236            raise NotImplementedError(f"The theoretical expectation of V2 for {field.name} fields is not implemented. If you know an expression, please get in touch.") from None
237        
238    elif isinstance(field, DataField):
239        try:
240            return field._V2(us, dus, verbose=verbose)
241        except AttributeError:
242            if field.first_der is None:
243                field.get_first_der()
244            if field.second_der is None:
245                field.get_second_der()
246                
247            stat = np.zeros_like(us)
248            curv = general_curvature(field, 1)
249            for ii in tqdm(np.arange(us.shape[0]), disable=not verbose):
250                this_mask = (us[ii] + dus[ii]/2. > field.field) & (us[ii] - dus[ii]/2. <= field.field)
251                stat[ii] = np.mean((curv*this_mask/dus[ii])[field.mask])
252            return _MF_prefactor(field.dim, 2) * stat
253        
254    else:
255        raise TypeError(f"The field must be either TheoryField or DataField (or a subclass).")
256        
257
258def V3(field, us, edges=False, verbose=True):
259    """Compute the third Minkowski Functional, $V_3$, normalized by the volume of the space.
260
261    Parameters
262    ----------
263    field : DataField or TheoryField
264        Field on which to compute $V_3$, which can be a theoretical field or a data field. Its dimension must be at least 3.
265        
266    us : np.array
267        The thresholds where $V_3$ is computed.
268        
269    edges : bool, optional
270        If False (default), the given `us` is assumed to be an array of uniformly distributed thresholds, which are taken as the central values of the bins.
271        If True, input `us` is assumed to be a monotonically increasing array of bin edges, including the rightmost edge, allowing for non-uniform distributions of thresholds. 
272        In the latter case, the effective thresholds are the central value of the given bins.
273        
274    verbose : bool, optional
275        If True (default), progress bars are shown for the computations on data.
276    
277    Returns
278    -------
279    V3 : np.array()
280        The values of the second Minkowski Functional at the given thresholds.
281    
282    """
283    if field.dim<3:
284        raise ValueError(f"V3 is defined for fields with at least 3 dimensions, but this field is {field.dim}D.")
285        
286    us = np.atleast_1d(us)
287    us, dus = define_ubins(us, edges)
288    
289    if isinstance(field, TheoryField):
290        us = subsample_us(us, dus)
291        try:
292            return np.nanmean(field.V3(us), axis=1)
293        except AttributeError:
294            raise NotImplementedError(f"The theoretical expectation of V3 for {field.name} fields is not implemented. If you know an expression, please get in touch.") from None
295        
296    elif isinstance(field, DataField):
297        try:
298            return field._V3(us, dus, verbose=verbose)
299        except AttributeError:
300            if field.first_der is None:
301                field.get_first_der()
302            if field.second_der is None:
303                field.get_second_der()
304                
305            stat = np.zeros_like(us)
306            curv = general_curvature(field, 2)
307            for ii in tqdm(np.arange(us.shape[0]), disable=not verbose):
308                this_mask = (us[ii] + dus[ii]/2. > field.field) & (us[ii] - dus[ii]/2. <= field.field)
309                stat[ii] = np.mean((curv*this_mask/dus[ii])[field.mask])
310            return _MF_prefactor(field.dim, 3) * stat
311        
312    else:
313        raise TypeError(f"The field must be either TheoryField or DataField (or a subclass).")
314
315
316__all__ = ["V0", "V1", "V2", "V3", "general_curvature"]
317
318__docformat__ = "numpy"