power_space.py 11 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

Martin Reinecke's avatar
Martin Reinecke committed
21
from d2o import distributed_data_object
22

23
from nifty.spaces.space import Space
theos's avatar
theos committed
24
25


Theo Steininger's avatar
Theo Steininger committed
26
class PowerSpace(Space):
Theo Steininger's avatar
Theo Steininger committed
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
    """ 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*
Martin Reinecke's avatar
Martin Reinecke committed
44
        Boundaries between the power spectrum bins.
Martin Reinecke's avatar
Martin Reinecke committed
45
46
        (If binbounds has n entries, there will be n+1 bins, the first bin
        starting at -inf, the last bin ending at +inf.)
Theo Steininger's avatar
Theo Steininger committed
47
48
49
50
51
52
53
54
        (default : None)
        if binbounds == None :
            Calculates the bounds from the kindex while applying the
            logarithmic and nbin keywords.

    Attributes
    ----------
    pindex : distributed_data_object
55
56
        This holds the information which pixel of the harmonic partner gets
        mapped to which power bin
Theo Steininger's avatar
Theo Steininger committed
57
    kindex : numpy.ndarray
58
        Sorted array of all k-modes.
Theo Steininger's avatar
Theo Steininger committed
59
60
61
    rho : numpy.ndarray
        The amount of k-modes that get mapped to one power bin is given by
        rho.
62
63
64
65
66
67
68
69
    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.
Martin Reinecke's avatar
step 4    
Martin Reinecke committed
70
    binbounds :  {list, array-like}
Martin Reinecke's avatar
Martin Reinecke committed
71
        Boundaries between the power spectrum bins
Theo Steininger's avatar
Theo Steininger committed
72
73
74
75
76
77
78

    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.

    """
79

80
81
    # ---Overwritten properties and methods---

82
    def __init__(self, harmonic_partner,
83
                 distribution_strategy='not',
84
                 logarithmic=False, nbin=None, binbounds=None):
Martin Reinecke's avatar
Martin Reinecke committed
85
        super(PowerSpace, self).__init__()
Martin Reinecke's avatar
step 1    
Martin Reinecke committed
86
        self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_k_array']
87

Martin Reinecke's avatar
Martin Reinecke committed
88
89
90
        if not (isinstance(harmonic_partner, Space) and \
                harmonic_partner.harmonic):
            raise ValueError("harmonic_partner must be a harmonic space.")
91
        self._harmonic_partner = harmonic_partner
Martin Reinecke's avatar
step x    
Martin Reinecke committed
92

Martin Reinecke's avatar
Martin Reinecke committed
93
94
        self._k_array = self._harmonic_partner.get_distance_array(distribution_strategy)
        self._compute_binbounds (binbounds, logarithmic, nbin)
Martin Reinecke's avatar
step x    
Martin Reinecke committed
95

Martin Reinecke's avatar
Martin Reinecke committed
96
        self._do_binning()
97

Theo Steininger's avatar
Theo Steininger committed
98
    def pre_cast(self, x, axes):
Theo Steininger's avatar
Theo Steininger committed
99
100
101
102
103
104
        """ 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.

105
106
107
        Parameters
        ----------
        x : {array-like, function array-like -> array-like}
Theo Steininger's avatar
Theo Steininger committed
108
109
110
111
112
113
            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.

114
115
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
116
117
118
        array-like
            discretized power spectrum

119
        """
Theo Steininger's avatar
Theo Steininger committed
120

Martin Reinecke's avatar
step x    
Martin Reinecke committed
121
        return x(self.kindex) if callable(x) else x
122

123
124
    # ---Mandatory properties and methods---

125
126
    def __repr__(self):
        return ("PowerSpace(harmonic_partner=%r, distribution_strategy=%r, "
Martin Reinecke's avatar
Martin Reinecke committed
127
                "binbounds=%r)"
128
                % (self.harmonic_partner, self.pindex.distribution_strategy,
Martin Reinecke's avatar
Martin Reinecke committed
129
                   self._binbounds))
130

131
132
133
    @property
    def harmonic(self):
        return True
134

135
136
    @property
    def shape(self):
137
        return self.kindex.shape
138

139
140
141
142
143
144
145
    @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
146
        return float(reduce(lambda x, y: x*y, self.pindex.shape))
147
148

    def copy(self):
149
        distribution_strategy = self.pindex.distribution_strategy
150
        return self.__class__(harmonic_partner=self.harmonic_partner,
151
                              distribution_strategy=distribution_strategy,
Martin Reinecke's avatar
step 3    
Martin Reinecke committed
152
                              binbounds=self._binbounds)
153

154
    def weight(self, x, power=1, axes=None, inplace=False):
Jait Dixit's avatar
Jait Dixit committed
155
156
        reshaper = [1, ] * len(x.shape)
        # we know len(axes) is always 1
157
158
        reshaper[axes[0]] = self.shape[0]

159
        weight = self.rho.reshape(reshaper)
160
        if power != 1:
161
            weight = weight ** np.float(power)
162
163
164
165
166
167

        if inplace:
            x *= weight
            result_x = x
        else:
            result_x = x*weight
168
169
170

        return result_x

171
    def get_distance_array(self, distribution_strategy):
Martin Reinecke's avatar
Martin Reinecke committed
172
        return distributed_data_object(
Martin Reinecke's avatar
Martin Reinecke committed
173
                                self.kindex, dtype=np.float64,
174
                                distribution_strategy=distribution_strategy)
theos's avatar
theos committed
175

176
    def get_fft_smoothing_kernel_function(self, sigma):
177
        raise NotImplementedError(
178
            "There is no fft smoothing function for PowerSpace.")
theos's avatar
theos committed
179

180
181
182
    # ---Added properties and methods---

    @property
183
    def harmonic_partner(self):
Theo Steininger's avatar
Theo Steininger committed
184
        """ Returns the Space of which this is the power space.
185
186
        """
        return self._harmonic_partner
187

Martin Reinecke's avatar
step 3    
Martin Reinecke committed
188
189
190
    @property
    def binbounds(self):
        return self._binbounds
191
192
193

    @property
    def pindex(self):
194
        """ A distributed_data_object having the shape of the harmonic partner
Theo Steininger's avatar
Theo Steininger committed
195
196
        space containing the indices of the power bin a pixel belongs to.
        """
197
198
199
200
        return self._pindex

    @property
    def kindex(self):
Theo Steininger's avatar
Theo Steininger committed
201
202
        """ Sorted array of all k-modes.
        """
203
204
205
206
        return self._kindex

    @property
    def rho(self):
Theo Steininger's avatar
Theo Steininger committed
207
208
        """Degeneracy factor of the individual k-vectors.
        """
209
        return self._rho
210

211
212
    @property
    def k_array(self):
Theo Steininger's avatar
Theo Steininger committed
213
214
215
        """ An array containing distances to the grid center (i.e. zero-mode)
        for every k-mode in the grid of the harmonic partner space.
        """
216
        return self._k_array
217
218
219
220

    # ---Serialization---

    def _to_hdf5(self, hdf5_group):
Jait Dixit's avatar
Jait Dixit committed
221
        hdf5_group['kindex'] = self.kindex
222
        hdf5_group['rho'] = self.rho
Martin Reinecke's avatar
step 3    
Martin Reinecke committed
223
        hdf5_group.attrs['binbounds'] = str(self._binbounds)
224

225
        #MR FIXME: why not "return None" as happens everywhere else?
226
        return {
227
            'harmonic_partner': self.harmonic_partner,
228
229
230
231
232
            'pindex': self.pindex,
            'k_array': self.k_array
        }

    @classmethod
Theo Steininger's avatar
Theo Steininger committed
233
    def _from_hdf5(cls, hdf5_group, repository):
Jait Dixit's avatar
Jait Dixit committed
234
235
236
237
        # make an empty PowerSpace object
        new_ps = EmptyPowerSpace()
        # reset class
        new_ps.__class__ = cls
Jait Dixit's avatar
Jait Dixit committed
238
        # call instructor so that classes are properly setup
Martin Reinecke's avatar
Martin Reinecke committed
239
        super(PowerSpace, new_ps).__init__()
Jait Dixit's avatar
Jait Dixit committed
240
        # set all values
Theo Steininger's avatar
Theo Steininger committed
241
242
        new_ps._harmonic_partner = repository.get('harmonic_partner',
                                                  hdf5_group)
Theo Steininger's avatar
Theo Steininger committed
243

Martin Reinecke's avatar
step 3    
Martin Reinecke committed
244
        exec("new_ps._binbounds = " + hdf5_group.attrs['binbounds'])
Jait Dixit's avatar
Jait Dixit committed
245

Theo Steininger's avatar
Theo Steininger committed
246
        new_ps._pindex = repository.get('pindex', hdf5_group)
Jait Dixit's avatar
Jait Dixit committed
247
248
        new_ps._kindex = hdf5_group['kindex'][:]
        new_ps._rho = hdf5_group['rho'][:]
Theo Steininger's avatar
Theo Steininger committed
249
        new_ps._k_array = repository.get('k_array', hdf5_group)
Martin Reinecke's avatar
step 1    
Martin Reinecke committed
250
        new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_k_array']
Jait Dixit's avatar
Jait Dixit committed
251
252
253

        return new_ps

Martin Reinecke's avatar
Martin Reinecke committed
254
255
256
    def _do_binning(self):
        """ Computes pindex, kindex and rho according to self._binbounds.
        """
Martin Reinecke's avatar
Martin Reinecke committed
257
258
259
260
261
262
        # prepare the pindex object
        self._pindex = distributed_data_object(
            global_shape=self._k_array.shape,
            dtype=np.int,
            distribution_strategy=self._k_array.distribution_strategy)

Martin Reinecke's avatar
Martin Reinecke committed
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
        # if no binning is requested, compute the "natural" binning, i.e. there
        # is one bin for every distinct k-vector length
        if self._binbounds is None:
            self._kindex = self._k_array.unique()
            # store the local pindex data in the global_pindex d2o
            self._pindex.set_local_data(
                np.searchsorted(self._kindex, self._k_array.get_local_data()))
            self._rho = self._pindex.bincount().get_full_data()

        # use the provided binbounds
        else:
            self._pindex.set_local_data(np.searchsorted(
                self._binbounds, self._k_array.get_local_data()))
            self._rho = self._pindex.bincount().get_full_data()
            self._kindex = self._pindex.bincount(
                weights=self._k_array).get_full_data()/self._rho

    def _compute_binbounds (self, binbounds, logarithmic, nbin):
        try:
            self._binbounds = np.sort(tuple(np.array(binbounds)))
        except(TypeError):
            self._binbounds = None

        if(self._binbounds is None):
            try:
                logarithmic = bool(logarithmic)
            except(TypeError):
                logarithmic = False
            try:
                nbin = int(nbin)
            except(TypeError):
                nbin = None

            if not logarithmic and nbin is None:
                return # no binbounds, use "natural" binning

            # equal binning
            # MR FIXME: this needs to improve
            kindex = self._k_array.unique()
            if (logarithmic):
                k = np.r_[0, np.log(kindex[1:])]
            else:
                k = kindex
            dk = np.max(k[2:] - k[1:-1])  # minimum dk to avoid empty bins
            if(nbin is None):
                nbin = int((k[-1] - 0.5 * (k[2] + k[1])) /
                           dk - 0.5)  # maximal nbin
            else:
                nbin = min(int(nbin), int(
                    (k[-1] - 0.5 * (k[2] + k[1])) / dk + 2.5))
                dk = (k[-1] - 0.5 * (k[2] + k[1])) / (nbin - 2.5)
            self._binbounds = np.r_[0.5 * (3 * k[1] - k[2]),
                              0.5 * (k[1] + k[2]) + dk * np.arange(nbin - 2)]
            if(logarithmic):
                self._binbounds = np.exp(self._binbounds)

Jait Dixit's avatar
Jait Dixit committed
319
320
321
322

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