energy_and_model_tests.py 3.13 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
def _get_acceptable_location(op, loc, lin):
Martin Reinecke's avatar
Martin Reinecke committed
30
    if not np.isfinite(lin.val.sum()):
Martin Reinecke's avatar
Martin Reinecke committed
31
32
33
34
        raise ValueError('Initial value must be finite')
    dir = from_random("normal", loc.domain)
    dirder = lin.jac(dir)
    if dirder.norm() == 0:
Martin Reinecke's avatar
Martin Reinecke committed
35
        dir = dir * (lin.val.norm()*1e-5)
Martin Reinecke's avatar
Martin Reinecke committed
36
    else:
Martin Reinecke's avatar
Martin Reinecke committed
37
        dir = dir * (lin.val.norm()*1e-5/dirder.norm())
Martin Reinecke's avatar
Martin Reinecke committed
38
39
40
41
42
43
44
45
46
47
48
49
50
51
    # 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
52

Martin Reinecke's avatar
Martin Reinecke committed
53
def _check_consistency(op, loc, tol, ntries, do_metric):
Martin Reinecke's avatar
Martin Reinecke committed
54
55
56
    for _ in range(ntries):
        lin = op(Linearization.make_var(loc))
        loc2, lin2 = _get_acceptable_location(op, loc, lin)
Martin Reinecke's avatar
Martin Reinecke committed
57
        dir = loc2-loc
Martin Reinecke's avatar
Martin Reinecke committed
58
59
60
61
62
63
        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
Martin Reinecke's avatar
Martin Reinecke committed
64
            numgrad = (lin2.val-lin.val)/dirnorm
Martin Reinecke's avatar
Martin Reinecke committed
65
            xtol = tol * dirder.norm() / np.sqrt(dirder.size)
Martin Reinecke's avatar
Martin Reinecke committed
66
67
68
69
70
71
            cond = (abs(numgrad-dirder) <= xtol).all()
            if do_metric:
                dgrad = linmid.metric(dir)/dirnorm
                dgrad2 = (lin2.gradient-lin.gradient)/dirnorm
                cond = cond and (abs(dgrad-dgrad2) <= xtol).all()
            if cond:
Martin Reinecke's avatar
Martin Reinecke committed
72
73
74
                break
            dir = dir*0.5
            dirnorm *= 0.5
Martin Reinecke's avatar
Martin Reinecke committed
75
            loc2, lin2 = locmid, linmid
Martin Reinecke's avatar
Martin Reinecke committed
76
77
78
        else:
            raise ValueError("gradient and value seem inconsistent")
        loc = locnext
Martin Reinecke's avatar
Martin Reinecke committed
79
80


Martin Reinecke's avatar
Martin Reinecke committed
81
82
83
84
def check_value_gradient_consistency(op, loc, tol=1e-8, ntries=100):
    _check_consistency(op, loc, tol, ntries, False)


Martin Reinecke's avatar
Martin Reinecke committed
85
def check_value_gradient_metric_consistency(op, loc, tol=1e-8, ntries=100):
Martin Reinecke's avatar
Martin Reinecke committed
86
    _check_consistency(op, loc, tol, ntries, True)