adjust_variances.py 3.58 KB
Newer Older
Philipp Arras's avatar
Cleanup  
Philipp Arras 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
17

Philipp Arras's avatar
Philipp Arras committed
18 19 20 21
from ..minimization.energy_adapter import EnergyAdapter
from ..multi_domain import MultiDomain
from ..multi_field import MultiField
from ..operators.distributors import PowerDistributor
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
22
from ..operators.energy_operators import Hamiltonian, InverseGammaLikelihood
Philipp Frank's avatar
Philipp Frank committed
23
from ..operators.scaling_operator import ScalingOperator
Martin Reinecke's avatar
Martin Reinecke committed
24
from ..operators.simple_linear_operators import ducktape
Philipp Frank's avatar
Philipp Frank committed
25

Philipp Arras's avatar
Cleanup  
Philipp Arras committed
26

Philipp Arras's avatar
Philipp Arras committed
27 28 29 30 31
def make_adjust_variances(a,
                          xi,
                          position,
                          samples=[],
                          scaling=None,
32
                          ic_samp=None):
33
    """Creates a Hamiltonian for constant likelihood optimizations.
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
34 35 36 37

    Constructs a Hamiltonian to solve constant likelihood optimizations of the
    form phi = a * xi under the constraint that phi remains constant.

Philipp Frank's avatar
Philipp Frank committed
38 39 40 41 42 43 44 45
    Parameters
    ----------
    a : Operator
        Operator which gives the amplitude when evaluated at a position
    xi : Operator
        Operator which gives the excitation when evaluated at a position
    postion : Field, MultiField
        Position of the whole problem
Philipp Arras's avatar
Philipp Arras committed
46 47
    samples : Field, MultiField
        Residual samples of the whole problem
Philipp Frank's avatar
Philipp Frank committed
48 49
    scaling : Float
        Optional rescaling of the Likelihood
50 51
    ic_samp : Controller
        Iteration Controller for Hamiltonian
Philipp Frank's avatar
Philipp Frank committed
52 53 54

    Returns
    -------
55 56
    Hamiltonian
        A Hamiltonian that can be used for further minimization
Philipp Frank's avatar
Philipp Frank committed
57 58
    """

Philipp Arras's avatar
Cleanup  
Philipp Arras committed
59
    d = a*xi
Philipp Frank's avatar
Philipp Frank committed
60 61
    d = (d.conjugate()*d).real
    n = len(samples)
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
62
    if n > 0:
Philipp Frank's avatar
Philipp Frank committed
63 64
        d_eval = 0.
        for i in range(n):
65
            d_eval = d_eval + d.force(position + samples[i])
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
66
        d_eval = d_eval/n
Philipp Frank's avatar
Philipp Frank committed
67
    else:
68
        d_eval = d.force(position)
Philipp Frank's avatar
Philipp Frank committed
69 70 71

    x = (a.conjugate()*a).real
    if scaling is not None:
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
72
        x = ScalingOperator(scaling, x.target)(x)
Philipp Frank's avatar
Philipp Frank committed
73

74
    return Hamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
Philipp Arras's avatar
Philipp Arras committed
75 76


77
def do_adjust_variances(position,
Philipp Arras's avatar
Philipp Arras committed
78
                        amplitude_operator,
79 80 81
                        minimizer,
                        xi_key='xi',
                        samples=[]):
82

83
    h_space = position[xi_key].domain[0]
Philipp Arras's avatar
Philipp Arras committed
84 85
    pd = PowerDistributor(h_space, amplitude_operator.target[0])
    a = pd(amplitude_operator)
Martin Reinecke's avatar
Martin Reinecke committed
86
    xi = ducktape(None, position.domain, xi_key)
Philipp Arras's avatar
Philipp Arras committed
87

88
    ham = make_adjust_variances(a, xi, position, samples=samples)
Philipp Arras's avatar
Philipp Arras committed
89 90

    # Minimize
Philipp Arras's avatar
Philipp Arras committed
91 92
    e = EnergyAdapter(
        position.extract(a.domain), ham, constants=[], want_metric=True)
Philipp Arras's avatar
Philipp Arras committed
93 94 95
    e, _ = minimizer(e)

    # Update position
Philipp Arras's avatar
Philipp Arras committed
96
    s_h_old = (a*xi).force(position)
Philipp Arras's avatar
Philipp Arras committed
97 98

    position = position.to_dict()
99
    position[xi_key] = s_h_old/a(e.position)
Philipp Arras's avatar
Philipp Arras committed
100 101 102
    position = MultiField.from_dict(position)
    position = MultiField.union([position, e.position])

Philipp Arras's avatar
Philipp Arras committed
103
    s_h_new = (a*xi).force(position)
Philipp Arras's avatar
Philipp Arras committed
104 105 106 107 108 109

    import numpy as np
    # TODO Move this into the tests
    np.testing.assert_allclose(s_h_new.to_global_data(),
                               s_h_old.to_global_data())
    return position