power_space.py 9.6 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
# Copyright(C) 2013-2019 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Philipp Arras committed
17

18
import numpy as np
Philipp Arras's avatar
Philipp Arras committed
19

Martin Reinecke's avatar
Martin Reinecke committed
20
from .. import dobj
Philipp Arras's avatar
Philipp Arras committed
21
from .structured_domain import StructuredDomain
Theo Steininger's avatar
Theo Steininger committed
22 23


Martin Reinecke's avatar
Martin Reinecke committed
24
class PowerSpace(StructuredDomain):
25
    """Represents non-equidistantly binned spaces for power spectra.
Theo Steininger's avatar
Theo Steininger committed
26

Martin Reinecke's avatar
Martin Reinecke committed
27
    A power space is the result of a projection of a harmonic domain where
Martin Reinecke's avatar
Martin Reinecke committed
28 29
    k-modes of equal length get mapped to one power index.

Theo Steininger's avatar
Theo Steininger committed
30 31
    Parameters
    ----------
Martin Reinecke's avatar
Martin Reinecke committed
32
    harmonic_partner : StructuredDomain
Martin Reinecke's avatar
Martin Reinecke committed
33
        The harmonic domain of which this is the power space.
Philipp Arras's avatar
Philipp Arras committed
34 35 36 37 38
    binbounds : None, or tuple of float
        By default (binbounds=None):
            There are as many bins as there are distinct k-vector lengths in
            the harmonic partner space.
            The `binbounds` property of the PowerSpace will be None.
Martin Reinecke's avatar
Martin Reinecke committed
39
        else:
Philipp Arras's avatar
Philipp Arras committed
40
            The bin bounds requested for this PowerSpace. The array
Martin Reinecke's avatar
Martin Reinecke committed
41 42
            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
Martin Reinecke's avatar
Martin Reinecke committed
43 44 45
            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
46
    """
47

48
    _powerIndexCache = {}
Martin Reinecke's avatar
Martin Reinecke committed
49
    _needed_for_hash = ["_harmonic_partner", "_binbounds"]
50

Martin Reinecke's avatar
Martin Reinecke committed
51
    @staticmethod
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
52
    def linear_binbounds(nbin, first_bound, last_bound):
Martin Reinecke's avatar
Martin Reinecke committed
53 54
        """Produces linearly spaced bin bounds.

Martin Reinecke's avatar
Martin Reinecke committed
55 56
        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
57
        nbin : int
Martin Reinecke's avatar
Martin Reinecke committed
58
            the number of bins
Martin Reinecke's avatar
Martin Reinecke committed
59
        first_bound, last_bound : float
Martin Reinecke's avatar
Martin Reinecke committed
60 61 62
            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.
Martin Reinecke's avatar
Martin Reinecke committed
63 64 65 66 67 68 69 70

        Returns
        -------
        numpy.ndarray(numpy.float64)
            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.
Martin Reinecke's avatar
Martin Reinecke committed
71 72
        """
        nbin = int(nbin)
73 74
        if nbin < 3:
            raise ValueError("nbin must be at least 3")
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
75
        return np.linspace(float(first_bound), float(last_bound), nbin-1)
Martin Reinecke's avatar
Martin Reinecke committed
76 77

    @staticmethod
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
78
    def logarithmic_binbounds(nbin, first_bound, last_bound):
Martin Reinecke's avatar
Martin Reinecke committed
79 80
        """Produces logarithmically spaced bin bounds.

Martin Reinecke's avatar
Martin Reinecke committed
81 82
        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
83
        nbin : int
Martin Reinecke's avatar
Martin Reinecke committed
84
            the number of bins
Martin Reinecke's avatar
Martin Reinecke committed
85
        first_bound, last_bound : float
Martin Reinecke's avatar
Martin Reinecke committed
86 87 88
            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.
Martin Reinecke's avatar
Martin Reinecke committed
89 90 91 92 93 94 95 96

        Returns
        -------
        numpy.ndarray(numpy.float64)
            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
97
        """
Martin Reinecke's avatar
Martin Reinecke committed
98
        nbin = int(nbin)
99 100
        if nbin < 3:
            raise ValueError("nbin must be at least 3")
Martin Reinecke's avatar
Martin Reinecke committed
101 102 103
        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
104

105 106
    @staticmethod
    def useful_binbounds(space, logarithmic, nbin=None):
Martin Reinecke's avatar
Martin Reinecke committed
107 108
        """Produces bin bounds suitable for a given domain.

Martin Reinecke's avatar
Martin Reinecke committed
109 110
        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
111 112 113 114
        space : StructuredDomain
            the domain for which the binbounds will be computed.
        logarithmic : bool
            If True bins will have equal size in linear space; otherwise they
Martin Reinecke's avatar
Martin Reinecke committed
115
            will have equal size in logarithmic space.
Martin Reinecke's avatar
Martin Reinecke committed
116 117 118
        nbin : int, optional
            the number of bins
            If None, the highest possible number of bins will be used
Martin Reinecke's avatar
Martin Reinecke committed
119 120 121 122 123 124 125 126

        Returns
        -------
        numpy.ndarray(numpy.float64)
            Binbounds array with `nbin-1` entries, if `nbin` is
            supplied, or the maximum number of entries that does not produce
            empty bins, if `nbin` is not supplied.
            The first and last bin boundary are inferred from `space`.
Martin Reinecke's avatar
Martin Reinecke committed
127
        """
Martin Reinecke's avatar
Martin Reinecke committed
128
        if not (isinstance(space, StructuredDomain) and space.harmonic):
129 130 131 132 133
            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)
134
        dists = space.get_unique_k_lengths()
135 136 137 138 139 140 141 142 143 144 145 146
        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
147 148
        if nbin < 3:
            raise ValueError("nbin must be at least 3")
149 150 151 152 153 154 155
        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
156
    def __init__(self, harmonic_partner, binbounds=None):
Martin Reinecke's avatar
Martin Reinecke committed
157
        if not (isinstance(harmonic_partner, StructuredDomain) and
Martin Reinecke's avatar
Martin Reinecke committed
158 159
                harmonic_partner.harmonic):
            raise ValueError("harmonic_partner must be a harmonic space.")
Martin Reinecke's avatar
Martin Reinecke committed
160
        if harmonic_partner.scalar_dvol is None:
Martin Reinecke's avatar
Martin Reinecke committed
161 162
            raise ValueError("harmonic partner must have "
                             "scalar volume factors")
163
        self._harmonic_partner = harmonic_partner
Martin Reinecke's avatar
Martin Reinecke committed
164
        pdvol = harmonic_partner.scalar_dvol
165

Martin Reinecke's avatar
Martin Reinecke committed
166 167
        if binbounds is not None:
            binbounds = tuple(binbounds)
Philipp Arras's avatar
Fixup  
Philipp Arras committed
168
            if min(binbounds) < 0:
169
                raise ValueError('Negative binbounds encountered')
170

Martin Reinecke's avatar
Martin Reinecke committed
171
        key = (harmonic_partner, binbounds)
172
        if self._powerIndexCache.get(key) is None:
173
            k_length_array = self.harmonic_partner.get_k_length_array()
Martin Reinecke's avatar
Martin Reinecke committed
174 175 176 177 178
            if binbounds is None:
                tmp = harmonic_partner.get_unique_k_lengths()
                tbb = 0.5*(tmp[:-1]+tmp[1:])
            else:
                tbb = binbounds
Martin Reinecke's avatar
Martin Reinecke committed
179
            locdat = np.searchsorted(tbb, k_length_array.local_data)
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
180
            temp_pindex = dobj.from_local_data(
181 182
                k_length_array.val.shape, locdat,
                dobj.distaxis(k_length_array.val))
Martin Reinecke's avatar
Martin Reinecke committed
183
            nbin = len(tbb)+1
Martin Reinecke's avatar
Martin Reinecke committed
184 185
            temp_rho = np.bincount(dobj.local_data(temp_pindex).ravel(),
                                   minlength=nbin)
Martin Reinecke's avatar
Martin Reinecke committed
186
            temp_rho = dobj.np_allreduce_sum(temp_rho)
187 188
            if (temp_rho == 0).any():
                raise ValueError("empty bins detected")
Martin Reinecke's avatar
Martin Reinecke committed
189 190 191
            # The explicit conversion to float64 is necessary because bincount
            # sometimes returns its result as an integer array, even when
            # floating-point weights are present ...
192 193
            temp_k_lengths = np.bincount(
                dobj.local_data(temp_pindex).ravel(),
Martin Reinecke's avatar
Martin Reinecke committed
194
                weights=k_length_array.local_data.ravel(),
Martin Reinecke's avatar
Martin Reinecke committed
195
                minlength=nbin).astype(np.float64, copy=False)
Martin Reinecke's avatar
Martin Reinecke committed
196
            temp_k_lengths = dobj.np_allreduce_sum(temp_k_lengths) / temp_rho
Martin Reinecke's avatar
Martin Reinecke committed
197 198
            temp_k_lengths.flags.writeable = False
            dobj.lock(temp_pindex)
Martin Reinecke's avatar
Martin Reinecke committed
199
            temp_dvol = temp_rho*pdvol
Martin Reinecke's avatar
Martin Reinecke committed
200
            temp_dvol.flags.writeable = False
Martin Reinecke's avatar
Martin Reinecke committed
201 202
            self._powerIndexCache[key] = (binbounds, temp_pindex,
                                          temp_k_lengths, temp_dvol)
203

Martin Reinecke's avatar
Martin Reinecke committed
204
        (self._binbounds, self._pindex, self._k_lengths, self._dvol) = \
205 206
            self._powerIndexCache[key]

207
    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
208 209
        return ("PowerSpace(harmonic_partner={}, binbounds={})"
                .format(self.harmonic_partner, self._binbounds))
210

211 212
    @property
    def harmonic(self):
Martin Reinecke's avatar
Martin Reinecke committed
213
        """bool : Always False for this class."""
214
        return False
215

216 217
    @property
    def shape(self):
Martin Reinecke's avatar
Martin Reinecke committed
218
        return self.k_lengths.shape
219

220
    @property
Martin Reinecke's avatar
Martin Reinecke committed
221
    def size(self):
222 223
        return self.shape[0]

Martin Reinecke's avatar
Martin Reinecke committed
224
    @property
225
    def scalar_dvol(self):
Martin Reinecke's avatar
Martin Reinecke committed
226 227
        return None

Martin Reinecke's avatar
Martin Reinecke committed
228
    @property
Martin Reinecke's avatar
Martin Reinecke committed
229 230
    def dvol(self):
        return self._dvol
231

232
    @property
233
    def harmonic_partner(self):
Martin Reinecke's avatar
Martin Reinecke committed
234
        """StructuredDomain : the harmonic domain associated with `self`."""
235
        return self._harmonic_partner
236 237

    @property
Martin Reinecke's avatar
Martin Reinecke committed
238
    def binbounds(self):
Martin Reinecke's avatar
Martin Reinecke committed
239 240 241 242
        """None or tuple of float : inner bin boundaries

        The boundaries between bins, starting with the right boundary of the
        first bin, up to the left boundary of the last bin.
Martin Reinecke's avatar
Martin Reinecke committed
243 244

        `None` is used to indicate natural binning.
Martin Reinecke's avatar
Martin Reinecke committed
245
        """
Martin Reinecke's avatar
Martin Reinecke committed
246
        return self._binbounds
247 248 249

    @property
    def pindex(self):
Martin Reinecke's avatar
Martin Reinecke committed
250 251 252
        """data_object : bin indices

        Bin index for every pixel in the harmonic partner.
Theo Steininger's avatar
Theo Steininger committed
253
        """
254 255 256
        return self._pindex

    @property
Martin Reinecke's avatar
Martin Reinecke committed
257
    def k_lengths(self):
Martin Reinecke's avatar
Martin Reinecke committed
258
        """numpy.ndarray(float) : k-vector length for each bin."""
Martin Reinecke's avatar
Martin Reinecke committed
259
        return self._k_lengths