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