Scheduled maintenance on Monday 2019-06-24 between 10:00-11:00 CEST

...
 
Commits (30)
......@@ -75,7 +75,7 @@ from .sugar import *
from .plot import Plot
from .library.smooth_linear_amplitude import (
SLAmplitude, LinearSLAmplitude, CepstrumOperator)
SLAmplitude, CepstrumOperator)
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator,
......@@ -83,7 +83,8 @@ from .library.dynamic_operator import (dynamic_operator,
from .library.light_cone_operator import LightConeOperator
from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.correlated_fields import CorrelatedField, MfCorrelatedField
from .library.correlated_fields import (CorrelatedField, MfCorrelatedField,
MfPartiallyCorrelatedField)
from .library.adjust_variances import (make_adjust_variances_hamiltonian,
do_adjust_variances)
from .library.gridder import GridderMaker
......
......@@ -15,6 +15,7 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from functools import reduce
from operator import mul
......@@ -24,8 +25,14 @@ from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape
from ..library.inverse_gamma_operator import InverseGammaOperator
from ..operators.value_inserter import ValueInserter
from ..sugar import makeOp
from ..field import Field
def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
def CorrelatedField(target, amplitude_operator, keys=['xi', "zero_mode"],
codomain=None, za=None, zq=None):
"""Constructs an operator which turns a white Gaussian excitation field
into a correlated field.
......@@ -69,14 +76,36 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol**-0.5
separate_zero_mode_prior = za is not None and zq is not None
if separate_zero_mode_prior:
za, zq = float(za), float(zq)
else:
if za is not None or zq is not None:
raise ValueError("za and zq have to be given together")
# When doubling the resolution of `tgt` the value of the highest k-mode
# will scale with a square root. `vol` cancels this effect such that the
# same power spectrum can be used for the spaces with the same volume,
# different resolutions and the same object in them.
return ht(vol*A*ducktape(h_space, None, name))
if not separate_zero_mode_prior:
return ht(vol*A*ducktape(h_space, None, keys[0]))
else:
zero_mode = ValueInserter(A.target, (0,)*len(A.target.shape))
zero_mode = (
zero_mode @
InverseGammaOperator(zero_mode.domain, za, zq).ducktape(keys[1]))
mask = np.ones(A.target.shape)
mask[(0,)*len(A.target.shape)] = 0.
mask = makeOp(Field.from_global_data(A.target, mask))
def MfCorrelatedField(target, amplitudes, name='xi'):
return ht(mask @ (vol*A*ducktape(h_space, None, keys[0])) + zero_mode)
# def MfCorrelatedField(target, amplitudes, name='xi'):
def MfCorrelatedField(target, amplitudes, keys=['xi', "zero_mode"],
za=None, zq=None):
"""Constructs an operator which turns white Gaussian excitation fields
into a correlated field defined on a DomainTuple with two entries and two
separate correlation structures.
......@@ -118,17 +147,93 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
ht = ht2 @ ht1
psp = [aa.target[0] for aa in amplitudes]
pd0 = PowerDistributor(hsp, psp[0], 0)
pd1 = PowerDistributor(pd0.domain, psp[1], 1)
pd0 = PowerDistributor(hsp, power_space=psp[0], space=0)
pd1 = PowerDistributor(pd0.domain, power_space=psp[1], space=1)
pd = pd0 @ pd1
dd0 = ContractionOperator(pd.domain, 1).adjoint
dd1 = ContractionOperator(pd.domain, 0).adjoint
d = [dd0, dd1]
blowup0 = ContractionOperator(pd.domain, 1).adjoint
blowup1 = ContractionOperator(pd.domain, 0).adjoint
d = [blowup0, blowup1]
a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
a = [BUp @ amplitudes[ii] for ii, BUp in enumerate(d)]
a = reduce(mul, a)
A = pd @ a
# For `vol` see comment in `CorrelatedField`
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ducktape(hsp, None, name))
separate_zero_mode_prior = za is not None and zq is not None
if separate_zero_mode_prior:
za, zq = float(za), float(zq)
else:
if za is not None or zq is not None:
raise ValueError("za and zq have to be given together")
# When doubling the resolution of `tgt` the value of the highest k-mode
# will scale with a square root. `vol` cancels this effect such that the
# same power spectrum can be used for the spaces with the same volume,
# different resolutions and the same object in them.
if not separate_zero_mode_prior:
return ht(vol*A*ducktape(hsp, None, keys[0]))
else:
zero_mode = ValueInserter(A.target, (0,)*len(A.target.shape))
zero_mode = (
zero_mode @
InverseGammaOperator(zero_mode.domain, za, zq).ducktape(keys[1]))
mask = np.ones(A.target.shape)
mask[(0,)*len(A.target.shape)] = 0.
mask = makeOp(Field.from_global_data(A.target, mask))
return ht(mask @ (vol*A*ducktape(hsp, None, keys[0])) + zero_mode)
def MfPartiallyCorrelatedField(target, energy_amplitude, name='xi_u'):
"""Constructs an operator which turns white Gaussian excitation fields
into a correlated field defined on a DomainTuple with two entries.
On the first domain, a white correlation structure is assumed. On
the second domain, the correlation structures given by energy_amplitude
are used.
This operator may be used as a model for multi-frequency reconstructions
with correlation structure only in the energy direction.
Parameters
----------
target : DomainTuple or tuple of Domain
Target of the operator. Must contain exactly two spaces.
It is assumed that the second space is the energy domain, so it must
be a one-dimensional RGSpace.
energy_amplitude: Operator
amplitude operator for the energy correlation structure
name : string
:class:`MultiField` key for xi-field.
Returns
-------
Operator
Correlated field
Notes
-----
In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore
the operator constructed by this method will output a correlated field
with *periodic* boundary conditions. If a non-periodic field is needed,
one needs to combine this operator with a :class:`FieldZeroPadder` or even
two (one for the energy and one for the spatial subdomain)
"""
tgt = DomainTuple.make(target)
if len(tgt) != 2:
raise ValueError
h_space = DomainTuple.make([dom.get_default_codomain() for dom in tgt])
ht1 = HarmonicTransformOperator(h_space, target=tgt[0], space=0)
ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
ht = ht2 @ ht1
p_space = energy_amplitude.target[0]
power_distributor = PowerDistributor(h_space[-1], p_space)
A = power_distributor(energy_amplitude)
dd = ContractionOperator(h_space, 0).adjoint
return ht((dd @ A)*ducktape(h_space, None, name))
......@@ -26,6 +26,8 @@ from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..sugar import makeOp
from ..library.inverse_gamma_operator import InverseGammaOperator
from ..operators.value_inserter import ValueInserter
def _ceps_kernel(k, a, k0):
......@@ -113,74 +115,6 @@ def CepstrumOperator(target, a, k0):
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
'''Operator for parametrizing smooth amplitudes (square roots of power
spectra).
The general guideline for setting up generative models in IFT is to
transform the problem into the eigenbase of the prior and formulate the
generative model in this base. This is done here for the case of an
amplitude which is smooth and has a linear component (both on
double-logarithmic scale).
This function assembles an :class:`Operator` which maps two a-priori white
Gaussian random fields to a smooth amplitude which is composed out of
a linear and a smooth component.
On double-logarithmic scale, i.e. both x and y-axis on logarithmic scale,
the output of the generated operator is:
AmplitudeOperator = 0.5*(smooth_component + linear_component)
This is then exponentiated and exponentially binned (in this order).
The prior on the linear component is parametrized by four real numbers,
being expected value and prior variance on the slope and the y-intercept
of the linear function.
The prior on the smooth component is parametrized by two real numbers: the
strength and the cutoff of the smoothness prior
(see :class:`CepstrumOperator`).
Parameters
----------
n_pix : int
Number of pixels of the space in which the .
target : PowerSpace
Target of the Operator.
a : float
Strength of smoothness prior (see :class:`CepstrumOperator`).
k0 : float
Cutoff of smothness prior in quefrency space (see
:class:`CepstrumOperator`).
sm : float
Expected exponent of power law.
sv : float
Prior standard deviation of exponent of power law.
im : float
Expected y-intercept of power law. This is the value at t_0 of the
LogRGSpace (see :class:`ExpTransform`).
iv : float
Prior standard deviation of y-intercept of power law.
Returns
-------
Operator
Operator which is defined on the space of white excitations fields and
which returns on its target a power spectrum which consists out of a
smooth and a linear part.
'''
return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0, sm=sm,
sv=sv, im=im, iv=iv, keys=keys).exp()
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
keys=['tau', 'phi']):
'''LinearOperator for parametrizing smooth log-amplitudes (square roots of
power spectra).
Logarithm of SLAmplitude
See documentation of SLAmplitude for more details
'''
if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)):
raise TypeError
......@@ -198,14 +132,16 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
# Linear component
sl = SlopeOperator(dom)
mean = np.array([sm, im + sm*dom.t_0[0]])
# mean = np.array([sm, im + sm*dom.t_0[0]])
mean = np.array([sm, im])
sig = np.array([sv, iv])
mean = Field.from_global_data(sl.domain, mean)
sig = Field.from_global_data(sl.domain, sig)
# print(mean.val, sig.val)
linear = sl @ Adder(mean) @ makeOp(sig).ducktape(keys[1])
# Combine linear and smooth component
loglog_ampl = 0.5*(smooth + linear)
# Go from loglog-space to linear-linear-space
return et @ loglog_ampl
return et @ loglog_ampl.exp()