Added numerical core routines for Jacobian multiply.

pyramid/numcore/phase_mag_real.pyx

@@ -92,8 +92,8 @@ def get_jacobi_core(
 @cython.wraparound(False)
 def get_jacobi_core(
     unsigned int v_dim, unsigned int u_dim,
-    np.ndarray v_phi, np.ndarray u_phi,
-    np.ndarray jacobi):
+    double[:, :] v_phi, double[:, :] u_phi,
+    double[:, :] jacobi):
     '''DOCSTRING!'''
     # TODO: Docstring!!!
     cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
@@ -105,3 +105,62 @@ def get_jacobi_core(
             # v_dim*u_dim columns for the v-component (note the minus!):
             jacobi[:, u_dim*v_dim+i+u_dim*j] = \
                 -np.reshape(v_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def multiply_jacobi_core(
+    unsigned int v_dim, unsigned int u_dim,
+    double[:, :] v_phi, double[:, :] u_phi,
+    double[:] vector,
+    double[:] result):
+
+    cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, ri, u, v, siz
+    cdef double v0, v1
+    siz = u_dim * v_dim
+
+    s = 0
+    for i in range(u_dim):
+        for j in range(v_dim):
+            u_min = (u_dim - 1) - i
+            v_min = (v_dim - 1) - j
+
+            ri = 0
+            for v in range(v_min, v_min + v_dim):
+                for u in range(u_min, u_min + u_dim):
+                    result[ri] += vector[s] * u_phi[v, u]
+                    result[ri] -= vector[s + siz] * v_phi[v, u]
+                    ri += 1
+            s += 1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def multiply_jacobi_core2(
+    int v_dim, int u_dim,
+    double[:, :] v_phi, double[:, :] u_phi,
+    double[:] vector,
+    double[:] result):
+
+    cdef int s1, s2, i, j, u_min, u_max, v_min, v_max, ri, u, v, siz, j_min, j_max, i_min, i_max
+
+    siz = u_dim * v_dim
+
+    for v in range(2 * v_dim - 1):
+        for u in range(2 * u_dim - 1):
+            i_min = max(0, (u_dim - 1) - u)
+            i_max = min(u_dim, ((2 * u_dim) - 1) - u)
+            for i in range(i_min, i_max):
+                s1 = i * u_dim
+                s2 = s1 + siz
+
+                j_min = max(0, (v_dim - 1) - v)
+                j_max = min(v_dim, ((2 * v_dim) - 1) - v)
+
+                u_min = (u_dim - 1) - i
+                v_min = (v_dim - 1) - j_min
+                ri = (v - ((v_dim - 1) - j_min)) * u_dim + (u - u_min)
+                for j in range(j_min, j_max):
+                    result[ri] += vector[s1 + j] * u_phi[v, u]
+                    result[ri] -= vector[s2 + j] * v_phi[v, u]
+                    ri += u_dim

pyramid/phasemapper.py

@@ -99,6 +99,22 @@ class Kernel:
             result += vector[s+size]*-self.v[v_min:v_max, u_min:u_max].reshape(-1)  # v        
         return result
 
+    def multiply_jacobi_core(self, vector):
+        v_dim, u_dim = self.dim
+        size = v_dim * u_dim
+        result = np.zeros(size)
+        nc.multiply_jacobi_core(
+            v_dim, u_dim, self.v, self.u, vector, result)
+        return result
+
+    def multiply_jacobi_core2(self, vector):
+        v_dim, u_dim = self.dim
+        size = v_dim * u_dim
+        result = np.zeros(size)
+        nc.multiply_jacobi_core2(
+            v_dim, u_dim, self.v, self.u, vector, result)
+        return result
+
     def multiply_jacobi_T(self, vector):
         # TODO: Docstring!!!
         # vector: v_dim*u_dim elements for u_mag and v_dim*u_dim elements for v_mag