energy_and_model_tests.py 3.72 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
Martin Reinecke committed
23
from ..linearization import Linearization
24

Martin Reinecke's avatar
Martin Reinecke committed
25
__all__ = ["check_value_gradient_consistency",
Martin Reinecke's avatar
Martin Reinecke committed
26
           "check_value_gradient_metric_consistency"]
27

Philipp Arras's avatar
Philipp Arras committed
28

Martin Reinecke's avatar
Martin Reinecke committed
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
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
53

Martin Reinecke's avatar
Martin Reinecke committed
54
def check_value_gradient_consistency(op, loc, tol=1e-8, ntries=100):
Martin Reinecke's avatar
Martin Reinecke committed
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
    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
Martin Reinecke's avatar
Martin Reinecke committed
77
78
79


def check_value_gradient_metric_consistency(op, loc, tol=1e-8, ntries=100):
Martin Reinecke's avatar
Martin Reinecke committed
80
81
82
83
84
85
86
87
88
89
90
91
92
    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
93
            dgrad2 = (lin2.gradient-lin.gradient)/dirnorm
Martin Reinecke's avatar
Martin Reinecke committed
94
95
            xtol = tol * dirder.norm() / np.sqrt(dirder.size)
            if ((abs(numgrad-dirder) <= xtol).all() and
Martin Reinecke's avatar
Martin Reinecke committed
96
97
                    (abs(dgrad-dgrad2) <= xtol).all()):
                break
Martin Reinecke's avatar
Martin Reinecke committed
98
99
100
101
102
103
104
            dir = dir*0.5
            dirnorm *= 0.5
            loc2 = locmid
            lin2 = linmid
        else:
            raise ValueError("gradient and value seem inconsistent")
        loc = locnext