added memorization of the residual in WienerFilterEnergy

demos/critical_filtering.py

@@ -11,18 +11,16 @@ rank = comm.rank
 np.random.seed(42)
 
 
-def plot_parameters(m,t,t_true, t_real, t_d):
+def plot_parameters(m,t,p, p_d):
 
     x = log(t.domain[0].kindex)
     m = fft.adjoint_times(m)
-    m_data = m.val.get_full_data().real
-    t_data = t.val.get_full_data().real
-    t_d_data = t_d.val.get_full_data().real
-    t_true_data = t_true.val.get_full_data().real
-    t_real_data = t_real.val.get_full_data().real
-    pl.plot([go.Heatmap(z=m_data)], filename='map.html')
-    pl.plot([go.Scatter(x=x,y=t_data), go.Scatter(x=x ,y=t_true_data),
-             go.Scatter(x=x, y=t_real_data), go.Scatter(x=x, y=t_d_data)], filename="t.html")
+    m = m.val.get_full_data().real
+    t = t.val.get_full_data().real
+    p = p.val.get_full_data().real
+    p_d = p_d.val.get_full_data().real
+    pl.plot([go.Heatmap(z=m)], filename='map.html')
+    pl.plot([go.Scatter(x=x,y=t), go.Scatter(x=x ,y=p), go.Scatter(x=x, y=p_d)], filename="t.html")
 
 
 class AdjointFFTResponse(LinearOperator):
@@ -63,11 +61,11 @@ if __name__ == "__main__":
     h_space = fft.target[0]
 
     # Setting up power space
-    p_space = PowerSpace(h_space, logarithmic=False,
+    p_space = PowerSpace(h_space, logarithmic=True,
                          distribution_strategy=distribution_strategy)#, nbin=5)
 
     # Choosing the prior correlation structure and defining correlation operator
-    p_spec = (lambda k: (.05 / (k + 1) ** 3))
+    p_spec = (lambda k: (.5 / (k + 1) ** 3))
     S = create_power_operator(h_space, power_spectrum=p_spec,
                               distribution_strategy=distribution_strategy)
 
@@ -87,7 +85,7 @@ if __name__ == "__main__":
     #Adding a harmonic transformation to the instrument
     R = AdjointFFTResponse(fft, Instrument)
 
-    noise = .1
+    noise = 1.
     N = DiagonalOperator(s_space, diagonal=noise, bare=True)
     n = Field.from_random(domain=s_space,
                           random_type='normal',
@@ -112,8 +110,8 @@ if __name__ == "__main__":
         print (x, iteration)
 
 
-    minimizer1 = RelaxedNewton(convergence_tolerance=0,
-                              convergence_level=1,
+    minimizer1 = RelaxedNewton(convergence_tolerance=10e-2,
+                              convergence_level=2,
                               iteration_limit=3,
                               callback=convergence_measure)
     minimizer2 = VL_BFGS(convergence_tolerance=0,
@@ -125,8 +123,8 @@ if __name__ == "__main__":
     flat_power = Field(p_space,val=10e-8)
     m0 = flat_power.power_synthesize(real_signal=True)
 
-    # t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
-    t0 = data_power- 1.
+    t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
+
 
 
     for i in range(500):
@@ -142,42 +140,15 @@ if __name__ == "__main__":
         m0 = map_energy.position
         D0 = map_energy.curvature
         # Initializing the power energy with updated parameters
-        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=100., samples=5)
+        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=100., samples=3)
 
         (power_energy, convergence) = minimizer1(power_energy)
 
+
         # Setting new power spectrum
         t0.val  = power_energy.position.val.real
 
         # Plotting current estimate
-        plot_parameters(m0,t0,log(sp),realized_power, data_power)
-
-    # Transforming fields to position space for plotting
-
-    ss = fft.adjoint_times(sh)
-    m = fft.adjoint_times(map_energy.position)
-
+        plot_parameters(m0,t0,log(sp), data_power)
 
-    # Plotting
-
-    d_data = d.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=d_data)], filename='data.html')
-
-    tt_data = power_energy.position.val.get_full_data().real
-    t_data = log(sp**2).val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Scatter(y=t_data),go.Scatter(y=tt_data)], filename="t.html")
-    ss_data = ss.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=ss_data)], filename='ss.html')
-
-    sh_data = sh.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=sh_data)], filename='sh.html')
-
-
-    m_data = m.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=m_data)], filename='map.html')
 

demos/laplace_testing.py

-from nifty import *
-
-import numpy as np
-from nifty import Field,\
-                  EndomorphicOperator,\
-                  PowerSpace
-import plotly.offline as pl
-import plotly.graph_objs as go
-
-import numpy as np
-from nifty import Field, \
-    EndomorphicOperator, \
-    PowerSpace
-
-
-class TestEnergy(Energy):
-    def __init__(self, position, Op):
-        super(TestEnergy, self).__init__(position)
-        self.Op = Op
-
-    def at(self, position):
-        return self.__class__(position=position, Op=self.Op)
-
-    @property
-    def value(self):
-        return 0.5 * self.position.dot(self.Op(self.position))
-
-    @property
-    def gradient(self):
-        return self.Op(self.position)
-
-    @property
-    def curvature(self):
-        curv = CurvOp(self.Op)
-        return curv
-
-class CurvOp(InvertibleOperatorMixin, EndomorphicOperator):
-    def __init__(self, Op,inverter=None, preconditioner=None):
-        self.Op = Op
-        self._domain = Op.domain
-        super(CurvOp, self).__init__(inverter=inverter,
-                                                     preconditioner=preconditioner)
-
-    def _times(self, x, spaces):
-        return self.Op(x)
-
-if __name__ == "__main__":
-
-    distribution_strategy = 'not'
-
-    # Set up position space
-    s_space = RGSpace([128,128])
-    # s_space = HPSpace(32)
-
-    # Define harmonic transformation and associated harmonic space
-    fft = FFTOperator(s_space)
-    h_space = fft.target[0]
-
-    # Setting up power space
-    p_space = PowerSpace(h_space, logarithmic=False,
-                         distribution_strategy=distribution_strategy, nbin=200)
-
-    # Choosing the prior correlation structure and defining correlation operator
-    pow_spec = (lambda k: (.05 / (k + 1) ** 2))
-    # t = Field(p_space, val=pow_spec)
-    t= Field.from_random("normal", domain=p_space)
-    lap = LaplaceOperator(p_space, logarithmic=True)
-    T = SmoothnessOperator(p_space,sigma=1., logarithmic=True)
-    test_energy = TestEnergy(t,T)
-
-    def convergence_measure(a_energy, iteration): # returns current energy
-        x = a_energy.value
-        print (x, iteration)
-    minimizer1 = VL_BFGS(convergence_tolerance=0,
-                       iteration_limit=10,
-                       callback=convergence_measure,
-                       max_history_length=10)
-
-
-    def explicify(op, domain):
-        space = domain
-        d = space.dim
-        res = np.zeros((d, d))
-        for i in range(d):
-            x = np.zeros(d)
-            x[i] = 1.
-            f = Field(space, val=x)
-            res[:, i] = op(f).val
-        return res
-    A = explicify(lap.times, p_space)
-    B = explicify(lap.adjoint_times, p_space)
-    test_energy,convergence = minimizer1(test_energy)
-    data = test_energy.position.val.get_full_data()
-    pl.plot([go.Scatter(x=(p_space.kindex)[1:], y=data[1:])], filename="t.html")
-    tt = Field.from_random("normal", domain=t.domain)
-    print "adjointness"
-    print t.dot(lap(tt))
-    print tt.dot(lap.adjoint_times(t))
-    print "log kindex"
-    aa = Field(p_space, val=p_space.kindex.copy())
-    aa.val[0] = 1
-
-    print lap(log(aa)**2).val
-    print "######################"
-    print test_energy.position.val
\ No newline at end of file

demos/train_example.py

-
-from nifty import *
-import plotly.offline as pl
-import plotly.graph_objs as go
-
-from mpi4py import MPI
-comm = MPI.COMM_WORLD
-rank = comm.rank
-
-np.random.seed(62)
-
-class NonlinearResponse(LinearOperator):
-    def __init__(self, FFT, Instrument, function, derivative, default_spaces=None):
-        super(NonlinearResponse, self).__init__(default_spaces)
-        self._domain = FFT.target
-        self._target = Instrument.target
-        self.function = function
-        self.derivative = derivative
-        self.I = Instrument
-        self.FFT = FFT
-
-
-    def _times(self, x, spaces=None):
-        return self.I(self.function(self.FFT.adjoint_times(x)))
-
-    def _adjoint_times(self, x, spaces=None):
-        return self.FFT(self.function(self.I.adjoint_times(x)))
-
-    def derived_times(self, x, position):
-        position_derivative = self.derivative(self.FFT.adjoint_times(position))
-        return self.I(position_derivative * self.FFT.adjoint_times(x))
-
-    def derived_adjoint_times(self, x, position):
-        position_derivative = self.derivative(self.FFT.adjoint_times(position))
-        return self.FFT(position_derivative * self.I.adjoint_times(x))
-
-    @property
-    def domain(self):
-        return self._domain
-
-    @property
-    def target(self):
-        return self._target
-
-    @property
-    def unitary(self):
-        return False
-
-def plot_parameters(m,t, t_real):
-    m = fft.adjoint_times(m)
-    x =  log(t.domain[0].kindex[1:])
-    m_data = m.val.get_full_data().real
-    t_data = t.val.get_full_data().real
-    t_real_data = t_real.val.get_full_data().real
-    pl.plot([go.Scatter(y=m_data)], filename='map.html')
-    pl.plot([go.Scatter(x = x,y=t_data),
-             go.Scatter(x= x, y=t_real_data[1:])], filename="t.html")
-class AdjointFFTResponse(LinearOperator):
-    def __init__(self, FFT, R, default_spaces=None):
-        super(AdjointFFTResponse, self).__init__(default_spaces)
-        self._domain = FFT.target
-        self._target = R.target
-        self.R = R
-        self.FFT = FFT
-
-    def _times(self, x, spaces=None):
-        return self.R(self.FFT.adjoint_times(x))
-
-    def _adjoint_times(self, x, spaces=None):
-        return self.FFT(self.R.adjoint_times(x))
-    @property
-    def domain(self):
-        return self._domain
-
-    @property
-    def target(self):
-        return self._target
-
-    @property
-    def unitary(self):
-        return False
-
-if __name__ == "__main__":
-
-    distribution_strategy = 'not'
-    full_data = np.genfromtxt("train_data.csv", delimiter = ',')
-    d = full_data.T[2]
-    d[0] = 0.
-    d -= d.mean()
-    d[0] = 0.
-    # Set up position space
-    s_space = RGSpace([len(d)])
-    # s_space = HPSpace(32)
-
-    d = Field(s_space, val=d)
-
-    # Define harmonic transformation and associated harmonic space
-    fft = FFTOperator(s_space)
-    h_space = fft.target[0]
-    p_space = PowerSpace(h_space, logarithmic=False,
-                         distribution_strategy=distribution_strategy)#, nbin=50)
-
-    # Choosing the measurement instrument
-#    Instrument = SmoothingOperator(s_space, sigma=0.01)
-    Instrument = DiagonalOperator(s_space, diagonal=1.)
-#    Instrument._diagonal.val[200:400, 200:400] = 0
-
-
-    #Adding a harmonic transformation to the instrument
-    R = AdjointFFTResponse(fft, Instrument)
-    noise = .1
-    N = DiagonalOperator(s_space, diagonal=noise, bare=True)
-
-    d_data = d.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Scatter(y=d_data)], filename='data.html')
-
-    # Choosing the minimization strategy
-
-    def convergence_measure(a_energy, iteration): # returns current energy
-        x = a_energy.value
-        print (x, iteration)
-
-    # minimizer1 = SteepestDescent(convergence_tolerance=0,
-    #                            iteration_limit=50,
-    #                            callback=convergence_measure)
-
-    minimizer1 = RelaxedNewton(convergence_tolerance=0,
-                              convergence_level=1,
-                              iteration_limit=2,
-                              callback=convergence_measure)
-    minimizer2 = RelaxedNewton(convergence_tolerance=0,
-                              convergence_level=1,
-                              iteration_limit=30,
-                              callback=convergence_measure)
-
-    minimizer3 = VL_BFGS(convergence_tolerance=0,
-                       iteration_limit=50,
-                       callback=convergence_measure,
-                       max_history_length=10)
-
-
-
-    # Setting starting position
-    flat_power = Field(p_space,val=10e-8)
-    m0 = flat_power.power_synthesize(real_signal=True)
-
-    t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
-    t0 = Field(p_space,val=-13.)
-    # t0 = log(sp.copy()**2)
-    S0 = create_power_operator(h_space, power_spectrum=exp(t0),
-                               distribution_strategy=distribution_strategy)
-
-    data_power  = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
-                                          nbin=p_space.config["nbin"]))
-    for i in range(500):
-        S0 = create_power_operator(h_space, power_spectrum=exp(t0),
-                                   distribution_strategy=distribution_strategy)
-
-        # Initializing the  nonlinear Wiener Filter energy
-        map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
-        # Minimization with chosen minimizer
-        map_energy = map_energy.analytic_solution()
-
-        # Updating parameters for correlation structure reconstruction
-        m0 = map_energy.position
-        D0 = map_energy.curvature
-        # Initializing the power energy with updated parameters
-        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=100000., samples=20)
-
-        (power_energy, convergence) = minimizer3(power_energy)
-
-        # Setting new power spectrum
-        t0.val = power_energy.position.val.real
-        # t0.val[-1] = t0.val[-2]
-        # Plotting current estimate
-        plot_parameters(m0, t0, data_power)
-
-    # Transforming fields to position space for plotting
-
-    ss = fft.adjoint_times(sh)
-    m = fft.adjoint_times(map_energy.position)
-
-
-
-    # Plotting
-
-    d_data = d.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=d_data)], filename='data.html')
-
-    tt_data = power_energy.position.val.get_full_data().real
-    t_data = log(sp**2).val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Scatter(y=t_data),go.Scatter(y=tt_data)], filename="t.html")
-    ss_data = ss.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=ss_data)], filename='ss.html')
-
-    sh_data = sh.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=sh_data)], filename='sh.html')
-
-
-    m_data = m.val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=m_data)], filename='map.html')
-
-    f_m_data = function(m).val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=f_m_data)], filename='f_map.html')
-    f_ss_data = function(ss).val.get_full_data().real
-    if rank == 0:
-        pl.plot([go.Heatmap(z=f_ss_data)], filename='f_ss.html')

demos/wiener_filter_advanced.py

@@ -42,7 +42,7 @@ if __name__ == "__main__":
     distribution_strategy = 'not'
 
     # Set up position space
-    s_space = RGSpace([512,512])
+    s_space = RGSpace([128,128])
     # s_space = HPSpace(32)
 
     # Define harmonic transformation and associated harmonic space
@@ -108,9 +108,29 @@ if __name__ == "__main__":
     # Initializing the Wiener Filter energy
     energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S)
 
+
     # Solving the problem analytically
     solution = energy.analytic_solution()
+    sample_variance = Field(sh.power_analyze(logarithmic=False).domain,val=0. + 0j)
+    sample_mean = Field(sh.domain,val=0. + 0j)
+    n_samples = 200
+    for i in range(n_samples):
+        sample = sugar.generate_posterior_sample(solution.position,solution.curvature)
+        sample_variance += sample.power_analyze(logarithmic=False)
+        sample_mean += sample
+    sample_variance = sample_variance/n_samples - (sample_mean/n_samples).power_analyze(logarithmic=False)
+
+    D = HarmonicPropagatorOperator(S=S, N=N, R=Instrument)
+
+
+    class DiagonalProber(DiagonalProberMixin, Prober):
+        pass
+
 
+    diagProber = DiagonalProber(domain=sh.domain)
+    diagProber(D)
+    diag = diagProber.diagonal
+    pr_diag = diag.power_analyze(logarithmic=False)**0.5
     # Solving the problem with chosen minimization strategy
     (energy, convergence) = minimizer(energy)
 

nifty/library/energy_library/wiener_filter_energy.py

 from nifty.energies.energy import Energy
+from nifty.energies.memoization import memo
 from nifty.library.operator_library import WienerFilterCurvature
 
 class WienerFilterEnergy(Energy):
@@ -33,16 +34,17 @@ class WienerFilterEnergy(Energy):
 
     @property
     def value(self):
+        residual = self._residual()
         energy = 0.5 * self.position.dot(self.S.inverse_times(self.position))
-        energy += 0.5 * (self.d - self.R(self.position)).dot(
-            self.N.inverse_times(self.d - self.R(self.position)))
+        energy += 0.5 * (residual).dot(self.N.inverse_times(residual))
         return energy.real
 
     @property
     def gradient(self):
+        residual = self._residual()
         gradient = self.S.inverse_times(self.position)
         gradient -= self.R.adjoint_times(
-                    self.N.inverse_times(self.d - self.R(self.position)))
+                    self.N.inverse_times(residual))
         return gradient
 
     @property
@@ -50,6 +52,12 @@ class WienerFilterEnergy(Energy):
         curvature = WienerFilterCurvature(R=self.R, N=self.N, S=self.S)
         return curvature
 
+    @memo
+    def _residual(self):
+        residual = self.d - self.R(self.position)
+        return residual
+
+
     def analytic_solution(self):
         D_inverse = self.curvature
         j = self.R.adjoint_times(self.N.inverse_times(self.d))