lm_space.py 4.9 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

csongor's avatar
csongor committed
19 20
from __future__ import division
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
21
from .structured_domain import StructuredDomain
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
22
from ..field import Field, exp
Martin Reinecke's avatar
fixes  
Martin Reinecke committed
23
from .. import dobj
Theo Steininger's avatar
Theo Steininger committed
24

csongor's avatar
csongor committed
25

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

Martin Reinecke's avatar
Martin Reinecke committed
29 30
    Its harmonic partner spaces are :class:`~nifty4.domains.hp_space.HPSpace`
    and :class:`~nifty4.domains.gl_space.GLSpace`.
Martin Reinecke's avatar
Martin Reinecke committed
31 32 33 34

    Parameters
    ----------
    lmax : int
Martin Reinecke's avatar
Martin Reinecke committed
35 36 37 38 39 40 41 42 43
        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
44 45
    """

46
    def __init__(self, lmax, mmax=None):
Martin Reinecke's avatar
Martin Reinecke committed
47
        super(LMSpace, self).__init__()
48 49 50 51 52 53 54 55 56
        self._needed_for_hash += ["_lmax", "_mmax"]
        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
57

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

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

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

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

77
    def scalar_dvol(self):
78
        return 1.
csongor's avatar
csongor committed
79

80
    def get_k_length_array(self):
81 82
        lmax = self._lmax
        mmax = self._mmax
Martin Reinecke's avatar
Martin Reinecke committed
83
        ldist = np.empty((self.size,), dtype=np.float64)
Martin Reinecke's avatar
stage1  
Martin Reinecke committed
84 85 86 87 88
        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
89 90 91
        for m in range(1, mmax+1):
            ldist[idx:idx+2*(lmax+1-m)] = tmp[2*m:]
            idx += 2*(lmax+1-m)
Martin Reinecke's avatar
fixes  
Martin Reinecke committed
92
        return Field((self,), dobj.from_global_data(ldist))
93

94
    def get_unique_k_lengths(self):
Martin Reinecke's avatar
Martin Reinecke committed
95 96
        return np.arange(self.lmax+1, dtype=np.float64)

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

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):
Martin Reinecke's avatar
Martin Reinecke committed
130 131 132 133 134 135 136
        """Returns a :class:`~nifty4.domains.gl_space.GLSpace` object, which is
        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
Martin Reinecke committed
139 140 141 142
        from .. import GLSpace
        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
Martin Reinecke committed
151 152
        from .. import GLSpace, HPSpace
        if not isinstance(codomain, (GLSpace, HPSpace)):
153
            raise TypeError("codomain must be a GLSpace or HPSpace.")