Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
ift
NIFTy
Commits
1167f4fa
Commit
1167f4fa
authored
Jan 07, 2019
by
Martin Reinecke
Browse files
Merge remote-tracking branch 'origin/NIFTy_5' into move_dynamic_prior
parents
7beb5b80
05d4cd29
Changes
113
Hide whitespace changes
Inline
Side-by-side
.gitlab-ci.yml
View file @
1167f4fa
...
...
@@ -34,32 +34,24 @@ build_docker_from_cache:
-
docker build -t $CONTAINER_TEST_IMAGE .
-
docker push $CONTAINER_TEST_IMAGE
test
_python2_with_coverage
:
test
:
stage
:
test
variables
:
OMPI_MCA_btl_vader_single_copy_mechanism
:
none
script
:
-
mpiexec -n 2 --bind-to none pytest -q test
-
pytest -q --cov=nifty5 test
-
mpiexec -n 2 --bind-to none pytest
-3
-q test
-
pytest
-3
-q --cov=nifty5 test
-
>
python -m coverage report --omit "*plotting*,*distributed_do*"
python
3
-m coverage report --omit "*plotting*,*distributed_do*"
-
>
python -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
python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
pages
:
stage
:
release
before_script
:
-
ls
script
:
-
python setup.py install --user -f
-
python
3
setup.py install --user -f
-
sh docs/generate.sh
-
mv docs/build/ public/
artifacts
:
...
...
@@ -69,7 +61,6 @@ pages:
-
NIFTy_4
before_script
:
-
python setup.py install --user -f
-
python3 setup.py install --user -f
run_ipynb
:
...
...
@@ -80,7 +71,6 @@ run_ipynb:
run_getting_started_1
:
stage
:
demo_runs
script
:
-
python 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
artifacts
:
...
...
@@ -90,7 +80,6 @@ run_getting_started_1:
run_getting_started_2
:
stage
:
demo_runs
script
:
-
python 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
artifacts
:
...
...
@@ -100,7 +89,6 @@ run_getting_started_2:
run_getting_started_3
:
stage
:
demo_runs
script
:
-
python demos/getting_started_3.py
-
python3 demos/getting_started_3.py
artifacts
:
paths
:
...
...
@@ -109,7 +97,6 @@ run_getting_started_3:
run_bernoulli
:
stage
:
demo_runs
script
:
-
python demos/bernoulli_demo.py
-
python3 demos/bernoulli_demo.py
artifacts
:
paths
:
...
...
@@ -118,7 +105,6 @@ run_bernoulli:
run_curve_fitting
:
stage
:
demo_runs
script
:
-
python demos/polynomial_fit.py
-
python3 demos/polynomial_fit.py
artifacts
:
paths
:
...
...
Dockerfile
View file @
1167f4fa
...
...
@@ -5,27 +5,23 @@ RUN apt-get update && apt-get install -y \
git \
# Packages needed for NIFTy
libfftw3-dev \
python python-pip python-dev python-future python-scipy cython \
python3 python3-pip python3-dev python3-future python3-scipy cython3 \
# Documentation build dependencies
python-sphinx python-sphinx-rtd-theme python-numpydoc \
python
3
-sphinx python
3
-sphinx-rtd-theme python
3
-numpydoc \
# Testing dependencies
python-nose python-coverage python-parameterized python-pytest python-pytest-cov \
python3-nose python3-coverage python3-parameterized python3-pytest python3-pytest-cov \
python3-coverage python3-parameterized python3-pytest python3-pytest-cov \
# Optional NIFTy dependencies
openmpi-bin libopenmpi-dev
python-mpi4py
python3-mpi4py \
openmpi-bin libopenmpi-dev python3-mpi4py \
# Packages needed for NIFTy
&& pip install pyfftw \
&& pip3 install pyfftw \
# 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 \
# Testing dependencies
&& rm -rf /var/lib/apt/lists/*
# Needed for demos to be running
RUN
apt-get update
&&
apt-get
install
-y
python-matplotlib
python3-matplotlib
\
&&
python3
-m
pip
install
--upgrade
pip
&&
python3
-m
pip
install
jupyter
&&
python
-m
pip
install
--upgrade
pip
&&
python
-m
pip
install
jupyter
\
RUN
apt-get update
&&
apt-get
install
-y
python3-matplotlib
\
&&
python3
-m
pip
install
--upgrade
pip
&&
python3
-m
pip
install
jupyter
\
&&
rm
-rf
/var/lib/apt/lists/
*
# Set matplotlib backend
...
...
README.md
View file @
1167f4fa
...
...
@@ -37,7 +37,7 @@ Installation
### 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/
)
-
[
pyFFTW
](
https://pypi.python.org/pypi/pyFFTW
)
...
...
@@ -61,8 +61,8 @@ distributions, the "apt" lines will need slight changes.
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_5
sudo apt-get install git libfftw3-dev python
3
python
3
-pip python
3
-dev
pip
3
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
using a pip-installed
`pyFFTW`
package. Unfortunately, some distributions are
...
...
@@ -71,35 +71,27 @@ with the installed `FFTW3` libraries.)
Plotting support is added via:
pip install --user matplotlib
pip
3
install --user matplotlib
Support for spherical harmonic transforms is added via:
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
pip
3
install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
MPI support is added via:
sudo apt-get install openmpi-bin libopenmpi-dev
pip 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`
pip3 install --user mpi4py
### Running the tests
In oder t
o run the tests
one needs two
additional packages:
T
o 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
following command in the repository root:
nosetests -x --with-coverage --cover-html --cover-package
=nifty5
pytest-3 --cov
=nifty5
test
### First Steps
...
...
@@ -108,7 +100,7 @@ For a quick start, you can browse through the [informal
introduction
](
http://ift.pages.mpcdf.de/NIFTy/code.html
)
or
dive into NIFTy by running one of the demonstrations, e.g.:
python demos/getting_started_1.py
python
3
demos/getting_started_1.py
### Acknowledgement
...
...
demos/bernoulli_demo.py
View file @
1167f4fa
...
...
@@ -11,33 +11,37 @@
# 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-201
8
Max-Planck-Society
# Copyright(C) 2013-201
9
Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
#####################################################################
# 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
nifty5
as
ift
if
__name__
==
'__main__'
:
# FIXME ABOUT THIS CODE
np
.
random
.
seed
(
41
)
# Set up the position space of the signal
#
# # One dimensional regular grid with uniform exposure
# position_space = ift.RGSpace(1024)
# exposure = np.ones(position_space.shape)
# Two-dimensional regular grid with inhomogeneous exposure
position_space
=
ift
.
RGSpace
([
512
,
512
])
# Sphere with uniform exposure
# position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and transform
mode
=
2
if
mode
==
0
:
# One-dimensional regular grid
position_space
=
ift
.
RGSpace
(
1024
)
elif
mode
==
1
:
# Two-dimensional regular grid
position_space
=
ift
.
RGSpace
([
512
,
512
])
else
:
# Sphere
position_space
=
ift
.
HPSpace
(
128
)
# Define harmonic space and transform
harmonic_space
=
position_space
.
get_default_codomain
()
HT
=
ift
.
HarmonicTransformOperator
(
harmonic_space
,
position_space
)
...
...
@@ -45,15 +49,13 @@ if __name__ == '__main__':
# Define power spectrum and amplitudes
def
sqrtpspec
(
k
):
return
1.
/
(
20.
+
k
**
2
)
return
1.
/
(
20.
+
k
**
2
)
A
=
ift
.
create_power_operator
(
harmonic_space
,
sqrtpspec
)
# Set up a sky model
sky
=
ift
.
positive_tanh
(
HT
(
A
))
# Set up a sky model and instrumental response
sky
=
ift
.
sigmoid
(
HT
(
A
))
GR
=
ift
.
GeometryRemover
(
position_space
)
# Set up instrumental response
R
=
GR
# Generate mock data
...
...
@@ -66,8 +68,8 @@ if __name__ == '__main__':
# Compute likelihood and Hamiltonian
position
=
ift
.
from_random
(
'normal'
,
harmonic_space
)
likelihood
=
ift
.
BernoulliEnergy
(
data
)(
p
)
ic_newton
=
ift
.
DeltaEnergyController
(
name
=
'Newton'
,
iteration_limit
=
100
,
tol_rel_deltaE
=
1e-8
)
ic_newton
=
ift
.
DeltaEnergyController
(
name
=
'Newton'
,
iteration_limit
=
100
,
tol_rel_deltaE
=
1e-8
)
minimizer
=
ift
.
NewtonCG
(
ic_newton
)
ic_sampling
=
ift
.
GradientNormController
(
iteration_limit
=
100
)
...
...
@@ -83,5 +85,4 @@ if __name__ == '__main__':
plot
.
add
(
reconstruction
,
title
=
'reconstruction'
)
plot
.
add
(
GR
.
adjoint_times
(
data
),
title
=
'data'
)
plot
.
add
(
sky
(
mock_position
),
title
=
'truth'
)
plot
.
output
(
nx
=
3
,
xsize
=
16
,
ysize
=
5
,
title
=
"results"
,
name
=
"bernoulli.png"
)
plot
.
output
(
nx
=
3
,
xsize
=
16
,
ysize
=
9
,
title
=
"results"
,
name
=
"bernoulli.png"
)
demos/getting_started_1.py
View file @
1167f4fa
...
...
@@ -11,31 +11,40 @@
# 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-201
8
Max-Planck-Society
# Copyright(C) 2013-201
9
Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
###############################################################################
# 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
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
)
for
i
in
range
(
4
):
for
j
in
range
(
4
):
if
(
i
+
j
)
%
2
==
0
:
mask
[
i
*
128
//
4
:(
i
+
1
)
*
128
//
4
,
j
*
128
//
4
:(
j
+
1
)
*
128
//
4
]
=
0
if
(
i
+
j
)
%
2
==
0
:
mask
[
i
*
128
//
4
:(
i
+
1
)
*
128
//
4
,
j
*
128
//
4
:(
j
+
1
)
*
128
//
4
]
=
0
return
mask
def
make_random_mask
():
# Random mask for spherical mode
mask
=
ift
.
from_random
(
'pm1'
,
position_space
)
mask
=
(
mask
+
1
)
/
2
mask
=
(
mask
+
1
)
/
2
return
mask
.
to_global_data
()
def
mask_to_nan
(
mask
,
field
):
# Set masked pixels to nan for plotting
masked_data
=
field
.
local_data
.
copy
()
masked_data
[
mask
.
local_data
==
0
]
=
np
.
nan
return
ift
.
from_local_data
(
field
.
domain
,
masked_data
)
...
...
@@ -43,46 +52,68 @@ def mask_to_nan(mask, field):
if
__name__
==
'__main__'
:
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
if
mode
==
0
:
# One
dimensional regular grid
# One
-
dimensional regular grid
position_space
=
ift
.
RGSpace
([
1024
])
mask
=
np
.
ones
(
position_space
.
shape
)
elif
mode
==
1
:
# Two
dimensional regular grid with che
ss
mask
# Two
-
dimensional regular grid with che
ckerboard
mask
position_space
=
ift
.
RGSpace
([
128
,
128
])
mask
=
make_che
ss
_mask
(
position_space
)
mask
=
make_che
ckerboard
_mask
(
position_space
)
else
:
# Sphere with half of its
location
s randomly masked
# Sphere with half of its
pixel
s randomly masked
position_space
=
ift
.
HPSpace
(
128
)
mask
=
make_random_mask
()
# Specify harmonic space corresponding to signal
harmonic_space
=
position_space
.
get_default_codomain
()
# Harmonic transform from harmonic space to position space
HT
=
ift
.
HarmonicTransformOperator
(
harmonic_space
,
target
=
position_space
)
# Set correlation
structur
e with a power spectrum
and build
#
prior correlation covariance
# Set
prior
correlation
covarianc
e with a power spectrum
leading to
#
homogeneous and isotropic statistics
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
)
# Mapping to (higher dimensional) harmonic space
PD
=
ift
.
PowerDistributor
(
harmonic_space
,
power_space
)
# Apply the mapping
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 is the prior field covariance
# Build instrument response consisting of a discretization, mask
# and harmonic transformaion
# Data is defined on a geometry-free space, thus the geometry is removed
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
.
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
))
# or with @
R
=
GR
@
Mask
@
HT
data_space
=
GR
.
target
# Set the noise covariance
# Set the noise covariance
N
noise
=
5.
N
=
ift
.
ScalingOperator
(
noise
,
data_space
)
...
...
@@ -91,29 +122,30 @@ if __name__ == '__main__':
MOCK_NOISE
=
N
.
draw_sample
()
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
))
D_inv
=
R
.
adjoint
(
N
.
inverse
(
R
))
+
S
.
inverse
# Make it invertible
# Make D_inv invertible (via Conjugate Gradient)
IC
=
ift
.
GradientNormController
(
iteration_limit
=
500
,
tol_abs_gradnorm
=
1e-3
)
D
=
ift
.
InversionEnabler
(
D_inv
,
IC
,
approximation
=
S
.
inverse
).
inverse
# WIENER FILTER
#
Calculate
WIENER FILTER
solution
m
=
D
(
j
)
# P
LOTTING
# P
lotting
rg
=
isinstance
(
position_space
,
ift
.
RGSpace
)
plot
=
ift
.
Plot
()
if
rg
and
len
(
position_space
.
shape
)
==
1
:
plot
.
add
([
HT
(
MOCK_SIGNAL
),
GR
.
adjoint
(
data
),
HT
(
m
)],
label
=
[
'Mock signal'
,
'Data'
,
'Reconstruction'
],
alpha
=
[
1
,
.
3
,
1
])
plot
.
add
(
mask_to_nan
(
mask
,
HT
(
m
-
MOCK_SIGNAL
)),
title
=
'Residuals'
)
plot
.
add
(
[
HT
(
MOCK_SIGNAL
),
GR
.
adjoint
(
data
),
HT
(
m
)],
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"
)
else
:
plot
.
add
(
HT
(
MOCK_SIGNAL
),
title
=
'Mock Signal'
)
plot
.
add
(
mask_to_nan
(
mask
,
(
GR
(
Mask
)).
adjoint
(
data
)),
title
=
'Data'
)
plot
.
add
(
mask_to_nan
(
mask
,
(
GR
(
Mask
)).
adjoint
(
data
)),
title
=
'Data'
)
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"
)
demos/getting_started_2.py
View file @
1167f4fa
...
...
@@ -11,16 +11,23 @@
# 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-201
8
Max-Planck-Society
# Copyright(C) 2013-201
9
Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
###############################################################################
# Log-normal field reconstruction from Poissonian data with inhomogenous
# exposure (in case for 2D mode)
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
###############################################################################
import
nifty5
as
ift
import
numpy
as
np
import
nifty5
as
ift
def
get_2D_exposure
():
def
exposure_2d
():
# Structured exposure for 2D mode
x_shape
,
y_shape
=
position_space
.
shape
exposure
=
np
.
ones
(
position_space
.
shape
)
...
...
@@ -31,72 +38,73 @@ def get_2D_exposure():
exposure
[:,
x_shape
*
4
//
5
:
x_shape
]
*=
.
1
exposure
[:,
x_shape
//
2
:
x_shape
*
3
//
2
]
*=
3.
exposure
=
ift
.
Field
.
from_global_data
(
position_space
,
exposure
)
return
exposure
return
ift
.
Field
.
from_global_data
(
position_space
,
exposure
)
if
__name__
==
'__main__'
:
# FIXME
description of the tutorial
# FIXME
All random seeds to 42
np
.
random
.
seed
(
41
)
# Set up the position space of the signal
#
# # One dimensional regular grid with uniform exposure
# position_space = ift.RGSpace(1024)
# exposure = ift.Field.full(position_space, 1.)
# Two-dimensional regular grid with inhomogeneous exposure
position_space
=
ift
.
RGSpace
([
512
,
512
])
exposure
=
get_2D_exposure
()
# Sphere with uniform exposure
# position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and transform
# Choose space on which the signal field is defined
mode
=
2
if
mode
==
0
:
# One-dimensional regular grid with uniform exposure
position_space
=
ift
.
RGSpace
(
1024
)
exposure
=
ift
.
Field
.
full
(
position_space
,
1.
)
elif
mode
==
1
:
# Two-dimensional regular grid with inhomogeneous exposure
position_space
=
ift
.
RGSpace
([
512
,
512
])
exposure
=
exposure_2d
()
else
:
# Sphere with uniform exposure
position_space
=
ift
.
HPSpace
(
128
)
exposure
=
ift
.
Field
.
full
(
position_space
,
1.
)
# Define harmonic space and harmonic transform
harmonic_space
=
position_space
.
get_default_codomain
()
HT
=
ift
.
HarmonicTransformOperator
(
harmonic_space
,
position_space
)
# Domain on which the field's degrees of freedom are defined
domain
=
ift
.
DomainTuple
.
make
(
harmonic_space
)
position
=
ift
.
from_random
(
'normal'
,
domain
)
# Define
power spectrum and amplitudes
# Define
amplitude (square root of power spectrum)
def
sqrtpspec
(
k
):
return
1.
/
(
20.
+
k
**
2
)
return
1.
/
(
20.
+
k
**
2
)
p_space
=
ift
.
PowerSpace
(
harmonic_space
)
pd
=
ift
.
PowerDistributor
(
harmonic_space
,
p_space
)
a
=
ift
.
PS_field
(
p_space
,
sqrtpspec
)
A
=
pd
(
a
)
#
Set up a
sky model
#
Define
sky model
sky
=
ift
.
exp
(
HT
(
ift
.
makeOp
(
A
)))
M
=
ift
.
DiagonalOperator
(
exposure
)
GR
=
ift
.
GeometryRemover
(
position_space
)
#
Set up
instrumental response
#
Define
instrumental response
R
=
GR
(
M
)
# Generate mock data
# Generate mock data
and define likelihood operator
d_space
=
R
.
target
[
0
]
lamb
=
R
(
sky
)
mock_position
=
ift
.
from_random
(
'normal'
,
domain
)
data
=
lamb
(
mock_position
)
data
=
np
.
random
.
poisson
(
data
.
to_global_data
().
astype
(
np
.
float64
))
data
=
ift
.
Field
.
from_global_data
(
d_space
,
data
)
# Compute likelihood and Hamiltonian
ic_newton
=
ift
.
DeltaEnergyController
(
name
=
'Newton'
,
iteration_limit
=
100
,
tol_rel_deltaE
=
1e-8
)
likelihood
=
ift
.
PoissonianEnergy
(
data
)(
lamb
)
# Settings for minimization
ic_newton
=
ift
.
DeltaEnergyController
(
name
=
'Newton'
,
iteration_limit
=
100
,
tol_rel_deltaE
=
1e-8
)
minimizer
=
ift
.
NewtonCG
(
ic_newton
)
#
M
inimiz
e
the Hamiltonian
#
Compute MAP solution by m
inimiz
ing
the
information
Hamiltonian
H
=
ift
.
Hamiltonian
(
likelihood
)
H
=
ift
.
EnergyAdapter
(
position
,
H
,
want_metric
=
True
)
initial_position
=
ift
.
from_random
(
'normal'
,
domain
)
H
=
ift
.
EnergyAdapter
(
initial_position
,
H
,
want_metric
=
True
)
H
,
convergence
=
minimizer
(
H
)
# Plot
results
# Plot
ting
signal
=
sky
(
mock_position
)
reconst
=
sky
(
H
.
position
)
plot
=
ift
.
Plot
()
...
...
@@ -104,4 +112,4 @@ if __name__ == '__main__':
plot
.
add
(
GR
.
adjoint
(
data
),
title
=
'Data'
)
plot
.
add
(
reconst
,
title
=
'Reconstruction'
)
plot
.
add
(
reconst
-
signal
,
title
=
'Residuals'
)
plot
.
output
(
name
=
'getting_started_2.p
ng
'
,
xsize
=
16
,
ysize
=
16
)
plot
.
output
(
name
=
'getting_started_2.p
df
'
,
xsize
=
16
,
ysize
=
16
)
demos/getting_started_3.py
View file @
1167f4fa
...
...
@@ -11,104 +11,137 @@
# 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-201
8
Max-Planck-Society
# Copyright(C) 2013-201
9
Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
############################################################
# Non-linear tomography
# The data is integrated lines of sight
# Random lines (set mode=0), radial lines (mode=1)
#############################################################
import
nifty5
as
ift
import
numpy
as
np
import
nifty5
as
ift
def
random_los
(
n_los
):
starts
=
list
(
np
.
random
.
uniform
(
0
,
1
,
(
n_los
,
2
)).
T
)
ends
=
list
(
0.5
+
0
*
np
.
random
.
uniform
(
0
,
1
,
(
n_los
,
2
)).
T
)
return
starts
,
ends
def
get_random_LOS
(
n_los
):
def
radial_los
(
n_los
):
starts
=
list
(
np
.
random
.
uniform
(
0
,
1
,
(
n_los
,
2
)).
T
)
ends
=
list
(
np
.
random
.
uniform
(
0
,
1
,
(
n_los
,
2
)).
T
)
return
starts
,
ends
if
__name__
==