adjust_variances.py 3.73 KB
Newer Older
Philipp Arras's avatar
Cleanup  
Philipp Arras committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
# 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.

from __future__ import absolute_import, division, print_function

from ..compat import *
Philipp Arras's avatar
Philipp Arras committed
22 23 24 25
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
26
from ..operators.energy_operators import Hamiltonian, InverseGammaLikelihood
Philipp Frank's avatar
Philipp Frank committed
27
from ..operators.scaling_operator import ScalingOperator
Philipp Arras's avatar
Philipp Arras committed
28
from ..operators.simple_linear_operators import FieldAdapter
Philipp Frank's avatar
Philipp Frank committed
29

Philipp Arras's avatar
Cleanup  
Philipp Arras committed
30

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

    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
42 43 44 45 46 47 48 49
    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
50 51
    samples : Field, MultiField
        Residual samples of the whole problem
Philipp Frank's avatar
Philipp Frank committed
52 53
    scaling : Float
        Optional rescaling of the Likelihood
54 55
    ic_samp : Controller
        Iteration Controller for Hamiltonian
Philipp Frank's avatar
Philipp Frank committed
56 57 58

    Returns
    -------
59 60
    Hamiltonian
        A Hamiltonian that can be used for further minimization
Philipp Frank's avatar
Philipp Frank committed
61 62
    """

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

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

78
    return Hamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
Philipp Arras's avatar
Philipp Arras committed
79 80


81 82 83 84 85
def do_adjust_variances(position,
                        amplitude_model,
                        minimizer,
                        xi_key='xi',
                        samples=[]):
86

87
    h_space = position[xi_key].domain[0]
Philipp Arras's avatar
Philipp Arras committed
88 89
    pd = PowerDistributor(h_space, amplitude_model.target[0])
    a = pd(amplitude_model)
Philipp Arras's avatar
Philipp Arras committed
90
    xi = FieldAdapter(h_space, xi_key)
Philipp Arras's avatar
Philipp Arras committed
91

92
    ham = make_adjust_variances(a, xi, position, samples=samples)
Philipp Arras's avatar
Philipp Arras committed
93 94

    # Minimize
Philipp Arras's avatar
Philipp Arras committed
95 96
    e = EnergyAdapter(
        position.extract(a.domain), ham, constants=[], want_metric=True)
Philipp Arras's avatar
Philipp Arras committed
97 98 99
    e, _ = minimizer(e)

    # Update position
Philipp Arras's avatar
Philipp Arras committed
100
    s_h_old = (a*xi).force(position)
Philipp Arras's avatar
Philipp Arras committed
101 102 103 104 105 106

    position = position.to_dict()
    position['xi'] = s_h_old/a(e.position)
    position = MultiField.from_dict(position)
    position = MultiField.union([position, e.position])

Philipp Arras's avatar
Philipp Arras committed
107
    s_h_new = (a*xi).force(position)
Philipp Arras's avatar
Philipp Arras committed
108 109 110 111 112 113

    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