power_space.py 10.6 KB
 Theo Steininger committed Apr 13, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 # 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 .  theos committed Jul 02, 2016 18   theos committed Jul 04, 2016 19 20 import numpy as np  theos committed Oct 27, 2016 21 22 import d2o  theos committed Aug 25, 2016 23 24 from power_index_factory import PowerIndexFactory  theos committed Jul 23, 2016 25 from nifty.spaces.space import Space  theos committed Aug 25, 2016 26 from nifty.spaces.rg_space import RGSpace  theos committed Jul 02, 2016 27 28   Theo Steininger committed Dec 05, 2016 29 class PowerSpace(Space):  Jait Dixit committed Nov 15, 2016 30   theos committed Aug 25, 2016 31 32  # ---Overwritten properties and methods---  Reimar H Leike committed May 10, 2017 33  def __init__(self, harmonic_partner=RGSpace((1,)),  theos committed Sep 13, 2016 34  distribution_strategy='not',  Reimar H Leike committed May 10, 2017 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68  logarithmic=False, nbin=None, binbounds=None): """Sets the attributes for a PowerSpace class instance. Parameters ---------- harmonic_partner : Space The harmonic Space of which this is the power space. distribution_strategy : str *optional* The distribution strategy of a d2o-object represeting a field over this PowerSpace. (default : 'not') logarithmic : bool *optional* True if logarithmic binning should be used. (default : False) nbin : {int, None} *optional* The number of bins this space has. (default : None) if nbin == None : It takes the nbin from its harmonic_partner binbounds : {list, array} *optional* Array-like inner boundaries of the used bins of the default indices. (default : None) if binbounds == None : Calculates the bounds from the kindex and corrects for logartihmic scale Notes ----- A power space is the result of a projection of a harmonic space where multiple k-modes get mapped to one power index. This can be regarded as a response operator :math:R going from harmonic space to power space. An array giving this map is stored in pindex (array which says in which power box a k-mode gets projected) An array for the adjoint of :math:R is given by kindex, which is an array of arrays stating which k-mode got mapped to a power index The a right-inverse to :math:R is given by the pundex which is an array giving one k-mode that maps to a power bin for every power bin. The amount of k-modes that get mapped to one power bin is given by rho. This is :math:RR^\dagger in the language of this projection operator Returns ------- None. """ #FIXME: default probably not working for log and normal scale  Martin Reinecke committed Apr 28, 2017 69  super(PowerSpace, self).__init__()  theos committed Oct 17, 2016 70 71  self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex', '_k_array']  theos committed Aug 25, 2016 72   Reimar H Leike committed May 10, 2017 73  if not isinstance(harmonic_partner, Space):  theos committed Oct 06, 2016 74  raise ValueError(  Reimar H Leike committed May 10, 2017 75 76  "harmonic_partner must be a Space.") if not harmonic_partner.harmonic:  theos committed Oct 06, 2016 77  raise ValueError(  Reimar H Leike committed May 10, 2017 78 79  "harmonic_partner must be a harmonic space.") self._harmonic_partner = harmonic_partner  theos committed Aug 25, 2016 80   Jait Dixit committed Dec 03, 2016 81  power_index = PowerIndexFactory.get_power_index(  Reimar H Leike committed May 10, 2017 82  domain=self.harmonic_partner,  Jait Dixit committed Dec 03, 2016 83  distribution_strategy=distribution_strategy,  Reimar H Leike committed May 10, 2017 84  logarithmic=logarithmic,  Jait Dixit committed Dec 03, 2016 85 86  nbin=nbin, binbounds=binbounds)  theos committed Aug 25, 2016 87 88  config = power_index['config']  Reimar H Leike committed May 10, 2017 89  self._logarithmic = config['logarithmic']  theos committed Aug 25, 2016 90 91 92 93 94 95  self._nbin = config['nbin'] self._binbounds = config['binbounds'] self._pindex = power_index['pindex'] self._kindex = power_index['kindex'] self._rho = power_index['rho']  theos committed Aug 26, 2016 96 97  self._pundex = power_index['pundex'] self._k_array = power_index['k_array']  theos committed Aug 25, 2016 98   Theo Steininger committed May 11, 2017 99  def pre_cast(self, x, axes):  Reimar H Leike committed May 10, 2017 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115  """Casts power spectra to discretized power spectra. This function takes an array or a function. If it is an array it does nothing, otherwise it intepretes the function as power spectrum and evaluates it at every k-mode. Parameters ---------- x : {array-like, function array-like -> array-like} power spectrum given either in discretized form or implicitly as a function axes : {tuple, int} *optional* does nothing (default : None) Returns ------- array-like : discretized power spectrum """  theos committed Sep 13, 2016 116 117 118 119 120  if callable(x): return x(self.kindex) else: return x  theos committed Aug 25, 2016 121 122 123 124 125  # ---Mandatory properties and methods--- @property def harmonic(self): return True  126   theos committed Jul 22, 2016 127 128  @property def shape(self):  theos committed Aug 25, 2016 129  return self.kindex.shape  theos committed Jul 22, 2016 130   theos committed Jul 22, 2016 131 132 133 134 135 136 137  @property def dim(self): return self.shape[0] @property def total_volume(self): # every power-pixel has a volume of 1  Jait Dixit committed Feb 12, 2017 138  return float(reduce(lambda x, y: x*y, self.pindex.shape))  theos committed Aug 25, 2016 139 140  def copy(self):  theos committed Sep 02, 2016 141  distribution_strategy = self.pindex.distribution_strategy  Reimar H Leike committed May 10, 2017 142  return self.__class__(harmonic_partner=self.harmonic_partner,  theos committed Sep 02, 2016 143  distribution_strategy=distribution_strategy,  Reimar H Leike committed May 10, 2017 144  logarithmic=self.logarithmic,  theos committed Aug 25, 2016 145  nbin=self.nbin,  Martin Reinecke committed Apr 28, 2017 146  binbounds=self.binbounds)  theos committed Jul 22, 2016 147   theos committed Jul 23, 2016 148  def weight(self, x, power=1, axes=None, inplace=False):  Jait Dixit committed Feb 12, 2017 149 150  reshaper = [1, ] * len(x.shape) # we know len(axes) is always 1  theos committed Jul 22, 2016 151 152  reshaper[axes[0]] = self.shape[0]  theos committed Aug 25, 2016 153  weight = self.rho.reshape(reshaper)  theos committed Jul 22, 2016 154 155  if power != 1: weight = weight ** power  theos committed Jul 23, 2016 156 157 158 159 160 161  if inplace: x *= weight result_x = x else: result_x = x*weight  theos committed Jul 22, 2016 162 163 164  return result_x  theos committed Sep 30, 2016 165  def get_distance_array(self, distribution_strategy):  theos committed Oct 27, 2016 166  result = d2o.distributed_data_object(  Martin Reinecke committed Apr 20, 2017 167  self.kindex, dtype=np.float64,  theos committed Oct 27, 2016 168 169  distribution_strategy=distribution_strategy) return result  theos committed Sep 15, 2016 170   theos committed Oct 27, 2016 171  def get_fft_smoothing_kernel_function(self, sigma):  theos committed Oct 06, 2016 172  raise NotImplementedError(  theos committed Oct 27, 2016 173  "There is no fft smoothing function for PowerSpace.")  theos committed Sep 15, 2016 174   theos committed Aug 25, 2016 175 176 177  # ---Added properties and methods--- @property  Reimar H Leike committed May 10, 2017 178 179 180 181 182 183 184  def harmonic_partner(self): """Returns the Space of which this is the power space. Returns ------- Space : The harmonic Space of which this is the power space. """ return self._harmonic_partner  theos committed Aug 25, 2016 185 186  @property  Reimar H Leike committed May 10, 2017 187 188 189 190 191 192 193  def logarithmic(self): """Returns a True if logarithmic binning is used. Returns ------- Bool : True if for this PowerSpace logarithmic binning is used. """ return self._logarithmic  theos committed Aug 25, 2016 194 195 196  @property def nbin(self):  Reimar H Leike committed May 10, 2017 197 198 199 200 201  """Returns the number of power bins. Returns ------- int : The number of bins this space has. """  theos committed Aug 25, 2016 202 203 204 205  return self._nbin @property def binbounds(self):  Reimar H Leike committed May 10, 2017 206 207 208 209 210 211 212 213  """ Inner boundaries of the used bins of the default indices. Returns ------- {list, array} : the inner boundaries of the used bins in the used scale, as they were set in __init__ or computed. """ # FIXME check wether this returns something sensible if 'None' was set in __init__  theos committed Aug 25, 2016 214 215 216 217  return self._binbounds @property def pindex(self):  Reimar H Leike committed May 10, 2017 218 219 220 221 222  """Index of the Fourier grid points that belong to a specific power index Returns ------- distributed_data_object : Index of the Fourier grid points in a distributed_data_object. """  theos committed Aug 25, 2016 223 224 225 226  return self._pindex @property def kindex(self):  Reimar H Leike committed May 10, 2017 227 228 229 230 231  """Array of all k-vector lengths. Returns ------- ndarray : Array which states for each k-mode which power index it maps to (adjoint to pindex) """  theos committed Aug 25, 2016 232 233 234 235  return self._kindex @property def rho(self):  Reimar H Leike committed May 10, 2017 236 237 238 239  """Degeneracy factor of the individual k-vectors. ndarray : Array stating how many k-modes are mapped to one power index for every power index """  theos committed Aug 25, 2016 240  return self._rho  theos committed Jul 22, 2016 241   theos committed Aug 26, 2016 242 243  @property def pundex(self):  Reimar H Leike committed May 10, 2017 244 245 246 247 248  """List of one k-mode per power bin which is in the bin. Returns ------- array-like : An array for which the n-th entry is an example one k-mode which belongs to the n-th power bin """  theos committed Aug 26, 2016 249 250 251 252  return self._pundex @property def k_array(self):  Reimar H Leike committed May 10, 2017 253 254 255 256 257 258  """This contains distances to zero for every k-mode of the harmonic partner. Returns ------- array-like : An array containing distances to the zero mode for every k-mode of the harmonic partner. """  theos committed Aug 26, 2016 259  return self._k_array  Jait Dixit committed Nov 15, 2016 260 261 262 263  # ---Serialization--- def _to_hdf5(self, hdf5_group):  Jait Dixit committed Dec 03, 2016 264 265 266  hdf5_group['kindex'] = self.kindex hdf5_group['rho'] = self.rho hdf5_group['pundex'] = self.pundex  Reimar H Leike committed May 10, 2017 267  hdf5_group['logarithmic'] = self.logarithmic  Theo Steininger committed Dec 05, 2016 268 269 270  # Store nbin as string, since it can be None hdf5_group.attrs['nbin'] = str(self.nbin) hdf5_group.attrs['binbounds'] = str(self.binbounds)  Jait Dixit committed Nov 15, 2016 271 272  return {  Reimar H Leike committed May 10, 2017 273  'harmonic_partner': self.harmonic_partner,  Jait Dixit committed Nov 15, 2016 274 275 276 277 278  'pindex': self.pindex, 'k_array': self.k_array } @classmethod  Theo Steininger committed Dec 05, 2016 279  def _from_hdf5(cls, hdf5_group, repository):  Jait Dixit committed Dec 03, 2016 280 281 282 283  # make an empty PowerSpace object new_ps = EmptyPowerSpace() # reset class new_ps.__class__ = cls  Jait Dixit committed Feb 12, 2017 284  # call instructor so that classes are properly setup  Martin Reinecke committed Apr 28, 2017 285  super(PowerSpace, new_ps).__init__()  Jait Dixit committed Dec 03, 2016 286  # set all values  Reimar H Leike committed May 10, 2017 287 288  new_ps._harmonic_partner = repository.get('harmonic_partner', hdf5_group) new_ps._logarithmic = hdf5_group['logarithmic'][()]  Theo Steininger committed Dec 05, 2016 289 290  exec('new_ps._nbin = ' + hdf5_group.attrs['nbin']) exec('new_ps._binbounds = ' + hdf5_group.attrs['binbounds'])  Jait Dixit committed Dec 03, 2016 291   Theo Steininger committed Dec 05, 2016 292  new_ps._pindex = repository.get('pindex', hdf5_group)  Jait Dixit committed Dec 03, 2016 293 294 295  new_ps._kindex = hdf5_group['kindex'][:] new_ps._rho = hdf5_group['rho'][:] new_ps._pundex = hdf5_group['pundex'][:]  Theo Steininger committed Dec 05, 2016 296  new_ps._k_array = repository.get('k_array', hdf5_group)  Jait Dixit committed Feb 12, 2017 297  new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',  Theo Steininger committed May 01, 2017 298  '_k_array']  Jait Dixit committed Dec 03, 2016 299 300 301 302 303 304 305  return new_ps class EmptyPowerSpace(PowerSpace): def __init__(self): pass