Commit 20577a56 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'switch_to_ducc' of gitlab.mpcdf.mpg.de:ift/nifty into switch_to_ducc

parents 98b214c2 fce04fb8
Pipeline #77095 passed with stages
in 25 minutes and 38 seconds
......@@ -9,6 +9,21 @@ now uses the DUCC package (<https://gitlab.mpcdf.mpg.de/mtr/ducc)>,
which is their successor.
Naming of operator tests
------------------------
The implementation tests for nonlinear operators are now available in
`ift.extra.check_operator()` and for linear operators
`ift.extra.check_linear_operator()`.
MetricGaussianKL interface
--------------------------
Users do not instanciate `MetricGaussianKL` by its constructor anymore. Rather
`MetricGaussianKL.make()` shall be used.
Changes since NIFTy 5
=====================
......@@ -70,6 +85,19 @@ print(met)
print(met.draw_sample())
```
New approach for sampling complex numbers
=========================================
When calling draw_sample_with_dtype with a complex dtype,
the variance is now used for the imaginary part and real part separately.
This is done in order to be consistent with the Hamiltonian.
Note that by this,
```
np.std(ift.from_random(domain, 'normal', dtype=np.complex128).val)
````
does not give 1, but sqrt(2) as a result.
MPI parallelisation over samples in MetricGaussianKL
----------------------------------------------------
......@@ -85,6 +113,7 @@ the generation of reproducible random numbers in the presence of MPI parallelism
and leads to cleaner code overall. Please see the documentation of
`nifty7.random` for details.
Interface Change for from_random and OuterProduct
-------------------------------------------------
......
......@@ -131,7 +131,7 @@ def main():
# Draw new samples to approximate the KL five times
for i in range(5):
# Draw new samples and minimize KL
KL = ift.MetricGaussianKL(mean, H, N_samples)
KL = ift.MetricGaussianKL.make(mean, H, N_samples)
KL, convergence = minimizer(KL)
mean = KL.position
......@@ -144,7 +144,7 @@ def main():
name=filename.format("loop_{:02d}".format(i)))
# Draw posterior samples
KL = ift.MetricGaussianKL(mean, H, N_samples)
KL = ift.MetricGaussianKL.make(mean, H, N_samples)
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(signal(sample + KL.position))
......
......@@ -152,10 +152,8 @@
"sigmas = [1.0, 0.5, 0.1]\n",
"\n",
"for i in range(3):\n",
" op = ift.library.correlated_fields._LognormalMomentMatching(mean=mean,\n",
" sig=sigmas[i],\n",
" key='foo',\n",
" N_copies=0)\n",
" op = ift.LognormalTransform(mean=mean, sigma=sigmas[i],\n",
" key='foo', N_copies=0)\n",
" op_samples = np.array(\n",
" [op(s).val for s in [ift.from_random(op.domain) for i in range(10000)]])\n",
"\n",
......
......@@ -131,7 +131,7 @@ def main():
for i in range(10):
# Draw new samples and minimize KL
KL = ift.MetricGaussianKL(mean, H, N_samples)
KL = ift.MetricGaussianKL.make(mean, H, N_samples)
KL, convergence = minimizer(KL)
mean = KL.position
......@@ -157,7 +157,7 @@ def main():
name=filename.format("loop_{:02d}".format(i)))
# Done, draw posterior samples
KL = ift.MetricGaussianKL(mean, H, N_samples)
KL = ift.MetricGaussianKL.make(mean, H, N_samples)
sc = ift.StatCalculator()
scA1 = ift.StatCalculator()
scA2 = ift.StatCalculator()
......
......@@ -34,6 +34,7 @@ from matplotlib.colors import LogNorm
import nifty7 as ift
def main():
dom = ift.UnstructuredDomain(1)
scale = 10
......@@ -90,7 +91,7 @@ def main():
plt.figure(figsize=[12, 8])
for ii in range(15):
if ii % 3 == 0:
mgkl = ift.MetricGaussianKL(pos, ham, 40)
mgkl = ift.MetricGaussianKL.make(pos, ham, 40)
plt.cla()
plt.imshow(z.T, origin='lower', norm=LogNorm(), vmin=1e-3,
......
......@@ -97,7 +97,7 @@ def main():
p_space = ift.UnstructuredDomain(N_params)
params = ift.full(p_space, 0.)
R = PolynomialResponse(p_space, x)
ift.extra.consistency_check(R)
ift.extra.check_linear_operator(R)
d_space = R.target
d = ift.makeField(d_space, y)
......
......@@ -52,6 +52,7 @@ from .operators.energy_operators import (
BernoulliEnergy, StandardHamiltonian, AveragedEnergy, QuadraticFormOperator,
Squared2NormOperator, StudentTEnergy, VariableCovarianceGaussianEnergy)
from .operators.convolution_operators import FuncConvolutionOperator
from .operators.normal_operators import NormalTransform, LognormalTransform
from .probing import probe_with_posterior_samples, probe_diagonal, \
StatCalculator, approximation2endo
......
This diff is collapsed.
......@@ -136,6 +136,8 @@ class Field(Operator):
The domain of the output random Field.
dtype : type
The datatype of the output random Field.
If the datatype is complex, each real and imaginary part
have variance 1
Returns
-------
......
......@@ -37,48 +37,11 @@ from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.linear_operator import LinearOperator
from ..operators.operator import Operator
from ..operators.simple_linear_operators import ducktape
from ..operators.normal_operators import NormalTransform, LognormalTransform
from ..probing import StatCalculator
from ..sugar import full, makeDomain, makeField, makeOp
def _reshaper(x, N):
x = np.asfarray(x)
if x.shape in [(), (1,)]:
return np.full(N, x) if N != 0 else x.reshape(())
elif x.shape == (N,):
return x
else:
raise TypeError("Shape of parameters cannot be interpreted")
def _lognormal_moments(mean, sig, N=0):
if N == 0:
mean, sig = np.asfarray(mean), np.asfarray(sig)
else:
mean, sig = (_reshaper(param, N) for param in (mean, sig))
if not np.all(mean > 0):
raise ValueError("mean must be greater 0; got {!r}".format(mean))
if not np.all(sig > 0):
raise ValueError("sig must be greater 0; got {!r}".format(sig))
logsig = np.sqrt(np.log1p((sig/mean)**2))
logmean = np.log(mean) - logsig**2/2
return logmean, logsig
def _normal(mean, sig, key, N=0):
if N == 0:
domain = DomainTuple.scalar_domain()
mean, sig = np.asfarray(mean), np.asfarray(sig)
return Adder(makeField(domain, mean)) @ (
sig * ducktape(domain, None, key))
domain = UnstructuredDomain(N)
mean, sig = (_reshaper(param, N) for param in (mean, sig))
return Adder(makeField(domain, mean)) @ (DiagonalOperator(
makeField(domain, sig)) @ ducktape(domain, None, key))
def _log_k_lengths(pspace):
"""Log(k_lengths) without zeromode"""
return np.log(pspace.k_lengths[1:])
......@@ -120,29 +83,6 @@ def _total_fluctuation_realized(samples):
return np.sqrt(res if np.isscalar(res) else res.val)
class _LognormalMomentMatching(Operator):
def __init__(self, mean, sig, key, N_copies):
key = str(key)
logmean, logsig = _lognormal_moments(mean, sig, N_copies)
self._mean = mean
self._sig = sig
op = _normal(logmean, logsig, key, N_copies).ptw("exp")
self._domain, self._target = op.domain, op.target
self.apply = op.apply
self._repr_str = f"_LognormalMomentMatching: " + op.__repr__()
@property
def mean(self):
return self._mean
@property
def std(self):
return self._sig
def __repr__(self):
return self._repr_str
class _SlopeRemover(EndomorphicOperator):
def __init__(self, domain, space=0):
self._domain = makeDomain(domain)
......@@ -441,10 +381,8 @@ class CorrelatedFieldMaker:
elif len(dofdex) != total_N:
raise ValueError("length of dofdex needs to match total_N")
N = max(dofdex) + 1 if total_N > 0 else 0
zm = _LognormalMomentMatching(offset_std_mean,
offset_std_std,
prefix + 'zeromode',
N)
zm = LognormalTransform(offset_std_mean, offset_std_std,
prefix + 'zeromode', N)
if total_N > 0:
zm = _Distributor(dofdex, zm.target, UnstructuredDomain(total_N)) @ zm
return CorrelatedFieldMaker(offset_mean, zm, prefix, total_N)
......@@ -532,17 +470,15 @@ class CorrelatedFieldMaker:
prefix = str(prefix)
# assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)
fluct = _LognormalMomentMatching(fluctuations_mean,
fluctuations_stddev,
self._prefix + prefix + 'fluctuations',
N)
flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev,
self._prefix + prefix + 'flexibility',
N)
asp = _LognormalMomentMatching(asperity_mean, asperity_stddev,
self._prefix + prefix + 'asperity', N)
avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
self._prefix + prefix + 'loglogavgslope', N)
fluct = LognormalTransform(fluctuations_mean, fluctuations_stddev,
self._prefix + prefix + 'fluctuations', N)
flex = LognormalTransform(flexibility_mean, flexibility_stddev,
self._prefix + prefix + 'flexibility', N)
asp = LognormalTransform(asperity_mean, asperity_stddev,
self._prefix + prefix + 'asperity', N)
avgsl = NormalTransform(loglogavgslope_mean, loglogavgslope_stddev,
self._prefix + prefix + 'loglogavgslope', N)
amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
self._azm, target_subdomain[-1].total_volume,
self._prefix + prefix + 'spectrum', dofdex)
......
This diff is collapsed.
......@@ -112,6 +112,8 @@ class MultiField(Operator):
The domain of the output random Field.
dtype : type
The datatype of the output random Field.
If the datatype is complex, each real an imaginary part have
variance 1.
Returns
-------
......@@ -248,6 +250,10 @@ class MultiField(Operator):
return MultiField(subset,
tuple(self[key] for key in subset.keys()))
def extract_by_keys(self, keys):
dom = MultiDomain.make({kk: vv for kk, vv in self.domain.items() if kk in keys})
return self.extract(dom)
def extract_part(self, subset):
if subset is self._domain:
return self
......
......@@ -17,6 +17,7 @@
from ..multi_domain import MultiDomain
from ..multi_field import MultiField
from ..utilities import indent
from .endomorphic_operator import EndomorphicOperator
from .linear_operator import LinearOperator
......@@ -46,7 +47,6 @@ class BlockDiagonalOperator(EndomorphicOperator):
else:
raise TypeError("LinearOperator expected")
def apply(self, x, mode):
self._check_input(x, mode)
val = tuple(op.apply(v, mode=mode) if op is not None else v
......@@ -80,3 +80,7 @@ class BlockDiagonalOperator(EndomorphicOperator):
res = {key: SumOperator.make([v1, v2], [selfneg, opneg])
for key, v1, v2 in zip(self._domain.keys(), self._ops, op._ops)}
return BlockDiagonalOperator(self._domain, res)
def __repr__(self):
s = "\n".join(f'{kk}: {self._ops[ii]}' for ii, kk in enumerate(self.domain.keys()))
return 'BlockDiagonalOperator:\n' + indent(s)
......@@ -58,7 +58,14 @@ class ChainOperator(LinearOperator):
fct = 1.
opsnew = []
lastdom = ops[-1].domain
dtype = None
for op in ops:
from .sampling_enabler import SamplingDtypeSetter
if isinstance(op, SamplingDtypeSetter) and isinstance(op._op, ScalingOperator):
if dtype is not None:
raise NotImplementedError
dtype = op._dtype
op = op._op
if (isinstance(op, ScalingOperator) and op._factor.imag == 0):
fct *= op._factor.real
else:
......@@ -72,7 +79,10 @@ class ChainOperator(LinearOperator):
break
if fct != 1 or len(opsnew) == 0:
# have to add the scaling operator at the end
opsnew.append(ScalingOperator(lastdom, fct))
op = ScalingOperator(lastdom, fct)
if dtype is not None:
op = SamplingDtypeSetter(op, dtype)
opsnew.append(op)
ops = opsnew
# combine DiagonalOperators where possible
opsnew = []
......@@ -142,7 +152,6 @@ class ChainOperator(LinearOperator):
from ..multi_domain import MultiDomain
if not isinstance(self._domain, MultiDomain):
return None, self
newop = None
for op in reversed(self._ops):
c_inp, t_op = op.simplify_for_constant_input(c_inp)
......
......@@ -48,6 +48,10 @@ def _check_sampling_dtype(domain, dtypes):
raise TypeError
def _iscomplex(dtype):
return np.issubdtype(dtype, np.complexfloating)
def _field_to_dtype(field):
if isinstance(field, Field):
dt = field.dtype
......@@ -127,10 +131,10 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
The covariance is assumed to be diagonal.
.. math ::
E(s,D) = - \\log G(s, D) = 0.5 (s)^\\dagger D^{-1} (s) + 0.5 tr log(D),
E(s,D) = - \\log G(s, C) = 0.5 (s)^\\dagger C (s) - 0.5 tr log(C),
an information energy for a Gaussian distribution with residual s and
diagonal covariance D.
inverse diagonal covariance C.
The domain of this energy will be a MultiDomain with two keys,
the target will be the scalar domain.
......@@ -139,10 +143,10 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
domain : Domain, DomainTuple, tuple of Domain
domain of the residual and domain of the covariance diagonal.
residual : key
residual_key : key
Residual key of the Gaussian.
inverse_covariance : key
inverse_covariance_key : key
Inverse covariance diagonal key of the Gaussian.
sampling_dtype : np.dtype
......@@ -156,7 +160,7 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
self._domain = MultiDomain.make({self._kr: dom, self._ki: dom})
self._dt = {self._kr: sampling_dtype, self._ki: np.float64}
_check_sampling_dtype(self._domain, self._dt)
self._cplx = np.issubdtype(sampling_dtype, np.complexfloating)
self._cplx = _iscomplex(sampling_dtype)
def apply(self, x):
self._check_input(x)
......@@ -171,6 +175,46 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
met = MultiField.from_dict({self._kr: i.val, self._ki: met**(-2)})
return res.add_metric(SamplingDtypeSetter(makeOp(met), self._dt))
def _simplify_for_constant_input_nontrivial(self, c_inp):
from .simplify_for_const import ConstantEnergyOperator
assert len(c_inp.keys()) == 1
key = c_inp.keys()[0]
assert key in self._domain.keys()
cst = c_inp[key]
if key == self._kr:
res = _SpecialGammaEnergy(cst).ducktape(self._ki)
else:
dt = self._dt[self._kr]
res = GaussianEnergy(inverse_covariance=makeOp(cst),
sampling_dtype=dt).ducktape(self._kr)
trlog = cst.log().sum().val_rw()
if not _iscomplex(dt):
trlog /= 2
res = res + ConstantEnergyOperator(-trlog)
res = res + ConstantEnergyOperator(0.)
assert res.target is self.target
return None, res
class _SpecialGammaEnergy(EnergyOperator):
def __init__(self, residual):
self._domain = DomainTuple.make(residual.domain)
self._resi = residual
self._cplx = _iscomplex(self._resi.dtype)
self._scale = ScalingOperator(self._domain, 1 if self._cplx else .5)
def apply(self, x):
self._check_input(x)
r = self._resi
if self._cplx:
res = 0.5*(r*x.real).vdot(r).real - x.ptw("log").sum()
else:
res = 0.5*((r*x).vdot(r) - x.ptw("log").sum())
if not x.want_metric:
return res
met = makeOp((self._scale(x.val))**(-2))
return res.add_metric(SamplingDtypeSetter(met, self._resi.dtype))
class GaussianEnergy(EnergyOperator):
"""Computes a negative-log Gaussian.
......@@ -192,6 +236,13 @@ class GaussianEnergy(EnergyOperator):
domain : Domain, DomainTuple, tuple of Domain or MultiDomain
Operator domain. By default it is inferred from `mean` or
`covariance` if specified
sampling_dtype : type
Here one can specify whether the distribution is a complex Gaussian or
not. Note that for a complex Gaussian the inverse_covariance is
.. math ::
(<ff^dagger>)^{-1}_P(f)/2,
where the additional factor of 2 is necessary because the
domain of s has double as many dimensions as in the real case.
Note
----
......@@ -225,14 +276,13 @@ class GaussianEnergy(EnergyOperator):
if sampling_dtype != _field_to_dtype(self._mean):
raise ValueError("Sampling dtype and mean not compatible")
self._icov = inverse_covariance
if inverse_covariance is None:
self._op = Squared2NormOperator(self._domain).scale(0.5)
self._met = ScalingOperator(self._domain, 1)
self._trivial_invcov = True
else:
self._op = QuadraticFormOperator(inverse_covariance)
self._met = inverse_covariance
self._trivial_invcov = False
if sampling_dtype is not None:
self._met = SamplingDtypeSetter(self._met, sampling_dtype)
......@@ -440,11 +490,9 @@ class StandardHamiltonian(EnergyOperator):
`<https://arxiv.org/abs/1812.04403>`_
"""
def __init__(self, lh, ic_samp=None, _c_inp=None, prior_dtype=np.float64):
def __init__(self, lh, ic_samp=None, prior_dtype=np.float64):
self._lh = lh
self._prior = GaussianEnergy(domain=lh.domain, sampling_dtype=prior_dtype)
if _c_inp is not None:
_, self._prior = self._prior.simplify_for_constant_input(_c_inp)
self._ic_samp = ic_samp
self._domain = lh.domain
......@@ -462,7 +510,7 @@ class StandardHamiltonian(EnergyOperator):
def _simplify_for_constant_input_nontrivial(self, c_inp):
out, lh1 = self._lh.simplify_for_constant_input(c_inp)
return out, StandardHamiltonian(lh1, self._ic_samp, _c_inp=c_inp)
return out, StandardHamiltonian(lh1, self._ic_samp)
class AveragedEnergy(EnergyOperator):
......
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.unstructured_domain import UnstructuredDomain
from ..operators.operator import Operator
from ..operators.adder import Adder
from ..operators.simple_linear_operators import ducktape
from ..operators.diagonal_operator import DiagonalOperator
from ..sugar import makeField
from ..utilities import value_reshaper, lognormal_moments
def NormalTransform(mean, sigma, key, N_copies=0):
"""Opchain that transforms standard normally distributed values to
normally distributed values with given mean an standard deviation.
Parameters:
-----------
mean : float
Mean of the field
sigma : float
Standard deviation of the field
key : string
Name of the operators domain (Multidomain)
N_copies : integer
If == 0, target will be a scalar field.
If >= 1, target will be an
:class:`~nifty7.unstructured_domain.UnstructuredDomain`.
"""
if N_copies == 0:
domain = DomainTuple.scalar_domain()
mean, sigma = np.asfarray(mean), np.asfarray(sigma)
mean_adder = Adder(makeField(domain, mean))
return mean_adder @ (sigma * ducktape(domain, None, key))
domain = UnstructuredDomain(N_copies)
mean, sigma = (value_reshaper(param, N_copies) for param in (mean, sigma))
mean_adder = Adder(makeField(domain, mean))
sigma_op = DiagonalOperator(makeField(domain, sigma))
return mean_adder @ sigma_op @ ducktape(domain, None, key)
def LognormalTransform(mean, sigma, key, N_copies):
"""Opchain that transforms standard normally distributed values to
log-normally distributed values with given mean an standard deviation.
Parameters:
-----------
mean : float
Mean of the field
sigma : float
Standard deviation of the field
key : string
Name of the domain
N_copies : integer
If == 0, target will be a scalar field.
If >= 1, target will be an
:class:`~nifty7.unstructured_domain.UnstructuredDomain`.
"""
logmean, logsigma = lognormal_moments(mean, sigma, N_copies)
return NormalTransform(logmean, logsigma, key, N_copies).ptw("exp")
......@@ -273,33 +273,49 @@ class Operator(metaclass=NiftyMeta):
def simplify_for_constant_input(self, c_inp):
from .energy_operators import EnergyOperator
from .simplify_for_const import ConstantEnergyOperator, ConstantOperator
if c_inp is None:
from ..multi_field import MultiField
from ..domain_tuple import DomainTuple
from ..sugar import makeDomain
if c_inp is None or (isinstance(c_inp, MultiField) and len(c_inp.keys()) == 0):
return None, self
dom = c_inp.domain
if isinstance(dom, MultiDomain) and len(dom) == 0:
return None, self
# Convention: If c_inp is MultiField, it needs to be defined on a
# subdomain of self._domain
if isinstance(self.domain, MultiDomain):
assert isinstance(c_inp.domain, MultiDomain)
if set(c_inp.keys()) > set(self.domain.keys()):
assert isinstance(dom, MultiDomain)
if not set(c_inp.keys()) <= set(self.domain.keys()):
raise ValueError
if c_inp.domain is self.domain:
if dom is self.domain:
if isinstance(self, DomainTuple):
raise RuntimeError
if isinstance(self, EnergyOperator):
op = ConstantEnergyOperator(self.domain, self(c_inp))
op = ConstantEnergyOperator(self(c_inp))
else:
op = ConstantOperator(self.domain, self(c_inp))
op = ConstantOperator(self.domain, self(c_inp))
return op(c_inp), op
if not isinstance(c_inp.domain, MultiDomain):
op = ConstantOperator(self(c_inp))
return None, op
if not isinstance(dom, MultiDomain):
raise RuntimeError
return self._simplify_for_constant_input_nontrivial(c_inp)
c_out, op = self._simplify_for_constant_input_nontrivial(c_inp)
vardom = makeDomain({kk: vv for kk, vv in self.domain.items()
if kk not in c_inp.keys()})
assert op.domain is vardom
assert op.target is self.target
assert isinstance(op, Operator)
if c_out is not None:
assert isinstance(c_out, MultiField)
assert len(set(c_out.keys()) & self.domain.keys()) == 0
assert set(c_out.keys()) <= set(c_inp.keys())
return c_out, op
def _simplify_for_constant_input_nontrivial(self, c_inp):