...
  View open merge request
Commits (47)
# 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-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
###############################################################################
# Log-normal field reconstruction from Poissonian data with inhomogenous
# exposure (in case for 2D mode)
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
###############################################################################
import sys
import numpy as np
import nifty5 as ift
def InvGammaFromMeanVar(domain, mean, variance):
base = mean**2 / variance + 1.
return ift.InverseGammaOperator(domain, alpha=base + 1., q=mean * base)
def exposure_2d():
# Structured exposure for 2D mode
x_shape, y_shape = sky_target.shape
exposure = np.ones(sky_target.shape)
exposure[x_shape//3:x_shape//2, :] *= 2.
exposure[x_shape*4//5:x_shape, :] *= .1
exposure[x_shape//2:x_shape*3//2, :] *= 3.
exposure[:, x_shape//3:x_shape//2] *= 2.
exposure[:, x_shape*4//5:x_shape] *= .1
exposure[:, x_shape//2:x_shape*3//2] *= 3.
return ift.Field.from_global_data(sky_target, exposure)
def cd(input):
return input.cast_domain(plot_domain)
if __name__ == '__main__':
np.random.seed(42)
if len(sys.argv) == 2:
mode = int(sys.argv[1])
else:
mode = 1
if not mode in (0, 1):
raise ValueError("Only modes 0 and 1 defined")
filename = f"demo-MultiCorrelatedField-mode_{mode}" + "-{}.png"
# Two-dimensional signal on RGs with inhomogeneous exposure
N = 128
position_space = ift.RGSpace(N)
energy_space = ift.RGSpace(N)
sky_target = ift.DomainTuple.make((position_space, energy_space))
plot_domain = ift.RGSpace((N, N))
exposure = 10. * exposure_2d()
# Define harmonic space and harmonic transform
harmonic_space_p = position_space.get_default_codomain()
harmonic_space_e = energy_space.get_default_codomain()
power_space_p = ift.PowerSpace(harmonic_space_p)
power_space_e = ift.PowerSpace(harmonic_space_e)
amp_p_dict = {
'target': power_space_p,
'n_pix': 64, # 64 spectral bins
'keys': ['tau_p', 'phi_p', 'zm_p'],
# Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
'k0': 0.4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum
'sm': -4., # Expected exponent of power law
'sv': 1., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 2., # Prior standard deviation of -||-
# zero-mode prior
#'zm': 2.5e-4, # Prior mean
#'zv': 1.0e-4 # Prior variance
'zm': None, # Prior mean
'zv': None # Prior variance
}
amp_e_dict = {
'target': power_space_e,
'n_pix': 64, # 64 spectral bins
'keys': ['tau_e', 'phi_e', 'zm_e'],
# Spectral smoothness (affects Gaussian process part)
'a': 3, # relatively high variance of spectral curvature
'k0': 0.4, # quefrency mode below which cepstrum flattens
# Power-law part of spectrum
'sm': -3., # Expected exponent of power law
'sv': 1., # Prior standard deviation of -||-
'im': 0., # Expected y-intercept of power law
'iv': 1., # Prior standard deviation of -||-
# zero-mode prior
#'zm': 3.0e-3, # Prior mean
#'zv': 2.0e-3 # Prior variance
'zm': None, # Prior mean
'zv': None # Prior variance
}
amp_p = ift.SLAmplitude(**amp_p_dict)
amp_e = ift.SLAmplitude(**amp_e_dict)
# Define sky operator
# za=4, zq=5 leads to an amplitude median of one
if mode == 1:
amplitudes = [amp_p, amp_e]
elif mode == 0:
amplitudes = [None, amp_e]
# Diffuse global scaling prior parameters
diff_gs_params = {'gs_alpha': 0.5, 'gs_q': 0.001}
log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes,
**diff_gs_params)
sky = log_diffuse.exp()
M = ift.DiagonalOperator(exposure)
GR = ift.GeometryRemover(M.target)
# Define instrumental response
R = GR(M)
# Generate mock data and define likelihood operator
lamb = R(sky)
mock_position = ift.from_random('normal', sky.domain)
data = lamb(mock_position)
signal = sky(mock_position)
# DEBUG
mock_log_diffuse_vals = log_diffuse(mock_position).to_global_data()
print("std(mock_log_diffuse):", mock_log_diffuse_vals.std())
data = np.random.poisson(data.to_global_data().astype(np.float64))
data = ift.Field.from_global_data(R.target, data)
likelihood = ift.PoissonianEnergy(data)(lamb)
# Settings for minimization
ic_sampling = ift.AbsDeltaEnergyController(name='Sampling',
iteration_limit=200,
deltaE=0.5)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
ic_newton = ift.AbsDeltaEnergyController(name='Newton',
iteration_limit=30,
deltaE=0.5)
minimizer = ift.NewtonCG(ic_newton)
pos = ift.from_random('normal', sky.domain) * 0.01
iteration = 0
# Minimize
for i in range(10):
KL = ift.MetricGaussianKL(pos, H, 10)
KL, convergence = minimizer(KL)
pos = KL.position
iteration += 1
# Plotting
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(sky(pos + sample))
reconst = sc.mean
plot = ift.Plot()
plot.add(cd(signal), title='Signal')
plot.add(cd(GR.adjoint(data)), title='Data')
plot.add(cd(reconst), title='Posterior Mean')
plot.add(cd(reconst - signal), title='Error')
plot.output(xsize=12,
ysize=10,
name=filename.format(f'{iteration:02d}-reconstr'))
p_p = ift.PowerDistributor(harmonic_space_p, power_space_p) @ amp_p
p_e = ift.PowerDistributor(harmonic_space_e, power_space_e) @ amp_e
if mode == 1:
amplitude_specs = [(p_p, 'position'), (p_e, 'energy')]
elif mode == 0:
amplitude_specs = [(p_e, 'energy')]
for p, n in amplitude_specs:
sample_specs = [
ift.power_analyze(p.force(pos + s)) for s in KL.samples
]
sc = ift.StatCalculator()
for spec in sample_specs:
sc.add(spec)
plot = ift.Plot()
plot.add(
sample_specs +
[sc.mean, ift.power_analyze(p.force(mock_position))],
linewidth=[1.] * len(sample_specs) + [3., 3.],
label=[None] * len(sample_specs) +
['Posterior Mean', 'Ground Truth'],
title=f"Sampled Posterior {n} spectrum")
plot.output(xsize=12,
ysize=10,
name=filename.format(f'{iteration:02d}-power_{n}'))
print("Saved results as '{}'.".format(filename.format('*')))
# Print InverseGamma prior fits
print("Variance Prior fit:")
key = 'xi_pspec_scaling'
op = ift.InverseGammaOperator(sky.domain[key], diff_gs_params['gs_alpha'],
diff_gs_params['gs_q'])
print(f"{key}:")
print(f"\tmock: {float(op(mock_position[key]).local_data):+1.2e}")
print(f"\t rec: {float(op(pos[key]).local_data):+1.2e}\n")
......@@ -73,7 +73,11 @@ if __name__ == '__main__':
'sm': -5, # preferred power-law slope
'sv': .5, # low variance of power-law slope
'im': 0, # y-intercept mean, in-/decrease for more/less contrast
'iv': .3 # y-intercept variance
'iv': .3, # y-intercept variance
# Zero mode prior (InverseGammaOperator)
'zm': 1.0, # Prior mean
'zv': 0.01 # Prior variance
}
A = ift.SLAmplitude(**dct)
......
......@@ -44,7 +44,7 @@ from .operators.symmetrizing_operator import SymmetrizingOperator
from .operators.block_diagonal_operator import BlockDiagonalOperator
from .operators.outer_product_operator import OuterProduct
from .operators.simple_linear_operators import (
VdotOperator, ConjugationOperator, Realizer,
VdotOperator, ConjugationOperator, WeightApplier, Realizer,
FieldAdapter, ducktape, GeometryRemover, NullOperator,
MatrixProductOperator, PartialExtractor)
from .operators.value_inserter import ValueInserter
......@@ -53,6 +53,7 @@ from .operators.energy_operators import (
BernoulliEnergy, StandardHamiltonian, AveragedEnergy, QuadraticFormOperator,
Squared2NormOperator)
from .operators.convolution_operators import FuncConvolutionOperator
from .operators.operator_normalizer import OperatorNormalizer
from .probing import probe_with_posterior_samples, probe_diagonal, \
StatCalculator, approximation2endo
......@@ -85,7 +86,7 @@ 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, MultiCorrelatedField
from .library.adjust_variances import (make_adjust_variances_hamiltonian,
do_adjust_variances)
from .library.gridder import GridderMaker
......
......@@ -17,12 +17,18 @@
from functools import reduce
from operator import mul
import numpy as np
from ..field import Field
from ..domain_tuple import DomainTuple
from ..operators.contraction_operator import ContractionOperator
from ..domains.power_space import PowerSpace
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape
from ..operators.simple_linear_operators import ducktape, VdotOperator
from ..operators.scaling_operator import ScalingOperator
from ..operators.contraction_operator import ContractionOperator
from ..operators.operator_normalizer import OperatorNormalizer
from ..library.inverse_gamma_operator import InverseGammaOperator
def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
......@@ -73,25 +79,44 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
# 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))
return ht(vol * A * ducktape(h_space, None, name))
def MfCorrelatedField(target, amplitudes, name='xi'):
def MultiCorrelatedField(target, amplitudes, gs_alpha, gs_q, name='xi'):
"""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.
into a (partially) correlated field defined on a DomainTuple with separate
correlation structures per subdomain.
Uses of this operator include modeling diffuse emissions with correlation
structure in the spatial and energy dimensions and point sources with
correlation structure in the energy and time dimensions, but not in the
spatial dimension. Generalizations to arbitray combinations of domains
with and without correlation structure are possible.
The field is defined as
`MultiCorrelatedField` = HT( A * xi * volume_factor * global_scaling )
This operator may be used as a model for multi-frequency reconstructions
with a correlation structure in both spatial and energy direction.
and takes a list of Amplitude operators from which to construct A
and the parameters of the global scaling prior, provided by a
:class:`InverseGammaOperator`.
Parameters
----------
target : Domain, DomainTuple or tuple of Domain
Target of the operator. Must contain exactly two spaces.
amplitudes: iterable of Operator
List of two amplitude operators.
Target of the operator.
amplitudes: List of Operators
List of amplitude operators. Must contain one per target subdomain.
If a subdomain without correlations is wanted,
pass `None` for this subdomain.
name : string
:class:`MultiField` key for xi-field.
gs_alpha : float
`alpha` parameter of the global_scaling prior
See :class:`InverseGammaOperator` for interpretation.
gs_q : float
`q` parameter of the global_scaling prior
See :class:`InverseGammaOperator` for interpretation.
Returns
-------
......@@ -107,28 +132,80 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
two (one for the energy and one for the spatial subdomain)
"""
tgt = DomainTuple.make(target)
if len(tgt) != 2:
if len(tgt) < 2:
raise ValueError
if len(amplitudes) != 2:
if len(tgt) != len(amplitudes):
raise ValueError
n_subdoms = len(tgt)
hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
ht = ht2 @ ht1
psp = [aa.target[0] for aa in amplitudes]
pd0 = PowerDistributor(hsp, psp[0], 0)
pd1 = PowerDistributor(pd0.domain, psp[1], 1)
pd = pd0 @ pd1
dd0 = ContractionOperator(pd.domain, 1).adjoint
dd1 = ContractionOperator(pd.domain, 0).adjoint
d = [dd0, dd1]
a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
a = reduce(mul, a)
# Assemble PowerDistributor stack
def gen_pd(pd_target, i):
if amplitudes[i] is not None:
return PowerDistributor(pd_target, amplitudes[i].target[0], i)
else:
return PowerDistributor(pd_target, PowerSpace(pd_target[i]), i)
pd = gen_pd(hsp, 0)
for i in range(1, n_subdoms - 1):
pd = pd @ gen_pd(pd.domain, i)
pd = pd @ gen_pd(pd.domain, n_subdoms - 1)
# Spread amplitudes into their respective constant dimensions
# and calculate their outer product
a_i = []
fct = 1.
for i in range(n_subdoms):
# Assemble spreader stack for this amplitude
spr = ScalingOperator(1., pd.domain)
for j in range(n_subdoms):
if j < i:
spr = spr @ ContractionOperator(spr.domain, 0).adjoint
if j > i:
spr = spr @ ContractionOperator(spr.domain, 1).adjoint
# Normalize amplitudes and apply spreader stack
if amplitudes[i] is not None:
# Dimension with correlation structure
normed_amplitude = OperatorNormalizer(amplitudes[i], p=2)
a_i.append(spr @ normed_amplitude)
else:
# No correlation structure in this dimension
# -> constant amplitude
# Still got to normalize, tough:
# Since the field contains only ones, 2-norm = sqrt(1-norm)
amplitude_integral = np.sqrt(Field.full(spr.domain, 1.).integrate())
fct /= amplitude_integral
fct = ScalingOperator(fct, pd.domain)
a = fct @ reduce(mul, a_i)
# Assemble global power spectrum
A = pd @ a
# For `vol` see comment in `CorrelatedField`
# Volume factors
# One (pixel volume)**(-0.5) included, for explanaition see
# respective comment in `CorrelatedField`.
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ducktape(hsp, None, name))
# Generate global power spectrum scaling operator
vdot_ig = VdotOperator(Field.full(A.target, 1.))
inv_gamma_op = InverseGammaOperator(vdot_ig.adjoint.domain, gs_alpha, gs_q)
global_amplitude = vdot_ig.adjoint @ inv_gamma_op.sqrt()
A_times_xi = vol * A * ducktape(hsp, None, name) * \
global_amplitude.ducktape(name + '_pspec_scaling')
# Assemble global HT
ht = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
for i in range(1, n_subdoms - 1):
htn = HarmonicTransformOperator(ht.target, target=tgt[i], space=i)
ht = htn @ ht
htn = HarmonicTransformOperator(ht.target,
target=tgt[-1],
space=n_subdoms - 1)
ht = htn @ ht
return ht(A_times_xi)
......@@ -25,7 +25,9 @@ from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..operators.simple_linear_operators import VdotOperator
from ..sugar import makeOp
from ..library.inverse_gamma_operator import InverseGammaOperator
def _ceps_kernel(k, a, k0):
......@@ -74,9 +76,9 @@ def CepstrumOperator(target, a, k0):
Larger values result in a weaker constraining prior.
"""
a = float(a)
target = DomainTuple.make(target)
if a <= 0:
raise ValueError
target = DomainTuple.make(target)
if len(target) > 1 or target[0].harmonic:
raise TypeError
if isinstance(k0, (float, int)):
......@@ -85,16 +87,16 @@ def CepstrumOperator(target, a, k0):
k0 = np.array(k0)
if len(k0) != len(target.shape):
raise ValueError
if np.any(np.array(k0) <= 0):
if np.any(k0 <= 0):
raise ValueError
qht = QHTOperator(target)
dom = qht.domain[0]
sym = SymmetrizingOperator(target)
# Compute cepstrum field
dim = len(dom.shape)
dom = qht.domain[0]
shape = dom.shape
dim = len(shape)
q_array = dom.get_k_array()
# Fill all non-zero modes
no_zero_modes = (slice(1, None),)*dim
......@@ -112,7 +114,8 @@ def CepstrumOperator(target, a, k0):
return sym @ qht @ makeOp(cepstrum.sqrt())
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
keys=['tau', 'phi', 'zm']):
'''Operator for parametrizing smooth amplitudes (square roots of power
spectra).
......@@ -130,6 +133,7 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
the output of the generated operator is:
AmplitudeOperator = 0.5*(smooth_component + linear_component)
+ zero_mode
This is then exponentiated and exponentially binned (in this order).
......@@ -161,6 +165,12 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
LogRGSpace (see :class:`ExpTransform`).
iv : float
Prior standard deviation of y-intercept of power law.
zm : float or None
Prior mean value of the zero mode.
If None is given, zero mode is set to zero.
zv : float or None
Prior variance of the zero mode.
If None is given, zero mode is set to zero.
Returns
-------
......@@ -169,12 +179,13 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
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()
return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0,
sm=sm, sv=sv, im=im, iv=iv, zm=zm, zv=zv,
keys=keys).exp()
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
keys=['tau', 'phi']):
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, zm, zv,
keys=['tau', 'phi', 'zm']):
'''LinearOperator for parametrizing smooth log-amplitudes (square roots of
power spectra).
......@@ -186,26 +197,50 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
a, k0 = float(a), float(k0)
sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
if sv <= 0 or iv <= 0:
if sv <= 0 or iv < 0:
raise ValueError
if zm is not None and zv is not None:
zero_mode = True
zm, zv = float(zm), float(zv)
if zm <= 0 or zv <= 0:
raise ValueError
if len(keys) < 3:
raise ValueError("Need a key for the zero mode domain")
else:
zero_mode = False
et = ExpTransform(target, n_pix)
dom = et.domain[0]
# Smooth component
dct = {'a': a, 'k0': k0}
smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
# Linear component
sl = SlopeOperator(dom)
mean = np.array([sm, im + sm*dom.t_0[0]])
sig = np.array([sv, iv])
mean = Field.from_global_data(sl.domain, mean)
sig = Field.from_global_data(sl.domain, sig)
linear = sl @ Adder(mean) @ makeOp(sig).ducktape(keys[1])
sigma = np.array([sv, iv])
sigma = Field.from_global_data(sl.domain, sigma)
linear = sl @ Adder(mean) @ makeOp(sigma).ducktape(keys[1])
# Smooth deviations from linear component
dct = {'a': a, 'k0': k0}
zm_mask = np.ones(dom.shape)
zm_mask[(0,) * len(dom.shape)] = 0.
zm_mask = Field.from_global_data(dom, zm_mask)
zm_mask_op = makeOp(zm_mask)
smooth = zm_mask_op @ CepstrumOperator(dom, **dct).ducktape(keys[0])
# Combine linear and smooth component
# Combine components
loglog_ampl = 0.5*(smooth + linear)
# Add optional zero mode prior
if zero_mode:
zm_mask_inv = Field.full(dom, 1.) - zm_mask
vdot = VdotOperator(zm_mask_inv)
b = zm**2 / zv + 1.
log_ig = InverseGammaOperator(vdot.target, alpha=b+1., q=zm*b).log()
zero_mode = vdot.adjoint @ log_ig.ducktape(keys[2])
loglog_ampl = loglog_ampl + zero_mode
# Go from loglog-space to linear-linear-space
return et @ loglog_ampl
# 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-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..field import Field
from .simple_linear_operators import VdotOperator, WeightApplier
def OperatorNormalizer(op, p=1):
"""Wraps `op` into a shell that ensures the fields returned by
(shell @ op) will have a p-norm of one.
Implements an operator analogously to
OperatorNormalizer(op) = op / |op|_p
where the integration is implemented as an operator chain.
Parameters
----------
op : Operator
p : int
Which p-norm to normalize under
"""
ones_field = Field.full(op.target, 1.)
vdot = VdotOperator(ones_field)
wa = WeightApplier(op.target, spaces=None, power=1)
# integrate as an operator
op_integrated_inv = vdot.adjoint @ (vdot @ wa @ (op)**p)**(-1/p)
return op * op_integrated_inv
......@@ -90,7 +90,7 @@ def testModelLibrary(space, seed):
np.random.seed(seed)
domain = ift.PowerSpace(space.get_default_codomain())
model = ift.SLAmplitude(target=domain, n_pix=4, a=.5, k0=2, sm=3, sv=1.5,
im=1.75, iv=1.3)
im=1.75, iv=1.3, zm=1., zv=0.5)
assert_(isinstance(model, ift.Operator))
S = ift.ScalingOperator(1., model.domain)
pos = S.draw_sample()
......@@ -102,7 +102,9 @@ def testModelLibrary(space, seed):
ift.extra.check_jacobian_consistency(model2, pos, ntries=20)
domtup = ift.DomainTuple.make((space, space))
model3 = ift.MfCorrelatedField(domtup, [model, model])
model3 = ift.MultiCorrelatedField(domtup, [model, model],
gs_alpha=0.5,
gs_q=2.0)
S = ift.ScalingOperator(1., model3.domain)
pos = S.draw_sample()
ift.extra.check_jacobian_consistency(model3, pos, ntries=20)
......
# ToDo Development
## LinearSLAmplitude:
- zero mode as log(InverseGammaOperator)
- use uninformative prior
## Demo:
- plot spectrum samples
## Tests:
- Statistics of fields drawn from sub-spectrum are the same as
statistics of fields drawn from global spectrum, sliced to the
sub-dimension
# Finally
- reintroduce tests for the MF model