power_space.py 10.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# NIFTy
# Copyright (C) 2017  Theo Steininger
#
# Author: Theo Steininger
#
# 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/>.
theos's avatar
theos committed
18

theos's avatar
theos committed
19
20
import numpy as np

21
22
import d2o

23
24
from power_index_factory import PowerIndexFactory

25
from nifty.spaces.space import Space
26
from nifty.spaces.rg_space import RGSpace
theos's avatar
theos committed
27
28


Theo Steininger's avatar
Theo Steininger committed
29
class PowerSpace(Space):
30

31
32
    # ---Overwritten properties and methods---

33
    def __init__(self, harmonic_partner=RGSpace((1,)),
34
                 distribution_strategy='not',
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
                 logarithmic=False, nbin=None, binbounds=None):
        """Sets the attributes for a PowerSpace class instance.

            Parameters
            ----------
            harmonic_partner : Space
                The harmonic Space of which this is the power space.
            distribution_strategy : str *optional*
                The distribution strategy of a d2o-object represeting a field over this PowerSpace.
                (default : 'not')
            logarithmic : bool *optional*
                True if logarithmic binning should be used.
                (default : False)
            nbin : {int, None} *optional*
                The number of bins this space has.
                (default : None) if nbin == None : It takes the nbin from its harmonic_partner
            binbounds :  {list, array} *optional*
                Array-like inner boundaries of the used bins of the default
                indices.
                (default : None) if binbounds == None : Calculates the bounds from the kindex and corrects for logartihmic scale 
            Notes
            -----
            A power space is the result of a projection of a harmonic space where multiple k-modes get mapped to one power index.
            This can be regarded as a response operator :math:`R` going from harmonic space to power space. 
            An array giving this map is stored in pindex (array which says in which power box a k-mode gets projected)
            An array for the adjoint of :math:`R` is given by kindex, which is an array of arrays stating which k-mode got mapped to a power index
            The a right-inverse to :math:`R` is given by the pundex which is an array giving one k-mode that maps to a power bin for every power bin.
            The amount of k-modes that get mapped to one power bin is given by rho. This is :math:`RR^\dagger` in the language of this projection operator            
            Returns
            -------
            None.

        """
        #FIXME: default probably not working for log and normal scale
Martin Reinecke's avatar
Martin Reinecke committed
69
        super(PowerSpace, self).__init__()
70
71
        self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
                                  '_k_array']
72

73
        if not isinstance(harmonic_partner, Space):
74
            raise ValueError(
75
76
                "harmonic_partner must be a Space.")
        if not harmonic_partner.harmonic:
77
            raise ValueError(
78
79
                "harmonic_partner must be a harmonic space.")
        self._harmonic_partner = harmonic_partner
80

Jait Dixit's avatar
Jait Dixit committed
81
        power_index = PowerIndexFactory.get_power_index(
82
                        domain=self.harmonic_partner,
Jait Dixit's avatar
Jait Dixit committed
83
                        distribution_strategy=distribution_strategy,
84
                        logarithmic=logarithmic,
Jait Dixit's avatar
Jait Dixit committed
85
86
                        nbin=nbin,
                        binbounds=binbounds)
87
88

        config = power_index['config']
89
        self._logarithmic = config['logarithmic']
90
91
92
93
94
95
        self._nbin = config['nbin']
        self._binbounds = config['binbounds']

        self._pindex = power_index['pindex']
        self._kindex = power_index['kindex']
        self._rho = power_index['rho']
96
97
        self._pundex = power_index['pundex']
        self._k_array = power_index['k_array']
98

Theo Steininger's avatar
Theo Steininger committed
99
    def pre_cast(self, x, axes):
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
        """Casts power spectra to discretized power spectra.
        
        This function takes an array or a function. If it is an array it does nothing,
        otherwise it intepretes the function as power spectrum and evaluates it at every
        k-mode.
        Parameters
        ----------
        x : {array-like, function array-like -> array-like}
            power spectrum given either in discretized form or implicitly as a function
        axes : {tuple, int} *optional*
            does nothing
            (default : None)
        Returns
        -------
        array-like : discretized power spectrum
        """
116
117
118
119
120
        if callable(x):
            return x(self.kindex)
        else:
            return x

121
122
123
124
125
    # ---Mandatory properties and methods---

    @property
    def harmonic(self):
        return True
126

127
128
    @property
    def shape(self):
129
        return self.kindex.shape
130

131
132
133
134
135
136
137
    @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
138
        return float(reduce(lambda x, y: x*y, self.pindex.shape))
139
140

    def copy(self):
141
        distribution_strategy = self.pindex.distribution_strategy
142
        return self.__class__(harmonic_partner=self.harmonic_partner,
143
                              distribution_strategy=distribution_strategy,
144
                              logarithmic=self.logarithmic,
145
                              nbin=self.nbin,
Martin Reinecke's avatar
Martin Reinecke committed
146
                              binbounds=self.binbounds)
147

148
    def weight(self, x, power=1, axes=None, inplace=False):
Jait Dixit's avatar
Jait Dixit committed
149
150
        reshaper = [1, ] * len(x.shape)
        # we know len(axes) is always 1
151
152
        reshaper[axes[0]] = self.shape[0]

153
        weight = self.rho.reshape(reshaper)
154
155
        if power != 1:
            weight = weight ** power
156
157
158
159
160
161

        if inplace:
            x *= weight
            result_x = x
        else:
            result_x = x*weight
162
163
164

        return result_x

165
    def get_distance_array(self, distribution_strategy):
166
        result = d2o.distributed_data_object(
Martin Reinecke's avatar
Martin Reinecke committed
167
                                self.kindex, dtype=np.float64,
168
169
                                distribution_strategy=distribution_strategy)
        return result
theos's avatar
theos committed
170

171
    def get_fft_smoothing_kernel_function(self, sigma):
172
        raise NotImplementedError(
173
            "There is no fft smoothing function for PowerSpace.")
theos's avatar
theos committed
174

175
176
177
    # ---Added properties and methods---

    @property
178
179
180
181
182
183
184
    def harmonic_partner(self):
        """Returns the Space of which this is the power space.
        Returns
        -------
        Space : The harmonic Space of which this is the power space.
        """
        return self._harmonic_partner
185
186

    @property
187
188
189
190
191
192
193
    def logarithmic(self):
        """Returns a True if logarithmic binning is used.
        Returns
        -------
        Bool : True if for this PowerSpace logarithmic binning is used.
        """
        return self._logarithmic
194
195
196

    @property
    def nbin(self):
197
198
199
200
201
        """Returns the number of power bins.
        Returns
        -------
        int : The number of bins this space has.
        """
202
203
204
205
        return self._nbin

    @property
    def binbounds(self):
206
207
208
209
210
211
212
213
        """ Inner boundaries of the used bins of the default
                indices.
        Returns
        -------
        {list, array} : the inner boundaries of the used bins in the used scale, as they were
        set in __init__ or computed.
        """
        # FIXME check wether this returns something sensible if 'None' was set in __init__
214
215
216
217
        return self._binbounds

    @property
    def pindex(self):
218
219
220
221
222
    """Index of the Fourier grid points that belong to a specific power index
    Returns
    -------
        distributed_data_object : Index of the Fourier grid points in a distributed_data_object.
    """
223
224
225
226
        return self._pindex

    @property
    def kindex(self):
227
228
229
230
231
    """Array of all k-vector lengths.
    Returns
    -------
        ndarray : Array which states for each k-mode which power index it maps to (adjoint to pindex)
    """
232
233
234
235
        return self._kindex

    @property
    def rho(self):
236
237
238
239
    """Degeneracy factor of the individual k-vectors.
    
    ndarray : Array stating how many k-modes are mapped to one power index for every power index
    """
240
        return self._rho
241

242
243
    @property
    def pundex(self):
244
245
246
247
248
    """List of one k-mode per power bin which is in the bin.
    Returns
    -------
    array-like : An array for which the n-th entry is an example one k-mode which belongs to the n-th power bin
    """
249
250
251
252
        return self._pundex

    @property
    def k_array(self):
253
254
255
256
257
258
    """This contains distances to zero for every k-mode of the harmonic partner.
    
    Returns
    -------
    array-like : An array containing distances to the zero mode for every k-mode of the harmonic partner.
    """
259
        return self._k_array
260
261
262
263

    # ---Serialization---

    def _to_hdf5(self, hdf5_group):
Jait Dixit's avatar
Jait Dixit committed
264
265
266
        hdf5_group['kindex'] = self.kindex
        hdf5_group['rho'] = self.rho
        hdf5_group['pundex'] = self.pundex
267
        hdf5_group['logarithmic'] = self.logarithmic
Theo Steininger's avatar
Theo Steininger committed
268
269
270
        # Store nbin as string, since it can be None
        hdf5_group.attrs['nbin'] = str(self.nbin)
        hdf5_group.attrs['binbounds'] = str(self.binbounds)
271
272

        return {
273
            'harmonic_partner': self.harmonic_partner,
274
275
276
277
278
            'pindex': self.pindex,
            'k_array': self.k_array
        }

    @classmethod
Theo Steininger's avatar
Theo Steininger committed
279
    def _from_hdf5(cls, hdf5_group, repository):
Jait Dixit's avatar
Jait Dixit committed
280
281
282
283
        # make an empty PowerSpace object
        new_ps = EmptyPowerSpace()
        # reset class
        new_ps.__class__ = cls
Jait Dixit's avatar
Jait Dixit committed
284
        # call instructor so that classes are properly setup
Martin Reinecke's avatar
Martin Reinecke committed
285
        super(PowerSpace, new_ps).__init__()
Jait Dixit's avatar
Jait Dixit committed
286
        # set all values
287
288
        new_ps._harmonic_partner = repository.get('harmonic_partner', hdf5_group)
        new_ps._logarithmic = hdf5_group['logarithmic'][()]
Theo Steininger's avatar
Theo Steininger committed
289
290
        exec('new_ps._nbin = ' + hdf5_group.attrs['nbin'])
        exec('new_ps._binbounds = ' + hdf5_group.attrs['binbounds'])
Jait Dixit's avatar
Jait Dixit committed
291

Theo Steininger's avatar
Theo Steininger committed
292
        new_ps._pindex = repository.get('pindex', hdf5_group)
Jait Dixit's avatar
Jait Dixit committed
293
294
295
        new_ps._kindex = hdf5_group['kindex'][:]
        new_ps._rho = hdf5_group['rho'][:]
        new_ps._pundex = hdf5_group['pundex'][:]
Theo Steininger's avatar
Theo Steininger committed
296
        new_ps._k_array = repository.get('k_array', hdf5_group)
Jait Dixit's avatar
Jait Dixit committed
297
        new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
298
                                    '_k_array']
Jait Dixit's avatar
Jait Dixit committed
299
300
301
302
303
304
305

        return new_ps


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