test_minimizers.py 7.73 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
# 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Martin Reinecke's avatar
Martin Reinecke committed
17

Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
18
import unittest
Martin Reinecke's avatar
changes    
Martin Reinecke committed
19
from itertools import product
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
20
from test.common import expand
21
22
23

import nifty5 as ift
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
24
from unittest import SkipTest
25
from numpy.testing import assert_allclose, assert_equal
Matevz, Sraml (sraml)'s avatar
Matevz, Sraml (sraml) committed
26

27
28
IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)

Martin Reinecke's avatar
changes    
Martin Reinecke committed
29
spaces = [ift.RGSpace([1024], distances=0.123), ift.HPSpace(32)]
30

Martin Reinecke's avatar
Martin Reinecke committed
31
32
minimizers = ['ift.VL_BFGS(IC)',
              'ift.NonlinearCG(IC, "Polak-Ribiere")',
33
              # 'ift.NonlinearCG(IC, "Hestenes-Stiefel"),
Martin Reinecke's avatar
Martin Reinecke committed
34
35
              'ift.NonlinearCG(IC, "Fletcher-Reeves")',
              'ift.NonlinearCG(IC, "5.49")',
Martin Reinecke's avatar
fix    
Martin Reinecke committed
36
              'ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=1000)',
37
38
              'ift.L_BFGS(IC)',
              'ift.NewtonCG(IC)']
Martin Reinecke's avatar
Martin Reinecke committed
39

Martin Reinecke's avatar
Martin Reinecke committed
40
newton_minimizers = ['ift.RelaxedNewton(IC)']
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
41
42
quadratic_only_minimizers = ['ift.ConjugateGradient(IC)',
                             'ift.ScipyCG(tol=1e-5, maxiter=300)']
Martin Reinecke's avatar
Martin Reinecke committed
43
slow_minimizers = ['ift.SteepestDescent(IC)']
44

45

Martin Reinecke's avatar
changes    
Martin Reinecke committed
46
class Test_Minimizers(unittest.TestCase):
47

Martin Reinecke's avatar
Martin Reinecke committed
48
49
    @expand(product(minimizers + newton_minimizers +
                    quadratic_only_minimizers + slow_minimizers, spaces))
50
    def test_quadratic_minimization(self, minimizer, space):
51
        np.random.seed(42)
Martin Reinecke's avatar
changes    
Martin Reinecke committed
52
53
54
        starting_point = ift.Field.from_random('normal', domain=space)*10
        covariance_diagonal = ift.Field.from_random(
                                  'uniform', domain=space) + 0.5
55
        covariance = ift.DiagonalOperator(covariance_diagonal)
Martin Reinecke's avatar
step 1    
Martin Reinecke committed
56
        required_result = ift.full(space, 1.)
57

Martin Reinecke's avatar
Martin Reinecke committed
58
        try:
Martin Reinecke's avatar
Martin Reinecke committed
59
            minimizer = eval(minimizer)
Martin Reinecke's avatar
Martin Reinecke committed
60
61
            energy = ift.QuadraticEnergy(A=covariance, b=required_result,
                                         position=starting_point)
62

Martin Reinecke's avatar
Martin Reinecke committed
63
64
65
66
67
            (energy, convergence) = minimizer(energy)
        except NotImplementedError:
            raise SkipTest

        assert_equal(convergence, IC.CONVERGED)
68
69
        assert_allclose(energy.position.local_data,
                        1./covariance_diagonal.local_data,
Martin Reinecke's avatar
changes    
Martin Reinecke committed
70
                        rtol=1e-3, atol=1e-3)
Martin Reinecke's avatar
Martin Reinecke committed
71

Martin Reinecke's avatar
Martin Reinecke committed
72
    @expand(product(minimizers+newton_minimizers))
73
    def test_rosenbrock(self, minimizer):
Martin Reinecke's avatar
Martin Reinecke committed
74
75
76
77
78
        try:
            from scipy.optimize import rosen, rosen_der, rosen_hess_prod
        except ImportError:
            raise SkipTest
        np.random.seed(42)
Martin Reinecke's avatar
Martin Reinecke committed
79
        space = ift.DomainTuple.make(ift.UnstructuredDomain((2,)))
Martin Reinecke's avatar
Martin Reinecke committed
80
81
82
83
84
85
86
87
        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):
88
                return rosen(self._position.to_global_data_rw())
Martin Reinecke's avatar
Martin Reinecke committed
89
90
91

            @property
            def gradient(self):
92
                inp = self._position.to_global_data_rw()
Martin Reinecke's avatar
Martin Reinecke committed
93
94
95
                out = ift.Field.from_global_data(space, rosen_der(inp))
                return out

Martin Reinecke's avatar
Martin Reinecke committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
            @property
            def metric(self):
                class RBCurv(ift.EndomorphicOperator):
                    def __init__(self, loc):
                        self._loc = loc.to_global_data_rw()
                        self._capability = self.TIMES
                        self._domain = space

                    def apply(self, x, mode):
                        self._check_input(x, mode)
                        inp = x.to_global_data_rw()
                        out = ift.Field.from_global_data(
                            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)

Martin Reinecke's avatar
Martin Reinecke committed
115
116
117
118
119
            def apply_metric(self, x):
                inp = x.to_global_data_rw()
                pos = self._position.to_global_data_rw()
                return ift.Field.from_global_data(
                    space, rosen_hess_prod(pos, inp))
Martin Reinecke's avatar
Martin Reinecke committed
120
121

        try:
Martin Reinecke's avatar
Martin Reinecke committed
122
            minimizer = eval(minimizer)
Martin Reinecke's avatar
Martin Reinecke committed
123
124
125
126
127
128
129
            energy = RBEnergy(position=starting_point)

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

        assert_equal(convergence, IC.CONVERGED)
130
        assert_allclose(energy.position.local_data, 1.,
Martin Reinecke's avatar
Martin Reinecke committed
131
                        rtol=1e-3, atol=1e-3)
132

133
134
    @expand(product(minimizers+slow_minimizers))
    def test_gauss(self, minimizer):
135
        space = ift.UnstructuredDomain((1,))
136
        starting_point = ift.Field.full(space, 3.)
137
138
139
140
141
142
143
144
145
146
147
148
149

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

            @property
            def value(self):
                x = self.position.to_global_data()[0]
                return -np.exp(-(x**2))

            @property
            def gradient(self):
                x = self.position.to_global_data()[0]
150
151
                return ift.Field.full(self.position.domain,
                                      2*x*np.exp(-(x**2)))
152

Martin Reinecke's avatar
Martin Reinecke committed
153
154
155
            def apply_metric(self, x):
                p = self.position.to_global_data()[0]
                v = (2 - 4*p*p)*np.exp(-p**2)
156
                return ift.DiagonalOperator(
Martin Reinecke's avatar
Martin Reinecke committed
157
                    ift.Field.full(self.position.domain, v))(x)
158
159

        try:
Martin Reinecke's avatar
Martin Reinecke committed
160
            minimizer = eval(minimizer)
161
162
163
164
165
166
167
            energy = ExpEnergy(position=starting_point)

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

        assert_equal(convergence, IC.CONVERGED)
168
        assert_allclose(energy.position.local_data, 0.,
169
                        atol=1e-3)
170

Martin Reinecke's avatar
Martin Reinecke committed
171
    @expand(product(minimizers+newton_minimizers+slow_minimizers))
172
173
    def test_cosh(self, minimizer):
        space = ift.UnstructuredDomain((1,))
174
        starting_point = ift.Field.full(space, 3.)
175
176
177
178
179
180
181
182
183
184
185
186
187

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

            @property
            def value(self):
                x = self.position.to_global_data()[0]
                return np.cosh(x)

            @property
            def gradient(self):
                x = self.position.to_global_data()[0]
188
                return ift.Field.full(self.position.domain, np.sinh(x))
189

Martin Reinecke's avatar
Martin Reinecke committed
190
191
192
193
194
195
196
            @property
            def metric(self):
                x = self.position.to_global_data()[0]
                v = np.cosh(x)
                return ift.DiagonalOperator(
                    ift.Field.full(self.position.domain, v))

Martin Reinecke's avatar
Martin Reinecke committed
197
198
199
            def apply_metric(self, x):
                p = self.position.to_global_data()[0]
                v = np.cosh(p)
200
                return ift.DiagonalOperator(
Martin Reinecke's avatar
Martin Reinecke committed
201
                    ift.Field.full(self.position.domain, v))(x)
202
203

        try:
Martin Reinecke's avatar
Martin Reinecke committed
204
            minimizer = eval(minimizer)
205
206
207
208
209
210
211
            energy = CoshEnergy(position=starting_point)

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

        assert_equal(convergence, IC.CONVERGED)
212
        assert_allclose(energy.position.local_data, 0., atol=1e-3)