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