test_minimizers.py 8.48 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
# 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
from numpy.testing import assert_allclose, assert_equal

Martin Reinecke's avatar
Martin Reinecke committed
22
import nifty7 as ift
Philipp Arras's avatar
Philipp Arras committed
23

24
from .common import setup_function, teardown_function
Philipp Arras's avatar
Philipp Arras committed
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

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
44
45
    'ift.ConjugateGradient(IC)',
    'ift.minimization.scipy_minimizer._ScipyCG(tol=1e-5, maxiter=300)'
Philipp Arras's avatar
Philipp Arras committed
46
]
47
slow_minimizers = ['ift.SteepestDescent(IC)', 'ift.ADVIOptimizer(10, resample=False)']
Philipp Arras's avatar
Philipp Arras committed
48
49
50
51
52
53


@pmp('minimizer', minimizers + newton_minimizers + quadratic_only_minimizers +
     slow_minimizers)
@pmp('space', spaces)
def test_quadratic_minimization(minimizer, space):
Martin Reinecke's avatar
merge    
Martin Reinecke committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
    with ift.random.Context(98765):
        starting_point = ift.Field.from_random(domain=space, random_type='normal') * 10
        covariance_diagonal = ift.Field.from_random(domain=space, random_type='uniform') + 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:
            pytest.skip()
    
        assert_equal(convergence, IC.CONVERGED)
        assert_allclose(
            energy.position.val,
            1./covariance_diagonal.val,
            rtol=1e-3,
            atol=1e-3)
Philipp Arras's avatar
Philipp Arras committed
75

Martin Reinecke's avatar
Martin Reinecke committed
76

77
78
79
80
@pmp('space', spaces)
def test_WF_curvature(space):
    required_result = ift.full(space, 1.)

81
    s = ift.Field.from_random(domain=space, random_type='uniform') + 0.5
82
    S = ift.DiagonalOperator(s)
83
    r = ift.Field.from_random(domain=space, random_type='uniform')
84
    R = ift.DiagonalOperator(r)
85
    n = ift.Field.from_random(domain=space, random_type='uniform') + 0.5
86
87
    N = ift.DiagonalOperator(n)
    all_diag = 1./s + r**2/n
Martin Reinecke's avatar
Martin Reinecke committed
88
    curv = ift.WienerFilterCurvature(R, N, S, iteration_controller=IC,
Philipp Arras's avatar
Philipp Arras committed
89
90
91
                                     iteration_controller_sampling=IC,
                                     data_sampling_dtype=np.float64,
                                     prior_sampling_dtype=np.float64)
92
93
    m = curv.inverse(required_result)
    assert_allclose(
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
94
95
        m.val,
        1./all_diag.val,
96
97
        rtol=1e-3,
        atol=1e-3)
Philipp Arras's avatar
Philipp Arras committed
98
99
    curv.draw_sample()
    curv.draw_sample(from_inverse=True)
100

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


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

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

        @property
        def value(self):
Martin Reinecke's avatar
Martin Reinecke committed
136
            return rosen(self._position.val_rw())
Philipp Arras's avatar
Philipp Arras committed
137
138
139

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

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

                def apply(self, x, mode):
                    self._check_input(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
154
                    inp = x.val_rw()
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
155
                    out = ift.Field(
Philipp Arras's avatar
Philipp Arras committed
156
157
158
159
160
161
162
163
                        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
164
165
            inp = x.val_rw()
            pos = self._position.val_rw()
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
166
            return ift.Field(space, rosen_hess_prod(pos, inp))
Philipp Arras's avatar
Philipp Arras committed
167
168
169
170
171
172
173

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

        (energy, convergence) = minimizer(energy)
    except NotImplementedError:
Philipp Arras's avatar
Philipp Arras committed
174
        pytest.skip()
Philipp Arras's avatar
Philipp Arras committed
175
176

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


@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
191
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
192
193
194
195
            return -np.exp(-(x**2))

        @property
        def gradient(self):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
196
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
197
198
199
            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
200
            p = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
201
202
203
204
205
206
207
208
209
210
            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:
Philipp Arras's avatar
Philipp Arras committed
211
        pytest.skip()
Philipp Arras's avatar
Philipp Arras committed
212
213

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


@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
228
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
229
230
231
232
            return np.cosh(x)

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

        @property
        def metric(self):
Martin Reinecke's avatar
stage2    
Martin Reinecke committed
238
            x = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
239
240
241
242
243
            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
244
            p = self.position.val[0]
Philipp Arras's avatar
Philipp Arras committed
245
246
247
248
249
250
251
252
253
254
            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:
Philipp Arras's avatar
Philipp Arras committed
255
        pytest.skip()
Philipp Arras's avatar
Philipp Arras committed
256
257

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