Fixed problems with the transposed operations of the forward model which did

not represent the correct matrizes (hopefully, now they do).

>>> PACKAGE
forwardmodel, kernel, phasemapper:
        b_0 is now again inherent part of the kernel.
projector:
        Corrected .jac_T_dot, now delivers product for the correct matrix.
>>> SCRIPTS
simple_reconstruction:
        Does simple compliance tests for 1 and 2 projections and tests a simple
        solver via scipy.sparse.linalg.cg for the reconstruction.

pyramid/forwardmodel.py

@@ -33,10 +33,9 @@ class ForwardModel:
         assert isinstance(kernel, Kernel), 'A Kernel object has to be provided!'
         self._kernel = kernel
 
-    def __init__(self, projectors, kernel, b_0):
+    def __init__(self, projectors, kernel):
         # TODO: Docstring!
         self.kernel = kernel
-        self.b_0 = b_0
         self.a = kernel.a
         self.dim_uv = kernel.dim_uv
         self.projectors = projectors
@@ -45,7 +44,7 @@ class ForwardModel:
         # TODO: Docstring!
 #        print 'FWD Model - __call__ -  input: ', len(x)
         result = [self.kernel.jac_dot(projector.jac_dot(x)) for projector in self.projectors]
-        result = self.b_0 * np.reshape(result, -1)
+        result = np.reshape(result, -1)
 #        print 'FWD Model - __call__ -  output:', len(result)
         return result
 

pyramid/kernel.py

@@ -76,7 +76,7 @@ class Kernel(object):
 
     '''# TODO: Can be used for several PhaseMappers via the fft arguments or via calling!
     
-    def __init__(self, a, dim_uv, numcore=True, geometry='disc'):
+    def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
         '''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
 
         Parameters
@@ -109,7 +109,7 @@ class Kernel(object):
         self.numcore = numcore
         self.geometry = geometry
         # Calculate kernel (single pixel phase):
-        coeff = -a**2 / (2*PHI_0)
+        coeff = -b_0 * a**2 / (2*PHI_0)
         v_dim, u_dim = dim_uv
         u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
         v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)

pyramid/phasemapper.py

@@ -46,7 +46,7 @@ class PMAdapterFM(PhaseMapper):
         assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
         self.a = a
         self.projector = projector
-        self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, geometry), b_0)
+        self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, b_0, geometry))
 
     def __call__(self, mag_data):
         assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
@@ -221,7 +221,7 @@ class PMConvolve(PhaseMapper):
         self.projector = projector
         self.b_0 = b_0
         self.threshold = threshold
-        self.kernel = Kernel(a, projector.dim_uv, geometry)
+        self.kernel = Kernel(a, projector.dim_uv, b_0, geometry)
 
     def __call__(self, mag_data):
         # Docstring!
@@ -235,7 +235,7 @@ class PMConvolve(PhaseMapper):
         u_phase = np.fft.irfftn(u_mag_fft * kernel.u_fft, kernel.dim_fft)[kernel.slice_fft].copy()
         v_phase = np.fft.irfftn(v_mag_fft * kernel.v_fft, kernel.dim_fft)[kernel.slice_fft].copy()
         # Return the result:
-        return PhaseMap(self.a, self.b_0*(u_phase-v_phase))
+        return PhaseMap(self.a, u_phase-v_phase)
 
 
 class PMReal(PhaseMapper):
@@ -271,7 +271,7 @@ class PMReal(PhaseMapper):
         self.projector = projector
         self.b_0 = b_0
         self.threshold = threshold
-        self.kernel = Kernel(a, projector.dim_uv, geometry)
+        self.kernel = Kernel(a, projector.dim_uv, b_0, geometry)
         self.numcore = numcore
 
     def __call__(self, mag_data):
@@ -281,8 +281,8 @@ class PMReal(PhaseMapper):
         threshold = self.threshold
         u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+dim_uv)
         # Create kernel (lookup-tables for the phase of one pixel):
-        u_phi = self.b_0 * self.kernel.u
-        v_phi = self.b_0 * self.kernel.v
+        u_phi = self.kernel.u
+        v_phi = self.kernel.v
         # Calculation of the phase:
         phase = np.zeros(dim_uv)
         if self.numcore:

pyramid/projector.py

@@ -85,18 +85,18 @@ class Projector(object):
         self.log.info('Calling _vector_field_projection_T')
         size_2d, size_3d = self.size_2d, self.size_3d
         result = np.zeros(3*size_3d)
-        # Go over all possible component projections (z, y, x) to (u, v):
-        if self.coeff[0][0] != 0:  # x to u
+        # Go over all possible component projections (u, v) to (z, y, x):
+        if self.coeff[0][0] != 0:  # u to x
             result[:size_3d] += self.coeff[0][0] * self.weight.T.dot(vector[:size_2d])
-        if self.coeff[0][1] != 0:  # y to u
-            result[:size_3d] += self.coeff[0][1] * self.weight.T.dot(vector[size_2d:])
-        if self.coeff[0][2] != 0:  # z to u
-            result[size_3d:2*size_3d] += self.coeff[0][2] * self.weight.T.dot(vector[:size_2d])
-        if self.coeff[1][0] != 0:  # x to v
-            result[size_3d:2*size_3d] += self.coeff[1][0] * self.weight.T.dot(vector[size_2d:])
-        if self.coeff[1][1] != 0:  # y to v
-            result[2*size_3d:] += self.coeff[1][1] * self.weight.T.dot(vector[:size_2d])
-        if self.coeff[1][2] != 0:  # z to v
+        if self.coeff[0][1] != 0:  # u to y
+            result[size_3d:2*size_3d] += self.coeff[0][1] * self.weight.T.dot(vector[:size_2d])
+        if self.coeff[0][2] != 0:  # u to z
+            result[2*size_3d:] += self.coeff[0][2] * self.weight.T.dot(vector[:size_2d])
+        if self.coeff[1][0] != 0:  # v to x
+            result[:size_3d] += self.coeff[1][0] * self.weight.T.dot(vector[size_2d:])
+        if self.coeff[1][1] != 0:  # v to y
+            result[size_3d:2*size_3d] += self.coeff[1][1] * self.weight.T.dot(vector[size_2d:])
+        if self.coeff[1][2] != 0:  # v to z
             result[2*size_3d:] += self.coeff[1][2] * self.weight.T.dot(vector[size_2d:])
         return result
 

scripts/simple_reconstruction.py

@@ -13,25 +13,167 @@ from pyramid.magdata import MagData
 from pyramid.projector import YTiltProjector
 from pyramid.phasemapper import PMConvolve
 from pyramid.datacollection import DataCollection
-import pyramid.optimizer as opt
 
 from pyramid.kernel import Kernel
 from pyramid.forwardmodel import ForwardModel
 from pyramid.costfunction import Costfunction
 
+from scipy.sparse.linalg import cg
+
+
+###################################################################################################
+print('--Compliance test for one projection')
 
 a = 1.
 b_0 = 1000.
-dim = (3, 3, 3)
+dim = (2, 2, 2)
 count = 1
 
+magnitude = np.zeros((3,)+dim)
+magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+
+mag_data = MagData(a, magnitude)
+
+tilts = np.linspace(0, 2*pi, num=count, endpoint=False)
+projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
+phase_maps = [pm(mag_data) for pm in phasemappers]
+
+dim_uv = dim[1:3]
+
+data_collection = DataCollection(a, dim_uv, b_0)
+
+[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+
+data = data_collection
+y = data.phase_vec
+kern = Kernel(data.a, data.dim_uv, data.b_0)
+F = ForwardModel(data.projectors, kern)
+C = Costfunction(y, F)
+
+size_2d = np.prod(dim[1]*dim[2])
+size_3d = np.prod(dim)
+
+P = np.array([data.projectors[0](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+K = np.array([kern(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+F_mult = K.dot(P)
+F_direct = np.array([F(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+np.testing.assert_almost_equal(F_direct, F_mult)
+
+P_jac = np.array([data.projectors[0].jac_dot(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+K_jac = np.array([kern.jac_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+F_jac_mult = K.dot(P)
+F_jac_direct = np.array([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+np.testing.assert_almost_equal(F_jac_direct, F_jac_mult)
+np.testing.assert_almost_equal(F_direct, F_jac_direct)
+
+P_t = np.array([data.projectors[0].jac_T_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+K_t = np.array([kern.jac_T_dot(np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
+F_t_mult = P_t.dot(K_t)
+F_t_direct = np.array([F.jac_T_dot(None, np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
+
+P_t_ref = P.T
+K_t_ref = K.T
+F_t_ref = F_mult.T
+np.testing.assert_almost_equal(P_t, P_t_ref)
+np.testing.assert_almost_equal(K_t, K_t_ref)
+np.testing.assert_almost_equal(F_t_mult, F_t_ref)
+np.testing.assert_almost_equal(F_t_direct, F_t_ref)
+
 ###################################################################################################
-print('--Generating input phase_maps')
+print('--Compliance test for two projections')
 
+a = 1.
+b_0 = 1000.
+dim = (2, 2, 2)
+count = 2
 
 magnitude = np.zeros((3,)+dim)
 magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
 
+mag_data = MagData(a, magnitude)
+
+tilts = np.linspace(0, 2*pi, num=count, endpoint=False)
+projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
+phase_maps = [pm(mag_data) for pm in phasemappers]
+
+dim_uv = dim[1:3]
+
+data_collection = DataCollection(a, dim_uv, b_0)
+
+[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+
+data = data_collection
+y = data.phase_vec
+kern = Kernel(data.a, data.dim_uv, data.b_0)
+F = ForwardModel(data.projectors, kern)
+C = Costfunction(y, F)
+
+size_2d = np.prod(dim[1]*dim[2])
+size_3d = np.prod(dim)
+
+P0 = np.array([data.projectors[0](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+P1 = np.array([data.projectors[1](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+P = np.vstack((P0, P1))
+K = np.array([kern(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+F_mult0 = K.dot(P0)
+F_mult1 = K.dot(P1)
+F_mult = np.vstack((F_mult0, F_mult1))
+F_direct = np.array([F(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+np.testing.assert_almost_equal(F_direct, F_mult)
+
+P_jac0 = np.array([data.projectors[0].jac_dot(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+P_jac1 = np.array([data.projectors[1].jac_dot(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+P = np.vstack((P0, P1))
+K_jac = np.array([kern.jac_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+F_jac_mult0 = K.dot(P_jac0)
+F_jac_mult1 = K.dot(P_jac1)
+F_jac_mult = np.vstack((F_jac_mult0, F_jac_mult1))
+F_jac_direct = np.array([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+np.testing.assert_almost_equal(F_jac_direct, F_jac_mult)
+np.testing.assert_almost_equal(F_direct, F_jac_direct)
+
+P_t0 = np.array([data.projectors[0].jac_T_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+P_t1 = np.array([data.projectors[1].jac_T_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+P_t = np.hstack((P_t0, P_t1))
+K_t = np.array([kern.jac_T_dot(np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
+F_t_mult0 = P_t0.dot(K_t)
+F_t_mult1 = P_t1.dot(K_t)
+F_t_mult = np.hstack((F_t_mult0, F_t_mult1))
+F_t_direct = np.array([F.jac_T_dot(None, np.eye(count*size_2d)[:, i]) for i in range(count*size_2d)]).T
+
+P_t_ref = P.T
+K_t_ref = K.T
+F_t_ref = F_mult.T
+np.testing.assert_almost_equal(P_t, P_t_ref)
+np.testing.assert_almost_equal(K_t, K_t_ref)
+np.testing.assert_almost_equal(F_t_mult, F_t_ref)
+np.testing.assert_almost_equal(F_t_direct, F_t_ref)
+
+
+
+
+
+
+
+
+###################################################################################################
+print('--STARTING RECONSTRUCTION')
+
+
+
+###################################################################################################
+print('--Generating input phase_maps')
+
+a = 1.
+b_0 = 1000.
+dim = (9, 9, 9)
+count = 8
+
+magnitude = np.zeros((3,)+dim)
+magnitude[0, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+
 mag_data = MagData(a, magnitude)
 mag_data.quiver_plot3d()
 
@@ -40,32 +182,67 @@ projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
 phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
 phase_maps = [pm(mag_data) for pm in phasemappers]
 
-#[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
-#                         for i, phase_map in enumerate(phase_maps)]
-phase_maps[0].display_phase()
+[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
+                         for i, phase_map in enumerate(phase_maps)]
 
 ###################################################################################################
 print('--Setting up data collection')
 
 dim_uv = dim[1:3]
 
+size_2d = np.prod(dim_uv)
+size_3d = np.prod(dim)
+
+
 data_collection = DataCollection(a, dim_uv, b_0)
 
 [data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
 
+data = data_collection
+y = data.phase_vec
+kern = Kernel(data.a, data.dim_uv, data.b_0)
+F = ForwardModel(data.projectors, kern)
+C = Costfunction(y, F)
+
 ###################################################################################################
-print('--Test optimizer')
+print('--Test simple solver')
 
-first_guess = MagData(a, np.zeros((3,)+dim))
+M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+lam = 1*10. ** -10
+MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
 
-first_guess.magnitude[1, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1
-#first_guess.magnitude[0, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = -1
+A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
+
+b = F.jac_T_dot(None, y)
+
+b_test = np.asarray((M.T.dot(y)).T)[0]
+
+x_f = cg(A, b)[0]
+
+mag_data_rec = MagData(a, np.zeros((3,)+dim))
+
+mag_data_rec.mag_vec = x_f
+
+mag_data_rec.quiver_plot3d()
+
+phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
+[phase_map.display_phase(title=u'Tilt series (rec.) $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
+                         for i, phase_map in enumerate(phase_maps_rec)]
 
-first_guess.quiver_plot3d()
 
-phase_guess = PMConvolve(first_guess.a, projectors[0], b_0)(first_guess)
 
-phase_guess.display_phase()
+
+#
+#first_guess = MagData(a, np.zeros((3,)+dim))
+#
+#first_guess.magnitude[1, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1
+#first_guess.magnitude[0, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = -1
+#
+#first_guess.quiver_plot3d()
+#
+#phase_guess = PMConvolve(first_guess.a, projectors[0], b_0)(first_guess)
+#
+#phase_guess.display_phase()
 
 #mag_opt = opt.optimize_cg(data_collection, first_guess)
 #
@@ -76,42 +253,42 @@ phase_guess.display_phase()
 #phase_opt.display_phase()
 
 
-###################################################################################################
-print('--Further testing')
-
-data = data_collection
-mag_0 = first_guess
-x_0 = first_guess.mag_vec
-y = data.phase_vec
-kern = Kernel(data.a, data.dim_uv)
-F = ForwardModel(data.projectors, kern, data.b_0)
-C = Costfunction(y, F)
-
-size_3d = np.prod(dim)
-
-cost = C(first_guess.mag_vec)
-cost_grad = C.jac(first_guess.mag_vec)
-cost_grad_del_x = cost_grad.reshape((3, 3, 3, 3))[0, ...]
-cost_grad_del_y = cost_grad.reshape((3, 3, 3, 3))[1, ...]
-cost_grad_del_z = cost_grad.reshape((3, 3, 3, 3))[2, ...]
-
-
-
-
-x_t = np.asmatrix(mag_data.mag_vec).T
-y = np.asmatrix(y).T
-K = np.array([F.jac_dot(x_0, np.eye(81)[:, i]) for i in range(81)]).T
-K = np.asmatrix(K)
-lam = 10. ** -10
-KTK = K.T * K + lam * np.asmatrix(np.eye(81))
-print lam, 
-#print pylab.cond(KTK),
-x_f = KTK.I * K.T * y
-print (np.asarray(K * x_f - y) ** 2).sum(),
-print (np.asarray(K * x_t - y) ** 2).sum()
-print x_f
-x_rec = np.asarray(x_f).reshape(3,3,3,3)
-print x_rec[0, ...]
+####################################################################################################
+#print('--Further testing')
+#
+#data = data_collection
+#mag_0 = first_guess
+#x_0 = first_guess.mag_vec
+#y = data.phase_vec
+#kern = Kernel(data.a, data.dim_uv)
+#F = ForwardModel(data.projectors, kern, data.b_0)
+#C = Costfunction(y, F)
+#
+#size_3d = np.prod(dim)
+#
+#cost = C(first_guess.mag_vec)
+#cost_grad = C.jac(first_guess.mag_vec)
+#cost_grad_del_x = cost_grad.reshape((3, 3, 3, 3))[0, ...]
+#cost_grad_del_y = cost_grad.reshape((3, 3, 3, 3))[1, ...]
+#cost_grad_del_z = cost_grad.reshape((3, 3, 3, 3))[2, ...]
+#
+#
+#
+#
+#x_t = np.asmatrix(mag_data.mag_vec).T
+#y = np.asmatrix(y).T
+#K = np.array([F.jac_dot(x_0, np.eye(81)[:, i]) for i in range(81)]).T
+#K = np.asmatrix(K)
+#lam = 10. ** -10
+#KTK = K.T * K + lam * np.asmatrix(np.eye(81))
+#print lam, 
+##print pylab.cond(KTK),
+#x_f = KTK.I * K.T * y
+#print (np.asarray(K * x_f - y) ** 2).sum(),
+#print (np.asarray(K * x_t - y) ** 2).sum()
+#print x_f
+#x_rec = np.asarray(x_f).reshape(3,3,3,3)
+#print x_rec[0, ...]
 #KTK = K.T * K
 #u,s,v = pylab.svd(KTK, full_matrices=1)
 #si = np.zeros_like(s)
@@ -123,48 +300,34 @@ print x_rec[0, ...]
 #
 #x_f = KTKI * K.T * y
 #print x_f
-#
-#
-
-
-###################################################################################################
-print('--Compliance test')
-
 
-# test F
 
-# test F.jac
 
-# test F.jac_dot
-
-#F_evaluate = F(mag_data.mag_vec)
-#F_jacobi = F.jac_dot(None, mag_data.mag_vec)
-#np.testing.assert_equal(F_evaluate, F_jacobi)
-#
-#F_reverse = F.jac_T_dot(None, F_jacobi)
-#np.testing.assert_equal(F_reverse, mag_data.mag_vec)
 
 
 
-from scipy.sparse.linalg import cg
 
-K = np.asmatrix([F.jac_dot(x_0, np.eye(81)[:, i]) for i in range(81)]).T
-lam = 10. ** -10
-KTK = K.T * K + lam * np.asmatrix(np.eye(81))
 
-A = KTK#np.array([F(np.eye(81)[:, i]) for i in range(81)])
-
-b = F.jac_T_dot(None, y)
 
-b_test = np.asarray((K.T * y).T)[0]
 
-x_f = cg(A, b_test)[0]
-
-mag_data_rec = MagData(a, np.zeros((3,)+dim))
-
-mag_data_rec.mag_vec = x_f
 
-mag_data_rec.quiver_plot3d()
+#K = np.asmatrix([F.jac_dot(None, np.eye(24)[:, i]) for i in range(24)]).T
+#lam = 10. ** -10
+#KTK = K.T * K + lam * np.asmatrix(np.eye(24))
+#
+#A = KTK#np.array([F(np.eye(81)[:, i]) for i in range(81)])
+#
+#b = F.jac_T_dot(None, y)
+#
+#b_test = np.asarray((K.T.dot(y)).T)[0]
+#
+#x_f = cg(A, b_test)[0]
+#
+#mag_data_rec = MagData(a, np.zeros((3,)+dim))
+#
+#mag_data_rec.mag_vec = x_f
+#
+#mag_data_rec.quiver_plot3d()