diff --git a/demos/getting_started_4_CorrelatedFields.ipynb b/demos/getting_started_4_CorrelatedFields.ipynb
index a969bddd88908fc6e405ca718db0c073f2194375..cf28e26ca4ea5e86698068856bdac6bb2b00af62 100644
--- a/demos/getting_started_4_CorrelatedFields.ipynb
+++ b/demos/getting_started_4_CorrelatedFields.ipynb
@@ -87,21 +87,19 @@
"main_sample = None\n",
"def init_model(m_pars, fl_pars):\n",
" global main_sample\n",
- " offset_mean = m_pars[\"offset_mean\"]\n",
- " offset_std = m_pars[\"offset_std\"]\n",
- " cf = ift.CorrelatedFieldMaker.make(m_pars[\"prefix\"])\n",
+ " cf = ift.CorrelatedFieldMaker(m_pars[\"prefix\"])\n",
+ " cf.set_amplitude_total_offset(m_pars[\"offset_mean\"], m_pars[\"offset_std\"])\n",
" cf.add_fluctuations(**fl_pars)\n",
- " field = cf.finalize(offset_mean, offset_std, prior_info=0)\n",
+ " field = cf.finalize(prior_info=0)\n",
" main_sample = ift.from_random(field.domain)\n",
" print(\"model domain keys:\", field.domain.keys())\n",
" \n",
"# Helper: field and spectrum from parameter dictionaries + plotting\n",
"def eval_model(m_pars, fl_pars, title=None, samples=None):\n",
- " offset_mean = m_pars[\"offset_mean\"]\n",
- " offset_std = m_pars[\"offset_std\"]\n",
- " cf = ift.CorrelatedFieldMaker.make(m_pars[\"prefix\"])\n",
+ " cf = ift.CorrelatedFieldMaker(m_pars[\"prefix\"])\n",
+ " cf.set_amplitude_total_offset(m_pars[\"offset_mean\"], m_pars[\"offset_std\"])\n",
" cf.add_fluctuations(**fl_pars)\n",
- " field = cf.finalize(offset_mean, offset_std, prior_info=0)\n",
+ " field = cf.finalize(prior_info=0)\n",
" spectrum = cf.amplitude\n",
" if samples is None:\n",
" samples = [main_sample]\n",
@@ -466,6 +464,18 @@
"display_name": "Python 3",
"language": "python",
"name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.2"
}
},
"nbformat": 4,
diff --git a/demos/getting_started_5_mf.py b/demos/getting_started_5_mf.py
index 3a2dbe213c2d1b12fd9d97934fee4d205ec893e1..e6bf3c2a67dc688a4bc78394ccbe65aecc042922 100644
--- a/demos/getting_started_5_mf.py
+++ b/demos/getting_started_5_mf.py
@@ -73,14 +73,17 @@ def main():
sp2 = ift.RGSpace(npix2)
# Set up signal model
- cfmaker = ift.CorrelatedFieldMaker.make('')
- 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(0., (1e-2, 1e-6))
-
- pspec1 = cfmaker.normalized_amplitudes[0]**2
- pspec2 = cfmaker.normalized_amplitudes[1]**2
+ cfmaker = ift.CorrelatedFieldMaker('')
+ 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')
+ cfmaker.set_amplitude_total_offset(0., (1e-2, 1e-6))
+ correlated_field = cfmaker.finalize()
+
+ normalized_amp = cfmaker.normalized_amplitudes()
+ pspec1 = normalized_amp[0]**2
+ pspec2 = normalized_amp[1]**2
DC = SingleDomain(correlated_field.target, position_space)
# Apply a nonlinearity
diff --git a/src/__init__.py b/src/__init__.py
index 827ed3c267817bc30e542d9402532156aec6adf8..ba33e119d986be40801ed2e121cf454bb1fdf094 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -88,7 +88,6 @@ from .library.adjust_variances import (make_adjust_variances_hamiltonian,
do_adjust_variances)
from .library.nft import Gridder, FinuFFT
from .library.correlated_fields import CorrelatedFieldMaker
-from .library.correlated_fields import moment_slice_to_average
from .library.correlated_fields_simple import SimpleCorrelatedField
from . import extra
diff --git a/src/library/correlated_fields.py b/src/library/correlated_fields.py
index b2673dc013e2006c07fb83b177a6b4e4ddeb3474..3183525be5843799850216627370f06b451dc63e 100644
--- a/src/library/correlated_fields.py
+++ b/src/library/correlated_fields.py
@@ -24,6 +24,7 @@ from operator import mul
import numpy as np
from .. import utilities
+from ..logger import logger
from ..domain_tuple import DomainTuple
from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
@@ -85,21 +86,6 @@ def _total_fluctuation_realized(samples):
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):
def __init__(self, domain, space=0):
self._domain = makeDomain(domain)
@@ -428,26 +414,11 @@ class CorrelatedFieldMaker:
See the methods :func:`make`, :func:`add_fluctuations`
and :func:`finalize` for further usage information."""
- def __init__(self, prefix, total_N):
- self._azm = None
- 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.
+ def __init__(self, prefix, total_N=0):
+ """Instantiate a CorrelatedFieldMaker object.
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 to the names of the domains of the cf operator to be made.
This determines the names of the operator domain.
@@ -456,16 +427,14 @@ class CorrelatedFieldMaker:
If not 0, the first entry of the operators target will be an
:class:`~nifty.domains.unstructured_domain.UnstructuredDomain`
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,
target_subdomain,
@@ -667,7 +636,45 @@ class CorrelatedFieldMaker:
self._a.append(amp)
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
operator.
@@ -677,20 +684,6 @@ class CorrelatedFieldMaker:
How many prior samples to draw for property verification 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)
if self._total_N > 0:
hspace = makeDomain(
@@ -705,7 +698,7 @@ class CorrelatedFieldMaker:
amp_space = 0
expander = ContractionOperator(hspace, spaces=spaces).adjoint
- azm = expander @ self._azm
+ azm = expander @ self.azm
ht = HarmonicTransformOperator(hspace,
self._target_subdomains[0][amp_space],
@@ -714,26 +707,12 @@ class CorrelatedFieldMaker:
ht = HarmonicTransformOperator(ht.target,
self._target_subdomains[i][amp_space],
space=spaces[i]) @ ht
- a = []
+ a = list(self.normalized_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:])
- pp = self._a[ii].target[amp_space]
+ pp = a[ii].target[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)
op = ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
@@ -775,11 +754,43 @@ class CorrelatedFieldMaker:
logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
@property
- def normalized_amplitudes(self):
- # TODO: rename
- """Returns the amplitude operators used in the model"""
+ def fluctuations(self):
+ """Returns the added fluctuations operators used in the model"""
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
def amplitude(self):
if len(self._a) > 1:
@@ -787,22 +798,12 @@ class CorrelatedFieldMaker:
' no unique set of amplitudes exist because only the',
' relative scale is determined.')
raise NotImplementedError(s)
- a_target = self._a[0].target
- 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)))
+ normal_amp = self.normalized_amplitudes()[0]
- zm_normalization = zm_unmask @ (zm_mask @ azm_expander(self._azm.ptw("reciprocal")))
-
- expand = ContractionOperator(a_target, len(a_target) - 1).adjoint
- return (zm_normalization * self._a[0]) * (expand @ self.amplitude_total_offset)
+ expand = ContractionOperator(
+ normal_amp.target, len(normal_amp.target) - 1
+ ).adjoint
+ return normal_amp * (expand @ self.azm)
@property
def power_spectrum(self):
@@ -810,8 +811,31 @@ class CorrelatedFieldMaker:
@property
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
+ @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}"
+ raise ValueError(msg.format(fluctuations_slice_mean))
+ from ..sugar import from_random
+ scm = 1.
+ for a in self._a:
+ op = a.fluctuation_amplitude * self.azm.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))
+
@property
def total_fluctuation(self):
"""Returns operator which acts on prior or posterior samples"""
@@ -821,9 +845,9 @@ class CorrelatedFieldMaker:
return self.average_fluctuation(0)
q = 1.
for a in self._a:
- fl = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
+ fl = a.fluctuation_amplitude*self.azm.ptw("reciprocal")
q = q*(Adder(full(fl.target, 1.)) @ fl**2)
- return (Adder(full(q.target, -1.)) @ q).ptw("sqrt")*self._azm
+ return (Adder(full(q.target, -1.)) @ q).ptw("sqrt")*self.azm
def slice_fluctuation(self, space):
"""Returns operator which acts on prior or posterior samples"""
@@ -835,12 +859,12 @@ class CorrelatedFieldMaker:
return self.average_fluctuation(0)
q = 1.
for j in range(len(self._a)):
- fl = self._a[j].fluctuation_amplitude*self._azm.ptw("reciprocal")
+ fl = self._a[j].fluctuation_amplitude*self.azm.ptw("reciprocal")
if j == space:
q = q*fl**2
else:
q = q*(Adder(full(fl.target, 1.)) @ fl**2)
- return q.ptw("sqrt")*self._azm
+ return q.ptw("sqrt")*self.azm
def average_fluctuation(self, space):
"""Returns operator which acts on prior or posterior samples"""
diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py
index 8eba29c8384bba1126a6cdb3f3609e6beda1a012..ba4c1f350ab4eecd9e74734a1200692966c1b6d2 100644
--- a/test/test_operators/test_correlated_fields.py
+++ b/test/test_operators/test_correlated_fields.py
@@ -71,7 +71,7 @@ def test_init(total_N, asperity, flexibility, ind, matern):
pytest.skip()
cfg = 1, 1
for dofdex in ([None], [None, [0]], [None, [0, 0], [0, 1], [1, 1]])[total_N]:
- cfm = ift.CorrelatedFieldMaker.make('', total_N)
+ cfm = ift.CorrelatedFieldMaker('', total_N)
cfm.add_fluctuations(ift.RGSpace(4), cfg, flexibility, asperity, (-2, 0.1))
if matern:
if total_N == 0:
@@ -81,7 +81,8 @@ def test_init(total_N, asperity, flexibility, ind, matern):
cfm.add_fluctuations_matern(ift.RGSpace(4), *(3*[cfg]))
else:
cfm.add_fluctuations(ift.RGSpace(4), *(4*[cfg]), index=ind)
- cfm.finalize(0, cfg, dofdex=dofdex, prior_info=0)
+ cfm.set_amplitude_total_offset(0, cfg, dofdex=dofdex)
+ cfm.finalize(prior_info=0)
@pmp('sspace', spaces)
@@ -94,12 +95,13 @@ def testAmplitudesInvariants(sspace, N):
astds = 0.2, 1.2
offset_std_mean = 1.3
- fa = ift.CorrelatedFieldMaker.make('', N)
+ fa = ift.CorrelatedFieldMaker('', N)
fa.add_fluctuations(sspace, (astds[0], 1e-2), (1.1, 2.), (2.1, .5), (-2, 1.),
'spatial', dofdex=dofdex2)
fa.add_fluctuations(fsspace, (astds[1], 1e-2), (3.1, 1.), (.5, .1), (-4, 1.),
'freq', dofdex=dofdex3)
- op = fa.finalize(1.2, (offset_std_mean, 1e-2), dofdex=dofdex1, prior_info=0)
+ fa.set_amplitude_total_offset(1.2, (offset_std_mean, 1e-2), dofdex=dofdex1)
+ op = fa.finalize(prior_info=0)
samples = [ift.from_random(op.domain) for _ in range(100)]
tot_flm, _ = _stats(fa.total_fluctuation, samples)
@@ -125,25 +127,22 @@ def testAmplitudesInvariants(sspace, N):
assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5)
assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5)
- fa = ift.CorrelatedFieldMaker.make('', N)
+ fa = ift.CorrelatedFieldMaker('', N)
+ fa.set_amplitude_total_offset(0., (offset_std_mean, .1), dofdex=dofdex1)
fa.add_fluctuations(fsspace, (astds[1], 1.), (3.1, 1.), (.5, .1), (-4, 1.), 'freq',
dofdex=dofdex3)
m = 3.
- n_zm = max(dofdex1) + 1 if N > 0 else 0
- zm = ift.LognormalTransform(*(offset_std_mean, .1), 'zeromode', n_zm)
- if N > 0:
- zm = ift.library.correlated_fields._Distributor(dofdex1, zm.target, ift.UnstructuredDomain(N)) @ zm
- x = ift.moment_slice_to_average(fa.normalized_amplitudes, zm, m)
+ x = fa.moment_slice_to_average(m)
fa.add_fluctuations(sspace, (x, 1.5), (1.1, 2.), (2.1, .5), (-2, 1.), 'spatial', 0,
dofdex=dofdex2)
- op = fa.finalize(0., (offset_std_mean, .1), dofdex=dofdex1, prior_info=0)
+ op = fa.finalize(prior_info=0)
em, estd = _stats(fa.slice_fluctuation(0), samples)
assert_allclose(m, em, rtol=0.5)
assert op.target[-2] == sspace
assert op.target[-1] == fsspace
- for ampl in fa.normalized_amplitudes:
+ for ampl in fa.normalized_amplitudes():
ift.extra.check_operator(ampl, 0.1*ift.from_random(ampl.domain), ntries=10)
ift.extra.check_operator(op, 0.1*ift.from_random(op.domain), ntries=10)
@@ -183,7 +182,7 @@ def test_complicated_vs_simple(seed, domain, without):
asperity,
loglogavgslope
)
- cfm = ift.CorrelatedFieldMaker.make(prefix)
+ cfm = ift.CorrelatedFieldMaker(prefix)
if asperity is not None and flexibility is None:
with pytest.raises(ValueError):
scf = ift.SimpleCorrelatedField(*scf_args, prefix=prefix,
@@ -196,8 +195,9 @@ def test_complicated_vs_simple(seed, domain, without):
harmonic_partner=hspace)
cfm.add_fluctuations(*add_fluct_args, prefix='',
harmonic_partner=hspace)
+ cfm.set_amplitude_total_offset(offset_mean, offset_std)
inp = ift.from_random(scf.domain)
- op1 = cfm.finalize(offset_mean, offset_std, prior_info=0)
+ op1 = cfm.finalize(prior_info=0)
assert scf.domain is op1.domain
ift.extra.assert_allclose(scf(inp), op1(inp))
ift.extra.check_operator(scf, inp, ntries=10)