lm_space.py 5.36 KB
 Theo Steininger committed Apr 13, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 # NIFTy # Copyright (C) 2017 Theo Steininger # # Author: Theo Steininger # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see .  csongor committed Jul 04, 2016 19 20 21 22 from __future__ import division import numpy as np  theos committed Jul 23, 2016 23 from nifty.spaces.space import Space  theos committed Jul 23, 2016 24   Jait Dixit committed Sep 17, 2016 25 26 from d2o import arange  csongor committed Jul 04, 2016 27   Theo Steininger committed Dec 05, 2016 28 class LMSpace(Space):  csongor committed Jul 04, 2016 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51  """ .. __ .. / / .. / / __ ____ ___ .. / / / _ _ | .. / /_ / / / / / / .. \___/ /__/ /__/ /__/ space class NIFTY subclass for spherical harmonics components, for representations of fields on the two-sphere. Parameters ---------- lmax : int Maximum :math:\ell-value up to which the spherical harmonics coefficients are to be used. See Also -------- hp_space : A class for the HEALPix discretization of the sphere [#]_. gl_space : A class for the Gauss-Legendre discretization of the sphere [#]_.  Theo Steininger committed May 11, 2017 52 53 54 55 56  Raises ------ ValueError If given lmax is negative.  Reimar H Leike committed May 10, 2017 57 58  Notes -----  Theo Steininger committed May 11, 2017 59  This implementation implicitly sets the mmax parameter to lmax.  Reimar H Leike committed May 10, 2017 60   csongor committed Jul 04, 2016 61 62 63 64 65 66 67 68 69 70  References ---------- .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere", *ApJ* 622..759G. .. [#] M. Reinecke and D. Sverre Seljebotn, 2013, "Libsharp - spherical harmonic transforms revisited"; arXiv:1303.4945 _ """  Martin Reinecke committed Apr 28, 2017 71 72  def __init__(self, lmax): super(LMSpace, self).__init__()  csongor committed Sep 02, 2016 73  self._lmax = self._parse_lmax(lmax)  csongor committed Jul 04, 2016 74   Theo Steininger committed Feb 25, 2017 75 76  def hermitian_decomposition(self, x, axes=None, preserve_gaussian_variance=False):  Theo Steininger committed Jan 11, 2017 77 78 79 80 81 82  hermitian_part = x.copy_empty() anti_hermitian_part = x.copy_empty() hermitian_part[:] = x.real anti_hermitian_part[:] = x.imag return (hermitian_part, anti_hermitian_part)  Theo Steininger committed Apr 21, 2017 83  # ---Mandatory properties and methods---  csongor committed Sep 02, 2016 84 85 86 87  @property def harmonic(self): return True  csongor committed Jul 04, 2016 88 89  @property  theos committed Sep 02, 2016 90 91  def shape(self): return (self.dim, )  csongor committed Jul 04, 2016 92 93  @property  theos committed Sep 02, 2016 94 95  def dim(self): l = self.lmax  Martin Reinecke committed Apr 20, 2017 96  # the LMSpace consists of the full triangle (including -m's!),  theos committed Sep 08, 2016 97 98  # minus two little triangles if mmax < lmax # dim = (((2*(l+1)-1)+1)**2/4 - 2 * (l-m)(l-m+1)/2  Theo Steininger committed Jan 24, 2017 99 100 101  # dim = np.int((l+1)**2 - (l-m)*(l-m+1.)) # We fix l == m return np.int((l+1)**2)  csongor committed Jul 04, 2016 102   theos committed Sep 02, 2016 103 104 105  @property def total_volume(self): # the individual pixels have a fixed volume of 1.  Martin Reinecke committed Apr 28, 2017 106  return np.float64(self.dim)  csongor committed Jul 04, 2016 107   theos committed Sep 02, 2016 108  def copy(self):  Martin Reinecke committed Apr 28, 2017 109  return self.__class__(lmax=self.lmax)  csongor committed Jul 04, 2016 110   theos committed Sep 02, 2016 111 112 113  def weight(self, x, power=1, axes=None, inplace=False): if inplace: return x  csongor committed Jul 04, 2016 114  else:  theos committed Sep 02, 2016 115  return x.copy()  csongor committed Sep 02, 2016 116   theos committed Oct 12, 2016 117 118 119 120 121  def get_distance_array(self, distribution_strategy): dists = arange(start=0, stop=self.shape[0], distribution_strategy=distribution_strategy) dists = dists.apply_scalar_function(  Theo Steininger committed Apr 21, 2017 122  lambda x: self._distance_array_helper(x, self.lmax),  Martin Reinecke committed Apr 28, 2017 123  dtype=np.float64)  theos committed Oct 12, 2016 124 125 126  return dists  Theo Steininger committed Apr 21, 2017 127 128 129 130 131 132 133 134  @staticmethod def _distance_array_helper(index_array, lmax): u = 2*lmax + 1 index_half = (index_array+np.minimum(lmax, index_array)+1)//2 m = (np.ceil((u - np.sqrt(u*u - 8*(index_half - lmax)))/2)).astype(int) res = (index_half - m*(u - m)//2).astype(np.float64) return res  theos committed Oct 27, 2016 135  def get_fft_smoothing_kernel_function(self, sigma):  Reimar H Leike committed May 10, 2017 136  # FIXME why x(x+1) ? add reference to paper!  theos committed Oct 12, 2016 137 138  return lambda x: np.exp(-0.5 * x * (x + 1) * sigma**2)  csongor committed Sep 02, 2016 139 140 141 142  # ---Added properties and methods--- @property def lmax(self):  Theo Steininger committed May 11, 2017 143 144  """ Returns the maximal :math:l value of any spherical harmonics :math:Y_{lm} that is represented in this Space.  145  """  csongor committed Sep 02, 2016 146 147 148 149  return self._lmax @property def mmax(self):  Theo Steininger committed May 11, 2017 150 151 152 153 154  """ Returns the maximal :math:m value of any spherical harmonic :math:Y_{lm} that is represented in this Space. As :math:m goes from :math:-l to :math:l for every :math:l this just returns the same as lmax.  155 156  See Also --------  Theo Steininger committed May 11, 2017 157 158 159  lmax : Returns the maximal :math:l-value of the spherical harmonics being used.  160  """  Theo Steininger committed May 11, 2017 161   Jait Dixit committed Sep 21, 2016 162  return self._lmax  csongor committed Sep 02, 2016 163 164 165  def _parse_lmax(self, lmax): lmax = np.int(lmax)  Martin Reinecke committed Apr 20, 2017 166 167  if lmax < 0: raise ValueError("lmax must be >=0.")  csongor committed Sep 02, 2016 168  return lmax  Jait Dixit committed Dec 03, 2016 169 170 171 172 173 174 175 176  # ---Serialization--- def _to_hdf5(self, hdf5_group): hdf5_group['lmax'] = self.lmax return None @classmethod  Theo Steininger committed Dec 05, 2016 177  def _from_hdf5(cls, hdf5_group, repository):  Jait Dixit committed Dec 03, 2016 178 179 180 181  result = cls( lmax=hdf5_group['lmax'][()], ) return result