gl_space.py 7.11 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 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/>.

csongor's avatar
csongor committed
19
20
from __future__ import division

Jait Dixit's avatar
Jait Dixit committed
21
import itertools
csongor's avatar
csongor committed
22
23
import numpy as np

24
from nifty.spaces.space import Space
25
from nifty.config import dependency_injector as gdi
26

Theo Steininger's avatar
Theo Steininger committed
27
28
29
pyHealpix = gdi.get('pyHealpix')


Theo Steininger's avatar
Theo Steininger committed
30
class GLSpace(Space):
csongor's avatar
csongor committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
    """
        ..                 __
        ..               /  /
        ..     ____ __  /  /
        ..   /   _   / /  /
        ..  /  /_/  / /  /_
        ..  \___   /  \___/  space class
        .. /______/

        NIFTY subclass for Gauss-Legendre pixelizations [#]_ of the two-sphere.

        Parameters
        ----------
        nlat : int
            Number of latitudinal bins, or rings.
46
47
        nlon : int, *optional*
            Number of longitudinal bins (default: ``2*nlat - 1``).
csongor's avatar
csongor committed
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

        See Also
        --------
        hp_space : A class for the HEALPix discretization of the sphere [#]_.
        lm_space : A class for spherical harmonic components.

        References
        ----------
        .. [#] M. Reinecke and D. Sverre Seljebotn, 2013, "Libsharp - spherical
               harmonic transforms revisited";
               `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
        .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
               High-Resolution Discretization and Fast Analysis of Data
               Distributed on the Sphere", *ApJ* 622..759G.

    """

65
66
    # ---Overwritten properties and methods---

Martin Reinecke's avatar
Martin Reinecke committed
67
    def __init__(self, nlat, nlon=None):
csongor's avatar
csongor committed
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
        """
            Sets the attributes for a gl_space class instance.

            Parameters
            ----------
            nlat : int
                Number of latitudinal bins, or rings.
            nlon : int, *optional*
                Number of longitudinal bins (default: ``2*nlat - 1``).

            Returns
            -------
            None

            Raises
            ------
            ValueError
85
                If input `nlat` or `nlon` is invalid.
csongor's avatar
csongor committed
86
87

        """
Theo Steininger's avatar
Theo Steininger committed
88
89
90
        if 'pyHealpix' not in gdi:
            raise ImportError(
                "The module pyHealpix is needed but not available.")
91

Martin Reinecke's avatar
Martin Reinecke committed
92
        super(GLSpace, self).__init__()
csongor's avatar
csongor committed
93

94
95
        self._nlat = self._parse_nlat(nlat)
        self._nlon = self._parse_nlon(nlon)
csongor's avatar
csongor committed
96

97
    # ---Mandatory properties and methods---
csongor's avatar
csongor committed
98

99
100
    @property
    def harmonic(self):
101
102
103
104
105
106
        """True if this can be regarded as a harmonic space.
        Returns
        -------
        bool : False
            Always returns False as the GLSpace cannot be regarded as harmonic space.
        """
107
        return False
csongor's avatar
csongor committed
108
109
110

    @property
    def shape(self):
111
        return (np.int((self.nlat * self.nlon)),)
csongor's avatar
csongor committed
112

113
    @property
114
    def dim(self):
115
        return np.int((self.nlat * self.nlon))
116
117
118
119

    @property
    def total_volume(self):
        return 4 * np.pi
120

121
122
    def copy(self):
        return self.__class__(nlat=self.nlat,
Martin Reinecke's avatar
Martin Reinecke committed
123
                              nlon=self.nlon)
124

Jait Dixit's avatar
Jait Dixit committed
125
    def weight(self, x, power=1, axes=None, inplace=False):
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
        """ Weights a field living on this space with a specified amount of volume-weights.

        Weights hereby refer to integration weights, as they appear in discretized integrals.
        Per default, this function mutliplies each bin of the field x by its volume, which lets
        it behave like a density (top form). However, different powers of the volume can be applied
        with the power parameter. If only certain axes are specified via the axes parameter,
        the weights are only applied with respect to these dimensions, yielding an object that
        behaves like a lower degree form.
        Parameters
        ----------
        x : Field
            A field with this space as domain to be weighted.
        power : int, *optional*
            The power to which the volume-weight is raised.
            (default: 1).
        axes : {int, tuple}, *optional*
            Specifies for which axes the weights should be applied.
            (default: None).
            If axes==None:
                weighting is applied with respect to all axes
        inplace : bool, *optional*
            If this is True, the weighting is done on the values of x,
            if it is False, x is not modified and this method returns a 
            weighted copy of x
            (default: False).

        Returns
        -------
        Field
            A weighted version of x, with volume-weights raised to power.
            
        """
        
159
160
        nlon = self.nlon
        nlat = self.nlat
Theo Steininger's avatar
Theo Steininger committed
161
        vol = pyHealpix.GL_weights(nlat, nlon) ** power
162
        weight = np.array(list(itertools.chain.from_iterable(
Theo Steininger's avatar
Theo Steininger committed
163
                          itertools.repeat(x, nlon) for x in vol)))
Jait Dixit's avatar
Jait Dixit committed
164
165
166

        if axes is not None:
            # reshape the weight array to match the input shape
167
            new_shape = np.ones(len(x.shape), dtype=np.int)
168
169
            # we know len(axes) is always 1
            new_shape[axes[0]] = len(weight)
Jait Dixit's avatar
Jait Dixit committed
170
171
172
173
174
            weight = weight.reshape(new_shape)

        if inplace:
            x *= weight
            result_x = x
csongor's avatar
csongor committed
175
        else:
Jait Dixit's avatar
Jait Dixit committed
176
            result_x = x * weight
csongor's avatar
csongor committed
177

Jait Dixit's avatar
Jait Dixit committed
178
        return result_x
179

180
    def get_distance_array(self, distribution_strategy):
181
182
183
184
185
186
187
        """This should not be used, it just raises an error when called.
        
        Raises
        ------
        NotImplementedError
            Always when called.
        """
Theo Steininger's avatar
Theo Steininger committed
188
        raise NotImplementedError
189

190
    def get_fft_smoothing_kernel_function(self, sigma):
191
192
193
194
195
196
197
        """This should not be used, it just raises an error when called.
        
        Raises
        ------
        NotImplementedError
            Always when called.
        """
Theo Steininger's avatar
Theo Steininger committed
198
        raise NotImplementedError
199

200
201
202
203
204
205
206
207
208
209
210
211
    # ---Added properties and methods---

    @property
    def nlat(self):
        return self._nlat

    @property
    def nlon(self):
        return self._nlon

    def _parse_nlat(self, nlat):
        nlat = int(nlat)
212
213
214
        if nlat < 1:
            raise ValueError(
                "nlat must be a positive number.")
215
216
217
218
219
220
221
        return nlat

    def _parse_nlon(self, nlon):
        if nlon is None:
            nlon = 2 * self.nlat - 1
        else:
            nlon = int(nlon)
222
223
            if nlon < 1:
                raise ValueError("nlon must be a positive number.")
224
        return nlon
225
226
227
228

    # ---Serialization---

    def _to_hdf5(self, hdf5_group):
Jait Dixit's avatar
Jait Dixit committed
229
230
231
        hdf5_group['nlat'] = self.nlat
        hdf5_group['nlon'] = self.nlon

232
233
234
        return None

    @classmethod
Theo Steininger's avatar
Theo Steininger committed
235
    def _from_hdf5(cls, hdf5_group, repository):
236
        result = cls(
Jait Dixit's avatar
Jait Dixit committed
237
238
239
240
            nlat=hdf5_group['nlat'][()],
            nlon=hdf5_group['nlon'][()],
            )

241
        return result