Include SpinGaussian and exclude EuclideanGaussian

pynkowski/theory/__init__.py

@@ -2,6 +2,7 @@
 
 - [`base_th`](theory/base_th.html): The base class for theoretical fields, to be used as the base for the other fields.
 - [`gaussian`](theory/gaussian.html): Isotropic Gaussian fields such as CMB temperature or initial density field.
+- [`spingaussian`](theory/spingaussian.html): Isotropic Gaussian fields in the SO(3) formalism, such as the CMB polarization.
 - [`chi2`](theory/chi2.html): Isotropic $\chi^2$ fields, such as the modulus of the CMB polarization.
 
 - There is also a general utilities submodule called [`utils_th`](theory/utils_th.html).
@@ -10,7 +11,7 @@
 '''
 
 from .base_th import TheoryField
-from .gaussian import SphericalGaussian, EuclideanGaussian, Gaussian
+from .gaussian import SphericalGaussian, Gaussian #, EuclideanGaussian
 from .spingaussian import SpinGaussian
 from .chi2 import SphericalChi2, Chi2
 

pynkowski/theory/gaussian.py

@@ -327,124 +327,124 @@ def maxima_dist(self, us):
  k2 = self.mu**2./self.C2
  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.))
 
-class EuclideanGaussian(Gaussian):
- """Class for Isotropic Gaussian fields defined in an Euclidean space (such as a 2D plane or a 3D cube).
+# class EuclideanGaussian(Gaussian):
+#  """Class for Isotropic Gaussian fields defined in an Euclidean space (such as a 2D plane or a 3D cube). 
 
- Parameters
- ----------
- dim : int
- Dimension of the space where the field is defined.
+#  Parameters
+#  ----------
+#  dim : int
+#  Dimension of the space where the field is defined.
 
- Pk : np.array or function
- 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.
+#  Pk : np.array or function
+#  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.
 
- k : np.array
- If `Pk` is an array, wavenumbers corresponding to the power spectrum. If `Pk` is a function, `k` is the array `[kmin, kmax]`.
+#  k : np.array
+#  If `Pk` is an array, wavenumbers corresponding to the power spectrum. If `Pk` is a function, `k` is the array `[kmin, kmax]`.
 
- normalise : bool, optional
- If `True`, normalise the field to unit variance.
- Default : True
+#  normalise : bool, optional
+#  If `True`, normalise the field to unit variance.
+#  Default : True
 
- volume : float, optional
- Hausdorff measure of the space where the field is defined (area if `dim=2`, volume if `dim=3`, and so on).
- Default : 1.
+#  volume : float, optional
+#  Hausdorff measure of the space where the field is defined (area if `dim=2`, volume if `dim=3`, and so on).
+#  Default : 1.
 
- Attributes
- ----------
- dim : int
- Dimension of the space where the field is defined.
+#  Attributes
+#  ----------
+#  dim : int
+#  Dimension of the space where the field is defined.
 
- Pk : function
- Power spectrum of the field, to be evaluated at wavenumbers `k`.
+#  Pk : function
+#  Power spectrum of the field, to be evaluated at wavenumbers `k`.
 
- klim : np.array
- Minimum and maximum wavenumbers to be considered for the power spectrum `Pk`.
+#  klim : np.array
+#  Minimum and maximum wavenumbers to be considered for the power spectrum `Pk`.
 
- name : str
- Name of the field, `"Euclidean Isotropic Gaussian"` by default.
+#  name : str
+#  Name of the field, `"Euclidean Isotropic Gaussian"` by default.
 
- sigma : float
- The standard deviation of the field.
+#  sigma : float
+#  The standard deviation of the field.
 
- mu : float
- The derivative of the covariance function at the origin, times $-2$. Equal to the variance of the first derivatives of the field.
+#  mu : float
+#  The derivative of the covariance function at the origin, times $-2$. Equal to the variance of the first derivatives of the field.
 
- nu : float
- 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).
+#  nu : float
+#  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).
 
- lkc_ambient : np.array
- The values for the Lipschitz–Killing Curvatures of the ambient space.
+#  lkc_ambient : np.array
+#  The values for the Lipschitz–Killing Curvatures of the ambient space.
 
- """ 
- def __init__(self, dim, Pk, k, normalise=True, volume=1.):
- if type(Pk) is np.ndarray:
- self.Pk = interp1d(k, Pk)
- self.klim = np.array([k.min(), k.max()])
- else:
- assert callable(Pk), 'The power spectrum must be a function or an array'
- assert k.shape == (2,), 'The power spectrum must be an array with two elements: `kmin, kmax`'
- self.klim = k
- self.Pk = Pk
+#  """ 
+#  def __init__(self, dim, Pk, k, normalise=True, volume=1.):
+#  if type(Pk) is np.ndarray:
+#  self.Pk = interp1d(k, Pk)
+#  self.klim = np.array([k.min(), k.max()])
+#  else:
+#  assert callable(Pk), 'The power spectrum must be a function or an array'
+#  assert k.shape == (2,), 'The power spectrum must be an array with two elements: `kmin, kmax`'
+#  self.klim = k
+#  self.Pk = Pk
 
- if dim != 3:
- 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.')
- # Are these formulae valid also for other dimensions?
- k2Pk = lambda k: self.Pk(k)*k**2.
- sigma = np.sqrt(quad(k2Pk, self.klim[0], self.klim[1])[0] / (2.*np.pi**2.)) #* volume
-
- k4Pk = lambda k: self.Pk(k)*k**4.
- mu = quad(k4Pk, self.klim[0], self.klim[1])[0] / (6.*np.pi**2.) #* volume
-
- k6Pk = lambda k: self.Pk(k)*k**6.
- nu = quad(k6Pk, self.klim[0], self.klim[1])[0] / (120.*np.pi**2.) #* volume
-
- if normalise:
- mu = mu/sigma**2.
- nu = nu/sigma**2.
- self.Pk = lambda k: self.Pk(k)/sigma**2.
- sigma = 1.
+#  if dim != 3:
+#  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.')
+#  # Are these formulae valid also for other dimensions?
+#  k2Pk = lambda k: self.Pk(k)*k**2.
+#  sigma = np.sqrt(quad(k2Pk, self.klim[0], self.klim[1])[0] / (2.*np.pi**2.)) #* volume
+
+#  k4Pk = lambda k: self.Pk(k)*k**4.
+#  mu = quad(k4Pk, self.klim[0], self.klim[1])[0] / (6.*np.pi**2.) #* volume
+
+#  k6Pk = lambda k: self.Pk(k)*k**6.
+#  nu = quad(k6Pk, self.klim[0], self.klim[1])[0] / (120.*np.pi**2.) #* volume
+
+#  if normalise:
+#  mu = mu/sigma**2.
+#  nu = nu/sigma**2.
+#  self.Pk = lambda k: self.Pk(k)/sigma**2.
+#  sigma = 1.
 
- lkc_ambient = np.zeros(dim+1)
- lkc_ambient[-1] = volume
- super().__init__(dim=dim, sigma=sigma, mu=mu, nu=nu, lkc_ambient=lkc_ambient)
- self.name = 'Euclidean Isotropic Gaussian'
- self._warn_non_euclidean = False
+#  lkc_ambient = np.zeros(dim+1)
+#  lkc_ambient[-1] = volume
+#  super().__init__(dim=dim, sigma=sigma, mu=mu, nu=nu, lkc_ambient=lkc_ambient)
+#  self.name = 'Euclidean Isotropic Gaussian'
+#  self._warn_non_euclidean = False
 
- def maxima_total(self):
- """Compute the expected values of local maxima of the field.
+#  def maxima_total(self):
+#  """Compute the expected values of local maxima of the field.
 
- Returns
- -------
- number_total : float
- Expected value of the number of local maxima.
+#  Returns
+#  -------
+#  number_total : float
+#  Expected value of the number of local maxima.
 
- """
- rho1 = -0.5*self.mu
- 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]
+#  """
+#  rho1 = -0.5*self.mu
+#  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]
 
- def maxima_dist(self, us):
- """Compute the expected distribution of local maxima of the field, as a function of threshold.
+#  def maxima_dist(self, us):
+#  """Compute the expected distribution of local maxima of the field, as a function of threshold.
 
- Parameters
- ----------
- us : np.array
- The thresholds considered for the computation.
+#  Parameters
+#  ----------
+#  us : np.array
+#  The thresholds considered for the computation.
 
- Returns
- -------
- density : np.array
- Density distribution of the local maxima.
+#  Returns
+#  -------
+#  density : np.array
+#  Density distribution of the local maxima.
 
- """
- rho1 = -0.5*self.mu
- κ = -rho1 / np.sqrt(self.nu)
- assert κ <= (self.dim + 2.)/self.dim, '`0.5*mu/sqrt(nu)` must be $≤ (dim+2)/dim$'
- 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.)
+#  """
+#  rho1 = -0.5*self.mu
+#  κ = -rho1 / np.sqrt(self.nu)
+#  assert κ <= (self.dim + 2.)/self.dim, '`0.5*mu/sqrt(nu)` must be $≤ (dim+2)/dim$'
+#  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.)
 
 
 
 
-__all__ = ["SphericalGaussian", "EuclideanGaussian", "Gaussian"]
+__all__ = ["SphericalGaussian", "Gaussian"] #, "EuclideanGaussian", 
 
 __docformat__ = "numpy"