lm_space.py 4.93 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
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2018 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15
16
17
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

19
from __future__ import absolute_import, division, print_function
Philipp Arras's avatar
Philipp Arras committed
20

csongor's avatar
csongor committed
21
import numpy as np
Philipp Arras's avatar
Philipp Arras committed
22
23

from ..compat import *
24
from ..field import Field
Philipp Arras's avatar
Philipp Arras committed
25
from .structured_domain import StructuredDomain
Theo Steininger's avatar
Theo Steininger committed
26

csongor's avatar
csongor committed
27

Martin Reinecke's avatar
Martin Reinecke committed
28
class LMSpace(StructuredDomain):
Martin Reinecke's avatar
Martin Reinecke committed
29
30
    """NIFTy subclass for sets of spherical harmonic coefficients.

Philipp Arras's avatar
Philipp Arras committed
31
32
    Its harmonic partner spaces are :class:`~nifty5.domains.hp_space.HPSpace`
    and :class:`~nifty5.domains.gl_space.GLSpace`.
Martin Reinecke's avatar
Martin Reinecke committed
33
34
35
36

    Parameters
    ----------
    lmax : int
Martin Reinecke's avatar
Martin Reinecke committed
37
38
39
40
41
42
43
44
45
        The maximum :math:`l` value of any spherical harmonic coefficient
        :math:`a_{lm}` that is represented by this object.
        Must be :math:`\ge 0`.

    mmax : int, optional
        The maximum :math:`m` value of any spherical harmonic coefficient
        :math:`a_{lm}` that is represented by this object.
        If not supplied, it is set to `lmax`.
        Must be :math:`\ge 0` and :math:`\le` `lmax`.
csongor's avatar
csongor committed
46
47
    """

Martin Reinecke's avatar
Martin Reinecke committed
48
49
    _needed_for_hash = ["_lmax", "_mmax"]

50
51
52
53
54
55
56
57
58
    def __init__(self, lmax, mmax=None):
        self._lmax = np.int(lmax)
        if self._lmax < 0:
            raise ValueError("lmax must be >=0.")
        if mmax is None:
            mmax = self._lmax
        self._mmax = np.int(mmax)
        if self._mmax < 0 or self._mmax > self._lmax:
            raise ValueError("mmax must be >=0 and <=lmax.")
csongor's avatar
csongor committed
59

60
    def __repr__(self):
61
        return ("LMSpace(lmax=%r, mmax=%r)" % (self.lmax, self.mmax))
62

csongor's avatar
csongor committed
63
64
65
    @property
    def harmonic(self):
        return True
csongor's avatar
csongor committed
66
67

    @property
68
    def shape(self):
Martin Reinecke's avatar
Martin Reinecke committed
69
        return (self.size, )
csongor's avatar
csongor committed
70
71

    @property
Martin Reinecke's avatar
Martin Reinecke committed
72
    def size(self):
73
74
        l = self._lmax
        m = self._mmax
75
        # the LMSpace consists of the full triangle (including -m's!),
Theo Steininger's avatar
Theo Steininger committed
76
        # minus two little triangles if mmax < lmax
77
        return (l+1)**2 - (l-m)*(l-m+1)
csongor's avatar
csongor committed
78

Martin Reinecke's avatar
Martin Reinecke committed
79
    @property
80
    def scalar_dvol(self):
81
        return 1.
csongor's avatar
csongor committed
82

83
    def get_k_length_array(self):
84
85
        lmax = self._lmax
        mmax = self._mmax
Martin Reinecke's avatar
Martin Reinecke committed
86
        ldist = np.empty((self.size,), dtype=np.float64)
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
87
88
89
90
91
        ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
        tmp = np.empty((2*lmax+2), dtype=np.float64)
        tmp[0::2] = np.arange(lmax+1)
        tmp[1::2] = np.arange(lmax+1)
        idx = lmax+1
92
93
94
        for m in range(1, mmax+1):
            ldist[idx:idx+2*(lmax+1-m)] = tmp[2*m:]
            idx += 2*(lmax+1-m)
95
        return Field.from_global_data(self, ldist)
96

97
    def get_unique_k_lengths(self):
Martin Reinecke's avatar
Martin Reinecke committed
98
99
        return np.arange(self.lmax+1, dtype=np.float64)

Martin Reinecke's avatar
Martin Reinecke committed
100
101
    @staticmethod
    def _kernel(x, sigma):
102
103
104
        # cf. "All-sky convolution for polarimetry experiments"
        # by Challinor et al.
        # http://arxiv.org/abs/astro-ph/0008228
105
        from ..sugar import exp
Martin Reinecke's avatar
Martin Reinecke committed
106
        return exp((x+1.) * x * (-0.5*sigma*sigma))
Martin Reinecke's avatar
Martin Reinecke committed
107

108
    def get_fft_smoothing_kernel_function(self, sigma):
Martin Reinecke's avatar
Martin Reinecke committed
109
        return lambda x: self._kernel(x, sigma)
110

csongor's avatar
csongor committed
111
112
    @property
    def lmax(self):
Martin Reinecke's avatar
Martin Reinecke committed
113
114
115
116
        """int : maximum allowed :math:`l`

        The maximum :math:`l` value of any spherical harmonic
        coefficient :math:`a_{lm}` that is represented in this domain.
117
        """
csongor's avatar
csongor committed
118
119
120
121
        return self._lmax

    @property
    def mmax(self):
Martin Reinecke's avatar
Martin Reinecke committed
122
123
124
125
        """int : maximum allowed :math:`m`

        The maximum :math:`m` value of any spherical harmonic
        coefficient :math:`a_{lm}` that is represented in this domain.
126
        """
127
        return self._mmax
Martin Reinecke's avatar
Martin Reinecke committed
128
129

    def get_default_codomain(self):
Philipp Arras's avatar
Philipp Arras committed
130
        """Returns a :class:`~nifty5.domains.gl_space.GLSpace` object, which is
Martin Reinecke's avatar
Martin Reinecke committed
131
132
133
134
135
136
        capable of storing an accurate representation of data residing on
        `self`.

        Returns
        -------
        GLSpace
Martin Reinecke's avatar
Martin Reinecke committed
137
            The partner domain
Martin Reinecke's avatar
Martin Reinecke committed
138
        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
139
        from ..domains.gl_space import GLSpace
Martin Reinecke's avatar
Martin Reinecke committed
140
141
142
        return GLSpace(self.lmax+1, self.mmax*2+1)

    def check_codomain(self, codomain):
Martin Reinecke's avatar
Martin Reinecke committed
143
144
145
146
147
148
149
150
        """Raises `TypeError` if `codomain` is not a matching partner domain
        for `self`.

        Notes
        -----
        This function only checks whether `codomain` is of type
        :class:`GLSpace` or :class:`HPSpace`.
        """
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
151
152
        from ..domains.gl_space import GLSpace
        from ..domains.hp_space import HPSpace
Martin Reinecke's avatar
Martin Reinecke committed
153
        if not isinstance(codomain, (GLSpace, HPSpace)):
154
            raise TypeError("codomain must be a GLSpace or HPSpace.")