test_noise.py 3.74 KB
Newer Older
Philipp Arras's avatar
Philipp Arras committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 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-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

import unittest
import nifty4 as ift
import numpy as np
from itertools import product
from test.common import expand
from numpy.testing import assert_allclose


# TODO Add also other space types


class Noise_Energy_Tests(unittest.TestCase):
    @expand(product([ift.RGSpace(64, distances=.789),
                     ift.RGSpace([32, 32], distances=.789)],
33
34
35
36
                    [ift.library.Exponential, ift.library.Linear],
                    [23, 131, 32]))
    def testNoise(self, space, nonlinearity, seed):
        np.random.seed(seed)
Philipp Arras's avatar
Philipp Arras committed
37
38
        f = nonlinearity()
        dim = len(space.shape)
39
40
        hspace = space.get_default_codomain()
        ht = ift.HarmonicTransformOperator(hspace, target=space)
Philipp Arras's avatar
Philipp Arras committed
41
42
43
44
45
46
47
48
        binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
        pspace = ift.PowerSpace(hspace, binbounds=binbounds)
        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
        xi = ift.Field.from_random(domain=hspace, random_type='normal')

        def pspec(k): return 1 / (1 + k**2)**dim
        tau = ift.PS_field(pspace, pspec)
        A = P.adjoint_times(ift.sqrt(tau))
49
50
51
52
53
54
        var = 1.
        n = ift.Field.from_random(domain=space, random_type='normal', std=np.sqrt(var))
        var = ift.Field(n.domain, val=var)
        N = ift.DiagonalOperator(var)
        eta0 = ift.log(var)
        s = ht(xi * A)
Philipp Arras's avatar
Philipp Arras committed
55
56
57
58
59
60
61
62
63
64
        diag = ift.Field.ones(space) * 10
        R = ift.DiagonalOperator(diag)
        diag = ift.Field.ones(space)
        d = R(f(s)) + n

        alpha = ift.Field(d.domain, val=2.)
        q = ift.Field(d.domain, val=1e-5)

        direction = ift.Field.from_random('normal', d.domain)
        direction /= np.sqrt(direction.var())
65
        eps = 1e-8
Philipp Arras's avatar
Philipp Arras committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
        eta1 = eta0 + eps * direction

        IC = ift.GradientNormController(
            name='IC',
            verbose=False,
            iteration_limit=100,
            tol_abs_gradnorm=1e-5)
        inverter = ift.ConjugateGradient(IC)

        S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
        D = ift.library.NonlinearWienerFilterEnergy(
            position=xi,
            d=d,
            Instrument=R,
            nonlinearity=f,
81
            ht=ht,
Philipp Arras's avatar
Philipp Arras committed
82
83
84
85
            power=A,
            N=N,
            S=S,
            inverter=inverter).curvature
86
87
        Nsamples = 10
        sample_list = [D.generate_posterior_sample() + xi for i in range(Nsamples)]
Philipp Arras's avatar
Philipp Arras committed
88
89
90
91

        energy0 = ift.library.NoiseEnergy(
            position=eta0, d=d, m=xi, D=D, t=tau, Instrument=R,
            alpha=alpha, q=q, Projection=P, nonlinearity=f,
92
            ht=ht, sample_list=sample_list)
Philipp Arras's avatar
Philipp Arras committed
93
94
95
        energy1 = ift.library.NoiseEnergy(
            position=eta1, d=d, m=xi, D=D, t=tau, Instrument=R,
            alpha=alpha, q=q, Projection=P, nonlinearity=f,
96
            ht=ht, sample_list=sample_list)
Philipp Arras's avatar
Philipp Arras committed
97
98
99

        a = (energy1.value - energy0.value) / eps
        b = energy0.gradient.vdot(direction)
100
        tol = 1e-5
Philipp Arras's avatar
Philipp Arras committed
101
        assert_allclose(a, b, rtol=tol, atol=tol)