adjust_variances.py 3.6 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
from ..minimization.energy_adapter import EnergyAdapter
from ..multi_field import MultiField
from ..operators.distributors import PowerDistributor
21
from ..operators.energy_operators import StandardHamiltonian, InverseGammaLikelihood
Philipp Frank's avatar
Philipp Frank committed
22
from ..operators.scaling_operator import ScalingOperator
Martin Reinecke's avatar
Martin Reinecke committed
23
from ..operators.simple_linear_operators import ducktape
Philipp Frank's avatar
Philipp Frank committed
24

Philipp Arras's avatar
Cleanup  
Philipp Arras committed
25

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

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

Philipp Arras's avatar
Philipp Arras committed
37 38
    FIXME xi is white.

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

    Returns
    -------
Martin Reinecke's avatar
Martin Reinecke committed
56
    StandardHamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
57
        A Hamiltonian that can be used for further minimization.
Philipp Frank's avatar
Philipp Frank committed
58 59
    """

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

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

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


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

87
    h_space = position[xi_key].domain[0]
Philipp Arras's avatar
Philipp Arras committed
88 89
    pd = PowerDistributor(h_space, amplitude_operator.target[0])
    a = pd(amplitude_operator)
Martin Reinecke's avatar
Martin Reinecke committed
90
    xi = ducktape(None, position.domain, 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

    position = position.to_dict()
103
    position[xi_key] = s_h_old/a(e.position)
Philipp Arras's avatar
Philipp Arras committed
104 105 106
    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