critical filter kind of working, still some issues

demos/critical_filtering.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_true, t_real):
+    m = fft.adjoint_times(m)
+    m_data = m.val.get_full_data().real
+    t_data = t.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(y=t_data), go.Scatter(y=t_true_data),
+             go.Scatter(y=t_real_data)], filename="t.html")
+
+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 = True,
+                         distribution_strategy=distribution_strategy)
+
+    # Choosing the prior correlation structure and defining correlation operator
+    pow_spec = (lambda k: (.05 / (k + 1) ** 3))
+    S = create_power_operator(h_space, power_spectrum=pow_spec,
+                              distribution_strategy=distribution_strategy)
+
+    # Drawing a sample sh from the prior distribution in harmonic space
+    sp = Field(p_space,  val=lambda z: pow_spec(z)**(1./2),
+               distribution_strategy=distribution_strategy)
+    sh = sp.power_synthesize(real_signal=True)
+
+
+    # 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
+
+    #   Choosing nonlinearity
+
+    # log-normal model:
+
+    function = exp
+    derivative = exp
+
+    # tan-normal model
+
+    # def function(x):
+    #     return 0.5 * tanh(x) + 0.5
+    # def derivative(x):
+    #     return 0.5*(1 - tanh(x)**2)
+
+    # no nonlinearity, Wiener Filter
+
+    # def function(x):
+    #     return x
+    # def derivative(x):
+    #     return 1
+
+    # small quadratic pertubarion
+
+    # def function(x):
+    #     return 0.5*x**2 + x
+    # def derivative(x):
+    #     return x + 1
+
+    # def function(x):
+    #     return 0.9*x**4 +0.2*x**2 + x
+    # def derivative(x):
+    #     return 0.9*4*x**3 + 0.4*x +1
+    #
+
+    #Adding a harmonic transformation to the instrument
+    R = NonlinearResponse(fft, Instrument, function, derivative)
+    noise = .01
+    N = DiagonalOperator(s_space, diagonal=noise, bare=True)
+    n = Field.from_random(domain=s_space,
+                          random_type='normal',
+                          std=sqrt(noise),
+                          mean=0)
+
+    # Creating the mock data
+    d = R(sh) + n
+    realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"])**2)
+    d_data = d.val.get_full_data().real
+    if rank == 0:
+        pl.plot([go.Heatmap(z=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=5,
+                              callback=convergence_measure)
+    # minimizer2 = RelaxedNewton(convergence_tolerance=0,
+    #                           convergence_level=1,
+    #                           iteration_limit=2,
+    #                           callback=convergence_measure)
+    #
+    # minimizer1 = VL_BFGS(convergence_tolerance=0,
+    #                    iteration_limit=5,
+    #                    callback=convergence_measure,
+    #                    max_history_length=3)
+
+
+
+    # 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=-8)
+    # t0 = log(sp.copy()**2)
+    S0 = create_power_operator(h_space, power_spectrum=exp(t0),
+                               distribution_strategy=distribution_strategy)
+
+
+
+    for i in range(100):
+        S0 = create_power_operator(h_space, power_spectrum=exp(t0),
+                              distribution_strategy=distribution_strategy)
+
+        # Initializing the  nonlinear Wiener Filter energy
+        map_energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
+        # Minimization with chosen minimizer
+        (map_energy, convergence) = minimizer1(map_energy)
+        # 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=10., samples=3)
+        (power_energy, convergence) = minimizer1(power_energy)
+        # Setting new power spectrum
+        t0 = power_energy.position
+        plot_parameters(m0,t0,log(sp**2),realized_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/nonlinear_wiener_filter.py

@@ -54,7 +54,7 @@ if __name__ == "__main__":
     distribution_strategy = 'not'
 
     # Set up position space
-    s_space = RGSpace([256,256])
+    s_space = RGSpace([100,100])
     # s_space = HPSpace(32)
 
     # Define harmonic transformation and associated harmonic space
@@ -65,7 +65,7 @@ if __name__ == "__main__":
     p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
 
     # Choosing the prior correlation structure and defining correlation operator
-    pow_spec = (lambda k: (0.42 / (k + 1) ** 3))
+    pow_spec = (lambda k: (1. / (k + 1) ** 3))
     S = create_power_operator(h_space, power_spectrum=pow_spec,
                               distribution_strategy=distribution_strategy)
 
@@ -81,23 +81,23 @@ if __name__ == "__main__":
 #    Instrument._diagonal.val[200:400, 200:400] = 0
 
     #   Choosing nonlinearity
-
+    #
     # function = exp
     # derivative = exp
     def function(x):
         return 0.5 * tanh(x) + 0.5
     def derivative(x):
-        return 0.5*(1 - tanh(x)**2)
+        return 0.5*(1. - tanh(x)**2)
     # def function(x):
-    #     return x
+    #     return 20*x
     # def derivative(x):
-    #     return 1
-    #
+    #     return 20.
+
     # def function(x):
-    #     return 0.5*x**2 + x
+    #     return 0.1*0.5*x**2 + x
     # def derivative(x):
-    #     return x + 1
-    #
+    #     return 0.1*x + 1
+
     # def function(x):
     #     return 0.9*x**4 +0.2*x**2 + x
     # def derivative(x):
@@ -105,14 +105,14 @@ if __name__ == "__main__":
     #
 
     #Adding a harmonic transformation to the instrument
-    R = NonlinearResponse(fft, Instrument, function, derivative)
-    noise = 0.01
+
+    noise = 10.
     N = DiagonalOperator(s_space, diagonal=noise, bare=True)
     n = Field.from_random(domain=s_space,
                           random_type='normal',
                           std=sqrt(noise),
                           mean=0)
-
+    R = NonlinearResponse(fft, Instrument, function, derivative)
     # Creating the mock data
     d = R(sh) + n
 
@@ -128,7 +128,7 @@ if __name__ == "__main__":
 #
     minimizer = RelaxedNewton(convergence_tolerance=0,
                               convergence_level=1,
-                              iteration_limit=10,
+                              iteration_limit=3,
                               callback=convergence_measure)
 
     # minimizer = VL_BFGS(convergence_tolerance=0,
@@ -136,10 +136,11 @@ if __name__ == "__main__":
     #                    callback=convergence_measure,
     #                    max_history_length=3)
     #
-    #
+
 
     # Setting starting position
-    m0 = Field.from_random(random_type="normal",domain = h_space, std=0.001)
+    flat_power = Field(p_space,val=10e-8)
+    m0 = flat_power.power_synthesize(real_signal=True)
 
     # Initializing the Wiener Filter energy
     energy = NonlinearWienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S)

nifty/library/energy_library/critical_power_energy.py

@@ -22,19 +22,20 @@ class CriticalPowerEnergy(Energy):
         The prior signal covariance in harmonic space.
     """
 
-    def __init__(self, position, m, D=None, alpha =1.0, q=0, sigma=0, w=None, samples=3):
+    def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0., w=None, samples=3):
         super(CriticalPowerEnergy, self).__init__(position = position)
         self.m = m
         self.D = D
         self.samples = samples
         self.sigma = sigma
-        self.alpha = alpha
-        self.q = q
-        self.T = SmoothnessOperator(domain=self.position.domain, sigma=self.sigma)
-        self.rho = self.position.domain.rho
-        if w is None:
+        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.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 + 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,
@@ -45,7 +46,7 @@ class CriticalPowerEnergy(Energy):
 
     @property
     def value(self):
-        energy = self.theta.sum()
+        energy = exp(-self.position).dot(self.q + self.w / 2.)
         energy += self.position.dot(self.alpha - 1 + self.rho / 2.)
         energy += 0.5 * self.position.dot(self.T(self.position))
         return energy.real
@@ -55,6 +56,7 @@ class CriticalPowerEnergy(Energy):
         gradient = - self.theta
         gradient += self.alpha - 1 + self.rho / 2.
         gradient += self.T(self.position)
+        gradient.val[0] = 0.
         return gradient
 
     @property
@@ -63,15 +65,19 @@ class CriticalPowerEnergy(Energy):
         return curvature
 
     def _calculate_w(self, m, D, samples):
-        w = Field(domain=self.position.domain, val=0)
+        w = Field(domain=self.position.domain, val=0. , dtype=m.dtype)
         if D is not None:
             for i in range(samples):
                 posterior_sample = generate_posterior_sample(m, D)
-                projected_sample =posterior_sample.project_power(domain=self.position.domain)
+                projected_sample = posterior_sample.project_power(
+                    logarithmic=self.position.domain[0].config["logarithmic"],
+                        decompose_power=False)
                 w += projected_sample
             w /= float(samples)
         else:
-            w = m.project_power(domain=self.position.domain)
+            w = m.project_power(
+                    logarithmic=self.position.domain[0].config["logarithmic"],
+                        decompose_power=False)
 
         return w
 

nifty/library/operator_library/critical_power_curvature.py

@@ -7,8 +7,8 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
 
         self.theta = theta
         self.T = T
-        if preconditioner is None:
-            preconditioner = self.T.times
+        # if preconditioner is None:
+        #     preconditioner = self.T.times
         self._domain = self.T.domain
         super(CriticalPowerCurvature, self).__init__(inverter=inverter,
                                                  preconditioner=preconditioner)

nifty/library/operator_library/laplace_operator.py

@@ -34,8 +34,8 @@ class LaplaceOperator(EndomorphicOperator):
         return False
     def _times(self, t, spaces):
         if t.val.distribution_strategy != 'not':
-            self.logger.warn("distribution_strategy should be 'not' but is %s"
-                             % t.val.distribution_strategy)
+            # self.logger.warn("distribution_strategy should be 'not' but is %s"
+            #                  % t.val.distribution_strategy)
             array = t.val.get_full_data()
         else:
             array = t.val.get_local_data(copy=False)
@@ -46,8 +46,8 @@ class LaplaceOperator(EndomorphicOperator):
         return Field(self.domain, val=ret)
     def _adjoint_times(self, t, spaces):
         if t.val.distribution_strategy != 'not':
-            self.logger.warn("distribution_strategy should be 'not' but is %s"
-                             % t.val.distribution_strategy)
+            # self.logger.warn("distribution_strategy should be 'not' but is %s"
+            #                  % t.val.distribution_strategy)
             array = t.val.get_full_data()
         else:
             array = t.val.get_local_data(copy=False)

nifty/library/operator_library/nonlinear_wiener_filter_curvature.py

@@ -9,8 +9,8 @@ class NonlinearWienerFilterCurvature(InvertibleOperatorMixin, EndomorphicOperato
         self.N = N
         self.S = S
         self.position = position
-        if preconditioner is None:
-            preconditioner = self.S.times
+        # if preconditioner is None:
+        #     preconditioner = self.S.times
         self._domain = self.S.domain
         super(NonlinearWienerFilterCurvature, self).__init__(inverter=inverter,
                                                  preconditioner=preconditioner)

nifty/sugar.py

@@ -25,13 +25,15 @@ from nifty import PowerSpace,\
 __all__ = ['create_power_operator']
 
 
-def create_power_operator(domain, power_spectrum, dtype=None,
+def create_power_operator(domain, power_spectrum, power_domain=None, dtype=None,
                           distribution_strategy='not'):
     if not domain.harmonic:
         fft = FFTOperator(domain)
         domain = fft.target[0]
-
-    power_domain = PowerSpace(domain,
+    if isinstance(power_spectrum, Field):
+        power_domain = power_spectrum.domain
+    elif power_domain is None:
+        power_domain = PowerSpace(domain,
                               distribution_strategy=distribution_strategy)
 
     fp = Field(power_domain,
@@ -49,15 +51,17 @@ def generate_posterior_sample(mean, covariance):
     S = covariance.S
     R = covariance.R
     N = covariance.N
-    power = S.diagonal()
+    power = sqrt(S.diagonal().power_analyze())
+    mock_signal = power.power_synthesize(real_signal=True)
+
+
     noise = N.diagonal().val
-    mock_signal = Field.from_random(random_type="normal", domain=S.domain,
-                                    std = sqrt(power), dtype = power.dtype)
+
     mock_noise = Field.from_random(random_type="normal", domain=N.domain,
                                    std = sqrt(noise), dtype = noise.dtype)
-    mock_data = R(mock_signal) + mock_noise
+    mock_data = R.derived_times(mock_signal, mean) + mock_noise
 
-    mock_j = R.adjoint_times(N.inverse_times(mock_data))
+    mock_j = R.derived_adjoint_times(N.inverse_times(mock_data), mean)
     mock_m = covariance.inverse_times(mock_j)
     sample = mock_signal - mock_m + mean
     return sample