From f710fe76e031dd959ae11874c9c4cf85cc17ee96 Mon Sep 17 00:00:00 2001 From: DavidMetzIMT Date: Sat, 7 May 2022 13:26:35 +0200 Subject: [PATCH] Issue setup small fix + optional meas_pattern (#45) - fix check/init stimulation and permitivity - solve, solve_eit, compute_jac and compute_b_matrix allow topass None ex_mat and None perm so we need to set them in each scope of those methods - notes: they are afterwards not memory in forward (only in eitbase) - rename _get in check fro better readability - add option to pass meas_pattern - need to be tested and documented - using of kwargs - typing error and some formatting - fix boolean comparison of ndarray - make dafualt val of ex_mat dependet of n_el - Add tests - add the default call (without arg) of solve, solve_eit, compute_jac, compute _b - Add test for optional meas_pattern - in fem.py improve the error str --- pyeit/eit/base.py | 15 +++- pyeit/eit/fem.py | 187 +++++++++++++++++++++++++++++++++++++--------- tests/test_fem.py | 46 +++++++++++- 3 files changed, 205 insertions(+), 43 deletions(-) diff --git a/pyeit/eit/base.py b/pyeit/eit/base.py index 6d66348..aa81303 100644 --- a/pyeit/eit/base.py +++ b/pyeit/eit/base.py @@ -63,8 +63,11 @@ def __init__( Notes ----- - parser is required for your code to be compatible with - (a) simulation data set or (b) FMMU data set + - parser is required for your code to be compatible with + (a) simulation data set or (b) FMMU data set + - To pass a custom measurement pattern use the kwarg meas_pattern + pay attenteion that the meas_pattern should be an nd.array of shape + (n_exc, n_meas_per_exc, 2). If not TypeError will be raised. """ if ex_mat is None: ex_mat = eit_scan_lines(len(el_pos), 8) @@ -73,6 +76,7 @@ def __init__( # build forward solver self.fwd = Forward(mesh, el_pos) + self.meas_pattern = kwargs.pop("meas_pattern", None) # solving mesh structure self.mesh = mesh @@ -222,6 +226,7 @@ def _compute_jac_matrix( perm=perm if perm is not None else self.perm, parser=self.parser, normalize=self.jac_normalized and allow_jac_norm, + meas_pattern=self.meas_pattern, ) def _compute_b_matrix(self) -> np.ndarray: @@ -234,7 +239,11 @@ def _compute_b_matrix(self) -> np.ndarray: BP matrix """ return self.fwd.compute_b_matrix( - ex_mat=self.ex_mat, step=self.step, perm=self.perm, parser=self.parser + ex_mat=self.ex_mat, + step=self.step, + perm=self.perm, + parser=self.parser, + meas_pattern=self.meas_pattern, ) def _check_solver_is_ready(self) -> None: diff --git a/pyeit/eit/fem.py b/pyeit/eit/fem.py index 4fa49d5..86bf675 100644 --- a/pyeit/eit/fem.py +++ b/pyeit/eit/fem.py @@ -109,6 +109,8 @@ def solve( (e.g. [0,7] or np.array([0,7]) or ex_mat[0].ravel). In that case a simplified version of `f` with shape (n_pts,) """ + ex_mat = self._check_ex_mat(ex_mat) # check/init stimulation + perm = self._check_perm(perm) # check/init permitivity f = self._compute_potential_distribution(ex_mat=ex_mat, perm=perm) # case ex_line has been passed instead of ex_mat # we return simplified version of f with shape (n_pts,) @@ -122,6 +124,7 @@ def solve_eit( step: int = 1, perm: Union[int, float, np.ndarray] = None, parser: Union[str, list[str]] = None, + **kwargs, ) -> FwdResult: """ EIT simulation, generate forward v measurement @@ -148,31 +151,20 @@ def solve_eit( Foward results comprising v: np.ndarray simulated boundary voltage measurements; shape(n_exc, n_el) - """ - f = self._compute_potential_distribution(ex_mat=ex_mat, perm=perm) - # boundary measurements, subtract_row-voltages on electrodes - diff_op = voltage_meter(ex_mat, n_el=self.n_el, step=step, parser=parser) - return FwdResult(v=self._get_boundary_voltages(f, diff_op)) - def _get_boundary_voltages(self, f: np.ndarray, diff_op: np.ndarray) -> np.ndarray: + Note + ---- + To pass a custom measurement pattern use the kwarg meas_pattern + pay attenteion that the meas_pattern should be an nd.array of shape + (n_exc, n_meas_per_exc, 2). If not TypeError will be raised. """ - Compute boundary voltages from potential distribution - - Parameters - ---------- - f : np.ndarray - potential on nodes ; shape (n_exc, n_pts) - diff_op : np.ndarray - measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2) + ex_mat = self._check_ex_mat(ex_mat) # check/init stimulation + perm = self._check_perm(perm) # check/init permitivity + f = self._compute_potential_distribution(ex_mat, perm) + # boundary measurements, subtract_row-voltages on electrodes + diff_op = self._build_meas_pattern(ex_mat, self.n_el, step, parser, **kwargs) - Returns - ------- - np.ndarray - simulated boundary voltage measurements; shape(n_exc, n_el) - """ - f_el = f[:, self.el_pos] - v = subtract_row(f_el, diff_op) - return np.hstack(v) + return FwdResult(v=self._get_boundary_voltages(f, diff_op)) def compute_jac( self, @@ -181,6 +173,7 @@ def compute_jac( perm: Union[int, float, np.ndarray] = None, parser: Union[str, list[str]] = None, normalize: bool = False, + **kwargs, ) -> np.ndarray: """ Compute the Jacobian matrix @@ -213,11 +206,15 @@ def compute_jac( ----- - initial boundary voltage meas. extimation v0 can be accessed after computation through call fwd.v0 + - To pass a custom measurement pattern use the kwarg meas_pattern + pay attenteion that the meas_pattern should be an nd.array of shape + (n_exc, n_meas_per_exc, 2). If not TypeError will be raised. """ - f = self._compute_potential_distribution( - ex_mat=ex_mat, perm=perm, memory_4_jac=True - ) + ex_mat = self._check_ex_mat(ex_mat) # check/init stimulation + perm = self._check_perm(perm) # check/init permitivity + + f = self._compute_potential_distribution(ex_mat, perm, memory_4_jac=True) # Build Jacobian matrix column wise (element wise) # Je = Re*Ke*Ve = (nex3) * (3x3) * (3x1) @@ -235,7 +232,8 @@ def jac_init(jac, k): self._r_matrix = None # clear memory self._ke = None # clear memory - diff_op = voltage_meter(ex_mat, n_el=self.n_el, step=step, parser=parser) + diff_op = self._build_meas_pattern(ex_mat, self.n_el, step, parser, **kwargs) + jac = subtract_row(jac_i, diff_op) self.v0 = self._get_boundary_voltages(f, diff_op) jac = np.vstack(jac) @@ -250,6 +248,7 @@ def compute_b_matrix( step: int = 1, perm: Union[int, float, np.ndarray] = None, parser: Union[str, list[str]] = None, + **kwargs, ) -> np.ndarray: """ Compute back-projection mappings (smear matrix) @@ -270,18 +269,28 @@ def compute_b_matrix( parser: Union[str, list[str]], optional see voltage_meter for more details, by default `None`. + Returns ------- np.ndarray back-projection mappings (smear matrix); shape(n_exc, n_pts, 1), dtype= bool + + Note + ---- + To pass a custom measurement pattern use the kwarg meas_pattern + pay attenteion that the meas_pattern should be an nd.array of shape + (n_exc, n_meas_per_exc, 2). If not TypeError will be raised. """ - f = self._compute_potential_distribution(ex_mat=ex_mat, perm=perm) + ex_mat = self._check_ex_mat(ex_mat) # check/init stimulation + perm = self._check_perm(perm) # check/init permitivity + + f = self._compute_potential_distribution(ex_mat, perm) f_el = f[:, self.el_pos] # build bp projection matrix # 1. we can either smear at the center of elements, using # >> fe = np.mean(f[:, self.tri], axis=1) # 2. or, simply smear at the nodes using f - diff_op = voltage_meter(ex_mat, n_el=self.n_el, step=step, parser=parser) + diff_op = self._build_meas_pattern(ex_mat, self.n_el, step, parser, **kwargs) # set new to `False` to get smear-computation from ChabaneAmaury b_matrix = smear(f, f_el, diff_op, new=True) return np.vstack(b_matrix) @@ -300,6 +309,29 @@ def set_ref_el(self, val: int = None) -> None: val if val is not None and val not in self.el_pos else max(self.el_pos) + 1 ) + ############################################################################ + # Intern methods + ############################################################################ + def _get_boundary_voltages(self, f: np.ndarray, diff_op: np.ndarray) -> np.ndarray: + """ + Compute boundary voltages from potential distribution + + Parameters + ---------- + f : np.ndarray + potential on nodes ; shape (n_exc, n_pts) + diff_op : np.ndarray + measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2) + + Returns + ------- + np.ndarray + simulated boundary voltage measurements; shape(n_exc, n_el) + """ + f_el = f[:, self.el_pos] + v = subtract_row(f_el, diff_op) + return np.hstack(v) + def _compute_potential_distribution( self, ex_mat: np.ndarray, perm: np.ndarray, memory_4_jac: bool = False ) -> np.ndarray: @@ -327,9 +359,6 @@ def _compute_potential_distribution( potential on nodes ; shape (n_exc, n_pts) """ - ex_mat = self._get_ex_mat(ex_mat) # check/init stimulation - perm = self._get_perm(perm) # check/init permitivity - # 1. calculate local stiffness matrix (on each element) ke = calculate_ke(self.pts, self.tri) # 2. assemble to global K @@ -350,7 +379,91 @@ def _compute_potential_distribution( .reshape(b.shape[0:2]) ) - def _get_perm(self, perm: Union[int, float, np.ndarray] = None) -> np.ndarray: + def _build_meas_pattern( + self, + ex_mat: np.ndarray, + n_el: int = 16, + step: int = 1, + parser: Union[str, list[str]] = None, + **kwargs, + ) -> np.ndarray: + """ + Build the measurement pattern (subtract_row-voltage pairs [N, M]) + for all excitations on boundary electrodes. + + Note + ---- + To pass a custom measurement pattern use the kwarg meas_pattern + pay attenteion that the meas_pattern should be an nd.array of shape + (n_exc, n_meas_per_exc, 2). If not TypeError will be raised. + + Parameters + ---------- + ex_mat : np.ndarray + Nx2 array, [positive electrode, negative electrode]. ; shape (n_exc, 2) + (see "voltage_meter") + n_el : int, optional + number of total electrodes, by default 16 + (see "voltage_meter") + step : int, optional + measurement method, by default 1 + (see "voltage_meter") + parser : Union[str, list[str]], optional + parsing the format of each frame in measurement/file, by default None + (see "voltage_meter") + + Returns + ------- + np.ndarray + measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2) + + """ + meas_pattern = self._check_meas_pattern(ex_mat.shape[0], **kwargs) + if meas_pattern is not None: + return meas_pattern + return voltage_meter(ex_mat, n_el, step, parser) + + def _check_meas_pattern( + self, n_exc: int, meas_pattern: np.ndarray = None + ) -> np.ndarray: + """ + Check measurement pattern + + Parameters + ---------- + 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), by default None + if None (no meas_pattern has been passed) None is returned + + Returns + ------- + np.ndarray + measurements pattern / subtract_row pairs [N, M]; shape (n_exc, n_meas_per_exc, 2) + + Raises + ------ + TypeError + raised if meas_pattern is not a nd.array of shape (n_exc, : , 2) + """ + + if meas_pattern is None: + return None + + if not isinstance(meas_pattern, np.ndarray): + raise TypeError( + f"Wrong type of {meas_pattern=}, expected an ndarray; shape ({n_exc}, n_meas_per_exc, 2)" + ) + # test shape is something like (n_exc, :, 2) + if meas_pattern.ndim != 3 or meas_pattern.shape[::2] != (n_exc, 2): + raise TypeError( + f"Wrong shape of {meas_pattern=}: {meas_pattern.shape=}, expected an ndarray; shape ({n_exc}, n_meas_per_exc, 2)" + ) + + return meas_pattern + + def _check_perm(self, perm: Union[int, float, np.ndarray] = None) -> np.ndarray: """ Check/init the permittivity on element @@ -384,7 +497,7 @@ def _get_perm(self, perm: Union[int, float, np.ndarray] = None) -> np.ndarray: ) return perm - def _get_ex_mat(self, ex_mat: np.ndarray = None) -> np.ndarray: + def _check_ex_mat(self, ex_mat: np.ndarray = None) -> np.ndarray: """ Check/init stimulation @@ -392,8 +505,8 @@ def _get_ex_mat(self, ex_mat: np.ndarray = None) -> np.ndarray: ---------- ex_mat : np.ndarray, optional stimulation/excitation matrix, of shape (n_exc, 2), by default `None`. - If `None` initialize stimulation matrix for 16 electrode and - apposition mode (see function `eit_scan_lines(16, 8)`) + If `None` initialize stimulation matrix for n_el electrode and + adjacent mode (see function `eit_scan_lines`) If single stimulation (ex_line) is passed only a list of length 2 and np.ndarray of size 2 will be treated. @@ -409,8 +522,8 @@ def _get_ex_mat(self, ex_mat: np.ndarray = None) -> np.ndarray: or np.ndarray of shape (n_exc, 2) """ if ex_mat is None: - # initialize the scan lines for 16 electrodes (default: apposition) - ex_mat = eit_scan_lines(16, 8) + # initialize the scan lines for 16 electrodes (default: adjacent) + ex_mat = eit_scan_lines(self.n_el, 1) elif isinstance(ex_mat, list) and len(ex_mat) == 2: # case ex_line has been passed instead of ex_mat ex_mat = np.array([ex_mat]).reshape((1, 2)) # build a 2D array diff --git a/tests/test_fem.py b/tests/test_fem.py index d292a89..7eb5824 100644 --- a/tests/test_fem.py +++ b/tests/test_fem.py @@ -145,6 +145,9 @@ def test_solve(self): f= fwd.solve(ex_mat, perm=mesh["perm"]) self.assertTrue(np.allclose(f, f_truth)) + # test without passing any argument + f= fwd.solve() + self.assertTrue(isinstance(f, np.ndarray)) def test_solve_eit(self): """test solve using a simple mesh structure""" @@ -156,14 +159,21 @@ def test_solve_eit(self): fwd.set_ref_el(mesh["ref"]) ex_mat = np.array([[0, 1], [1, 0]]) # include voltage differences on driving electrodes - fwd = fwd.solve_eit(ex_mat, parser="meas_current") + res = fwd.solve_eit(ex_mat, parser="meas_current") vdiff_truth = f_truth[el_pos[1]] - f_truth[el_pos[0]] v_truth = vdiff_truth * np.array([1, -1, -1, 1]) - self.assertTrue(np.allclose(fwd.v, v_truth)) + self.assertTrue(np.allclose(res.v, v_truth)) + # test without passing any argument + res= fwd.solve_eit() + self.assertTrue(isinstance(res.v, np.ndarray)) + + # test passing meas_pattern + res= fwd.solve_eit(ex_mat, meas_pattern=np.array([[[0, 1]],[[1, 0]]])) + self.assertTrue(isinstance(res.v, np.ndarray)) def test_compute_jac(self): - """test solve using a simple mesh structure""" + """test compute_jac using a simple mesh structure""" #TODO @ liubenyuan please checkt this test # compute_jac return jac with the "subtract_row part" here you wanted to test only the jac_i get from old solve method # the jac_i_truth correspond to the jac_truth! @@ -179,6 +189,36 @@ def test_compute_jac(self): jac = fwd.compute_jac(ex_mat, perm=mesh["perm"], parser="meas_current" ) self.assertTrue(np.allclose(jac, jac_truth)) + # test without passing any argument + jac = fwd.compute_jac() + self.assertTrue(isinstance(jac, np.ndarray)) + + # test passing meas_pattern + jac = fwd.compute_jac(ex_mat, meas_pattern=np.array([[[0, 1]]])) + self.assertTrue(isinstance(jac, np.ndarray)) + + def test_compute_b_matrix(self): + """test compute_jac using a simple mesh structure""" + + mesh, el_pos = _mesh_obj() + b_truth = np.array([]) # TODO @ liubenyuan + + # testing solve + ex_mat = np.array([[0, 1]]) + fwd = pyeit.eit.fem.Forward(mesh, el_pos) + # fix ref to be exactly the one in mesh + fwd.set_ref_el(mesh["ref"]) + b = fwd.compute_b_matrix(ex_mat, perm=mesh["perm"], parser="meas_current" ) + + # self.assertTrue(np.allclose(b, b_truth)) # TODO @ liubenyuan + # test without passing any argument + b = fwd.compute_b_matrix() + self.assertTrue(isinstance(b, np.ndarray)) + + # test passing meas_pattern + b = fwd.compute_b_matrix(ex_mat, meas_pattern=np.array([[[0, 1]]])) + self.assertTrue(isinstance(b, np.ndarray)) + if __name__ == "__main__":