adjust_variances.py 4.45 KB
Newer Older
Philipp Arras's avatar
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
Philipp Arras committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Philipp Arras committed
17

Philipp Arras's avatar
Philipp Arras committed
18 19
from ..minimization.energy_adapter import EnergyAdapter
from ..multi_field import MultiField
20 21
from ..operators.energy_operators import (InverseGammaLikelihood,
                                          StandardHamiltonian)
22
from ..operators.scaling_operator import ScalingOperator
23
from ..operators.simple_linear_operators import ducktape
24

Philipp Arras's avatar
Philipp Arras committed
25

26 27 28 29 30 31
def make_adjust_variances_hamiltonian(a,
                                      xi,
                                      position,
                                      samples=[],
                                      scaling=None,
                                      ic_samp=None):
32
    """Creates a Hamiltonian for constant likelihood optimizations.
Philipp Arras's avatar
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.

37 38
    xi is desired to be a Gaussian white Field, thus variations that are
    more easily represented by a should be absorbed in a.
Philipp Arras's avatar
Philipp Arras committed
39

40 41 42
    Parameters
    ----------
    a : Operator
43
        Gives the amplitude when evaluated at position.
44
    xi : Operator
45 46
        Field Adapter selecting a part of position.
        xi is desired to be a Gaussian white Field.
47
    position : Field, MultiField
48
        Contains the initial values for the operators a and xi, to be adjusted
Philipp Arras's avatar
Philipp Arras committed
49
    samples : Field, MultiField
50
        Residual samples of position.
51
    scaling : Float
Philipp Arras's avatar
Philipp Arras committed
52
        Optional rescaling of the Likelihood.
53
    ic_samp : Controller
Philipp Arras's avatar
Philipp Arras committed
54
        Iteration Controller for Hamiltonian.
55 56 57

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

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

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

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


81
def do_adjust_variances(position, A, minimizer, xi_key='xi', samples=[]):
82
    """Adjusts the variance of xi_key to be represented by amplitude_operator.
Martin Reinecke's avatar
Martin Reinecke committed
83

84
    Solves a constant likelihood optimization of the
85 86
    form phi = A * position[xi_key] under the constraint that phi remains
    constant.
87

Martin Reinecke's avatar
Martin Reinecke committed
88
    The field indexed by xi_key is desired to be a Gaussian white Field,
89 90
    thus variations that are more easily represented by A will be absorbed in
    A.
91 92 93 94 95 96

    Parameters
    ----------
    position : Field, MultiField
        Contains the initial values for amplitude_operator and the key xi_key,
        to be adjusted.
97
    A : Operator
98 99 100 101 102 103 104 105 106 107 108 109 110
        Gives the amplitude when evaluated at position.
    minimizer : Minimizer
        Used to solve the optimization problem.
    xi_key : String
        Key of the Field containing undesired variations. This Field is
        contained in position.
    samples : Field, MultiField, optional
        Residual samples of position. If samples are supplied then phi remains
        only approximately constant. Default: [].

    Returns
    -------
    MultiField
111
        The new position after variances have been adjusted.
112
    """
113
    xi = ducktape(None, position.domain, xi_key)
114
    ham = make_adjust_variances_hamiltonian(A, xi, position, samples=samples)
Philipp Arras's avatar
Philipp Arras committed
115 116

    # Minimize
Philipp Arras's avatar
Philipp Arras committed
117
    e = EnergyAdapter(
118
        position.extract(A.domain), ham, constants=[], want_metric=True)
Philipp Arras's avatar
Philipp Arras committed
119 120 121
    e, _ = minimizer(e)

    # Update position
122
    s_h_old = (A*xi).force(position)
Philipp Arras's avatar
Philipp Arras committed
123 124

    position = position.to_dict()
125
    position[xi_key] = s_h_old/A(e.position)
Philipp Arras's avatar
Philipp Arras committed
126
    position = MultiField.from_dict(position)
127
    return MultiField.union([position, e.position])