power_space.py 10.4 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.
theos's avatar
theos committed
18

theos's avatar
theos committed
19
20
import numpy as np

21
22
import d2o

23
24
from power_index_factory import PowerIndexFactory

25
from nifty.spaces.space import Space
theos's avatar
theos committed
26
27


Theo Steininger's avatar
Theo Steininger committed
28
class PowerSpace(Space):
Theo Steininger's avatar
Theo Steininger committed
29
30
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
    """ 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
56
57
        This holds the information which pixel of the harmonic partner gets
        mapped to which power bin
Theo Steininger's avatar
Theo Steininger committed
58
    kindex : numpy.ndarray
59
        Sorted array of all k-modes.
Theo Steininger's avatar
Theo Steininger committed
60
    pundex : numpy.ndarray
61
        Flat index of the first occurence of a k-vector with length==kindex[n]
62
        in the k_array.
Theo Steininger's avatar
Theo Steininger committed
63
64
65
    rho : numpy.ndarray
        The amount of k-modes that get mapped to one power bin is given by
        rho.
66
67
68
69
70
71
72
73
    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.
74
75
76
    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
77
78
79
80
81
82
83

    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.

    """
84

85
86
    # ---Overwritten properties and methods---

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

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

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

109
        self._config = power_index['config']
110
111
112
113

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

117
118
119
120
121
        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
122
    def pre_cast(self, x, axes):
Theo Steininger's avatar
Theo Steininger committed
123
124
125
126
127
128
        """ 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.

129
130
131
        Parameters
        ----------
        x : {array-like, function array-like -> array-like}
Theo Steininger's avatar
Theo Steininger committed
132
133
134
135
136
137
            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.

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

143
        """
Theo Steininger's avatar
Theo Steininger committed
144

145
146
147
148
149
        if callable(x):
            return x(self.kindex)
        else:
            return x

150
151
    # ---Mandatory properties and methods---

152
153
154
155
156
157
158
    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']))

159
160
161
    @property
    def harmonic(self):
        return True
162

163
164
    @property
    def shape(self):
165
        return self.kindex.shape
166

167
168
169
170
171
172
173
    @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
174
        return float(reduce(lambda x, y: x*y, self.pindex.shape))
175
176

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

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

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

        if inplace:
            x *= weight
            result_x = x
        else:
            result_x = x*weight
198
199
200

        return result_x

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

207
    def get_fft_smoothing_kernel_function(self, sigma):
208
        raise NotImplementedError(
209
            "There is no fft smoothing function for PowerSpace.")
theos's avatar
theos committed
210

211
212
213
    # ---Added properties and methods---

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

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

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

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

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

245
246
    @property
    def pundex(self):
Theo Steininger's avatar
Theo Steininger committed
247
248
249
250
        """ 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.
        """
251
252
253
254
        return self._pundex

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

    # ---Serialization---

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

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

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

290
291
292
293
        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
294

Theo Steininger's avatar
Theo Steininger committed
295
        new_ps._pindex = repository.get('pindex', hdf5_group)
Jait Dixit's avatar
Jait Dixit committed
296
297
298
        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
299
        new_ps._k_array = repository.get('k_array', hdf5_group)
Jait Dixit's avatar
Jait Dixit committed
300
        new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
301
                                    '_k_array']
Jait Dixit's avatar
Jait Dixit committed
302
303
304
305
306
307
308

        return new_ps


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