Commit 3e343cd4 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'modular_amplitudes' into 'NIFTy_5'

Rewrite power spectrum priors

See merge request ift/nifty-dev!120
parents 309dd635 67b317d9
......@@ -29,10 +29,11 @@ def get_random_LOS(n_los):
if __name__ == '__main__':
# FIXME description of the tutorial
np.random.seed(42)
np.seterr(all='raise')
position_space = ift.RGSpace([128, 128])
# Setting up an amplitude model
A = ift.AmplitudeModel(position_space, 16, 1, 10, -4., 1, 0., 1.)
A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -4., 1, 1., 1.)
# Building the model for a correlated signal
harmonic_space = position_space.get_default_codomain()
......@@ -40,10 +41,11 @@ if __name__ == '__main__':
power_space = A.target[0]
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
correlated_field = ht(
power_distributor(A)*ift.FieldAdapter(harmonic_space, "xi"))
vol = harmonic_space.scalar_dvol
vol = ift.ScalingOperator(vol ** (-0.5),harmonic_space)
correlated_field = ht(vol(power_distributor(A))*ift.FieldAdapter(harmonic_space, "xi"))
# alternatively to the block above one can do:
# correlated_field = ift.CorrelatedField(position_space, A)
#correlated_field = ift.CorrelatedField(position_space, A)
# apply some nonlinearity
signal = ift.positive_tanh(correlated_field)
......
......@@ -34,6 +34,7 @@ from .operators.inversion_enabler import InversionEnabler
from .operators.laplace_operator import LaplaceOperator
from .operators.linear_operator import LinearOperator
from .operators.mask_operator import MaskOperator
from .operators.offset_operator import OffsetOperator
from .operators.qht_operator import QHTOperator
from .operators.regridding_operator import RegriddingOperator
from .operators.sampling_enabler import SamplingEnabler
......@@ -73,7 +74,7 @@ from .minimization.kl_energy import KL_Energy
from .sugar import *
from .plot import Plot
from .library.amplitude_model import AmplitudeModel, SmoothAmplitudeModel
from .library.amplitude_model import AmplitudeModel
from .library.inverse_gamma_model import InverseGammaModel
from .library.los_response import LOSResponse
......
......@@ -91,22 +91,22 @@ class LogRGSpace(StructuredDomain):
return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True)
def get_k_length_array(self):
ib = dobj.ibegin_from_shape(self._shape)
res = np.arange(self.local_shape[0], dtype=np.float64) + ib[0]
res = np.minimum(res, self.shape[0] - res)*self.bindistances[0]
if len(self.shape) == 1:
return Field.from_local_data(self, res)
res *= res
for i in range(1, len(self.shape)):
tmp = np.arange(self.local_shape[i], dtype=np.float64) + ib[i]
tmp = np.minimum(tmp, self.shape[i] - tmp)*self.bindistances[i]
tmp *= tmp
res = np.add.outer(res, tmp)
return Field.from_local_data(self, np.sqrt(res))
def get_expk_length_array(self):
# FIXME This is a hack! Only for plotting. Seems not to be the final version.
out = exp(self.get_k_length_array()).to_global_data_rw()
out[1:] = out[:-1]
out[0] = 0
return Field.from_global_data(self, out)
if not self.harmonic:
raise NotImplementedError
ks = self.get_k_array()
return Field.from_global_data(self, np.linalg.norm(ks, axis=0))
def get_k_array(self):
ndim = len(self.shape)
k_array = np.zeros((ndim,) + self.shape)
dist = self.bindistances
for i in range(ndim):
ks = np.zeros(self.shape[i])
ks[1:] = np.minimum(self.shape[i] - 1 - np.arange(self.shape[i]-1), np.arange(self.shape[i]-1)) * dist[i]
if self.harmonic:
ks[0] = np.nan
else:
ks[0] = -np.inf
ks[1:] += self.t_0[i]
k_array[i] += ks.reshape((1,)*i + (self.shape[i],) + (1,)*(ndim-i-1))
return k_array
......@@ -23,13 +23,11 @@ import numpy as np
from ..compat import *
from ..domains.power_space import PowerSpace
from ..field import Field
from ..multi_domain import MultiDomain
from ..operators.operator import Operator
from ..sugar import makeOp, sqrt
def _ceps_kernel(dof_space, k, a, k0):
return a**2/(1+(k/(k0*dof_space.bindistances[0]))**2)**2
return a**2/(1 + (k/k0)**2)**2
def create_cepstrum_amplitude_field(domain, cepstrum):
......@@ -45,22 +43,12 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
"""
dim = len(domain.shape)
dist = domain.bindistances
shape = domain.shape
# Prepare q_array
q_array = np.zeros((dim,) + shape)
if dim == 1:
ks = domain.get_k_length_array().to_global_data()
q_array = np.array([ks])
else:
for i in range(dim):
ks = np.minimum(shape[i] - np.arange(shape[i]) +
1, np.arange(shape[i])) * dist[i]
q_array[i] += ks.reshape((1,)*i + (shape[i],) + (1,)*(dim-i-1))
q_array = domain.get_k_array()
# Fill cepstrum field (all non-zero modes)
no_zero_modes = (slice(1, None),) * dim
no_zero_modes = (slice(1, None),)*dim
ks = q_array[(slice(None),) + no_zero_modes]
cepstrum_field = np.zeros(shape)
cepstrum_field[no_zero_modes] = cepstrum(ks)
......@@ -78,81 +66,55 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
return Field.from_global_data(domain, cepstrum_field)
class AmplitudeModel(Operator):
def CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
'''
Computes a smooth power spectrum.
Output lives in PowerSpace.
Parameters
----------
Npixdof : #pix in dof_space
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0
'''
from ..operators.qht_operator import QHTOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
qht = QHTOperator(target=logk_space)
dof_space = qht.domain[0]
sym = SymmetrizingOperator(logk_space)
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
res = sym(qht(makeOp(sqrt(cepstrum))))
if not zero_mode:
shp = res.target.shape
mask = np.ones(shp)
mask[(0,)*len(shp)] = 0.
mask = makeOp(Field.from_global_data(res.target, mask))
res = mask(res)
return res
def SlopeModel(logk_space, sm, sv, im, iv):
'''
Parameters
----------
sm, sv : slope_mean = expected exponent of power law (e.g. -4),
slope_variance (default=1)
im, iv : y-intercept_mean, y-intercept_variance of power_slope
im, iv : y-intercept_mean, y-intercept_std of power_slope
'''
def __init__(self, s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi']):
from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
h_space = s_space.get_default_codomain()
self._exp_transform = ExpTransform(PowerSpace(h_space), Npixdof)
logk_space = self._exp_transform.domain[0]
qht = QHTOperator(target=logk_space)
dof_space = qht.domain[0]
sym = SymmetrizingOperator(logk_space)
phi_mean = np.array([sm, im])
phi_sig = np.array([sv, iv])
self._slope = SlopeOperator(logk_space, phi_sig)
self._norm_phi_mean = Field.from_global_data(self._slope.domain,
phi_mean/phi_sig)
self._domain = MultiDomain.make({keys[0]: dof_space,
keys[1]: self._slope.domain})
self._target = self._exp_transform.target
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
self._smooth_op = sym(qht(makeOp(sqrt(cepstrum))))
self._keys = tuple(keys)
self._qht = qht
self._ceps = makeOp(sqrt(cepstrum))
def apply(self, x):
self._check_input(x)
smooth_spec = self._smooth_op(x[self._keys[0]])
phi = x[self._keys[1]] + self._norm_phi_mean
linear_spec = self._slope(phi)
loglog_spec = smooth_spec + linear_spec
return self._exp_transform((0.5*loglog_spec).exp())
@property
def qht(self):
return self._qht
@property
def ceps(self):
return self._ceps
@property
def norm_phi_mean(self):
return self._norm_phi_mean
class SmoothAmplitudeModel(Operator):
from ..operators.slope_operator import SlopeOperator
from ..operators.offset_operator import OffsetOperator
phi_mean = np.array([sm, im + sm*logk_space.t_0[0]])
phi_sig = np.array([sv, iv])
slope = SlopeOperator(logk_space)
phi_mean = Field.from_global_data(slope.domain, phi_mean)
phi_sig = Field.from_global_data(slope.domain, phi_sig)
return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
def AmplitudeModel(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi'], zero_mode=True):
'''
Computes a smooth power spectrum.
Output lives in PowerSpace.
......@@ -165,39 +127,26 @@ class SmoothAmplitudeModel(Operator):
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0
'''
def __init__(self, s_space, Npixdof, ceps_a, ceps_k, name='tau'):
from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
h_space = s_space.get_default_codomain()
self._exp_transform = ExpTransform(PowerSpace(h_space), Npixdof)
logk_space = self._exp_transform.domain[0]
qht = QHTOperator(target=logk_space)
dof_space = qht.domain[0]
self._domain = MultiDomain.make({name: dof_space})
self._target = self._exp_transform.target
sm, sv : slope_mean = expected exponent of power law (e.g. -4),
slope_variance (default=1)
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
im, iv : y-intercept_mean, y-intercept_variance of power_slope
'''
self._smooth_op = qht(makeOp(sqrt(cepstrum)))
self._key = name
from ..operators.exp_transform import ExpTransform
from ..operators.simple_linear_operators import FieldAdapter
from ..operators.scaling_operator import ScalingOperator
self._qht = qht
self._ceps = makeOp(sqrt(cepstrum))
h_space = s_space.get_default_codomain()
et = ExpTransform(PowerSpace(h_space), Npixdof)
logk_space = et.domain[0]
def apply(self, x):
self._check_input(x)
smooth_spec = self._smooth_op(x[self._key])
loglog_spec = smooth_spec
return self._exp_transform((0.5*loglog_spec).exp())
smooth = CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
linear = SlopeModel(logk_space, sm, sv, im, iv)
@property
def qht(self):
return self._qht
fa_smooth = FieldAdapter(smooth.domain, keys[0])
fa_linear = FieldAdapter(linear.domain, keys[1])
@property
def ceps(self):
return self._ceps
fac = ScalingOperator(0.5, smooth.target)
return et((fac(smooth(fa_smooth) + linear(fa_linear))).exp())
......@@ -25,6 +25,7 @@ from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import FieldAdapter
from ..operators.scaling_operator import ScalingOperator
def CorrelatedField(s_space, amplitude_model, name='xi'):
......@@ -45,7 +46,9 @@ def CorrelatedField(s_space, amplitude_model, name='xi'):
p_space = amplitude_model.target[0]
power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_model)
return ht(A*FieldAdapter(h_space, name))
vol = h_space.scalar_dvol
vol = ScalingOperator(vol**(-0.5), h_space)
return ht(vol(A)*FieldAdapter(h_space, name))
def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
......
# 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/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from __future__ import absolute_import, division, print_function
from ..compat import *
from .operator import Operator
class OffsetOperator(Operator):
def __init__(self, field):
self._field = field
self._domain = self._target = field.domain
def apply(self, x):
self._check_input(x)
return x + self._field
......@@ -108,6 +108,7 @@ class _FunctionApplier(Operator):
return getattr(x, self._funcname)()
# FIXME Is this class used except in _OpChain? Can it be merged?
class _CombinedOperator(Operator):
def __init__(self, ops, _callingfrommake=False):
if not _callingfrommake:
......
......@@ -24,7 +24,6 @@ from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..multi_domain import MultiDomain
from ..multi_field import MultiField
from ..sugar import full
from .endomorphic_operator import EndomorphicOperator
from .linear_operator import LinearOperator
......
......@@ -46,43 +46,27 @@ class SlopeOperator(LinearOperator):
sigmas : np.array, shape=(2,)
The slope variance and the y-intercept variance.
"""
def __init__(self, target, sigmas):
def __init__(self, target):
if not isinstance(target, LogRGSpace):
raise TypeError
if len(target.shape) != 1:
raise ValueError("Slope Operator only works for ndim == 1")
self._domain = DomainTuple.make(UnstructuredDomain((2,)))
self._target = DomainTuple.make(target)
self._capability = self.TIMES | self.ADJOINT_TIMES
self._sigmas = sigmas
self.ndim = len(self.target[0].shape)
self.pos = np.zeros((self.ndim,) + self.target[0].shape)
if self.ndim == 1:
self.pos[0] = self.target[0].get_k_length_array().to_global_data()
else:
shape = self.target[0].shape
for i in range(self.ndim):
rng = np.arange(target.shape[i])
tmp = np.minimum(
rng, target.shape[i]+1-rng) * target.bindistances[i]
self.pos[i] += tmp.reshape(
(1,)*i + (shape[i],) + (1,)*(self.ndim-i-1))
pos = self.target[0].get_k_array() - self.target[0].t_0[0]
self._pos = pos[0, 1:]
def apply(self, x, mode):
self._check_input(x, mode)
# Times
inp = x.to_global_data()
if mode == self.TIMES:
inp = x.to_global_data()
res = self._sigmas[-1] * inp[-1]
for i in range(self.ndim):
res += self._sigmas[i] * inp[i] * self.pos[i]
return Field.from_global_data(self.target, res)
# Adjoint times
res = np.zeros(self.domain[0].shape, dtype=x.dtype)
xglob = x.to_global_data()
res[-1] = np.sum(xglob) * self._sigmas[-1]
for i in range(self.ndim):
res[i] = np.sum(self.pos[i] * xglob) * self._sigmas[i]
return Field.from_global_data(self.domain, res)
res = np.empty(self.target.shape, dtype=x.dtype)
res[0] = 0
res[1:] = inp[1] + inp[0]*self._pos
else:
res = np.array(
[np.sum(self._pos*inp[1:]),
np.sum(inp[1:])], dtype=x.dtype)
return Field.from_global_data(self._tgt(mode), res)
......@@ -73,8 +73,7 @@ class Consistency_Tests(unittest.TestCase):
def testSlopeOperator(self, args, dtype):
tmp = ift.ExpTransform(ift.PowerSpace(args[0]), args[1], args[2])
tgt = tmp.domain[0]
sig = np.array([0.3, 0.13])
op = ift.SlopeOperator(tgt, sig)
op = ift.SlopeOperator(tgt)
ift.extra.consistency_check(op, dtype, dtype)
@expand(product(_h_spaces + _p_spaces + _pow_spaces,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment