energy_and_model_tests.py 7.25 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 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-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

19
20
from __future__ import absolute_import, division, print_function
from ..compat import *
21
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
22
from ..sugar import from_random
Martin Reinecke's avatar
tweak    
Martin Reinecke committed
23
24
from ..minimization.energy import Energy
from ..models.model import Model
Martin Reinecke's avatar
Martin Reinecke committed
25
from ..linearization import Linearization
26

Martin Reinecke's avatar
Martin Reinecke committed
27
__all__ = ["check_value_gradient_consistency",
Martin Reinecke's avatar
Martin Reinecke committed
28
29
30
           "check_value_gradient_metric_consistency",
           "check_value_gradient_metric_consistency2",
           "check_value_gradient_consistency2"]
31

Philipp Arras's avatar
Philipp Arras committed
32

33
34
35
36
37
def _get_acceptable_model(M):
    val = M.value
    if not np.isfinite(val.sum()):
        raise ValueError('Initial Model value must be finite')
    dir = from_random("normal", M.position.domain)
38
    dirder = M.jacobian(dir)
39
    if dirder.norm() == 0:
40
        dir = dir * val.norm() * 1e-5
41
    else:
42
        dir = dir * val.norm() * (1e-5/dirder.norm())
Philipp Arras's avatar
Philipp Arras committed
43
    # Find a step length that leads to a "reasonable" Model
44
45
46
47
48
49
50
    for i in range(50):
        try:
            M2 = M.at(M.position+dir)
            if np.isfinite(M2.value.sum()) and abs(M2.value.sum()) < 1e20:
                break
        except FloatingPointError:
            pass
51
        dir = dir*0.5
52
53
54
55
    else:
        raise ValueError("could not find a reasonable initial step")
    return M2

56

Martin Reinecke's avatar
Martin Reinecke committed
57
def _get_acceptable_energy(E):
58
59
    val = E.value
    if not np.isfinite(val):
60
        raise ValueError('Initial Energy must be finite')
61
    dir = from_random("normal", E.position.domain)
62
    dirder = E.gradient.vdot(dir)
63
    dir = dir * (np.abs(val)/np.abs(dirder)*1e-5)
Philipp Arras's avatar
Philipp Arras committed
64
    # Find a step length that leads to a "reasonable" energy
Martin Reinecke's avatar
Martin Reinecke committed
65
66
67
68
69
70
71
    for i in range(50):
        try:
            E2 = E.at(E.position+dir)
            if np.isfinite(E2.value) and abs(E2.value) < 1e20:
                break
        except FloatingPointError:
            pass
72
        dir = dir*0.5
Martin Reinecke's avatar
Martin Reinecke committed
73
74
75
76
    else:
        raise ValueError("could not find a reasonable initial step")
    return E2

Martin Reinecke's avatar
Martin Reinecke committed
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
def _get_acceptable_location(op, loc, lin):
    val = lin.val
    if not np.isfinite(val.sum()):
        raise ValueError('Initial value must be finite')
    dir = from_random("normal", loc.domain)
    dirder = lin.jac(dir)
    if dirder.norm() == 0:
        dir = dir * val.norm() * 1e-5
    else:
        dir = dir * val.norm() * (1e-5/dirder.norm())
    # Find a step length that leads to a "reasonable" location
    for i in range(50):
        try:
            loc2 = loc+dir
            lin2 = op(Linearization.make_var(loc2))
            if np.isfinite(lin2.val.sum()) and abs(lin2.val.sum()) < 1e20:
                break
        except FloatingPointError:
            pass
        dir = dir*0.5
    else:
        raise ValueError("could not find a reasonable initial step")
    return loc2, lin2

Martin Reinecke's avatar
Martin Reinecke committed
101

Martin Reinecke's avatar
Martin Reinecke committed
102
def check_value_gradient_consistency(E, tol=1e-8, ntries=100):
103
    for _ in range(ntries):
104
105
106
107
        if isinstance(E, Energy):
            E2 = _get_acceptable_energy(E)
        else:
            E2 = _get_acceptable_model(E)
108
        val = E.value
Martin Reinecke's avatar
Martin Reinecke committed
109
        dir = E2.position - E.position
Martin Reinecke's avatar
tweak    
Martin Reinecke committed
110
        Enext = E2
Martin Reinecke's avatar
Martin Reinecke committed
111
        dirnorm = dir.norm()
112
        for i in range(50):
Martin Reinecke's avatar
Martin Reinecke committed
113
            Emid = E.at(E.position + 0.5*dir)
114
115
116
            if isinstance(E, Energy):
                dirder = Emid.gradient.vdot(dir)/dirnorm
            else:
117
                dirder = Emid.jacobian(dir)/dirnorm
118
            numgrad = (E2.value-val)/dirnorm
119
            if isinstance(E, Model):
120
                xtol = tol * dirder.norm() / np.sqrt(dirder.size)
121
                if (abs(numgrad-dirder) <= xtol).all():
122
123
124
                    break
            else:
                xtol = tol*Emid.gradient_norm
125
                if abs(numgrad-dirder) <= xtol:
126
                    break
127
            dir = dir*0.5
Martin Reinecke's avatar
Martin Reinecke committed
128
            dirnorm *= 0.5
Martin Reinecke's avatar
Martin Reinecke committed
129
            E2 = Emid
130
        else:
Martin Reinecke's avatar
Martin Reinecke committed
131
            raise ValueError("gradient and value seem inconsistent")
Martin Reinecke's avatar
tweak    
Martin Reinecke committed
132
        E = Enext
Martin Reinecke's avatar
Martin Reinecke committed
133

Martin Reinecke's avatar
Martin Reinecke committed
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
def check_value_gradient_consistency2(op, loc, tol=1e-8, ntries=100):
    for _ in range(ntries):
        lin = op(Linearization.make_var(loc))
        loc2, lin2 = _get_acceptable_location(op, loc, lin)
        val = lin.val
        dir = loc2 - loc
        locnext = loc2
        dirnorm = dir.norm()
        for i in range(50):
            locmid = loc + 0.5*dir
            linmid = op(Linearization.make_var(locmid))
            dirder = linmid.jac(dir)/dirnorm
            numgrad = (lin2.val-val)/dirnorm
            xtol = tol * dirder.norm() / np.sqrt(dirder.size)
            if (abs(numgrad-dirder) <= xtol).all():
                break
            dir = dir*0.5
            dirnorm *= 0.5
            loc2 = locmid
            lin2 = linmid
        else:
            raise ValueError("gradient and value seem inconsistent")
        loc = locnext
def check_value_gradient_metric_consistency2(op, loc, tol=1e-8, ntries=100):
    for _ in range(ntries):
        lin = op(Linearization.make_var(loc))
        loc2, lin2 = _get_acceptable_location(op, loc, lin)
        val = lin.val
        dir = loc2 - loc
        locnext = loc2
        dirnorm = dir.norm()
        for i in range(50):
            locmid = loc + 0.5*dir
            linmid = op(Linearization.make_var(locmid))
            dirder = linmid.jac(dir)/dirnorm
            numgrad = (lin2.val-val)/dirnorm
            dgrad = linmid.metric(dir)/dirnorm
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
171
            dgrad2 = (lin2.gradient-lin.gradient)/dirnorm
Martin Reinecke's avatar
Martin Reinecke committed
172
173
            xtol = tol * dirder.norm() / np.sqrt(dirder.size)
            if ((abs(numgrad-dirder) <= xtol).all() and
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
174
                (abs(dgrad-dgrad2) <= xtol).all()):
Martin Reinecke's avatar
Martin Reinecke committed
175
176
177
178
179
180
181
182
                    break
            dir = dir*0.5
            dirnorm *= 0.5
            loc2 = locmid
            lin2 = linmid
        else:
            raise ValueError("gradient and value seem inconsistent")
        loc = locnext
Martin Reinecke's avatar
Martin Reinecke committed
183

Martin Reinecke's avatar
Martin Reinecke committed
184
def check_value_gradient_metric_consistency(E, tol=1e-8, ntries=100):
185
    if isinstance(E, Model):
Martin Reinecke's avatar
Martin Reinecke committed
186
        raise ValueError('Models have no metric, thus it cannot be tested.')
Martin Reinecke's avatar
Martin Reinecke committed
187
188
    for _ in range(ntries):
        E2 = _get_acceptable_energy(E)
189
        val = E.value
Martin Reinecke's avatar
Martin Reinecke committed
190
        dir = E2.position - E.position
Martin Reinecke's avatar
tweak    
Martin Reinecke committed
191
        Enext = E2
Martin Reinecke's avatar
Martin Reinecke committed
192
        dirnorm = dir.norm()
193
        for i in range(50):
194
195
            Emid = E.at(E.position + 0.5*dir)
            dirder = Emid.gradient.vdot(dir)/dirnorm
Martin Reinecke's avatar
Martin Reinecke committed
196
            dgrad = Emid.metric(dir)/dirnorm
Martin Reinecke's avatar
Martin Reinecke committed
197
198
199
            xtol = tol*Emid.gradient_norm
            if abs((E2.value-val)/dirnorm - dirder) < xtol and \
               (abs((E2.gradient-E.gradient)/dirnorm-dgrad) < xtol).all():
200
                break
201
            dir = dir*0.5
Martin Reinecke's avatar
Martin Reinecke committed
202
            dirnorm *= 0.5
203
            E2 = Emid
204
        else:
Martin Reinecke's avatar
Martin Reinecke committed
205
            raise ValueError("gradient, value and metric seem inconsistent")
Martin Reinecke's avatar
tweak    
Martin Reinecke committed
206
        E = Enext