Commit 0e644a88 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'support_fixed_zeromode' into 'NIFTy_7'

correlated_fields_simple.py: Support fixed zero mode

See merge request !618
parents 0eff4f8e 8371ff35
Pipeline #101733 passed with stages
in 26 minutes and 25 seconds
......@@ -59,11 +59,11 @@ class GLSpace(StructuredDomain):
@property
def shape(self):
return (np.int((self.nlat * self.nlon)),)
return (int((self.nlat * self.nlon)),)
@property
def size(self):
return np.int((self.nlat * self.nlon))
return int((self.nlat * self.nlon))
@property
def scalar_dvol(self):
......
......@@ -53,7 +53,7 @@ class HPSpace(StructuredDomain):
@property
def size(self):
return np.int(12 * self.nside * self.nside)
return int(12 * self.nside * self.nside)
@property
def scalar_dvol(self):
......
......@@ -131,7 +131,7 @@ class LMSpace(StructuredDomain):
op = HarmonicTransformOperator(lm0, gl)
kernel_lm = op.adjoint_times(kernel_sphere.weight(1)).val
# evaluate the k lengths of the harmonic space
k_lengths = self.get_k_length_array().val.astype(np.int)
k_lengths = self.get_k_length_array().val.astype(int)
return Field.from_raw(self, kernel_lm[k_lengths])
@property
......
......@@ -259,7 +259,7 @@ class Field(Operator):
if np.isscalar(wgt):
fct *= wgt
else:
new_shape = np.ones(len(self.shape), dtype=np.int)
new_shape = np.ones(len(self.shape), dtype=int)
new_shape[self._domain.axes[ind][0]:
self._domain.axes[ind][-1]+1] = wgt.shape
wgt = wgt.reshape(new_shape)
......
......@@ -642,9 +642,15 @@ class CorrelatedFieldMaker:
----------
offset_mean : float
Mean offset from zero of the correlated field to be made.
offset_std : tuple of float or instance of :class:`~nifty7.operators.operator.Operator` acting on scalar domain
offset_std : tuple of float, instance of \
:class:`~nifty7.operators.operator.Operator` acting on scalar \
domain, scalar or None
Mean standard deviation and standard deviation of the standard
deviation of the offset. No, this is not a word duplication.
The option to specify `None` or equivalently a scalar value of `1.`
only really makes sense for one dimensional amplitude spectral.
Take special care if using this option for multi-dimensional
amplitude spectra that this is really what you want.
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
......@@ -658,7 +664,9 @@ class CorrelatedFieldMaker:
logger.warning("Overwriting the previous mean offset and zero-mode")
self._offset_mean = offset_mean
if isinstance(offset_std, Operator):
if offset_std is None or (np.isscalar(offset_std) and offset_std == 1.):
self._azm = 1.
elif isinstance(offset_std, Operator):
self._azm = offset_std
else:
if dofdex is None:
......@@ -667,7 +675,11 @@ class CorrelatedFieldMaker:
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
te = (
"`offset_std` of invalid type and/or shape"
f"; expected a 2D tuple of floats; got '{offset_std!r}'"
)
raise TypeError(te)
zm = LognormalTransform(*offset_std, self._prefix + 'zeromode', N)
if self._total_N > 0:
zm = _Distributor(dofdex, zm.target, UnstructuredDomain(self._total_N)) @ zm
......@@ -704,10 +716,10 @@ class CorrelatedFieldMaker:
self._target_subdomains[i][amp_space],
space=spaces[i]) @ ht
expander = ContractionOperator(hspace, spaces=spaces).adjoint
azm = expander @ self.azm
a = list(self.get_normalized_amplitudes())
if np.isscalar(self.azm):
a = list(self.fluctuations)
else:
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]
......@@ -715,7 +727,12 @@ class CorrelatedFieldMaker:
a[ii] = co.adjoint @ pd @ a[ii]
corr = reduce(mul, a)
xi = ducktape(hspace, None, self._prefix + 'xi')
op = ht(azm * corr * xi)
if np.isscalar(self.azm):
op = ht(corr * xi)
else:
expander = ContractionOperator(hspace, spaces=spaces).adjoint
azm = expander @ self.azm
op = ht(azm * corr * xi)
if self._offset_mean is not None:
offset = self._offset_mean
......@@ -772,7 +789,15 @@ class CorrelatedFieldMaker:
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.
Notes
-----
In the case of no zero-mode, i.e. an assumed zero-mode of unity, this
call is equivalent to the `fluctuations` property.
"""
if np.isscalar(self.azm):
return self.fluctuations
normal_amp = []
for amp in self._a:
a_target = amp.target
......@@ -808,10 +833,13 @@ class CorrelatedFieldMaker:
raise NotImplementedError(s)
normal_amp = self.get_normalized_amplitudes()[0]
expand = ContractionOperator(
normal_amp.target, len(normal_amp.target) - 1
).adjoint
return normal_amp * (expand @ self.azm)
if np.isscalar(self.azm):
return normal_amp
else:
expand = ContractionOperator(
normal_amp.target, len(normal_amp.target) - 1
).adjoint
return normal_amp * (expand @ self.azm)
@property
def power_spectrum(self):
......
......@@ -59,17 +59,16 @@ def SimpleCorrelatedField(
else:
target.check_codomain(harmonic_partner)
harmonic_partner.check_codomain(target)
for kk in [offset_std, fluctuations, loglogavgslope]:
for kk in (fluctuations, loglogavgslope):
if len(kk) != 2:
raise TypeError
for kk in [flexibility, asperity]:
for kk in (offset_std, flexibility, asperity):
if not (kk is None or len(kk) == 2):
raise TypeError
if flexibility is None and asperity is not None:
raise ValueError
fluct = LognormalTransform(*fluctuations, prefix + 'fluctuations', 0)
avgsl = NormalTransform(*loglogavgslope, prefix + 'loglogavgslope', 0)
zm = LognormalTransform(*offset_std, prefix + 'zeromode', 0)
pspace = PowerSpace(harmonic_partner)
twolog = _TwoLogIntegrations(pspace)
......@@ -107,8 +106,11 @@ def SimpleCorrelatedField(
maskzm = np.ones(pspace.shape)
maskzm[0] = 0
maskzm = makeOp(makeField(pspace, maskzm))
insert = ValueInserter(pspace, (0,))
a = (maskzm @ ((ps_expander @ fluct)*a)) + insert(zm)
a = (maskzm @ ((ps_expander @ fluct)*a))
if offset_std is not None:
zm = LognormalTransform(*offset_std, prefix + 'zeromode', 0)
insert = ValueInserter(pspace, (0,))
a = a + insert(zm)
a = a.scale(target.total_volume)
ht = HarmonicTransformOperator(harmonic_partner, target)
......
......@@ -64,10 +64,10 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
c_first = np.ceil(start[:, i]+direction*dmin)
c_first = np.where(direction > 0., c_first, c_first-1.)
c_first = (c_first-start[:, i])/dirx
pos1 = np.asarray((start[:, i]+dmin*direction), dtype=np.int)
pos1 = np.asarray((start[:, i]+dmin*direction), dtype=int)
pos1 = np.sum(pos1*inc)
cdist = np.empty(0, dtype=np.float64)
add = np.empty(0, dtype=np.int)
add = np.empty(0, dtype=int)
for j in range(ndim):
if direction[j] != 0:
step = inc[j] if direction[j] > 0 else -inc[j]
......
......@@ -138,8 +138,8 @@ class LineSearch(metaclass=NiftyMeta):
max_zoom_iterations=100):
self.preferred_initial_step_size = preferred_initial_step_size
self.c1 = np.float(c1)
self.c2 = np.float(c2)
self.c1 = float(c1)
self.c2 = float(c2)
self.max_step_size = max_step_size
self.max_iterations = int(max_iterations)
self.max_zoom_iterations = int(max_zoom_iterations)
......
......@@ -97,8 +97,8 @@ class DOFDistributor(LinearOperator):
firstaxis = self._target.axes[self._space][0]
lastaxis = self._target.axes[self._space][-1]
arrshape = self._target.shape
presize = np.prod(arrshape[0:firstaxis], dtype=np.int)
postsize = np.prod(arrshape[lastaxis+1:], dtype=np.int)
presize = np.prod(arrshape[0:firstaxis], dtype=int)
postsize = np.prod(arrshape[lastaxis+1:], dtype=int)
self._hshape = (presize, self._domain[self._space].shape[0], postsize)
self._pshape = (presize, self._dofdex.size, postsize)
......
......@@ -66,7 +66,7 @@ class RegriddingOperator(LinearOperator):
self._frac = [None] * ndim
for d in range(ndim):
tmp = np.arange(new_shape[d])*(newdist[d]/dom.distances[d])
self._bindex[d] = np.minimum(dom.shape[d]-2, tmp.astype(np.int))
self._bindex[d] = np.minimum(dom.shape[d]-2, tmp.astype(int))
self._frac[d] = tmp-self._bindex[d]
def apply(self, x, mode):
......
......@@ -31,7 +31,7 @@ SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES]
@pmp('domain', SPACE_COMBINATIONS)
@pmp('attribute_desired_type',
[['domain', ift.DomainTuple], ['val', np.ndarray],
['shape', tuple], ['size', (np.int, np.int64)]])
['shape', tuple], ['size', (int, np.int64)]])
def test_return_types(domain, attribute_desired_type):
attribute = attribute_desired_type[0]
desired_type = attribute_desired_type[1]
......@@ -286,7 +286,7 @@ def test_stdfunc():
f = ift.Field.full(s, 27)
assert_equal(f.val, 27)
assert_equal(f.shape, (200,))
assert_equal(f.dtype, np.int)
assert_equal(f.dtype, int)
fx = ift.full(f.domain, 0)
assert_equal(f.dtype, fx.dtype)
assert_equal(f.shape, fx.shape)
......
......@@ -62,11 +62,12 @@ def testDistributor(dofdex, seed):
@pmp('total_N', [0, 1, 2])
@pmp('offset_std', [None, (1, 1)])
@pmp('asperity', [None, (1, 1)])
@pmp('flexibility', [None, (1, 1)])
@pmp('ind', [None, 1])
@pmp('matern', [True, False])
def test_init(total_N, asperity, flexibility, ind, matern):
def test_init(total_N, offset_std, asperity, flexibility, ind, matern):
if flexibility is None and asperity is not None:
pytest.skip()
cfg = 1, 1
......@@ -81,10 +82,39 @@ 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.set_amplitude_total_offset(0, cfg, dofdex=dofdex)
cfm.set_amplitude_total_offset(0, offset_std, dofdex=dofdex)
cfm.finalize(prior_info=0)
@pmp('sspace', spaces)
@pmp('asperity', [None, (1, 1)])
@pmp('flexibility', [None, (1, 1)])
@pmp('matern', [True, False])
def test_constant_zero_mode(sspace, asperity, flexibility, matern):
if flexibility is None and asperity is not None:
pytest.skip()
cfg = 1, 1
cfm = ift.CorrelatedFieldMaker('')
if matern:
cfm.add_fluctuations_matern(sspace, *(3 * [cfg]))
else:
cfm.add_fluctuations(sspace, *(4 * [cfg]))
cfm.set_amplitude_total_offset(0, None)
cf = cfm.finalize(prior_info=0)
r = ift.from_random(cf.domain).to_dict()
r_xi = np.copy(r["xi"].val)
r_xi[0] = 1.
r["xi"] = ift.Field(r["xi"].domain, r_xi)
r = ift.MultiField.from_dict(r)
cf_r = cf(r)
rtol = 1e-7
if isinstance(sspace, (ift.HPSpace, ift.GLSpace)):
rtol = 1e-2
assert_allclose(cf_r.s_integrate(), sspace.total_volume, rtol=rtol)
@pmp('sspace', spaces)
@pmp('N', [0, 2])
def testAmplitudesInvariants(sspace, N):
......@@ -149,7 +179,7 @@ def testAmplitudesInvariants(sspace, N):
@pmp('seed', [42, 31])
@pmp('domain', spaces)
@pmp('without', (('asperity', ), ('flexibility', ), ('flexibility', 'asperity')))
@pmp('without', (('offset_std', ), ('asperity', ), ('flexibility', ), ('flexibility', 'asperity')))
def test_complicated_vs_simple(seed, domain, without):
with ift.random.Context(seed):
offset_mean = _rand()
......
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