Commit 4a9980ba authored by Philipp Frank's avatar Philipp Frank
Browse files

Merge branch 'move_dynamic_prior' of https://gitlab.mpcdf.mpg.de/ift/nifty-dev...

Merge branch 'move_dynamic_prior' of https://gitlab.mpcdf.mpg.de/ift/nifty-dev into move_dynamic_prior
parents b0c150f9 5ceea44c
......@@ -34,18 +34,22 @@ build_docker_from_cache:
- docker build -t $CONTAINER_TEST_IMAGE .
- docker push $CONTAINER_TEST_IMAGE
test:
test_serial:
stage: test
variables:
OMPI_MCA_btl_vader_single_copy_mechanism: none
script:
- mpiexec -n 2 --bind-to none pytest-3 -q test
- pytest-3 -q --cov=nifty5 test
- >
python3 -m coverage report --omit "*plotting*,*distributed_do*"
- >
python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
test_mpi:
stage: test
variables:
OMPI_MCA_btl_vader_single_copy_mechanism: none
script:
- mpiexec -n 2 --bind-to none pytest-3 -q test
pages:
stage: release
before_script:
......
......@@ -9,7 +9,7 @@ RUN apt-get update && apt-get install -y \
# Documentation build dependencies
python3-sphinx python3-sphinx-rtd-theme python3-numpydoc \
# Testing dependencies
python3-coverage python3-parameterized python3-pytest python3-pytest-cov \
python3-coverage python3-pytest python3-pytest-cov \
# Optional NIFTy dependencies
openmpi-bin libopenmpi-dev python3-mpi4py \
# Packages needed for NIFTy
......
......@@ -28,14 +28,14 @@ from ..field import Field
from ..operators.linear_operator import LinearOperator
def _gaussian_error_function(x):
return 0.5/erfc(x*np.sqrt(2.))
def _gaussian_sf(x):
return 0.5*erfc(x/np.sqrt(2.))
def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
ndim = start.shape[0]
nlos = start.shape[1]
inc = np.full(len(shp), 1)
inc = np.full(len(shp), 1, dtype=np.int64)
for i in range(-2, -len(shp)-1, -1):
inc[i] = inc[i+1]*shp[i+1]
......@@ -59,7 +59,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
dmin += 1e-7
dmax -= 1e-7
if dmin >= dmax: # no intersection
out[i] = (np.full(0, 0), np.full(0, 0.))
out[i] = (np.full(0, 0, dtype=np.int64), np.full(0, 0.))
continue
# determine coordinates of first cell crossing
c_first = np.ceil(start[:, i]+dir*dmin)
......@@ -75,7 +75,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
tmp = np.arange(start=c_first[j], stop=dmax,
step=abs(1./dir[j]))
cdist = np.append(cdist, tmp)
add = np.append(add, np.full(len(tmp), step))
add = np.append(add, np.full(len(tmp), step, dtype=np.int64))
idx = np.argsort(cdist)
cdist = cdist[idx]
add = add[idx]
......@@ -85,21 +85,19 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
cdist *= corfac
wgt = np.diff(cdist)
mdist = 0.5*(cdist[:-1]+cdist[1:])
wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], erf)
wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], sig[i], erf)
add = np.append(pos1, add)
add = np.cumsum(add)
out[i] = (add, wgt)
return out
def apply_erf(wgt, dist, lo, mid, hi, erf):
def apply_erf(wgt, dist, lo, mid, hi, sig, erf):
wgt = wgt.copy()
mask = dist > hi
wgt[mask] = 0.
mask = (dist > mid) & (dist <= hi)
wgt[mask] *= erf((dist[mask]-mid)/(hi-mid))
mask = (dist <= mid) & (dist > lo)
wgt[mask] *= erf((dist[mask]-mid)/(mid-lo))
mask = (dist > lo) & (dist <= hi)
wgt[mask] *= erf((-1/dist[mask]+1/mid)/sig)
return wgt
......@@ -119,17 +117,32 @@ class LOSResponse(LinearOperator):
of sight. The first dimension must have as many entries as `domain`
has dimensions. The second dimensions must be identical for both arrays
and indicated the total number of lines of sight.
sigmas_low, sigmas_up : numpy.ndarray(float) (optional)
For expert use. If unsure, leave blank.
sigmas: numpy.ndarray(float) (optional)
If this is not None, the inverse of the lengths of the LOSs are assumed
to be Gaussian distributed with these sigmas. The start point will
remain the same, but the endpoint is assumed to be unknown.
This is a typical statistical model for astrophysical parallaxes.
The LOS response then returns the expected integral
over the input given that the length of the LOS is unknown and
therefore the result is averaged over different endpoints.
default: None
truncation: float (optional)
Use only if the sigmas keyword argument is used!
This truncates the probability of the endpoint lying more sigmas away
than the truncation. Used to speed up computation and to avoid negative
distances. It should hold that `1./(1./length-sigma*truncation)>0`
for all lengths of the LOSs and all corresponding sigma of sigmas.
If unsure, leave blank.
default: 3.
Notes
-----
`starts, `ends`, `sigmas_low`, and `sigmas_up` have to be identical on
`starts, `ends`, `sigmas`, and `truncation` have to be identical on
every calling MPI task (i.e. the full LOS information has to be provided on
every task).
"""
def __init__(self, domain, starts, ends, sigmas_low=None, sigmas_up=None):
def __init__(self, domain, starts, ends, sigmas=None, truncation=3.):
self._domain = DomainTuple.make(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES
......@@ -141,20 +154,15 @@ class LOSResponse(LinearOperator):
starts = np.array(starts)
nlos = starts.shape[1]
ends = np.array(ends)
if sigmas_low is None:
sigmas_low = np.zeros(nlos, dtype=np.float32)
if sigmas_up is None:
sigmas_up = np.zeros(nlos, dtype=np.float32)
sigmas_low = np.array(sigmas_low)
sigmas_up = np.array(sigmas_up)
if sigmas is None:
sigmas = np.zeros(nlos, dtype=np.float32)
sigmas = np.array(sigmas)
if starts.shape[0] != ndim:
raise TypeError("dimension mismatch")
if nlos != sigmas_low.shape[0]:
if nlos != sigmas.shape[0]:
raise TypeError("dimension mismatch")
if starts.shape != ends.shape:
raise TypeError("dimension mismatch")
if sigmas_low.shape != sigmas_up.shape:
raise TypeError("dimension mismatch")
self._local_shape = dobj.local_shape(self.domain[0].shape)
local_zero_point = (np.array(
......@@ -164,7 +172,11 @@ class LOSResponse(LinearOperator):
diffs = ends-starts
difflen = np.linalg.norm(diffs, axis=0)
diffs /= difflen
real_ends = ends + sigmas_up*diffs
real_distances = 1./(1./difflen - truncation*sigmas)
if np.any(real_distances < 0):
raise ValueError("parallax error truncation to high: "
"getting negative distances")
real_ends = starts + diffs*real_distances
lzp = local_zero_point.reshape((-1, 1))
dist = np.array(self.domain[0].distances).reshape((-1, 1))
localized_pixel_starts = (starts-lzp)/dist + 0.5
......@@ -175,8 +187,11 @@ class LOSResponse(LinearOperator):
localized_pixel_ends,
self._local_shape,
np.array(self.domain[0].distances),
difflen-sigmas_low, difflen, difflen+sigmas_up,
_gaussian_error_function)
1./(1./difflen+truncation*sigmas),
difflen,
1./(1./difflen-truncation*sigmas),
sigmas,
_gaussian_sf)
boxsz = 16
nlos = len(w_i)
......@@ -202,7 +217,7 @@ class LOSResponse(LinearOperator):
tmp += (fullidx[j]//boxsz)*fct
fct *= self._local_shape[j]
tmp += cnt/float(nlos)
tmp += iarr[ofs:ofs+nval]/float(nlos*npix)
tmp += iarr[ofs:ofs+nval]/(float(nlos)*float(npix))
pri[ofs:ofs+nval] = tmp
ofs += nval
cnt += 1
......
......@@ -15,20 +15,12 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from builtins import str
import pytest
import numpy as np
from parameterized import parameterized
np.seterr(all='raise', under='ignore')
def list2fixture(lst):
@pytest.fixture(params=lst)
def myfixture(request):
return request.param
def _custom_name_func(testcase_func, param_num, param):
return "{}_{}".format(
testcase_func.__name__,
parameterized.to_safe_name("_".join(str(x) for x in param.args)),
)
def expand(*args, **kwargs):
return parameterized.expand(*args, func=_custom_name_func, **kwargs)
return myfixture
# 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-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import unittest
from itertools import product
from test.common import expand
import nifty5 as ift
import numpy as np
class Energy_Tests(unittest.TestCase):
def make_operator(self, **kwargs):
np.random.seed(kwargs['seed'])
S = ift.ScalingOperator(1., kwargs['space'])
s = S.draw_sample()
return ift.MultiField.from_dict({kwargs['space_key']: s})
@expand(product(
[ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]
))
def testGaussian(self, space, seed):
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
energy = ift.GaussianEnergy(domain=space)
ift.extra.check_value_gradient_consistency(energy, op)
# @expand(product(
# [ift.GLSpace(15),
# ift.RGSpace(64, distances=.789),
# ift.RGSpace([32, 32], distances=.789)],
# [4, 78, 23]
# ))
# def testQuadratic(self, type1, space, seed):
# np.random.seed(seed)
# S = ift.ScalingOperator(1., space)
# s = [S.draw_sample() for _ in range(3)]
# energy = ift.QuadraticEnergy(s[0], ift.makeOp(s[1]), s[2])
# ift.extra.check_value_gradient_consistency(energy)
@expand(
product([
ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)
], [4, 78, 23]))
def testInverseGammaLikelihood(self, space, seed):
op = self.make_operator(space_key='s1', space=space, seed=seed)['s1']
op = op.exp()
d = np.random.normal(10, size=space.shape)**2
d = ift.Field.from_global_data(space, d)
energy = ift.InverseGammaLikelihood(d)
ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
@expand(product(
[ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]
))
def testPoissonian(self, space, seed):
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
op = op.exp()
d = np.random.poisson(120, size=space.shape)
d = ift.Field.from_global_data(space, d)
energy = ift.PoissonianEnergy(d)
ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
@expand(product(
[ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]
))
def testHamiltonian_and_KL(self, space, seed):
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
op = op.exp()
lh = ift.GaussianEnergy(domain=space)
hamiltonian = ift.Hamiltonian(lh)
ift.extra.check_value_gradient_consistency(hamiltonian, op)
S = ift.ScalingOperator(1., space)
samps = [S.draw_sample() for i in range(3)]
kl = ift.SampledKullbachLeiblerDivergence(hamiltonian, samps)
ift.extra.check_value_gradient_consistency(kl, op)
@expand(product(
[ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]
))
def testBernoulli(self, space, seed):
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
op = op.sigmoid()
d = np.random.binomial(1, 0.1, size=space.shape)
d = ift.Field.from_global_data(space, d)
energy = ift.BernoulliEnergy(d)
ift.extra.check_value_gradient_consistency(energy, op, tol=1e-6)
# 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-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
import pytest
import nifty5 as ift
from itertools import product
# Currently it is not possible to parametrize fixtures. But this will
# hopefully be fixed in the future.
# https://docs.pytest.org/en/latest/proposals/parametrize_with_fixtures.html
SPACES = [
ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)
]
SEEDS = [4, 78, 23]
PARAMS = product(SEEDS, SPACES)
@pytest.fixture(params=PARAMS)
def field(request):
np.random.seed(request.param[0])
S = ift.ScalingOperator(1., request.param[1])
s = S.draw_sample()
return ift.MultiField.from_dict({'s1': s})['s1']
def test_gaussian(field):
energy = ift.GaussianEnergy(domain=field.domain)
ift.extra.check_value_gradient_consistency(energy, field)
def test_inverse_gamma(field):
field = field.exp()
space = field.domain
d = np.random.normal(10, size=space.shape)**2
d = ift.Field.from_global_data(space, d)
energy = ift.InverseGammaLikelihood(d)
ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7)
def testPoissonian(field):
field = field.exp()
space = field.domain
d = np.random.poisson(120, size=space.shape)
d = ift.Field.from_global_data(space, d)
energy = ift.PoissonianEnergy(d)
ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7)
def test_hamiltonian_and_KL(field):
field = field.exp()
space = field.domain
lh = ift.GaussianEnergy(domain=space)
hamiltonian = ift.Hamiltonian(lh)
ift.extra.check_value_gradient_consistency(hamiltonian, field)
S = ift.ScalingOperator(1., space)
samps = [S.draw_sample() for i in range(3)]
kl = ift.SampledKullbachLeiblerDivergence(hamiltonian, samps)
ift.extra.check_value_gradient_consistency(kl, field)
def test_bernoulli(field):
field = field.sigmoid()
space = field.domain
d = np.random.binomial(1, 0.1, size=space.shape)
d = ift.Field.from_global_data(space, d)
energy = ift.BernoulliEnergy(d)
ift.extra.check_value_gradient_consistency(energy, field, tol=1e-6)
This diff is collapsed.
......@@ -15,61 +15,67 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import unittest
from itertools import product
from test.common import expand
import numpy as np
import pytest
import nifty5 as ift
import numpy as np
def _flat_PS(k):
return np.ones_like(k)
class Energy_Tests(unittest.TestCase):
@expand(product([ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
["tanh", "exp", ""],
[1, 1e-2, 1e2],
[4, 78, 23]))
def testGaussianEnergy(self, space, nonlinearity, noise, seed):
np.random.seed(seed)
dim = len(space.shape)
hspace = space.get_default_codomain()
ht = ift.HarmonicTransformOperator(hspace, target=space)
binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
pspace = ift.PowerSpace(hspace, binbounds=binbounds)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
xi0_var = ift.Linearization.make_var(xi0)
pmp = pytest.mark.parametrize
def pspec(k): return 1 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
A = Dist(ift.sqrt(pspec))
N = ift.ScalingOperator(noise, space)
n = N.draw_sample()
s = ht(ift.makeOp(A)(xi0_var))
R = ift.ScalingOperator(10., space)
def d_model():
if nonlinearity == "":
return R(ht(ift.makeOp(A)))
else:
tmp = ht(ift.makeOp(A))
nonlin = getattr(tmp, nonlinearity)()
return R(nonlin)
@pmp('space', [
ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)
])
@pmp('nonlinearity', ["tanh", "exp", ""])
@pmp('noise', [1, 1e-2, 1e2])
@pmp('seed', [4, 78, 23])
def test_gaussian_energy(space, nonlinearity, noise, seed):
np.random.seed(seed)
dim = len(space.shape)
hspace = space.get_default_codomain()
ht = ift.HarmonicTransformOperator(hspace, target=space)
binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
pspace = ift.PowerSpace(hspace, binbounds=binbounds)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
# FIXME Needed?
xi0_var = ift.Linearization.make_var(xi0)
d = d_model()(xi0) + n
def pspec(k):
return 1/(1 + k**2)**dim
if noise == 1:
N = None
pspec = ift.PS_field(pspace, pspec)
A = Dist(ift.sqrt(pspec))
N = ift.ScalingOperator(noise, space)
n = N.draw_sample()
# FIXME Needed?
s = ht(ift.makeOp(A)(xi0_var))
R = ift.ScalingOperator(10., space)
energy = ift.GaussianEnergy(d, N)(d_model())
def d_model():
if nonlinearity == "":
ift.extra.check_value_gradient_metric_consistency(
energy, xi0, ntries=10)
return R(ht(ift.makeOp(A)))
else:
ift.extra.check_value_gradient_consistency(
energy, xi0, ntries=10, tol=5e-8)
tmp = ht(ift.makeOp(A))
nonlin = getattr(tmp, nonlinearity)()
return R(nonlin)
d = d_model()(xi0) + n
if noise == 1:
N = None
energy = ift.GaussianEnergy(d, N)(d_model())
if nonlinearity == "":
ift.extra.check_value_gradient_metric_consistency(
energy, xi0, ntries=10)
else:
ift.extra.check_value_gradient_consistency(
energy, xi0, ntries=10, tol=5e-8)
# 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-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import unittest
from itertools import product
from test.common import expand
import nifty5 as ift
import numpy as np
from unittest import SkipTest
from numpy.testing import assert_allclose, assert_equal
IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
spaces = [ift.RGSpace([1024], distances=0.123), ift.HPSpace(32)]
minimizers = ['ift.VL_BFGS(IC)',
'ift.NonlinearCG(IC, "Polak-Ribiere")',
# 'ift.NonlinearCG(IC, "Hestenes-Stiefel"),
'ift.NonlinearCG(IC, "Fletcher-Reeves")',
'ift.NonlinearCG(IC, "5.49")',
'ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=1000)',
'ift.L_BFGS(IC)',
'ift.NewtonCG(IC)']
newton_minimizers = ['ift.RelaxedNewton(IC)']
quadratic_only_minimizers = ['ift.ConjugateGradient(IC)',
'ift.ScipyCG(tol=1e-5, maxiter=300)']
slow_minimizers = ['ift.SteepestDescent(IC)']
class Test_Minimizers(unittest.TestCase):
@expand(product(minimizers + newton_minimizers +
quadratic_only_minimizers + slow_minimizers, spaces))
def test_quadratic_minimization(self, minimizer, space):
np.random.seed(42)
starting_point = ift.Field.from_random('normal', domain=space)*10
covariance_diagonal = ift.Field.from_random(