Commit 46ece3be authored by Philipp Arras's avatar Philipp Arras

Merge branch 'simplifymodel' into 'NIFTy_7'

Amplitude -> Power spectrum

See merge request !558
parents 69313966 43761bb6
Pipeline #78770 passed with stages
in 13 minutes and 16 seconds
......@@ -13,6 +13,12 @@ Furthermore, it is now possible to disable the asperity and the flexibility
together with the asperity in the correlated field model. Note that disabling
only the flexibility is not possible.
Additionally, the parameters `flexibility`, `asperity` and most importantly
`loglogavgslope` refer to the power spectrum instead of the amplitude now.
For existing codes that means that both values in the tuple `loglogavgslope`
and `flexibility` need to be doubled. The transformation of the `asperity`
parameter is nontrivial.
SimpleCorrelatedField
---------------------
......
......@@ -63,20 +63,20 @@ def main():
'offset_std': (1e-3, 1e-6),
# Amplitude of field fluctuations
'fluctuations': (2., 1.), # 1.0, 1e-2
'fluctuations': (2., 1.), # 1.0, 1e-2
# Exponent of power law power spectrum component
'loglogavgslope': (-2., 0.5), # -3.0, 0.5
'loglogavgslope': (-4., 1), # -6.0, 1
# Amplitude of integrated Wiener process power spectrum component
'flexibility': (2.5, 1.), # 1.0, 0.5
'flexibility': (5, 2.), # 2.0, 1.0
# How ragged the integrated Wiener process component is
'asperity': (0.5, 0.5) # 0.1, 0.5
}
correlated_field = ift.SimpleCorrelatedField(position_space, **args)
A = correlated_field.amplitude
pspec = correlated_field.power_spectrum
# Apply a nonlinearity
signal = ift.sigmoid(correlated_field)
......@@ -113,7 +113,7 @@ def main():
plot = ift.Plot()
plot.add(signal(mock_position), title='Ground Truth')
plot.add(R.adjoint_times(data), title='Data')
plot.add([A.force(mock_position)], title='Power Spectrum')
plot.add([pspec.force(mock_position)], title='Power Spectrum')
plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename.format("setup"))
# number of samples used to estimate the KL
......@@ -129,7 +129,7 @@ def main():
# Plot current reconstruction
plot = ift.Plot()
plot.add(signal(KL.position), title="Latent mean")
plot.add([A.force(KL.position + ss) for ss in KL.samples],
plot.add([pspec.force(KL.position + ss) for ss in KL.samples],
title="Samples power spectrum")
plot.output(ny=1, ysize=6, xsize=16,
name=filename.format("loop_{:02d}".format(i)))
......@@ -144,13 +144,13 @@ def main():
plot.add(sc.mean, title="Posterior Mean")
plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
powers = [A.force(s + KL.position) for s in KL.samples]
powers = [pspec.force(s + KL.position) for s in KL.samples]
sc = ift.StatCalculator()
for pp in powers:
sc.add(pp)
plot.add(
powers + [A.force(mock_position),
A.force(KL.position), sc.mean],
powers + [pspec.force(mock_position),
pspec.force(KL.position), sc.mean],
title="Sampled Posterior Power Spectrum",
linewidth=[1.]*len(powers) + [3., 3., 3.],
label=[None]*len(powers) + ['Ground truth', 'Posterior latent mean', 'Posterior mean'])
......
......@@ -279,7 +279,7 @@
"metadata": {},
"outputs": [],
"source": [
"vary_parameter('loglogavgslope', [(-3., 1e-16), (-1., 1e-16), (1., 1e-16)], samples_vary_in='xi')"
"vary_parameter('loglogavgslope', [(-6., 1e-16), (-2., 1e-16), (2., 1e-16)], samples_vary_in='xi')"
]
},
{
......@@ -295,8 +295,8 @@
"metadata": {},
"outputs": [],
"source": [
"vary_parameter('loglogavgslope', [(-1., 0.01), (-1., 0.1), (-1., 1.0)], samples_vary_in='loglogavgslope')\n",
"cf_x_fluct_pars['loglogavgslope'] = (-1., 1e-16)"
"vary_parameter('loglogavgslope', [(-2., 0.02), (-2., 0.2), (-2., 2.0)], samples_vary_in='loglogavgslope')\n",
"cf_x_fluct_pars['loglogavgslope'] = (-2., 1e-16)"
]
},
{
......@@ -322,7 +322,7 @@
"metadata": {},
"outputs": [],
"source": [
"vary_parameter('flexibility', [(0.1, 1e-16), (1.0, 1e-16), (3.0, 1e-16)], samples_vary_in='spectrum')"
"vary_parameter('flexibility', [(0.4, 1e-16), (4.0, 1e-16), (12.0, 1e-16)], samples_vary_in='spectrum')"
]
},
{
......@@ -338,8 +338,8 @@
"metadata": {},
"outputs": [],
"source": [
"vary_parameter('flexibility', [(2., 0.01), (2., 0.1), (2., 1.)], samples_vary_in='flexibility')\n",
"cf_x_fluct_pars['flexibility'] = (2., 1e-16)"
"vary_parameter('flexibility', [(4., 0.02), (4., 0.2), (4., 2.)], samples_vary_in='flexibility')\n",
"cf_x_fluct_pars['flexibility'] = (4., 1e-16)"
]
},
{
......
......@@ -74,13 +74,13 @@ def main():
# Set up signal model
cfmaker = ift.CorrelatedFieldMaker.make(0., (1e-2, 1e-6), '')
cfmaker.add_fluctuations(sp1, (0.1, 1e-2), (1, .1), (.01, .5), (-2, 1.), 'amp1')
cfmaker.add_fluctuations(sp2, (0.1, 1e-2), (1, .1), (.01, .5),
(-1.5, .5), 'amp2')
cfmaker.add_fluctuations(sp1, (0.1, 1e-2), (2, .2), (.01, .5), (-4, 2.), 'amp1')
cfmaker.add_fluctuations(sp2, (0.1, 1e-2), (2, .2), (.01, .5),
(-3, 1), 'amp2')
correlated_field = cfmaker.finalize()
A1 = cfmaker.normalized_amplitudes[0]
A2 = cfmaker.normalized_amplitudes[1]
pspec1 = cfmaker.normalized_amplitudes[0]**2
pspec2 = cfmaker.normalized_amplitudes[1]**2
DC = SingleDomain(correlated_field.target, position_space)
# Apply a nonlinearity
......@@ -103,8 +103,8 @@ def main():
plot = ift.Plot()
plot.add(signal(mock_position), title='Ground Truth')
plot.add(R.adjoint_times(data), title='Data')
plot.add([A1.force(mock_position)], title='Power Spectrum 1')
plot.add([A2.force(mock_position)], title='Power Spectrum 2')
plot.add([pspec1.force(mock_position)], title='Power Spectrum 1')
plot.add([pspec2.force(mock_position)], title='Power Spectrum 2')
plot.output(ny=2, nx=2, xsize=10, ysize=10, name=filename.format("setup"))
# Minimization parameters
......@@ -139,11 +139,11 @@ def main():
plot = ift.Plot()
plot.add(signal(mock_position), title="ground truth")
plot.add(signal(KL.position), title="reconstruction")
plot.add([A1.force(KL.position),
A1.force(mock_position)],
plot.add([pspec1.force(KL.position),
pspec1.force(mock_position)],
title="power1")
plot.add([A2.force(KL.position),
A2.force(mock_position)],
plot.add([pspec2.force(KL.position),
pspec2.force(mock_position)],
title="power2")
plot.add((ic_newton.history, ic_sampling.history,
minimizer.inversion_history),
......@@ -165,8 +165,8 @@ def main():
powers2 = []
for sample in KL.samples:
sc.add(signal(sample + KL.position))
p1 = A1.force(sample + KL.position)
p2 = A2.force(sample + KL.position)
p1 = pspec1.force(sample + KL.position)
p2 = pspec2.force(sample + KL.position)
scA1.add(p1)
powers1.append(p1)
scA2.add(p2)
......@@ -178,12 +178,12 @@ def main():
plot.add(sc.mean, title="Posterior Mean")
plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
powers1 = [A1.force(s + KL.position) for s in KL.samples]
powers2 = [A2.force(s + KL.position) for s in KL.samples]
plot.add(powers1 + [scA1.mean, A1.force(mock_position)],
powers1 = [pspec1.force(s + KL.position) for s in KL.samples]
powers2 = [pspec2.force(s + KL.position) for s in KL.samples]
plot.add(powers1 + [scA1.mean, pspec1.force(mock_position)],
title="Sampled Posterior Power Spectrum 1",
linewidth=[1.]*len(powers1) + [3., 3.])
plot.add(powers2 + [scA2.mean, A2.force(mock_position)],
plot.add(powers2 + [scA2.mean, pspec2.force(mock_position)],
title="Sampled Posterior Power Spectrum 2",
linewidth=[1.]*len(powers2) + [3., 3.])
plot.output(ny=2, nx=2, xsize=15, ysize=15, name=filename_res)
......
......@@ -88,23 +88,23 @@ class _SlopeRemover(EndomorphicOperator):
self._domain = makeDomain(domain)
assert isinstance(self._domain[space], PowerSpace)
logkl = _relative_log_k_lengths(self._domain[space])
self._sc = logkl/float(logkl[-1])
sc = logkl/float(logkl[-1])
self._space = space
axis = self._domain.axes[space][0]
self._last = (slice(None),)*axis + (-1,) + (None,)
self._extender = (None,)*(axis) + (slice(None),) + (None,)*(self._domain.axes[-1][-1]-axis)
extender = (None,)*(axis) + (slice(None),) + (None,)*(self._domain.axes[-1][-1]-axis)
self._sc = sc[extender]
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
x = x.val
if mode == self.TIMES:
res = x - x[self._last]*self._sc[self._extender]
x = x.val
res = x - x[self._last]*self._sc
else:
res = x.copy()
res[self._last] -= (x*self._sc[self._extender]).sum(
axis=self._space, keepdims=True)
res = x.val_rw()
res[self._last] -= (res*self._sc).sum(axis=self._space, keepdims=True)
return makeField(self._tgt(mode), res)
......@@ -163,16 +163,15 @@ class _Normalization(Operator):
mode_multiplicity = pd.adjoint(full(pd.target, 1.)).val_rw()
zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,)
mode_multiplicity[zero_mode] = 0
self._mode_multiplicity = makeOp(makeField(self._domain, mode_multiplicity))
self._specsum = _SpecialSum(self._domain, space)
multipl = makeOp(makeField(self._domain, mode_multiplicity))
self._specsum = _SpecialSum(self._domain, space) @ multipl
def apply(self, x):
self._check_input(x)
amp = x.ptw("exp")
spec = amp**2
spec = x.ptw("exp")
# FIXME This normalizes also the zeromode which is supposed to be left
# untouched by this operator
return self._specsum(self._mode_multiplicity(spec))**(-0.5)*amp
return (self._specsum(spec).reciprocal()*spec).sqrt()
class _SpecialSum(EndomorphicOperator):
......@@ -412,12 +411,11 @@ class CorrelatedFieldMaker:
on which they apply.
The parameters `fluctuations`, `flexibility`, `asperity` and
`loglogavgslope` configure the power spectrum model ("amplitude")
used on the target field subdomain `target_subdomain`.
It is assembled as the sum of a power law component
(linear slope in log-log power-frequency-space),
a smooth varying component (integrated Wiener process) and
a ragged component (un-integrated Wiener process).
`loglogavgslope` configure the power spectrum model used on the target
field subdomain `target_subdomain`. It is assembled as the sum of a
power law component (linear slope in log-log power-frequency-space), a
smooth varying component (integrated Wiener process) and a ragged
component (un-integrated Wiener process).
Multiple calls to `add_fluctuations` are possible, in which case
the constructed field will have the outer product of the individual
......@@ -610,7 +608,7 @@ class CorrelatedFieldMaker:
@property
def normalized_amplitudes(self):
"""Returns the power spectrum operators used in the model"""
"""Returns the amplitude operators used in the model"""
return self._a
@property
......@@ -624,6 +622,10 @@ class CorrelatedFieldMaker:
expand = ContractionOperator(dom, len(dom)-1).adjoint
return self._a[0]*(expand @ self.amplitude_total_offset)
@property
def power_spectrum(self):
return self.amplitude**2
@property
def amplitude_total_offset(self):
return self._azm
......
......@@ -71,21 +71,23 @@ class SimpleCorrelatedField(Operator):
if flexibility is not None:
flex = LognormalTransform(*flexibility, prefix + 'flexibility', 0)
dom = twolog.domain[0]
vflex = np.zeros(dom.shape)
vflex = np.empty(dom.shape)
vflex[0] = vflex[1] = np.sqrt(_log_vol(pspace))
vflex = makeOp(makeField(dom, vflex))
sig_flex = vflex @ expander @ flex
xi = ducktape(dom, None, prefix + 'spectrum')
shift = np.ones(dom.shape)
shift = np.empty(dom.shape)
shift[0] = _log_vol(pspace)**2 / 12.
shift[1] = 1
shift = makeField(dom, shift)
if asperity is None:
asp = makeOp(shift.ptw("sqrt")) @ (xi*sig_flex)
else:
asp = LognormalTransform(*asperity, prefix + 'asperity', 0)
vasp = np.zeros(dom.shape)
vasp = np.empty(dom.shape)
vasp[0] = 1
vasp[1] = 0
vasp = makeOp(makeField(dom, vasp))
sig_asp = vasp @ expander @ asp
asp = xi*sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt")
......@@ -112,3 +114,8 @@ class SimpleCorrelatedField(Operator):
def amplitude(self):
"""Analoguous to :func:`~nifty7.library.correlated_fields.CorrelatedFieldMaker.amplitude`."""
return self._a
@property
def power_spectrum(self):
"""Analoguous to :func:`~nifty7.library.correlated_fields.CorrelatedFieldMaker.power_spectrum`."""
return self.amplitude**2
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