Error Analysis

import numpy as np
from typing import Callable, Tuple


class ErrorAnalyzer:
    """Comprehensive error analysis toolkit."""

    @staticmethod
    def condition_number_matrix(A: np.ndarray) -> float:
        """Compute 2-norm condition number of a matrix."""
        return np.linalg.cond(A, ord=2)

    @staticmethod
    def relative_error(computed: float, exact: float) -> float:
        """Compute relative error."""
        if exact == 0:
            return abs(computed)
        return abs(computed - exact) / abs(exact)

    @staticmethod
    def finite_difference_error(f: Callable, df: Callable, x: float,
                                h_range: Tuple[float, float] = (1e-16, 1)):
        """Analyze truncation vs roundoff error in finite differences."""
        h_values = np.logspace(np.log10(h_range[0]), np.log10(h_range[1]), 200)
        errors = []

        exact = df(x)
        for h in h_values:
            # Central difference
            approx = (f(x + h) - f(x - h)) / (2 * h)
            error = abs(approx - exact)
            errors.append(error)

        errors = np.array(errors)
        optimal_h = h_values[np.argmin(errors)]

        return h_values, errors, optimal_h

    @staticmethod
    def demonstrate_cancellation():
        """Demonstrate catastrophic cancellation."""
        print("\n--- Catastrophic Cancellation Demo ---")
        x = 1.0
        h_values = np.logspace(-1, -16, 16)

        print(f"Computing derivative of f(x) = sin(x) at x = 1")
        exact = np.cos(1.0)

        for h in h_values:
            # Forward difference (suffers from cancellation)
            fwd = (np.sin(1 + h) - np.sin(1)) / h
            # Central difference (better but still has issues)
            cen = (np.sin(1 + h) - np.sin(1 - h)) / (2 * h)

            fwd_err = abs(fwd - exact)
            cen_err = abs(cen - exact)

            print(f"  h={h:.0e}: Forward err={fwd_err:.2e}, "
                  f"Central err={cen_err:.2e}")

    @staticmethod
    def welford_online(data: np.ndarray) -> Tuple[float, float]:
        """Welford's online algorithm for mean and variance."""
        n = 0
        mean = 0.0
        M2 = 0.0

        for x in data:
            n += 1
            delta = x - mean
            mean += delta / n
            delta2 = x - mean
            M2 += delta * delta2

        variance = M2 / n if n > 1 else 0.0
        return mean, variance

    @staticmethod
    def naive_mean_variance(data: np.ndarray) -> Tuple[float, float]:
        """Naive two-pass mean and variance (unstable)."""
        n = len(data)
        mean = np.mean(data)
        variance = np.mean((data - mean) ** 2)
        return mean, variance

    @staticmethod
    def demonstrate_stability():
        """Compare stable vs unstable algorithms."""
        print("\n--- Algorithm Stability Demo ---")

        # Generate data with large mean and small variance
        np.random.seed(42)
        data = np.random.uniform(1e7, 1e7 + 1, 1_000_000)
        true_mean = np.mean(data)
        true_var = 1/12  # Variance of uniform(0,1)

        # Naive approach
        naive_mean, naive_var = ErrorAnalyzer.naive_mean_variance(data)

        # Welford
        welford_mean, welford_var = ErrorAnalyzer.welford_online(data)

        print(f"True mean:       {true_mean:.10f}")
        print(f"Naive mean:      {naive_mean:.10f}")
        print(f"Welford mean:    {welford_mean:.10f}")
        print(f"Mean error (naive):   {abs(naive_mean - true_mean):.2e}")
        print(f"Mean error (Welford): {abs(welford_mean - true_mean):.2e}")

        # Variance comparison
        naive_var_pass1 = np.mean(data**2) - np.mean(data)**2  # One-pass (unstable)
        print(f"\nNaive one-pass var: {naive_var_pass1:.6f}")
        print(f"Naive two-pass var: {naive_var:.6f}")
        print(f"Welford var:        {welford_var:.6f}")
        print(f"True var:           {true_var:.6f}")

    @staticmethod
    def demonstrate_condition_number():
        """Show effect of condition number on linear system solving."""
        print("\n--- Condition Number Effect on Linear Systems ---")

        # Well-conditioned
        A1 = np.array([[1, 0], [0, 2]])
        b1 = np.array([1, 1])
        x1 = np.linalg.solve(A1, b1)

        # Ill-conditioned
        A2 = np.array([[1, 1], [1, 1.0001]])
        b2 = np.array([2, 2.0001])
        x2 = np.linalg.solve(A2, b2)

        print(f"Well-conditioned (κ={np.linalg.cond(A1):.1f}):")
        print(f"  Solution: {x1}")

        print(f"Ill-conditioned (κ={np.linalg.cond(A2):.0f}):")
        print(f"  Solution: {x2}")

        # Perturbation experiment
        print(f"\nPerturbation experiment:")
        delta_b = np.array([1e-10, -1e-10])

        x1_pert = np.linalg.solve(A1, b1 + delta_b)
        x2_pert = np.linalg.solve(A2, b2 + delta_b)

        print(f"  Well-conditioned: Δx = {np.linalg.norm(x1_pert - x1):.2e}")
        print(f"  Ill-conditioned:  Δx = {np.linalg.norm(x2_pert - x2):.2e}")

    @staticmethod
    def demonstrate_backward_error():
        """Show backward error interpretation."""
        print("\n--- Backward Error Demonstration ---")

        # Solve Ax = b with known solution
        A = np.array([[1, 2], [3, 4]], dtype=float)
        x_true = np.array([1, 1])
        b = A @ x_true

        # Solve with perturbed matrix
        delta_A = np.array([[1e-10, 0], [0, 1e-10]])
        A_pert = A + delta_A
        x_computed = np.linalg.solve(A_pert, b)

        # Backward error: what perturbation to b would give this answer?
        delta_b = b - A @ x_computed
        backward_error = np.linalg.norm(delta_b) / np.linalg.norm(b)

        print(f"True solution:     {x_true}")
        print(f"Computed solution: {x_computed}")
        print(f"Forward error:     {np.linalg.norm(x_computed - x_true):.2e}")
        print(f"Backward error:    {backward_error:.2e}")
        print(f"Condition number:  {np.linalg.cond(A):.1f}")


def full_error_analysis():
    """Run comprehensive error analysis demonstrations."""
    print("=" * 60)
    print("NUMERICAL ERROR ANALYSIS")
    print("=" * 60)

    analyzer = ErrorAnalyzer()

    # 1. Condition numbers
    print("\n--- Condition Number Examples ---")
    matrices = {
        "Identity": np.eye(3),
        "Hilbert (5×5)": np.linalg.inv(np.vander(np.arange(5), increasing=True)),
        "Random 3×3": np.random.randn(3, 3),
        "Near-singular": np.array([[1, 1], [1, 1.0001]]),
    }

    for name, A in matrices.items():
        kappa = analyzer.condition_number_matrix(A)
        print(f"  {name:<20}: κ = {kappa:.4e}")

    # 2. Finite difference error
    print("\n--- Finite Difference Error Analysis ---")
    f = np.sin
    df = np.cos
    x = 1.0

    h_vals, errors, optimal_h = analyzer.finite_difference_error(f, df, x)
    print(f"  Optimal step size h: {optimal_h:.2e}")
    print(f"  Minimum error:       {np.min(errors):.2e}")
    print(f"  Theoretical optimal: {np.sqrt(np.finfo(float).eps):.2e}")

    # 3. Cancellation
    analyzer.demonstrate_cancellation()

    # 4. Stability comparison
    analyzer.demonstrate_stability()

    # 5. Condition number effect
    analyzer.demonstrate_condition_number()

    # 6. Backward error
    analyzer.demonstrate_backward_error()


full_error_analysis()