Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Issue setup small fix + optional meas_pattern #45

Merged
merged 10 commits into from
May 7, 2022
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