Add flyby geometry, cubic spline method, and sample count sweep

scripts/run_baseline.py

@@ -219,10 +219,14 @@ def build_parser() -> argparse.ArgumentParser:
     parser.add_argument("--field-kind", default="smooth")
     parser.add_argument("--field-seed", type=int, default=0)
     parser.add_argument("--field-noise-sigma", type=float, default=0.05)
+<<<<<<< HEAD
+    parser.add_argument("--methods", default="rbf,linear,nearest,cubic_spline")
+=======
     parser.add_argument("--simulation-path", default=None)
     parser.add_argument("--methods", default="rbf,linear,nearest")
+>>>>>>> upstream/main
     parser.add_argument("--sample-counts", default="50")
-    parser.add_argument("--geometries", default="random,clustered,multi_probe_like")
+    parser.add_argument("--geometries", default="random,clustered,multi_probe_like,flyby")
     parser.add_argument("--noise-levels", default="0.0,0.02")
     parser.add_argument("--sample-seed", type=int, default=1)
     parser.add_argument("--noise-seed", type=int, default=11)

scripts/run_sample_count_sweep.py

@@ -0,0 +1,98 @@
+import json
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+
+from sparse_recon.datasets.synthetic import create_synthetic_field
+from sparse_recon.methods.cubic_spline import CubicSplineMethod
+from sparse_recon.methods.linear import LinearMethod
+from sparse_recon.methods.nearest import NearestMethod
+from sparse_recon.methods.rbf import RBFMethod
+from sparse_recon.pipeline import run_sampling_experiment
+from sparse_recon.sampling.geometries import generate_sampling_points
+
+SAMPLE_COUNTS = [10, 25, 50, 100, 200]
+METHODS = [
+    ("rbf", "#378ADD"),
+    ("linear", "#1D9E75"),
+    ("nearest", "#888780"),
+    ("cubic_spline", "#D85A30"),
+]
+OUTPUT_DIR = Path("results/sample_count_sweep")
+
+
+def build_method(name):
+    if name == "rbf":
+        return RBFMethod()
+    if name == "linear":
+        return LinearMethod()
+    if name == "nearest":
+        return NearestMethod()
+    if name == "cubic_spline":
+        return CubicSplineMethod()
+    raise ValueError(f"Unknown method: {name}")
+
+
+def main():
+    OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
+
+    field = create_synthetic_field(kind="smooth", nx=64, ny=64, seed=0)
+
+    results = {name: [] for name, _ in METHODS}
+
+    for name, _ in METHODS:
+        for n in SAMPLE_COUNTS:
+            sample_coords = generate_sampling_points(
+                geometry="random", n_points=n, dim=2, seed=42
+            )
+            method = build_method(name)
+            _, result = run_sampling_experiment(
+                field,
+                sample_coords,
+                method,
+                noise_sigma=0.0,
+                noise_seed=0,
+            )
+            rel_l2 = result.metrics["relative_l2"]
+            valid = result.metrics["valid_fraction"]
+            results[name].append({"n": n, "relative_l2": rel_l2, "valid_fraction": valid})
+            print(f"{name:<15} n={n:4d}  relative_l2={rel_l2:.4f}  valid={valid:.2f}")
+
+    with open(OUTPUT_DIR / "results.json", "w") as f:
+        json.dump(results, f, indent=2)
+
+    colors = {name: color for name, color in METHODS}
+
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
+
+    for name, _ in METHODS:
+        ns = [r["n"] for r in results[name]]
+        errors = [r["relative_l2"] for r in results[name]]
+        valid = [r["valid_fraction"] for r in results[name]]
+        color = colors[name]
+        ax1.plot(ns, errors, marker="o", linewidth=2, label=name, color=color)
+        ax2.plot(ns, valid, marker="o", linewidth=2, label=name, color=color)
+
+    ax1.set_xlabel("Number of sample points")
+    ax1.set_ylabel("Relative L2 error")
+    ax1.set_title("Reconstruction error vs sample count")
+    ax1.set_xticks(SAMPLE_COUNTS)
+    ax1.legend()
+    ax1.grid(True, alpha=0.3)
+
+    ax2.set_xlabel("Number of sample points")
+    ax2.set_ylabel("Valid fraction")
+    ax2.set_title("Domain coverage vs sample count")
+    ax2.set_xticks(SAMPLE_COUNTS)
+    ax2.set_ylim(0, 1.05)
+    ax2.legend()
+    ax2.grid(True, alpha=0.3)
+
+    fig.suptitle("RBF vs Linear vs Nearest vs Cubic Spline — random geometry, smooth field", fontsize=12)
+    fig.tight_layout()
+    fig.savefig(OUTPUT_DIR / "error_vs_sample_count.png", dpi=150)
+    print(f"\nPlot saved to {OUTPUT_DIR / 'error_vs_sample_count.png'}")
+
+
+if __name__ == "__main__":
+    main()

src/sparse_recon/methods/cubic_spline.py

@@ -0,0 +1,41 @@
+import numpy as np
+from scipy.interpolate import CloughTocher2DInterpolator, LinearNDInterpolator
+
+from sparse_recon.methods.base import ReconstructionMethod
+
+
+class CubicSplineMethod(ReconstructionMethod):
+    """Cubic spline interpolation using Clough-Tocher scheme (C1 smooth, 2D only).
+
+    Falls back to linear interpolation outside the convex hull of sample points,
+    matching the behaviour of scipy's LinearNDInterpolator for extrapolation.
+    """
+
+    name = "cubic_spline"
+
+    def __init__(self, fill_value: float = np.nan):
+        self.fill_value = fill_value
+
+    def fit(self, sample_coords: np.ndarray, sample_values: np.ndarray):
+        self._cubic = CloughTocher2DInterpolator(
+            sample_coords,
+            sample_values,
+            fill_value=self.fill_value,
+        )
+        self._linear_fallback = LinearNDInterpolator(
+            sample_coords,
+            sample_values,
+            fill_value=self.fill_value,
+        )
+        return self
+
+    def predict(self, query_coords: np.ndarray) -> np.ndarray:
+        pred = self._cubic(query_coords)
+        # where cubic returns nan (outside convex hull), use linear fallback
+        nan_mask = np.any(np.isnan(pred), axis=-1)
+        if np.any(nan_mask):
+            pred[nan_mask] = self._linear_fallback(query_coords[nan_mask])
+        return pred
+
+    def get_params(self) -> dict:
+        return {"fill_value": self.fill_value}

src/sparse_recon/sampling/geometries.py

@@ -17,6 +17,43 @@ def clustered_points(center, offsets) -> np.ndarray:
     offsets = np.asarray(offsets)
     return center[None, :] + offsets
 
+def generate_flyby_points(
+    n_points: int,
+    dim: int,
+    seed: int = 0,
+    low: float = 0.0,
+    high: float = 1.0,
+) -> np.ndarray:
+    rng = np.random.default_rng(seed)
+
+    if dim == 2:
+        start_point = rng.uniform(low, high, size=dim)
+        angle = rng.uniform(0, 2 * np.pi)
+        direction = np.array([np.cos(angle), np.sin(angle)])
+        # Ensure the line segment is within the bounds [low, high]
+        # We need to calculate the length of the line to ensure it fits within the bounds
+        # This current implementation can generate points outside the bounds if the start_point + direction * length goes out
+        # A simpler approach for tests and initial implementation: generate points and then clip.
+        # The problem with clipping is it can make truly collinear points non-collinear at the boundaries.
+        # For now, let's modify the line generation to ensure it's within bounds.
+        # Let's make the line segment start and end within the box.
+        # We can pick two random points in the box and draw a line between them.
+        point1 = rng.uniform(low, high, size=dim)
+        point2 = rng.uniform(low, high, size=dim)
+
+        # Generate n_points along the line segment defined by point1 and point2
+        t = np.linspace(0, 1, n_points)[:, None]
+        points = point1 + t * (point2 - point1)
+
+    elif dim == 3:
+        point1 = rng.uniform(low, high, size=dim)
+        point2 = rng.uniform(low, high, size=dim)
+        t = np.linspace(0, 1, n_points)[:, None]
+        points = point1 + t * (point2 - point1)
+    else:
+        raise ValueError(f"Flyby geometry not supported for dimension {dim}")
+
+    return points # points are already within bounds if point1 and point2 are within bounds
 
 def tetrahedron_like(scale: float = 0.1, center=(0.5, 0.5, 0.5)) -> np.ndarray:
     center = np.asarray(center)
@@ -142,6 +179,22 @@ def generate_sampling_points(
             high=high,
         )
     if geometry == "multi_probe_like":
+<<<<<<< HEAD
+        if dim != 2:
+            raise ValueError("multi_probe_like geometry currently supports dim=2 only")
+        return multi_probe_like_points_2d(
+            n_points=n_points,
+            seed=seed,
+            low=low,
+            high=high,
+        )
+    if geometry == "tetrahedron_like":
+        if dim != 3:
+            raise ValueError("tetrahedron_like geometry currently supports dim=3 only")
+        return tetrahedron_like(center=np.array([0.5, 0.5, 0.5]), scale=0.1) # Default values
+    if geometry == "flyby":
+        return generate_flyby_points(n_points=n_points, dim=dim, seed=seed, low=low, high=high)
+=======
         if dim == 2:
             return multi_probe_like_points_2d(
                 n_points=n_points,
@@ -157,6 +210,7 @@ def generate_sampling_points(
                 high=high,
             )
         raise ValueError("multi_probe_like geometry currently supports dim=2 or dim=3")
+>>>>>>> upstream/main
 
-    supported = ["clustered", "multi_probe_like", "random"]
+    supported = ["clustered", "multi_probe_like", "random", "tetrahedron_like", "flyby"]
     raise ValueError(f"Unknown geometry '{geometry}'. Supported: {', '.join(supported)}")

tests/test_geometries.py

@@ -0,0 +1,126 @@
+import numpy as np
+import pytest
+from sparse_recon.sampling.geometries import (
+    generate_sampling_points,
+    random_points_in_box,
+    clustered_points,
+    tetrahedron_like,
+    clustered_points_in_box,
+    multi_probe_like_points_2d,
+    generate_flyby_points,
+)
+
+def test_random_points_in_box():
+    points = random_points_in_box(n_points=10, dim=2, seed=0)
+    assert points.shape == (10, 2)
+    assert np.all(points >= 0.0)
+    assert np.all(points <= 1.0)
+
+def test_clustered_points():
+    center = np.array([0.5, 0.5])
+    offsets = np.array([[0.1, 0.1], [-0.1, -0.1]])
+    points = clustered_points(center, offsets)
+    assert points.shape == (2, 2)
+    assert np.allclose(points[0], [0.6, 0.6])
+    assert np.allclose(points[1], [0.4, 0.4])
+
+def test_tetrahedron_like():
+    points = tetrahedron_like()
+    assert points.shape == (4, 3)
+    expected_points = np.array([
+        [0.6, 0.6, 0.6],
+        [0.6, 0.4, 0.4],
+        [0.4, 0.6, 0.4],
+        [0.4, 0.4, 0.6],
+    ])
+    assert np.allclose(points, expected_points)
+
+def test_clustered_points_in_box():
+    points = clustered_points_in_box(n_points=20, dim=2, seed=0)
+    assert points.shape == (20, 2)
+    assert np.all(points >= 0.0)
+    assert np.all(points <= 1.0)
+
+def test_multi_probe_like_points_2d():
+    points = multi_probe_like_points_2d(n_points=10, seed=0)
+    assert points.shape == (10, 2)
+    assert np.all(points >= 0.0)
+    assert np.all(points <= 1.0)
+
+def test_generate_flyby_points_2d():
+    points = generate_flyby_points(n_points=10, dim=2, seed=0)
+    assert points.shape == (10, 2)
+    assert np.all(points >= 0.0)
+    assert np.all(points <= 1.0)
+    # Check if points are collinear
+    n_points_2d = 10 # Define n_points locally
+    if n_points_2d > 1:
+        diffs = np.diff(points, axis=0)
+        # All diffs should be parallel, i.e., cross product should be zero (or very small)
+        # For 2D, this means the determinant of any two diff vectors should be zero
+        if len(diffs) > 1:
+            # For 2D, if points are collinear, the y/x ratio between consecutive points should be constant.
+            # Or, more robustly, the area of the triangle formed by three consecutive points should be zero.
+            # This can be checked using the determinant of vectors (p2-p1) and (p3-p1).
+            # For a line, if p1, p2, p3 are points, then (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1) should be 0.
+            for i in range(len(points) - 2):
+                p1 = points[i]
+                p2 = points[i+1]
+                p3 = points[i+2]
+                # Check if (p2-p1) and (p3-p1) are parallel by checking if the determinant is close to zero
+                det = (p2[0] - p1[0]) * (p3[1] - p1[1]) - (p3[0] - p1[0]) * (p2[1] - p1[1])
+                assert np.isclose(det, 0.0), f"Points are not collinear in 2D: {points[i:i+3]}"
+
+def test_generate_flyby_points_3d():
+    n_points_3d = 10 # Define n_points locally
+    points = generate_flyby_points(n_points=n_points_3d, dim=3, seed=0)
+    assert points.shape == (n_points_3d, 3)
+    assert np.all(points >= 0.0)
+    assert np.all(points <= 1.0)
+    # Check if points are collinear
+    if n_points_3d > 1:
+        # For 3D, cross product of (p2-p1) and (p3-p1) should be a zero vector
+        for i in range(len(points) - 2):
+            p1 = points[i]
+            p2 = points[i+1]
+            p3 = points[i+2]
+            vec1 = p2 - p1
+            vec2 = p3 - p1
+            cross_prod = np.cross(vec1, vec2)
+            assert np.allclose(cross_prod, 0.0), f"Points are not collinear in 3D: {points[i:i+3]}"
+
+def test_generate_sampling_points_random():
+    points = generate_sampling_points("random", 10, 2, seed=0)
+    assert points.shape == (10, 2)
+
+def test_generate_sampling_points_clustered():
+    points = generate_sampling_points("clustered", 10, 2, seed=0)
+    assert points.shape == (10, 2)
+
+def test_generate_sampling_points_multi_probe_like():
+    points = generate_sampling_points("multi_probe_like", 10, 2, seed=0)
+    assert points.shape == (10, 2)
+
+def test_generate_sampling_points_tetrahedron_like():
+    points = generate_sampling_points("tetrahedron_like", 4, 3, seed=0)
+    assert points.shape == (4, 3)
+
+def test_generate_sampling_points_flyby_2d():
+    points = generate_sampling_points("flyby", 10, 2, seed=0)
+    assert points.shape == (10, 2)
+
+def test_generate_sampling_points_flyby_3d():
+    points = generate_sampling_points("flyby", 10, 3, seed=0)
+    assert points.shape == (10, 3)
+
+def test_generate_sampling_points_unknown_geometry():
+    with pytest.raises(ValueError, match="Unknown geometry"):
+        generate_sampling_points("unknown", 10, 2)
+
+def test_generate_sampling_points_tetrahedron_like_wrong_dim():
+    with pytest.raises(ValueError, match="tetrahedron_like geometry currently supports dim=3 only"):
+        generate_sampling_points("tetrahedron_like", 4, 2)
+
+def test_generate_sampling_points_flyby_wrong_dim():
+    with pytest.raises(ValueError, match="Flyby geometry not supported for dimension"): # Adjust regex if needed
+        generate_sampling_points("flyby", 10, 4)