power_space.py 8.55 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
Martin Reinecke's avatar
Martin Reinecke committed
23
from ... import dobj
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
    """ 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
33
34
35
36
37
38
39
40
41
42
43
44
    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
45
46
47
48
        (default : None)

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

    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.

    """
73

74
75
    _powerIndexCache = {}

76
77
    # ---Overwritten properties and methods---

Martin Reinecke's avatar
Martin Reinecke committed
78
    @staticmethod
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
79
    def linear_binbounds(nbin, first_bound, last_bound):
Martin Reinecke's avatar
Martin Reinecke committed
80
81
82
83
84
85
86
87
88
89
90
91
        """
        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
92
93
        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
94
95

    @staticmethod
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
96
    def logarithmic_binbounds(nbin, first_bound, last_bound):
Martin Reinecke's avatar
Martin Reinecke committed
97
98
99
100
101
102
103
104
105
106
107
108
        """
        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
109
        nbin = int(nbin)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
110
        assert nbin >= 3, "nbin must be at least 3"
Martin Reinecke's avatar
Martin Reinecke committed
111
112
113
        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
114

115
116
117
118
119
120
121
122
    @staticmethod
    def useful_binbounds(space, logarithmic, nbin=None):
        if not (isinstance(space, Space) and space.harmonic):
            raise ValueError("first argument must be a harmonic space.")
        if logarithmic is None and nbin is None:
            return None
        nbin = None if nbin is None else int(nbin)
        logarithmic = bool(logarithmic)
123
        dists = space.get_unique_k_lengths()
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
        if len(dists) < 3:
            raise ValueError("Space does not have enough unique k lengths")
        lbound = 0.5*(dists[0]+dists[1])
        rbound = 0.5*(dists[-2]+dists[-1])
        dists[0] = lbound
        dists[-1] = rbound
        if logarithmic:
            dists = np.log(dists)
        binsz_min = np.max(np.diff(dists))
        nbin_max = int((dists[-1]-dists[0])/binsz_min)+2
        if nbin is None:
            nbin = nbin_max
        assert nbin >= 3, "nbin must be at least 3"
        if nbin > nbin_max:
            raise ValueError("nbin is too large")
        if logarithmic:
            return PowerSpace.logarithmic_binbounds(nbin, lbound, rbound)
        else:
            return PowerSpace.linear_binbounds(nbin, lbound, rbound)

Martin Reinecke's avatar
Martin Reinecke committed
144
    def __init__(self, harmonic_partner, binbounds=None):
Martin Reinecke's avatar
Martin Reinecke committed
145
        super(PowerSpace, self).__init__()
146
        self._needed_for_hash += ['_harmonic_partner', '_binbounds']
147

Martin Reinecke's avatar
Martin Reinecke committed
148
149
150
        if not (isinstance(harmonic_partner, Space) and
                harmonic_partner.harmonic):
            raise ValueError("harmonic_partner must be a harmonic space.")
151
        self._harmonic_partner = harmonic_partner
152

Martin Reinecke's avatar
Martin Reinecke committed
153
154
        if binbounds is not None:
            binbounds = tuple(binbounds)
155

Martin Reinecke's avatar
Martin Reinecke committed
156
        key = (harmonic_partner, binbounds)
157
        if self._powerIndexCache.get(key) is None:
158
            k_length_array = self.harmonic_partner.get_k_length_array()
159
            temp_pindex = self._compute_pindex(
160
                                harmonic_partner=self.harmonic_partner,
161
                                k_length_array=k_length_array,
Martin Reinecke's avatar
Martin Reinecke committed
162
                                binbounds=binbounds)
Martin Reinecke's avatar
Martin Reinecke committed
163
            temp_rho = dobj.bincount(temp_pindex.ravel())
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
164
            assert not np.any(temp_rho == 0), "empty bins detected"
Martin Reinecke's avatar
Martin Reinecke committed
165
            temp_k_lengths = np.bincount(temp_pindex.ravel(),
166
                                      weights=k_length_array.ravel()) \
167
                / temp_rho
Martin Reinecke's avatar
Martin Reinecke committed
168
            self._powerIndexCache[key] = (binbounds,
169
                                          temp_pindex,
Martin Reinecke's avatar
Martin Reinecke committed
170
                                          temp_k_lengths,
171
172
                                          temp_rho)

Martin Reinecke's avatar
Martin Reinecke committed
173
        (self._binbounds, self._pindex, self._k_lengths, self._rho) = \
174
175
            self._powerIndexCache[key]

176
    @staticmethod
177
    def _compute_pindex(harmonic_partner, k_length_array, binbounds):
178
        if binbounds is None:
179
            tmp = harmonic_partner.get_unique_k_lengths()
180
            binbounds = 0.5*(tmp[:-1]+tmp[1:])
181
        return np.searchsorted(binbounds, k_length_array)
182

183
184
    # ---Mandatory properties and methods---

185
    def __repr__(self):
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
186
187
        return ("PowerSpace(harmonic_partner=%r, binbounds=%r)"
                % (self.harmonic_partner, self._binbounds))
188

189
190
191
    @property
    def harmonic(self):
        return True
192

193
194
    @property
    def shape(self):
Martin Reinecke's avatar
Martin Reinecke committed
195
        return self.k_lengths.shape
196

197
198
199
200
    @property
    def dim(self):
        return self.shape[0]

201
    def scalar_dvol(self):
202
        return 1.
203

204
    def get_k_length_array(self):
Martin Reinecke's avatar
Martin Reinecke committed
205
        return self.k_lengths.copy()
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
Martin Reinecke's avatar
Martin Reinecke committed
220
221
    def binbounds(self):
        return self._binbounds
222
223
224

    @property
    def pindex(self):
Martin Reinecke's avatar
Martin Reinecke committed
225
        """ A data object having the shape of the harmonic partner
Theo Steininger's avatar
Theo Steininger committed
226
227
        space containing the indices of the power bin a pixel belongs to.
        """
228
229
230
        return self._pindex

    @property
Martin Reinecke's avatar
Martin Reinecke committed
231
    def k_lengths(self):
Theo Steininger's avatar
Theo Steininger committed
232
233
        """ Sorted array of all k-modes.
        """
Martin Reinecke's avatar
Martin Reinecke committed
234
        return self._k_lengths
235
236
237

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