power_conversion.py 10.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2015 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## 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/>.
21

Marco Selig's avatar
Marco Selig committed
22
import numpy as np
23
from numpy import pi
24
from nifty.config import about
25
from nifty.field import Field
26 27 28 29 30 31 32 33 34 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 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 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
from nifty.nifty_simple_math import sqrt, exp, log


def power_backward_conversion_lm(k_space, p, mean=None):
    """
        This function is designed to convert a theoretical/statistical power
        spectrum of a log-normal field to the theoretical power spectrum of
        the underlying Gaussian field.
        The function only works for power spectra defined for lm_spaces

        Parameters
        ----------
        k_space : nifty.rg_space,
            a regular grid space with the attribute `Fourier = True`
        p : np.array,
            the power spectrum of the log-normal field.
            Needs to have the same number of entries as
            `k_space.get_power_indices()[0]`
        mean : float, *optional*
            specifies the mean of the log-normal field. If `mean` is not
            specified the function will use the monopole of the power spectrum.
            If it is specified the function will NOT use the monopole of the
            spectrum. (default: None)
            WARNING: a mean that is too low can violate positive definiteness
            of the log-normal field. In this case the function produces an
            error.

        Returns
        -------
        mean : float,
            the recovered mean of the underlying Gaussian distribution.
        p1 : np.array,
            the power spectrum of the underlying Gaussian field, where the
            monopole has been set to zero. Eventual monopole power has been
            shifted to the mean.

        References
        ----------
        .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
            `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
    """

    p = np.copy(p)
    if(mean is not None):
        p[0] = 4*pi*mean**2

    klen = k_space.get_power_indices()[0]
    C_0_Omega = Field(k_space,val=0)
    C_0_Omega.val[:len(klen)] = p*sqrt(2*klen+1)/sqrt(4*pi)
    C_0_Omega = C_0_Omega.transform()

    if(np.any(C_0_Omega.val<0.)):
        raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
        return None

    lC = log(C_0_Omega)

    Z = lC.transform()

    spec = Z.val[:len(klen)]

    mean = (spec[0]-0.5*sqrt(4*pi)*log((p*(2*klen+1)/(4*pi)).sum()))/sqrt(4*pi)

    spec[0] = 0.

    spec = spec*sqrt(4*pi)/sqrt(2*klen+1)

    spec = np.real(spec)

    if(np.any(spec<0.)):
        spec = spec*(spec>0.)
        about.warnings.cprint("WARNING: negative modes set to zero.")

    return mean.real,spec


def power_forward_conversion_lm(k_space,p,mean=0):
    """
        This function is designed to convert a theoretical/statistical power
        spectrum of a Gaussian field to the theoretical power spectrum of
        the exponentiated field.
        The function only works for power spectra defined for lm_spaces

        Parameters
        ----------
        k_space : nifty.rg_space,
            a regular grid space with the attribute `Fourier = True`
        p : np.array,
            the power spectrum of the Gaussian field.
            Needs to have the same number of entries as
            `k_space.get_power_indices()[0]`
        m : float, *optional*
            specifies the mean of the Gaussian field (default: 0).

        Returns
        -------
        p1 : np.array,
            the power spectrum of the exponentiated Gaussian field.

        References
        ----------
        .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
            `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
    """
    m = mean
    klen = k_space.get_power_indices()[0]
    C_0_Omega = Field(k_space,val=0)
    C_0_Omega.val[:len(klen)] = p*sqrt(2*klen+1)/sqrt(4*pi)
    C_0_Omega = C_0_Omega.transform()

    C_0_0 = (p*(2*klen+1)/(4*pi)).sum()

    exC = exp(C_0_Omega+C_0_0+2*m)

    Z = exC.transform()

    spec = Z.val[:len(klen)]

    spec = spec*sqrt(4*pi)/sqrt(2*klen+1)

    spec = np.real(spec)

    if(np.any(spec<0.)):
        spec = spec*(spec>0.)
        about.warnings.cprint("WARNING: negative modes set to zero.")

    return spec
Marco Selig's avatar
Marco Selig committed
153 154


155
def power_backward_conversion_rg(k_space, p, mean=None, bare=True):
Marco Selig's avatar
Marco Selig committed
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
    """
        This function is designed to convert a theoretical/statistical power
        spectrum of a log-normal field to the theoretical power spectrum of
        the underlying Gaussian field.
        The function only works for power spectra defined for rg_spaces

        Parameters
        ----------
        k_space : nifty.rg_space,
            a regular grid space with the attribute `Fourier = True`
        p : np.array,
            the power spectrum of the log-normal field.
            Needs to have the same number of entries as
            `k_space.get_power_indices()[0]`
        mean : float, *optional*
            specifies the mean of the log-normal field. If `mean` is not
            specified the function will use the monopole of the power spectrum.
            If it is specified the function will NOT use the monopole of the
            spectrum (default: None).
            WARNING: a mean that is too low can violate positive definiteness
            of the log-normal field. In this case the function produces an
            error.
        bare : bool, *optional*
            whether `p` is the bare power spectrum or not (default: True).

        Returns
        -------
        mean : float,
            the recovered mean of the underlying Gaussian distribution.
        p1 : np.array,
            the power spectrum of the underlying Gaussian field, where the
            monopole has been set to zero. Eventual monopole power has been
            shifted to the mean.

        References
        ----------
192 193
        .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter
               power spectrum";
Marco Selig's avatar
Marco Selig committed
194 195
            `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
    """
196

197
    pindex = k_space.power_indices['pindex']
198 199
    weight = k_space.get_weight()

200
    monopole_index = pindex.argmin()
201 202

    # Cast the supplied spectrum
203
    spec = k_space.enforce_power(p)
204 205 206 207
    # Now we mimick the weightning behaviour of
    # spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
    # by appliying the weight from the k_space
    if bare:
208
        spec *= weight
209

210 211 212 213 214
    #TODO: Does this realy set the mean to the monopole as promised in the docs? -> Check!
    if mean is None:
        mean = 0.
    else:
        spec[0] = 0.
215 216 217

    p_val = pindex.apply_scalar_function(lambda x: spec[x],
                                         dtype=spec.dtype.type)
218
    power_field = Field(k_space, val=p_val, zerocenter=True).transform()
219 220
    power_field += (mean**2)

221 222
    if power_field.min() < 0:
        raise ValueError(about._errors.cstring(
223 224
            "ERROR: spectrum or mean incompatible with positive " +
            "definiteness. \n Try increasing the mean."))
225

226
    log_of_power_field = power_field.apply_scalar_function(np.log,
227 228 229
                                                           inplace=True)
    power_spectrum_1 = log_of_power_field.power()**(0.5)
    power_spectrum_1[0] = log_of_power_field.transform()[monopole_index]
230

231 232 233
    power_spectrum_0 = k_space.calc_weight(p_val).sum() + (mean**2)
    power_spectrum_0 = np.log(power_spectrum_0)
    power_spectrum_0 *= (0.5 / weight)
Marco Selig's avatar
Marco Selig committed
234

235
    log_mean = weight * (power_spectrum_1[0] - power_spectrum_0)
Marco Selig's avatar
Marco Selig committed
236

237
    power_spectrum_1[0] = 0.
Marco Selig's avatar
Marco Selig committed
238

239 240 241 242 243
    # Mimik
    # power_operator(k_space,spec=power_spectrum_1,bare=False).\
    #  get_power(bare=True).real
    if bare:
        power_spectrum_1 /= weight
Marco Selig's avatar
Marco Selig committed
244

245
    return log_mean.real, power_spectrum_1.real
Marco Selig's avatar
Marco Selig committed
246 247


248
def power_forward_conversion_rg(k_space, p, mean=0, bare=True):
Marco Selig's avatar
Marco Selig committed
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274
    """
        This function is designed to convert a theoretical/statistical power
        spectrum of a Gaussian field to the theoretical power spectrum of
        the exponentiated field.
        The function only works for power spectra defined for rg_spaces

        Parameters
        ----------
        k_space : nifty.rg_space,
            a regular grid space with the attribute `Fourier = True`
        p : np.array,
            the power spectrum of the Gaussian field.
            Needs to have the same number of entries as
            `k_space.get_power_indices()[0]`
        mean : float, *optional*
            specifies the mean of the Gaussian field (default: 0).
        bare : bool, *optional*
            whether `p` is the bare power spectrum or not (default: True).

        Returns
        -------
        p1 : np.array,
            the power spectrum of the exponentiated Gaussian field.

        References
        ----------
275 276
        .. [#] M. Greiner and T.A. Ensslin,
            "Log-transforming the matter power spectrum";
Marco Selig's avatar
Marco Selig committed
277 278 279
            `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
    """

280
    pindex = k_space.power_indices['pindex']
281 282
    weight = k_space.get_weight()
    # Cast the supplied spectrum
283
    spec = k_space.enforce_power(p)
284 285 286 287
    # Now we mimick the weightning behaviour of
    # spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
    # by appliying the weight from the k_space
    if bare:
288
        spec *= weight
289 290 291 292 293

    S_val = pindex.apply_scalar_function(lambda x: spec[x],
                                         dtype=spec.dtype.type)

    # S_x is a field
294
    S_x = Field(k_space, val=S_val, zerocenter=True).transform()
295 296 297 298
    # s_0 is a scalar
    s_0 = k_space.calc_weight(S_val, power=1).sum()

    # Calculate the new power_field
299 300
    S_x += s_0
    S_x += 2*mean
301

302 303
    power_field = S_x.apply_scalar_function(np.exp, inplace=True)

304
    new_spec = power_field.power()**(0.5)
Marco Selig's avatar
Marco Selig committed
305

306 307 308 309
    # Mimik
    # power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real
    if bare:
        new_spec /= weight
Marco Selig's avatar
Marco Selig committed
310

311
    return new_spec.real