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