Upload New File

demos/paper_demos/reconstruction.py

+# Program to generate figures of article "Information theory for fields" 
+# by Torsten Ensslin, Annalen der Physik, submitted to special edition 
+# "Physics of Information" in April 2018
+
+import numpy as np
+import nifty4 as ift
+import matplotlib.pyplot as plt
+
+
+if __name__ == "__main__":
+    np.random.seed(20)
+
+   # Set up physical constants
+    nu = 1.        # excitation field level
+    kappa = 10.    # diffusion constant
+    eps = 1e-8     # small number to tame zero mode
+    sigma_n = 0.2  # noise level
+    sigma_n2 = sigma_n**2
+    L = 1.         # Total length of interval or volume the field lives on
+    nprobes = 1000  # Number of probes for uncertainty quantification
+
+    # Define resolution (pixels per dimension)
+    N_pixels = 1024
+
+    # Define two data gaps
+    N1a = int(0.6*N_pixels)
+    N1b = int(0.64*N_pixels)
+    N2a = int(0.67*N_pixels)
+    N2b = int(0.8*N_pixels)
+
+    # Set up derived constants
+    amp = nu/(2*kappa)  # spectral normalization
+    pow_spec = (lambda k: amp / (eps + k**2))
+    lambda2 = 2*kappa*sigma_n2/nu  # resulting correlation length squared
+    lambda1 = np.sqrt(lambda2)
+    pixel_width = L/N_pixels
+    x = np.arange(0,1,pixel_width)
+
+    # Set up the geometry
+    s_domain = ift.RGSpace([N_pixels], distances=pixel_width)
+    h_domain = s_domain.get_default_codomain()
+    p_domain = ift.PowerSpace(h_domain)
+    HT = ift.HarmonicTransformOperator(h_domain, s_domain)
+    aHT = HT.adjoint
+
+    # Create mock signal
+    Phi_h = ift.create_power_operator(h_domain, power_spectrum=pow_spec)
+    phi_h = Phi_h.draw_sample()
+    # remove zero mode as this is part of T_0
+    phi_h.val[0] = 0
+    phi = HT(phi_h)
+
+    # Setting up an exemplary response
+    GeoRem = ift.GeometryRemover(s_domain)
+    GeoAdd = GeoRem.adjoint
+    d_domain = GeoRem.target[0]
+    mask = np.ones(d_domain.shape)
+    mask[N1a:N1b] = 0.
+    mask[N2a:N2b] = 0.
+    mask = ift.Field.from_global_data(d_domain, mask)
+    Mask = ift.DiagonalOperator(mask)
+    R0 = Mask*GeoRem
+    R = R0*HT
+
+    # Linear measurement scenario
+    N = ift.ScalingOperator(sigma_n2, d_domain)  # Noise covariance
+    n = Mask(N.draw_sample())  # seting the noise to zero in masked region
+    d = R(phi_h) + n
+
+    # Wiener filter
+    j = R.adjoint_times(N.inverse_times(d))
+    IC = ift.GradientNormController(name="inverter", iteration_limit=500,
+                                    tol_abs_gradnorm=1e-3)
+    inverter = ift.ConjugateGradient(controller=IC)
+    D = (ift.SandwichOperator(R, N.inverse) + Phi_h.inverse).inverse
+    D = ift.InversionEnabler(D, inverter, approximation=Phi_h)
+    m = HT(D(j))
+
+    # Uncertainty
+    D = ift.SandwichOperator(aHT, D)  # real space propagator
+    Dhat = ift.probe_with_posterior_samples(D.inverse, None, nprobes=nprobes)[1]
+    sig = ift.sqrt(Dhat)
+
+    # Plotting
+    x_mod = x*mask.val+2*(1-mask.val)
+    plt.rcParams["text.usetex"] =True
+    plt.fill_between(x, m.val-sig.val, m.val+sig.val, color='pink',
+                      alpha=None)
+    plt.plot(x,phi.val, label=r"$\varphi$", color='black')
+    plt.scatter(x_mod,d.val, label=r'$d$', s=1, color='blue', alpha=0.5)
+    plt.plot(x,m.val, label=r'$m$', color='red')
+    plt.xlim([0,L])
+    plt.ylim([-1,1])
+    #plt.xlabel('coordinate')
+    #plt.ylabel('field value')
+    plt.title('Wiener filter reconstruction')
+    plt.legend()
+    #plt.show()
+    plt.savefig('Wiener_filter.pdf')
+    plt.close()
+
+
+    # Nonlinear measurement scenario
+    rho_0 = pixel_width
+    R = ift.ScalingOperator(rho_0, d_domain)*R
+    aR = R.adjoint
+
+    def Response(phi_h):
+        return R(aHT(ift.exp(3*HT(phi_h))))
+
+    lam = Response(phi_h)
+    data = ift.Field(d_domain, val=np.random.poisson(lam.val),
+                     dtype=np.float64, copy=True)
+
+    # defining energy, gradient, curvature
+    class Hamiltonian(ift.Energy):
+        def __init__(self, position, inverter):
+            super(Hamiltonian, self).__init__(position=position)
+            self._inverter = inverter
+            lam = Response(position)
+            rho = 3*ift.exp(3*HT(position))
+            Rho = ift.DiagonalOperator(rho)
+            eps = 1e-100  # to catch harmless 0/0 where data is zero
+            weight = (data+eps)/(lam**2+eps)
+            W = ift.DiagonalOperator(weight)
+
+            prior_energy = 0.5*position.vdot(Phi_h.inverse_adjoint_times(position))
+            likel_energy = lam.sum()-data.vdot(ift.log(lam+eps))
+
+            prior_grad = Phi_h.inverse_adjoint_times(position)
+            likel_grad = aHT(rho*HT(aR((lam-data)/(lam+eps))))
+
+            prior_curv = Phi_h.inverse
+    #        j = R.adjoint((data-lam)/(lam+eps))
+    #        tmp = aHT(3*rho*HT(j))
+    #        like1_curv = ift.DiagonalOperator(tmp)  # dangerous term
+            R1 = R0*Rho*HT
+            aR1 = R1.adjoint
+            like2_curv = ift.SandwichOperator(R1, W)  # = aR1*W*R1
+
+            self._curv = prior_curv + like2_curv
+            self._grad = prior_grad + likel_grad
+            self._grad.lock()
+            self._value = prior_energy + likel_energy
+
+        def at(self, position):
+            return self.__class__(position, self._inverter)
+
+        @property
+        def value(self):
+            return self._value
+
+        @property
+        def gradient(self):
+            return self._grad
+
+        @property
+        def curvature(self):
+            return ift.InversionEnabler(self._curv, self._inverter,
+                                        approximation=Phi_h.inverse)
+
+
+    # initial guess
+    psi0 = ift.Field.full(h_domain, 1e-7)
+    energy = Hamiltonian(psi0, inverter)
+    IC1 = ift.GradientNormController(name="IC1", iteration_limit=200,
+                                     tol_abs_gradnorm=1e-4)
+    minimizer = ift.RelaxedNewton(IC1)
+    energy = minimizer(energy)[0]
+
+    var = ift.probe_with_posterior_samples(energy.curvature, HT, nprobes)[1]
+    sig = ift.sqrt(var)
+
+    m = HT(energy.position)
+    phi = HT(phi_h)
+    plt.rcParams["text.usetex"] = True
+    c1 = np.exp(3*(m.val-sig.val))
+    c2 = np.exp(3*(m.val+sig.val))
+    plt.fill_between(x, c1, c2, color='pink',
+                      alpha=None)
+    plt.plot(x,np.exp(3*phi.val), label=r"$e^{3\varphi}$", color='black')
+    plt.scatter(x_mod,data.val, label=r'$d$', s=1, color='blue', alpha=0.5)
+    plt.plot(x,np.exp(3*m.val), label=r'$e^{3\varphi_\mathrm{cl}}$', color='red')
+    plt.xlim([0,L])
+    plt.ylim([-0.1,7.5])
+    #plt.xlabel('coordinate')
+    #plt.ylabel('field value')
+    plt.title('Poisson log-normal reconstruction')
+    plt.legend()
+    #plt.show()
+    plt.savefig('Poisson.pdf')
+    exit()
\ No newline at end of file