Commit d8b20aad authored by Gordian Edenhofer's avatar Gordian Edenhofer
Browse files

correlated_fields.py: No un-normalized amplitudes

Adapt the density_estimator accordingly.
parent ddfde472
......@@ -32,18 +32,21 @@ import nifty7 as ift
def density_estimator(
domain, exposure=1., pad=1., cf_fluctuations_sane_defaults=None
domain, exposure=1., pad=1., cf_fluctuations=None, cf_azm=None
):
from functools import reduce
cf_azm_sane_default = (0., (1e-2, 1e-6))
cf_fluctuations_sane_default = {
"scale": (1.0, 0.5),
"cutoff": (10.0, 5.0),
"loglogslope": (-12.0, 6.0)
}
domain = ift.DomainTuple.make(domain)
dom_scaling = 1. + np.broadcast_to(pad, (len(domain.axes), ))
if cf_fluctuations_sane_defaults is None:
cf_fluctuations_sane_defaults = {
"scale": (1.0, 0.5),
"cutoff": (10.0, 5.0),
"loglogslope": (-12.0, 6.0)
}
if cf_fluctuations is None:
cf_fluctuations = cf_fluctuations_sane_default
if cf_azm is None:
cf_azm = cf_azm_sane_default
domain_padded = []
for d_scl, d in zip(dom_scaling, domain):
......@@ -63,11 +66,13 @@ def density_estimator(
prefix = "de" # density estimator
cfmaker = ift.CorrelatedFieldMaker(prefix)
for i, d in enumerate(domain_padded):
cfmaker.add_fluctuations_matern(
d, **cf_fluctuations_sane_defaults, prefix=f"ax{i}"
)
# cfmaker.set_amplitude_total_offset(0., (1e-2, 1e-6))
correlated_field = cfmaker.finalize(0, normalize=False)
if isinstance(cf_fluctuations, (list, tuple)):
cf_fl = cf_fluctuations[i]
else:
cf_fl = cf_fluctuations
cfmaker.add_fluctuations_matern(d, **cf_fl, prefix=f"ax{i}")
cfmaker.set_amplitude_total_offset(*cf_azm)
correlated_field = cfmaker.finalize(0)
domain_shape = tuple(d.shape for d in domain)
slc = ift.SliceOperator(correlated_field.target, domain_shape)
......
......@@ -676,7 +676,7 @@ class CorrelatedFieldMaker:
zm = _Distributor(dofdex, zm.target, UnstructuredDomain(self._total_N)) @ zm
self._azm = zm
def finalize(self, prior_info=100, normalize=True):
def finalize(self, prior_info=100):
"""Finishes model construction process and returns the constructed
operator.
......@@ -707,36 +707,18 @@ class CorrelatedFieldMaker:
self._target_subdomains[i][amp_space],
space=spaces[i]) @ ht
if normalize:
expander = ContractionOperator(hspace, spaces=spaces).adjoint
azm = expander @ self.azm
expander = ContractionOperator(hspace, spaces=spaces).adjoint
azm = expander @ self.azm
a = list(self.get_normalized_amplitudes())
else:
a = list(self.fluctuations)
a = list(self.get_normalized_amplitudes())
for ii in range(n_amplitudes):
co = ContractionOperator(hspace, spaces[:ii] + spaces[ii + 1:])
pp = a[ii].target[amp_space]
pd = PowerDistributor(co.target, pp, amp_space)
a[ii] = co.adjoint @ pd @ a[ii]
corr = reduce(mul, a)
if normalize:
xi = ducktape(hspace, None, self._prefix + 'xi')
op = ht(azm * corr * xi)
else:
if self._total_N > 0:
re = (
"initializing a correlated field without an amplitude total"
" offset with `total_N` larger than zero is not supported"
)
raise RuntimeError(re)
m = np.zeros(hspace.shape, dtype=bool)
m.flat[0] = True
mask_zm = MaskOperator(Field(hspace, m))
xi = ducktape(mask_zm.adjoint.domain, None, self._prefix + "xi")
xi = mask_zm.adjoint @ xi
op = ht(corr * xi)
xi = ducktape(hspace, None, self._prefix + 'xi')
op = ht(azm * corr * xi)
if self._offset_mean is not None:
offset = self._offset_mean
......
Supports Markdown
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