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

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

Martin Reinecke's avatar
Martin Reinecke committed
23
from ...spaces.space import Space
Martin Reinecke's avatar
Martin Reinecke committed
24
from functools import reduce
theos's avatar
theos committed
25
26


Theo Steininger's avatar
Theo Steininger committed
27
class PowerSpace(Space):
Theo Steininger's avatar
Theo Steininger committed
28
29
30
31
32
33
34
    """ NIFTY class for spaces of power spectra.

    Parameters
    ----------
    harmonic_partner : Space
        The harmonic Space of which this is the power space.
    logarithmic : bool *optional*
Martin Reinecke's avatar
Martin Reinecke committed
35
        True if logarithmic binning should be used (default : None).
Theo Steininger's avatar
Theo Steininger committed
36
37
38
39
40
    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
41
42
43
        Boundaries between the power spectrum bins.
        (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
44
        (default : None)
Martin Reinecke's avatar
Martin Reinecke committed
45
        if binbounds == None:
Theo Steininger's avatar
Theo Steininger committed
46
47
            Calculates the bounds from the kindex while applying the
            logarithmic and nbin keywords.
Martin Reinecke's avatar
Martin Reinecke committed
48
49
50
51
52
    Note: if "bindounds" is not None, both "logarithmic" and "nbin" must be
        None!
    Note: if "binbounds", "logarithmic", and "nbin" are all None, then
        "natural" binning is performed, i.e. there will be one bin for every
        distinct k-vector length.
Theo Steininger's avatar
Theo Steininger committed
53
54
55

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

    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.

    """
81

82
83
    _powerIndexCache = {}

84
85
    # ---Overwritten properties and methods---

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
86
    def __init__(self, harmonic_partner,
Martin Reinecke's avatar
Martin Reinecke committed
87
                 logarithmic=None, nbin=None, binbounds=None):
Martin Reinecke's avatar
Martin Reinecke committed
88
        super(PowerSpace, self).__init__()
Martin Reinecke's avatar
Martin Reinecke committed
89
        self._ignore_for_hash += ['_pindex', '_kindex', '_rho']
90

Martin Reinecke's avatar
Martin Reinecke committed
91
92
93
        if not (isinstance(harmonic_partner, Space) and
                harmonic_partner.harmonic):
            raise ValueError("harmonic_partner must be a harmonic space.")
94
        self._harmonic_partner = harmonic_partner
95

Martin Reinecke's avatar
Martin Reinecke committed
96
97
98
99
        # sanity check
        if binbounds is not None and not(nbin is None and logarithmic is None):
            raise ValueError(
                "if binbounds is defined, nbin and logarithmic must be None")
100

Martin Reinecke's avatar
Martin Reinecke committed
101
102
        if binbounds is not None:
            binbounds = tuple(binbounds)
103

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
104
        key = (harmonic_partner, logarithmic, nbin, binbounds)
105
106
        if self._powerIndexCache.get(key) is None:
            distance_array = \
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
107
                self.harmonic_partner.get_distance_array()
108
109
110
111
112
113
            temp_binbounds = self._compute_binbounds(
                                  harmonic_partner=self.harmonic_partner,
                                  logarithmic=logarithmic,
                                  nbin=nbin,
                                  binbounds=binbounds)
            temp_pindex = self._compute_pindex(
114
                                harmonic_partner=self.harmonic_partner,
115
                                distance_array=distance_array,
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
116
117
118
119
120
121
122
                                binbounds=temp_binbounds)
            #temp_rho = temp_pindex.bincount().get_full_data()
            temp_rho = np.bincount(temp_pindex.flatten())
            #temp_kindex = \
            #    (temp_pindex.bincount(weights=distance_array).get_full_data() /
            #     temp_rho)
            temp_kindex = np.bincount(temp_pindex.flatten(),weights=distance_array.flatten())/temp_rho
123
124
125
126
127
128
129
130
            self._powerIndexCache[key] = (temp_binbounds,
                                          temp_pindex,
                                          temp_kindex,
                                          temp_rho)

        (self._binbounds, self._pindex, self._kindex, self._rho) = \
            self._powerIndexCache[key]

131
    @staticmethod
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
132
    def _compute_binbounds(harmonic_partner,
133
                           logarithmic, nbin, binbounds):
134

Martin Reinecke's avatar
Martin Reinecke committed
135
        if logarithmic is None and nbin is None and binbounds is None:
136
            result = None
Martin Reinecke's avatar
Martin Reinecke committed
137
138
139
140
141
142
143
144
145
146
147
        else:
            if binbounds is not None:
                bb = np.sort(np.array(binbounds))
            else:
                if logarithmic is not None:
                    logarithmic = bool(logarithmic)
                if nbin is not None:
                    nbin = int(nbin)

                # equidistant binning (linear or log)
                # MR FIXME: this needs to improve
148
                kindex = harmonic_partner.get_unique_distances()
Martin Reinecke's avatar
Martin Reinecke committed
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
                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)
                bb = np.r_[0.5 * (3 * k[1] - k[2]),
                           0.5 * (k[1] + k[2]) + dk * np.arange(nbin-2)]
                if(logarithmic):
                    bb = np.exp(bb)
165
166
            result = tuple(bb)
        return result
167

168
    @staticmethod
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
169
    def _compute_pindex(harmonic_partner, distance_array, binbounds):
170
        if binbounds is None:
171
            binbounds = harmonic_partner.get_natural_binbounds()
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
172
        return np.searchsorted(binbounds, distance_array)
173

Theo Steininger's avatar
Theo Steininger committed
174
    def pre_cast(self, x, axes):
Theo Steininger's avatar
Theo Steininger committed
175
176
177
178
179
180
        """ 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.

181
182
183
        Parameters
        ----------
        x : {array-like, function array-like -> array-like}
Theo Steininger's avatar
Theo Steininger committed
184
185
186
187
188
189
            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.

190
191
        Returns
        -------
Theo Steininger's avatar
Theo Steininger committed
192
193
194
        array-like
            discretized power spectrum

195
        """
Theo Steininger's avatar
Theo Steininger committed
196

Martin Reinecke's avatar
Martin Reinecke committed
197
        return x(self.kindex) if callable(x) else x
198

199
200
    # ---Mandatory properties and methods---

201
    def __repr__(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
202
203
        return ("PowerSpace(harmonic_partner=%r, binbounds=%r)"
                % (self.harmonic_partner, self._binbounds))
204

205
206
207
    @property
    def harmonic(self):
        return True
208

209
210
    @property
    def shape(self):
211
        return self.kindex.shape
212

213
214
215
216
217
218
219
    @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
220
        return float(reduce(lambda x, y: x*y, self.pindex.shape))
221
222

    def copy(self):
223
        return self.__class__(harmonic_partner=self.harmonic_partner,
Martin Reinecke's avatar
Martin Reinecke committed
224
                              binbounds=self._binbounds)
225

Martin Reinecke's avatar
Martin Reinecke committed
226
    def weight(self, x, power, axes, inplace=False):
Jait Dixit's avatar
Jait Dixit committed
227
228
        reshaper = [1, ] * len(x.shape)
        # we know len(axes) is always 1
229
230
        reshaper[axes[0]] = self.shape[0]

231
        weight = self.rho.reshape(reshaper)
232
        if power != 1:
233
            weight = weight ** np.float(power)
234
235
236
237
238
239

        if inplace:
            x *= weight
            result_x = x
        else:
            result_x = x*weight
240
241
242

        return result_x

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
243
244
    def get_distance_array(self):
        return self.kindex.copy()
theos's avatar
theos committed
245

246
    def get_fft_smoothing_kernel_function(self, sigma):
247
        raise NotImplementedError(
248
            "There is no fft smoothing function for PowerSpace.")
theos's avatar
theos committed
249

250
251
252
    # ---Added properties and methods---

    @property
253
    def harmonic_partner(self):
Theo Steininger's avatar
Theo Steininger committed
254
        """ Returns the Space of which this is the power space.
255
256
        """
        return self._harmonic_partner
257
258

    @property
Martin Reinecke's avatar
Martin Reinecke committed
259
260
    def binbounds(self):
        return self._binbounds
261
262
263

    @property
    def pindex(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
264
        """ A numpy.ndarray having the shape of the harmonic partner
Theo Steininger's avatar
Theo Steininger committed
265
266
        space containing the indices of the power bin a pixel belongs to.
        """
267
268
269
270
        return self._pindex

    @property
    def kindex(self):
Theo Steininger's avatar
Theo Steininger committed
271
272
        """ Sorted array of all k-modes.
        """
273
274
275
276
        return self._kindex

    @property
    def rho(self):
Theo Steininger's avatar
Theo Steininger committed
277
278
        """Degeneracy factor of the individual k-vectors.
        """
279
        return self._rho