test_minimizers.py 7.34 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/>.
#
Martin Reinecke's avatar
Martin Reinecke committed
14
# Copyright(C) 2013-2018 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
15
16
17
18
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

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

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

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

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

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

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

46

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

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

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

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

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

73
74
    @expand(product(minimizers+newton_minimizers))
    def test_rosenbrock(self, minimizer):
Martin Reinecke's avatar
Martin Reinecke committed
75
76
77
78
79
        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
80
        space = ift.DomainTuple.make(ift.UnstructuredDomain((2,)))
Martin Reinecke's avatar
Martin Reinecke committed
81
82
83
84
85
86
87
88
        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):
89
                return rosen(self._position.to_global_data_rw())
Martin Reinecke's avatar
Martin Reinecke committed
90
91
92

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

            @property
Martin Reinecke's avatar
Martin Reinecke committed
98
            def metric(self):
Martin Reinecke's avatar
Martin Reinecke committed
99
100
                class RBCurv(ift.EndomorphicOperator):
                    def __init__(self, loc):
101
                        self._loc = loc.to_global_data_rw()
Martin Reinecke's avatar
Martin Reinecke committed
102
                        self._capability = self.TIMES
103
                        self._domain = space
Martin Reinecke's avatar
Martin Reinecke committed
104
105
106

                    def apply(self, x, mode):
                        self._check_input(x, mode)
107
                        inp = x.to_global_data_rw()
Martin Reinecke's avatar
Martin Reinecke committed
108
109
110
                        out = ift.Field.from_global_data(
                            space, rosen_hess_prod(self._loc.copy(), inp))
                        return out
Martin Reinecke's avatar
Martin Reinecke committed
111

Martin Reinecke's avatar
Martin Reinecke committed
112
113
                t1 = ift.GradientNormController(tol_abs_gradnorm=1e-5,
                                                iteration_limit=1000)
114
                return ift.InversionEnabler(RBCurv(self._position), t1)
Martin Reinecke's avatar
Martin Reinecke committed
115
116

        try:
Martin Reinecke's avatar
Martin Reinecke committed
117
            minimizer = eval(minimizer)
Martin Reinecke's avatar
Martin Reinecke committed
118
119
120
121
122
123
124
            energy = RBEnergy(position=starting_point)

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

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

128
129
    @expand(product(minimizers+slow_minimizers))
    def test_gauss(self, minimizer):
130
        space = ift.UnstructuredDomain((1,))
131
        starting_point = ift.Field.full(space, 3.)
132
133
134
135
136
137
138
139
140
141
142
143
144

        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]
145
146
                return ift.Field.full(self.position.domain,
                                      2*x*np.exp(-(x**2)))
147
148

            @property
Martin Reinecke's avatar
Martin Reinecke committed
149
            def metric(self):
150
151
152
                x = self.position.to_global_data()[0]
                v = (2 - 4*x*x)*np.exp(-x**2)
                return ift.DiagonalOperator(
153
                    ift.Field.full(self.position.domain, v))
154
155

        try:
Martin Reinecke's avatar
Martin Reinecke committed
156
            minimizer = eval(minimizer)
157
158
159
160
161
162
163
            energy = ExpEnergy(position=starting_point)

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

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

    @expand(product(minimizers+newton_minimizers+slow_minimizers))
    def test_cosh(self, minimizer):
        space = ift.UnstructuredDomain((1,))
170
        starting_point = ift.Field.full(space, 3.)
171
172
173
174
175
176
177
178
179
180
181
182
183

        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]
184
                return ift.Field.full(self.position.domain, np.sinh(x))
185
186

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

        try:
Martin Reinecke's avatar
Martin Reinecke committed
194
            minimizer = eval(minimizer)
195
196
197
198
199
200
201
            energy = CoshEnergy(position=starting_point)

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

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