test_minimizers.py 8 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
# 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.

from unittest import SkipTest

import numpy as np
import pytest
from numpy.testing import assert_allclose, assert_equal

Martin Reinecke's avatar
5->6    
Martin Reinecke committed
24
import nifty6 as ift
25
from .common import setup_function, teardown_function
Philipp Arras's avatar
Philipp Arras committed
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

pmp = pytest.mark.parametrize
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 = [
Martin Reinecke's avatar
Martin Reinecke committed
45
46
    'ift.ConjugateGradient(IC)',
    'ift.minimization.scipy_minimizer._ScipyCG(tol=1e-5, maxiter=300)'
Philipp Arras's avatar
Philipp Arras committed
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
]
slow_minimizers = ['ift.SteepestDescent(IC)']


@pmp('minimizer', minimizers + newton_minimizers + quadratic_only_minimizers +
     slow_minimizers)
@pmp('space', spaces)
def test_quadratic_minimization(minimizer, space):
    starting_point = ift.Field.from_random('normal', domain=space)*10
    covariance_diagonal = ift.Field.from_random('uniform', domain=space) + 0.5
    covariance = ift.DiagonalOperator(covariance_diagonal)
    required_result = ift.full(space, 1.)

    try:
        minimizer = eval(minimizer)
        energy = ift.QuadraticEnergy(
            A=covariance, b=required_result, position=starting_point)

        (energy, convergence) = minimizer(energy)
    except NotImplementedError:
        raise SkipTest

    assert_equal(convergence, IC.CONVERGED)
    assert_allclose(
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
71
72
        energy.position.val,
        1./covariance_diagonal.val,
Philipp Arras's avatar
Philipp Arras committed
73
74
75
        rtol=1e-3,
        atol=1e-3)

Martin Reinecke's avatar
Martin Reinecke committed
76

77
78
79
80
81
82
83
84
85
86
87
@pmp('space', spaces)
def test_WF_curvature(space):
    required_result = ift.full(space, 1.)

    s = ift.Field.from_random('uniform', domain=space) + 0.5
    S = ift.DiagonalOperator(s)
    r = ift.Field.from_random('uniform', domain=space)
    R = ift.DiagonalOperator(r)
    n = ift.Field.from_random('uniform', domain=space) + 0.5
    N = ift.DiagonalOperator(n)
    all_diag = 1./s + r**2/n
Martin Reinecke's avatar
Martin Reinecke committed
88
89
    curv = ift.WienerFilterCurvature(R, N, S, iteration_controller=IC,
                                     iteration_controller_sampling=IC)
90
91
    m = curv.inverse(required_result)
    assert_allclose(
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
92
93
        m.val,
        1./all_diag.val,
94
95
96
97
98
        rtol=1e-3,
        atol=1e-3)
    curv.draw_sample()
    curv.draw_sample(from_inverse=True)

99
100
101
102
103
    if len(space.shape) == 1:
        R = ift.ValueInserter(space, [0])
        n = ift.from_random('uniform', R.domain) + 0.5
        N = ift.DiagonalOperator(n)
        all_diag = 1./s + R(1/n)
Martin Reinecke's avatar
Martin Reinecke committed
104
105
106
        curv = ift.WienerFilterCurvature(R.adjoint, N, S,
                                         iteration_controller=IC,
                                         iteration_controller_sampling=IC)
107
108
        m = curv.inverse(required_result)
        assert_allclose(
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
109
110
            m.val,
            1./all_diag.val,
111
112
113
114
115
116
            rtol=1e-3,
            atol=1e-3)
        curv.draw_sample()
        curv.draw_sample(from_inverse=True)


Philipp Arras's avatar
Philipp Arras committed
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
@pmp('minimizer', minimizers + newton_minimizers)
def test_rosenbrock(minimizer):
    try:
        from scipy.optimize import rosen, rosen_der, rosen_hess_prod
    except ImportError:
        raise SkipTest
    space = ift.DomainTuple.make(ift.UnstructuredDomain((2,)))
    starting_point = ift.Field.from_random('normal', domain=space)*10

    class RBEnergy(ift.Energy):
        def __init__(self, position):
            super(RBEnergy, self).__init__(position)

        @property
        def value(self):
Martin Reinecke's avatar
Martin Reinecke committed
132
            return rosen(self._position.val_rw())
Philipp Arras's avatar
Philipp Arras committed
133
134
135

        @property
        def gradient(self):
Martin Reinecke's avatar
Martin Reinecke committed
136
            inp = self._position.val_rw()
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
137
            out = ift.Field(space, rosen_der(inp))
Philipp Arras's avatar
Philipp Arras committed
138
139
140
141
142
143
            return out

        @property
        def metric(self):
            class RBCurv(ift.EndomorphicOperator):
                def __init__(self, loc):
Martin Reinecke's avatar
Martin Reinecke committed
144
                    self._loc = loc.val_rw()
Philipp Arras's avatar
Philipp Arras committed
145
146
147
148
149
                    self._capability = self.TIMES
                    self._domain = space

                def apply(self, x, mode):
                    self._check_input(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
150
                    inp = x.val_rw()
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
151
                    out = ift.Field(
Philipp Arras's avatar
Philipp Arras committed
152
153
154
155
156
157
158
159
                        space, rosen_hess_prod(self._loc.copy(), inp))
                    return out

            t1 = ift.GradientNormController(
                tol_abs_gradnorm=1e-5, iteration_limit=1000)
            return ift.InversionEnabler(RBCurv(self._position), t1)

        def apply_metric(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
160
161
            inp = x.val_rw()
            pos = self._position.val_rw()
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
162
            return ift.Field(space, rosen_hess_prod(pos, inp))
Philipp Arras's avatar
Philipp Arras committed
163
164
165
166
167
168
169
170
171
172

    try:
        minimizer = eval(minimizer)
        energy = RBEnergy(position=starting_point)

        (energy, convergence) = minimizer(energy)
    except NotImplementedError:
        raise SkipTest

    assert_equal(convergence, IC.CONVERGED)
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
173
    assert_allclose(energy.position.val, 1., rtol=1e-3, atol=1e-3)
Philipp Arras's avatar
Philipp Arras committed
174
175
176
177
178
179
180
181
182
183
184
185
186


@pmp('minimizer', minimizers + slow_minimizers)
def test_gauss(minimizer):
    space = ift.UnstructuredDomain((1,))
    starting_point = ift.Field.full(space, 3.)

    class ExpEnergy(ift.Energy):
        def __init__(self, position):
            super(ExpEnergy, self).__init__(position)

        @property
        def value(self):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
187
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
188
189
190
191
            return -np.exp(-(x**2))

        @property
        def gradient(self):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
192
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
193
194
195
            return ift.Field.full(self.position.domain, 2*x*np.exp(-(x**2)))

        def apply_metric(self, x):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
196
            p = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
197
198
199
200
201
202
203
204
205
206
207
208
209
            v = (2 - 4*p*p)*np.exp(-p**2)
            return ift.DiagonalOperator(
                ift.Field.full(self.position.domain, v))(x)

    try:
        minimizer = eval(minimizer)
        energy = ExpEnergy(position=starting_point)

        (energy, convergence) = minimizer(energy)
    except NotImplementedError:
        raise SkipTest

    assert_equal(convergence, IC.CONVERGED)
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
210
    assert_allclose(energy.position.val, 0., atol=1e-3)
Philipp Arras's avatar
Philipp Arras committed
211
212
213
214
215
216
217
218
219
220
221
222
223


@pmp('minimizer', minimizers + newton_minimizers + slow_minimizers)
def test_cosh(minimizer):
    space = ift.UnstructuredDomain((1,))
    starting_point = ift.Field.full(space, 3.)

    class CoshEnergy(ift.Energy):
        def __init__(self, position):
            super(CoshEnergy, self).__init__(position)

        @property
        def value(self):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
224
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
225
226
227
228
            return np.cosh(x)

        @property
        def gradient(self):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
229
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
230
231
232
233
            return ift.Field.full(self.position.domain, np.sinh(x))

        @property
        def metric(self):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
234
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
235
236
237
238
239
            v = np.cosh(x)
            return ift.DiagonalOperator(
                ift.Field.full(self.position.domain, v))

        def apply_metric(self, x):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
240
            p = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
241
242
243
244
245
246
247
248
249
250
251
252
253
            v = np.cosh(p)
            return ift.DiagonalOperator(
                ift.Field.full(self.position.domain, v))(x)

    try:
        minimizer = eval(minimizer)
        energy = CoshEnergy(position=starting_point)

        (energy, convergence) = minimizer(energy)
    except NotImplementedError:
        raise SkipTest

    assert_equal(convergence, IC.CONVERGED)
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
254
    assert_allclose(energy.position.val, 0., atol=1e-3)