test_noise.py 3.52 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
        var = 1.
Philipp Arras's avatar
Philipp Arras committed
50
51
52
53
        n = ift.Field.from_random(
            domain=space,
            random_type='normal',
            std=np.sqrt(var))
54
55
56
57
        var = ift.Field(n.domain, val=var)
        N = ift.DiagonalOperator(var)
        eta0 = ift.log(var)
        s = ht(xi * A)
Martin Reinecke's avatar
Martin Reinecke committed
58
        R = ift.ScalingOperator(10., space)
Philipp Arras's avatar
Philipp Arras committed
59
60
61
62
63
64
65
        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())
66
        eps = 1e-8
Philipp Arras's avatar
Philipp Arras committed
67
68
69
70
71
72
73
74
75
76
77
78
79
        eta1 = eta0 + eps * direction

        IC = ift.GradientNormController(
            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,
80
            ht=ht,
Philipp Arras's avatar
Philipp Arras committed
81
82
83
84
            power=A,
            N=N,
            S=S,
            inverter=inverter).curvature
85
        Nsamples = 10
Philipp Arras's avatar
Philipp Arras committed
86
87
88
        xi_sample_list = [
            D.generate_posterior_sample() +
            xi for i in range(Nsamples)]
Philipp Arras's avatar
Philipp Arras committed
89
90

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

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