Allow Inhomogeneous Measurement Patterns (Accelerate) (#88)

* Inhomogeneous measurement pattern
* Optimize the calculation of smearing matrix B and the Jacobian

pyeit/eit/fem.py

@@ -180,12 +180,6 @@ def _check_mesh_protocol_compatibility(
                 stacklevel=2,
             )
 
-        if m_n_el < p_n_el:
-            raise ValueError(
-                f"Protocol is not compatible with mesh :\
-The mesh use {m_n_el} electrodes, and the protocol use only {p_n_el} electrodes "
-            )
-
     def solve_eit(
         self,
         perm: Optional[Union[int, float, complex, np.ndarray]] = None,
@@ -239,37 +233,35 @@ def compute_jac(
             Jacobian matrix, initial boundary voltage meas. extimation v0
 
         """
-        # update k if necessary and calculate r=inv(k)
+        # update k if necessary and calculate r=inv(k), dense matrix, slow
         self.assemble_pde(perm)
-        r_el = la.inv(self.kg.toarray())[self.mesh.el_pos]
-
-        # calculate v, jac per excitation pattern (ex_line)
-        _jac = np.zeros(
-            (self.protocol.n_exc, self.protocol.n_meas, self.mesh.n_elems),
-            dtype=self.mesh.dtype,
-        )
-        # v = np.zeros((self.protocol.n_exc, self.protocol.n_meas), dtype=self.mesh.dtype)
-
+        r_mat = la.inv(self.kg.toarray())[self.mesh.el_pos]
+        r_el = np.full((self.protocol.ex_mat.shape[0],) + r_mat.shape, r_mat)
+        # nodes potential
         f = self.solve_vectorized(self.protocol.ex_mat)
-        r_el = np.full((self.protocol.ex_mat.shape[0],) + r_el.shape, r_el)
-        v = subtract_row_vectorized(f[:, self.mesh.el_pos], self.protocol.meas_mat)
+        f_el = f[:, self.mesh.el_pos]
+        # build measurements and node resistance
+        v = subtract_row_vectorized(f_el, self.protocol.meas_mat)
         ri = subtract_row_vectorized(r_el, self.protocol.meas_mat)
-        # Build Jacobian matrix column wise (element wise)
-        #    Je = Re*Ke*Ve = (nex3) * (3x3) * (3x1)
-        for i in range(self.protocol.ex_mat.shape[0]):
-            for e, ijk in enumerate(self.mesh.element):
-                _jac[i, :, e] = np.dot(np.dot(ri[i][:, ijk], self.se[e]), f[i][ijk])
+        v0 = v.reshape(-1)
 
+        ## calculate v, jac per excitation (ex_line)
+        # _jac = np.zeros((self.protocol.n_meas, self.mesh.n_elems), dtype=self.mesh.dtype)
         # for i, ex_line in enumerate(self.protocol.ex_mat):
         #     f = self.solve(ex_line)
         #     v[i] = subtract_row(f[self.mesh.el_pos], self.protocol.meas_mat[i])
         #     ri = subtract_row(r_el, self.protocol.meas_mat[i])
         #     for (e, ijk) in enumerate(self.mesh.element):
         #         _jac[i, :, e] = np.dot(np.dot(ri[:, ijk], self.se[e]), f[ijk])
+        # jac = np.concatenate(_jac)
 
-        # measurement protocol
-        jac = np.concatenate(_jac)
-        v0 = v.reshape(-1)
+        # Build Jacobian matrix element wise (column wise)
+        # Je = Re*Ke*Ve = (n_measx3) * (3x3) * (3x1)
+        jac = np.zeros((self.protocol.n_meas, self.mesh.n_elems), dtype=self.mesh.dtype)
+        indices = self.protocol.meas_mat[:, 2]
+        f_n = f[indices]  # replica of potential on nodes of difference excitations
+        for e, ijk in enumerate(self.mesh.element):
+            jac[:, e] = np.sum(np.dot(ri[:, ijk], self.se[e]) * f_n[:, ijk], axis=1)
 
         # Jacobian normalization: divide each row of J (J[i]) by abs(v0[i])
         if normalize:
@@ -296,6 +288,16 @@ def compute_b_matrix(
             back-projection mappings (smear matrix); shape(n_exc, n_pts, 1), dtype= bool
         """
         self.assemble_pde(perm)
+        f = self.solve_vectorized(self.protocol.ex_mat)
+        f_el = f[:, self.mesh.el_pos]
+        return _smear_nd(f, f_el, self.protocol.meas_mat)
+
+    def compute_b_matrix_iter(
+        self,
+        perm: Optional[Union[int, float, complex, np.ndarray]] = None,
+    ):
+        """Compute back-projection mappings (smear matrix) [obsolete]"""
+        self.assemble_pde(perm)
         b_mat = np.zeros((self.protocol.n_exc, self.protocol.n_meas, self.mesh.n_nodes))
 
         for i in range(self.protocol.n_exc):
@@ -313,7 +315,7 @@ def compute_b_matrix(
 
 def _smear(f: np.ndarray, fb: np.ndarray, pairs: np.ndarray):
     """
-    Build smear matrix B for bp for one exitation
+    Build smear matrix B for bp for one exitation [obsolete]
 
     used for the smear matrix computation by @ChabaneAmaury
 
@@ -337,11 +339,9 @@ def _smear(f: np.ndarray, fb: np.ndarray, pairs: np.ndarray):
     return (f_min < f) & (f <= f_max)
 
 
-def smear_nd(
-    f: np.ndarray, fb: np.ndarray, meas_pattern: np.ndarray, new: bool = False
-) -> np.ndarray:
+def _smear_nd(f: np.ndarray, fb: np.ndarray, meas_pattern: np.ndarray) -> np.ndarray:
     """
-    Build smear matrix B for bp
+    Build smear matrix B for bp (vectorized version using exc_idx from meas_pattern)
 
     Parameters
     ----------
@@ -350,40 +350,26 @@ def smear_nd(
     fb: np.ndarray
         potential on adjacent electrodes; shape (n_exc, n_el)
     meas_pattern: np.ndarray
-        electrodes numbering pairs; shape (n_exc, n_meas, 2)
-    new : bool, optional
-        flag to use new matrices based computation, by default False.
-        If `False` to smear-computation from ChabaneAmaury is used
+        electrodes numbering pairs; shape (n_meas_tot, 3)
 
     Returns
     -------
     np.ndarray
-        back-projection (smear) matrix; shape (n_exc, n_meas, n_pts), dtype= bool
+        back-projection (smear) matrix; shape (n_meas_tot, n_pts), dtype= bool
     """
-    if new:
-        # new implementation not much faster! :(
-        idx_meas_0 = meas_pattern[:, :, 0]
-        idx_meas_1 = meas_pattern[:, :, 1]
-        n_exc = meas_pattern.shape[0]  # number of excitations
-        n_meas = meas_pattern.shape[1]  # number of measurements per excitations
-        n_pts = f.shape[1]  # number of nodes
-        idx_exc = np.ones_like(idx_meas_0, dtype=int) * np.arange(n_exc).reshape(
-            n_exc, 1
-        )
-        f_min = np.minimum(fb[idx_exc, idx_meas_0], fb[idx_exc, idx_meas_1])
-        f_max = np.maximum(fb[idx_exc, idx_meas_0], fb[idx_exc, idx_meas_1])
-        # contruct matrices of shapes (n_exc, n_meas, n_pts) for comparison
-        f_min = np.repeat(f_min[:, :, np.newaxis], n_pts, axis=2)
-        f_max = np.repeat(f_max[:, :, np.newaxis], n_pts, axis=2)
-        f0 = np.repeat(f[:, :, np.newaxis], n_meas, axis=2)
-        f0 = f0.swapaxes(1, 2)
-        return np.array((f_min < f0) & (f0 <= f_max))
-    else:
-        # Replacing the below code by a faster implementation in Numpy
-        def b_matrix_init(k):
-            return _smear(f[k], fb[k], meas_pattern[k])
+    n = meas_pattern[:, 0]
+    m = meas_pattern[:, 1]
+    exc_id = meas_pattern[:, 2]
+    # (n_meas_tot,) voltages on electrodes
+    f_min = np.minimum(fb[exc_id, n], fb[exc_id, m])
+    f_max = np.maximum(fb[exc_id, n], fb[exc_id, m])
+    # contruct matrix of shapes (n_meas_tot, n_pts) for comparison
+    n_pts = f.shape[1]
+    f_min = np.repeat(f_min[:, np.newaxis], n_pts, axis=1)
+    f_max = np.repeat(f_max[:, np.newaxis], n_pts, axis=1)
+    f_pts = f[exc_id]  # voltages on nodes of all excitations
 
-        return np.array(list(map(b_matrix_init, np.arange(f.shape[0]))))
+    return np.array((f_min < f_pts) & (f_pts <= f_max))
 
 
 def subtract_row(v: np.ndarray, meas_pattern: np.ndarray):
@@ -411,24 +397,24 @@ def subtract_row(v: np.ndarray, meas_pattern: np.ndarray):
 def subtract_row_vectorized(v: np.ndarray, meas_pattern: np.ndarray):
     """
     Build the voltage differences on axis=1 using the meas_pattern.
-    v_diff[k] = v[i, :] - v[j, :]
+    v_diff[k] = v[exc_id, i] - v[exc_id, j]
 
     New implementation 33% less computation time
 
     Parameters
     ----------
     v: np.ndarray
-        Nx1 boundary measurements vector or NxM matrix; shape (n_exc,n_el,1)
+        (n_exc, n_el) boundary measurements or (n_exc, (n_el, n_element)) nodes resistance
     meas_pattern: np.ndarray
-        Nx2 subtract_row pairs; shape (n_exc, n_meas, 2)
+        Nx2 subtract_row pairs; shape (n_meas_tot, 3)
 
     Returns
     -------
     np.ndarray
         difference measurements v_diff
     """
-    idx = np.indices(meas_pattern.shape[:-1])[0]
-    return v[idx, meas_pattern[:, :, 0]] - v[idx, meas_pattern[:, :, 1]]
+    idx = meas_pattern[:, 2]
+    return v[idx, meas_pattern[:, 0]] - v[idx, meas_pattern[:, 1]]
 
 
 def assemble(

pyeit/eit/protocol.py

@@ -21,7 +21,7 @@ class PyEITProtocol:
     ex_mat: np.ndarray
         excitation matrix (pairwise)
     meas_mat: np.ndarray
-        measurement matrix (differential pairs)
+        measurement matrix (differential pairs), support inhomogeneous number of measurements per excitation pair.
     keep_ba: np.ndarray
         boolean array index for keeping measurements
     """
@@ -81,25 +81,22 @@ def _check_meas_mat(self, meas_mat: np.ndarray) -> np.ndarray:
         n_exc : int
             number of excitations/stimulations
         meas_pattern : np.ndarray, optional
-           measurements pattern / subtract_row pairs [N, M] to check; shape (n_exc, n_meas_per_exc, 2)
+           measurements pattern / subtract_row pairs [N, M]; shape (n_meas_tot, 3)
 
         Returns
         -------
         np.ndarray
-            measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2)
+            measurements pattern / subtract_row pairs [N, M]; shape (n_meas_tot, 3)
 
         Raises
         ------
         TypeError
-            raised if meas_pattern is not a np.ndarray of shape (n_exc, : , 2)
+            raised if meas_pattern is not a np.ndarray of shape (n_meas_tot, 3)
         """
         if not isinstance(meas_mat, np.ndarray):
             raise TypeError(f"Wrong type of {type(meas_mat)=}, expected an ndarray;")
-        # test shape is something like (n_exc, :, 2)
-        if meas_mat.ndim != 3 or meas_mat.shape[::2] != (self.n_exc, 2):
-            raise TypeError(
-                f"Wrong shape of {meas_mat.shape=}, should be ({self.n_exc}, n_meas_per_exc, 2);"
-            )
+        if meas_mat.ndim != 2 or meas_mat.shape[-1] != 3:
+            raise TypeError(f"{meas_mat.shape=} must be (n_meas_tot, 3);")
 
         return meas_mat
 
@@ -126,29 +123,23 @@ def n_meas(self) -> int:
         Returns
         -------
         int
-            number of measurements per excitations
+            total amount of measurements (n_meas_tot)
         """
-        return self.meas_mat.shape[1]
-
-    @property
-    def n_meas_tot(self) -> int:
-        """
-        Returns
-        -------
-        int
-            total amount of measurements
-        """
-        return self.n_meas * self.n_exc
+        return self.meas_mat.shape[0]
 
     @property
     def n_el(self) -> int:
         """
         Returns
         -------
         int
-            number of electrodes used in the excitation and
+            infer the number of electrodes used in the excitation and measurements patterns,
+            where the electrodes are numbered [0, n_el-1].
         """
-        return int(max(max(self.ex_mat.flatten()), max(self.meas_mat.flatten()))) + 1
+        return (
+            int(max(max(self.ex_mat.flatten()), max(self.meas_mat[:, :-1].flatten())))
+            + 1
+        )
 
 
 def create(
@@ -206,8 +197,8 @@ def build_meas_pattern_std(
     parser: Union[str, List[str]] = "std",
 ) -> Tuple[np.ndarray, np.ndarray]:
     """
-    Build the measurement pattern (subtract_row-voltage pairs [N, M])
-    for all excitations on boundary electrodes.
+    Build the measurement pattern (subtract_row-voltage pairs [N, M]) for all excitations on boundary electrodes.
+    The excitation index (exc_id) are also recorded for computing subtract_row_vectorized and smear_nd.
 
     we direct operate on measurements or Jacobian on electrodes,
     so, we can use LOCAL index in this module, do not require el_pos.
@@ -241,30 +232,33 @@ def build_meas_pattern_std(
     Returns
     -------
     diff_op: np.ndarray
-        measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2)
+        measurements pattern / subtract_row pairs, and the excitation indice;
+        shape (n_meas_tot, 3), for each row, it represents [Ni, Mi, exc_id]
     keep_ba: np.ndarray
-        (n_exc*n_meas_per_exc,) boolean array
+        (n_meas_tot,) boolean array
     """
     if not isinstance(parser, list):  # transform parser into list
         parser = [parser]
     meas_current = "meas_current" in parser
     fmmu_rotate = any(p in ("fmmu", "rotate_meas") for p in parser)
 
     diff_op, keep_ba = [], []
-    for ex_line in ex_mat:
-        a, b = ex_line[0], ex_line[1]
+    for exc_id, exc_line in enumerate(ex_mat):
+        a, b = exc_line[0], exc_line[1]
         i0 = a if fmmu_rotate else 0
+        # build [[m, n, idx]_i] array
         m = (i0 + np.arange(n_el)) % n_el
         n = (m + step) % n_el
-        meas_pattern = np.vstack([n, m]).T
+        idx = exc_id * np.ones(n_el)
+        meas_pattern = np.vstack([n, m, idx]).T
 
         diff_keep = np.logical_and.reduce((m != a, m != b, n != a, n != b))
         keep_ba.append(diff_keep)
         if not meas_current:
             meas_pattern = meas_pattern[diff_keep]
-        diff_op.append(meas_pattern)
+        diff_op.append(meas_pattern.astype(int))
 
-    return np.array(diff_op), np.array(keep_ba).ravel()
+    return np.vstack(diff_op), np.array(keep_ba).ravel()
 
 
 def build_exc_pattern_std(n_el: int = 16, dist: int = 1) -> np.ndarray:

tests/test_fem.py

@@ -115,7 +115,7 @@ def test_meas_pattern(self):
                 protocol = _protocol_obj(ex_mat, n_el, 1, parser)
                 fwd = pyeit.eit.fem.EITForward(mesh, protocol)
                 diff_truth = _meas_pattern(ex_line, n_el, 1, parser)
-                diff = fwd.protocol.meas_mat[0]
+                diff = fwd.protocol.meas_mat[:, :-1]
 
                 assert np.allclose(diff, diff_truth)
 
@@ -132,18 +132,25 @@ def test_subtract_row(self):
         self.assertTrue(np.allclose(vd.ravel(), vd_truth.ravel()))
 
     def test_subtract_row_vectorized(self):
-        """calculate f[diff_op[0]] - f[diff_op[1]]"""
-        n_exe = 10
+        """calculate f[exc_id, diff_op[0]] - f[exc_id, diff_op[1]]"""
         n_el = 16
-        n_rows = 3
-        v = np.full((n_rows, n_el), np.random.randn(n_el))
-        diff_pairs = np.array(
-            [[np.random.permutation(n_el)[:2] for _ in range(n_exe)]] * n_rows
-        )
-        vd_truth = np.zeros((n_rows, 10))
-        for i in range(vd_truth.shape[0]):
-            vd_truth[i] = np.array([v[i, d[0]] - v[i, d[1]] for d in diff_pairs[i]])
-        vd = pyeit.eit.fem.subtract_row_vectorized(v, diff_pairs)
+        n_meas = 16
+        n_exc = 3
+        n_meas_tot = n_exc * n_meas
+        v = np.full((n_exc, n_el), np.random.randn(n_el))
+        # build measurement pattern, [m, n, exc_id] per row
+        diff_pairs = []
+        for i in range(n_exc):
+            for _ in range(n_meas):
+                diff_pairs.append(np.hstack([np.random.permutation(n_el)[:2], i]))
+        meas_pattern = np.vstack(diff_pairs)
+        print(meas_pattern)
+        # calculate ground truth
+        vd_truth = np.zeros((n_meas_tot,))
+        for i in range(n_meas_tot):
+            v_exc = v[meas_pattern[i, 2]]
+            vd_truth[i] = v_exc[meas_pattern[i, 0]] - v_exc[meas_pattern[i, 1]]
+        vd = pyeit.eit.fem.subtract_row_vectorized(v, meas_pattern)
 
         self.assertTrue(vd_truth.size == vd.size)
         self.assertTrue(np.allclose(vd.ravel(), vd_truth.ravel()))