Skip to content

Commit

Permalink
Issue setup small fix + optional meas_pattern (#45)
Browse files Browse the repository at this point in the history
- 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
  • Loading branch information
davmetz authored May 7, 2022
1 parent 8c4e228 commit f710fe7
Show file tree
Hide file tree
Showing 3 changed files with 205 additions and 43 deletions.
15 changes: 12 additions & 3 deletions pyeit/eit/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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:
Expand Down
187 changes: 150 additions & 37 deletions pyeit/eit/fem.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,)
Expand All @@ -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
Expand All @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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:
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -384,16 +497,16 @@ 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
Parameters
----------
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.
Expand All @@ -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
Expand Down
Loading

0 comments on commit f710fe7

Please sign in to comment.