...
 
Commits (105)
......@@ -26,7 +26,6 @@ test_min:
test_mpi:
stage: test
script:
- ci/install_pyHealpix.sh
- ci/install_mpi4py.sh
- nosetests -vv
- nosetests3 -vv
......@@ -34,7 +33,6 @@ test_mpi:
test_mpi_fftw:
stage: test
script:
- ci/install_pyHealpix.sh
- ci/install_mpi4py.sh
- ci/install_pyfftw.sh
- nosetests -vv
......@@ -43,14 +41,15 @@ test_mpi_fftw:
test_mpi_fftw_hdf5:
stage: test
script:
- ci/install_pyHealpix.sh
- ci/install_mpi4py.sh
- ci/install_pyfftw.sh
- ci/install_h5py.sh
- mpiexec --allow-run-as-root -n 2 nosetests -x
- mpiexec --allow-run-as-root -n 2 nosetests3 -x
- mpiexec --allow-run-as-root -n 4 nosetests -x
- mpiexec --allow-run-as-root -n 4 nosetests3 -x
- nosetests -vv
- nosetests3 -vv
#- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 2 nosetests -x
#- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests -x
#- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 2 nosetests3 -x
#- OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests3 -x
- nosetests -x --with-coverage --cover-package=nifty --cover-branches
- >
coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }'
......
......@@ -14,10 +14,7 @@ ADD ci/requirements_extras.txt /tmp/requirements_extras.txt
RUN pip install --upgrade -r /tmp/requirements_extras.txt
# install pyHealpix, pyfftw and h5py
ADD ci/install_pyHealpix.sh /tmp/install_pyHealpix.sh
RUN cd /tmp && chmod +x install_pyHealpix.sh && ./install_pyHealpix.sh
# install pyfftw and h5py
ADD ci/install_mpi4py.sh /tmp/install_mpi4py.sh
RUN cd /tmp && chmod +x install_mpi4py.sh && ./install_mpi4py.sh
......
NIFTY - Numerical Information Field Theory
==========================================
[![build status](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/master/build.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/master)
[![coverage report](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/master/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/master)
[![build status](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_3/build.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_3)
[![coverage report](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_3/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_3)
**NIFTY** project homepage:
[http://www.mpa-garching.mpg.de/ift/nifty/](http://www.mpa-garching.mpg.de/ift/nifty/)
......@@ -91,8 +91,7 @@ Starting with a fresh Ubuntu installation move to a folder like
- Install pyHealpix:
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
(cd pyHealpix && autoreconf -i && ./configure --prefix=$HOME/.local --enable-openmp --enable-native-optimizations && make -j4 install)
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
- Finally, NIFTy:
......@@ -123,14 +122,7 @@ may cause trouble.
- Install pyHealpix:
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
cd pyHealpix
autoreconf -i && ./configure --prefix=`python-config --prefix` --enable-openmp --enable-native-optimizations && make -j4 && sudo make install
cd ..
(The third command installs the package system-wide. User-specific
installation would be preferrable, but we haven't found a simple recipe yet
how to determine the installation prefix ...)
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
- Install NIFTy:
......@@ -164,7 +156,7 @@ Please, acknowledge the use of NIFTY in your publication(s) by using a
phrase such as the following:
> *"Some of the results in this publication have been derived using the
> NIFTY package [Selig et al., 2013]."*
> NIFTY package [Steininger et al., 2017]."*
### References
......@@ -180,7 +172,5 @@ The NIFTY package is licensed under the
**NIFTY** project homepage:
[](http://www.mpa-garching.mpg.de/ift/nifty/)
[1] Selig et al., "NIFTY - Numerical Information Field Theory - a
versatile Python library for signal inference", [A&A, vol. 554, id.
A26](http://dx.doi.org/10.1051/0004-6361/201321236), 2013;
[arXiv:1301.4499](http://www.arxiv.org/abs/1301.4499)
[1] Steininger et al., "NIFTy 3 - Numerical Information Field Theory - A Python framework for multicomponent signal inference on HPC clusters", 2017, submitted to PLOS One;
[arXiv:1708.01073](https://arxiv.org/abs/1708.01073)
......@@ -2,3 +2,4 @@
apt-get install -y libhdf5-10 libhdf5-dev libhdf5-openmpi-10 libhdf5-openmpi-dev hdf5-tools
CC=mpicc HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/openmpi HDF5_MPI="ON" pip install --no-binary=h5py h5py
CC=mpicc HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/openmpi HDF5_MPI="ON" pip3 install --no-binary=h5py h5py
#!/bin/bash
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
(cd pyHealpix && autoreconf -i && ./configure --enable-openmp && make -j4 install)
(cd pyHealpix && autoreconf -i && PYTHON=python3 ./configure --enable-openmp && make -j4 install)
rm -rf pyHealpix
......@@ -8,3 +8,4 @@ coverage
git+https://gitlab.mpcdf.mpg.de/ift/mpi_dummy.git
git+https://gitlab.mpcdf.mpg.de/ift/keepers.git
git+https://gitlab.mpcdf.mpg.de/ift/D2O.git
git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
import numpy as np
from nifty import (VL_BFGS, DiagonalOperator, FFTOperator, Field,
LinearOperator, PowerSpace, RelaxedNewton, RGSpace,
SteepestDescent, create_power_operator, exp, log, sqrt)
import nifty as ift
from nifty.library.critical_filter import CriticalPowerEnergy
from nifty.library.wiener_filter import WienerFilterEnergy
......@@ -12,23 +10,24 @@ from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(42)
np.random.seed(44)
def plot_parameters(m, t, p, p_d):
def plot_parameters(m, t, p, p_sig,p_d):
x = log(t.domain[0].kindex)
x = np.log(t.domain[0].kindex)
m = fft.adjoint_times(m)
m = m.val.get_full_data().real
t = t.val.get_full_data().real
p = p.val.get_full_data().real
pd_sig = p_sig.val.get_full_data()
p_d = p_d.val.get_full_data().real
pl.plot([go.Heatmap(z=m)], filename='map.html', auto_open=False)
pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p),
go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
go.Scatter(x=x, y=p_d),go.Scatter(x=x, y=pd_sig)], filename="t.html", auto_open=False)
class AdjointFFTResponse(LinearOperator):
class AdjointFFTResponse(ift.LinearOperator):
def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces)
self._domain = FFT.target
......@@ -60,50 +59,55 @@ if __name__ == "__main__":
distribution_strategy = 'not'
# Set up position space
s_space = RGSpace([128, 128])
# s_space = HPSpace(32)
dist = 1/128. *0.1
s_space = ift.RGSpace([128, 128], distances=[dist,dist])
# s_space = ift.HPSpace(32)
# Define harmonic transformation and associated harmonic space
fft = FFTOperator(s_space)
fft = ift.FFTOperator(s_space)
h_space = fft.target[0]
# Set up power space
p_space = PowerSpace(h_space, logarithmic=True,
distribution_strategy=distribution_strategy)
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True),
distribution_strategy=distribution_strategy)
# Choose the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
# p_spec = (lambda k: 1)
S = ift.create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
# Draw a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sp = ift.Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
# Choose the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument = ift.SmoothingOperator(s_space, sigma=0.01)
Instrument = ift.DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0
# Add a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
std=sqrt(noise),
mean=0)
noise = 10.
ndiag = ift.Field(s_space, noise).weight(1)
N = ift.DiagonalOperator(s_space, ndiag)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=np.sqrt(noise),
mean=0)
# Create mock data
d = R(sh) + n
# The information source
j = R.adjoint_times(N.inverse_times(d))
realized_power = log(sh.power_analyze(binbounds=p_space.binbounds))
data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
realized_power = ift.log(sh.power_analyze(binbounds=p_space.binbounds))
data_power = ift.log(fft(d).power_analyze(binbounds=p_space.binbounds))
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False)
......@@ -113,28 +117,25 @@ if __name__ == "__main__":
x = a_energy.value
print(x, iteration)
minimizer1 = RelaxedNewton(convergence_tolerance=1e-4,
convergence_level=1,
iteration_limit=5,
callback=convergence_measure)
minimizer2 = VL_BFGS(convergence_tolerance=1e-10,
convergence_level=1,
iteration_limit=30,
callback=convergence_measure,
max_history_length=20)
minimizer3 = SteepestDescent(convergence_tolerance=1e-4,
iteration_limit=100,
callback=convergence_measure)
IC1 = ift.GradientNormController(iteration_limit=5,
tol_abs_gradnorm=0.1)
minimizer1 = ift.RelaxedNewton(IC1)
IC2 = ift.GradientNormController(iteration_limit=30,
tol_abs_gradnorm=0.1)
minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
IC3 = ift.GradientNormController(iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer3 = ift.SteepestDescent(IC3)
# Set starting position
flat_power = Field(p_space, val=1e-8)
flat_power = ift.Field(p_space, val=1e-8)
m0 = flat_power.power_synthesize(real_signal=True)
t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
# t0 = ift.Field(p_space, val=np.log(1./(1+p_space.kindex)**2))
t0 = ift.Field(p_space, val=-5)
for i in range(500):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
S0 = ift.create_power_operator(h_space, power_spectrum=ift.exp(t0),
distribution_strategy=distribution_strategy)
# Initialize non-linear Wiener Filter energy
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
......@@ -143,14 +144,15 @@ if __name__ == "__main__":
m0 = D0.inverse_times(j)
# Initialize power energy with updated parameters
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0,
smoothness_prior=10., samples=3)
smoothness_prior=1., samples=5)
(power_energy, convergence) = minimizer2(power_energy)
(power_energy, convergence) = minimizer1(power_energy)
# Set new power spectrum
t0.val = power_energy.position.val.real
# Plot current estimate
print(i)
if i % 5 == 0:
plot_parameters(m0, t0, log(sp), data_power)
if i % 1 == 0:
plot_parameters(sh, t0, ift.log(sp), ift.log(sh.power_analyze(binbounds=p_space.binbounds)),data_power)
print ift.log(sh.power_analyze(binbounds=p_space.binbounds)).val - t0.val
# -*- coding: utf-8 -*-
from nifty import *
import nifty as ift
import numpy as np
if __name__ == "__main__":
nifty_configuration['default_distribution_strategy'] = 'fftw'
ift.nifty_configuration['default_distribution_strategy'] = 'not'
# Setting up parameters |\label{code:wf_parameters}|
correlation_length = 1. # Typical distance over which the field is correlated
......@@ -20,53 +21,53 @@ if __name__ == "__main__":
L = 2. # Total side-length of the domain
N_pixels = 128 # Grid resolution (pixels per axis)
#signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space = HPSpace(16)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
fft = FFTOperator(harmonic_space, target=signal_space, target_dtype=np.float)
power_space = PowerSpace(harmonic_space)
signal_space = ift.HPSpace(16)
harmonic_space = ift.FFTOperator.get_default_codomain(signal_space)
fft = ift.FFTOperator(harmonic_space, target=signal_space, target_dtype=np.float)
power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}|
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = Field(power_space, val=power_spectrum)
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.Field(power_space, val=power_spectrum)
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response
mask = Field(signal_space, val=1.)
mask = ift.Field(signal_space, val=1.)
N10 = int(N_pixels/10)
#mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
noise = Field.from_random(domain=data_domain, random_type='normal',
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
data = R(ift.exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
# Wiener filter
m0 = Field(harmonic_space, val=0.j)
energy = library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S)
m0 = ift.Field(harmonic_space, val=0.j)
energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S)
minimizer1 = VL_BFGS(convergence_tolerance=1e-5,
iteration_limit=3000,
#callback=convergence_measure,
max_history_length=20)
minimizer2 = RelaxedNewton(convergence_tolerance=1e-5,
iteration_limit=10,
#callback=convergence_measure
)
minimizer3 = SteepestDescent(convergence_tolerance=1e-5, iteration_limit=1000)
IC1 = ift.GradientNormController(iteration_limit=30,
tol_abs_gradnorm=0.1)
minimizer1 = ift.VL_BFGS(IC1, max_history_length=20)
IC2 = ift.GradientNormController(iteration_limit=5,
tol_abs_gradnorm=0.1)
minimizer2 = ift.RelaxedNewton(IC2)
IC3 = ift.GradientNormController(iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer3 = ift.SteepestDescent(IC3)
# me1 = minimizer1(energy)
# me2 = minimizer2(energy)
me2 = minimizer2(energy)
# me3 = minimizer3(energy)
# m1 = fft(me1[0].position)
# m2 = fft(me2[0].position)
m2 = fft(me2[0].position)
# m3 = fft(me3[0].position)
#
......@@ -81,19 +82,19 @@ if __name__ == "__main__":
#Plotting #|\label{code:wf_plotting}|
#plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter = plotting.HealpixPlotter(color_map=plotting.colormaps.PlankCmap())
plotter = ift.plotting.HealpixPlotter(color_map=ift.plotting.colormaps.PlankCmap())
plotter.figure.xaxis = plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = plotting.Axis(label='Pixel Index')
plotter.figure.xaxis = ift.plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = ift.plotting.Axis(label='Pixel Index')
plotter.plot.zmax = 5; plotter.plot.zmin = -5
#plotter.plot.zmax = 5; plotter.plot.zmin = -5
## plotter(variance, path = 'variance.html')
# #plotter.plot.zmin = exp(mock_signal.min());
# plotter(mock_signal.real, path='mock_signal.html')
# plotter(Field(signal_space, val=np.log(data.val.get_full_data().real).reshape(signal_space.shape)),
# path = 'log_of_data.html')
#plotter.plot.zmin = np.exp(mock_signal.min());
plotter(mock_signal.real, path='mock_signal.html')
plotter(ift.Field(signal_space, val=np.log(data.val.get_full_data().real).reshape(signal_space.shape)),
path = 'log_of_data.html')
#
# plotter(m1.real, path='m_LBFGS.html')
# plotter(m2.real, path='m_Newton.html')
plotter(m2.real, path='m_Newton.html')
# plotter(m3.real, path='m_SteepestDescent.html')
#
......@@ -7,7 +7,7 @@ from nifty import plotting
from keepers import Repository
if __name__ == "__main__":
ift.nifty_configuration['default_distribution_strategy'] = 'fftw'
ift.nifty_configuration['default_distribution_strategy'] = 'not'
signal_to_noise = 1.5 # The signal to noise ratio
......@@ -93,8 +93,8 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=(0, 1, 0, 1))
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise,
bare=True)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0)
......@@ -110,7 +110,7 @@ if __name__ == "__main__":
# Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober): pass
proby = Proby((signal_space_1, signal_space_2), probe_count=100)
proby = Proby((signal_space_1, signal_space_2), probe_count=1)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
# sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1))
......@@ -130,7 +130,7 @@ if __name__ == "__main__":
plotter.plot.zmin = 0.
plotter.plot.zmax = 3.
sm = ift.SmoothingOperator.make(plot_space, sigma=0.03)
sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03)
plotter(ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), path='uncertainty.html')
plotter.plot.zmin = np.real(mock_signal.min());
......
......@@ -7,7 +7,7 @@ from keepers import Repository
if __name__ == "__main__":
ift.nifty_configuration['default_distribution_strategy'] = 'fftw'
ift.nifty_configuration['default_distribution_strategy'] = 'not'
# Setting up parameters |\label{code:wf_parameters}|
correlation_length_scale = 1. # Typical distance over which the field is correlated
......@@ -41,7 +41,8 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise #|\label{code:wf_mock_data}|
......@@ -54,10 +55,10 @@ if __name__ == "__main__":
# Probing the uncertainty |\label{code:wf_uncertainty_probing}|
class Proby(ift.DiagonalProberMixin, ift.Prober): pass
proby = Proby(signal_space, probe_count=800)
proby = Proby(signal_space, probe_count=1)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) #|\label{code:wf_variance_fft_wrap}|
sm = ift.SmoothingOperator.make(signal_space, sigma=0.03)
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03)
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) #|\label{code:wf_variance_weighting}|
repo = Repository('repo_800.h5')
......
......@@ -8,7 +8,7 @@ from nifty.library import WienerFilterCurvature
if __name__ == "__main__":
nifty_configuration['default_distribution_strategy'] = 'fftw'
nifty_configuration['default_distribution_strategy'] = 'not'
nifty_configuration['harmonic_rg_base'] = 'real'
# Setting up variable parameters
......@@ -33,7 +33,7 @@ if __name__ == "__main__":
# Total side-length of the domain
L = 2.
# Grid resolution (pixels per axis)
N_pixels = 512
N_pixels = 128
signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
......@@ -54,9 +54,8 @@ if __name__ == "__main__":
data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
N = DiagonalOperator(data_domain,
diagonal=mock_signal.var()/signal_to_noise,
bare=True)
ndiag = Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = DiagonalOperator(data_domain, ndiag)
noise = Field.from_random(domain=data_domain,
random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
......
import d2o
from nifty import *
import plotly.offline as pl
......@@ -8,7 +11,7 @@ from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(42)
# d2o.random.seed(42)
class AdjointFFTResponse(LinearOperator):
......@@ -40,7 +43,8 @@ class AdjointFFTResponse(LinearOperator):
if __name__ == "__main__":
distribution_strategy = 'not'
nifty_configuration['default_distribution_strategy'] = 'not'
nifty_configuration['harmonic_rg_base'] = 'real'
# Set up position space
s_space = RGSpace([128, 128])
......@@ -51,30 +55,29 @@ if __name__ == "__main__":
h_space = fft.target[0]
# Setting up power space
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
p_space = PowerSpace(h_space)
# Choosing the prior correlation structure and defining
# correlation operator
p_spec = (lambda k: (42 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
S = create_power_operator(h_space, power_spectrum=p_spec)
# Drawing a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sp = Field(p_space, val=p_spec)
sh = sp.power_synthesize(real_signal=True)
ss = fft.adjoint_times(sh)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.05)
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
Instrument._diagonal.val[64:80, 32:80] = 0.
# Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1.
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
ndiag = Field(s_space, val=ss.var()/signal_to_noise).weight()
N = DiagonalOperator(s_space, diagonal=ndiag)
n = Field.from_random(domain=s_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
......@@ -94,9 +97,17 @@ if __name__ == "__main__":
# iteration_limit=50,
# callback=convergence_measure)
minimizer = RelaxedNewton(convergence_tolerance=0,
iteration_limit=1,
callback=convergence_measure)
controller = GradientNormController(iteration_limit=50,
callback=convergence_measure)
minimizer = VL_BFGS(controller=controller)
controller = GradientNormController(iteration_limit=1,
callback=convergence_measure)
minimizer = RelaxedNewton(controller=controller)
#
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
......@@ -104,9 +115,6 @@ if __name__ == "__main__":
# max_history_length=3)
#
inverter = ConjugateGradient(convergence_level=3,
convergence_tolerance=1e-5,
preconditioner=None)
# Setting starting position
m0 = Field(h_space, val=.0)
......@@ -116,14 +124,25 @@ if __name__ == "__main__":
# Solving the problem analytically
m0 = D0.inverse_times(j)
energy, convergence = minimizer(energy)
m = energy.position
D = energy.curvature
m = minimizer(energy)[0].position
plotter = plotting.RG2DPlotter()
plotter(ss, path='signal.html')
plotter(fft.inverse_times(m), path='m.html')
sample_variance = Field(sh.domain, val=0.)
sample_mean = Field(sh.domain, val=0.)
sample_variance = Field(s_space, val=0.)
sample_mean = Field(s_space, val=0.)
# sampling the uncertainty map
n_samples = 10
n_samples = 50
for i in range(n_samples):
sample = fft(sugar.generate_posterior_sample(0., D0))
sample = fft.adjoint_times(sugar.generate_posterior_sample(m, D))
sample_variance += sample**2
sample_mean += sample
variance = (sample_variance - sample_mean**2)/n_samples
sample_mean /= n_samples
sample_variance /= n_samples
variance = (sample_variance - sample_mean**2)
Required packages:
sphinx
numpydoc
sphinx_rtd_theme
#Two step creation of webpages:
sphinx-apidoc -l -e -d 3 -o docs/source/mod/ nifty/ nifty/plotting/ nifty/spaces/power_space/power_indices.py nifty/spaces/power_space/power_index_factory.py nifty/config/ nifty/basic_arithmetics.py nifty/nifty_meta.py nifty/random.py nifty/version.py nifty/field_types/ nifty/operators/fft_operator/transformations/rg_transforms.py
#creates all .rst files neccesary for ModuleIndex excluding helper modules
sphinx-build -b html docs/source/ docs/build/
#generates html file and build directory
Two step creation of webpages:
sphinx-apidoc -l -e -d 3 -o sphinx/source/mod/ nifty/ nifty/plotting/ nifty/spaces/power_space/power_indices.py nifty/spaces/power_space/power_index_factory.py nifty/config/ nifty/basic_arithmetics.py nifty/nifty_meta.py nifty/random.py nifty/version.py nifty/field_types/ nifty/operators/fft_operator/transformations/rg_transforms.py
creates all .rst files neccesary for ModuleIndex excluding helper modules
sphinx-build -b html sphinx/source/ sphinx/build/
generates html filel amd build directory
......@@ -30,6 +30,8 @@ from .config import dependency_injector,\
nifty_configuration,\
d2o_configuration
logger.logger.setLevel(nifty_configuration['loglevel'])
from d2o import distributed_data_object, d2o_librarian
from .field import Field
......@@ -44,6 +46,8 @@ from .field_types import *
from .energies import *
from .memoization import memo
from .minimization import *
from .spaces import *
......
......@@ -27,116 +27,114 @@ __all__ = ['cos', 'sin', 'cosh', 'sinh', 'tan', 'tanh', 'arccos', 'arcsin',
'conjugate', 'clipped_exp', 'limited_exp', 'limited_exp_deriv']
def _math_helper(x, function):
if isinstance(x, Field):
result_val = x.val.apply_scalar_function(function)
result = x.copy_empty(dtype=result_val.dtype)
result.val = result_val
elif isinstance(x, distributed_data_object):
result = x.apply_scalar_function(function, inplace=False)
def _math_helper(x, function, out):
if not isinstance(x, Field):
raise TypeError("This function only accepts Field objects.")
if out is not None:
if not isinstance(out, Field) or x.domain!=out.domain:
raise ValueError("Bad 'out' argument")
function(x.val, out=out.val)
return out
else:
result = function(np.asarray(x))
return Field(domain=x.domain, val=function(x.val))
def cos(x, out=None):
return _math_helper(x, np.cos, out)
def sin(x, out=None):
return _math_helper(x, np.sin, out)
return result
def cosh(x, out=None):
return _math_helper(x, np.cosh, out)
def cos(x):
return _math_helper(x, np.cos)
def sinh(x, out=None):
return _math_helper(x, np.sinh, out)
def sin(x):
return _math_helper(x, np.sin)
def tan(x, out=None):
return _math_helper(x, np.tan, out)
def cosh(x):
return _math_helper(x, np.cosh)
def tanh(x, out=None):
return _math_helper(x, np.tanh, out)
def sinh(x):
return _math_helper(x, np.sinh)
def arccos(x, out=None):
return _math_helper(x, np.arccos, out)
def tan(x):
return _math_helper(x, np.tan)
def arcsin(x, out=None):
return _math_helper(x, np.arcsin, out)
def tanh(x):
return _math_helper(x, np.tanh)
def arccosh(x, out=None):
return _math_helper(x, np.arccosh, out)
def arccos(x):
return _math_helper(x, np.arccos)
def arcsinh(x, out=None):
return _math_helper(x, np.arcsinh, out)
def arcsin(x):
return _math_helper(x, np.arcsin)
def arctan(x, out=None):
return _math_helper(x, np.arctan, out)
def arccosh(x):
return _math_helper(x, np.arccosh)
def arctanh(x, out=None):
return _math_helper(x, np.arctanh, out)
def arcsinh(x):
return _math_helper(x, np.arcsinh)
def sqrt(x, out=None):
return _math_helper(x, np.sqrt, out)
def arctan(x):
return _math_helper(x, np.arctan)
def exp(x, out=None):
return _math_helper(x, np.exp, out)
def arctanh(x):
return _math_helper(x, np.arctanh)
def log(x, out=None):
return _math_helper(x, np.log, out)
def sqrt(x):
return _math_helper(x, np.sqrt)
def conjugate(x, out=None):
return _math_helper(x, np.conjugate, out)
def exp(x):
return _math_helper(x, np.exp)
def conj(x, out=None):
return _math_helper(x, np.conj, out)
def clipped_exp(x):
return _math_helper(x, lambda z: np.exp(np.minimum(200, z)))
def clipped_exp(x, out=None):
return _math_helper(x, lambda z: np.exp(np.minimum(200, z)), out)
def limited_exp(x, out=None):
return _math_helper(x, _limited_exp_helper, out)
def limited_exp(x):
return _math_helper(x, _limited_exp_helper)
def _limited_exp_helper(x):
thr = 200.
mask = x>thr
mask = x > thr
if np.count_nonzero(mask) == 0:
return np.exp(x)
result = ((1.-thr) + x)*np.exp(thr)
result[~mask] = np.exp(x[~mask])
return result
def limited_exp_deriv(x):
return _math_helper(x, _limited_exp_deriv_helper)
def limited_exp_deriv(x, out=None):
return _math_helper(x, _limited_exp_deriv_helper, out)
def _limited_exp_deriv_helper(x):
thr = 200.
mask = x>thr
mask = x > thr
if np.count_nonzero(mask) == 0:
return np.exp(x)
result = np.empty_like(x)
result[mask] = np.exp(thr)
result[~mask] = np.exp(x[~mask])
return result
def log(x, base=None):
result = _math_helper(x, np.log)
if base is not None:
result = result/log(base)
return result
def conjugate(x):
return _math_helper(x, np.conjugate)
def conj(x):
return _math_helper(x, np.conjugate)
......@@ -76,12 +76,26 @@ variable_harmonic_rg_base = keepers.Variable(
lambda z: z in ['real', 'complex'],
genus='str')
variable_threads = keepers.Variable(
'threads',
[1],
lambda z: np.int(abs(z)) == z,
genus='int')
variable_loglevel = keepers.Variable(
'loglevel',
[10],
lambda z: np.int(z) == z and 0 <= z <= 50,
genus='int')
nifty_configuration = keepers.get_Configuration(
name='NIFTy',
variables=[variable_fft_module,
variable_default_field_dtype,
variable_default_distribution_strategy,
variable_harmonic_rg_base],
variable_harmonic_rg_base,
variable_threads,
variable_loglevel],
file_name='NIFTy.conf',
search_paths=[os.path.expanduser('~') + "/.config/nifty/",
os.path.expanduser('~') + "/.config/",
......
......@@ -25,8 +25,8 @@ from keepers import Loggable,\
from future.utils import with_metaclass
class DomainObject(with_metaclass(NiftyMeta, type('NewBase',
(Versionable, Loggable, object), {}))):
class DomainObject(with_metaclass(
NiftyMeta, type('NewBase', (Versionable, Loggable, object), {}))):
"""The abstract class that can be used as a domain for a field.
This holds all the information and functionality a field needs to know
......@@ -43,8 +43,7 @@ class DomainObject(with_metaclass(NiftyMeta, type('NewBase',
"""
def __init__(self):
# _global_id is used in the Versioning module from keepers
self._ignore_for_hash = ['_global_id']
self._needed_for_hash = []
@abc.abstractmethod
def __repr__(self):
......@@ -53,11 +52,8 @@ class DomainObject(with_metaclass(NiftyMeta, type('NewBase',
def __hash__(self):
# Extract the identifying parts from the vars(self) dict.
result_hash = 0
for key in sorted(vars(self).keys()):
item = vars(self)[key]
if key in self._ignore_for_hash or key == '_ignore_for_hash':
continue
result_hash ^= item.__hash__() ^ int(hash(key)//117)
for key in self._needed_for_hash:
result_hash ^= hash(vars(self)[key])
return result_hash
def __eq__(self, x):
......@@ -74,18 +70,14 @@ class DomainObject(with_metaclass(NiftyMeta, type('NewBase',
True if `self` and x describe the same manifold.
"""
if isinstance(x, type(self)):
for key in list(vars(self).keys()):
item1 = vars(self)[key]
if key in self._ignore_for_hash or key == '_ignore_for_hash':
continue
item2 = vars(x)[key]
if item1 != item2:
return False
if self is x: # shortcut for simple case
return True
else:
if not isinstance(x, type(self)):
return False
for key in self._needed_for_hash:
if vars(self)[key] != vars(x)[key]:
return False
return True
def __ne__(self, x):
return not self.__eq__(x)
......
......@@ -17,5 +17,5 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from .energy import Energy
from .quadratic_energy import QuadraticEnergy
from .line_energy import LineEnergy
from .memoization import memo
......@@ -17,12 +17,14 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..nifty_meta import NiftyMeta
from ..memoization import memo
from keepers import Loggable
from future.utils import with_metaclass
class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {}))):
class Energy(with_metaclass(NiftyMeta,
type('NewBase', (Loggable, object), {}))):
""" Provides the functional used by minimization schemes.
The Energy object is an implementation of a scalar function including its
......@@ -41,7 +43,7 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
value : np.float
The value of the energy functional at given `position`.
gradient : Field
The gradient at given `position` in parameter direction.
The gradient at given `position`.
curvature : LinearOperator, callable
A positive semi-definite operator or function describing the curvature
of the potential at the given `position`.
......@@ -64,12 +66,18 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
"""
def __init__(self, position):
def __init__(self, position, gradient=None, curvature=None):
super(Energy, self).__init__()
self._cache = {}
self._position = position.copy()
def at(self, position):
self._cache = {}
if gradient is not None:
self._cache['gradient'] = gradient
if curvature is not None:
self._cache['curvature'] = curvature
def at(self, position, gradient=None, curvature=None):
""" Initializes and returns a new Energy object at the new position.
Parameters
......@@ -83,8 +91,9 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
Energy object at new position.
"""
return self.__class__(position)
return self.__class__(position,
gradient=gradient,
curvature=curvature)
@property
def position(self):
......@@ -97,6 +106,7 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
return self._position
@property
@memo
def value(self):
"""
The value of the energy functional at given `position`.
......@@ -106,15 +116,37 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
raise NotImplementedError
@property
@memo
def gradient(self):
"""
The gradient at given `position` in parameter direction.
The gradient at given `position`.
"""
raise NotImplementedError
@property
@memo
def gradient_norm(self):
"""
The length of the gradient at given `position`.
"""
return self.gradient.norm()
@property
@memo
def gradient_infnorm(self):
"""
The infinity norm of the gradient at given `position`.
"""
return abs(self.gradient).max()
@property
@memo
def curvature(self):
"""
A positive semi-definite operator or function describing the curvature
......
......@@ -18,8 +18,13 @@
from __future__ import print_function
from keepers import Loggable
from ..nifty_meta import NiftyMeta
from future.utils import with_metaclass
class LineEnergy(object):
class LineEnergy((with_metaclass(NiftyMeta,
type('NewBase', (Loggable, object), {})))):
""" Evaluates an underlying Energy along a certain line direction.
Given an Energy class and a line direction, its position is parametrized by
......@@ -74,7 +79,7 @@ class LineEnergy(object):
self._line_position = float(line_position)
self._line_direction = line_direction
if self._line_position==float(offset):
if self._line_position == float(offset):
self.energy = energy
else:
pos = energy.position \
......@@ -116,6 +121,6 @@ class LineEnergy(object):
def directional_derivative(self):
res = self.energy.gradient.vdot(self.line_direction)
if abs(res.imag) / max(abs(res.real), 1.) > 1e-12:
print ("directional derivative has non-negligible "
"imaginary part:", res)
self.logger.warn("directional derivative has non-negligible "
"imaginary part:", res)
return res.real
from .energy import Energy
from ..memoization import memo
class QuadraticEnergy(Energy):
"""The Energy for a quadratic form.
The most important aspect of this energy is that its curvature must be
position-independent.
"""
def __init__(self, position, A, b, gradient=None, curvature=None):
super(QuadraticEnergy, self).__init__(position=position,
gradient=gradient,
curvature=curvature)
self._A = A
self._b = b
if gradient is not None:
self._Ax = gradient + self._b
else:
self._Ax = self._A(self.position)
def at(self, position, gradient=None, curvature=None):
return self.__class__(position=position, A=self._A, b=self._b,
gradient=gradient, curvature=curvature)
@property
@memo
def value(self):
return 0.5*self.position.vdot(self._Ax) - self._b.vdot(self.position)
@property
@memo
def gradient(self):
return self._Ax - self._b
@property
def curvature(self):
return self._A
......@@ -676,7 +676,7 @@ class Field(Loggable, Versionable, object):
result_list[0].val.get_axes_local_distribution_strategy(
result_list[0].domain_axes[power_space_index])
if pindex.distribution_strategy is not local_distribution_strategy:
if pindex.distribution_strategy != local_distribution_strategy:
raise AttributeError(
"The distribution_strategy of pindex does not fit the "
"slice_local distribution strategy of the synthesized field.")
......@@ -814,7 +814,7 @@ class Field(Loggable, Versionable, object):
try:
return int(reduce(lambda x, y: x * y, dim_tuple))
except TypeError:
return 0
return 1
@property
def dof(self):
......@@ -842,7 +842,7 @@ class Field(Loggable, Versionable, object):
@property
def real(self):
""" The real part of the field (data is not copied).
""" The real part of the field.
"""
real_part = self.val.real
result = self.copy_empty(dtype=real_part.dtype)
......@@ -851,11 +851,11 @@ class Field(Loggable, Versionable, object):
@property
def imag(self):
""" The imaginary part of the field (data is not copied).
""" The imaginary part of the field.
"""
real_part = self.val.imag
result = self.copy_empty(dtype=real_part.dtype)
result.set_val(new_val=real_part, copy=False)
imag_part = self.val.imag
result = self.copy_empty(dtype=imag_part.dtype)
result.set_val(new_val=imag_part, copy=False)
return result
# ---Special unary/binary operations---
......@@ -1072,7 +1072,7 @@ class Field(Loggable, Versionable, object):
new_field.set_val(new_val=new_val, copy=False)
return new_field
def vdot(self, x=None, spaces=None, bare=False):
def vdot(self, x=None, spaces=None):
""" Computes the volume-factor-aware dot product of 'self' with x.
Parameters
......@@ -1081,11 +1081,8 @@ class Field(Loggable, Versionable, object):
The domain of x must contain `self.domain`
spaces : tuple of ints
If the domain of `self` and `x` are not the same, `spaces` specfies
the mapping.
bare : boolean
If true, no volume factors will be included in the computation.
If the domain of `self` and `x` are not the same, `spaces`
specifies the mapping.
Returns
-------
......@@ -1097,10 +1094,7 @@ class Field(Loggable, Versionable, object):
"the NIFTy field class")
# Compute the dot respecting the fact of discrete/continuous spaces
if bare:
y = self
else:
y = self.weight(power=1)
y = self.weight(power=1)
if spaces is None:
x_val = x.get_val(copy=False)
......@@ -1590,7 +1584,7 @@ class Field(Loggable, Versionable, object):
new_field.dtype = np.dtype(hdf5_group.attrs['dtype'])
new_field.distribution_strategy =\
hdf5_group.attrs['distribution_strategy']
str(hdf5_group.attrs['distribution_strategy'])
return new_field
......
......@@ -24,6 +24,8 @@ class FieldArray(FieldType):
def __init__(self, shape):
super(FieldArray, self).__init__()
self._needed_for_hash += ["_shape"]
try:
self._shape = tuple([int(i) for i in shape])
except TypeError:
......
from ...operators.endomorphic_operator import EndomorphicOperator
from ...operators.invertible_operator_mixin import InvertibleOperatorMixin
from ...operators.diagonal_operator import DiagonalOperator
from ...minimization import ConjugateGradient
class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
......@@ -21,17 +22,16 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
# ---Overwritten properties and methods---
def __init__(self, theta, T, inverter=None, preconditioner=None, **kwargs):
def __init__(self, theta, T, inverter=None, **kwargs):
self.theta = DiagonalOperator(theta.domain, diagonal=theta)
self.theta = DiagonalOperator(theta.domain, diagonal=theta, copy=False)
self.T = T
if preconditioner is None:
preconditioner = self.theta.inverse_times
if inverter is None:
inverter = ConjugateGradient(
preconditioner=self.theta.inverse_times)
self._domain = self.theta.domain
super(CriticalPowerCurvature, self).__init__(
inverter=inverter,
preconditioner=preconditioner,
**kwargs)
inverter=inverter, **kwargs)
def _times(self, x, spaces):
return self.T(x) + self.theta(x)
......
<
from ...energies.energy import Energy
from ...operators.smoothness_operator import SmoothnessOperator
from ...operators.diagonal_operator import DiagonalOperator
from ...operators.power_projection_operator import PowerProjection
from . import CriticalPowerCurvature
from ...energies.memoization import memo
from ...memoization import memo
from ...minimization import ConjugateGradient
from ...sugar import generate_posterior_sample
from ... import Field, exp
......@@ -54,8 +57,11 @@ class CriticalPowerEnergy(Energy):
# ---Overwritten properties and methods---
def __init__(self, position, m, D=None, alpha=1.0, q=0.,
smoothness_prior=0., logarithmic=True, samples=3, w=None):
super(CriticalPowerEnergy, self).__init__(position=position)
smoothness_prior=0., logarithmic=True, samples=3, w=None,
inverter=None, gradient=None, curvature=None):
super(CriticalPowerEnergy, self).__init__(position=position,
gradient=gradient,
curvature=curvature)
self.m = m
self.D = D
self.samples = samples
......@@ -65,7 +71,21 @@ class CriticalPowerEnergy(Energy):
strength=smoothness_prior,
logarithmic=logarithmic)
self.rho = self.position.domain[0].rho
self.P = PowerProjection(domain=self.m.domain,
target=self.position.domain)
self._w = w if w is not None else None
if inverter is None:
preconditioner = DiagonalOperator(self._theta.domain,
diagonal=self._theta,
copy=False)
inverter = ConjugateGradient(preconditioner=preconditioner)
self._inverter = inverter
self.one = Field(self.position.domain, val=1.)
self.constants = self.one/2. + self.alpha - 1
@property
def inverter(self):
return self._inverter
# ---Mandatory properties and methods---
......@@ -73,28 +93,31 @@ class CriticalPowerEnergy(Energy):
return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
q=self.q, smoothness_prior=self.smoothness_prior,
logarithmic=self.logarithmic,
w=self.w, samples=self.samples)
w=self.w, samples=self.samples,
inverter=self.inverter)
@property
@memo
def value(self):
energy = self._theta.sum()
energy += self.position.vdot(self._rho_prime, bare=True)
energy = self.one.vdot(self._theta)
energy += self.position.vdot(self.constants)
energy += 0.5 * self.position.vdot(self._Tt)
return energy.real
@property
@memo
def gradient(self):
gradient = -self._theta.weight(-1)
gradient += (self._rho_prime).weight(-1)
gradient = -self._theta
gradient += (self.constants)
gradient += self._Tt
gradient.val = gradient.val.real
return gradient
@property
@memo
def curvature(self):
curvature = CriticalPowerCurvature(theta=self._theta.weight(-1),
T=self.T)