Electrode Models & Forward Problem

AbstractElectrodeModel

How current enters and voltage is measured on the boundary (continuum, point, gap or complete electrode model). A ForwardModel turns an electrode model on a discretization into injection and measurement matrices.

ContinuumModel()

Continuum model: every boundary dof acts as an electrode. Current patterns are boundary current densities (one value per boundary dof, disc.boundary_dofs order), voltages are the boundary values of the potential.

PointElectrodeModel(points; measure = points)

Point electrodes: current I_ℓ enters at the boundary dof closest to points[ℓ], voltages are read at the boundary dofs closest to measure. Injection and measurement points may differ.

GapModel(electrodes; measure = electrodes)

Gap model: current I_ℓ is uniformly distributed over electrode electrodes[ℓ] (a vector of boundary facets in the back end's representation, e.g. FacetIndex for Ferrite), zero in the gaps. The voltage of a measurement electrode is the mean potential on it. Injection and measurement electrodes may differ (measure).

CompleteElectrodeModel(electrodes, z; measure = 1:length(electrodes))

Complete electrode model with contact impedances z (scalar or one per electrode). The electrode voltages U are unknowns of the forward problem (grounded by Σ U_ℓ = 0); measure selects the electrodes whose voltages are measured (e.g. only electrodes that carry no current).

angular_electrodes(disc, L; coverage = 0.5, offset = 0.0, center = nothing, angles = nothing)

L electrodes (sets of boundary facets in the back end's representation) around center, centred at the angles offset + 2π(ℓ-1)/L or at angles, covering the fraction coverage of the boundary. Back end contract.

electrode_length(disc, electrode)

Length (area in 3D) of an electrode. Back end contract.

AbstractForwardModel

A discretized EIT forward problem: system matrix A(σ) = A₀ + Σₐ σₐ ∂A/∂σₐ, injection matrix P, measurement matrix Q, null space and grounding of the current-driven problem, and the structure of the voltage-driven (Dirichlet) problem.

ForwardModel(disc, model::AbstractElectrodeModel)

Discrete forward model of the electrode model model on the discretization disc.

Fields: n (system size; n_u for most models, n_u + L for the complete electrode model), n_u, n_σ, A (system matrix, updated by system_matrix!), A₀ (constant part of its stored values or nothing), tensor (ConductivityTensor), P (injection, n × n_inject), Q (measurement, n_measure × n), nullspace, grounding (the functional w with wᵀx = 0 for the current-driven solution), measure_weights (weights of the mean removed from measured voltages), dirichlet_dofs, free_dofs, E (Dirichlet expansion), C (current representation), angles (angular positions of the injection sites, used by trigonometric_patterns).

See forward_neumann and forward_dirichlet.

system_matrix!(fm, σ)

Update fm.A = A₀ + L(σ) (and its Dirichlet block) for the conductivity coefficients σ; returns fm.A. Non-positive coefficients throw an InfeasibleConductivityError, which the optimizers treat as a rejected trial point.

n_inject(fm), n_measure(fm), n_control(fm)

Number of injection sites (columns of P), measurements (rows of Q) and Dirichlet controls (columns of E; equals n_inject).

trigonometric_patterns(fm, K)

n_inject × 2K matrix of trigonometric patterns cos(kθ), sin(kθ), k = 1..K, at the angular positions of the injection sites, shifted so that each pattern injects zero net current (for the continuum model: ∫ g ds = 0). Also usable as voltage patterns.

forward_neumann(fm, σ, currents; solver = DirectSolver()) -> (voltages, X)

Current-driven forward problem: solve A(σ) X = P currents (grounded) and measure voltages = Q X. currents is a vector or an n_inject × s matrix of patterns; X holds the full states (n × s).

forward_dirichlet(fm, σ, voltages; solver = DirectSolver()) -> (currents, X)

Voltage-driven forward problem: prescribe voltages (n_control × s) on the Dirichlet dofs (boundary nodes, electrode nodes or CEM electrode voltages), solve for the rest and return the currents in the representation of the injection (so that it inverts forward_neumann for consistent data).

reground(V, w)
reground(fm, V, grounding)

Shift every column of V by a constant so that wᵀv = 0 (V - 1 (wᵀV)/(wᵀ1)). Groundings differ only by such shifts, so data can be moved between them and back without loss.

With a forward model: grounding = :integral uses fm.measure_weights (zero boundary mean for the continuum model), :nodal uses unit weights (zero sum of the measured values). For the electrode models both are the plain mean over the measured electrodes.

pattern_svd(disc, fm, currents, voltages; metric = :euclidean, noise = nothing, reference = nothing)
pattern_svd(currents, voltages, Mi, Mv; noise = nothing, reference = nothing)

New pairs of boundary data from measured pairs (currents, voltages) (n_inject × s, n_measure × s): currents * C, voltages * C for the s × s matrix C that makes the new currents orthonormal in the inner product Mi and the new voltages orthogonal in Mv, sorted by decreasing singular value (the norms of the new voltages). The voltages are first regrounded consistently with the metric (Mv-weighted mean zero). Returns (; currents, voltages, values, Mi, Mv, combination, projection, noise, noise_levels, reference) with the combination matrix C (new pairs = old pairs · C) and the measurement modes projection (one row per mode, orthonormal in the Mv⁻¹ inner product, ordered like the patterns): projection * (voltages - reference) is the diagonal matrix of values. For a two-sided truncation of the data, keep the leading rows as a ProjectedMisfit.

With a reference (the voltages of a reference conductivity for the same currents, or, with a forward model, the reference conductivity itself), the SVD is taken of the difference voltages - reference (Isaacson's distinguishability: the new currents are the patterns that best distinguish the unknown from the reference), values are the singular values of the difference, and the returned voltages and reference are the measured and reference data rotated by the same C, so that the pairs remain valid data for any objective.

With a noise model of the measured voltages (GaussianNoise, RelativeGaussianNoise; the currents are taken as exact), noise is the noise of the new voltages, a GaussianNoise with one standard deviation per entry (the rotation mixes patterns of different noise; use it for discrepancy_target), and noise_levels the expected Mv-norm of the noise in every new pattern, the scale to compare values with (see truncate_patterns). Both are nothing without a noise model.

For relative noise the absolute data have a signal-to-noise ratio of about 1/δ in every pattern, so truncation at the noise level removes nothing: it is the difference to a reference whose singular values decay below the noise.

metric:

• :euclidean: Mi = Mv = I (nodal/electrode values; matches the nodal-sum ground),
• :L2: the L²(Γ) inner products, i.e. the boundary mass matrix for the continuum model (current densities and voltages), and for gap and complete electrode models |e_ℓ| for electrode voltages and 1/|e_ℓ| for electrode currents (their densities I_ℓ/|e_ℓ|). Mesh independent; not defined for point electrodes.
• a tuple (Mi, Mv) of symmetric positive definite matrices.

On uniform boundary meshes the two built-in metrics give the same patterns up to scaling. For the continuum model with the L² metric the values approximate the singular values of the Neumann-to-Dirichlet map on the span of the input patterns.

truncate_patterns(p; τ = 2, measurements = nothing)
truncate_patterns(p, K; measurements = nothing)

The leading pairs of a pattern_svd result p: the first K, or those before the first pattern whose singular value drops to τ times its noise level (needs pattern_svd(...; noise); at least one pair is kept). The trailing pairs mostly carry noise; discarding them regularises in data space, without an assumption on the conductivity. Directions that carry only noise do not have singular values below the noise level but at it (the SVD of noisy data has a noise floor), so τ must lie clearly above 1. Returns a named tuple with the same fields, restricted to the retained pairs. With measurements = M, only the leading M measurement modes are kept in projection; use them as misfit = ProjectedMisfit(t.projection) for K M residuals.