Commit 5036a58e authored by Gordian Edenhofer's avatar Gordian Edenhofer
Browse files

Set zm in correlated field via dedicated function

parent 9947fa6d
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Notebook showcasing the NIFTy 6 Correlated Field model # Notebook showcasing the NIFTy 6 Correlated Field model
**Skip to `Parameter Showcases` for the meat/veggies ;)** **Skip to `Parameter Showcases` for the meat/veggies ;)**
The field model roughly works like this: The field model roughly works like this:
`f = HT( A * zero_mode * xi ) + offset` `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. `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 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. a representation of the field in harmonic space.
It is then transformed into the target real space and a offset added. 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 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. and an integrated Wiener process whose amplitude and roughness can be set.
## Setup code ## Setup code
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import nifty7 as ift import nifty7 as ift
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
ift.random.push_sseq_from_seed(43) ift.random.push_sseq_from_seed(43)
n_pix = 256 n_pix = 256
x_space = ift.RGSpace(n_pix) x_space = ift.RGSpace(n_pix)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Plotting routine # Plotting routine
def plot(fields, spectra, title=None): def plot(fields, spectra, title=None):
# Plotting preparation is normally handled by nifty7.Plot # Plotting preparation is normally handled by nifty7.Plot
# It is done manually here to be able to tweak details # It is done manually here to be able to tweak details
# Fields are assumed to have identical domains # Fields are assumed to have identical domains
fig = plt.figure(tight_layout=True, figsize=(12, 3.5)) fig = plt.figure(tight_layout=True, figsize=(12, 3.5))
if title is not None: if title is not None:
fig.suptitle(title, fontsize=14) fig.suptitle(title, fontsize=14)
# Field # Field
ax1 = fig.add_subplot(1, 2, 1) ax1 = fig.add_subplot(1, 2, 1)
ax1.axhline(y=0., color='k', linestyle='--', alpha=0.25) ax1.axhline(y=0., color='k', linestyle='--', alpha=0.25)
for field in fields: for field in fields:
dom = field.domain[0] dom = field.domain[0]
xcoord = np.arange(dom.shape[0]) * dom.distances[0] xcoord = np.arange(dom.shape[0]) * dom.distances[0]
ax1.plot(xcoord, field.val) ax1.plot(xcoord, field.val)
ax1.set_xlim(xcoord[0], xcoord[-1]) ax1.set_xlim(xcoord[0], xcoord[-1])
ax1.set_ylim(-5., 5.) ax1.set_ylim(-5., 5.)
ax1.set_xlabel('x') ax1.set_xlabel('x')
ax1.set_ylabel('f(x)') ax1.set_ylabel('f(x)')
ax1.set_title('Field realizations') ax1.set_title('Field realizations')
# Spectrum # Spectrum
ax2 = fig.add_subplot(1, 2, 2) ax2 = fig.add_subplot(1, 2, 2)
for spectrum in spectra: for spectrum in spectra:
xcoord = spectrum.domain[0].k_lengths xcoord = spectrum.domain[0].k_lengths
ycoord = spectrum.val_rw() ycoord = spectrum.val_rw()
ycoord[0] = ycoord[1] ycoord[0] = ycoord[1]
ax2.plot(xcoord, ycoord) ax2.plot(xcoord, ycoord)
ax2.set_ylim(1e-6, 10.) ax2.set_ylim(1e-6, 10.)
ax2.set_xscale('log') ax2.set_xscale('log')
ax2.set_yscale('log') ax2.set_yscale('log')
ax2.set_xlabel('k') ax2.set_xlabel('k')
ax2.set_ylabel('p(k)') ax2.set_ylabel('p(k)')
ax2.set_title('Power Spectrum') ax2.set_title('Power Spectrum')
fig.align_labels() fig.align_labels()
plt.show() plt.show()
# Helper: draw main sample # Helper: draw main sample
main_sample = None main_sample = None
def init_model(m_pars, fl_pars): def init_model(m_pars, fl_pars):
global main_sample global main_sample
offset_mean = m_pars["offset_mean"] cf = ift.CorrelatedFieldMaker(m_pars["prefix"])
offset_std = m_pars["offset_std"] cf.set_amplitude_total_offset(m_pars["offset_mean"], m_pars["offset_std"])
cf = ift.CorrelatedFieldMaker.make(m_pars["prefix"])
cf.add_fluctuations(**fl_pars) cf.add_fluctuations(**fl_pars)
field = cf.finalize(offset_mean, offset_std, prior_info=0) field = cf.finalize(prior_info=0)
main_sample = ift.from_random(field.domain) main_sample = ift.from_random(field.domain)
print("model domain keys:", field.domain.keys()) print("model domain keys:", field.domain.keys())
# Helper: field and spectrum from parameter dictionaries + plotting # Helper: field and spectrum from parameter dictionaries + plotting
def eval_model(m_pars, fl_pars, title=None, samples=None): def eval_model(m_pars, fl_pars, title=None, samples=None):
offset_mean = m_pars["offset_mean"] cf = ift.CorrelatedFieldMaker(m_pars["prefix"])
offset_std = m_pars["offset_std"] cf.set_amplitude_total_offset(m_pars["offset_mean"], m_pars["offset_std"])
cf = ift.CorrelatedFieldMaker.make(m_pars["prefix"])
cf.add_fluctuations(**fl_pars) cf.add_fluctuations(**fl_pars)
field = cf.finalize(offset_mean, offset_std, prior_info=0) field = cf.finalize(prior_info=0)
spectrum = cf.amplitude spectrum = cf.amplitude
if samples is None: if samples is None:
samples = [main_sample] samples = [main_sample]
field_realizations = [field(s) for s in samples] field_realizations = [field(s) for s in samples]
spectrum_realizations = [spectrum.force(s) for s in samples] spectrum_realizations = [spectrum.force(s) for s in samples]
plot(field_realizations, spectrum_realizations, title) plot(field_realizations, spectrum_realizations, title)
def gen_samples(key_to_vary): def gen_samples(key_to_vary):
if key_to_vary is None: if key_to_vary is None:
return [main_sample] return [main_sample]
dct = main_sample.to_dict() dct = main_sample.to_dict()
subdom_to_vary = dct.pop(key_to_vary).domain subdom_to_vary = dct.pop(key_to_vary).domain
samples = [] samples = []
for i in range(8): for i in range(8):
d = dct.copy() d = dct.copy()
d[key_to_vary] = ift.from_random(subdom_to_vary) d[key_to_vary] = ift.from_random(subdom_to_vary)
samples.append(ift.MultiField.from_dict(d)) samples.append(ift.MultiField.from_dict(d))
return samples return samples
def vary_parameter(parameter_key, values, samples_vary_in=None): def vary_parameter(parameter_key, values, samples_vary_in=None):
s = gen_samples(samples_vary_in) s = gen_samples(samples_vary_in)
for v in values: for v in values:
if parameter_key in cf_make_pars.keys(): if parameter_key in cf_make_pars.keys():
m_pars = {**cf_make_pars, parameter_key: v} m_pars = {**cf_make_pars, parameter_key: v}
eval_model(m_pars, cf_x_fluct_pars, f"{parameter_key} = {v}", s) eval_model(m_pars, cf_x_fluct_pars, f"{parameter_key} = {v}", s)
else: else:
fl_pars = {**cf_x_fluct_pars, parameter_key: v} fl_pars = {**cf_x_fluct_pars, parameter_key: v}
eval_model(cf_make_pars, fl_pars, f"{parameter_key} = {v}", s) eval_model(cf_make_pars, fl_pars, f"{parameter_key} = {v}", s)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Before the Action: The Moment-Matched Log-Normal Distribution ## Before the Action: The Moment-Matched Log-Normal Distribution
Many properties of the correlated field are modelled as being lognormally distributed. 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). 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, 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!) here are a few example histograms: (observe the x-axis!)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(13, 3.5)) fig = plt.figure(figsize=(13, 3.5))
mean = 1.0 mean = 1.0
sigmas = [1.0, 0.5, 0.1] sigmas = [1.0, 0.5, 0.1]
for i in range(3): for i in range(3):
op = ift.LognormalTransform(mean=mean, sigma=sigmas[i], op = ift.LognormalTransform(mean=mean, sigma=sigmas[i],
key='foo', N_copies=0) key='foo', N_copies=0)
op_samples = np.array( op_samples = np.array(
[op(s).val for s in [ift.from_random(op.domain) for i in range(10000)]]) [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 = fig.add_subplot(1, 3, i + 1)
ax.hist(op_samples, bins=50) ax.hist(op_samples, bins=50)
ax.set_title(f"mean = {mean}, sigma = {sigmas[i]}") ax.set_title(f"mean = {mean}, sigma = {sigmas[i]}")
ax.set_xlabel('x') ax.set_xlabel('x')
del op_samples del op_samples
plt.show() plt.show()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## The Neutral Field ## The Neutral Field
To demonstrate the effect of all parameters, first a 'neutral' set of parameters 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 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 on the generated field realizations and the underlying power spectrum from which
they were drawn. they were drawn.
As a neutral field, a model with a white power spectrum and vanishing spectral power was chosen. As a neutral field, a model with a white power spectrum and vanishing spectral power was chosen.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Neutral model parameters yielding a quasi-constant field # Neutral model parameters yielding a quasi-constant field
cf_make_pars = { cf_make_pars = {
'offset_mean': 0., 'offset_mean': 0.,
'offset_std': (1e-3, 1e-16), 'offset_std': (1e-3, 1e-16),
'prefix': '' 'prefix': ''
} }
cf_x_fluct_pars = { cf_x_fluct_pars = {
'target_subdomain': x_space, 'target_subdomain': x_space,
'fluctuations': (1e-3, 1e-16), 'fluctuations': (1e-3, 1e-16),
'flexibility': (1e-3, 1e-16), 'flexibility': (1e-3, 1e-16),
'asperity': (1e-3, 1e-16), 'asperity': (1e-3, 1e-16),
'loglogavgslope': (0., 1e-16) 'loglogavgslope': (0., 1e-16)
} }
init_model(cf_make_pars, cf_x_fluct_pars) init_model(cf_make_pars, cf_x_fluct_pars)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Show neutral field # Show neutral field
eval_model(cf_make_pars, cf_x_fluct_pars, "Neutral Field") eval_model(cf_make_pars, cf_x_fluct_pars, "Neutral Field")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Parameter Showcases # Parameter Showcases
## The `fluctuations` parameters of `add_fluctuations()` ## The `fluctuations` parameters of `add_fluctuations()`
determine the **amplitude of variations along the field dimension** determine the **amplitude of variations along the field dimension**
for which `add_fluctuations` is called. for which `add_fluctuations` is called.
`fluctuations[0]` set the _average_ amplitude of the fields fluctuations along the given dimension,\ `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. `fluctuations[1]` sets the width and shape of the amplitude distribution.
The amplitude is modelled as being log-normally distributed, The amplitude is modelled as being log-normally distributed,
see `The Moment-Matched Log-Normal Distribution` above for details. see `The Moment-Matched Log-Normal Distribution` above for details.
#### `fluctuations` mean: #### `fluctuations` mean:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('fluctuations', [(0.05, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi') vary_parameter('fluctuations', [(0.05, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
#### `fluctuations` std: #### `fluctuations` std:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('fluctuations', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='fluctuations') vary_parameter('fluctuations', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='fluctuations')
cf_x_fluct_pars['fluctuations'] = (1., 1e-16) cf_x_fluct_pars['fluctuations'] = (1., 1e-16)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## The `loglogavgslope` parameters of `add_fluctuations()` ## The `loglogavgslope` parameters of `add_fluctuations()`
determine **the slope of the loglog-linear (power law) component of the power spectrum**. determine **the slope of the loglog-linear (power law) component of the power spectrum**.
The slope is modelled to be normally distributed. The slope is modelled to be normally distributed.
#### `loglogavgslope` mean: #### `loglogavgslope` mean:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('loglogavgslope', [(-6., 1e-16), (-2., 1e-16), (2., 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: %% Cell type:markdown id: tags:
#### `loglogavgslope` std: #### `loglogavgslope` std:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('loglogavgslope', [(-2., 0.02), (-2., 0.2), (-2., 2.0)], samples_vary_in='loglogavgslope') vary_parameter('loglogavgslope', [(-2., 0.02), (-2., 0.2), (-2., 2.0)], samples_vary_in='loglogavgslope')
cf_x_fluct_pars['loglogavgslope'] = (-2., 1e-16) cf_x_fluct_pars['loglogavgslope'] = (-2., 1e-16)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## The `flexibility` parameters of `add_fluctuations()` ## The `flexibility` parameters of `add_fluctuations()`
determine **the amplitude of the integrated Wiener process component of the power spectrum** determine **the amplitude of the integrated Wiener process component of the power spectrum**
(how strong the power spectrum varies besides the power-law). (how strong the power spectrum varies besides the power-law).
`flexibility[0]` sets the _average_ amplitude of the i.g.p. component,\ `flexibility[0]` sets the _average_ amplitude of the i.g.p. component,\
`flexibility[1]` sets how much the amplitude can vary.\ `flexibility[1]` sets how much the amplitude can vary.\
These two parameters feed into a moment-matched log-normal distribution model, These two parameters feed into a moment-matched log-normal distribution model,
see above for a demo of its behavior. see above for a demo of its behavior.
#### `flexibility` mean: #### `flexibility` mean:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('flexibility', [(0.4, 1e-16), (4.0, 1e-16), (12.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: %% Cell type:markdown id: tags:
#### `flexibility` std: #### `flexibility` std:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('flexibility', [(4., 0.02), (4., 0.2), (4., 2.)], samples_vary_in='flexibility') vary_parameter('flexibility', [(4., 0.02), (4., 0.2), (4., 2.)], samples_vary_in='flexibility')
cf_x_fluct_pars['flexibility'] = (4., 1e-16) cf_x_fluct_pars['flexibility'] = (4., 1e-16)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## The `asperity` parameters of `add_fluctuations()` ## The `asperity` parameters of `add_fluctuations()`
`asperity` determines **how rough the integrated Wiener process component of the power spectrum is**. `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.\ `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, These two parameters feed into a moment-matched log-normal distribution model,
see above for a demo of its behavior. see above for a demo of its behavior.
#### `asperity` mean: #### `asperity` mean:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('asperity', [(0.001, 1e-16), (1.0, 1e-16), (5., 1e-16)], samples_vary_in='spectrum') vary_parameter('asperity', [(0.001, 1e-16), (1.0, 1e-16), (5., 1e-16)], samples_vary_in='spectrum')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
#### `asperity` std: #### `asperity` std:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('asperity', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='asperity') vary_parameter('asperity', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='asperity')
cf_x_fluct_pars['asperity'] = (1., 1e-16) cf_x_fluct_pars['asperity'] = (1., 1e-16)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## The `offset_mean` parameter of `CorrelatedFieldMaker.make()` ## The `offset_mean` parameter of `CorrelatedFieldMaker.make()`
The `offset_mean` parameter defines a global additive offset on the field realizations. 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. 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: %% Cell type:code id: tags:
``` python ``` python
# Reset model to neutral # Reset model to neutral
cf_x_fluct_pars['fluctuations'] = (1e-3, 1e-16) cf_x_fluct_pars['fluctuations'] = (1e-3, 1e-16)
cf_x_fluct_pars['flexibility'] = (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['asperity'] = (1e-3, 1e-16)
cf_x_fluct_pars['loglogavgslope'] = (1e-3, 1e-16) cf_x_fluct_pars['loglogavgslope'] = (1e-3, 1e-16)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('offset_mean', [3., 0., -2.]) vary_parameter('offset_mean', [3., 0., -2.])
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## The `offset_std` parameters of `CorrelatedFieldMaker.make()` ## The `offset_std` parameters of `CorrelatedFieldMaker.make()`
Variation of the global offset of the field are modelled as being log-normally distributed. 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. 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[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. The `offset_std[1]` parameters defines the with and shape of the offset variation distribution.
#### `offset_std` mean: #### `offset_std` mean:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('offset_std', [(1e-16, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi') vary_parameter('offset_std', [(1e-16, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
#### `offset_std` std: #### `offset_std` std:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
vary_parameter('offset_std', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='zeromode') vary_parameter('offset_std', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='zeromode')
``` ```
......
...@@ -73,14 +73,17 @@ def main(): ...@@ -73,14 +73,17 @@ def main():
sp2 = ift.RGSpace(npix2) sp2 = ift.RGSpace(npix2)
# Set up signal model # Set up signal model
cfmaker = ift.CorrelatedFieldMaker.make('') cfmaker = ift.CorrelatedFieldMaker('')
cfmaker.add_fluctuations(sp1, (0.1, 1e-2), (2, .2), (.01, .5), (-4, 2.), 'amp1') cfmaker.add_fluctuations(sp1, (0.1, 1e-2), (2, .2), (.01, .5), (-4, 2.),
cfmaker.add_fluctuations(sp2, (0.1, 1e-2), (2, .2), (.01, .5), 'amp1')
(-3, 1), 'amp2') cfmaker.add_fluctuations(sp2, (0.1, 1e-2), (2, .2), (.01, .5), (-3, 1),
correlated_field = cfmaker.finalize(0., (1e-2, 1e-6)) 'amp2')
cfmaker.set_amplitude_total_offset(0., (1e-2, 1e-6))
pspec1 = cfmaker.normalized_amplitudes[0]**2 correlated_field = cfmaker.finalize()
pspec2 = cfmaker.normalized_amplitudes[1]**2
normalized_amp = cfmaker.normalized_amplitudes()
pspec1 = normalized_amp[0]**2
pspec2 = normalized_amp[1]**2
DC = SingleDomain(correlated_field.target, position_space) DC = SingleDomain(correlated_field.target, position_space)
# Apply a nonlinearity # Apply a nonlinearity
......
...@@ -88,7 +88,6 @@ from .library.adjust_variances import (make_adjust_variances_hamiltonian, ...@@ -88,7 +88,6 @@ from .library.adjust_variances import (make_adjust_variances_hamiltonian,
do_adjust_variances) do_adjust_variances)
from .library.nft import Gridder, FinuFFT from .library.nft import Gridder, FinuFFT
from .library.correlated_fields import CorrelatedFieldMaker from .library.correlated_fields import CorrelatedFieldMaker
from .library.correlated_fields import moment_slice_to_average
from .library.correlated_fields_simple import SimpleCorrelatedField from .library.correlated_fields_simple import SimpleCorrelatedField
from . import extra from . import extra
......
...@@ -24,6 +24,7 @@ from operator import mul ...@@ -24,6 +24,7 @@ from operator import mul
import numpy as np import numpy as np
from .. import utilities from .. import utilities
from ..logger import logger
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from ..domains.power_space import PowerSpace from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain from ..domains.unstructured_domain import UnstructuredDomain
...@@ -85,21 +86,6 @@ def _total_fluctuation_realized(samples): ...@@ -85,21 +86,6 @@ def _total_fluctuation_realized(samples):
return np.sqrt(res if np.isscalar(res) else res.val) return np.sqrt(res if np.isscalar(res) else res.val)
def moment_slice_to_average(amp, zm, fluctuations_slice_mean, nsamples=1000):
fluctuations_slice_mean = float(fluctuations_slice_mean)
if not fluctuations_slice_mean > 0:
msg = "fluctuations_slice_mean must be greater zero; got {!r}"
raise ValueError(msg.format(fluctuations_slice_mean))
from ..sugar import from_random
scm = 1.
for a in amp:
op = a.fluctuation_amplitude*zm.ptw("reciprocal")
res = np.array([op(from_random(op.domain, 'normal')).val
for _ in range(nsamples)])
scm *= res**2 + 1.
return fluctuations_slice_mean/np.mean(np.sqrt(scm))
class _SlopeRemover(EndomorphicOperator): class _SlopeRemover(EndomorphicOperator):
def __init__(self, domain, space=0): def __init__(self, domain, space=0):
self._domain = makeDomain(domain) self._domain = makeDomain(domain)
...@@ -428,26 +414,11 @@ class CorrelatedFieldMaker: ...@@ -428,26 +414,11 @@ class CorrelatedFieldMaker:
See the methods :func:`make`, :func:`add_fluctuations` See the methods :func:`make`, :func:`add_fluctuations`
and :func:`finalize` for further usage information.""" and :func:`finalize` for further usage information."""
def __init__(self, prefix, total_N): def __init__(self, prefix, total_N=0):
self._azm = None """Instantiate a CorrelatedFieldMaker object.
self._offset_mean = None
self._a = []
self._target_subdomains = []
self._prefix = prefix
self._total_N = total_N
@staticmethod
def make(prefix, total_N=0):
"""Returns a CorrelatedFieldMaker object.
Parameters Parameters
---------- ----------
offset_mean : float
Mean offset from zero of the correlated field to be made.
offset_std : tuple of float
Mean standard deviation and standard deviation of the standard
deviation of the offset. No, this is not a word duplication.
prefix : string prefix : string
Prefix to the names of the domains of the cf operator to be made. Prefix to the names of the domains of the cf operator to be made.
This determines the names of the operator domain. This determines the names of the operator domain.
...@@ -456,16 +427,14 @@ class CorrelatedFieldMaker: ...@@ -456,16 +427,14 @@ class CorrelatedFieldMaker:
If not 0, the first entry of the operators target will be an If not 0, the first entry of the operators target will be an
:class:`~nifty.domains.unstructured_domain.UnstructuredDomain` :class:`~nifty.domains.unstructured_domain.UnstructuredDomain`
with length `total_N`. with length `total_N`.
dofdex : np.array of integers, optional
An integer array specifying the zero mode models used if
total_N > 1. It needs to have length of total_N. If total_N=3 and
dofdex=[0,0,1], that means that two models for the zero mode are
instantiated; the first one is used for the first and second
field model and the second is used for the third field model.
*If not specified*, use the same zero mode model for all
constructed field models.
""" """
return CorrelatedFieldMaker(prefix, total_N) self._azm = None
self._offset_mean = None
self._a = []
self._target_subdomains = []
self._prefix = prefix
self._total_N = total_N
def add_fluctuations(self, def add_fluctuations(self,
target_subdomain, target_subdomain,
...@@ -667,7 +636,45 @@ class CorrelatedFieldMaker: ...@@ -667,7 +636,45 @@ class CorrelatedFieldMaker:
self._a.append(amp) self._a.append(amp)
self._target_subdomains.append(target_subdomain) self._target_subdomains.append(target_subdomain)
def finalize(self, offset_mean, offset_std, dofdex=None, prior_info=100): def set_amplitude_total_offset(self, offset_mean, offset_std, dofdex=None):
"""Sets the zero-mode for the combined amplitude operator
Parameters
----------
offset_mean : float
Mean offset from zero of the correlated field to be made.
offset_std : tuple of float
Mean standard deviation and standard deviation of the standard
deviation of the offset. No, this is not a word duplication.
dofdex : np.array of integers, optional
An integer array specifying the zero mode models used if
total_N > 1. It needs to have length of total_N. If total_N=3 and
dofdex=[0,0,1], that means that two models for the zero mode are
instantiated; the first one is used for the first and second
field model and the second is used for the third field model.
*If not specified*, use the same zero mode model for all
constructed field models.
"""
if self._offset_mean is not None and self._azm is not None:
logger.warning("Overwriting the previous mean offset and zero-mode")
self._offset_mean = offset_mean
if isinstance(offset_std, Operator):
self._azm = offset_std
else:
if dofdex is None:
dofdex = np.full(self._total_N, 0)
elif len(dofdex) != self._total_N:
raise ValueError("length of dofdex needs to match total_N")
N = max(dofdex) + 1 if self._total_N > 0 else 0
if len(offset_std) != 2:
raise TypeError
zm = LognormalTransform(*offset_std, self._prefix + 'zeromode', N)
if self._total_N > 0:
zm = _Distributor(dofdex, zm.target, UnstructuredDomain(self._total_N)) @ zm
self._azm = zm
def finalize(self, prior_info=100):
"""Finishes model construction process and returns the constructed """Finishes model construction process and returns the constructed
operator. operator.
...@@ -677,20 +684,6 @@ class CorrelatedFieldMaker: ...@@ -677,20 +684,6 @@ class CorrelatedFieldMaker:
How many prior samples to draw for property verification statistics How many prior samples to draw for property verification statistics
If zero, skips calculating and displaying statistics. If zero, skips calculating and displaying statistics.
""" """
if dofdex is None:
dofdex = np.full(self._total_N, 0)
elif len(dofdex) != self._total_N:
raise ValueError("length of dofdex needs to match total_N")
N = max(dofdex) + 1 if self._total_N > 0 else 0
# TODO: Allow for `offset_std` being an arbitrary operator
if len(offset_std) != 2:
raise TypeError
zm = LognormalTransform(*offset_std, self._prefix + 'zeromode', N)
if self._total_N > 0:
zm = _Distributor(dofdex, zm.target, UnstructuredDomain(self._total_N)) @ zm
self._azm = zm
self._offset_mean = offset_mean
n_amplitudes = len(self._a) n_amplitudes = len(self._a)
if self._total_N > 0: if self._total_N > 0:
hspace = makeDomain( hspace = makeDomain(
...@@ -705,7 +698,7 @@ class CorrelatedFieldMaker: ...@@ -705,7 +698,7 @@ class CorrelatedFieldMaker:
amp_space = 0 amp_space = 0
expander = ContractionOperator(hspace, spaces=spaces).adjoint expander = ContractionOperator(hspace, spaces=spaces).adjoint
azm = expander @ self._azm azm = expander @ self.azm
ht = HarmonicTransformOperator(hspace, ht = HarmonicTransformOperator(hspace,
self._target_subdomains[0][amp_space], self._target_subdomains[0][amp_space],
...@@ -714,26 +707,12 @@ class CorrelatedFieldMaker: ...@@ -714,26 +707,12 @@ class CorrelatedFieldMaker:
ht = HarmonicTransformOperator(ht.target, ht = HarmonicTransformOperator(ht.target,
self._target_subdomains[i][amp_space], self._target_subdomains[i][amp_space],
space=spaces[i]) @ ht space=spaces[i]) @ ht
a = [] a = list(self.normalized_amplitudes())
for ii in range(n_amplitudes): for ii in range(n_amplitudes):
a_target = self._a[ii].target
a_space = 0 if not hasattr(self._a[ii], "_space") else self._a[ii]._space
a_pp = self._a[ii].target[a_space]
assert isinstance(a_pp, PowerSpace)
azm_expander = ContractionOperator(a_target, spaces=a_space).adjoint
zm_unmask, zm_mask = [np.zeros(a_pp.shape) for _ in range(2)]
zm_mask[1:] = zm_unmask[0] = 1.
zm_mask = DiagonalOperator(makeField(a_pp, zm_mask), a_target, a_space)
zm_unmask = DiagonalOperator(makeField(a_pp, zm_unmask), a_target, a_space)
zm_unmask = Adder(zm_unmask(full(zm_unmask.domain, 1)))
zm_normalization = zm_unmask @ (zm_mask @ azm_expander(self._azm.ptw("reciprocal")))
co = ContractionOperator(hspace, spaces[:ii] + spaces[ii + 1:]) co = ContractionOperator(hspace, spaces[:ii] + spaces[ii + 1:])
pp = self._a[ii].target[amp_space] pp = a[ii].target[amp_space]
pd = PowerDistributor(co.target, pp, amp_space) pd = PowerDistributor(co.target, pp, amp_space)
a.append(co.adjoint @ pd @ (zm_normalization * self._a[ii])) a[ii] = co.adjoint @ pd @ a[ii]
corr = reduce(mul, a) corr = reduce(mul, a)
op = ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi')) op = ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
...@@ -775,11 +754,43 @@ class CorrelatedFieldMaker: ...@@ -775,11 +754,43 @@ class CorrelatedFieldMaker:
logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s)) logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
@property @property
def normalized_amplitudes(self): def fluctuations(self):
# TODO: rename """Returns the added fluctuations operators used in the model"""
"""Returns the amplitude operators used in the model"""
return self._a return self._a
def normalized_amplitudes(self):
"""Returns the normalized amplitude operators used in the final model
The amplitude operators are corrected for the otherwise degenerate
zero-mode. Their scales are only meaningful relative to one another.
Their absolute scale bares no information.
"""
normal_amp = []
for amp in self._a:
a_target = amp.target
a_space = 0 if not hasattr(amp, "_space") else amp._space
a_pp = amp.target[a_space]
assert isinstance(a_pp, PowerSpace)
azm_expander = ContractionOperator(
a_target, spaces=a_space
).adjoint
zm_unmask, zm_mask = [np.zeros(a_pp.shape) for _ in range(2)]
zm_mask[1:] = zm_unmask[0] = 1.
zm_mask = DiagonalOperator(
makeField(a_pp, zm_mask), a_target, a_space
)
zm_unmask = DiagonalOperator(
makeField(a_pp, zm_unmask), a_target, a_space
)
zm_unmask = Adder(zm_unmask(full(zm_unmask.domain, 1)))
zm_normalization = zm_unmask @ (
zm_mask @ azm_expander(self.azm.ptw("reciprocal"))
)
normal_amp.append(zm_normalization * amp)
return tuple(normal_amp)
@property @property
def amplitude(self): def amplitude(self):
if len(self._a) > 1: if len(self._a) > 1:
...@@ -787,22 +798,12 @@ class CorrelatedFieldMaker: ...@@ -787,22 +798,12 @@ class CorrelatedFieldMaker:
' no unique set of amplitudes exist because only the', ' no unique set of amplitudes exist because only the',
' relative scale is determined.') ' relative scale is determined.')
raise NotImplementedError(s) raise NotImplementedError(s)
a_target = self._a[0].target normal_amp = self.normalized_amplitudes()[0]
a_space = 0 if not hasattr(self._a[0], "_space") else self._a[0]._space
a_pp = self._a[0].target[a_space]
assert isinstance(a_pp, PowerSpace)
azm_expander = ContractionOperator(a_target, spaces=a_space).adjoint
zm_unmask, zm_mask = [np.zeros(a_pp.shape) for _ in range(2)]
zm_mask[1:] = zm_unmask[0] = 1.
zm_mask = DiagonalOperator(makeField(a_pp, zm_mask), a_target, a_space)
zm_unmask = DiagonalOperator(makeField(a_pp, zm_unmask), a_target, a_space)
zm_unmask = Adder(zm_unmask(full(zm_unmask.domain, 1)))
zm_normalization = zm_unmask @ (zm_mask @ azm_expander(self._azm.ptw("reciprocal"))) expand = ContractionOperator(
normal_amp.target, len(normal_amp.target) - 1
expand = ContractionOperator(a_target, len(a_target) - 1).adjoint ).adjoint
return (zm_normalization * self._a[0]) * (expand @ self.amplitude_total_offset) return normal_amp * (expand @ self.azm)
@property @property
def power_spectrum(self): def power_spectrum(self):
...@@ -810,8 +811,31 @@ class CorrelatedFieldMaker: ...@@ -810,8 +811,31 @@ class CorrelatedFieldMaker:
@property @property
def amplitude_total_offset(self): def amplitude_total_offset(self):
"""Returns the total offset of the amplitudes"""
if self._azm is None:
nie = "You need to set the `amplitude_total_offset` first"
raise NotImplementedError(nie)
return self._azm return self._azm
@property
def azm(self):
"""Alias for `amplitude_total_offset`"""
return self.amplitude_total_offset
def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
fluctuations_slice_mean = float(fluctuations_slice_mean)
if not fluctuations_slice_mean > 0:
msg = "fluctuations_slice_mean must be greater zero; got {!r}"