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

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)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
34
35
36
        #self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma)
        self.Laplace = LaplaceOperator(self.position.domain[0])
        self.roughness = self.Laplace(self.position)
37
38
39
        self.rho = self.position.domain[0].rho
        self.w = w
        if self.w is None:
40
            self.w = self._calculate_w(self.m, self.D, self.samples)
41
        self.theta = exp(-self.position) * (self.q + self.w / 2.)
42
43

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

    @property
    def value(self):
52
        energy = exp(-self.position).dot(self.q + self.w / 2.)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
53
54
        energy += self.position.dot(self.alpha - 1. + self.rho / 2.)
        energy += 0.5 * self.roughness.dot(self.roughness) / self.sigma**2
55
56
57
58
        return energy.real

    @property
    def gradient(self):
59
        gradient = - self.theta
Jakob Knollmueller's avatar
Jakob Knollmueller committed
60
61
        gradient += self.alpha - 1. + self.rho / 2.
        gradient += self.Laplace(self.roughness) / self.sigma **2
62
63
64
65
        return gradient

    @property
    def curvature(self):
Jakob Knollmueller's avatar
Jakob Knollmueller committed
66
67
        curvature = CriticalPowerCurvature(theta=self.theta, Laplace = self.Laplace,
                                           sigma = self.sigma)
68
69
70
        return curvature

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

88
        return w
89
90
91