add perm option of EIT solvers. update testsuite and default examples…

… for EIT solvers.

examples/eit_dynamic_bp.py

@@ -23,7 +23,6 @@
     mesh_obj = mesh.create(n_el, h0=0.1)
 
 """ 1. problem setup """
-
 anomaly = PyEITAnomaly_Circle(center=[0.5, 0.5], r=0.1, perm=10.0)
 mesh_new = mesh.set_perm(mesh_obj, anomaly=anomaly, background=1.0)
 

examples/eit_dynamic_greit.py

@@ -52,8 +52,8 @@
 
 """ 3. Construct using GREIT """
 eit = greit.GREIT(mesh_obj, protocol_obj)
-eit.setup(p=0.50, lamb=0.001)
-ds = eit.solve(v1, v0)
+eit.setup(p=0.50, lamb=0.01, perm=1, jac_normalized=True)
+ds = eit.solve(v1, v0, normalize=True)
 x, y, ds = eit.mask_value(ds, mask_value=np.NAN)
 
 # show alpha

examples/eit_dynamic_jac.py

@@ -49,7 +49,7 @@
 # (mostly) the shape and the electrode positions are not exactly the same
 # as in mesh generating the jac, then data must be normalized.
 eit = jac.JAC(mesh_obj, protocol_obj)
-eit.setup(p=0.5, lamb=0.01, method="kotre")
+eit.setup(p=0.5, lamb=0.01, method="kotre", perm=1, jac_normalized=True)
 ds = eit.solve(v1, v0, normalize=True)
 ds_n = sim2pts(pts, tri, np.real(ds))
 

examples/eit_dynamic_svd.py

@@ -50,7 +50,7 @@
 # (mostly) the shape and the electrode positions are not exactly the same
 # as in mesh generating the jac, then data must be normalized.
 eit = svd.SVD(mesh_obj, protocol_obj)
-eit.setup(n=35, method="svd")
+eit.setup(n=50, method="svd", perm=1, jac_normalized=True)
 ds = eit.solve(v1, v0, normalize=True)
 ds_n = sim2pts(pts, tri, np.real(ds))
 

pyeit/eit/__init__.py

@@ -13,10 +13,12 @@
 """
 from .bp import BP
 from .jac import JAC
+from .svd import SVD
 from .greit import GREIT
 
 __all__ = [
     "BP",
     "JAC",
+    "SVD",
     "GREIT",
 ]

pyeit/eit/bp.py

@@ -5,27 +5,32 @@
 # Distributed under the (new) BSD License. See LICENSE.txt for more info.
 from __future__ import division, absolute_import, print_function
 
+from typing import Union
 import numpy as np
 from .base import EitBase
 
 
 class BP(EitBase):
     """A naive inversion of (Euclidean) back projection."""
 
-    def setup(self, weight: str = "none") -> None:
+    def setup(
+        self, weight: str = "none", perm: Union[int, float, np.ndarray] = None
+    ) -> None:
         """
         Setup BP solver
 
         Parameters
         ----------
         weight : str, optional
             BP parameter, by default "none"
+        perm : Union[int, float, np.ndarray], optional
+            If perm is not None, a prior of perm distribution is used to build the smear matrix
         """
         self.params = {"weight": weight}
 
         # build the weighting matrix
         # BP: in node imaging, H is the smear matrix (transpose of B)
-        self.B = self.fwd.compute_b_matrix()
+        self.B = self.fwd.compute_b_matrix(perm=perm)
         self.H = self._compute_h(b_matrix=self.B)
         self.is_ready = True
 

pyeit/eit/fem.py

@@ -165,7 +165,9 @@ def solve_eit(
             simulated boundary voltage measurements; shape(n_exe*n_el,)
         """
         self.assemble_pde(perm)
-        v = np.zeros((self.protocol.n_exc, self.protocol.n_meas))
+        v = np.zeros(
+            (self.protocol.n_exc, self.protocol.n_meas), dtype=self.mesh.perm.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])
@@ -206,7 +208,9 @@ def compute_jac(
             (self.protocol.n_exc, self.protocol.n_meas, self.mesh.n_elems),
             dtype=self.mesh.perm.dtype,
         )
-        v = np.zeros((self.protocol.n_exc, self.protocol.n_meas))
+        v = np.zeros(
+            (self.protocol.n_exc, self.protocol.n_meas), dtype=self.mesh.perm.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])

pyeit/eit/greit.py

@@ -13,11 +13,9 @@
 # Distributed under the (new) BSD License. See LICENSE.txt for more info.
 from __future__ import absolute_import, division, print_function
 
-from typing import Tuple
-
+from typing import Tuple, Union
 import numpy as np
 import scipy.linalg as la
-
 from .base import EitBase
 from .interp2d import meshgrid, weight_sigmod
 
@@ -34,6 +32,7 @@ def setup(
         n: int = 32,
         s: float = 20.0,
         ratio: float = 0.1,
+        perm: Union[int, float, np.ndarray] = None,
         jac_normalized: bool = False,
     ) -> None:
         """
@@ -55,6 +54,8 @@ def setup(
             control the blur, by default 20.0
         ratio : float, optional
             desired ratio, by default 0.1
+        perm : Union[int, float, np.ndarray], optional
+            If perm is not None, a prior of perm distribution is used to build Jacobian
         jac_normalized : bool, optional
             normalize the jacobian using f0 computed from input perm, by
             default False
@@ -93,7 +94,7 @@ def setup(
         self.xg, self.yg, self.mask = meshgrid(self.mesh.node, n=n)
 
         w_mat = self._compute_grid_weights(self.xg, self.yg)
-        self.J, self.v0 = self.fwd.compute_jac(normalize=jac_normalized)
+        self.J, self.v0 = self.fwd.compute_jac(perm=perm, normalize=jac_normalized)
         self.H = self._compute_h(jac=self.J, w_mat=w_mat)
         self.is_ready = True
 

pyeit/eit/jac.py

@@ -8,10 +8,8 @@
 from __future__ import division, absolute_import, print_function
 
 from typing import Union
-
 import numpy as np
 import scipy.linalg as la
-
 from .base import EitBase
 
 
@@ -23,6 +21,7 @@ def setup(
         p: float = 0.20,
         lamb: float = 0.001,
         method: str = "kotre",
+        perm: Union[int, float, np.ndarray] = None,
         jac_normalized: bool = False,
     ) -> None:
         """
@@ -38,6 +37,8 @@ def setup(
             JAC parameters, by default 0.001
         method : str, optional
             regularization methods ("kotre", "lm", "dgn" ), by default "kotre"
+        perm : Union[int, float, np.ndarray], optional
+            If perm is not None, a prior of perm distribution is used to build jac
         jac_normalized : bool, optional
             normalize the jacobian using f0 computed from input perm, by
             default False
@@ -51,7 +52,7 @@ def setup(
         }
         # pre-compute H0 for dynamical imaging
         # H = (J.T*J + R)^(-1) * J.T
-        self.J, self.v0 = self.fwd.compute_jac(normalize=jac_normalized)
+        self.J, self.v0 = self.fwd.compute_jac(perm=perm, normalize=jac_normalized)
         self.H = self._compute_h(self.J, p, lamb, method)
         self.is_ready = True
 

pyeit/eit/svd.py

@@ -5,15 +5,22 @@
 # Distributed under the (new) BSD License. See LICENSE.txt for more info.
 from __future__ import division, absolute_import, print_function
 
+from typing import Union
 import numpy as np
-
 from .jac import JAC
 
 
 class SVD(JAC):
     """implementing a sensitivity-based EIT imaging class"""
 
-    def setup(self, n: int = 25, rcond: float = 1e-2, method: str = "svd") -> None:
+    def setup(
+        self,
+        n: int = 25,
+        rcond: float = 1e-2,
+        method: str = "svd",
+        perm: Union[int, float, np.ndarray] = None,
+        jac_normalized: bool = False,
+    ) -> None:
         """
         Setup of SVD solver, singular value decomposition based reconstruction.
 
@@ -27,9 +34,14 @@ def setup(self, n: int = 25, rcond: float = 1e-2, method: str = "svd") -> None:
             reconstruction method, by default "svd"
             'svd': SVD truncation,
             'pinv': pseudo inverse
+        perm : Union[int, float, np.ndarray], optional
+            If perm is not None, a prior of perm distribution is used to build jac
+        jac_normalized : bool, optional
+            normalize the jacobian using f0 computed from input perm, by
+            default False
         """
         # correct n_ord
-        self.J, self.v0 = self.fwd.compute_jac()
+        self.J, self.v0 = self.fwd.compute_jac(perm=perm, normalize=jac_normalized)
         nm, ne = self.J.shape
         n_ord = np.min([nm, ne, n])
 

tests/test_eit.py

@@ -1,36 +1,117 @@
 # test for eit
 import unittest
 import numpy as np
-import pyeit.eit
-from pyeit.eit.protocol import PyEITProtocol, build_meas_pattern_std
-from pyeit.mesh import PyEITMesh
 
+from pyeit.eit.fem import EITForward
+import pyeit.mesh as mesh
+from pyeit.mesh.wrapper import PyEITAnomaly_Circle
+import pyeit.eit.protocol as protocol
+import pyeit.eit
 
-def _mesh_obj():
-    """build a simple mesh, which is used in FMMU.CEM"""
-    node = np.array([[0, 0], [1, 1], [1, 2], [0, 1]])
-    element = np.array([[0, 1, 3], [1, 2, 3]])
-    perm = np.array([3.0, 1.0])  # assemble should not use perm.dtype
-    el_pos = np.array([1, 2])
 
-    return PyEITMesh(node=node, element=element, perm=perm, el_pos=el_pos, ref_node=3)
+def eit_loc_eval(ds: np.ndarray, mesh_obj: mesh.PyEITMesh, mode: str = "element"):
+    """
+    mode: string, default "element"
+        if mode = "element", ds are values on element
+        if mode = "node", ds are values on node (BP)
+    """
+    loc = np.argmax(np.abs(ds))
+    if mode == "node":
+        loc_xyz = mesh_obj.node[loc]
+    else:
+        loc_xyz = mesh_obj.elem_centers[loc]
 
+    ds_max = ds[loc]
+    ds_sign = np.sign(ds_max)
 
-def _protocol_obj(ex_mat, n_el, step_meas, parser_meas):
-    meas_mat = build_meas_pattern_std(ex_mat, n_el, step_meas, parser_meas)
-    return PyEITProtocol(ex_mat, meas_mat)
+    return loc_xyz, ds_max, ds_sign
 
 
 class TestFem(unittest.TestCase):
+    def setUp(self):
+        n_el = 16
+        self.n_el = n_el
+        self.mesh_obj = mesh.create(self.n_el, h0=0.1)
+
+        # set anomaly
+        self.anomaly = {"center": [0.5, 0.5], "r": 0.1, "perm": 10.0, "sign": True}
+        anomaly = PyEITAnomaly_Circle(
+            center=self.anomaly["center"],
+            r=self.anomaly["r"],
+            perm=self.anomaly["perm"],
+        )
+        self.mesh_new = mesh.set_perm(self.mesh_obj, anomaly=anomaly, background=1.0)
+        self.protocol_obj = protocol.create(
+            n_el, dist_exc=1, step_meas=1, parser_meas="std"
+        )
+
+        # calculate simulated data
+        self.fwd = EITForward(self.mesh_obj, self.protocol_obj)
+        self.v0 = self.fwd.solve_eit()
+        self.v1 = self.fwd.solve_eit(perm=self.mesh_new.perm)
+
     def test_bp(self):
         """test back projection"""
-        mesh = _mesh_obj()
-        ex_mat = np.array([[0, 1], [1, 0]])
-        protocol = _protocol_obj(ex_mat, mesh.n_el, 1, "meas_current")
-        solver = pyeit.eit.BP(mesh=mesh, protocol=protocol)
-        solver.setup()
+        eit = pyeit.eit.bp.BP(self.mesh_obj, self.protocol_obj)
+        # setup BP with no prior
+        eit.setup(weight="none", perm=1)
+        ds = 192.0 * eit.solve(self.v1, self.v0, normalize=False)
+
+        # evaluate
+        loc, ds_max, ds_sign = eit_loc_eval(ds, self.mesh_obj, mode="node")
+        # print(loc, ds_max, ds_sign)
+        loc = loc[: len(self.anomaly["center"])]
+        dist = np.linalg.norm(loc - self.anomaly["center"])
+
+        self.assertTrue(ds_sign == int(self.anomaly["sign"]))
+        self.assertTrue(dist < self.anomaly["r"])
+
+    def test_jac(self):
+        """test jac"""
+        eit = pyeit.eit.jac.JAC(self.mesh_obj, self.protocol_obj)
+        eit.setup(p=0.5, lamb=0.01, method="kotre", perm=1, jac_normalized=True)
+        ds = eit.solve(self.v1, self.v0, normalize=True)
+
+        # evaluate
+        loc, ds_max, ds_sign = eit_loc_eval(ds, self.mesh_obj, mode="element")
+        # print(loc, ds_max, ds_sign)
+        loc = loc[: len(self.anomaly["center"])]
+        dist = np.linalg.norm(loc - self.anomaly["center"])
+
+        self.assertTrue(ds_sign == int(self.anomaly["sign"]))
+        self.assertTrue(dist < self.anomaly["r"])
+
+    def test_svd(self):
+        """test SVD"""
+        eit = pyeit.eit.svd.SVD(self.mesh_obj, self.protocol_obj)
+        eit.setup(n=50, method="svd", perm=1, jac_normalized=True)
+        ds = eit.solve(self.v1, self.v0, normalize=True)
+        # evaluate
+        loc, ds_max, ds_sign = eit_loc_eval(ds, self.mesh_obj, mode="element")
+        # print(loc, ds_max, ds_sign)
+        loc = loc[: len(self.anomaly["center"])]
+        dist = np.linalg.norm(loc - self.anomaly["center"])
+
+        self.assertTrue(ds_sign == int(self.anomaly["sign"]))
+        # more tolerance on SVD
+        self.assertTrue(dist < 2 * self.anomaly["r"])
+
+    def test_greit(self):
+        """test GREIT"""
+        eit = pyeit.eit.greit.GREIT(self.mesh_obj, self.protocol_obj)
+        eit.setup(p=0.50, lamb=0.01, perm=1, jac_normalized=True)
+        ds = eit.solve(self.v1, self.v0, normalize=True)
+        x, y, ds = eit.mask_value(ds, mask_value=np.NAN)
+
+        # evaluate GREIT
+        loc = np.where(np.abs(ds) == np.nanmax(np.abs(ds)))
+        center = np.array([x[loc][0], y[loc][0]])
+        ds_sign = np.sign(ds[loc][0])
+        # print(loc, center, ds_sign)
+        dist = np.linalg.norm(center - self.anomaly["center"])
 
-        assert solver.B.shape[0] > 0
+        self.assertTrue(ds_sign == int(self.anomaly["sign"]))
+        self.assertTrue(dist < self.anomaly["r"])
 
 
 if __name__ == "__main__":