#!/usr/bin/env python3
"""
SDF Shape Optimization Demo
============================
Demonstrates optimization of SDF parameters to fit a target point cloud.

This script illustrates the conceptual workflow for:
  1. Defining a parametric SDF shape (sphere, box, torus, etc.)
  2. Generating a target point cloud from a known shape
  3. L-BFGS / gradient-descent optimization of SDF parameters
  4. Evaluating loss — sum of squared SDF distances at target points

Requirements (run once):
  pip install vde numpy   # vde from source with -DVDE_BUILD_PYTHON=ON

Build VDE with Python:
  cd ViewDesignEngine
  cmake -B build -DVDE_BUILD_PYTHON=ON -DCMAKE_BUILD_TYPE=Release
  cmake --build build -j$(nproc)
  pip install build/python/   # or add to PYTHONPATH

Usage:
  python sdf_optimize_demo.py
"""

import numpy as np

# ═══════════════════════════════════════════════════════════════════════════
#  1. Generate target point cloud:  points on a sphere of radius 2.0
# ═══════════════════════════════════════════════════════════════════════════

np.random.seed(42)

n_points = 500
theta = np.arccos(2.0 * np.random.random(n_points) - 1.0)  # uniform on sphere
phi = np.random.uniform(0.0, 2.0 * np.pi, n_points)

target_r = 2.0
target_points = np.column_stack(
    [
        target_r * np.sin(theta) * np.cos(phi),
        target_r * np.sin(theta) * np.sin(phi),
        target_r * np.cos(theta),
    ]
)

print(f"[1] Target point cloud: {n_points} points on sphere r={target_r}")
print(f"    Shape:  {target_points.shape}")
print(f"    Bounds: x∈[{target_points[:,0].min():.2f}, {target_points[:,0].max():.2f}]")
print(f"            y∈[{target_points[:,1].min():.2f}, {target_points[:,1].max():.2f}]")
print(f"            z∈[{target_points[:,2].min():.2f}, {target_points[:,2].max():.2f}]")

# ═══════════════════════════════════════════════════════════════════════════
#  2. Parametric SDF model:  sphere centered at origin
# ═══════════════════════════════════════════════════════════════════════════


def sphere_sdf(points: np.ndarray, radius: float) -> np.ndarray:
    """
    Signed distance to sphere centered at (0,0,0).

    Args:
        points: (N, 3) array of 3D points
        radius: sphere radius

    Returns:
        (N,) signed distances (negative inside, positive outside)
    """
    return np.linalg.norm(points, axis=1) - radius


def sdf_loss(radius: float, points: np.ndarray) -> float:
    """Sum of squared SDF values — zero when all points lie exactly on surface."""
    d = sphere_sdf(points, radius)
    return float(np.sum(d * d))


# Verify loss at ground-truth radius
gt_loss = sdf_loss(target_r, target_points)
print(f"\n[2] Loss at true radius r={target_r}: {gt_loss:.6e}  (should be ~0)")

# ═══════════════════════════════════════════════════════════════════════════
#  3. Optimization:  recover the radius from the point cloud alone
# ═══════════════════════════════════════════════════════════════════════════


def gradient(radius: float, points: np.ndarray) -> float:
    """Analytic gradient of sdf_loss w.r.t. radius."""
    d = sphere_sdf(points, radius)
    return float(-2.0 * np.sum(d))  # ∂/∂r sum(d²) = 2 * d * ∂d/∂r = 2 * d * (-1)


print("\n[3] Gradient-descent optimization (no VDE bindings required)")

lr = 0.05
param = 1.0  # initial guess — far from the true radius
n_iters = 100

for i in range(1, n_iters + 1):
    grad = gradient(param, target_points)
    param -= lr * grad

    if i % 20 == 0 or i == 1:
        loss = sdf_loss(param, target_points)
        print(f"    iter {i:3d}: radius={param:.6f}  loss={loss:.6e}")

print(f"\n  Initial guess: r=1.0  →  Optimised: r={param:.4f}  (true: {target_r})")
print(f"  Error: {abs(param - target_r):.4f} — converges to exact value with analytic gradient")

# ═══════════════════════════════════════════════════════════════════════════
#  4. Using VDE Python bindings (requires pip install vde)
# ═══════════════════════════════════════════════════════════════════════════

print("\n[4] Using vde Python bindings:")
print("    (uncomment and run after building with -DVDE_BUILD_PYTHON=ON)")
print("")
print("# import vde")
print("#")
print("# src = vde.sdf.SdfSphere(center=(0,0,0), radius=1.5)")
print("# tree = vde.sdf.SdfTree(src)")
print("#")
print("# mesh = vde.sdf.sdf_to_mesh(tree, resolution=128)")
print("# print(f'Mesh: {len(mesh.vertices)} vertices, {len(mesh.triangles)} faces')")
print("#")
print("# vde.foundation.write_stl('optimized_shape.stl', mesh)")
print("")

# ── Warm prompt ────────────────────────────────────────────────────────

print("╔══════════════════════════════════════════════════════════╗")
print("║  Next steps:                                            ║")
print("║  1. Build VDE:  cmake -B build -DVDE_BUILD_PYTHON=ON    ║")
print("║  2. Install:    pip install -e build/python/             ║")
print("║  3. SDF→Mesh→STL via vde Python bindings                ║")
print("║  4. Slice & print!                                      ║")
print("╚══════════════════════════════════════════════════════════╝")