diff --git a/scripts/run_baseline.py b/scripts/run_baseline.py index e2e71cd..1982696 100644 --- a/scripts/run_baseline.py +++ b/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) diff --git a/scripts/run_sample_count_sweep.py b/scripts/run_sample_count_sweep.py new file mode 100644 index 0000000..26d5d56 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/src/sparse_recon/methods/cubic_spline.py b/src/sparse_recon/methods/cubic_spline.py new file mode 100644 index 0000000..e6ec796 --- /dev/null +++ b/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} diff --git a/src/sparse_recon/sampling/geometries.py b/src/sparse_recon/sampling/geometries.py index 2b8bc97..efd4a61 100644 --- a/src/sparse_recon/sampling/geometries.py +++ b/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)}") diff --git a/tests/test_geometries.py b/tests/test_geometries.py new file mode 100644 index 0000000..672d566 --- /dev/null +++ b/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)