Commit 110a3ce5 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_4' into notebook_demo

parents 81cab5d6 9a5abf11
......@@ -17,10 +17,11 @@ test_min:
stage: test
script:
- pip install --user .
- nosetests -q --with-coverage --cover-package=nifty4 --cover-branches
- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests -q
- pip3 install --user .
- nosetests -q
- nosetests3 -q
- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests3 -q
- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests -q 2>/dev/null
- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests3 -q 2>/dev/null
- nosetests -q --with-coverage --cover-package=nifty4 --cover-branches --cover-erase
- >
coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }' || true
coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }'
......@@ -37,13 +37,11 @@ if __name__ == "__main__":
R = R*ift.DiagonalOperator(mask)
R = R*ht
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
#R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
# sensitivity=(mask,))
data_domain = R.target[0]
# Setting up the noise covariance and drawing a random noise realization
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.std()/signal_to_noise
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.DiagonalOperator(
ift.Field.full(data_domain, noise_amplitude**2))
noise = ift.Field.from_random(
......@@ -58,16 +56,15 @@ if __name__ == "__main__":
inverter = ift.ConjugateGradient(controller=ctrl)
energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R,
N, S, inverter=inverter)
# minimizer = ift.VL_BFGS(controller=ctrl2, max_history_length=20)
#minimizer = ift.VL_BFGS(controller=ctrl2, max_history_length=20)
minimizer = ift.RelaxedNewton(controller=ctrl2)
# minimizer = ift.SteepestDescent(controller=ctrl2)
#minimizer = ift.SteepestDescent(controller=ctrl2)
me = minimizer(energy)
m = ht(me[0].position)
# Plotting
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
plotdict = {"colormap": "Planck-like"}
ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict)
logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
ift.plot(ift.Field(signal_space, val=ift.dobj.from_global_data(logdata)),
......
......@@ -24,10 +24,10 @@ if __name__ == "__main__":
# Set up position space
# s_space = ift.RGSpace([1024])
s_space = ift.HPSpace(32)
h_space = s_space.get_default_codomain()
# Define harmonic transformation and associated harmonic space
FFT = ift.FFTOperator(s_space)
h_space = FFT.target[0]
HT = ift.HarmonicTransformOperator(h_space, target=s_space)
# Setting up power space
p_space = ift.PowerSpace(h_space,
......@@ -51,22 +51,28 @@ if __name__ == "__main__":
mask = ift.Field(s_space, val=ift.dobj.from_global_data(mask))
MaskOperator = ift.DiagonalOperator(mask)
InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0],
exposure=[1.0])
MeasurementOperator = InstrumentResponse*MaskOperator
R = ift.GeometryRemover(s_space)
R = R*MaskOperator
#R = R*HT
#R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
MeasurementOperator = R
d_space = MeasurementOperator.target
noise_covariance = ift.Field(d_space, val=noise_level**2).weight()
N = ift.DiagonalOperator(noise_covariance)
n = ift.Field.from_random(domain=d_space, random_type='normal',
std=noise_level)
Projection = ift.PowerProjectionOperator(domain=h_space,
power_space=p_space)
power = Projection.adjoint_times(ift.exp(0.5*log_p))
# Creating the mock data
true_sky = nonlinearity(FFT.adjoint_times(power*sh))
d = MeasurementOperator(true_sky) + n
true_sky = nonlinearity(HT(power*sh))
noiseless_data = MeasurementOperator(true_sky)
noise_amplitude = noiseless_data.val.std()*noise_level
N = ift.DiagonalOperator(
ift.Field.full(d_space, noise_amplitude**2))
n = ift.Field.from_random(
domain=d_space, random_type='normal',
std=noise_amplitude, mean=0)
# Creating the mock data
d = noiseless_data + n
m0 = ift.power_synthesize(ift.Field(p_space, val=1e-7))
t0 = ift.Field(p_space, val=-4.)
......@@ -84,7 +90,7 @@ if __name__ == "__main__":
for i in range(20):
power0 = Projection.adjoint_times(ift.exp(0.5*t0))
map0_energy = ift.library.NonlinearWienerFilterEnergy(
m0, d, MeasurementOperator, nonlinearity, FFT, power0, N, S,
m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S,
inverter=inverter)
# Minimization with chosen minimizer
......@@ -96,7 +102,7 @@ if __name__ == "__main__":
# Initializing the power energy with updated parameters
power0_energy = ift.library.NonlinearPowerEnergy(
position=t0, d=d, N=N, m=m0, D=D0, FFT=FFT,
position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Projection=Projection, sigma=1., samples=2, inverter=inverter)
......@@ -110,6 +116,6 @@ if __name__ == "__main__":
m0, t0 = adjust_zero_mode(m0, t0)
ift.plot(true_sky)
ift.plot(nonlinearity(FFT.adjoint_times(power0*m0)),
ift.plot(nonlinearity(HT(power0*m0)),
title='reconstructed_sky')
ift.plot(MeasurementOperator.adjoint_times(d))
......@@ -81,10 +81,10 @@ if __name__ == "__main__":
data_domain = R.target
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.std()/signal_to_noise
noise_amplitude = noiseless_data.val.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, noise_amplitude**2)
N = ift.DiagonalOperator(ndiag)
#ndiag = ift.Field.full(data_domain, noise_amplitude**2)
N = ift.ScalingOperator(noise_amplitude**2, data_domain)
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=noise_amplitude, mean=0)
......@@ -108,12 +108,5 @@ if __name__ == "__main__":
ift.plot(ift.Field(plot_space,val=data.val), name='data.png', **plotdict)
ift.plot(ift.Field(plot_space,val=m.val), name='map.png', **plotdict)
# sampling the uncertainty map
sample_variance = ift.Field.zeros(signal_domain)
sample_mean = ift.Field.zeros(signal_domain)
n_samples = 10
for i in range(n_samples):
sample = ht(wiener_curvature.generate_posterior_sample()) + m
sample_variance += sample**2
sample_mean += sample
variance = sample_variance/n_samples - (sample_mean/n_samples)**2
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 10)
ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict)
......@@ -42,7 +42,7 @@ if __name__ == "__main__":
data_domain = R.target[0]
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.std()/signal_to_noise
noise_amplitude = noiseless_data.val.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, noise_amplitude**2)
N = ift.DiagonalOperator(ndiag)
......
......@@ -42,7 +42,7 @@ if __name__ == "__main__":
noiseless_data = R(sh)
signal_to_noise = 1.
noise_amplitude = noiseless_data.std()/signal_to_noise
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2))
n = ift.Field.from_random(domain=s_space,
random_type='normal',
......
......@@ -58,7 +58,7 @@ if __name__ == "__main__":
data_domain = R.target[0]
noiseless_data = R(mock_signal)
noise_amplitude = noiseless_data.std()/signal_to_noise
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.DiagonalOperator(
ift.Field.full(data_domain, noise_amplitude**2))
noise = ift.Field.from_random(
......
......@@ -36,7 +36,7 @@ if __name__ == "__main__":
R = Instrument*ht
noiseless_data = R(sh)
signal_to_noise = 1.
noise_amplitude = noiseless_data.std()/signal_to_noise
noise_amplitude = noiseless_data.val.std()/signal_to_noise
N = ift.DiagonalOperator(ift.Field.full(s_space, noise_amplitude**2))
n = ift.Field.from_random(domain=s_space,
random_type='normal',
......@@ -68,12 +68,6 @@ if __name__ == "__main__":
sample_variance = ift.Field.zeros(s_space)
sample_mean = ift.Field.zeros(s_space)
n_samples = 50
for i in range(n_samples):
sample = ht(D.generate_posterior_sample() + m)
sample_variance += sample**2
sample_mean += sample
sample_mean /= n_samples
sample_variance /= n_samples
variance = sample_variance - sample_mean**2
mean, variance = ift.probe_with_posterior_samples(D, ht, 50)
ift.plot(variance, name="variance.png", colormap="Planck-like")
ift.plot(mean, name="mean.png", colormap="Planck-like")
Significant differences between NIFTy nightly and nifty2go
==========================================================
Significant differences between NIFTy3 and nifty4
=================================================
1) Field domains in nifty2go are stored in DomainTuple objects, which makes
1) Field domains in nifty4 are stored in DomainTuple objects, which makes
comparisons between domains and computation of axis indices etc. much simpler
and more efficient.
No impact on the user ... these objects are generated whenever needed and
have all necessary functions to make them look like tuples of spaces.
2) In nifty2go an operator's domain and target refer to the _full_ domains of
2) In nifty4 an operator's domain and target refer to the _full_ domains of
the input and output fields read/written by times(), adjoint_times() etc.
In NIFTy nightly, domain and target only refer to the (sub-)domain on
which the operator actually acts. This leads to complications like the need
for the "default_spaces" argument in the operator constructor and the
"spaces" keywords in all operator calls.
Advantages of the nifty2go approach:
Advantages of the nifty4 approach:
- less error-prone and easier to understand; less code overall
- operators have more knowledge and may tune themselves better
- difficulties with the design of ComposedOperator (issue 152) resolve
......@@ -25,13 +25,13 @@ Significant differences between NIFTy nightly and nifty2go
However, I have not found any such situation in the current code base, so
it appears to be rare.
3) nifty2go uses one of two different "data_object" modules for array
3) nifty4 uses one of two different "data_object" modules for array
storage instead of D2O.
A "data_object" module consists of a class called "data_object" which
provides a subset of the numpy.ndarray interface, plus a few additional
functions for manipulating these data objects.
If no MPI support is found on the system, or if a computation is run on a
single task, nifty2go automatically loads a minimalistic "data_object"
single task, nifty4 automatically loads a minimalistic "data_object"
module where the data_object class is simply identical to numpy.ndarray.
The support functions are mostly trivial as well.
If MPI is required, another module is loaded, which supports parallel
......@@ -53,11 +53,11 @@ Significant differences between NIFTy nightly and nifty2go
kindex -> k_lengths
(because this is not an index)
6) In nifty2go, PowerSpace is not a harmonic space.
6) In nifty4, PowerSpace is not a harmonic space.
7) In nifty2go, parallel probing should work (needs systematic testing)
7) In nifty4, parallel probing should work (needs systematic testing)
9) Many default arguments have been removed in nifty2go, wherever there is no
9) Many default arguments have been removed in nifty4, wherever there is no
sensible default (in my opinion). My personal impression is that this has
actually made the demos more readable, but I'm sure not everyone will agree
:)
......@@ -80,8 +80,16 @@ Significant differences between NIFTy nightly and nifty2go
more or less always needed).
14) A new approach is used for FFTs along axes that are distributed among
MPI tasks. As a consequence, nifty2go works well with the standard version
MPI tasks. As a consequence, nifty4 works well with the standard version
of pyfftw and does not need the MPI-enabled fork.
15) Arithmetic functions working on Fields have been moved from
basic_arithmetics.py to field.py.
16) Operators can be comined via "*", "+" and "-", resulting in new combined
operators.
17) Every operator has the properties ".adjoint" and ".inverse", which return
its adjoint and inverse, respectively.
18) Handling of volume factors has been changed completely.
......@@ -2,10 +2,9 @@ from .version import __version__
from . import dobj
from .domain_object import DomainObject
from .spaces.field_array import FieldArray
from .spaces.space import Space
from .spaces.domain import Domain
from .spaces.unstructured_domain import UnstructuredDomain
from .spaces.structured_domain import StructuredDomain
from .spaces.rg_space import RGSpace
from .spaces.lm_space import LMSpace
from .spaces.hp_space import HPSpace
......@@ -23,8 +22,7 @@ from .operators.harmonic_transform_operator import HarmonicTransformOperator
from .operators.fft_operator import FFTOperator
from .operators.fft_smoothing_operator import FFTSmoothingOperator
from .operators.direct_smoothing_operator import DirectSmoothingOperator
from .operators.response_operator import ResponseOperator
from .operators.response_operator import GeometryRemover
from .operators.geometry_remover import GeometryRemover
from .operators.laplace_operator import LaplaceOperator
from .operators.power_projection_operator import PowerProjectionOperator
from .operators.inversion_enabler import InversionEnabler
......@@ -34,6 +32,7 @@ from .field import Field, sqrt, exp, log
from .probing.prober import Prober
from .probing.diagonal_prober_mixin import DiagonalProberMixin
from .probing.trace_prober_mixin import TraceProberMixin
from .probing.utils import probe_with_posterior_samples
from .minimization.line_search import LineSearch
from .minimization.line_search_strong_wolfe import LineSearchStrongWolfe
......@@ -55,11 +54,11 @@ from .sugar import *
from .plotting.plot import plot
from . import library
__all__ = ["DomainObject", "FieldArray", "Space", "RGSpace", "LMSpace",
__all__ = ["Domain", "UnstructuredDomain", "StructuredDomain", "RGSpace", "LMSpace",
"HPSpace", "GLSpace", "DOFSpace", "PowerSpace", "DomainTuple",
"LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "FFTOperator", "FFTSmoothingOperator",
"DirectSmoothingOperator", "ResponseOperator", "LaplaceOperator",
"DirectSmoothingOperator", "LaplaceOperator",
"PowerProjectionOperator", "InversionEnabler",
"Field", "sqrt", "exp", "log",
"Prober", "DiagonalProberMixin", "TraceProberMixin"]
......@@ -17,7 +17,7 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from functools import reduce
from .domain_object import DomainObject
from .spaces.domain import Domain
class DomainTuple(object):
......@@ -28,12 +28,10 @@ class DomainTuple(object):
self._axtuple = self._get_axes_tuple()
shape_tuple = tuple(sp.shape for sp in self._dom)
self._shape = reduce(lambda x, y: x + y, shape_tuple, ())
self._dim = reduce(lambda x, y: x * y, self._shape, 1)
self._accdims = (1,)
self._size = reduce(lambda x, y: x * y, self._shape, 1)
prod = 1
for dom in self._dom:
prod *= dom.dim
self._accdims += (prod,)
prod *= dom.size
def _get_axes_tuple(self):
i = 0
......@@ -60,16 +58,16 @@ class DomainTuple(object):
def _parse_domain(domain):
if domain is None:
return ()
if isinstance(domain, DomainObject):
if isinstance(domain, Domain):
return (domain,)
if not isinstance(domain, tuple):
domain = tuple(domain)
for d in domain:
if not isinstance(d, DomainObject):
if not isinstance(d, Domain):
raise TypeError(
"Given object contains something that is not an "
"instance of DomainObject class.")
"instance of Domain class.")
return domain
def __getitem__(self, i):
......@@ -80,8 +78,8 @@ class DomainTuple(object):
return self._shape
@property
def dim(self):
return self._dim
def size(self):
return self._size
@property
def axes(self):
......
......@@ -35,7 +35,7 @@ class Field(object):
Parameters
----------
domain : None, DomainTuple, tuple of DomainObjects, or single DomainObject
domain : None, DomainTuple, tuple(Domain), or Domain
val : None, Field, data_object, or scalar
The values the array should contain after init. A scalar input will
......@@ -156,7 +156,7 @@ class Field(object):
'pm1', 'normal', 'uniform' are the supported arguments for this
method.
domain : DomainObject
domain : DomainTuple
The domain of the output random field
dtype : type
......@@ -201,7 +201,7 @@ class Field(object):
return self._domain.shape
@property
def dim(self):
def size(self):
""" Returns the total number of pixel-dimensions the field has.
Effectively, all values from shape are multiplied.
......@@ -211,13 +211,13 @@ class Field(object):
out : int
The dimension of the Field.
"""
return self._domain.dim
return self._domain.size
@property
def real(self):
""" The real part of the field (data is not copied)."""
if not np.issubdtype(self.dtype, np.complexfloating):
raise ValueError(".real called on a non-complex Field")
return self
return Field(self._domain, self.val.real)
@property
......
......@@ -74,32 +74,8 @@ class WienerFilterCurvature(EndomorphicOperator):
covariance.
"""
power = power_analyze(sqrt(self.S.diagonal()))
mock_signal = power_synthesize(power, real_signal=True)
noise = self.N.diagonal()
mock_noise = Field.from_random(random_type="normal",
domain=self.N.domain, dtype=noise.dtype)
mock_noise *= sqrt(noise)
mock_data = self.R(mock_signal) + mock_noise
mock_j = self.R.adjoint_times(self.N.inverse_times(mock_data))
mock_m = self.inverse_times(mock_j)
return mock_signal - mock_m
def generate_posterior_sample2(self):
power = self.S.diagonal()
mock_signal = Field.from_random(random_type="normal",
domain=self.S.domain,
dtype=power.dtype)
mock_signal *= sqrt(power)
noise = self.N.diagonal()
mock_noise = Field.from_random(random_type="normal",
domain=self.N.domain, dtype=noise.dtype)
mock_noise *= sqrt(noise)
mock_signal = self.S.generate_posterior_sample()
mock_noise = self.N.generate_posterior_sample()
mock_data = self.R(mock_signal) + mock_noise
......
......@@ -18,7 +18,7 @@
from __future__ import division
import numpy as np
from ..field import Field
from ..field import Field, sqrt
from ..domain_tuple import DomainTuple
from .endomorphic_operator import EndomorphicOperator
from .. import utilities
......@@ -140,3 +140,27 @@ class DiagonalOperator(EndomorphicOperator):
def adjoint(self):
return DiagonalOperator(self._diagonal.conjugate(), self._domain,
self._spaces)
def generate_posterior_sample(self):
""" Generates a posterior sample from a Gaussian distribution with
given mean and covariance.
This method generates samples by setting up the observation and
reconstruction of a mock signal in order to obtain residuals of the
right correlation which are added to the given mean.
Returns
-------
sample : Field
Returns the a sample from the Gaussian of given mean and
covariance.
"""
if self._spaces is not None:
raise ValueError("Cannot draw (yet) from this operator")
res = Field.from_random(random_type="normal",
domain=self._domain,
dtype=self._diagonal.dtype)
res *= sqrt(self._diagonal)
return res
......@@ -17,20 +17,16 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..field import Field
from ..spaces.field_array import FieldArray
from ..spaces.unstructured_domain import UnstructuredDomain
from ..domain_tuple import DomainTuple
from .linear_operator import LinearOperator
from .fft_smoothing_operator import FFTSmoothingOperator
from .diagonal_operator import DiagonalOperator
from .scaling_operator import ScalingOperator
import numpy as np
class GeometryRemover(LinearOperator):
def __init__(self, domain):
super(GeometryRemover, self).__init__()
self._domain = DomainTuple.make(domain)
target_list = [FieldArray(dom.shape) for dom in self._domain]
target_list = [UnstructuredDomain(dom.shape) for dom in self._domain]
self._target = DomainTuple.make(target_list)
@property
......@@ -50,38 +46,3 @@ class GeometryRemover(LinearOperator):
if mode == self.TIMES:
return Field(self._target, val=x.weight(1).val)
return Field(self._domain, val=x.val).weight(1)
def ResponseOperator(domain, sigma, sensitivity):
# sensitivity has units 1/field/volume and gives a measure of how much
# the instrument will excited when it is exposed to a certain field
# volume amplitude
domain = DomainTuple.make(domain)
ncomp = len(sensitivity)
if len(sigma) != ncomp or len(domain) != ncomp:
raise ValueError("length mismatch between sigma, sensitivity "
"and domain")
ncomp = len(sigma)
if ncomp == 0:
raise ValueError("Empty response operator not allowed")
kernel = None
sensi = None
for i in range(ncomp):
if sigma[i] > 0:
op = FFTSmoothingOperator(domain, sigma[i], space=i)
kernel = op if kernel is None else op*kernel
if np.isscalar(sensitivity[i]):
if sensitivity[i] != 1.:
op = ScalingOperator(sensitivity[i], domain)
sensi = op if sensi is None else op*sensi
elif isinstance(sensitivity[i], Field):
op = DiagonalOperator(sensitivity[i], domain=domain, spaces=i)
sensi = op if sensi is None else op*sensi
res = GeometryRemover(domain)
if sensi is not None:
res = res * sensi
if kernel is not None:
res = res * kernel
return res
......@@ -88,3 +88,23 @@ class ScalingOperator(EndomorphicOperator):
return self.TIMES | self.ADJOINT_TIMES
return (self.TIMES | self.ADJOINT_TIMES |
self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES)
def generate_posterior_sample(self):
""" Generates a posterior sample from a Gaussian distribution with
given mean and covariance.
This method generates samples by setting up the observation and
reconstruction of a mock signal in order to obtain residuals of the
right correlation which are added to the given mean.
Returns
-------
sample : Field
Returns the a sample from the Gaussian of given mean and
covariance.
"""
return Field.from_random(random_type="normal",
domain=self._domain,
std=np.sqrt(self._factor),
dtype=np.result_type(self._factor))
......@@ -16,33 +16,40 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import unittest
from numpy.testing import assert_allclose
import nifty4 as ift
from itertools import product
from test.common import expand
class ResponseOperator_Tests(unittest.TestCase):
spaces = [ift.RGSpace(128), ift.GLSpace(nlat=37)]
@expand(product(spaces, [0., 5., 1.], [0., 1., .33]))
def test_property(self, space, sigma, sensitivity):
if not isinstance(space, ift.RGSpace): # no smoothing supported
sigma = 0.
op = ift.ResponseOperator(space, sigma=[sigma],
sensitivity=[sensitivity])
if op.domain[0] != space:
raise TypeError