power_space.py 8.05 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 ...spaces.space import Space
Martin Reinecke's avatar
Martin Reinecke committed
22
from functools import reduce
theos's avatar
theos committed
23
24


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

    Parameters
    ----------
    harmonic_partner : Space
        The harmonic Space of which this is the power space.
Martin Reinecke's avatar
Martin Reinecke committed
32
33
34
35
36
37
38
39
40
41
42
43
    binbounds: None, or tuple/array/list of float
        if None:
            There will be as many bins as there are distinct k-vector lengths
            in the harmonic partner space.
            The "binbounds" property of the PowerSpace will also be None.

        else:
            the bin bounds requested for this PowerSpace. The array
            must be sorted and strictly ascending. The first entry is the right
            boundary of the first bin, and the last entry is the left boundary
            of the last bin, i.e. thee will be len(binbounds)+1 bins in total,
            with the first and last bins reaching to -+infinity, respectively.
Theo Steininger's avatar
Theo Steininger committed
44
45
46
47
        (default : None)

    Attributes
    ----------
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
48
    pindex : numpy.ndarray
49
50
        This holds the information which pixel of the harmonic partner gets
        mapped to which power bin
Theo Steininger's avatar
Theo Steininger committed
51
    kindex : numpy.ndarray
52
        Sorted array of all k-modes.
Theo Steininger's avatar
Theo Steininger committed
53
54
55
    rho : numpy.ndarray
        The amount of k-modes that get mapped to one power bin is given by
        rho.
56
57
58
    dim : np.int
        Total number of dimensionality, i.e. the number of pixels.
    harmonic : bool
Martin Reinecke's avatar
Martin Reinecke committed
59
        Always True for this space.
60
61
62
63
    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
64
65
66
    binbounds : tuple or None
        Boundaries between the power spectrum bins; None is used to indicate
        natural binning
Theo Steininger's avatar
Theo Steininger committed
67
68
69
70
71
72
73

    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.

    """
74

75
76
    _powerIndexCache = {}

77
78
    # ---Overwritten properties and methods---

Martin Reinecke's avatar
Martin Reinecke committed
79
    @staticmethod
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
80
    def linear_binbounds(nbin, first_bound, last_bound):
Martin Reinecke's avatar
Martin Reinecke committed
81
82
83
84
85
86
87
88
89
90
91
92
        """
        nbin: integer
            the number of bins
        first_bound, last_bound: float
            the k values for the right boundary of the first bin and the left
            boundary of the last bin, respectively. They are given in length
            units of the harmonic partner space.
        This will produce a binbounds array with nbin-1 entries with
        binbounds[0]=first_bound and binbounds[-1]=last_bound and the remaining
        values equidistantly spaced (in linear scale) between these two.
        """
        nbin = int(nbin)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
93
94
        assert nbin >= 3, "nbin must be at least 3"
        return np.linspace(float(first_bound), float(last_bound), nbin-1)
Martin Reinecke's avatar
Martin Reinecke committed
95
96

    @staticmethod
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
97
    def logarithmic_binbounds(nbin, first_bound, last_bound):
Martin Reinecke's avatar
Martin Reinecke committed
98
99
100
101
102
103
104
105
106
107
108
109
        """
        nbin: integer
            the number of bins
        first_bound, last_bound: float
            the k values for the right boundary of the first bin and the left
            boundary of the last bin, respectively. They are given in length
            units of the harmonic partner space.
        This will produce a binbounds array with nbin-1 entries with
        binbounds[0]=first_bound and binbounds[-1]=last_bound and the remaining
        values equidistantly spaced (in natural logarithmic scale)
        between these two.
        """
Martin Reinecke's avatar
Martin Reinecke committed
110
        nbin = int(nbin)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
111
        assert nbin >= 3, "nbin must be at least 3"
Martin Reinecke's avatar
Martin Reinecke committed
112
113
114
        return np.logspace(np.log(float(first_bound)),
                           np.log(float(last_bound)),
                           nbin-1, base=np.e)
Martin Reinecke's avatar
Martin Reinecke committed
115

Martin Reinecke's avatar
Martin Reinecke committed
116
    def __init__(self, harmonic_partner, binbounds=None):
Martin Reinecke's avatar
Martin Reinecke committed
117
        super(PowerSpace, self).__init__()
Martin Reinecke's avatar
Martin Reinecke committed
118
        self._ignore_for_hash += ['_pindex', '_kindex', '_rho']
119

Martin Reinecke's avatar
Martin Reinecke committed
120
121
122
        if not (isinstance(harmonic_partner, Space) and
                harmonic_partner.harmonic):
            raise ValueError("harmonic_partner must be a harmonic space.")
123
        self._harmonic_partner = harmonic_partner
124

Martin Reinecke's avatar
Martin Reinecke committed
125
126
        if binbounds is not None:
            binbounds = tuple(binbounds)
127

Martin Reinecke's avatar
Martin Reinecke committed
128
        key = (harmonic_partner, binbounds)
129
130
        if self._powerIndexCache.get(key) is None:
            distance_array = \
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
131
                self.harmonic_partner.get_distance_array()
132
            temp_pindex = self._compute_pindex(
133
                                harmonic_partner=self.harmonic_partner,
134
                                distance_array=distance_array,
Martin Reinecke's avatar
Martin Reinecke committed
135
                                binbounds=binbounds)
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
136
            temp_rho = np.bincount(temp_pindex.flatten())
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
137
            assert not np.any(temp_rho == 0), "empty bins detected"
Martin Reinecke's avatar
Martin Reinecke committed
138
            temp_kindex = np.bincount(temp_pindex.flatten(),
Martin Reinecke's avatar
Martin Reinecke committed
139
                                   weights=distance_array.flatten()) / temp_rho
Martin Reinecke's avatar
Martin Reinecke committed
140
            self._powerIndexCache[key] = (binbounds,
141
142
143
144
145
146
147
                                          temp_pindex,
                                          temp_kindex,
                                          temp_rho)

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

148
    @staticmethod
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
149
    def _compute_pindex(harmonic_partner, distance_array, binbounds):
150
        if binbounds is None:
151
            binbounds = harmonic_partner.get_natural_binbounds()
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
152
        return np.searchsorted(binbounds, distance_array)
153

154
155
    # ---Mandatory properties and methods---

156
    def __repr__(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
157
158
        return ("PowerSpace(harmonic_partner=%r, binbounds=%r)"
                % (self.harmonic_partner, self._binbounds))
159

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

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

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

    def copy(self):
178
        return self.__class__(harmonic_partner=self.harmonic_partner,
Martin Reinecke's avatar
Martin Reinecke committed
179
                              binbounds=self._binbounds)
180

Martin Reinecke's avatar
Martin Reinecke committed
181
    def weight(self, x, power, axes, inplace=False):
Jait Dixit's avatar
Jait Dixit committed
182
183
        reshaper = [1, ] * len(x.shape)
        # we know len(axes) is always 1
184
185
        reshaper[axes[0]] = self.shape[0]

186
        weight = self.rho.reshape(reshaper)
187
        if power != 1:
188
            weight = weight ** np.float(power)
189
190
191
192
193
194

        if inplace:
            x *= weight
            result_x = x
        else:
            result_x = x*weight
195
196
197

        return result_x

Martin Reinecke's avatar
stage1    
Martin Reinecke committed
198
199
    def get_distance_array(self):
        return self.kindex.copy()
theos's avatar
theos committed
200

201
    def get_fft_smoothing_kernel_function(self, sigma):
202
        raise NotImplementedError(
203
            "There is no fft smoothing function for PowerSpace.")
theos's avatar
theos committed
204

205
206
207
    # ---Added properties and methods---

    @property
208
    def harmonic_partner(self):
Theo Steininger's avatar
Theo Steininger committed
209
        """ Returns the Space of which this is the power space.
210
211
        """
        return self._harmonic_partner
212
213

    @property
Martin Reinecke's avatar
Martin Reinecke committed
214
215
    def binbounds(self):
        return self._binbounds
216
217
218

    @property
    def pindex(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
219
        """ A numpy.ndarray having the shape of the harmonic partner
Theo Steininger's avatar
Theo Steininger committed
220
221
        space containing the indices of the power bin a pixel belongs to.
        """
222
223
224
225
        return self._pindex

    @property
    def kindex(self):
Theo Steininger's avatar
Theo Steininger committed
226
227
        """ Sorted array of all k-modes.
        """
228
229
230
231
        return self._kindex

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