adjust_variances.py 4.84 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

27 28 29 30 31 32
def make_adjust_variances_hamiltonian(a,
                                      xi,
                                      position,
                                      samples=[],
                                      scaling=None,
                                      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.

38 39
    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
40

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

    Returns
    -------
Martin Reinecke's avatar
Martin Reinecke committed
59
    StandardHamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
60
        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
fixups  
Philipp Frank committed
64
    d = (d.conjugate()*d).real
Philipp Frank's avatar
Philipp Frank committed
65
    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

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


82
def do_adjust_variances(position,
Philipp Arras's avatar
Philipp Arras committed
83
                        amplitude_operator,
84 85 86
                        minimizer,
                        xi_key='xi',
                        samples=[]):
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
    """Adjusts the variance of xi_key to be represented by amplitude_operator.
    
    Solves a constant likelihood optimization of the
    form phi = amplitude_operator * position[xi_key] under the constraint that
    phi remains constant.

    The field indexed by xi_key is desired to be a Gaussian white Field, 
    thus variations that are more easily represented by amplitude_operator 
    will be absorbed in amplitude_operator.

    Parameters
    ----------
    position : Field, MultiField
        Contains the initial values for amplitude_operator and the key xi_key,
        to be adjusted.
    amplitude_operator : Operator
        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
        The new position after variances were adjusted.
    """

119

120
    h_space = position[xi_key].domain[0]
Philipp Arras's avatar
Philipp Arras committed
121 122
    pd = PowerDistributor(h_space, amplitude_operator.target[0])
    a = pd(amplitude_operator)
Martin Reinecke's avatar
Martin Reinecke committed
123
    xi = ducktape(None, position.domain, xi_key)
Philipp Arras's avatar
Philipp Arras committed
124

125
    ham = make_adjust_variances_hamiltonian(a, xi, position, samples=samples)
Philipp Arras's avatar
Philipp Arras committed
126 127

    # Minimize
Philipp Arras's avatar
Philipp Arras committed
128 129
    e = EnergyAdapter(
        position.extract(a.domain), ham, constants=[], want_metric=True)
Philipp Arras's avatar
Philipp Arras committed
130 131 132
    e, _ = minimizer(e)

    # Update position
Philipp Arras's avatar
Philipp Arras committed
133
    s_h_old = (a*xi).force(position)
Philipp Arras's avatar
Philipp Arras committed
134 135

    position = position.to_dict()
136
    position[xi_key] = s_h_old/a(e.position)
Philipp Arras's avatar
Philipp Arras committed
137 138 139 140
    position = MultiField.from_dict(position)
    position = MultiField.union([position, e.position])

    return position