Skip to content

GitLab

Commit 43761bb6 authored by Philipp Arras's avatar Philipp Arras
Browse files

Simplify correlated field model and fixups 

A = a / (sqrt(sum(a**2)))

-->

A = sqrt(pspec / sqrt(sum(pspec)))
parent 3e42a35d
Hide whitespace changes
Inline Side-by-side
ChangeLog.md
View file @ 43761bb6
......@@ -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
---------------------
......
demos/getting_started_3.py
View file @ 43761bb6
......@@ -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'])
......
demos/getting_started_4_CorrelatedFields.ipynb
View file @ 43761bb6
%% Cell type:markdown id: tags:
# Notebook showcasing the NIFTy 6 Correlated Field model
**Skip to `Parameter Showcases` for the meat/veggies ;)**
The field model roughly works like this:
`f = HT( A * zero_mode * xi ) + offset`
`A` is a spectral power field which is constructed from power spectra that hold on subdomains of the target domain.
It is scaled by a zero mode operator and then pointwise multiplied by a gaussian excitation field, yielding
a representation of the field in harmonic space.
It is then transformed into the target real space and a offset added.
The power spectra `A` is constructed of are in turn constructed as the sum of a power law component
and an integrated Wiener process whose amplitude and roughness can be set.
## Setup code
%% Cell type:code id: tags:
``` python
import nifty7 as ift
import matplotlib.pyplot as plt
import numpy as np
ift.random.push_sseq_from_seed(43)
n_pix = 256
x_space = ift.RGSpace(n_pix)
```
%% Cell type:code id: tags:
``` python
# Plotting routine
def plot(fields, spectra, title=None):
# Plotting preparation is normally handled by nifty7.Plot
# It is done manually here to be able to tweak details
# Fields are assumed to have identical domains
fig = plt.figure(tight_layout=True, figsize=(12, 3.5))
if title is not None:
fig.suptitle(title, fontsize=14)
# Field
ax1 = fig.add_subplot(1, 2, 1)
ax1.axhline(y=0., color='k', linestyle='--', alpha=0.25)
for field in fields:
dom = field.domain[0]
xcoord = np.arange(dom.shape[0]) * dom.distances[0]
ax1.plot(xcoord, field.val)
ax1.set_xlim(xcoord[0], xcoord[-1])
ax1.set_ylim(-5., 5.)
ax1.set_xlabel('x')
ax1.set_ylabel('f(x)')
ax1.set_title('Field realizations')
# Spectrum
ax2 = fig.add_subplot(1, 2, 2)
for spectrum in spectra:
xcoord = spectrum.domain[0].k_lengths
ycoord = spectrum.val_rw()
ycoord[0] = ycoord[1]
ax2.plot(xcoord, ycoord)
ax2.set_ylim(1e-6, 10.)
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.set_xlabel('k')
ax2.set_ylabel('p(k)')
ax2.set_title('Power Spectrum')
fig.align_labels()
plt.show()
# Helper: draw main sample
main_sample = None
def init_model(m_pars, fl_pars):
global main_sample
cf = ift.CorrelatedFieldMaker.make(**m_pars)
cf.add_fluctuations(**fl_pars)
field = cf.finalize(prior_info=0)
main_sample = ift.from_random(field.domain)
print("model domain keys:", field.domain.keys())
# Helper: field and spectrum from parameter dictionaries + plotting
def eval_model(m_pars, fl_pars, title=None, samples=None):
cf = ift.CorrelatedFieldMaker.make(**m_pars)
cf.add_fluctuations(**fl_pars)
field = cf.finalize(prior_info=0)
spectrum = cf.amplitude
if samples is None:
samples = [main_sample]
field_realizations = [field(s) for s in samples]
spectrum_realizations = [spectrum.force(s) for s in samples]
plot(field_realizations, spectrum_realizations, title)
def gen_samples(key_to_vary):
if key_to_vary is None:
return [main_sample]
dct = main_sample.to_dict()
subdom_to_vary = dct.pop(key_to_vary).domain
samples = []
for i in range(8):
d = dct.copy()
d[key_to_vary] = ift.from_random(subdom_to_vary)
samples.append(ift.MultiField.from_dict(d))
return samples
def vary_parameter(parameter_key, values, samples_vary_in=None):
s = gen_samples(samples_vary_in)
for v in values:
if parameter_key in cf_make_pars.keys():
m_pars = {**cf_make_pars, parameter_key: v}
eval_model(m_pars, cf_x_fluct_pars, f"{parameter_key} = {v}", s)
else:
fl_pars = {**cf_x_fluct_pars, parameter_key: v}
eval_model(cf_make_pars, fl_pars, f"{parameter_key} = {v}", s)
```
%% Cell type:markdown id: tags:
## Before the Action: The Moment-Matched Log-Normal Distribution
Many properties of the correlated field are modelled as being lognormally distributed.
The distribution models are parametrized via their means and standard-deviations (first and second position in tuple).
To get a feeling of how the ratio of the `mean` and `stddev` parameters influences the distribution shape,
here are a few example histograms: (observe the x-axis!)
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(13, 3.5))
mean = 1.0
sigmas = [1.0, 0.5, 0.1]
for i in range(3):
op = ift.LognormalTransform(mean=mean, sigma=sigmas[i],
key='foo', N_copies=0)
op_samples = np.array(
[op(s).val for s in [ift.from_random(op.domain) for i in range(10000)]])
ax = fig.add_subplot(1, 3, i + 1)
ax.hist(op_samples, bins=50)
ax.set_title(f"mean = {mean}, sigma = {sigmas[i]}")
ax.set_xlabel('x')
del op_samples
plt.show()
```
%% Cell type:markdown id: tags:
## The Neutral Field
To demonstrate the effect of all parameters, first a 'neutral' set of parameters
is defined which then are varied one by one, showing the effect of the variation
on the generated field realizations and the underlying power spectrum from which
they were drawn.
As a neutral field, a model with a white power spectrum and vanishing spectral power was chosen.
%% Cell type:code id: tags:
``` python
# Neutral model parameters yielding a quasi-constant field
cf_make_pars = {
'offset_mean': 0.,
'offset_std': (1e-3, 1e-16),
'prefix': ''
}
cf_x_fluct_pars = {
'target_subdomain': x_space,
'fluctuations': (1e-3, 1e-16),
'flexibility': (1e-3, 1e-16),
'asperity': (1e-3, 1e-16),
'loglogavgslope': (0., 1e-16)
}
init_model(cf_make_pars, cf_x_fluct_pars)
```
%% Cell type:code id: tags:
``` python
# Show neutral field
eval_model(cf_make_pars, cf_x_fluct_pars, "Neutral Field")
```
%% Cell type:markdown id: tags:
# Parameter Showcases
## The `fluctuations` parameters of `add_fluctuations()`
determine the **amplitude of variations along the field dimension**
for which `add_fluctuations` is called.
`fluctuations[0]` set the _average_ amplitude of the fields fluctuations along the given dimension,\
`fluctuations[1]` sets the width and shape of the amplitude distribution.
The amplitude is modelled as being log-normally distributed,
see `The Moment-Matched Log-Normal Distribution` above for details.
#### `fluctuations` mean:
%% Cell type:code id: tags:
``` python
vary_parameter('fluctuations', [(0.05, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi')
```
%% Cell type:markdown id: tags:
#### `fluctuations` std:
%% Cell type:code id: tags:
``` python
vary_parameter('fluctuations', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='fluctuations')
cf_x_fluct_pars['fluctuations'] = (1., 1e-16)
```
%% Cell type:markdown id: tags:
## The `loglogavgslope` parameters of `add_fluctuations()`
determine **the slope of the loglog-linear (power law) component of the power spectrum**.
The slope is modelled to be normally distributed.
#### `loglogavgslope` mean:
%% Cell type:code id: tags:
``` python
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')
```
%% Cell type:markdown id: tags:
#### `loglogavgslope` std:
%% Cell type:code id: tags:
``` python
vary_parameter('loglogavgslope', [(-1., 0.01), (-1., 0.1), (-1., 1.0)], samples_vary_in='loglogavgslope')
cf_x_fluct_pars['loglogavgslope'] = (-1., 1e-16)
vary_parameter('loglogavgslope', [(-2., 0.02), (-2., 0.2), (-2., 2.0)], samples_vary_in='loglogavgslope')
cf_x_fluct_pars['loglogavgslope'] = (-2., 1e-16)
```
%% Cell type:markdown id: tags:
## The `flexibility` parameters of `add_fluctuations()`
determine **the amplitude of the integrated Wiener process component of the power spectrum**
(how strong the power spectrum varies besides the power-law).
`flexibility[0]` sets the _average_ amplitude of the i.g.p. component,\
`flexibility[1]` sets how much the amplitude can vary.\
These two parameters feed into a moment-matched log-normal distribution model,
see above for a demo of its behavior.
#### `flexibility` mean:
%% Cell type:code id: tags:
``` python
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')
```
%% Cell type:markdown id: tags:
#### `flexibility` std:
%% Cell type:code id: tags:
``` python
vary_parameter('flexibility', [(2., 0.01), (2., 0.1), (2., 1.)], samples_vary_in='flexibility')
cf_x_fluct_pars['flexibility'] = (2., 1e-16)
vary_parameter('flexibility', [(4., 0.02), (4., 0.2), (4., 2.)], samples_vary_in='flexibility')
cf_x_fluct_pars['flexibility'] = (4., 1e-16)
```
%% Cell type:markdown id: tags:
## The `asperity` parameters of `add_fluctuations()`
`asperity` determines **how rough the integrated Wiener process component of the power spectrum is**.
`asperity[0]` sets the average roughness, `asperity[1]` sets how much the roughness can vary.\
These two parameters feed into a moment-matched log-normal distribution model,
see above for a demo of its behavior.
#### `asperity` mean:
%% Cell type:code id: tags:
``` python
vary_parameter('asperity', [(0.001, 1e-16), (1.0, 1e-16), (5., 1e-16)], samples_vary_in='spectrum')
```
%% Cell type:markdown id: tags:
#### `asperity` std:
%% Cell type:code id: tags:
``` python
vary_parameter('asperity', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='asperity')
cf_x_fluct_pars['asperity'] = (1., 1e-16)
```
%% Cell type:markdown id: tags:
## The `offset_mean` parameter of `CorrelatedFieldMaker.make()`
The `offset_mean` parameter defines a global additive offset on the field realizations.
If the field is used for a lognormal model `f = field.exp()`, this acts as a global signal magnitude offset.
%% Cell type:code id: tags:
``` python
# Reset model to neutral
cf_x_fluct_pars['fluctuations'] = (1e-3, 1e-16)
cf_x_fluct_pars['flexibility'] = (1e-3, 1e-16)
cf_x_fluct_pars['asperity'] = (1e-3, 1e-16)
cf_x_fluct_pars['loglogavgslope'] = (1e-3, 1e-16)
```
%% Cell type:code id: tags:
``` python
vary_parameter('offset_mean', [3., 0., -2.])
```
%% Cell type:markdown id: tags:
## The `offset_std` parameters of `CorrelatedFieldMaker.make()`
Variation of the global offset of the field are modelled as being log-normally distributed.
See `The Moment-Matched Log-Normal Distribution` above for details.
The `offset_std[0]` parameter sets how much NIFTy will vary the offset *on average*.\
The `offset_std[1]` parameters defines the with and shape of the offset variation distribution.
#### `offset_std` mean:
%% Cell type:code id: tags:
``` python
vary_parameter('offset_std', [(1e-16, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi')
```
%% Cell type:markdown id: tags:
#### `offset_std` std:
%% Cell type:code id: tags:
``` python
vary_parameter('offset_std', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='zeromode')
```
......
demos/getting_started_5_mf.py
View file @ 43761bb6
......@@ -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)
......
src/library/correlated_fields.py
View file @ 43761bb6
......@@ -168,11 +168,10 @@ class _Normalization(Operator):
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(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
......@@ -606,7 +604,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
......@@ -620,6 +618,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
......
src/library/correlated_fields_simple.py
View file @ 43761bb6
......@@ -114,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