Gbrunin/analyzer refactoring

.gitignore

@@ -54,6 +54,7 @@ coverage.xml
 .pytest_cache/
 cover/
 prof/
+.ruff_cache/
 
 # Translations
 *.mo

CONTRIBUTING.md

@@ -87,7 +87,7 @@ parsers' handling of orthogonal cells and periodic boundary conditions.
 ## Adding a new contact-angle method
 
 Subclass `BaseTrajectoryAnalyzer`
-([src/wetting_angle_kit/analysis/analyzer.py](src/wetting_angle_kit/analysis/analyzer.py))
+([src/wetting_angle_kit/analysis/_base.py](src/wetting_angle_kit/analysis/_base.py))
 and add an integration test in `tests/test_analysis/` that
 exercises the method on one of the fixture trajectories.
 

README.md

@@ -8,21 +8,54 @@
 [![License: BSD 3-Clause](https://img.shields.io/badge/License-BSD_3--Clause-blue.svg)](LICENSE)
 [![Documentation](https://img.shields.io/badge/docs-matgenix.github.io-blue)](https://matgenix.github.io/wetting-angle-kit)
 
-wetting-angle-kit provides modular tools to parse MD trajectories (LAMMPS dump, XYZ, ASE) and compute droplet contact angles using two complementary approaches:
+wetting-angle-kit parses MD trajectories (LAMMPS dump, XYZ, ASE) and computes the contact angle of a droplet sitting on a planar wall. The package follows the same conceptual recipe every method uses — extract the liquid-vapor interface from atom positions, decide where the wall plane sits, fit a geometric shape, read off the angle from the shape/wall intersection — but exposes each step as a swappable component so users can match the method to their system.
 
-1. Slicing Method (per-frame circle fit) – robust against transient shape changes.
-2. Binning Density Method – averages frames into a density field for a single representative angle.
+## How the methods are built
 
-The documentation is available [here](https://matgenix.github.io/wetting-angle-kit), you can find examples and tutorials.
+### Interface extraction: how do we turn atoms into a surface?
+
+The liquid-vapor interface isn't a sharp surface in an MD simulation — the density drops smoothly over ~1 Å. Two extraction strategies recover a clean set of interface points from the noisy atom cloud:
+
+The package exposes two orthogonal strategy axes for interface extraction. A `SpaceSampling` decides *where* density is evaluated; a `DensityEstimator` decides *how*. An `InterfaceExtractor` composes one of each.
+
+- **Sampling: `SpaceSampling.rays(...)`** emits a fan of rays from the droplet centre of mass; the interface along each ray is the half-density point of a 1D tanh fit. The fan layout is azimuthal slices in the `(x, z)` plane (for a per-slice fit) or a Fibonacci sphere of directions (for a whole-shape fit).
+- **Sampling: `SpaceSampling.grid(...)`** builds a fixed-cell grid in space; the interface is the iso-density contour at the half-bulk level, traced via marching squares (slicing mode) or marching cubes (whole mode). In slicing mode the grid iterates per slice (per azimuthal angle for spherical droplets, per axial step for cylinder droplets), so the slicing fit sees one `(s, z)` contour per slice and can expose per-slice asymmetry. Closer to the "average over many frames" intuition than ray fans; works well when atom statistics are limited per frame.
+- **Density: `DensityEstimator.gaussian(density_sigma=…)`** is a 3D Gaussian KDE (smooth, no per-cell Poisson noise).
+- **Density: `DensityEstimator.binning(bin_width=…)`** is a 3D top-hat histogram (cheap; bin_width required only for the rays sampling, where it sets the pointwise kernel size).
+
+Any sampling × any density is a valid extractor.
+
+### Surface fitting: what geometric shape do we fit to those points?
+
+- **Slicing fit** — independently fits an algebraic circle in each slice's `(x, z)` plane, then averages the per-slice contact angles. Good when the droplet might be slightly non-spherical: the per-slice scatter naturally reports a `±σ` band.
+- **Whole fit** — fits a single sphere (spherical droplet) or cylinder (cylindrical droplet) to the entire 3D interface shell. Uses the algebraic Taubin method, plus optional bootstrap resampling to put an uncertainty on the recovered angle.
+- **Coupled fit** (joint approach) — a 7-parameter (2D) or 9-parameter (3D) hyperbolic-tangent density model that solves "where is the interface", "where is the wall plane", and "what's the cap geometry" in one nonlinear least-squares fit on a density field. The per-cell density is computed by a pluggable `DensityEstimator` strategy: a top-hat histogram (`DensityEstimator.binning()`, default) or a 3D Gaussian KDE evaluated at the cell centres (`DensityEstimator.gaussian(density_sigma=…)`); the KDE variant trades a small constant cost for a smooth, Poisson-noise-free density. Statistically efficient when you pool many frames per batch.
+
+### Wall detection: where is the wall plane?
+
+The contact angle is measured at the cap–wall intersection, so the wall plane has to be located explicitly:
+
+- `min_plus_offset`: derive the wall from the interface itself (lowest interface point + offset). Works for slicing geometries and full-sphere ray fans, where the interface points reach the wall.
+- `from_atoms`: read the actual wall atom positions from the trajectory and place the wall at the mean of the top atomic layer. Most physically faithful when the simulation explicitly contains substrate atoms.
+- `explicit`: caller supplies the wall z directly — useful when the wall position is known a priori from the simulation setup.
+
+### Frame batching: per-frame angle or pooled batch?
+
+The `TemporalAggregator` groups trajectory frames into batches before fitting. `batch_size=1` runs the full pipeline once per frame (giving you an angle vs time curve); `batch_size=N` pools `N` frames together and fits one angle per pool (more atoms per fit → less noise, less time resolution); `batch_size=-1` pools everything into a single batch.
+
+## Two top-level entry points
+
+1. **`TrajectoryAnalyzer`** — composes the four strategies above (`InterfaceExtractor` × `SurfaceFitter` × `WallDetector` × `TemporalAggregator`). Use it when you want per-frame time resolution or when you want to mix-and-match approaches (e.g. ray-fan extractor + whole-fit + explicit wall + 5-frame batches).
+2. **`CoupledFit2DAnalyzer` / `CoupledFit3DAnalyzer`** — the joint-fit alternative. One robust angle per pooled batch via the hyperbolic-tangent density model. The per-cell density estimator is pluggable (`DensityEstimator.binning()` or `DensityEstimator.gaussian(...)`). Best when you have many frames and don't need per-frame time resolution.
+
+The documentation is available [here](https://matgenix.github.io/wetting-angle-kit), with worked examples and tutorials.
 
 ## Installation
 
 ### Prerequisites
 
-Before installing wetting-angle-kit, ensure you have the following prerequisites:
-
-1. **Python 3.10 or higher**: Make sure you have Python 3.10 or higher installed on your system.
-2. **Conda**: Ensure you have Conda installed. If not, you can install it from [here](https://docs.conda.io/en/latest/miniconda.html).
+Before installing wetting-angle-kit, ensure you have **Python 3.10 or higher**
+installed on your system.
 
 Core (only to analyse simple xyz trajectories):
 
@@ -45,45 +78,69 @@ pip install wetting-angle-kit[all]
 
 #### Install OVITO
 
-OVITO must be installed first in the conda environment and using the following Conda command:
+OVITO must be installed first using pip:
 
 ```sh
-conda install --strict-channel-priority -c https://conda.ovito.org -c conda-forge ovito=3.11.3
+pip install  ovito==3.11.3
 ```
 
 ## Quick Start
 
 
 ```python
 from wetting_angle_kit.analysis import (
-    BinningTrajectoryAnalyzer,
-    SlicingTrajectoryAnalyzer,
+    CoupledFit2DAnalyzer,
+    DensityEstimator,
+    InterfaceExtractor,
+    SpaceSampling,
+    SurfaceFitter,
+    TrajectoryAnalyzer,
+    WallDetector,
 )
+from wetting_angle_kit.analysis.temporal import TemporalAggregator
 from wetting_angle_kit.parsers import XYZParser, XYZWaterFinder
 
 trajectory_file = "trajectory.xyz"
 
-# Identify water oxygen atoms by neighbor count. ``particle_type_wall``
-# lists the symbols of the substrate atoms so they are excluded.
-finder = XYZWaterFinder(trajectory_file, particle_type_wall=["C"])
+# Identify water oxygen atoms by neighbour count.
+finder = XYZWaterFinder(trajectory_file)
 oxygen_ids = finder.get_water_oxygen_indices(frame_index=0)
 
 parser = XYZParser(trajectory_file)
 
-slicing = SlicingTrajectoryAnalyzer(
-    parser,
+# --- Composable pipeline (per-frame slicing-fit angles) ---
+slicing = TrajectoryAnalyzer(
+    parser=parser,
     atom_indices=oxygen_ids,
     droplet_geometry="spherical",
-    delta_gamma=5,
+    interface_extractor=InterfaceExtractor(
+        sampling=SpaceSampling.rays(
+            delta_azimuthal=5.0,  # 5° between slicing planes
+            delta_polar=8.0,
+        ),
+        density=DensityEstimator.gaussian(),
+    ),
+    surface_fitter=SurfaceFitter.slicing(surface_filter_offset=2.0),
+    wall_detector=WallDetector.min_plus_offset(offset=0.0),
+    temporal_aggregator=TemporalAggregator(batch_size=1),  # one angle per frame
 )
-results = slicing.analyze(frame_range=range(0, 50))
+results = slicing.analyze(range(0, 24))
 print(results.mean_angle, results.std_angle)
 
-binning = BinningTrajectoryAnalyzer(
-    parser,
+# --- Joint coupled-fit (one robust angle over a pooled batch) ---
+coupled_fit = CoupledFit2DAnalyzer(
+    parser=parser,
     atom_indices=oxygen_ids,
     droplet_geometry="spherical",
+    binning_params={
+        "xi_0": 0.0, "xi_f": 70.0, "bin_width_x": 2.0,
+        "zi_0": 0.0, "zi_f": 70.0, "bin_width_z": 2.0,
+    },
+    # Default: histogram density. Swap in `DensityEstimator.gaussian(
+    # density_sigma=2.5)` for a smooth Gaussian-KDE density field —
+    # useful on per-frame batches or sparse systems.
+    density_estimator=DensityEstimator.binning(),
 )
-results_binning = binning.analyze(frame_range=range(0, 200))
-print(results_binning.mean_angle, results_binning.std_angle)
+results_coupled_fit = coupled_fit.analyze(range(0, 200))
+print(results_coupled_fit.mean_angle, results_coupled_fit.std_angle)
 ```

docs/examples/coupled_fit_3d_ca.py

@@ -0,0 +1,82 @@
+"""Coupled-fit 3D contact-angle example.
+
+Runs the 3D coupled hyperbolic-tangent fit on a full ``(xi, yi, zi)``
+density grid via :class:`CoupledFit3DAnalyzer`. The nine-parameter
+model fits the spherical-cap interface and the wall plane
+simultaneously, recovering a single robust angle per pooled batch.
+
+Only spherical droplets are supported — cylindrical droplets reduce to
+the 2D coupled fit by translational symmetry.
+"""
+
+from wetting_angle_kit.analysis import (
+    CoupledFit3DAnalyzer,
+    DensityEstimator,
+)
+from wetting_angle_kit.analysis.temporal import TemporalAggregator
+from wetting_angle_kit.parsers import LammpsDumpParser, LammpsDumpWaterFinder
+
+# --- Step 1: Define the trajectory file ---
+filename = "../../tests/trajectories/traj_spherical_drop_4k.lammpstrj"
+
+# --- Step 2: Identify water-oxygen atoms ---
+wat_find = LammpsDumpWaterFinder(
+    filename,
+    oxygen_type=1,
+    hydrogen_type=2,
+)
+oxygen_indices = wat_find.get_water_oxygen_indices(frame_index=0)
+print("Number of water molecules:", len(oxygen_indices))
+
+# --- Step 3: Define the 3D grid ---
+# xi/yi are in the droplet-centred frame; zi is in the lab frame so
+# the wall position retains physical meaning.
+grid_params = {
+    "xi_0": -30.0,
+    "xi_f": 30.0,
+    "dx": 3.2,
+    "yi_0": -30.0,
+    "yi_f": 30.0,
+    "dy": 3.2,
+    "zi_0": 0.0,
+    "zi_f": 60.0,
+    "dz": 4.0,
+}
+
+# --- Step 4: Pick a density estimator ---
+# Top-hat histogram on the 3D sampling grid (default):
+estimator = DensityEstimator.binning()
+# Swap in the Gaussian KDE for smoother per-cell density:
+# estimator = DensityEstimator.gaussian(density_sigma=3.0)
+
+# --- Step 5: Build the analyzer ---
+analyzer = CoupledFit3DAnalyzer(
+    parser=LammpsDumpParser(filename),
+    atom_indices=oxygen_indices,
+    droplet_geometry="spherical",
+    grid_params=grid_params,
+    density_estimator=estimator,
+    # Pool all frames into a single batch — the 3D density needs more
+    # atoms than the 2D one for comparable per-cell noise.
+    temporal_aggregator=TemporalAggregator(batch_size=-1),
+)
+
+# --- Step 6: Run analysis ---
+n_frames = LammpsDumpParser(filename).frame_count()
+results = analyzer.analyze(range(0, n_frames))
+print("Mean contact angle (°):", results.mean_angle)
+
+# Per-batch detail:
+batch = results.batches[0]
+print(
+    f"Frames {batch.frames[0]}–{batch.frames[-1]}: "
+    f"angle = {batch.angle:.2f}°, "
+    f"R_eq = {batch.model_params['R_eq']:.2f} Å, "
+    f"z_wall = {batch.model_params['zi_0']:.2f} Å"
+)
+print(
+    f"Droplet centre: "
+    f"xi_c = {batch.model_params['xi_c']:.2f} Å, "
+    f"yi_c = {batch.model_params['yi_c']:.2f} Å, "
+    f"zi_c = {batch.model_params['zi_c']:.2f} Å"
+)

docs/examples/coupled_fit_ca.py

@@ -0,0 +1,81 @@
+"""Coupled-fit contact-angle example.
+
+Runs the coupled hyperbolic-tangent fit on a 2D density grid via
+:class:`CoupledFit2DAnalyzer`. The analyzer solves interface extraction,
+wall detection, and surface fitting together — one robust angle per pooled
+batch.
+
+Two density estimators are shown:
+
+- :meth:`DensityEstimator.binning` (the default) — top-hat histogram
+  with geometry-aware ``dV`` normalisation. Fast and exact; intrinsically
+  noisy at low per-cell counts.
+- :meth:`DensityEstimator.gaussian` — 3D Gaussian KDE on the cell
+  centres. Smooth density field with no per-cell Poisson noise; the
+  estimator of choice when running per-frame analyses or on systems
+  with low atom density per cell.
+"""
+
+from wetting_angle_kit.analysis import (
+    CoupledFit2DAnalyzer,
+    DensityEstimator,
+)
+from wetting_angle_kit.analysis.temporal import TemporalAggregator
+from wetting_angle_kit.parsers import LammpsDumpParser, LammpsDumpWaterFinder
+
+# --- Step 1: Define the trajectory file ---
+filename = "../../tests/trajectories/traj_spherical_drop_4k.lammpstrj"
+
+# --- Step 2: Identify water-oxygen atoms ---
+wat_find = LammpsDumpWaterFinder(
+    filename,
+    oxygen_type=1,
+    hydrogen_type=2,
+)
+oxygen_indices = wat_find.get_water_oxygen_indices(frame_index=0)
+print("Number of water molecules:", len(oxygen_indices))
+
+# --- Step 3: Define the grid ---
+grid_params = {
+    "xi_0": 0.0,
+    "xi_f": 70.0,
+    "dx": 2.0,
+    "zi_0": 0.0,
+    "zi_f": 70.0,
+    "dz": 2.0,
+}
+
+# --- Step 4: Pick a density estimator ---
+# Top-hat histogram on the sampling grid (default):
+estimator = DensityEstimator.binning()
+# Swap in the Gaussian KDE for smoother per-cell density. ``density_sigma``
+# is the Gaussian kernel width; 3 Å is a sensible default for
+# room-temperature water:
+# estimator = DensityEstimator.gaussian(density_sigma=2.5)
+
+# --- Step 5: Build the analyzer ---
+analyzer = CoupledFit2DAnalyzer(
+    parser=LammpsDumpParser(filename),
+    atom_indices=oxygen_indices,
+    droplet_geometry="spherical",
+    grid_params=grid_params,
+    density_estimator=estimator,
+    # Pool 10 frames per batch — the coupled fit benefits from
+    # statistics; ``batch_size=-1`` pools the entire trajectory.
+    temporal_aggregator=TemporalAggregator(batch_size=10),
+)
+
+# --- Step 6: Run analysis on a frame range ---
+# 20 frames at batch_size=10 gives two pooled batches.
+results = analyzer.analyze(range(0, 20))
+print("Mean contact angle (°):", results.mean_angle)
+print("Std across batches (°):", results.std_angle)
+
+# Per-batch detail:
+batch = results.batches[0]
+print(
+    f"Frames {batch.frames[0]}–{batch.frames[-1]}: "
+    f"angle = {batch.angle:.2f}°, "
+    f"R_eq = {batch.model_params['R_eq']:.2f} Å, "
+    f"z_wall = {batch.model_params['zi_0']:.2f} Å"
+)

docs/examples/parsing_trajectory_files.py

@@ -25,7 +25,7 @@
 )
 
 # --- Identify water oxygen indices for the first frame ---
-oxygen_indices = wat_find.get_water_oxygen_ids(frame_index=0)
+oxygen_indices = wat_find.get_water_oxygen_indices(frame_index=0)
 print(f"Number of water molecules: {len(oxygen_indices)}")
 
 # --- Initialize parser ---
@@ -90,5 +90,5 @@
 print("Total atoms loaded:", len(positions))
 
 # --- Extract subset of atoms (first 50) ---
-subset = xyz_parser.parse(frame_index=0, indices=list(range(50)))
+subset = xyz_parser.parse(frame_index=0, indices=list(range(24)))
 print("Subset (50 atoms) shape:", subset.shape)