adjust_variances.py 3.67 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
Martin Reinecke's avatar
Martin Reinecke committed
21 22
from ..operators.energy_operators import (StandardHamiltonian,
                                          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 Arras's avatar
Philipp Arras committed
38 39
    FIXME xi is white.

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
76 77
    return StandardHamiltonian(InverseGammaLikelihood(d_eval/2.)(x),
                               ic_samp=ic_samp)
Philipp Arras's avatar
Philipp Arras committed
78 79


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

89
    h_space = position[xi_key].domain[0]
Philipp Arras's avatar
Philipp Arras committed
90 91
    pd = PowerDistributor(h_space, amplitude_operator.target[0])
    a = pd(amplitude_operator)
Martin Reinecke's avatar
Martin Reinecke committed
92
    xi = ducktape(None, position.domain, xi_key)
Philipp Arras's avatar
Philipp Arras committed
93

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

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

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

    position = position.to_dict()
105
    position[xi_key] = s_h_old/a(e.position)
Philipp Arras's avatar
Philipp Arras committed
106 107 108
    position = MultiField.from_dict(position)
    position = MultiField.union([position, e.position])

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

    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