Commit f1caff21 authored by Martin Reinecke's avatar Martin Reinecke

update to NIFTY5

parents ffc3fa61 a198e78a
Pipeline #31304 passed with stages
in 1 minute and 23 seconds
......@@ -34,7 +34,7 @@ test_python2_with_coverage:
script:
- python setup.py install --user -f
- mpiexec -n 2 --bind-to none nosetests -q 2> /dev/null
- nosetests -q --with-coverage --cover-package=nifty4 --cover-erase
- nosetests -q --with-coverage --cover-package=nifty5 --cover-erase
- >
coverage report --omit "*plotting*,*distributed_do*"
- >
......@@ -57,4 +57,4 @@ pages:
paths:
- public
only:
- NIFTy_4
- NIFTy_5
......@@ -49,8 +49,8 @@ Optional dependencies:
### Sources
The current version of Nifty4 can be obtained by cloning the repository and
switching to the NIFTy_4 branch:
The current version of Nifty5 can be obtained by cloning the repository and
switching to the NIFTy_5 branch:
git clone https://gitlab.mpcdf.mpg.de/ift/NIFTy.git
......@@ -59,7 +59,7 @@ switching to the NIFTy_4 branch:
In the following, we assume a Debian-based distribution. For other
distributions, the "apt" lines will need slight changes.
NIFTy4 and its mandatory dependencies can be installed via:
NIFTy5 and its mandatory dependencies can be installed via:
sudo apt-get install git libfftw3-dev python python-pip python-dev
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_4
......@@ -99,7 +99,7 @@ In oder to run the tests one needs two additional packages:
Afterwards the tests (including a coverage report) can be run using the
following command in the repository root:
nosetests -x --with-coverage --cover-html --cover-package=nifty4
nosetests -x --with-coverage --cover-html --cover-package=nifty5
### First Steps
......
......@@ -140,7 +140,7 @@
"source": [
"import numpy as np\n",
"np.random.seed(40)\n",
"import nifty4 as ift\n",
"import nifty5 as ift\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
......
import nifty4 as ift
from nifty4.library.nonlinearities import Linear, Tanh, Exponential
import nifty5 as ift
from nifty5.library.nonlinearities import Linear, Tanh, Exponential
import numpy as np
np.random.seed(42)
......@@ -14,8 +14,8 @@ def adjust_zero_mode(m0, t0):
return (ift.Field.from_global_data(m0.domain, mtmp),
ift.Field.from_global_data(t0.domain, ttmp))
if __name__ == "__main__":
if __name__ == "__main__":
noise_level = 1.
p_spec = (lambda k: (.5 / (k + 1) ** 3))
......
import nifty4 as ift
from nifty4.library.nonlinearities import Exponential
import nifty5 as ift
from nifty5.library.nonlinearities import Exponential
import numpy as np
np.random.seed(42)
......@@ -14,8 +14,8 @@ def adjust_zero_mode(m0, t0):
return (ift.Field.from_global_data(m0.domain, mtmp),
ift.Field.from_global_data(t0.domain, ttmp))
if __name__ == "__main__":
if __name__ == "__main__":
noise_level = 1.
p_spec = (lambda k: (1. / (k + 1) ** 2))
......
import nifty4 as ift
from nifty4.library.nonlinearities import Linear, Exponential, Tanh
import nifty5 as ift
from nifty5.library.nonlinearities import Linear, Exponential, Tanh
import numpy as np
np.random.seed(20)
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
if __name__ == "__main__":
signal_to_noise = 0.5 # The signal to noise ratio
......@@ -74,9 +74,12 @@ if __name__ == "__main__":
# Wiener filter
j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
sampling_ctrl = ift.GradientNormController(name="sampling",
tol_abs_gradnorm=1e2)
inverter = ift.ConjugateGradient(controller=ctrl)
sampling_inverter = ift.ConjugateGradient(controller=sampling_ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter)
S=S, N=N, R=R, inverter=inverter, sampling_inverter=sampling_inverter)
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
......@@ -88,7 +91,7 @@ if __name__ == "__main__":
ift.plot(data.cast_domain(plot_space), name='data.png', **plotdict)
ift.plot(m.cast_domain(plot_space), name='map.png', **plotdict)
# sampling the uncertainty map
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 10)
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 50)
ift.plot(ift.sqrt(variance).cast_domain(plot_space),
name="uncertainty.png", **plotdict)
ift.plot((mean+m).cast_domain(plot_space),
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
......@@ -48,9 +48,12 @@ if __name__ == "__main__":
# Wiener filter
j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-2)
sampling_ctrl = ift.GradientNormController(name="sampling",
tol_abs_gradnorm=2e1)
inverter = ift.ConjugateGradient(controller=ctrl)
sampling_inverter = ift.ConjugateGradient(controller=sampling_ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter)
S=S, N=N, R=R, inverter=inverter, sampling_inverter=sampling_inverter)
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
......@@ -60,6 +63,6 @@ if __name__ == "__main__":
ift.plot(m, name="map.png", **plotdict)
# sampling the uncertainty map
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 5)
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 50)
ift.plot(ift.sqrt(variance), name="uncertainty.png", **plotdict)
ift.plot(mean+m, name="posterior_mean.png", **plotdict)
# Program to generate figures of article "Information theory for fields"
# by Torsten Ensslin, Annalen der Physik, submitted to special edition
# "Physics of Information" in April 2018
# "Physics of Information" in April 2018, arXiv:1804.03350
# to get exact figure of paper comment out lines marked by "COMMENT OUT"
import numpy as np
import nifty4 as ift
import nifty5 as ift
import matplotlib.pyplot as plt
......@@ -25,7 +26,8 @@ if __name__ == "__main__":
sigma_n = 0.2 # noise level
sigma_n2 = sigma_n**2
L = 1. # Total length of interval or volume the field lives on
nprobes = 1000 # Number of probes for uncertainty quantification
nprobes = 1000 # Number of probes for uncertainty quantification used in paper
nprobes = 10 # COMMENT OUT TO REPRODUCE PAPER FIGURE EXACTLY
# Define resolution (pixels per dimension)
N_pixels = 1024
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
import matplotlib.pyplot as plt
from nifty4.sugar import create_power_operator
from nifty5.sugar import create_power_operator
np.random.seed(42)
......@@ -36,48 +36,63 @@ D_inv = ift.SandwichOperator.make(R_p, N.inverse) + S.inverse
N_samps = 200
N_iter = 100
IC = ift.GradientNormController(tol_abs_gradnorm=1e-3, iteration_limit=N_iter)
m, samps = ift.library.generate_krylov_samples(D_inv, S, j, N_samps, IC)
m_x = sky(m)
tol = 1e-3
IC = ift.GradientNormController(tol_abs_gradnorm=tol, iteration_limit=N_iter)
inverter = ift.ConjugateGradient(IC)
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, inverter=inverter,
sampling_inverter=inverter)
m_xi = curv.inverse_times(j)
samps_long = [curv.draw_sample(from_inverse=True) for i in range(N_samps)]
tol = 1e2
IC = ift.GradientNormController(tol_abs_gradnorm=tol, iteration_limit=N_iter)
inverter = ift.ConjugateGradient(IC)
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, inverter=inverter)
samps_old = [curv.draw_sample(from_inverse=True) for i in range(N_samps)]
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, inverter=inverter,
sampling_inverter=inverter)
samps_short = [curv.draw_sample(from_inverse=True) for i in range(N_samps)]
# Compute mean
sc = ift.StatCalculator()
for samp in samps_long:
sc.add(samp)
m_x = sky(sc.mean + m_xi)
plt.plot(d.to_global_data(), '+', label="data", alpha=.5)
plt.plot(s_x.to_global_data(), label="original")
plt.plot(m_x.to_global_data(), label="reconstruction")
plt.legend()
plt.savefig('Krylov_reconstruction.png')
plt.savefig('reconstruction.png')
plt.close()
pltdict = {'alpha': .3, 'linewidth': .2}
for i in range(N_samps):
if i == 0:
plt.plot(sky(samps_old[i]).to_global_data(), color='b',
label='Traditional samples (residuals)',
plt.plot(sky(samps_short[i]).to_global_data(), color='b',
label='Short samples (residuals)',
**pltdict)
plt.plot(sky(samps[i]).to_global_data(), color='r',
label='Krylov samples (residuals)',
plt.plot(sky(samps_long[i]).to_global_data(), color='r',
label='Long samples (residuals)',
**pltdict)
else:
plt.plot(sky(samps_old[i]).to_global_data(), color='b', **pltdict)
plt.plot(sky(samps[i]).to_global_data(), color='r', **pltdict)
plt.plot(sky(samps_short[i]).to_global_data(), color='b', **pltdict)
plt.plot(sky(samps_long[i]).to_global_data(), color='r', **pltdict)
plt.plot((s_x - m_x).to_global_data(), color='k', label='signal - mean')
plt.legend()
plt.savefig('Krylov_samples_residuals.png')
plt.savefig('samples_residuals.png')
plt.close()
D_hat_old = ift.full(x_space, 0.).to_global_data()
D_hat_new = ift.full(x_space, 0.).to_global_data()
for i in range(N_samps):
D_hat_old += sky(samps_old[i]).to_global_data()**2
D_hat_new += sky(samps[i]).to_global_data()**2
plt.plot(np.sqrt(D_hat_old / N_samps), 'r--', label='Traditional uncertainty')
D_hat_old += sky(samps_short[i]).to_global_data()**2
D_hat_new += sky(samps_long[i]).to_global_data()**2
plt.plot(np.sqrt(D_hat_old / N_samps), 'r--', label='Short uncertainty')
plt.plot(-np.sqrt(D_hat_old / N_samps), 'r--')
plt.fill_between(range(len(D_hat_new)), -np.sqrt(D_hat_new / N_samps), np.sqrt(
D_hat_new / N_samps), facecolor='0.5', alpha=0.5,
label='Krylov uncertainty')
label='Long uncertainty')
plt.plot((s_x - m_x).to_global_data(), color='k', label='signal - mean')
plt.legend()
plt.savefig('Krylov_uncertainty.png')
plt.savefig('uncertainty.png')
plt.close()
import nifty4 as ift
import nifty5 as ift
import numpy as np
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
# TODO: MAKE RESPONSE MPI COMPATIBLE OR USE LOS RESPONSE INSTEAD
......@@ -43,6 +43,7 @@ class CustomResponse(ift.LinearOperator):
def capability(self):
return self.TIMES | self.ADJOINT_TIMES
if __name__ == "__main__":
np.random.seed(43)
# Set up physical constants
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
if __name__ == "__main__":
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
import numericalunits as nu
......@@ -11,6 +11,7 @@ def init_nu():
data = comm.bcast(data, root=0)
nu.reset_units(data[0])
if __name__ == "__main__":
# In MPI mode, the random seed for numericalunits must be set by hand
init_nu()
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
np.random.seed(42)
......@@ -54,7 +54,8 @@ if __name__ == "__main__":
# Initialize Wiener filter energy
energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S,
inverter=inverter)
inverter=inverter,
sampling_inverter=inverter)
energy, convergence = minimizer(energy)
m = energy.position
......
rm -rf docs/build docs/source/mod
sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty4
sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty5
sphinx-build -b html docs/source/ docs/build/
......@@ -7,7 +7,7 @@ NIFTy-related publications
author = {{Martin Reinecke, Theo Steininger, Marco Selig}},
title = {NIFTy -- Numerical Information Field TheorY},
url = {https://gitlab.mpcdf.mpg.de/ift/NIFTy},
version = {nifty4},
version = {nifty5},
date = {2018-04-05},
}
......
.. currentmodule:: nifty4
.. currentmodule:: nifty5
=============
Code Overview
......@@ -35,7 +35,7 @@ Domains
Abstract base class
-------------------
One of the fundamental building blocks of the NIFTy4 framework is the *domain*.
One of the fundamental building blocks of the NIFTy5 framework is the *domain*.
Its required capabilities are expressed by the abstract :class:`Domain` class.
A domain must be able to answer the following queries:
......@@ -151,7 +151,7 @@ performance benefits.
Linear Operators
================
A linear operator (represented by NIFTy4's abstract :class:`LinearOperator`
A linear operator (represented by NIFTy5's abstract :class:`LinearOperator`
class) can be interpreted as an (implicitly defined) matrix.
It can be applied to :class:`Field` instances, resulting in other :class:`Field`
instances that potentially live on other domains.
......@@ -215,7 +215,7 @@ Further operator classes provided by NIFTy are
Syntactic sugar
---------------
Nifty4 allows simple and intuitive construction of altered and combined
Nifty5 allows simple and intuitive construction of altered and combined
operators.
As an example, if ``A``, ``B`` and ``C`` are of type :class:`LinearOperator`
and ``f1`` and ``f2`` are of type :class:`Field`, writing::
......@@ -245,7 +245,7 @@ high-dimensional functions, which are often nonlinear.
Energy functionals
------------------
In NIFTy4 such functions are represented by objects of type :class:`Energy`.
In NIFTy5 such functions are represented by objects of type :class:`Energy`.
These hold the prescription how to calculate the function's
:attr:`~Energy.value`, :attr:`~Energy.gradient` and
(optionally) :attr:`~Energy.curvature` at any given :attr:`~Energy.position`
......@@ -256,9 +256,9 @@ linear operator objects.
Energies are classes that typically have to be provided by the user when
tackling new IFT problems.
Some examples of concrete energy classes delivered with NIFTy4 are
Some examples of concrete energy classes delivered with NIFTy5 are
:class:`QuadraticEnergy` (with position-independent curvature, mainly used with
conjugate gradient minimization) and :class:`~nifty4.library.WienerFilterEnergy`.
conjugate gradient minimization) and :class:`~nifty5.library.WienerFilterEnergy`.
Iteration control
......@@ -269,7 +269,7 @@ checking the quality of the current solution estimate and stopping once
it is sufficiently accurate. In case of numerical problems, the iteration needs
to be terminated as well, returning a suitable error description.
In NIFTy4, this functionality is encapsulated in the abstract
In NIFTy5, this functionality is encapsulated in the abstract
:class:`IterationController` class, which is provided with the initial energy
object before starting the minimization, and is updated with the improved
energy after every iteration. Based on this information, it can either continue
......
......@@ -73,7 +73,7 @@ source_suffix = '.rst'
master_doc = 'index'
# General information about the project.
project = u'NIFTy4'
project = u'NIFTy5'
copyright = u'2013-2018, Max-Planck-Society'
author = u'Theo Steininger / Martin Reinecke'
......
......@@ -5,7 +5,7 @@ Installation
In the following, we assume a Debian-based Linux distribution. For other
distributions, the "apt" lines will need slight changes.
NIFTy4 and its mandatory dependencies can be installed via::
NIFTy5 and its mandatory dependencies can be installed via::
sudo apt-get install git libfftw3-dev python python-pip python-dev
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_4
......
# 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-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from ..minimization.quadratic_energy import QuadraticEnergy
def generate_krylov_samples(D_inv, S, j, N_samps, controller):
"""
Generates inverse samples from a curvature D.
This algorithm iteratively generates samples from
a curvature D by applying conjugate gradient steps
and resampling the curvature in search direction.
Parameters
----------
D_inv : EndomorphicOperator
The curvature which will be the inverse of the covarianc
of the generated samples
S : EndomorphicOperator (from which one can sample)
A prior covariance operator which is used to generate prior
samples that are then iteratively updated
j : Field, optional
A Field to which the inverse of D_inv is applied. The solution
of this matrix inversion problem is a side product of generating
the samples.
If not supplied, it is sampled from the inverse prior.
N_samps : Int
How many samples to generate.
controller : IterationController
convergence controller for the conjugate gradient iteration
Returns
-------
(solution, samples) : A tuple of a field 'solution' and a list of fields
'samples'. The first entry of the tuple is the solution x to
D_inv(x) = j
and the second entry are a list of samples from D_inv.inverse
"""
# RL FIXME: make consistent with complex numbers
j = S.draw_sample(from_inverse=True) if j is None else j
energy = QuadraticEnergy(j.empty_copy().fill(0.), D_inv, j)
y = [S.draw_sample() for _ in range(N_samps)]
status = controller.start(energy)
if status != controller.CONTINUE:
return energy.position, y
r = energy.gradient
d = r.copy()
previous_gamma = r.vdot(r).real
if previous_gamma == 0:
return energy.position, y
while True:
q = energy.curvature(d)
ddotq = d.vdot(q).real
if ddotq == 0.:
logger.error("Error: ConjugateGradient: ddotq==0.")
return energy.position, y
alpha = previous_gamma/ddotq
if alpha < 0:
logger.error("Error: ConjugateGradient: alpha<0.")
return energy.position, y
for i in range(len(y)):
y[i] += (np.random.randn()*np.sqrt(ddotq) - y[i].vdot(q))/ddotq * d
q *= -alpha
r = r + q
energy = energy.at_with_grad(energy.position - alpha*d, r)
gamma = r.vdot(r).real
if gamma == 0:
return energy.position, y
status = controller.check(energy)
if status != controller.CONTINUE:
return energy.position, y
d *= max(0, gamma/previous_gamma)
d += r
previous_gamma = gamma
......@@ -4,6 +4,7 @@ from . import dobj
from .domains import *
from .domain_tuple import DomainTuple
from .field import Field
from .models import *
from .operators import *
from .probing.utils import probe_with_posterior_samples, probe_diagonal, \
StatCalculator
......@@ -20,7 +21,5 @@ from .logger import logger
from .multi import *
__all__ = ["__version__", "dobj", "DomainTuple"] + \
domains.__all__ + operators.__all__ + minimization.__all__ + \
["Field"] + sugar.__all__ + \
multi.__all__
# We deliberately don't set __all__ here, because we don't want people to do a
# "from nifty5 import *"; that would swamp the global namespace.
......@@ -211,6 +211,7 @@ class data_object(object):
def fill(self, value):
self._data.fill(value)
for op in ["__add__", "__radd__", "__iadd__",
"__sub__", "__rsub__", "__isub__",
"__mul__", "__rmul__", "__imul__",
......
......@@ -25,7 +25,7 @@ class GLSpace(StructuredDomain):
"""NIFTy subclass for Gauss-Legendre pixelizations of the two-sphere.
Its harmonic partner domain is the
:class:`~nifty4.domains.lm_space.LMSpace`.
:class:`~nifty5.domains.lm_space.LMSpace`.
Parameters
----------
......
......@@ -25,7 +25,7 @@ class HPSpace(StructuredDomain):
"""NIFTy subclass for HEALPix discretizations of the two-sphere.
Its harmonic partner domain is the
:class:`~nifty4.domains.lm_space.LMSpace`.
:class:`~nifty5.domains.lm_space.LMSpace`.
Parameters
----------
......
......@@ -25,8 +25,8 @@ from ..field import Field
class LMSpace(StructuredDomain):
"""NIFTy subclass for sets of spherical harmonic coefficients.
Its harmonic partner spaces are :class:`~nifty4.domains.hp_space.HPSpace`
and :class:`~nifty4.domains.gl_space.GLSpace`.
Its harmonic partner spaces are :class:`~nifty5.domains.hp_space.HPSpace`
and :class:`~nifty5.domains.gl_space.GLSpace`.
Parameters
----------
......@@ -130,7 +130,7 @@ class LMSpace(StructuredDomain):
return self._mmax
def get_default_codomain(self):
"""Returns a :class:`~nifty4.domains.gl_space.GLSpace` object, which is
"""Returns a :class:`~nifty5.domains.gl_space.GLSpace` object, which is
capable of storing an accurate representation of data residing on
`self`.
......
......@@ -21,7 +21,7 @@ from functools import reduce
class UnstructuredDomain(Domain):
"""A :class:`~nifty4.domains.domain.Domain` subclass for spaces with no
"""A :class:`~nifty5.domains.domain.Domain` subclass for spaces with no
associated geometry.
Typically used for data spaces.
......
......@@ -23,7 +23,6 @@ from . import utilities
from .domain_tuple import DomainTuple
from functools import reduce
from . import dobj
import sys
__all__ = ["Field", "sqrt", "exp", "log", "conjugate"]
......@@ -714,10 +713,10 @@ class Field(object):
self.local_data[()] = other.local_data[()]
def __repr__(self):
return "<nifty4.Field>"
return "<nifty5.Field>"
def __str__(self):
return "nifty4.Field instance\n- domain = " + \
return "nifty5.Field instance\n- domain = " + \
self._domain.__str__() + \
"\n- val = " + repr(self.val)
......@@ -730,6 +729,7 @@ class Field(object):
return False
return (self._val == other._val).all()
for op in ["__add__", "__radd__", "__iadd__",
"__sub__", "__rsub__", "__isub__",
"__mul__", "__rmul__", "__imul__",
......
from .wiener_filter_energy import WienerFilterEnergy
from .wiener_filter_curvature import WienerFilterCurvature
from .los_response import LOSResponse
from .noise_energy import NoiseEnergy
from .nonlinear_power_energy import NonlinearPowerEnergy
from .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
from .nonlinearities import Exponential, Linear, PositiveTanh, Tanh
from .poisson_energy import PoissonEnergy
from .nonlinearities import Exponential, Linear, Tanh, PositiveTanh
from .krylov_sampling import generate_krylov_samples
from .los_response import LOSResponse
from .unit_log_gauss import UnitLogGauss
from .poisson_log_likelihood import PoissonLogLikelihood
from .wiener_filter_curvature import WienerFilterCurvature
from .wiener_filter_energy import WienerFilterEnergy
......@@ -122,8 +122,10 @@ class LOSResponse(LinearOperator):
starts = np.array(starts)
nlos = starts.shape[1]
ends = np.array(ends)
assert starts.shape[0] == ndim, "dimension mismatch"