Commit c12019fb authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'NIFTy_5' into move_dynamic_prior

parents 7beb5b80 dd23c622
...@@ -34,32 +34,24 @@ build_docker_from_cache: ...@@ -34,32 +34,24 @@ build_docker_from_cache:
- docker build -t $CONTAINER_TEST_IMAGE . - docker build -t $CONTAINER_TEST_IMAGE .
- docker push $CONTAINER_TEST_IMAGE - docker push $CONTAINER_TEST_IMAGE
test_python2_with_coverage: test:
stage: test stage: test
variables: variables:
OMPI_MCA_btl_vader_single_copy_mechanism: none OMPI_MCA_btl_vader_single_copy_mechanism: none
script: script:
- mpiexec -n 2 --bind-to none pytest -q test - mpiexec -n 2 --bind-to none pytest-3 -q test
- pytest -q --cov=nifty5 test - pytest-3 -q --cov=nifty5 test
- > - >
python -m coverage report --omit "*plotting*,*distributed_do*" python3 -m coverage report --omit "*plotting*,*distributed_do*"
- > - >
python -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }' python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
test_python3:
stage: test
variables:
OMPI_MCA_btl_vader_single_copy_mechanism: none
script:
- pytest-3 -q
- mpiexec -n 2 --bind-to none pytest-3 -q
pages: pages:
stage: release stage: release
before_script: before_script:
- ls - ls
script: script:
- python setup.py install --user -f - python3 setup.py install --user -f
- sh docs/generate.sh - sh docs/generate.sh
- mv docs/build/ public/ - mv docs/build/ public/
artifacts: artifacts:
...@@ -69,7 +61,6 @@ pages: ...@@ -69,7 +61,6 @@ pages:
- NIFTy_4 - NIFTy_4
before_script: before_script:
- python setup.py install --user -f
- python3 setup.py install --user -f - python3 setup.py install --user -f
run_ipynb: run_ipynb:
...@@ -80,7 +71,6 @@ run_ipynb: ...@@ -80,7 +71,6 @@ run_ipynb:
run_getting_started_1: run_getting_started_1:
stage: demo_runs stage: demo_runs
script: script:
- python demos/getting_started_1.py
- python3 demos/getting_started_1.py - python3 demos/getting_started_1.py
- mpiexec -n 2 --bind-to none python3 demos/getting_started_1.py 2> /dev/null - mpiexec -n 2 --bind-to none python3 demos/getting_started_1.py 2> /dev/null
artifacts: artifacts:
...@@ -90,7 +80,6 @@ run_getting_started_1: ...@@ -90,7 +80,6 @@ run_getting_started_1:
run_getting_started_2: run_getting_started_2:
stage: demo_runs stage: demo_runs
script: script:
- python demos/getting_started_2.py
- python3 demos/getting_started_2.py - python3 demos/getting_started_2.py
- mpiexec -n 2 --bind-to none python3 demos/getting_started_2.py 2> /dev/null - mpiexec -n 2 --bind-to none python3 demos/getting_started_2.py 2> /dev/null
artifacts: artifacts:
...@@ -100,7 +89,6 @@ run_getting_started_2: ...@@ -100,7 +89,6 @@ run_getting_started_2:
run_getting_started_3: run_getting_started_3:
stage: demo_runs stage: demo_runs
script: script:
- python demos/getting_started_3.py
- python3 demos/getting_started_3.py - python3 demos/getting_started_3.py
artifacts: artifacts:
paths: paths:
...@@ -109,7 +97,6 @@ run_getting_started_3: ...@@ -109,7 +97,6 @@ run_getting_started_3:
run_bernoulli: run_bernoulli:
stage: demo_runs stage: demo_runs
script: script:
- python demos/bernoulli_demo.py
- python3 demos/bernoulli_demo.py - python3 demos/bernoulli_demo.py
artifacts: artifacts:
paths: paths:
...@@ -118,7 +105,6 @@ run_bernoulli: ...@@ -118,7 +105,6 @@ run_bernoulli:
run_curve_fitting: run_curve_fitting:
stage: demo_runs stage: demo_runs
script: script:
- python demos/polynomial_fit.py
- python3 demos/polynomial_fit.py - python3 demos/polynomial_fit.py
artifacts: artifacts:
paths: paths:
......
...@@ -5,27 +5,23 @@ RUN apt-get update && apt-get install -y \ ...@@ -5,27 +5,23 @@ RUN apt-get update && apt-get install -y \
git \ git \
# Packages needed for NIFTy # Packages needed for NIFTy
libfftw3-dev \ libfftw3-dev \
python python-pip python-dev python-future python-scipy cython \
python3 python3-pip python3-dev python3-future python3-scipy cython3 \ python3 python3-pip python3-dev python3-future python3-scipy cython3 \
# Documentation build dependencies # Documentation build dependencies
python-sphinx python-sphinx-rtd-theme python-numpydoc \ python3-sphinx python3-sphinx-rtd-theme python3-numpydoc \
# Testing dependencies # Testing dependencies
python-nose python-coverage python-parameterized python-pytest python-pytest-cov \ python3-coverage python3-parameterized python3-pytest python3-pytest-cov \
python3-nose python3-coverage python3-parameterized python3-pytest python3-pytest-cov \
# Optional NIFTy dependencies # Optional NIFTy dependencies
openmpi-bin libopenmpi-dev python-mpi4py python3-mpi4py \ openmpi-bin libopenmpi-dev python3-mpi4py \
# Packages needed for NIFTy # Packages needed for NIFTy
&& pip install pyfftw \
&& pip3 install pyfftw \ && pip3 install pyfftw \
# Optional NIFTy dependencies # Optional NIFTy dependencies
&& pip install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \ && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
# Testing dependencies # Testing dependencies
&& rm -rf /var/lib/apt/lists/* && rm -rf /var/lib/apt/lists/*
# Needed for demos to be running # Needed for demos to be running
RUN apt-get update && apt-get install -y python-matplotlib python3-matplotlib \ RUN apt-get update && apt-get install -y python3-matplotlib \
&& python3 -m pip install --upgrade pip && python3 -m pip install jupyter && python -m pip install --upgrade pip && python -m pip install jupyter \ && python3 -m pip install --upgrade pip && python3 -m pip install jupyter \
&& rm -rf /var/lib/apt/lists/* && rm -rf /var/lib/apt/lists/*
# Set matplotlib backend # Set matplotlib backend
......
...@@ -37,7 +37,7 @@ Installation ...@@ -37,7 +37,7 @@ Installation
### Requirements ### Requirements
- [Python](https://www.python.org/) (v2.7.x or 3.5.x) - [Python 3](https://www.python.org/) (3.5.x or later)
- [SciPy](https://www.scipy.org/) - [SciPy](https://www.scipy.org/)
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW) - [pyFFTW](https://pypi.python.org/pypi/pyFFTW)
...@@ -61,8 +61,8 @@ distributions, the "apt" lines will need slight changes. ...@@ -61,8 +61,8 @@ distributions, the "apt" lines will need slight changes.
NIFTy5 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 sudo apt-get install git libfftw3-dev python3 python3-pip python3-dev
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_5 pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_5
(Note: If you encounter problems related to `pyFFTW`, make sure that you are (Note: If you encounter problems related to `pyFFTW`, make sure that you are
using a pip-installed `pyFFTW` package. Unfortunately, some distributions are using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
...@@ -71,35 +71,27 @@ with the installed `FFTW3` libraries.) ...@@ -71,35 +71,27 @@ with the installed `FFTW3` libraries.)
Plotting support is added via: Plotting support is added via:
pip install --user matplotlib pip3 install --user matplotlib
Support for spherical harmonic transforms is added via: Support for spherical harmonic transforms is added via:
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
MPI support is added via: MPI support is added via:
sudo apt-get install openmpi-bin libopenmpi-dev sudo apt-get install openmpi-bin libopenmpi-dev
pip install --user mpi4py pip3 install --user mpi4py
### Installation for Python 3
If you want to run NIFTy with Python 3, you need to make the following changes
to the instructions above:
- in all `apt-get` commands, replace `python-*` by `python3-*`
- in all `pip` commands, replace `pip` by `pip3`
### Running the tests ### Running the tests
In oder to run the tests one needs two additional packages: To run the tests, additional packages are required:
pip install --user nose parameterized coverage sudo apt-get install python3-coverage python3-parameterized python3-pytest python3-pytest-cov
Afterwards the tests (including a coverage report) can be run using the Afterwards the tests (including a coverage report) can be run using the
following command in the repository root: following command in the repository root:
nosetests -x --with-coverage --cover-html --cover-package=nifty5 pytest-3 --cov=nifty5 test
### First Steps ### First Steps
...@@ -108,7 +100,7 @@ For a quick start, you can browse through the [informal ...@@ -108,7 +100,7 @@ For a quick start, you can browse through the [informal
introduction](http://ift.pages.mpcdf.de/NIFTy/code.html) or introduction](http://ift.pages.mpcdf.de/NIFTy/code.html) or
dive into NIFTy by running one of the demonstrations, e.g.: dive into NIFTy by running one of the demonstrations, e.g.:
python demos/getting_started_1.py python3 demos/getting_started_1.py
### Acknowledgement ### Acknowledgement
......
...@@ -11,33 +11,37 @@ ...@@ -11,33 +11,37 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2018 Max-Planck-Society # Copyright(C) 2013-2019 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
# and financially supported by the Studienstiftung des deutschen Volkes.
#####################################################################
# Bernoulli reconstruction
# Reconstruct an event probability field with values between 0 and 1
# from the observed events
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
#####################################################################
import nifty5 as ift
import numpy as np import numpy as np
import nifty5 as ift
if __name__ == '__main__': if __name__ == '__main__':
# FIXME ABOUT THIS CODE
np.random.seed(41) np.random.seed(41)
# Set up the position space of the signal # Set up the position space of the signal
# mode = 2
# # One dimensional regular grid with uniform exposure if mode == 0:
# position_space = ift.RGSpace(1024) # One-dimensional regular grid
# exposure = np.ones(position_space.shape) position_space = ift.RGSpace(1024)
elif mode == 1:
# Two-dimensional regular grid with inhomogeneous exposure # Two-dimensional regular grid
position_space = ift.RGSpace([512, 512]) position_space = ift.RGSpace([512, 512])
else:
# Sphere with uniform exposure # Sphere
# position_space = ift.HPSpace(128) position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
# Define harmonic space and transform
# Defining harmonic space and transform
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space) HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
...@@ -45,15 +49,13 @@ if __name__ == '__main__': ...@@ -45,15 +49,13 @@ if __name__ == '__main__':
# Define power spectrum and amplitudes # Define power spectrum and amplitudes
def sqrtpspec(k): def sqrtpspec(k):
return 1. / (20. + k**2) return 1./(20. + k**2)
A = ift.create_power_operator(harmonic_space, sqrtpspec) A = ift.create_power_operator(harmonic_space, sqrtpspec)
# Set up a sky model # Set up a sky operator and instrumental response
sky = ift.positive_tanh(HT(A)) sky = ift.sigmoid(HT(A))
GR = ift.GeometryRemover(position_space) GR = ift.GeometryRemover(position_space)
# Set up instrumental response
R = GR R = GR
# Generate mock data # Generate mock data
...@@ -66,8 +68,8 @@ if __name__ == '__main__': ...@@ -66,8 +68,8 @@ if __name__ == '__main__':
# Compute likelihood and Hamiltonian # Compute likelihood and Hamiltonian
position = ift.from_random('normal', harmonic_space) position = ift.from_random('normal', harmonic_space)
likelihood = ift.BernoulliEnergy(data)(p) likelihood = ift.BernoulliEnergy(data)(p)
ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100, ic_newton = ift.DeltaEnergyController(
tol_rel_deltaE=1e-8) name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
minimizer = ift.NewtonCG(ic_newton) minimizer = ift.NewtonCG(ic_newton)
ic_sampling = ift.GradientNormController(iteration_limit=100) ic_sampling = ift.GradientNormController(iteration_limit=100)
...@@ -83,5 +85,4 @@ if __name__ == '__main__': ...@@ -83,5 +85,4 @@ if __name__ == '__main__':
plot.add(reconstruction, title='reconstruction') plot.add(reconstruction, title='reconstruction')
plot.add(GR.adjoint_times(data), title='data') plot.add(GR.adjoint_times(data), title='data')
plot.add(sky(mock_position), title='truth') plot.add(sky(mock_position), title='truth')
plot.output(nx=3, xsize=16, ysize=5, title="results", plot.output(nx=3, xsize=16, ysize=9, title="results", name="bernoulli.png")
name="bernoulli.png")
...@@ -11,31 +11,40 @@ ...@@ -11,31 +11,40 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2018 Max-Planck-Society # Copyright(C) 2013-2019 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
# and financially supported by the Studienstiftung des deutschen Volkes.
###############################################################################
# Compute a Wiener filter solution with NIFTy
# Shows how measurement gaps are filled in
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
###############################################################################
import nifty5 as ift
import numpy as np import numpy as np
import nifty5 as ift
def make_chess_mask(position_space): def make_checkerboard_mask(position_space):
# Checkerboard mask for 2D mode
mask = np.ones(position_space.shape) mask = np.ones(position_space.shape)
for i in range(4): for i in range(4):
for j in range(4): for j in range(4):
if (i+j) % 2 == 0: if (i + j) % 2 == 0:
mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0 mask[i*128//4:(i + 1)*128//4, j*128//4:(j + 1)*128//4] = 0
return mask return mask
def make_random_mask(): def make_random_mask():
# Random mask for spherical mode
mask = ift.from_random('pm1', position_space) mask = ift.from_random('pm1', position_space)
mask = (mask+1)/2 mask = (mask + 1)/2
return mask.to_global_data() return mask.to_global_data()
def mask_to_nan(mask, field): def mask_to_nan(mask, field):
# Set masked pixels to nan for plotting
masked_data = field.local_data.copy() masked_data = field.local_data.copy()
masked_data[mask.local_data == 0] = np.nan masked_data[mask.local_data == 0] = np.nan
return ift.from_local_data(field.domain, masked_data) return ift.from_local_data(field.domain, masked_data)
...@@ -43,46 +52,68 @@ def mask_to_nan(mask, field): ...@@ -43,46 +52,68 @@ def mask_to_nan(mask, field):
if __name__ == '__main__': if __name__ == '__main__':
np.random.seed(42) np.random.seed(42)
# FIXME description of the tutorial
# Choose problem geometry and masking # Choose space on which the signal field is defined
mode = 1 mode = 1
if mode == 0: if mode == 0:
# One dimensional regular grid # One-dimensional regular grid
position_space = ift.RGSpace([1024]) position_space = ift.RGSpace([1024])
mask = np.ones(position_space.shape) mask = np.ones(position_space.shape)
elif mode == 1: elif mode == 1:
# Two dimensional regular grid with chess mask # Two-dimensional regular grid with checkerboard mask
position_space = ift.RGSpace([128, 128]) position_space = ift.RGSpace([128, 128])
mask = make_chess_mask(position_space) mask = make_checkerboard_mask(position_space)
else: else:
# Sphere with half of its locations randomly masked # Sphere with half of its pixels randomly masked
position_space = ift.HPSpace(128) position_space = ift.HPSpace(128)
mask = make_random_mask() mask = make_random_mask()
# Specify harmonic space corresponding to signal
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
# Harmonic transform from harmonic space to position space
HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space) HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
# Set correlation structure with a power spectrum and build # Set prior correlation covariance with a power spectrum leading to
# prior correlation covariance # homogeneous and isotropic statistics
def power_spectrum(k): def power_spectrum(k):
return 100. / (20.+k**3) return 100./(20. + k**3)
# 1D spectral space on which the power spectrum is defined
power_space = ift.PowerSpace(harmonic_space) power_space = ift.PowerSpace(harmonic_space)
# Mapping to (higher dimensional) harmonic space
PD = ift.PowerDistributor(harmonic_space, power_space) PD = ift.PowerDistributor(harmonic_space, power_space)
# Apply the mapping
prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum)) prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
# Insert the result into the diagonal of an harmonic space operator
S = ift.DiagonalOperator(prior_correlation_structure) S = ift.DiagonalOperator(prior_correlation_structure)
# S is the prior field covariance
# Build instrument response consisting of a discretization, mask # Build instrument response consisting of a discretization, mask
# and harmonic transformaion # and harmonic transformaion
# Data is defined on a geometry-free space, thus the geometry is removed
GR = ift.GeometryRemover(position_space) GR = ift.GeometryRemover(position_space)
# Masking operator to model that parts of the field have not been observed
mask = ift.Field.from_global_data(position_space, mask) mask = ift.Field.from_global_data(position_space, mask)
Mask = ift.DiagonalOperator(mask) Mask = ift.DiagonalOperator(mask)
# The response operator consists of
# - an harmonic transform (to get to image space)
# - the application of the mask
# - the removal of geometric information
# Operators can be composed either with parenthesis
R = GR(Mask(HT)) R = GR(Mask(HT))
# or with @
R = GR @ Mask @ HT
data_space = GR.target data_space = GR.target
# Set the noise covariance # Set the noise covariance N
noise = 5. noise = 5.
N = ift.ScalingOperator(noise, data_space) N = ift.ScalingOperator(noise, data_space)
...@@ -91,29 +122,30 @@ if __name__ == '__main__': ...@@ -91,29 +122,30 @@ if __name__ == '__main__':
MOCK_NOISE = N.draw_sample() MOCK_NOISE = N.draw_sample()
data = R(MOCK_SIGNAL) + MOCK_NOISE data = R(MOCK_SIGNAL) + MOCK_NOISE
# Build propagator D and information source j # Build inverse propagator D and information source j
D_inv = R.adjoint @ N.inverse @ R + S.inverse
j = R.adjoint_times(N.inverse_times(data)) j = R.adjoint_times(N.inverse_times(data))
D_inv = R.adjoint(N.inverse(R)) + S.inverse # Make D_inv invertible (via Conjugate Gradient)
# Make it invertible
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3) IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
# WIENER FILTER # Calculate WIENER FILTER solution
m = D(j) m = D(j)
# PLOTTING # Plotting
rg = isinstance(position_space, ift.RGSpace) rg = isinstance(position_space, ift.RGSpace)
plot = ift.Plot() plot = ift.Plot()
if rg and len(position_space.shape) == 1: if rg and len(position_space.shape) == 1:
plot.add([HT(MOCK_SIGNAL), GR.adjoint(data), HT(m)], plot.add(
label=['Mock signal', 'Data', 'Reconstruction'], [HT(MOCK_SIGNAL), GR.adjoint(data),
alpha=[1, .3, 1]) HT(m)],
plot.add(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals') label=['Mock signal', 'Data', 'Reconstruction'],
alpha=[1, .3, 1])
plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
plot.output(nx=2, ny=1, xsize=10, ysize=4, title="getting_started_1") plot.output(nx=2, ny=1, xsize=10, ysize=4, title="getting_started_1")
else: else:
plot.add(HT(MOCK_SIGNAL), title='Mock Signal') plot.add(HT(MOCK_SIGNAL), title='Mock Signal')
plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)), plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)), title='Data')
title='Data')
plot.add(HT(m), title='Reconstruction') plot.add(HT(m), title='Reconstruction')
plot.add(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals') plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
plot.output(nx=2, ny=2, xsize=10, ysize=10, title="getting_started_1") plot.output(nx=2, ny=2, xsize=10, ysize=10, title="getting_started_1")
...@@ -11,16 +11,23 @@ ...@@ -11,16 +11,23 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2018 Max-Planck-Society # Copyright(C) 2013-2019 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
# and financially supported by the Studienstiftung des deutschen Volkes.