Skip to content

GitLab

Commit 773e7b8c authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge remote-tracking branch 'origin/NIFTy_6' into NIFTy_7

parents 03765dbf ad0434e3
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
View file @ 773e7b8c
......@@ -13,7 +13,7 @@ stages:
build_docker_from_scratch:
only:
- schedules
image: docker:stable
image: docker:19.03.8
stage: build_docker
before_script:
- ls
......@@ -25,7 +25,7 @@ build_docker_from_scratch:
build_docker_from_cache:
except:
- schedules
image: docker:stable
image: docker:19.03.8
stage: build_docker
before_script:
- ls
......
ChangeLog
View file @ 773e7b8c
......@@ -16,6 +16,8 @@ In addition to the below changes, the following operators were introduced:
* PartialConjugate: Conjugates parts of a multi-field
* SliceOperator: Geometry preserving mask operator
* SplitOperator: Splits a single field into a multi-field
* MatrixProductOperator: Applies matrices (scipy.sparse, numpy) to fields
* IntegrationOperator: Integrates over subspaces of fields
FFT convention adjusted
=======================
......
src/library/correlated_fields.py
View file @ 773e7b8c
......@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
# Authors: Philipp Frank, Philipp Arras, Philipp Haim
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -25,7 +25,6 @@ from ..domain_tuple import DomainTuple
from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..linearization import Linearization
from ..logger import logger
from ..multi_field import MultiField
from ..operators.adder import Adder
......@@ -39,6 +38,7 @@ from ..operators.operator import Operator
from ..operators.simple_linear_operators import ducktape
from ..probing import StatCalculator
from ..sugar import full, makeDomain, makeField, makeOp
from .. import utilities
def _reshaper(x, N):
......@@ -244,10 +244,9 @@ class _SpecialSum(EndomorphicOperator):
class _Distributor(LinearOperator):
def __init__(self, dofdex, domain, target):
self._dofdex = dofdex
self._target = makeDomain(target)
self._domain = makeDomain(domain)
self._dofdex = np.array(dofdex)
self._target = DomainTuple.make(target)
self._domain = DomainTuple.make(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
......@@ -256,8 +255,8 @@ class _Distributor(LinearOperator):
if mode == self.TIMES:
res = x[self._dofdex]
else:
res = np.empty(self._tgt(mode).shape)
res[self._dofdex] = x
res = np.zeros(self._tgt(mode).shape, dtype=x.dtype)
res = utilities.special_add_at(res, 0, self._dofdex, x)
return makeField(self._tgt(mode), res)
......@@ -354,11 +353,46 @@ class _Amplitude(Operator):
class CorrelatedFieldMaker:
"""Constrution helper for hirarchical correlated field models.
The correlated field models are parametrized by creating
power spectrum operators ("amplitudes")
acting on their target subdomains
via calls to :func:`add_fluctuations`.
During creation of the :class:`CorrelatedFieldMaker` via
:func:`make`, a global offset from zero can be added to the
field to be created and an operator applying gaussian fluctuations
around this offset needs to be parametrized.
The resulting correlated field model operator has a
:class:`~nifty6.multi_domain.MultiDomain` as its domain and
expects its input values to be univariately gaussian.
The target of the constructed operator will be a
:class:`~nifty6.domain_tuple.DomainTuple`
containing the `target_subdomains` of the added fluctuations in the
order of the `add_fluctuations` calls.
Creation of the model operator is finished by calling the method
:func:`finalize`, which returns the configured operator.
An operator representing an array of correlated field models
can be constructed by setting the `total_N` parameter of
:func:`make`. It will have an
:class:`~nifty.domains.unstructucture_domain.UnstructureDomain`
of shape `(total_N,)` prepended to its target domain and represent
`total_N` correlated fields simulataneously.
The degree of information sharing between the correlated field
models can be configured via the `dofdex` parameters
of :func:`make` and :func:`add_fluctuations`.
See the methods :func:`make`, :func:`add_fluctuations`
and :func:`finalize` for further usage information."""
def __init__(self, offset_mean, offset_fluctuations_op, prefix, total_N):
if not isinstance(offset_fluctuations_op, Operator):
raise TypeError("offset_fluctuations_op needs to be an operator")
self._a = []
self._position_spaces = []
self._target_subdomains = []
self._offset_mean = offset_mean
self._azm = offset_fluctuations_op
......@@ -366,8 +400,7 @@ class CorrelatedFieldMaker:
self._total_N = total_N
@staticmethod
def make(offset_mean, offset_std_mean, offset_std_std, prefix,
total_N=0,
def make(offset_mean, offset_std_mean, offset_std_std, prefix, total_N=0,
dofdex=None):
"""Returns a CorrelatedFieldMaker object.
......@@ -376,15 +409,25 @@ class CorrelatedFieldMaker:
offset_mean : float
Mean offset from zero of the correlated field to be made.
offset_std_mean : float
Mean standard deviation of the offset value.
Mean standard deviation of the offset.
offset_std_std : float
Standard deviation of the stddev of the offset value.
Standard deviation of the stddev of the offset.
prefix : string
Prefix to the names of the domains of the cf operator to be made.
total_N : integer
?
dofdex : np.array
?
This determines the names of the operator domain.
total_N : integer, optional
Number of field models to create.
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
instanciated; 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 dofdex is None:
dofdex = np.full(total_N, 0)
......@@ -400,7 +443,7 @@ class CorrelatedFieldMaker:
return CorrelatedFieldMaker(offset_mean, zm, prefix, total_N)
def add_fluctuations(self,
position_space,
target_subdomain,
fluctuations_mean,
fluctuations_stddev,
flexibility_mean,
......@@ -413,11 +456,60 @@ class CorrelatedFieldMaker:
index=None,
dofdex=None,
harmonic_partner=None):
"""Function to add correlation structures to the field to be made.
Correlations are described by their power spectrum and the subdomain
on which they apply.
The parameters `fluctuations`, `flexibility`, `asperity` and
`loglogavgslope` configure the power spectrum model ("amplitude")
used on the target field subdomain `target_subdomain`.
It is assembled as the sum of a power law component
(linear slope in log-log power-frequency-space),
a smooth varying component (integrated wiener process) and
a ragged componenent (unintegrated wiener process).
Multiple calls to `add_fluctuations` are possible, in which case
the constructed field will have the outer product of the individual
power spectra as its global power spectrum.
Parameters
----------
target_subdomain : :class:`~nifty6.domain.Domain`, \
:class:`~nifty6.domain_tuple.DomainTuple`
Target subdomain on which the correlation structure defined
in this call should hold.
fluctuations_{mean,stddev} : float
Total spectral energy -> Amplitude of the fluctuations
flexibility_{mean,stddev} : float
Smooth variation speed of the power spectrum
asperity_{mean,stddev} : float
Strength of unsmoothed power spectrum variations
Used to accomodate single frequency peaks
loglogavgslope_{mean,stddev} : float
Power law component exponent
prefix : string
prefix of the power spectrum parameter domain names
index : int
Position target_subdomain in the final total domain of the
correlated field operator.
dofdex : np.array, optional
An integer array specifying the power spectrum 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 power spectrum models are
instanciated; the first one is used for the first and second
field model and the second one is used for the third field model.
*If not given*, use the same power spectrum model for all
constructed field models.
harmonic_partner : :class:`~nifty6.domain.Domain`, \
:class:`~nifty6.domain_tuple.DomainTuple`
In which harmonic space to define the power spectrum
"""
if harmonic_partner is None:
harmonic_partner = position_space.get_default_codomain()
harmonic_partner = target_subdomain.get_default_codomain()
else:
position_space.check_codomain(harmonic_partner)
harmonic_partner.check_codomain(position_space)
target_subdomain.check_codomain(harmonic_partner)
harmonic_partner.check_codomain(target_subdomain)
if dofdex is None:
dofdex = np.full(self._total_N, 0)
......@@ -426,12 +518,12 @@ class CorrelatedFieldMaker:
if self._total_N > 0:
N = max(dofdex) + 1
position_space = makeDomain((UnstructuredDomain(N), position_space))
target_subdomain = makeDomain((UnstructuredDomain(N), target_subdomain))
else:
N = 0
position_space = makeDomain(position_space)
target_subdomain = makeDomain(target_subdomain)
prefix = str(prefix)
# assert isinstance(position_space[space], (RGSpace, HPSpace, GLSpace)
# assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)
fluct = _LognormalMomentMatching(fluctuations_mean,
fluctuations_stddev,
......@@ -445,17 +537,26 @@ class CorrelatedFieldMaker:
avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
self._prefix + prefix + 'loglogavgslope', N)
amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
self._azm, position_space[-1].total_volume,
self._azm, target_subdomain[-1].total_volume,
self._prefix + prefix + 'spectrum', dofdex)
if index is not None:
self._a.insert(index, amp)
self._position_spaces.insert(index, position_space)
self._target_subdomains.insert(index, target_subdomain)
else:
self._a.append(amp)
self._position_spaces.append(position_space)
self._target_subdomains.append(target_subdomain)
def finalize(self, prior_info=100):
"""Finishes model construction process and returns the constructed
operator.
def _finalize_from_op(self):
Parameters
----------
prior_info : integer
How many prior samples to draw for property verification statistics
If zero, skips calculating and displaying statistics.
"""
n_amplitudes = len(self._a)
if self._total_N > 0:
hspace = makeDomain(
......@@ -473,11 +574,11 @@ class CorrelatedFieldMaker:
azm = expander @ self._azm
ht = HarmonicTransformOperator(hspace,
self._position_spaces[0][amp_space],
self._target_subdomains[0][amp_space],
space=spaces[0])
for i in range(1, n_amplitudes):
ht = HarmonicTransformOperator(ht.target,
self._position_spaces[i][amp_space],
self._target_subdomains[i][amp_space],
space=spaces[i]) @ ht
a = []
for ii in range(n_amplitudes):
......@@ -486,13 +587,8 @@ class CorrelatedFieldMaker:
pd = PowerDistributor(co.target, pp, amp_space)
a.append(co.adjoint @ pd @ self._a[ii])
corr = reduce(mul, a)
return ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
op = ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
def finalize(self, prior_info=100):
"""
offset vs zeromode: volume factor
"""
op = self._finalize_from_op()
if self._offset_mean is not None:
offset = self._offset_mean
# Deviations from this offset must not be considered here as they
......@@ -546,6 +642,7 @@ class CorrelatedFieldMaker:
@property
def normalized_amplitudes(self):
"""Returns the power spectrum operators used in the model"""
return self._a
@property
......
src/minimization/metric_gaussian_kl.py
View file @ 773e7b8c
......@@ -114,6 +114,11 @@ class MetricGaussianKL(Energy):
If not None, samples will be distributed as evenly as possible
across this communicator. If `mirror_samples` is set, then a sample and
its mirror image will always reside on the same task.
nanisinf : bool
If true, nan energies which can happen due to overflows in the forward
model are interpreted as inf. Thereby, the code does not crash on
these occaisions but rather the minimizer is told that the position it
has tried is not sensible.
_local_samples : None
Only a parameter for internal uses. Typically not to be set by users.
......@@ -131,7 +136,8 @@ class MetricGaussianKL(Energy):
def __init__(self, mean, hamiltonian, n_samples, constants=[],
point_estimates=[], mirror_samples=False,
napprox=0, comm=None, _local_samples=None):
napprox=0, comm=None, _local_samples=None,
nanisinf=False):
super(MetricGaussianKL, self).__init__(mean)
if not isinstance(hamiltonian, StandardHamiltonian):
......@@ -142,6 +148,7 @@ class MetricGaussianKL(Energy):
raise TypeError
self._constants = tuple(constants)
self._point_estimates = tuple(point_estimates)
self._mitigate_nans = nanisinf
if not isinstance(mirror_samples, bool):
raise TypeError
......@@ -195,6 +202,8 @@ class MetricGaussianKL(Energy):
v += tmp.val.val
g = g + tmp.gradient
self._val = _np_allreduce_sum(self._comm, v)[()] / self._n_eff_samples
if np.isnan(self._val) and self._mitigate_nans:
self._val = np.inf
self._grad = _allreduce_sum_field(self._comm, g) / self._n_eff_samples
self._metric = None
......@@ -202,7 +211,7 @@ class MetricGaussianKL(Energy):
return MetricGaussianKL(
position, self._hamiltonian, self._n_samples, self._constants,
self._point_estimates, self._mirror_samples, comm=self._comm,
_local_samples=self._local_samples)
_local_samples=self._local_samples, nanisinf=self._mitigate_nans)
@property
def value(self):
......
src/multi_field.py
View file @ 773e7b8c
......@@ -83,9 +83,9 @@ class MultiField(Operator):
def domain(self):
return self._domain
# @property
# def dtype(self):
# return {key: val.dtype for key, val in self._val.items()}
@property
def dtype(self):
return {key: val.dtype for key, val in self.items()}
def _transform(self, op):
return MultiField(self._domain, tuple(op(v) for v in self._val))
......
src/operators/outer_product_operator.py
View file @ 773e7b8c
......@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -35,7 +35,7 @@ class OuterProduct(LinearOperator):
self._domain = DomainTuple.make(domain)
self._field = field
self._target = DomainTuple.make(
tuple(sub_d for sub_d in field.domain._dom + domain._dom))
tuple(sub_d for sub_d in field.domain._dom + self._domain._dom))
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
......
src/operators/sampling_enabler.py
View file @ 773e7b8c
......@@ -124,6 +124,10 @@ class SamplingDtypeSetter(EndomorphicOperator):
need to conincide the with keys of the `MultiDomain`.
"""
def __init__(self, op, dtype):
if isinstance(op, SamplingDtypeSetter):
if op._dtype != dtype:
raise ValueError('Dtype for sampling already set to another dtype.')
op = op._op
if not isinstance(op, EndomorphicOperator):
raise TypeError
if not hasattr(op, 'draw_sample_with_dtype'):
......
src/plot.py
View file @ 773e7b8c
......@@ -405,7 +405,7 @@ def _plot2D(f, ax, **kwargs):
ax.set_ylabel(kwargs.pop("ylabel", ""))
dom = dom[x_space]
if not have_rgb:
cmap = kwargs.pop("colormap", plt.rcParams['image.cmap'])
cmap = kwargs.pop("cmap", plt.rcParams['image.cmap'])
if isinstance(dom, RGSpace):
nx, ny = dom.shape
......@@ -417,7 +417,7 @@ def _plot2D(f, ax, **kwargs):
else:
im = ax.imshow(
f.val.T, extent=[0, nx*dx, 0, ny*dy],
vmin=kwargs.get("zmin"), vmax=kwargs.get("zmax"),
vmin=kwargs.get("vmin"), vmax=kwargs.get("vmax"),
cmap=cmap, origin="lower", **norm, **aspect)
plt.colorbar(im)
_limit_xy(**kwargs)
......@@ -453,7 +453,7 @@ def _plot2D(f, ax, **kwargs):
if have_rgb:
plt.imshow(res, origin="lower")
else:
plt.imshow(res, vmin=kwargs.get("zmin"), vmax=kwargs.get("zmax"),
plt.imshow(res, vmin=kwargs.get("vmin"), vmax=kwargs.get("vmax"),
norm=norm.get('norm'), cmap=cmap, origin="lower")
plt.colorbar(orientation="horizontal")
return
......@@ -518,7 +518,7 @@ class Plot(object):
Label for the y axis.
[xyz]min, [xyz]max: float
Limits for the values to plot.
colormap: string
cmap: string
Color map to use for the plot (if it is a 2D plot).
linewidth: float or list of floats
Line width.
......
src/pointwise.py
View file @ 773e7b8c
......@@ -25,9 +25,13 @@ def _sqrt_helper(v):
def _sinc_helper(v):
tmp = np.sinc(v)
tmp2 = (np.cos(np.pi*v)-tmp)/v
return (tmp, np.where(v==0., 0, tmp2))
fv = np.sinc(v)
df = np.empty(v.shape, dtype=v.dtype)
sel = v != 0.
v = v[sel]
df[sel] = (np.cos(np.pi*v)-fv[sel])/v
df[~sel] = 0
return (fv, df)
def _expm1_helper(v):
......@@ -54,13 +58,13 @@ def _reciprocal_helper(v):
def _abs_helper(v):
if np.issubdtype(v.dtype, np.complexfloating):
raise TypeError("Argument must not be complex")
return (np.abs(v), np.where(v==0, np.nan, np.sign(v)))
return (np.abs(v), np.where(v == 0, np.nan, np.sign(v)))
def _sign_helper(v):
if np.issubdtype(v.dtype, np.complexfloating):
raise TypeError("Argument must not be complex")
return (np.sign(v), np.where(v==0, np.nan, 0))
return (np.sign(v), np.where(v == 0, np.nan, 0))
def _power_helper(v, expo):
......@@ -73,21 +77,21 @@ def _clip_helper(v, a_min, a_max):
tmp = np.clip(v, a_min, a_max)
tmp2 = np.ones(v.shape)
if a_min is not None:
tmp2 = np.where(tmp==a_min, 0., tmp2)
tmp2 = np.where(tmp == a_min, 0., tmp2)
if a_max is not None:
tmp2 = np.where(tmp==a_max, 0., tmp2)
tmp2 = np.where(tmp == a_max, 0., tmp2)
return (tmp, tmp2)
ptw_dict = {
"sqrt": (np.sqrt, _sqrt_helper),
"sin" : (np.sin, lambda v: (np.sin(v), np.cos(v))),
"cos" : (np.cos, lambda v: (np.cos(v), -np.sin(v))),
"tan" : (np.tan, lambda v: (np.tan(v), 1./np.cos(v)**2)),
"sin": (np.sin, lambda v: (np.sin(v), np.cos(v))),
"cos": (np.cos, lambda v: (np.cos(v), -np.sin(v))),
"tan": (np.tan, lambda v: (np.tan(v), 1./np.cos(v)**2)),
"sinc": (np.sinc, _sinc_helper),
"exp" : (np.exp, lambda v: (2*(np.exp(v),))),
"expm1" : (np.expm1, _expm1_helper),
"log" : (np.log, lambda v: (np.log(v), 1./v)),
"exp": (np.exp, lambda v: (2*(np.exp(v),))),
"expm1": (np.expm1, _expm1_helper),
"log": (np.log, lambda v: (np.log(v), 1./v)),
"log10": (np.log10, lambda v: (np.log10(v), (1./np.log(10.))/v)),
"log1p": (np.log1p, lambda v: (np.log1p(v), 1./(1.+v))),
"sinh": (np.sinh, lambda v: (np.sinh(v), np.cosh(v))),
......
src/sugar.py
View file @ 773e7b8c
......@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -20,21 +20,20 @@ from time import time
import numpy as np
from .logger import logger
from . import utilities
from . import pointwise, utilities
from .domain_tuple import DomainTuple
from .domains.power_space import PowerSpace
from .field import Field
from .logger import logger
from .multi_domain import MultiDomain
from .multi_field import MultiField
from .operators.block_diagonal_operator import BlockDiagonalOperator
from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import PowerDistributor
from .operators.operator import Operator
from .operators.sampling_enabler import SamplingDtypeSetter
from .operators.scaling_operator import ScalingOperator
from .plot import Plot
from . import pointwise
__all__ = ['PS_field', 'power_analyze', 'create_power_operator',
'create_harmonic_smoothing_operator', 'from_random',
......@@ -501,17 +500,26 @@ def calculate_position(operator, output):
if output.domain != operator.target:
raise TypeError
if isinstance(output, MultiField):
cov = 1e-3*max([vv.max() for vv in output.val.values()])**2
cov = 1e-3*max([np.max(np.abs(vv)) for vv in output.val.values()])**2
invcov = ScalingOperator(output.domain, cov).inverse
dtype = list(set([ff.dtype for ff in output.values()]))
if len(dtype) != 1:
raise ValueError('Only MultiFields with one dtype supported.')
dtype = dtype[0]
else:
cov = 1e-3*output.val.max()**2
cov = 1e-3*np.max(np.abs(output.val))**2
dtype = output.dtype
invcov = ScalingOperator(output.domain, cov).inverse
d = output + invcov.draw_sample_with_dtype(dtype=output.dtype, from_inverse=True)
invcov = SamplingDtypeSetter(invcov, output.dtype)
invcov = SamplingDtypeSetter(invcov, output.dtype)
d = output + invcov.draw_sample(from_inverse=True)
lh = GaussianEnergy(d, invcov) @ operator
H = StandardHamiltonian(
lh, ic_samp=GradientNormController(iteration_limit=200))
pos = 0.1 * from_random(operator.domain, 'normal')
minimizer = NewtonCG(GradientNormController(iteration_limit=10))
pos = 0.1*from_random(operator.domain)
minimizer = NewtonCG(GradientNormController(iteration_limit=10, name='findpos'))
for ii in range(3):
logger.info(f'Start iteration {ii+1}/3')
kl = MetricGaussianKL(pos, H, 3, mirror_samples=True)
kl, _ = minimizer(kl)
pos = kl.position
......
test/test_linearization.py
View file @ 773e7b8c
......@@ -57,14 +57,28 @@ def test_special_gradients():
'log', 'exp', 'sqrt', 'sin', 'cos', 'tan', 'sinc', 'sinh', 'cosh', 'tanh',
'absolute', 'reciprocal', 'sigmoid', 'log10', 'log1p', "expm1"
])
def test_actual_gradients(f):
@pmp('cplxpos', [True, False])
@pmp('cplxdir', [True, False])
@pmp('holomorphic', [True, False])
def test_actual_gradients(f, cplxpos, cplxdir, holomorphic):
if (cplxpos or cplxdir) and f in ['absolute']:
return
if holomorphic and f in ['absolute']:
# These function are not holomorphic
return
dom = ift.UnstructuredDomain((1,))
fld = ift.full(dom, 2.4)
eps = 1e-8
if cplxpos:
fld = fld + 0.21j
eps = 1e-7
if cplxdir:
eps *= 1j
if holomorphic:
eps *= (1+0.78j)
var0 = ift.Linearization.make_var(fld)
var1 = ift.Linearization.make_var(fld + eps)
f0 = var0.ptw(f).val.val
f1 = var1.ptw(f).val.val
df0 = (f1 - f0)/eps
df1 = _lin2grad(var0.ptw(f))
assert_allclose(df0, df1, rtol=100*eps)
assert_allclose(df0, df1, rtol=100*np.abs(eps))