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