critical_power_energy.py 3.28 KB
Newer Older
1
from nifty.energies.energy import Energy
2
from nifty.library.operator_library import CriticalPowerCurvature,\
3
                                            SmoothnessOperator
Jakob Knollmueller's avatar
Jakob Knollmueller committed
4

5
from nifty.sugar import generate_posterior_sample
6
from nifty import Field, exp
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

class CriticalPowerEnergy(Energy):
    """The Energy for the Gaussian lognormal case.

    It describes the situation of linear measurement  of a
    lognormal signal with Gaussian noise and Gaussain signal prior.

    Parameters
    ----------
    d : Field,
        the data.
    R : Operator,
        The nonlinear response operator, describtion of the measurement process.
    N : EndomorphicOperator,
        The noise covariance in data space.
    S : EndomorphicOperator,
        The prior signal covariance in harmonic space.
    """

26
    def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0., w=None, samples=3):
27
        super(CriticalPowerEnergy, self).__init__(position = position)
28
        self.m = m
29
30
31
        self.D = D
        self.samples = samples
        self.sigma = sigma
32
33
        self.alpha = Field(self.position.domain, val=alpha)
        self.q = Field(self.position.domain, val=q)
34
        self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma)
35
36
37
        self.rho = self.position.domain[0].rho
        self.w = w
        if self.w is None:
38
            self.w = self._calculate_w(self.m, self.D, self.samples)
39
        self.theta = (exp(-self.position) * (self.q + self.w / 2.))
40
41

    def at(self, position):
42
43
44
45
46
        return self.__class__(position, self.m, D=self.D,
                              alpha =self.alpha,
                              q=self.q,
                              sigma=self.sigma, w=self.w,
                              samples=self.samples)
47
48
49

    @property
    def value(self):
50
51
52
        energy = exp(-self.position).dot(self.q + self.w / 2., bare= True)
        energy += self.position.dot(self.alpha - 1. + self.rho / 2., bare=True)
        energy += 0.5 * self.position.dot(self.T(self.position))
53
54
55
56
        return energy.real

    @property
    def gradient(self):
57
58
59
60
        gradient = - self.theta.weight(-1)
        gradient += (self.alpha - 1. + self.rho / 2.).weight(-1)
        gradient +=  self.T(self.position)
        gradient.val = gradient.val.real
61
62
63
64
        return gradient

    @property
    def curvature(self):
65
        curvature = CriticalPowerCurvature(theta=self.theta.weight(-1), T=self.T)
66
67
68
        return curvature

    def _calculate_w(self, m, D, samples):
69
        w = Field(domain=self.position.domain, val=0. , dtype=m.dtype)
70
71
72
        if D is not None:
            for i in range(samples):
                posterior_sample = generate_posterior_sample(m, D)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
73
                projected_sample = posterior_sample.power_analyze(
74
                    logarithmic=self.position.domain[0].config["logarithmic"],
Jakob Knollmueller's avatar
Jakob Knollmueller committed
75
                    nbin= self.position.domain[0].config["nbin"],
76
                        decompose_power=False)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
77
                w += (projected_sample **2) * self.rho
78
79
            w /= float(samples)
        else:
Jakob Knollmueller's avatar
Jakob Knollmueller committed
80
            w = m.power_analyze(
81
                    logarithmic=self.position.domain[0].config["logarithmic"],
82
                        nbin=self.position.domain[0].config["nbin"],
83
                        decompose_power=False)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
84
85
            w **= 2
            w *= self.rho
86

87
        return w
88
89
90