An educational Python package implementing Direct Sampling Methods (DSM) and Iterative Direct Sampling Methods (IDSM) for elliptic inverse problems, with comprehensive Jupyter notebook tutorials.
The recently developed Iterative Direct Sampling Methods (IDSM) by Jin, Wang, Zou, and Ito represent a significant advance in solving inverse problems like Electrical Impedance Tomography (EIT). They combine the speed of direct methods with the accuracy of iterative ones. However, the mathematical depth of the original papers can be a barrier. This project bridges this gap by providing a clear, well-documented, open-source educational implementation.
Target Audience: Graduate students and researchers new to inverse problems or direct sampling methods.
Core Philosophy: Clarity over speed. Code is simple, well-commented, and modular. Mathematical steps in the papers are explicitly linked to lines of code.
The primary model is the Electrical Impedance Tomography (EIT) problem on an ellipse
Given boundary measurements (Cauchy data)
The framework extends to the generalized model with a zeroth-order potential term (Paper 1, Examples 2–3):
which covers diffuse optical tomography (DOT,
| Method | Reference | Key Idea |
|---|---|---|
| DSM (Direct Sampling) | Paper 1, Section 2.2 | Single-pass indicator via Green's function correlation |
| IDSM (Full Data) | Paper 1, Algorithm 3.2 | Iterative refinement with regularized DtN map + quasi-Newton |
| IDSM (Partial Data) | Paper 3, Algorithm 5.1 | Data completion + HR-DtN + stabilization-correction |
| IDSM (Parabolic) | Paper 2, Algorithm 4.1 | Time-segmented inversion of moving inclusions via backward adjoint |
| Phaseless DSM + DSM-DL | Paper 4, Section 3-5 | Corrected phaseless data (Eq. 3.11–3.12) and U-Net for inhomogeneous + impenetrable scatterers |
- Python 3.8+
- NumPy, SciPy, Matplotlib, MeshPy, scikit-fem
# Create conda environment (recommended)
conda create -n IDSM python=3.10
conda activate IDSM
# Install dependencies
pip install -r requirements.txt
# Verify
python -c "import sys; sys.path.insert(0,'.'); from src.idsm import run_idsm; print('OK')"The recommended way to get started is through the Jupyter notebooks, which provide step-by-step theory and code:
cd IDSM-easy-to-begin/notebooks
jupyter notebook 01_forward_problem.ipynb # Phase 1: FEM, mesh, forward solverThe six notebooks are designed to be followed in order:
| Notebook | Content | Approx. Time |
|---|---|---|
01_forward_problem.ipynb |
Weak form, P1 FEM, saddle-point system, noise model | ~1 min |
02_classical_dsm.ipynb |
Laplace-Beltrami, DSM indicator, limitations | ~2 min |
03_iterative_dsm.ipynb |
Regularized DtN map, IDSM Algorithm 3.2, DFP/BFG | ~8 min |
04_comparative_study.ipynb |
DSM vs IDSM, partial data, ablation, noise sweep | ~15 min |
05_parabolic_idsm.ipynb |
Parabolic IDSM (Algorithm 4.1), FreeFEM-default Section 5 port, saved 05_*.png tutorial figures |
~10–20 min |
06_phaseless_dsmdl.ipynb |
Phaseless DSM (Section 5.1, BIE Dirichlet/Neumann + VIE Lippmann-Schwinger), DSM-DL with U-Net (Section 5.2 polygon/MNIST/mixed-circle), plus reference-code Python port | ~1 min from cache; full retrain scripts/run_phaseless_full_repro.py ~1 h |
Alternatively, you can use the package API directly:
from src.mesh import generate_elliptic_mesh
from src.forward_solver import make_conductivity_example1, generate_cauchy_data
from src.idsm import run_idsm
# Generate mesh and ground truth
mesh = generate_elliptic_mesh(n_boundary=256)
sigma_true, u_true = make_conductivity_example1(mesh)
# Generate synthetic Cauchy data with 10% noise
sources = [lambda x, y: x, lambda x, y: y]
data = generate_cauchy_data(mesh, sigma_true, sources, noise_level=0.1)
# Run IDSM (Algorithm 3.2)
history = run_idsm(mesh, data, n_iter=22, sigma_range=0.01,
problem_type="conductivity", lowrank_method="BFG")
# Result: history['sigma_final'] is the reconstructed conductivity (P0 array)
print(f"Residual: {history['residuals'][0]:.4e} -> {history['residuals'][-1]:.4e}")IDSM/
src/
__init__.py # Package exports
config.py # Centralized hyperparameters (RuntimeConfig, etc.)
mesh.py # Mesh generation, boundary extraction, coarse mesh
fem.py # P1 FEM public API (delegates to scikit-fem backend)
fem_skfem.py # scikit-fem backed FEM assembly (default backend)
fem_legacy.py # Hand-written P1 FEM (retained for regression testing)
forward_solver.py # Forward PDE solver, Cauchy data generation
dsm.py # Classical DSM (Laplace-Beltrami, indicator)
idsm.py # Full-data IDSM (Algorithm 3.2)
idsm_partial.py # Partial-data IDSM (Algorithm 5.1)
idsm_parabolic.py # Parabolic IDSM (Algorithm 4.1, Paper 2)
phaseless_scattering.py # Helmholtz forward (LU Lippmann-Schwinger) + DSM index (Eq. 3.11)
phaseless_bie.py # BIE Dirichlet/Neumann solvers for Section 5.1 impenetrable obstacles
phaseless_dsmdl.py # DSM-DL: datasets, U-Net (paper Fig.1), TV+window-SSIM loss, training/eval
phaseless_paper_figs.py # Shared assemblers used by the script and Notebook 06 for Fig.6-10
utils.py # Visualization, distance, IoU metrics
notebooks/
01_forward_problem.ipynb # Phase 1: FEM, mesh, forward solver
02_classical_dsm.ipynb # Phase 2: DSM baseline
03_iterative_dsm.ipynb # Phase 3: IDSM with DtN map
04_comparative_study.ipynb # Phase 4: Full comparison
05_parabolic_idsm.ipynb # Phase 5: Parabolic IDSM (moving inclusions)
06_phaseless_dsmdl.ipynb # Phase 6: Phaseless DSM + DSM-DL (Helmholtz, U-Net, Fig.2-10)
tests/
test_mesh.py # Mesh area, boundary, coarsening
test_fem.py # Stiffness symmetry, mass, Neumann solver
test_fem_regression.py # skfem vs legacy numerical agreement
test_forward.py # Forward solver, noise model
test_dsm.py # Eigenvalues, indicator
test_idsm.py # Box constraints, residual decrease
test_idsm_partial.py # Data completion, HR-DtN, lambda
test_utils.py # Distance, IoU, grid projection
test_config.py # Configuration defaults and env vars
test_phaseless_scattering.py # DSM corrected data shape, noise model, peak diagnostics
test_phaseless_dsmdl.py # U-Net forward shape, loss finite, dataset shape
reference/
Example1.edp ... Example5.edp # FreeFEM reference code
figures/ # Generated figures from notebooks
parabolic/ # Paper Section 5 multi-frame reproduction figures
06_phaseless/ # Phaseless DSM + DSM-DL paper figures (Fig.2-10)
results/
parabolic/ # Reproducible NPZ caches consumed by figures/05_parabolic
phaseless/ # full_summary.json, full_comparison.md, checkpoints/
scripts/
run_all_examples.py # Regenerate Paper Section 5 NPZ caches with ThFine/Th/ThCoarse
plot_parabolic_paper_figures.py # Regenerate figures/05_parabolic from NPZ caches
run_nb04_figures.py # Regenerate NB04 comparative-study figures (figures/04_comparative/04_*.png)
run_phaseless_full_repro.py # Phase 6: paper-scale Helmholtz reproduction (Fig.2-10 + checkpoints)
run_phaseless_forward_generate.py # Phase 6: pure Python port of ISP_forward/DataMnist.m
run_phaseless_port_repro.py # Phase 6: pure Python port of DSMDL_phaseless/main.py
requirements.txt
README.md
generate_elliptic_mesh(n_boundary)-- Create unstructured triangular mesh on the ellipsegenerate_sampling_grid(n_grid)-- Regular grid for DSM indicator evaluationgenerate_coarse_mesh(target_triangles)-- Coarse mesh for Paper 3 stabilizationfine_to_coarse_p0/coarse_to_fine_p0-- P0 inter-mesh projection
FEM public API, backed by scikit-fem (default) or hand-written legacy implementation (IDSM_FEM_LEGACY=1).
-
assemble_stiffness_matrix(mesh, sigma)-- Stiffness matrix$K_{ij} = \int \sigma \nabla\phi_j \cdot \nabla\phi_i$ -
assemble_mass_matrix(mesh, coeff)-- Mass matrix$M_{ij} = \int q,\phi_j \phi_i$ -
assemble_boundary_mass_matrix(mesh)-- Boundary mass$M_\Gamma$ -
assemble_partial_boundary_mass_matrix(mesh, gamma_d_node_mask)-- Split$M_\Gamma \to M_D + M_N$ (Paper 3) -
compute_boundary_normal_flux(mesh, sigma, y)-- Boundary normal flux$\sigma,\partial y/\partial n$ -
compute_boundary_normal_derivative(mesh, z, sigma_bg)-- Generic$\sigma_0\partial z/\partial n$ (geometry-independent) -
solve_neumann_system(K, b, B)-- Saddle-point solve$[K,B; B^T,0]$ -
solve_robin_system(mesh, A, alpha, v)-- Robin BVP for DtN map
-
run_idsm(mesh, cauchy_data, ...)-- Main IDSM loop (Algorithm 3.2) -
apply_regularized_dtn(mesh, v, A_op, alpha)-- Double Robin BVP (Eq. 3.20) -
compute_p0_gradient(mesh, w_list, y_list)-- Gradient$\zeta_k$ (Eq. 3.17) -
LowRankPreconditioner-- DFP/BFG quasi-Newton (Eq. 3.14-3.15)
-
run_idsm_partial(mesh, cauchy_data, gamma_d_info, ...)-- Paper 3 Algorithm 5.1 -
define_accessible_boundary(mesh, theta_range)-- Define$\Gamma_D$ from angle range -
complete_data(y_data, y_current, mask)-- Data completion (Eq. 4.1) -
apply_hr_dtn(mesh, v, A_op, alpha_d, alpha_n, ...)-- Heterogeneous DtN (Eq. 4.2) -
StabilizedLowRankResolver-- Stabilization-correction scheme (Eq. 4.10-4.16)
-
run_idsm_parabolic(coarse_mesh, fine_mesh, cfg, c_func, v_func, solve_mesh=..., ...)-- Algorithm 4.1, time-segmented IDSM for moving inclusions.fine_meshis the data meshThFine,solve_meshis the P1 PDE meshTh, andcoarse_meshis the P0 coefficient meshThCoarse; the function auto-dispatches to σ/V double-field or U-recovery single-field path viacfg.model -
solve_forward_segment(...)-- Crank–Nicolson forward solver per segment (linear σ + V) -
solve_forward_segment_nonlinear(...)-- Newton + Crank–Nicolson solver for §5.3 nonlinear$|y|y\cdot U$ source -
iterate_segment_nonlinear(...)-- U-recovery inner loop (single-field P0 reconstruction projected to$[u_A, 2u_B]$ ) -
edp_cfg_example_5_{1..5}(noise=...)/paper_cfg_example_5_{1..5}(noise=...)-- Section 5 Examples 5.1–5.5 configurations. Both currently keep the FreeFEM/reference-program numerical defaults;paper_cfgsupplies paper-style noise choices whileedp_cfgpreserves each reference script default. These are not a strict paper-text parameter profile when the paper and.edpdefaults differ (e.g. damping, time step, and some low-rank choices). -
solve_adjoint_segment(...)-- Backward adjoint PDE (Eq. 4.6) with terminal condition -
synthesize_full_forward(...)-- Reference fine-mesh forward pass that produces noisy boundary measurements -
c_func_example_5_{1..5}/v_func_example_5_{1..5}/u_func_example_5_3-- Ground-truth coefficient fields used both for synthesizing data and for IoU evaluation -
trajectory_example_5_{1..5}/radius_example_5_{1,4,5}-- Inclusion centre/radius trajectories per example -
ParabolicConfig-- Parabolic IDSM hyperparameters (defaults matchparabolic_*.edp)
Notebook 05 produces tutorial figures directly under figures/05_parabolic/05_*.png. The FreeFEM-default multi-frame plots are generated from three-mesh Python caches:
# Main paper-noise caches, including Example 5.1 at 5% and 10% noise.
python scripts/run_all_examples.py --mode paper --include-ex51-n10 --mesh-mode edp --quiet
# Figure 5.3 and 5.4 read the FreeFEM-default cache names used by the plotter.
python scripts/run_all_examples.py --example 5.3 --mode edp --mesh-mode edp --quiet
python scripts/run_all_examples.py --example 5.4 --mode edp --mesh-mode edp --quiet
python scripts/plot_parabolic_paper_figures.pyThe run commands write numerical caches to results/parabolic/*.npz; the plot command reads those caches and refreshes figures/05_parabolic/*.png plus figures/05_parabolic/table1.txt. The NPZ files are not source code, but they are kept so the FreeFEM-default figures can be inspected without rerunning the full parabolic experiment every time.
-
PhaselessBatchSimulator-- Hybrid forward solver: medium scatterers go through a batched Lippmann-Schwinger LU solve using the paper Eq. (2.4) coefficient contrast$k^2(n-1)$ (torch.linalg.lu_factorcached per sample, reused across incidences), while sound-soft scatterers ($u=0$ on$\partial D$ , paper Eq. 2.3) are dispatched to the Dirichlet boundary-integral equation inphaseless_bie.py. Mixed-circle scenes use a Born-superposition approximation (medium VIE field + per-soft-scatterer BIE field). Bothcompute_dsm_inputs(labels, ...)(legacy complex-$n$ surrogate) andcompute_dsm_inputs_with_meta(labels, metas, ...)(strict BIE+VIE Born) yield either Eq. (3.11) phaseless or Eq. (3.1) phased DSM indices. -
run_example_dsm_paper(example_key, noise_level, n_incident, ...)-- Section 5.1 dispatcher: VIE for medium square / close squares / ring, and BIE Dirichlet for sound-soft squares. -
add_phaseless_noise,corrected_phaseless_data,compute_phaseless_dsm_indicator-- Building blocks for the paper's Eq. (3.11)-(3.12).
-
boundary_circle,boundary_square,stack_boundaries-- Boundary discretizations used by Section 5.1. -
solve_sound_soft-- Single-layer Dirichlet BIE$S,\varphi = -u_{\text{inc}}$ . -
solve_sound_hard-- Indirect single-layer Neumann BIE$(-I/2 + K^T),\varphi = -\partial_n u_{\text{inc}}$ . -
evaluate_total_at-- Reconstruct the total field at the receiver ring from the boundary density.
-
DatasetConfig,TrainingConfig-- Strict paper-protocol settings (wavelength, receiver geometry,$N_i$ ,$n_{\text{soft}}$ , batch / epochs / LR schedule). -
build_dataset/build_dataset_with_norm-- Generate the polygon, MNIST (official train/test split), and mixed-circle labels and the matching DSM indicator inputs. -
UNetDSMDL(in_channels, out_channels=1, base_channels=64, depth=4)-- U-Net mirroring paper Fig.1 (64-128-256-512-1024). -
dsmdl_loss-- Eq. (4.1)$\text{MSE} + 0.5,\text{TV}(X) + 0.5(1-\text{SSIM}_{11\times 11}(X,Y))$ . -
BalancedPolygonBatchSampler,make_polygon_dataloader-- 5+5 medium/soft per-batch sampler (paper §5.2.1.1). -
train_unet,predict,threshold_to_classes,pixel_accuracy,relative_l2_error-- Training loop and Eq. (5.3)/(5.4) metrics. -
save_model_checkpoint,load_model_checkpoint-- Persist trained networks for fast reuse.
- Shared assemblers for paper Fig.6–Fig.10. Both
scripts/run_phaseless_full_repro.pyandnotebooks/06_phaseless_dsmdl.ipynbimport from here, so the layout stays identical between the offline run and the inline notebook view. load_case_from_checkpoint(path, in_channels, device)-- Rehydrate a model + dataset metadata from.ptfor figure rendering.render_and_save_all(ckpt_dir, out_fig, seed, device)-- Convenience entry that regenerates all DSM-DL paper figures from checkpoints.
The strict run is driven by a single script and produces every paper output (Fig.2–Fig.10, Table 1, Table 2, Section 5.2.3 metrics, and the trained checkpoints):
# Strict paper-protocol reproduction. ~60 minutes on an A100.
conda run -n IDSM python scripts/run_phaseless_full_repro.pyThe reference-code path is also available in pure Python:
# Optional: regenerate forward data without MATLAB.
conda run -n IDSM python scripts/run_phaseless_forward_generate.py --n-samples 200
# Reproduce DSMDL_phaseless/main.py behavior (U_Net3Ab + legacy loss weights).
conda run -n IDSM python scripts/run_phaseless_port_repro.pyIt writes:
results/phaseless/full_summary.jsonandresults/phaseless/full_comparison.md(Table 1, Table 2, mixed-circle metrics vs paper),results/phaseless/checkpoints/{polygon_Ni1, polygon_Ni4, mnist_Ni4, mnist_Ni16, mixed_circle_Ni10}.pt,figures/06_phaseless/{full_fig2_4_examples1_3, full_fig5_ring_ni_sweep, fig6_polygon_recon, fig7_mnist_recon, fig8_chinese_recon, fig9_austria_recon, fig10_mixed_recon}.png.
Notebook 06 §6 then renders Fig.6–Fig.10 inline from those checkpoints in well under a minute, without retraining.
-
compute_dsm_indicator(mesh, cauchy_data, gamma, ...)-- DSM indicator$\eta(x)$ (Eq. 2.8) -
discretize_laplace_beltrami(mesh, gamma)--$(-\Delta_\Gamma)^\gamma$ eigendecomposition -
LaplaceBeltramiOperator-- Discrete$(-\Delta_\Gamma)^\gamma$ operator withapply()method -
compute_scattering_data(cauchy_data)-- Scattering$y_d^s = y_\emptyset - y_d$ -
compute_dsm_numerator(mesh, scatter, lb_op)-- Auxiliary PDE solve (Eq. 2.9) -
compute_dsm_denominator_integral(mesh, points)-- FreeFEM-style normalization (Eq. 2.10)
-
solve_forward(mesh, sigma, f_func)-- EIT forward solve:$\nabla\cdot(\sigma\nabla y)=0$ -
solve_forward_general(mesh, sigma, potential, f_func)-- Generalized with zeroth-order term -
generate_cauchy_data(mesh, sigma, sources, noise_level)-- Noisy Cauchy data pairs -
generate_cauchy_data_general(mesh, sigma, potential, sources, noise_level)-- DOT Cauchy data -
make_conductivity_example1(mesh)-- Example 1 (two insulating squares) -
make_conductivity_conductive(mesh)-- Conductive variant ($\sigma=3.0$ ) -
make_conductivity_single(mesh)-- Single circular inclusion -
make_double_example2(mesh)-- Example 2 (simultaneous$\sigma$ +$v$ , double type) -
make_potential_example3(mesh)-- Example 3 (potential-only, DOT)
compute_iou(u_true, u_recon, mesh)-- Area-matched Intersection over Uniondistance_to_boundary(mesh, points)-- Min distance from points to boundary edges
RuntimeConfig-- GPU/seed/backend settings (from environment variables)MeshConfig-- Mesh resolution parametersFullIDSMConfig-- Full-data IDSM hyperparameters (Algorithm 3.2)PartialIDSMConfig-- Partial-data IDSM hyperparameters (Algorithm 5.1)DoubleIDSMConfig-- Example 2 (double type) hyperparameters (FreeFEM Example2.edp)Notebook01Config...Notebook04Config-- Per-notebook configuration dataclasses
| Notebook | Phase | Content |
|---|---|---|
01_forward_problem.ipynb |
1 | Weak form, P1 FEM, saddle-point system, mesh convergence, noise model |
02_classical_dsm.ipynb |
2 | Laplace-Beltrami eigenproblem, DSM indicator, gamma parameter, limitations |
03_iterative_dsm.ipynb |
3 | Regularized DtN map, Algorithm 3.2, DFP/BFG, convergence analysis |
04_comparative_study.ipynb |
4 | DSM vs IDSM, partial data, damping factor, ablation, noise sweep |
05_parabolic_idsm.ipynb |
5 | Crank-Nicolson forward, backward adjoint (Eq. 4.6), Algorithm 4.1, FreeFEM-default live runs for Examples 5.1–5.5 (merging / mixed / nonlinear / fading / diminishing) |
06_phaseless_dsmdl.ipynb |
6 | Helmholtz forward, phaseless DSM (Eq. 3.11–3.12), BIE Dirichlet/Neumann for Section 5.1, U-Net DSM-DL (Eq. 4.1, Fig.1), Fig.2–10 reproduction with checkpoint cache |
All hyperparameters are centralized in src/config.py:
from IDSM.src.config import Notebook04Config, RuntimeConfig
cfg = Notebook04Config()
cfg.full.alpha = 1.0 # Robin regularization
cfg.full.n_iter = 22 # Iterations (FreeFEM: storeNum=22)
cfg.partial.alpha_d = 0.05 # Paper 3 Table 1
cfg.partial.alpha_n = 2.0GPU toggle (currently CPU-only due to sparse solver constraints):
runtime = RuntimeConfig(use_gpu=True) # Will warn and fall back to CPUcd IDSM
pytest tests/ -v # Default: scikit-fem backend (85 tests)
IDSM_FEM_LEGACY=1 pytest tests/ -v # Legacy hand-written FEM backendQ: Why scikit-fem instead of FEniCSx?
A: scikit-fem is a pure Python library (pip-installable, no compiled dependencies) that covers all P1 triangular FEM operations needed by this project. FEniCSx requires PETSc/MPI compilation and has heavier installation requirements, which creates unnecessary barriers for an educational package. Both are mathematically equivalent for the P1 case; our regression tests verify numerical agreement to machine precision.
Q: Can I switch back to the hand-written FEM?
A: Yes. Set IDSM_FEM_LEGACY=1 as an environment variable. The adapter layer in fem.py will delegate to fem_legacy.py instead of fem_skfem.py. All tests pass with both backends.
-
K. Ito, B. Jin, F. Wang, J. Zou, "Iterative direct sampling method for elliptic inverse problems with limited Cauchy data," SIAM J. Imaging Sci. 18(2), 2025. arXiv:2503.00423
-
B. Jin, F. Wang, J. Zou, "A direct sampling method for simultaneously recovering inhomogeneous inclusions of different nature," SIAM J. Sci. Comput. 43(3), 2021. arXiv:2005.05499
-
B. Jin, F. Wang, J. Zou, "A stable iterative direct sampling method for elliptic inverse problems with partial Cauchy data," J. Comput. Phys. 550, 2026. arXiv:2511.08171
-
B. Jin, F. Wang, J. Zou, "An iterative direct sampling method for reconstructing moving inhomogeneities in parabolic problems," 2025. arXiv:2511.08197
-
J. Ning, F. Han, J. Zou, "A direct sampling method and its integration with deep learning for inverse scattering problems with phaseless data," SIAM J. Sci. Comput., 2025. arXiv:2403.02584
FreeFEM reference code:
- Elliptic: github.com/RaulWangfr/IDSM-elliptic
- Parabolic: github.com/RaulWangfr/IDSM-parabolic
- Phaseless DSM-DL: no public reference repository identified; Phase 6 is a paper-faithful re-implementation.
This project is for educational purposes.