critical filtering seems to work

.gitignore

@@ -7,6 +7,7 @@
 *.egg
 *.egg-info
 *~
+*.DS_Store
 *.csv
 .git/
 .svn/

demos/critical_filtering.py

@@ -8,7 +8,7 @@ from mpi4py import MPI
 comm = MPI.COMM_WORLD
 rank = comm.rank
 
-np.random.seed(232)
+np.random.seed(42)
 
 
 def plot_parameters(m,t,t_true, t_real, t_d):
@@ -60,8 +60,8 @@ if __name__ == "__main__":
     h_space = fft.target[0]
 
     # Setting up power space
-    p_space = PowerSpace(h_space, logarithmic=False,
-                         distribution_strategy=distribution_strategy, nbin=128)
+    p_space = PowerSpace(h_space, logarithmic=True,
+                         distribution_strategy=distribution_strategy, nbin=120)
 
     # Choosing the prior correlation structure and defining correlation operator
     pow_spec = (lambda k: (.05 / (k + 1) ** 2))
@@ -83,7 +83,7 @@ if __name__ == "__main__":
     #Adding a harmonic transformation to the instrument
     R = AdjointFFTResponse(fft, Instrument)
 
-    noise = 1.
+    noise = 5.
     N = DiagonalOperator(s_space, diagonal=noise, bare=True)
     n = Field.from_random(domain=s_space,
                           random_type='normal',
@@ -95,7 +95,7 @@ if __name__ == "__main__":
 
     realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"],
                                           nbin=p_space.config["nbin"])**2)
-    data_power  = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
+    data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
                                           nbin=p_space.config["nbin"])**2)
     d_data = d.val.get_full_data().real
     if rank == 0:
@@ -110,24 +110,24 @@ if __name__ == "__main__":
 
     minimizer1 = RelaxedNewton(convergence_tolerance=0,
                               convergence_level=1,
-                              iteration_limit=2,
+                              iteration_limit=3,
                               callback=convergence_measure)
     minimizer2 = VL_BFGS(convergence_tolerance=0,
                        iteration_limit=50,
                        callback=convergence_measure,
-                       max_history_length=3)
+                       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 = data_power
+    # t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
+    t0 = data_power- 1.
     S0 = create_power_operator(h_space, power_spectrum=exp(t0),
                                distribution_strategy=distribution_strategy)
 
 
-    for i in range(100):
+    for i in range(500):
         S0 = create_power_operator(h_space, power_spectrum=exp(t0),
                               distribution_strategy=distribution_strategy)
 
@@ -140,14 +140,13 @@ 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=10., samples=3)
-        if i > 0:
-            (power_energy, convergence) = minimizer1(power_energy)
-        else:
-            (power_energy, convergence) = minimizer2(power_energy)
+        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=1000., samples=1)
+
+        (power_energy, convergence) = minimizer1(power_energy)
+
         # Setting new power spectrum
-        t0 = power_energy.position
-        t0.val[-1] = t0.val[-2]
+        t0.val  = power_energy.position.val.real
+        # t0.val[-1] = t0.val[-2]
         # Plotting current estimate
         plot_parameters(m0,t0,log(sp**2),realized_power, data_power)
 

demos/laplace_testing.py

@@ -58,23 +58,23 @@ if __name__ == "__main__":
 
     # Setting up power space
     p_space = PowerSpace(h_space, logarithmic=False,
-                         distribution_strategy=distribution_strategy, nbin=70)
+                         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)
-    T = SmoothnessOperator(p_space,sigma=1.)
+    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=1000,
+                       iteration_limit=10,
                        callback=convergence_measure,
-                       max_history_length=3)
+                       max_history_length=10)
 
 
     def explicify(op, domain):
@@ -91,7 +91,7 @@ if __name__ == "__main__":
     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=log(p_space.kindex)[1:], y=data[1:])], filename="t.html")
+    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))
@@ -100,6 +100,6 @@ if __name__ == "__main__":
     aa = Field(p_space, val=p_space.kindex.copy())
     aa.val[0] = 1
 
-    print lap(log(aa)).val
+    print lap(log(aa)**2).val
     print "######################"
     print test_energy.position.val
\ No newline at end of file

nifty/library/energy_library/critical_power_energy.py

 from nifty.energies.energy import Energy
 from nifty.library.operator_library import CriticalPowerCurvature,\
-                                            LaplaceOperator
+                                            SmoothnessOperator
 
 from nifty.sugar import generate_posterior_sample
 from nifty import Field, exp
@@ -31,14 +31,12 @@ class CriticalPowerEnergy(Energy):
         self.sigma = sigma
         self.alpha = Field(self.position.domain, val=alpha)
         self.q = Field(self.position.domain, val=q)
-        #self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma)
-        self.Laplace = LaplaceOperator(self.position.domain[0])
-        self.roughness = self.Laplace(self.position)
+        self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma)
         self.rho = self.position.domain[0].rho
         self.w = w
         if self.w is None:
             self.w = self._calculate_w(self.m, self.D, self.samples)
-        self.theta = exp(-self.position) * (self.q + self.w / 2.)
+        self.theta = (exp(-self.position) * (self.q + self.w / 2.))
 
     def at(self, position):
         return self.__class__(position, self.m, D=self.D,
@@ -49,22 +47,22 @@ class CriticalPowerEnergy(Energy):
 
     @property
     def value(self):
-        energy = exp(-self.position).dot(self.q + self.w / 2.)
-        energy += self.position.dot(self.alpha - 1. + self.rho / 2.)
-        energy += 0.5 * self.roughness.dot(self.roughness) / self.sigma**2
+        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))
         return energy.real
 
     @property
     def gradient(self):
-        gradient = - self.theta
-        gradient += self.alpha - 1. + self.rho / 2.
-        gradient += self.Laplace(self.roughness) / self.sigma **2
+        gradient = - self.theta.weight(-1)
+        gradient += (self.alpha - 1. + self.rho / 2.).weight(-1)
+        gradient +=  self.T(self.position)
+        gradient.val = gradient.val.real
         return gradient
 
     @property
     def curvature(self):
-        curvature = CriticalPowerCurvature(theta=self.theta, Laplace = self.Laplace,
-                                           sigma = self.sigma)
+        curvature = CriticalPowerCurvature(theta=self.theta.weight(-1), T=self.T)
         return curvature
 
     def _calculate_w(self, m, D, samples):
@@ -81,6 +79,7 @@ class CriticalPowerEnergy(Energy):
         else:
             w = m.power_analyze(
                     logarithmic=self.position.domain[0].config["logarithmic"],
+                        nbin=self.position.domain[0].config["nbin"],
                         decompose_power=False)
             w **= 2
             w *= self.rho

nifty/library/operator_library/critical_power_curvature.py

@@ -3,11 +3,10 @@ from nifty.operators import EndomorphicOperator,\
 
 
 class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
-    def __init__(self, theta, Laplace, sigma, inverter=None, preconditioner=None):
+    def __init__(self, theta, T, inverter=None, preconditioner=None):
 
         self.theta = theta
-        self.Laplace = Laplace
-        self.sigma = sigma
+        self.T = T
         # if preconditioner is None:
         #     preconditioner = self.T.times
         self._domain = self.theta.domain
@@ -28,5 +27,4 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
     # ---Added properties and methods---
 
     def _times(self, x, spaces):
-        return self.Laplace.adjoint_times(self.Laplace(x)) / self.sigma ** 2 \
-               + self.theta * x
+        return self.T(x) + self.theta * x

nifty/library/operator_library/smoothness_operator.py

@@ -6,7 +6,7 @@ from laplace_operator import LaplaceOperator
 
 class SmoothnessOperator(EndomorphicOperator):
 
-    def __init__(self, domain, sigma=1.,
+    def __init__(self, domain, sigma=1.,logarithmic = True,
                  default_spaces=None):
 
         super(SmoothnessOperator, self).__init__(default_spaces)
@@ -21,7 +21,7 @@ class SmoothnessOperator(EndomorphicOperator):
 
         self._sigma = sigma
 
-        self._Laplace = LaplaceOperator(domain=self._domain[0])
+        self._Laplace = LaplaceOperator(domain=self._domain[0], logarithmic=logarithmic)
 
     """
     SmoothnessOperator