power_space.py 10.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# 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 <http://www.gnu.org/licenses/>.
Theo Steininger's avatar
Theo Steininger committed
13
14
15
16
17
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
Theo Steininger's avatar
Theo Steininger committed
18

Martin Reinecke's avatar
fix    
Martin Reinecke committed
19
# from builtins import str
20
21
import numpy as np

22
23
import d2o

Martin Reinecke's avatar
Martin Reinecke committed
24
from .power_index_factory import PowerIndexFactory
25

Martin Reinecke's avatar
Martin Reinecke committed
26
from ..space import Space
Martin Reinecke's avatar
Martin Reinecke committed
27
from functools import reduce
Theo Steininger's avatar
Theo Steininger committed
28
29


Theo Steininger's avatar
Theo Steininger committed
30
class PowerSpace(Space):
Theo Steininger's avatar
Theo Steininger committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
    """ NIFTY class for spaces of power spectra.

    Parameters
    ----------
    harmonic_partner : Space
        The harmonic Space of which this is the power space.
    distribution_strategy : str *optional*
        The distribution strategy used for the distributed_data_objects
        derived from this PowerSpace, e.g. the pindex.
        (default : 'not')
    logarithmic : bool *optional*
        True if logarithmic binning should be used (default : False).
    nbin : {int, None} *optional*
        The number of bins that should be used for power spectrum binning
        (default : None).
        if nbin == None, then nbin is set to the length of kindex.
    binbounds :  {list, array-like} *optional*
        Array-like inner boundaries of the used bins of the default
        indices.
        (default : None)
        if binbounds == None :
            Calculates the bounds from the kindex while applying the
            logarithmic and nbin keywords.

    Attributes
    ----------
    pindex : distributed_data_object
58
59
        This holds the information which pixel of the harmonic partner gets
        mapped to which power bin
Theo Steininger's avatar
Theo Steininger committed
60
    kindex : numpy.ndarray
61
        Sorted array of all k-modes.
Theo Steininger's avatar
Theo Steininger committed
62
    pundex : numpy.ndarray
63
        Flat index of the first occurence of a k-vector with length==kindex[n]
64
        in the k_array.
Theo Steininger's avatar
Theo Steininger committed
65
66
67
    rho : numpy.ndarray
        The amount of k-modes that get mapped to one power bin is given by
        rho.
68
69
70
71
72
73
74
75
    dim : np.int
        Total number of dimensionality, i.e. the number of pixels.
    harmonic : bool
        Specifies whether the space is a signal or harmonic space.
    total_volume : np.float
        The total volume of the space.
    shape : tuple of np.ints
        The shape of the space's data array.
76
77
78
    config : {logarithmic, nbin, binbounds}
        Dictionary storing the values for `logarithmic`, `nbin`, and
        `binbounds` that were used during initialization.
Theo Steininger's avatar
Theo Steininger committed
79
80
81
82
83
84
85

    Notes
    -----
    A power space is the result of a projection of a harmonic space where
    k-modes of equal length get mapped to one power index.

    """
86

87
88
    # ---Overwritten properties and methods---

89
    def __init__(self, harmonic_partner,
90
                 distribution_strategy='not',
91
                 logarithmic=False, nbin=None, binbounds=None):
Martin Reinecke's avatar
Martin Reinecke committed
92
        super(PowerSpace, self).__init__()
93
94
        self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
                                  '_k_array']
95

96
        if not isinstance(harmonic_partner, Space):
97
            raise ValueError(
98
99
                "harmonic_partner must be a Space.")
        if not harmonic_partner.harmonic:
100
            raise ValueError(
101
102
                "harmonic_partner must be a harmonic space.")
        self._harmonic_partner = harmonic_partner
103

Jait Dixit's avatar
Jait Dixit committed
104
        power_index = PowerIndexFactory.get_power_index(
105
                        domain=self.harmonic_partner,
Jait Dixit's avatar
Jait Dixit committed
106
                        distribution_strategy=distribution_strategy,
107
                        logarithmic=logarithmic,
Jait Dixit's avatar
Jait Dixit committed
108
109
                        nbin=nbin,
                        binbounds=binbounds)
110

111
        self._config = power_index['config']
112
113
114
115

        self._pindex = power_index['pindex']
        self._kindex = power_index['kindex']
        self._rho = power_index['rho']
116
117
        self._pundex = power_index['pundex']
        self._k_array = power_index['k_array']
118

119
120
121
122
123
        if self.config['nbin'] is not None:
            if self.config['nbin'] > len(self.kindex):
                self.logger.warn("nbin was set to a value being larger than "
                                 "the length of kindex!")

Theo Steininger's avatar
Theo Steininger committed
124
    def pre_cast(self, x, axes):
Theo Steininger's avatar
Theo Steininger committed
125
126
127
128
129
130
        """ Casts power spectrum functions to discretized power spectra.

        This function takes an array or a function. If it is an array it does
        nothing, otherwise it interpretes the function as power spectrum and
        evaluates it at every k-mode.

131
132
133
        Parameters
        ----------
        x : {array-like, function array-like -> array-like}
Theo Steininger's avatar
Theo Steininger committed
134
135
136
137
138
139
            power spectrum given either in discretized form or implicitly as a
            function
        axes : tuple of ints
            Specifies the axes of x which correspond to this space. For
            explicifying the power spectrum function, this is ignored.

140
141
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
142
143
144
        array-like
            discretized power spectrum

145
        """
Theo Steininger's avatar
Theo Steininger committed
146

147
148
149
150
151
        if callable(x):
            return x(self.kindex)
        else:
            return x

152
153
    # ---Mandatory properties and methods---

154
155
156
157
158
159
160
    def __repr__(self):
        return ("PowerSpace(harmonic_partner=%r, distribution_strategy=%r, "
                "logarithmic=%r, nbin=%r, binbounds=%r)"
                % (self.harmonic_partner, self.pindex.distribution_strategy,
                   self.config['logarithmic'], self.config['nbin'],
                   self.config['binbounds']))

161
162
163
    @property
    def harmonic(self):
        return True
164

165
166
    @property
    def shape(self):
167
        return self.kindex.shape
168

169
170
171
172
173
174
175
    @property
    def dim(self):
        return self.shape[0]

    @property
    def total_volume(self):
        # every power-pixel has a volume of 1
Jait Dixit's avatar
Jait Dixit committed
176
        return float(reduce(lambda x, y: x*y, self.pindex.shape))
177
178

    def copy(self):
179
        distribution_strategy = self.pindex.distribution_strategy
180
        return self.__class__(harmonic_partner=self.harmonic_partner,
181
                              distribution_strategy=distribution_strategy,
Theo Steininger's avatar
Theo Steininger committed
182
183
184
                              logarithmic=self.config["logarithmic"],
                              nbin=self.config["nbin"],
                              binbounds=self.config["binbounds"])
185

186
    def weight(self, x, power=1, axes=None, inplace=False):
Jait Dixit's avatar
Jait Dixit committed
187
188
        reshaper = [1, ] * len(x.shape)
        # we know len(axes) is always 1
189
190
        reshaper[axes[0]] = self.shape[0]

191
        weight = self.rho.reshape(reshaper)
192
        if power != 1:
193
            weight = weight ** np.float(power)
194
195
196
197
198
199

        if inplace:
            x *= weight
            result_x = x
        else:
            result_x = x*weight
200
201
202

        return result_x

203
    def get_distance_array(self, distribution_strategy):
204
        result = d2o.distributed_data_object(
Martin Reinecke's avatar
Martin Reinecke committed
205
                                self.kindex, dtype=np.float64,
206
207
                                distribution_strategy=distribution_strategy)
        return result
208

209
    def get_fft_smoothing_kernel_function(self, sigma):
210
        raise NotImplementedError(
211
            "There is no fft smoothing function for PowerSpace.")
212

213
214
215
    # ---Added properties and methods---

    @property
216
    def harmonic_partner(self):
Theo Steininger's avatar
Theo Steininger committed
217
        """ Returns the Space of which this is the power space.
218
219
        """
        return self._harmonic_partner
220
221

    @property
222
223
224
    def config(self):
        """ Returns the configuration which was used for `logarithmic`, `nbin`
        and `binbounds` during initialization.
225
        """
226
        return self._config
227
228
229

    @property
    def pindex(self):
230
        """ A distributed_data_object having the shape of the harmonic partner
Theo Steininger's avatar
Theo Steininger committed
231
232
        space containing the indices of the power bin a pixel belongs to.
        """
233
234
235
236
        return self._pindex

    @property
    def kindex(self):
Theo Steininger's avatar
Theo Steininger committed
237
238
        """ Sorted array of all k-modes.
        """
239
240
241
242
        return self._kindex

    @property
    def rho(self):
Theo Steininger's avatar
Theo Steininger committed
243
244
        """Degeneracy factor of the individual k-vectors.
        """
245
        return self._rho
246

247
248
    @property
    def pundex(self):
Theo Steininger's avatar
Theo Steininger committed
249
250
251
252
        """ An array for which the n-th entry gives the flat index of the
        first occurence of a k-vector with length==kindex[n] in the
        k_array.
        """
253
254
255
256
        return self._pundex

    @property
    def k_array(self):
Theo Steininger's avatar
Theo Steininger committed
257
258
259
        """ An array containing distances to the grid center (i.e. zero-mode)
        for every k-mode in the grid of the harmonic partner space.
        """
260
        return self._k_array
261
262
263
264

    # ---Serialization---

    def _to_hdf5(self, hdf5_group):
Jait Dixit's avatar
Jait Dixit committed
265
        hdf5_group['kindex'] = self.kindex
266
267
        hdf5_group['rho'] = self.rho
        hdf5_group['pundex'] = self.pundex
Theo Steininger's avatar
Theo Steininger committed
268
        hdf5_group['logarithmic'] = self.config["logarithmic"]
Theo Steininger's avatar
Theo Steininger committed
269
        # Store nbin as string, since it can be None
270
271
        hdf5_group.attrs['nbin'] = str(self.config["nbin"])
        hdf5_group.attrs['binbounds'] = str(self.config["binbounds"])
272

273
        #MR FIXME: why not "return None" as happens everywhere else?
274
        return {
275
            'harmonic_partner': self.harmonic_partner,
276
277
278
279
280
            'pindex': self.pindex,
            'k_array': self.k_array
        }

    @classmethod
Theo Steininger's avatar
Theo Steininger committed
281
    def _from_hdf5(cls, hdf5_group, repository):
Jait Dixit's avatar
Jait Dixit committed
282
283
284
285
        # make an empty PowerSpace object
        new_ps = EmptyPowerSpace()
        # reset class
        new_ps.__class__ = cls
Jait Dixit's avatar
Jait Dixit committed
286
        # call instructor so that classes are properly setup
Martin Reinecke's avatar
Martin Reinecke committed
287
        super(PowerSpace, new_ps).__init__()
Jait Dixit's avatar
Jait Dixit committed
288
        # set all values
Theo Steininger's avatar
Theo Steininger committed
289
290
        new_ps._harmonic_partner = repository.get('harmonic_partner',
                                                  hdf5_group)
Theo Steininger's avatar
Theo Steininger committed
291

292
293
294
295
        new_ps._config = {}
        new_ps._config['logarithmic'] = hdf5_group['logarithmic'][()]
        exec("new_ps._config['nbin'] = " + hdf5_group.attrs['nbin'])
        exec("new_ps._config['binbounds'] = " + hdf5_group.attrs['binbounds'])
Jait Dixit's avatar
Jait Dixit committed
296

Theo Steininger's avatar
Theo Steininger committed
297
        new_ps._pindex = repository.get('pindex', hdf5_group)
Jait Dixit's avatar
Jait Dixit committed
298
299
300
        new_ps._kindex = hdf5_group['kindex'][:]
        new_ps._rho = hdf5_group['rho'][:]
        new_ps._pundex = hdf5_group['pundex'][:]
Theo Steininger's avatar
Theo Steininger committed
301
        new_ps._k_array = repository.get('k_array', hdf5_group)
Jait Dixit's avatar
Jait Dixit committed
302
        new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
303
                                    '_k_array']
Jait Dixit's avatar
Jait Dixit committed
304
305
306
307
308
309
310

        return new_ps


class EmptyPowerSpace(PowerSpace):
    def __init__(self):
        pass