fixed some bugs, got reconstruction to work.

pyramid/costfunction.py

@@ -53,8 +53,9 @@ class Costfunction(object):
         # Extract important information:
         self.y = data_set.phase_vec
         self.Se_inv = data_set.Se_inv
-        self.n = data_set.n * 3
+        self.n = data_set.m
         self.m = len(self.y)
+        self.chisq, self.chisq_a, self.chisq_m = None, None, None
         self.LOG.debug('Created '+str(self))
 
     def __repr__(self):
@@ -73,7 +74,9 @@ class Costfunction(object):
     def __call__(self, x):
         self.LOG.debug('Calling __call__')
         delta_y = self.fwd_model(x) - self.y
-        self.chisq = delta_y.dot(self.Se_inv.dot(delta_y)) + self.regularisator(x)
+        self.chisq_m = delta_y.dot(self.Se_inv.dot(delta_y))
+        self.chisq_a = self.regularisator(x)
+        self.chisq = self.chisq_m + self.chisq_a
         return self.chisq
 
     def jac(self, x):

pyramid/dataset.py

@@ -88,7 +88,7 @@ class DataSet(object):
             'Dimension has to be a tuple of length 3!'
         if mask is not None:
             assert mask.shape == dim, 'Mask dimensions must match!'
-            self.m = 3 * np.sum(self.mask)
+            self.m = 3 * np.sum(mask)
         else:
             self.m = 3 * np.prod(dim)
         self.a = a

pyramid/forwardmodel.py

@@ -59,6 +59,7 @@ class ForwardModel(object):
 
     def __call__(self, x):
         self.LOG.debug('Calling __call__')
+        self.mag_data.magnitude[:] = 0
         self.mag_data.set_vector(x, self.data_set.mask)
         # TODO: Multiprocessing
         result = np.zeros(self.n)
@@ -89,6 +90,7 @@ class ForwardModel(object):
 
         '''
         self.LOG.debug('Calling jac_dot')
+        self.mag_data.magnitude[:] = 0
         self.mag_data.set_vector(vector, self.data_set.mask)
         result = np.zeros(self.n)
         hp = self.hook_points
@@ -119,7 +121,7 @@ class ForwardModel(object):
         '''
         self.LOG.debug('Calling jac_T_dot')
 
-        result = np.zeros(self.m)
+        result = np.zeros(3*np.prod(self.data_set.dim))
         hp = self.hook_points
         for i, projector in enumerate(self.data_set.projectors):
             vec = vector[hp[i]:hp[i+1]]

pyramid/kernel.py

@@ -105,3 +105,38 @@ class Kernel(object):
         self.LOG.debug('Calling __str__')
         return 'Kernel(a=%s, dim_uv=%s, numcore=%s, geometry=%s)' % \
             (self.a, self.dim_uv, self.numcore, self.geometry)
+
+    def _multiply_jacobi(self, vector):
+        self.LOG.debug('Calling _multiply_jacobi')
+        v_dim, u_dim = self.dim_uv
+        assert len(vector) == 2 * self.size, \
+            'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.size)
+        result = np.zeros(self.size)
+        # Iterate over all contributing pixels (numbered consecutively)
+        for s in range(self.size):  # column-wise (two columns at a time, u- and v-component)
+            i = s % u_dim  # u-coordinate of current contributing pixel
+            j = int(s/u_dim)  # v-coordinate of current ccontributing pixel
+            u_min = (u_dim-1) - i  # u_dim-1: center of the kernel
+            u_max = (2*u_dim-1) - i  # = u_min + u_dim
+            v_min = (v_dim-1) - j  # v_dim-1: center of the kernel
+            v_max = (2*v_dim-1) - j  # = v_min + v_dim
+            result += vector[s] * self.u[v_min:v_max, u_min:u_max].reshape(-1)  # u
+            result -= vector[s+self.size] * self.v[v_min:v_max, u_min:u_max].reshape(-1)  # v
+        return result
+
+    def _multiply_jacobi_T(self, vector):
+        self.LOG.debug('Calling _multiply_jacobi_T')
+        v_dim, u_dim = self.dim_uv
+        result = np.zeros(2*self.size)
+        # Iterate over all contributing pixels (numbered consecutively):
+        for s in range(self.size):  # row-wise (two rows at a time, u- and v-component)
+            i = s % u_dim  # u-coordinate of current contributing pixel
+            j = int(s/u_dim)  # v-coordinate of current contributing pixel
+            u_min = (u_dim-1) - i  # u_dim-1: center of the kernel
+            u_max = (2*u_dim-1) - i  # = u_min + u_dim
+            v_min = (v_dim-1) - j  # v_dim-1: center of the kernel
+            v_max = (2*v_dim-1) - j  # = v_min + v_dim
+            result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1))  # u
+            result[s+self.size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))  # v
+        return result
+#

pyramid/magdata.py

@@ -83,7 +83,8 @@ class MagData(object):
     @mag_vec.setter
     def mag_vec(self, mag_vec):
         assert isinstance(mag_vec, np.ndarray), 'Vector has to be a numpy array!'
-        assert np.size(mag_vec) == 3*np.prod(self.dim), 'Vector has to match magnitude dimensions!'
+        assert np.size(mag_vec) == 3*np.prod(self.dim), \
+            'Vector has to match magnitude dimensions! {} {}'.format(mag_vec.shape, 3*np.prod(self.dim))
         self.magnitude = mag_vec.reshape((3,)+self.dim)
 
     def __init__(self, a, magnitude):
@@ -281,9 +282,13 @@ class MagData(object):
             Order is: first all `x`-, then all `y`-, then all `z`-components.
 
         '''
-        return np.reshape([self.magnitude[2][mask],
+        if mask is not None:
+            return np.reshape([self.magnitude[0][mask],
                            self.magnitude[1][mask],
-                           self.magnitude[0][mask]], -1)
+                           self.magnitude[2][mask]], -1)
+        else:
+            return self.mag_vec
+
 
     def set_vector(self, vector, mask=None):
         '''Set the magnetic components of the masked pixels to the values specified by `vector`.
@@ -304,9 +309,9 @@ class MagData(object):
         assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
         count = np.size(vector)/3
         if mask is not None:
-            self.magnitude[2][mask] = vector[:count]  # x-component
+            self.magnitude[0][mask] = vector[:count]  # x-component
             self.magnitude[1][mask] = vector[count:2*count]  # y-component
-            self.magnitude[0][mask] = vector[2*count:]  # z-component
+            self.magnitude[2][mask] = vector[2*count:]  # z-component
         else:
             self.mag_vec = vector
 

pyramid/phasemap.py

@@ -293,6 +293,7 @@ class PhaseMap(object):
         phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
         a = phase_file.a
         phase = phase_file.variables['phase'][:]
+        #phase = phase[28:36, 28:36]
         phase_file.close()
         return PhaseMap(a, phase)
 

pyramid/phasemapper.py

@@ -138,7 +138,9 @@ class PhaseMapperRDFC(PhaseMapper):
         assert len(vector) == self.m, \
             'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
         u_mag, v_mag = np.reshape(vector, (2,)+self.kernel.dim_uv)
-        return self._convolve(u_mag, v_mag)
+        result = self._convolve(u_mag, v_mag).reshape(-1)
+        return result
+        #return self.kernel._multiply_jacobi(vector)
 
     def jac_T_dot(self, vector):
         '''Calculate the product of the transposed Jacobi matrix with a given `vector`.
@@ -161,10 +163,10 @@ class PhaseMapperRDFC(PhaseMapper):
             'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
         # TODO: directly? ask Jörn again!
         result = np.zeros(self.m)
+        v_dim, u_dim = self.kernel.dim_uv
         nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1],
                                    self.kernel.u, self.kernel.v, vector, result)
-        return result
-
+        return self.kernel._multiply_jacobi_T(vector)
 
 class PhaseMapperRDRC(PhaseMapper):
 
@@ -317,8 +319,8 @@ class PhaseMapperRDRC(PhaseMapper):
                 u_max = (2*u_dim-1) - i  # = u_min + u_dim
                 v_min = (v_dim-1) - j  # v_dim-1: center of the kernel
                 v_max = (2*v_dim-1) - j  # = v_min + v_dim
-                result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1))
-                result[s+self.n] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))
+                result[s] = np.sum(vector*self.kernel.u[v_min:v_max, u_min:u_max].reshape(-1))
+                result[s+self.n] = np.sum(vector*-self.kernel.v[v_min:v_max, u_min:u_max].reshape(-1))
         return result
 
 

pyramid/reconstruction.py

@@ -119,11 +119,11 @@ def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
     fwd_model = ForwardModel(data)
     cost = Costfunction(data, regularisator)
     print cost(np.zeros(cost.n))
-    x_opt = jutil.cg.conj_grad_minimize(cost)
+    x_opt = jutil.cg.conj_grad_minimize(cost, max_iter=20)
     print cost(x_opt)
     # Create and return fitting MagData object:
     mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
-    mag_opt.mag_vec = x_opt
+    mag_opt.set_vector(x_opt, data.mask)
     return mag_opt
 
 
@@ -152,15 +152,38 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
     LOG.debug('Calling optimize_cg')
     if first_guess is None:
         first_guess = MagData(data.a, np.zeros((3,)+data.dim))
-    x_0 = first_guess.mag_vec
+    x_0 = first_guess.get_vector(data.mask)
     fwd_model = ForwardModel(data)
     cost = Costfunction(data, regularisator)
+
+#    proj = fwd_model.data_set.projectors[0]
+#    proj_jac1 = np.array([proj.jac_dot(np.eye(proj.m, 1, -i).squeeze()) for i in range(proj.m)])
+#    proj_jac2 = np.array([proj.jac_T_dot(np.eye(proj.n, 1, -i).squeeze()) for i in range(proj.n)])
+#    print jac1, jac2.T, abs(jac1-jac2.T).sum()
+#    print jac1.shape, jac2.shape
+
+#    pm = fwd_model.phase_mappers[proj.dim_uv]
+#    pm_jac1 = np.array([pm.jac_dot(np.eye(pm.m)[:, i]) for i in range(pm.m)])
+#    pm_jac2 = np.array([pm.jac_T_dot(np.eye(pm.n)[:, i]) for i in range(pm.n)])
+#    print jac1, jac2.T, abs(jac1-jac2.T).sum()
+#    print jac1.shape, jac2.shape
+
+
+ #   jac1 = np.array([fwd_model.jac_dot(x_0, np.eye(fwd_model.m)[:, i]) for i in range(fwd_model.m)])
+ #   jac2 = np.array([fwd_model.jac_T_dot(x_0, np.eye(fwd_model.n)[:, i]) for i in range(fwd_model.n)])
+ #   print proj_jac1.dot(pm_jac1)
+ #   print (pm_jac2.dot(proj_jac2)).T
+ #   print jac1
+#    print jac2.T
+#    print abs(jac1-jac2.T).sum()
+#    print jac1.shape, jac2.shape
+
     assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
-    result = jutil.minimizer.minimize(cost, x_0, options={"conv_rel":1e-20}, tol={"max_iteration":1})
+    result = jutil.minimizer.minimize(cost, x_0, options={"conv_rel":1e-2}, tol={"max_iteration":4})
     x_opt = result.x
     print cost(x_opt)
     mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
-    mag_opt.mag_vec = x_opt
+    mag_opt.set_vector(x_opt, data.mask)
     return mag_opt
 
 

pyramid/regularisator.py

@@ -110,28 +110,14 @@ class ZeroOrderRegularisator(Regularisator):
         self.LOG.debug('Created '+str(self))
 
 
-#class FirstOrderRegularisator(Regularisator):
-#    # TODO: Docstring!
-#
-#    def __init__(self, fwd_model, lam, x_a=None):
-#        size_3d = fwd_model.size_3d
-#        dim = fwd_model.dim
-#        converter = IndexConverter(dim)
-#        row = []
-#        col = []
-#        data = []
-#
-#        for i in range(size_3d):
-#            neighbours = converter.find_neighbour_ind(i)
-#
-#            Sa_inv = csr_matrix(coo_matrix(data, (rows, columns)), shape=(3*size_3d, 3*size_3d))
-#
-#        term2 = []
-#        for i in range(3):
-#            component = mag_data[i, ...]
-#            for j in range(3):
-#                if component.shape[j] > 1:
-#                    term2.append(np.diff(component, axis=j).reshape(-1))
-#
-#        super(FirstOrderRegularisator, self).__init__(Sa_inv, x_a)
-#        self.LOG.debug('Created '+str(self))
+class FirstOrderRegularisator(Regularisator):
+    # TODO: Docstring!
+
+    def __init__(self, mask, lam, x_a=None):
+        import jutil
+        D0 = jutil.diff.get_diff_operator(mask, 0, 3)
+        D1 = jutil.diff.get_diff_operator(mask, 1, 3)
+        D = jutil.operator.VStack([D0, D1])
+        norm = jutil.norms.WeightedL2Square(D)
+        super(FirstOrderRegularisator, self).__init__(norm, lam)
+        self.LOG.debug('Created '+str(self))