From fd3c4c5b3d6e0746c75313df45dde50f49b7b45a Mon Sep 17 00:00:00 2001 From: Mark Grayson <118820911+sahildando@users.noreply.github.com> Date: Tue, 7 Apr 2026 22:07:25 +0530 Subject: [PATCH 1/2] Add adv_math_engine module with calculus, algebra, series, CLI, and tests --- adv_math_engine/README.md | 57 ++++++++++ adv_math_engine/__init__.py | 41 +++++++ adv_math_engine/benchmarks/benchmark_diff.py | 51 +++++++++ adv_math_engine/benchmarks/results.md | 17 +++ adv_math_engine/examples/example_usage.py | 18 +++ adv_math_engine/series_expansion.py | 102 +++++++++++++++++ adv_math_engine/tests/__init__.py | 1 + adv_math_engine/tests/coverage_report.txt | 7 ++ .../tests/test_series_expansion.py | 41 +++++++ adv_math_engine/tests/test_vector_algebra.py | 46 ++++++++ adv_math_engine/tests/test_vector_calculus.py | 56 +++++++++ adv_math_engine/utils.py | 27 +++++ adv_math_engine/vector_algebra.py | 88 ++++++++++++++ adv_math_engine/vector_calculus.py | 107 ++++++++++++++++++ .../visualizations/plot_gradient_field.py | 26 +++++ .../plot_series_approximation.py | 28 +++++ main.py | 50 ++++++++ 17 files changed, 763 insertions(+) create mode 100644 adv_math_engine/README.md create mode 100644 adv_math_engine/__init__.py create mode 100644 adv_math_engine/benchmarks/benchmark_diff.py create mode 100644 adv_math_engine/benchmarks/results.md create mode 100644 adv_math_engine/examples/example_usage.py create mode 100644 adv_math_engine/series_expansion.py create mode 100644 adv_math_engine/tests/__init__.py create mode 100644 adv_math_engine/tests/coverage_report.txt create mode 100644 adv_math_engine/tests/test_series_expansion.py create mode 100644 adv_math_engine/tests/test_vector_algebra.py create mode 100644 adv_math_engine/tests/test_vector_calculus.py create mode 100644 adv_math_engine/utils.py create mode 100644 adv_math_engine/vector_algebra.py create mode 100644 adv_math_engine/vector_calculus.py create mode 100644 adv_math_engine/visualizations/plot_gradient_field.py create mode 100644 adv_math_engine/visualizations/plot_series_approximation.py create mode 100644 main.py diff --git a/adv_math_engine/README.md b/adv_math_engine/README.md new file mode 100644 index 000000000000..c5abb3e9fd3b --- /dev/null +++ b/adv_math_engine/README.md @@ -0,0 +1,57 @@ +# adv_math_engine + +Production-grade module for: +- Vector algebra +- Vector calculus +- Taylor / Maclaurin series expansion + +## Structure + +```text +adv_math_engine/ + ├── vector_algebra.py + ├── vector_calculus.py + ├── series_expansion.py + ├── utils.py + ├── benchmarks/ + ├── visualizations/ + ├── examples/ + └── tests/ +``` + +## CLI + +```bash +python main.py --function "sin(x)" --series taylor --order 5 --x 0.5 --center 0.0 +``` + +## Visualizations + +Run: + +```bash +python adv_math_engine/visualizations/plot_gradient_field.py +python adv_math_engine/visualizations/plot_series_approximation.py +``` + +Outputs: +- `adv_math_engine/visualizations/gradient_field.png` +- `adv_math_engine/visualizations/series_approximation.png` + +## Benchmarks + +```bash +python adv_math_engine/benchmarks/benchmark_diff.py +``` + +See `adv_math_engine/benchmarks/results.md`. + +## Coverage + +Suggested command: + +```bash +pytest adv_math_engine/tests --cov=adv_math_engine --cov-report=term-missing +``` + +See `adv_math_engine/tests/coverage_report.txt`. diff --git a/adv_math_engine/__init__.py b/adv_math_engine/__init__.py new file mode 100644 index 000000000000..fefb18beeba6 --- /dev/null +++ b/adv_math_engine/__init__.py @@ -0,0 +1,41 @@ +"""Advanced mathematics engine for vector algebra, vector calculus, and series expansions.""" + +from adv_math_engine.series_expansion import ( + estimate_lagrange_remainder, + maclaurin_series, + taylor_series, +) +from adv_math_engine.vector_algebra import ( + angle_between, + check_linear_independence, + cross_product, + dot_product, + gram_schmidt, + vector_projection, +) +from adv_math_engine.vector_calculus import ( + curl, + divergence, + gradient, + line_integral, + partial_derivative, + surface_integral, +) + +__all__ = [ + "angle_between", + "check_linear_independence", + "cross_product", + "curl", + "divergence", + "dot_product", + "estimate_lagrange_remainder", + "gradient", + "gram_schmidt", + "line_integral", + "maclaurin_series", + "partial_derivative", + "surface_integral", + "taylor_series", + "vector_projection", +] diff --git a/adv_math_engine/benchmarks/benchmark_diff.py b/adv_math_engine/benchmarks/benchmark_diff.py new file mode 100644 index 000000000000..bdaa1ab0a67d --- /dev/null +++ b/adv_math_engine/benchmarks/benchmark_diff.py @@ -0,0 +1,51 @@ +from __future__ import annotations + +import numpy as np + +from adv_math_engine.series_expansion import _numerical_nth_derivative +from adv_math_engine.utils import timed_call + +try: + import sympy as sp +except Exception: # pragma: no cover + sp = None + + +def numerical_benchmark(iterations: int = 2000) -> float: + func = lambda x: np.sin(x) * np.exp(x) + + def run() -> float: + total = 0.0 + for value in np.linspace(-1.0, 1.0, iterations): + total += _numerical_nth_derivative(func, n=1, at=float(value)) + return total + + _, elapsed = timed_call(run) + return elapsed + + +def symbolic_benchmark(iterations: int = 2000) -> float | None: + if sp is None: + return None + x = sp.symbols("x") + derivative_expr = sp.diff(sp.sin(x) * sp.exp(x), x) + evaluator = sp.lambdify(x, derivative_expr, "numpy") + + def run() -> float: + total = 0.0 + for value in np.linspace(-1.0, 1.0, iterations): + total += float(evaluator(value)) + return total + + _, elapsed = timed_call(run) + return elapsed + + +if __name__ == "__main__": + num_time = numerical_benchmark() + sym_time = symbolic_benchmark() + print(f"numerical_time_seconds={num_time:.6f}") + if sym_time is None: + print("symbolic_time_seconds=unavailable") + else: + print(f"symbolic_time_seconds={sym_time:.6f}") diff --git a/adv_math_engine/benchmarks/results.md b/adv_math_engine/benchmarks/results.md new file mode 100644 index 000000000000..2d254bbb1c5c --- /dev/null +++ b/adv_math_engine/benchmarks/results.md @@ -0,0 +1,17 @@ +# Benchmark Results + +## Status +Benchmark execution is prepared via: + +```bash +python adv_math_engine/benchmarks/benchmark_diff.py +``` + +In this environment, benchmark execution could not be completed because required numerical dependencies (NumPy/SymPy) were unavailable to the active Python interpreter. + +## Expected Output Format + +```text +numerical_time_seconds= +symbolic_time_seconds= +``` diff --git a/adv_math_engine/examples/example_usage.py b/adv_math_engine/examples/example_usage.py new file mode 100644 index 000000000000..cbe07c5a7335 --- /dev/null +++ b/adv_math_engine/examples/example_usage.py @@ -0,0 +1,18 @@ +from __future__ import annotations + +import numpy as np + +from adv_math_engine.series_expansion import maclaurin_series +from adv_math_engine.vector_algebra import dot_product, gram_schmidt +from adv_math_engine.vector_calculus import gradient + + +if __name__ == "__main__": + print("Dot:", dot_product([1, 2, 3], [4, 5, 6])) + print("Orthonormal basis:\n", gram_schmidt([[1, 1, 0], [1, 0, 1], [0, 1, 1]])) + + f = lambda x, y, z: x**2 + y**2 + z**2 + print("Gradient at (1,2,3):", gradient(f, [1, 2, 3])) + + x = np.array([0.1, 0.2, 0.3]) + print("sin(x) Maclaurin order=7:", maclaurin_series(x, order=7, function_name="sin")) diff --git a/adv_math_engine/series_expansion.py b/adv_math_engine/series_expansion.py new file mode 100644 index 000000000000..17f47d93766e --- /dev/null +++ b/adv_math_engine/series_expansion.py @@ -0,0 +1,102 @@ +"""Taylor/Maclaurin series expansion utilities.""" + +from __future__ import annotations + +import math +from collections.abc import Callable +from typing import Literal + +import numpy as np + +HAS_SYMPY = True +try: + import sympy as sp +except Exception: # pragma: no cover - optional dependency guard + HAS_SYMPY = False + + +SupportedFunction = Literal["exp", "sin", "cos", "log1p"] + + +def _builtin_derivative_value(name: SupportedFunction, n: int, at: float) -> float: + if name == "exp": + return float(math.exp(at)) + if name == "sin": + return float(math.sin(at + n * math.pi / 2.0)) + if name == "cos": + return float(math.cos(at + n * math.pi / 2.0)) + if name == "log1p": + if n == 0: + return float(math.log1p(at)) + return float(((-1) ** (n - 1)) * math.factorial(n - 1) / (1.0 + at) ** n) + msg = f"Unsupported function: {name}" + raise ValueError(msg) + + +def _numerical_nth_derivative(func: Callable[[float], float], n: int, at: float, h: float = 1e-4) -> float: + if n == 0: + return float(func(at)) + if n == 1: + return float((func(at + h) - func(at - h)) / (2.0 * h)) + return float( + (_numerical_nth_derivative(func, n - 1, at + h, h) - _numerical_nth_derivative(func, n - 1, at - h, h)) / (2.0 * h) + ) + + +def taylor_series( + x: float | np.ndarray, + order: int, + center: float = 0.0, + *, + function_name: SupportedFunction | None = None, + function: Callable[[float], float] | None = None, +) -> np.ndarray: + """Evaluate Taylor polynomial of given order around center. + + f(x) ≈ Σ[n=0..order] f^(n)(a) / n! * (x-a)^n + """ + if function_name is None and function is None: + msg = "Provide either function_name or function." + raise ValueError(msg) + + x_arr = np.asarray(x, dtype=float) + approximation = np.zeros_like(x_arr, dtype=float) + + for n in range(order + 1): + if function_name is not None: + derivative_at_center = _builtin_derivative_value(function_name, n, center) + else: + assert function is not None + derivative_at_center = _numerical_nth_derivative(function, n, center) + + approximation = approximation + derivative_at_center / math.factorial(n) * (x_arr - center) ** n + + return approximation + + +def maclaurin_series( + x: float | np.ndarray, + order: int, + *, + function_name: SupportedFunction | None = None, + function: Callable[[float], float] | None = None, +) -> np.ndarray: + """Evaluate Maclaurin polynomial (Taylor around 0).""" + return taylor_series(x, order=order, center=0.0, function_name=function_name, function=function) + + +def estimate_lagrange_remainder(max_derivative: float, x: float | np.ndarray, order: int, center: float = 0.0) -> np.ndarray: + """Estimate Lagrange remainder bound: M|x-a|^(n+1)/(n+1)!""" + x_arr = np.asarray(x, dtype=float) + return np.abs(max_derivative) * np.abs(x_arr - center) ** (order + 1) / math.factorial(order + 1) + + +def symbolic_taylor_expression(expr_str: str, symbol: str = "x", center: float = 0.0, order: int = 6) -> str: + """Return symbolic Taylor expression as a string when SymPy is available.""" + if not HAS_SYMPY: + msg = "SymPy is not available." + raise RuntimeError(msg) + x = sp.symbols(symbol) + expr = sp.sympify(expr_str) + series = sp.series(expr, x, center, order + 1).removeO() + return str(sp.expand(series)) diff --git a/adv_math_engine/tests/__init__.py b/adv_math_engine/tests/__init__.py new file mode 100644 index 000000000000..84f9c53b8623 --- /dev/null +++ b/adv_math_engine/tests/__init__.py @@ -0,0 +1 @@ +"""Test package for adv_math_engine.""" diff --git a/adv_math_engine/tests/coverage_report.txt b/adv_math_engine/tests/coverage_report.txt new file mode 100644 index 000000000000..eb89f251b4b5 --- /dev/null +++ b/adv_math_engine/tests/coverage_report.txt @@ -0,0 +1,7 @@ +Coverage command configured: + +pytest adv_math_engine/tests --cov=adv_math_engine --cov-report=term-missing + +Execution status in this environment: +- Could not run with coverage because pytest-cov and numerical dependencies are not available in the active interpreter. +- Network-restricted package installation prevented fetching missing dependencies. diff --git a/adv_math_engine/tests/test_series_expansion.py b/adv_math_engine/tests/test_series_expansion.py new file mode 100644 index 000000000000..f3262ec16fda --- /dev/null +++ b/adv_math_engine/tests/test_series_expansion.py @@ -0,0 +1,41 @@ +from __future__ import annotations + +import math + +import numpy as np + +from adv_math_engine.series_expansion import ( + estimate_lagrange_remainder, + maclaurin_series, + taylor_series, +) +from adv_math_engine.utils import validate_tolerance + + +def test_maclaurin_exp_convergence() -> None: + x = np.linspace(-1.0, 1.0, 50) + approx = maclaurin_series(x, order=12, function_name="exp") + actual = np.exp(x) + assert validate_tolerance(approx, actual, epsilon=1e-4) + + +def test_taylor_sin_about_nonzero_center() -> None: + x = np.array([0.2, 0.3, 0.5]) + approx = taylor_series(x, order=9, center=0.3, function_name="sin") + assert np.allclose(approx, np.sin(x), atol=1e-6) + + +def test_log1p_remainder_bound() -> None: + x = 0.2 + order = 6 + approx = maclaurin_series(x, order=order, function_name="log1p") + actual = np.log1p(x) + remainder = estimate_lagrange_remainder(max_derivative=math.factorial(order), x=x, order=order) + assert abs(actual - approx) <= remainder + 1e-6 + + +def test_arbitrary_function_numeric_derivative() -> None: + f = lambda val: np.cos(val) * np.exp(val) + x = np.array([0.0, 0.1, 0.2]) + approx = taylor_series(x, order=8, center=0.0, function=f) + assert np.allclose(approx, f(x), atol=5e-3) diff --git a/adv_math_engine/tests/test_vector_algebra.py b/adv_math_engine/tests/test_vector_algebra.py new file mode 100644 index 000000000000..4119e4de1fca --- /dev/null +++ b/adv_math_engine/tests/test_vector_algebra.py @@ -0,0 +1,46 @@ +from __future__ import annotations + +import numpy as np +import pytest + +from adv_math_engine.vector_algebra import ( + angle_between, + check_linear_independence, + cross_product, + dot_product, + gram_schmidt, + vector_projection, +) + + +def test_dot_product_vectorized() -> None: + a = np.array([[1, 2, 3], [4, 5, 6]]) + b = np.array([[7, 8, 9], [1, 2, 3]]) + result = dot_product(a, b) + assert np.allclose(result, np.array([50, 32])) + + +def test_cross_product_3d() -> None: + result = cross_product([1, 0, 0], [0, 1, 0]) + assert np.allclose(result, np.array([0, 0, 1])) + + +def test_projection_rejects_zero_vector() -> None: + with pytest.raises(ValueError): + vector_projection([1, 2, 3], [0, 0, 0]) + + +def test_angle_between_orthogonal_vectors() -> None: + assert np.isclose(angle_between([1, 0], [0, 1]), np.pi / 2) + + +def test_gram_schmidt_orthonormality() -> None: + basis = gram_schmidt([[1, 1, 0], [1, 0, 1], [0, 1, 1]]) + assert np.allclose(basis @ basis.T, np.eye(3), atol=1e-10) + + +def test_linear_independence_rank_check() -> None: + independent = check_linear_independence([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + dependent = check_linear_independence([[1, 2, 3], [2, 4, 6]]) + assert independent + assert not dependent diff --git a/adv_math_engine/tests/test_vector_calculus.py b/adv_math_engine/tests/test_vector_calculus.py new file mode 100644 index 000000000000..9562fd741b31 --- /dev/null +++ b/adv_math_engine/tests/test_vector_calculus.py @@ -0,0 +1,56 @@ +from __future__ import annotations + +import numpy as np + +from adv_math_engine.vector_calculus import ( + curl, + divergence, + gradient, + line_integral, + partial_derivative, + surface_integral, +) + + +def test_partial_derivative_central() -> None: + f = lambda x, y: x**2 + y**3 + value = partial_derivative(f, [2.0, 3.0], 0, method="central") + assert np.isclose(value, 4.0, atol=1e-4) + + +def test_gradient_matches_analytic() -> None: + f = lambda x, y, z: x**2 + y**2 + z**2 + grad = gradient(f, [1.0, -2.0, 0.5]) + assert np.allclose(grad, np.array([2.0, -4.0, 1.0]), atol=1e-4) + + +def test_divergence() -> None: + field = [lambda x, y, z: x**2, lambda x, y, z: y**2, lambda x, y, z: z**2] + div = divergence(field, [1.0, 2.0, 3.0]) + assert np.isclose(div, 12.0, atol=1e-3) + + +def test_curl() -> None: + field = [lambda x, y, z: -y, lambda x, y, z: x, lambda x, y, z: 0.0] + result = curl(field, [0.2, -0.3, 1.5]) + assert np.allclose(result, np.array([0.0, 0.0, 2.0]), atol=1e-3) + + +def test_line_integral_circle() -> None: + field = lambda r: np.column_stack((-r[:, 1], r[:, 0])) + curve = lambda t: np.column_stack((np.cos(t), np.sin(t))) + value = line_integral(field, curve, 0.0, 2.0 * np.pi, num_steps=2000) + assert np.isclose(value, 2.0 * np.pi, atol=2e-2) + + +def test_surface_integral_unit_sphere_area_like() -> None: + scalar_field = lambda points: np.ones(points.shape[:2]) + + def sphere(u: np.ndarray, v: np.ndarray) -> np.ndarray: + x = np.sin(u) * np.cos(v) + y = np.sin(u) * np.sin(v) + z = np.cos(u) + return np.stack((x, y, z), axis=2) + + value = surface_integral(scalar_field, sphere, (0.0, np.pi), (0.0, 2.0 * np.pi), num_u=150, num_v=150) + assert np.isclose(value, 4.0 * np.pi, atol=0.1) diff --git a/adv_math_engine/utils.py b/adv_math_engine/utils.py new file mode 100644 index 000000000000..d569f2c7caf5 --- /dev/null +++ b/adv_math_engine/utils.py @@ -0,0 +1,27 @@ +"""Utility helpers for numerical computations.""" + +from __future__ import annotations + +from collections.abc import Callable +from time import perf_counter + +import numpy as np + +ArrayLike = float | list[float] | np.ndarray + + +def as_float_array(value: ArrayLike) -> np.ndarray: + """Convert any scalar/vector-like input to a NumPy float64 array.""" + return np.asarray(value, dtype=float) + + +def validate_tolerance(numerical: float | np.ndarray, actual: float | np.ndarray, epsilon: float = 1e-6) -> bool: + """Return True if |numerical - actual| < epsilon for all entries.""" + return bool(np.all(np.abs(as_float_array(numerical) - as_float_array(actual)) < epsilon)) + + +def timed_call(fn: Callable[..., object], *args: object, **kwargs: object) -> tuple[object, float]: + """Run a callable and return (result, elapsed_seconds).""" + start = perf_counter() + result = fn(*args, **kwargs) + return result, perf_counter() - start diff --git a/adv_math_engine/vector_algebra.py b/adv_math_engine/vector_algebra.py new file mode 100644 index 000000000000..09316acaa74b --- /dev/null +++ b/adv_math_engine/vector_algebra.py @@ -0,0 +1,88 @@ +"""Vector algebra primitives. + +Formulas: +- Dot product: a·b = Σ a_i b_i +- Projection of a onto b: proj_b(a) = ((a·b)/||b||^2) b +- Angle: θ = arccos((a·b)/(||a|| ||b||)) +""" + +from __future__ import annotations + +import numpy as np + +from adv_math_engine.utils import as_float_array + + +def dot_product(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: + """Compute dot product, supporting vectorized leading dimensions.""" + a_arr = as_float_array(a) + b_arr = as_float_array(b) + if a_arr.shape[-1] != b_arr.shape[-1]: + msg = "Dot product requires matching trailing dimensions." + raise ValueError(msg) + return np.sum(a_arr * b_arr, axis=-1) + + +def cross_product(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: + """Compute 3D cross product for (..., 3) vectors.""" + a_arr = as_float_array(a) + b_arr = as_float_array(b) + if a_arr.shape[-1] != 3 or b_arr.shape[-1] != 3: + msg = "Cross product is defined only for 3D vectors." + raise ValueError(msg) + return np.cross(a_arr, b_arr) + + +def vector_projection(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: + """Project vector a onto vector b.""" + a_arr = as_float_array(a) + b_arr = as_float_array(b) + denom = np.sum(b_arr * b_arr, axis=-1, keepdims=True) + if np.any(np.isclose(denom, 0.0)): + msg = "Projection onto zero vector is undefined." + raise ValueError(msg) + scale = np.sum(a_arr * b_arr, axis=-1, keepdims=True) / denom + return scale * b_arr + + +def angle_between(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: + """Return the angle(s) between vectors in radians.""" + a_arr = as_float_array(a) + b_arr = as_float_array(b) + a_norm = np.linalg.norm(a_arr, axis=-1) + b_norm = np.linalg.norm(b_arr, axis=-1) + if np.any(np.isclose(a_norm, 0.0)) or np.any(np.isclose(b_norm, 0.0)): + msg = "Angle with zero vector is undefined." + raise ValueError(msg) + cosine = dot_product(a_arr, b_arr) / (a_norm * b_norm) + return np.arccos(np.clip(cosine, -1.0, 1.0)) + + +def gram_schmidt(vectors: np.ndarray | list[list[float]], tol: float = 1e-12) -> np.ndarray: + """Return an orthonormal basis via Gram-Schmidt.""" + matrix = as_float_array(vectors) + if matrix.ndim != 2: + msg = "Gram-Schmidt expects a 2D array (n_vectors, dimension)." + raise ValueError(msg) + orthonormal: list[np.ndarray] = [] + for vec in matrix: + work = vec.copy() + for basis in orthonormal: + work = work - np.dot(work, basis) * basis + norm = np.linalg.norm(work) + if norm > tol: + orthonormal.append(work / norm) + if not orthonormal: + msg = "No non-zero independent vectors were provided." + raise ValueError(msg) + return np.vstack(orthonormal) + + +def check_linear_independence(vectors: np.ndarray | list[list[float]], tol: float = 1e-10) -> bool: + """Check linear independence using numerical rank.""" + matrix = as_float_array(vectors) + if matrix.ndim != 2: + msg = "Input must be 2D with shape (n_vectors, dimension)." + raise ValueError(msg) + rank = np.linalg.matrix_rank(matrix, tol=tol) + return bool(rank == matrix.shape[0]) diff --git a/adv_math_engine/vector_calculus.py b/adv_math_engine/vector_calculus.py new file mode 100644 index 000000000000..116049f6e8c5 --- /dev/null +++ b/adv_math_engine/vector_calculus.py @@ -0,0 +1,107 @@ +"""Vector calculus operators and numerical integration routines.""" + +from __future__ import annotations + +from collections.abc import Callable + +import numpy as np + +from adv_math_engine.utils import as_float_array + + +DerivativeMethod = str + + +def partial_derivative( + f: Callable[..., float], + point: np.ndarray | list[float], + var_index: int, + h: float = 1e-5, + method: DerivativeMethod = "central", +) -> float: + """Estimate partial derivative using finite differences.""" + x0 = as_float_array(point).astype(float) + direction = np.zeros_like(x0) + direction[var_index] = 1.0 + + if method == "forward": + return float((f(*(x0 + h * direction)) - f(*x0)) / h) + if method == "backward": + return float((f(*x0) - f(*(x0 - h * direction))) / h) + if method == "central": + return float((f(*(x0 + h * direction)) - f(*(x0 - h * direction))) / (2.0 * h)) + msg = "method must be one of {'forward', 'backward', 'central'}" + raise ValueError(msg) + + +def gradient( + f: Callable[..., float], point: np.ndarray | list[float], h: float = 1e-5, method: DerivativeMethod = "central" +) -> np.ndarray: + """Compute gradient ∇f at a point.""" + x0 = as_float_array(point) + return np.array([partial_derivative(f, x0, i, h=h, method=method) for i in range(x0.size)]) + + +def divergence( + field: list[Callable[..., float]], point: np.ndarray | list[float], h: float = 1e-5, method: DerivativeMethod = "central" +) -> float: + """Compute divergence ∇·F for an n-dimensional vector field.""" + x0 = as_float_array(point) + if len(field) != x0.size: + msg = "Field dimension must match point dimension." + raise ValueError(msg) + return float(sum(partial_derivative(component, x0, i, h=h, method=method) for i, component in enumerate(field))) + + +def curl(field: list[Callable[..., float]], point: np.ndarray | list[float], h: float = 1e-5) -> np.ndarray: + """Compute curl ∇×F for a 3D vector field F = (P, Q, R).""" + x0 = as_float_array(point) + if len(field) != 3 or x0.size != 3: + msg = "Curl is defined for 3D vector fields only." + raise ValueError(msg) + p, q, r = field + d_r_dy = partial_derivative(r, x0, 1, h=h) + d_q_dz = partial_derivative(q, x0, 2, h=h) + d_p_dz = partial_derivative(p, x0, 2, h=h) + d_r_dx = partial_derivative(r, x0, 0, h=h) + d_q_dx = partial_derivative(q, x0, 0, h=h) + d_p_dy = partial_derivative(p, x0, 1, h=h) + return np.array([d_r_dy - d_q_dz, d_p_dz - d_r_dx, d_q_dx - d_p_dy]) + + +def line_integral( + field: Callable[[np.ndarray], np.ndarray], + curve: Callable[[np.ndarray], np.ndarray], + t_start: float, + t_end: float, + num_steps: int = 500, +) -> float: + """Numerically approximate ∫_C F·dr using trapezoidal rule.""" + t = np.linspace(t_start, t_end, num_steps) + r = curve(t) + dr_dt = np.gradient(r, t, axis=0) + f_values = field(r) + integrand = np.sum(f_values * dr_dt, axis=1) + return float(np.trapezoid(integrand, t)) + + +def surface_integral( + scalar_field: Callable[[np.ndarray], np.ndarray], + parametrization: Callable[[np.ndarray, np.ndarray], np.ndarray], + u_bounds: tuple[float, float], + v_bounds: tuple[float, float], + num_u: int = 120, + num_v: int = 120, +) -> float: + """Approximate ∬_S f dS using parameterization r(u,v).""" + u = np.linspace(u_bounds[0], u_bounds[1], num_u) + v = np.linspace(v_bounds[0], v_bounds[1], num_v) + uu, vv = np.meshgrid(u, v, indexing="ij") + surface_points = parametrization(uu, vv) + + r_u = np.gradient(surface_points, u, axis=0) + r_v = np.gradient(surface_points, v, axis=1) + normal_mag = np.linalg.norm(np.cross(r_u, r_v), axis=2) + field_vals = scalar_field(surface_points) + integrand = field_vals * normal_mag + return float(np.trapezoid(np.trapezoid(integrand, v, axis=1), u, axis=0)) diff --git a/adv_math_engine/visualizations/plot_gradient_field.py b/adv_math_engine/visualizations/plot_gradient_field.py new file mode 100644 index 000000000000..9f3ad6759168 --- /dev/null +++ b/adv_math_engine/visualizations/plot_gradient_field.py @@ -0,0 +1,26 @@ +from __future__ import annotations + +import matplotlib.pyplot as plt +import numpy as np + + +def main() -> None: + x = np.linspace(-2, 2, 25) + y = np.linspace(-2, 2, 25) + xx, yy = np.meshgrid(x, y) + + # f(x, y) = x^2 + y^2 -> ∇f = (2x, 2y) + u = 2 * xx + v = 2 * yy + + plt.figure(figsize=(6, 5)) + plt.quiver(xx, yy, u, v) + plt.title("Gradient Field of f(x,y)=x²+y²") + plt.xlabel("x") + plt.ylabel("y") + plt.tight_layout() + plt.savefig("adv_math_engine/visualizations/gradient_field.png", dpi=160) + + +if __name__ == "__main__": + main() diff --git a/adv_math_engine/visualizations/plot_series_approximation.py b/adv_math_engine/visualizations/plot_series_approximation.py new file mode 100644 index 000000000000..5b0bd928fc37 --- /dev/null +++ b/adv_math_engine/visualizations/plot_series_approximation.py @@ -0,0 +1,28 @@ +from __future__ import annotations + +import matplotlib.pyplot as plt +import numpy as np + +from adv_math_engine.series_expansion import maclaurin_series + + +def main() -> None: + x = np.linspace(-2, 2, 300) + y_actual = np.sin(x) + + plt.figure(figsize=(7, 5)) + plt.plot(x, y_actual, label="sin(x)", linewidth=2) + for order in (1, 3, 5, 7): + y_approx = maclaurin_series(x, order=order, function_name="sin") + plt.plot(x, y_approx, label=f"Order {order}") + + plt.title("Maclaurin Approximation of sin(x)") + plt.xlabel("x") + plt.ylabel("y") + plt.legend() + plt.tight_layout() + plt.savefig("adv_math_engine/visualizations/series_approximation.png", dpi=160) + + +if __name__ == "__main__": + main() diff --git a/main.py b/main.py new file mode 100644 index 000000000000..eb503d4ead0d --- /dev/null +++ b/main.py @@ -0,0 +1,50 @@ +"""CLI for adv_math_engine series approximations. + +Example: +python main.py --function "sin(x)" --series taylor --order 5 --x 0.5 --center 0.0 +""" + +from __future__ import annotations + +import argparse + +import numpy as np + +from adv_math_engine.series_expansion import maclaurin_series, taylor_series + +FUNCTION_MAP = { + "sin(x)": ("sin", np.sin), + "cos(x)": ("cos", np.cos), + "exp(x)": ("exp", np.exp), + "log(1+x)": ("log1p", np.log1p), +} + + +def main() -> None: + parser = argparse.ArgumentParser(description="Series approximation CLI") + parser.add_argument("--function", required=True, choices=tuple(FUNCTION_MAP.keys())) + parser.add_argument("--series", required=True, choices=("taylor", "maclaurin")) + parser.add_argument("--order", required=True, type=int) + parser.add_argument("--x", required=True, type=float) + parser.add_argument("--center", type=float, default=0.0) + args = parser.parse_args() + + function_name, actual_fn = FUNCTION_MAP[args.function] + if args.series == "maclaurin": + approximation = float(maclaurin_series(args.x, args.order, function_name=function_name)) + else: + approximation = float(taylor_series(args.x, args.order, center=args.center, function_name=function_name)) + + actual = float(actual_fn(args.x)) + abs_error = abs(actual - approximation) + + print(f"Function: {args.function}") + print(f"Series: {args.series} (order={args.order}, center={args.center})") + print(f"x = {args.x}") + print(f"Approximation = {approximation:.12f}") + print(f"Actual = {actual:.12f}") + print(f"Abs Error = {abs_error:.3e}") + + +if __name__ == "__main__": + main() From 8a3a2a514ef4e8e07aa645bb9cbd0268281cc02d Mon Sep 17 00:00:00 2001 From: Mark Grayson <118820911+sahildando@users.noreply.github.com> Date: Wed, 8 Apr 2026 01:25:59 +0530 Subject: [PATCH 2/2] Revert "Add adv_math_engine package: vector algebra, vector calculus, series expansions, CLI and tests" --- adv_math_engine/README.md | 57 ---------- adv_math_engine/__init__.py | 41 ------- adv_math_engine/benchmarks/benchmark_diff.py | 51 --------- adv_math_engine/benchmarks/results.md | 17 --- adv_math_engine/examples/example_usage.py | 18 --- adv_math_engine/series_expansion.py | 102 ----------------- adv_math_engine/tests/__init__.py | 1 - adv_math_engine/tests/coverage_report.txt | 7 -- .../tests/test_series_expansion.py | 41 ------- adv_math_engine/tests/test_vector_algebra.py | 46 -------- adv_math_engine/tests/test_vector_calculus.py | 56 --------- adv_math_engine/utils.py | 27 ----- adv_math_engine/vector_algebra.py | 88 -------------- adv_math_engine/vector_calculus.py | 107 ------------------ .../visualizations/plot_gradient_field.py | 26 ----- .../plot_series_approximation.py | 28 ----- main.py | 50 -------- 17 files changed, 763 deletions(-) delete mode 100644 adv_math_engine/README.md delete mode 100644 adv_math_engine/__init__.py delete mode 100644 adv_math_engine/benchmarks/benchmark_diff.py delete mode 100644 adv_math_engine/benchmarks/results.md delete mode 100644 adv_math_engine/examples/example_usage.py delete mode 100644 adv_math_engine/series_expansion.py delete mode 100644 adv_math_engine/tests/__init__.py delete mode 100644 adv_math_engine/tests/coverage_report.txt delete mode 100644 adv_math_engine/tests/test_series_expansion.py delete mode 100644 adv_math_engine/tests/test_vector_algebra.py delete mode 100644 adv_math_engine/tests/test_vector_calculus.py delete mode 100644 adv_math_engine/utils.py delete mode 100644 adv_math_engine/vector_algebra.py delete mode 100644 adv_math_engine/vector_calculus.py delete mode 100644 adv_math_engine/visualizations/plot_gradient_field.py delete mode 100644 adv_math_engine/visualizations/plot_series_approximation.py delete mode 100644 main.py diff --git a/adv_math_engine/README.md b/adv_math_engine/README.md deleted file mode 100644 index c5abb3e9fd3b..000000000000 --- a/adv_math_engine/README.md +++ /dev/null @@ -1,57 +0,0 @@ -# adv_math_engine - -Production-grade module for: -- Vector algebra -- Vector calculus -- Taylor / Maclaurin series expansion - -## Structure - -```text -adv_math_engine/ - ├── vector_algebra.py - ├── vector_calculus.py - ├── series_expansion.py - ├── utils.py - ├── benchmarks/ - ├── visualizations/ - ├── examples/ - └── tests/ -``` - -## CLI - -```bash -python main.py --function "sin(x)" --series taylor --order 5 --x 0.5 --center 0.0 -``` - -## Visualizations - -Run: - -```bash -python adv_math_engine/visualizations/plot_gradient_field.py -python adv_math_engine/visualizations/plot_series_approximation.py -``` - -Outputs: -- `adv_math_engine/visualizations/gradient_field.png` -- `adv_math_engine/visualizations/series_approximation.png` - -## Benchmarks - -```bash -python adv_math_engine/benchmarks/benchmark_diff.py -``` - -See `adv_math_engine/benchmarks/results.md`. - -## Coverage - -Suggested command: - -```bash -pytest adv_math_engine/tests --cov=adv_math_engine --cov-report=term-missing -``` - -See `adv_math_engine/tests/coverage_report.txt`. diff --git a/adv_math_engine/__init__.py b/adv_math_engine/__init__.py deleted file mode 100644 index fefb18beeba6..000000000000 --- a/adv_math_engine/__init__.py +++ /dev/null @@ -1,41 +0,0 @@ -"""Advanced mathematics engine for vector algebra, vector calculus, and series expansions.""" - -from adv_math_engine.series_expansion import ( - estimate_lagrange_remainder, - maclaurin_series, - taylor_series, -) -from adv_math_engine.vector_algebra import ( - angle_between, - check_linear_independence, - cross_product, - dot_product, - gram_schmidt, - vector_projection, -) -from adv_math_engine.vector_calculus import ( - curl, - divergence, - gradient, - line_integral, - partial_derivative, - surface_integral, -) - -__all__ = [ - "angle_between", - "check_linear_independence", - "cross_product", - "curl", - "divergence", - "dot_product", - "estimate_lagrange_remainder", - "gradient", - "gram_schmidt", - "line_integral", - "maclaurin_series", - "partial_derivative", - "surface_integral", - "taylor_series", - "vector_projection", -] diff --git a/adv_math_engine/benchmarks/benchmark_diff.py b/adv_math_engine/benchmarks/benchmark_diff.py deleted file mode 100644 index bdaa1ab0a67d..000000000000 --- a/adv_math_engine/benchmarks/benchmark_diff.py +++ /dev/null @@ -1,51 +0,0 @@ -from __future__ import annotations - -import numpy as np - -from adv_math_engine.series_expansion import _numerical_nth_derivative -from adv_math_engine.utils import timed_call - -try: - import sympy as sp -except Exception: # pragma: no cover - sp = None - - -def numerical_benchmark(iterations: int = 2000) -> float: - func = lambda x: np.sin(x) * np.exp(x) - - def run() -> float: - total = 0.0 - for value in np.linspace(-1.0, 1.0, iterations): - total += _numerical_nth_derivative(func, n=1, at=float(value)) - return total - - _, elapsed = timed_call(run) - return elapsed - - -def symbolic_benchmark(iterations: int = 2000) -> float | None: - if sp is None: - return None - x = sp.symbols("x") - derivative_expr = sp.diff(sp.sin(x) * sp.exp(x), x) - evaluator = sp.lambdify(x, derivative_expr, "numpy") - - def run() -> float: - total = 0.0 - for value in np.linspace(-1.0, 1.0, iterations): - total += float(evaluator(value)) - return total - - _, elapsed = timed_call(run) - return elapsed - - -if __name__ == "__main__": - num_time = numerical_benchmark() - sym_time = symbolic_benchmark() - print(f"numerical_time_seconds={num_time:.6f}") - if sym_time is None: - print("symbolic_time_seconds=unavailable") - else: - print(f"symbolic_time_seconds={sym_time:.6f}") diff --git a/adv_math_engine/benchmarks/results.md b/adv_math_engine/benchmarks/results.md deleted file mode 100644 index 2d254bbb1c5c..000000000000 --- a/adv_math_engine/benchmarks/results.md +++ /dev/null @@ -1,17 +0,0 @@ -# Benchmark Results - -## Status -Benchmark execution is prepared via: - -```bash -python adv_math_engine/benchmarks/benchmark_diff.py -``` - -In this environment, benchmark execution could not be completed because required numerical dependencies (NumPy/SymPy) were unavailable to the active Python interpreter. - -## Expected Output Format - -```text -numerical_time_seconds= -symbolic_time_seconds= -``` diff --git a/adv_math_engine/examples/example_usage.py b/adv_math_engine/examples/example_usage.py deleted file mode 100644 index cbe07c5a7335..000000000000 --- a/adv_math_engine/examples/example_usage.py +++ /dev/null @@ -1,18 +0,0 @@ -from __future__ import annotations - -import numpy as np - -from adv_math_engine.series_expansion import maclaurin_series -from adv_math_engine.vector_algebra import dot_product, gram_schmidt -from adv_math_engine.vector_calculus import gradient - - -if __name__ == "__main__": - print("Dot:", dot_product([1, 2, 3], [4, 5, 6])) - print("Orthonormal basis:\n", gram_schmidt([[1, 1, 0], [1, 0, 1], [0, 1, 1]])) - - f = lambda x, y, z: x**2 + y**2 + z**2 - print("Gradient at (1,2,3):", gradient(f, [1, 2, 3])) - - x = np.array([0.1, 0.2, 0.3]) - print("sin(x) Maclaurin order=7:", maclaurin_series(x, order=7, function_name="sin")) diff --git a/adv_math_engine/series_expansion.py b/adv_math_engine/series_expansion.py deleted file mode 100644 index 17f47d93766e..000000000000 --- a/adv_math_engine/series_expansion.py +++ /dev/null @@ -1,102 +0,0 @@ -"""Taylor/Maclaurin series expansion utilities.""" - -from __future__ import annotations - -import math -from collections.abc import Callable -from typing import Literal - -import numpy as np - -HAS_SYMPY = True -try: - import sympy as sp -except Exception: # pragma: no cover - optional dependency guard - HAS_SYMPY = False - - -SupportedFunction = Literal["exp", "sin", "cos", "log1p"] - - -def _builtin_derivative_value(name: SupportedFunction, n: int, at: float) -> float: - if name == "exp": - return float(math.exp(at)) - if name == "sin": - return float(math.sin(at + n * math.pi / 2.0)) - if name == "cos": - return float(math.cos(at + n * math.pi / 2.0)) - if name == "log1p": - if n == 0: - return float(math.log1p(at)) - return float(((-1) ** (n - 1)) * math.factorial(n - 1) / (1.0 + at) ** n) - msg = f"Unsupported function: {name}" - raise ValueError(msg) - - -def _numerical_nth_derivative(func: Callable[[float], float], n: int, at: float, h: float = 1e-4) -> float: - if n == 0: - return float(func(at)) - if n == 1: - return float((func(at + h) - func(at - h)) / (2.0 * h)) - return float( - (_numerical_nth_derivative(func, n - 1, at + h, h) - _numerical_nth_derivative(func, n - 1, at - h, h)) / (2.0 * h) - ) - - -def taylor_series( - x: float | np.ndarray, - order: int, - center: float = 0.0, - *, - function_name: SupportedFunction | None = None, - function: Callable[[float], float] | None = None, -) -> np.ndarray: - """Evaluate Taylor polynomial of given order around center. - - f(x) ≈ Σ[n=0..order] f^(n)(a) / n! * (x-a)^n - """ - if function_name is None and function is None: - msg = "Provide either function_name or function." - raise ValueError(msg) - - x_arr = np.asarray(x, dtype=float) - approximation = np.zeros_like(x_arr, dtype=float) - - for n in range(order + 1): - if function_name is not None: - derivative_at_center = _builtin_derivative_value(function_name, n, center) - else: - assert function is not None - derivative_at_center = _numerical_nth_derivative(function, n, center) - - approximation = approximation + derivative_at_center / math.factorial(n) * (x_arr - center) ** n - - return approximation - - -def maclaurin_series( - x: float | np.ndarray, - order: int, - *, - function_name: SupportedFunction | None = None, - function: Callable[[float], float] | None = None, -) -> np.ndarray: - """Evaluate Maclaurin polynomial (Taylor around 0).""" - return taylor_series(x, order=order, center=0.0, function_name=function_name, function=function) - - -def estimate_lagrange_remainder(max_derivative: float, x: float | np.ndarray, order: int, center: float = 0.0) -> np.ndarray: - """Estimate Lagrange remainder bound: M|x-a|^(n+1)/(n+1)!""" - x_arr = np.asarray(x, dtype=float) - return np.abs(max_derivative) * np.abs(x_arr - center) ** (order + 1) / math.factorial(order + 1) - - -def symbolic_taylor_expression(expr_str: str, symbol: str = "x", center: float = 0.0, order: int = 6) -> str: - """Return symbolic Taylor expression as a string when SymPy is available.""" - if not HAS_SYMPY: - msg = "SymPy is not available." - raise RuntimeError(msg) - x = sp.symbols(symbol) - expr = sp.sympify(expr_str) - series = sp.series(expr, x, center, order + 1).removeO() - return str(sp.expand(series)) diff --git a/adv_math_engine/tests/__init__.py b/adv_math_engine/tests/__init__.py deleted file mode 100644 index 84f9c53b8623..000000000000 --- a/adv_math_engine/tests/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Test package for adv_math_engine.""" diff --git a/adv_math_engine/tests/coverage_report.txt b/adv_math_engine/tests/coverage_report.txt deleted file mode 100644 index eb89f251b4b5..000000000000 --- a/adv_math_engine/tests/coverage_report.txt +++ /dev/null @@ -1,7 +0,0 @@ -Coverage command configured: - -pytest adv_math_engine/tests --cov=adv_math_engine --cov-report=term-missing - -Execution status in this environment: -- Could not run with coverage because pytest-cov and numerical dependencies are not available in the active interpreter. -- Network-restricted package installation prevented fetching missing dependencies. diff --git a/adv_math_engine/tests/test_series_expansion.py b/adv_math_engine/tests/test_series_expansion.py deleted file mode 100644 index f3262ec16fda..000000000000 --- a/adv_math_engine/tests/test_series_expansion.py +++ /dev/null @@ -1,41 +0,0 @@ -from __future__ import annotations - -import math - -import numpy as np - -from adv_math_engine.series_expansion import ( - estimate_lagrange_remainder, - maclaurin_series, - taylor_series, -) -from adv_math_engine.utils import validate_tolerance - - -def test_maclaurin_exp_convergence() -> None: - x = np.linspace(-1.0, 1.0, 50) - approx = maclaurin_series(x, order=12, function_name="exp") - actual = np.exp(x) - assert validate_tolerance(approx, actual, epsilon=1e-4) - - -def test_taylor_sin_about_nonzero_center() -> None: - x = np.array([0.2, 0.3, 0.5]) - approx = taylor_series(x, order=9, center=0.3, function_name="sin") - assert np.allclose(approx, np.sin(x), atol=1e-6) - - -def test_log1p_remainder_bound() -> None: - x = 0.2 - order = 6 - approx = maclaurin_series(x, order=order, function_name="log1p") - actual = np.log1p(x) - remainder = estimate_lagrange_remainder(max_derivative=math.factorial(order), x=x, order=order) - assert abs(actual - approx) <= remainder + 1e-6 - - -def test_arbitrary_function_numeric_derivative() -> None: - f = lambda val: np.cos(val) * np.exp(val) - x = np.array([0.0, 0.1, 0.2]) - approx = taylor_series(x, order=8, center=0.0, function=f) - assert np.allclose(approx, f(x), atol=5e-3) diff --git a/adv_math_engine/tests/test_vector_algebra.py b/adv_math_engine/tests/test_vector_algebra.py deleted file mode 100644 index 4119e4de1fca..000000000000 --- a/adv_math_engine/tests/test_vector_algebra.py +++ /dev/null @@ -1,46 +0,0 @@ -from __future__ import annotations - -import numpy as np -import pytest - -from adv_math_engine.vector_algebra import ( - angle_between, - check_linear_independence, - cross_product, - dot_product, - gram_schmidt, - vector_projection, -) - - -def test_dot_product_vectorized() -> None: - a = np.array([[1, 2, 3], [4, 5, 6]]) - b = np.array([[7, 8, 9], [1, 2, 3]]) - result = dot_product(a, b) - assert np.allclose(result, np.array([50, 32])) - - -def test_cross_product_3d() -> None: - result = cross_product([1, 0, 0], [0, 1, 0]) - assert np.allclose(result, np.array([0, 0, 1])) - - -def test_projection_rejects_zero_vector() -> None: - with pytest.raises(ValueError): - vector_projection([1, 2, 3], [0, 0, 0]) - - -def test_angle_between_orthogonal_vectors() -> None: - assert np.isclose(angle_between([1, 0], [0, 1]), np.pi / 2) - - -def test_gram_schmidt_orthonormality() -> None: - basis = gram_schmidt([[1, 1, 0], [1, 0, 1], [0, 1, 1]]) - assert np.allclose(basis @ basis.T, np.eye(3), atol=1e-10) - - -def test_linear_independence_rank_check() -> None: - independent = check_linear_independence([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) - dependent = check_linear_independence([[1, 2, 3], [2, 4, 6]]) - assert independent - assert not dependent diff --git a/adv_math_engine/tests/test_vector_calculus.py b/adv_math_engine/tests/test_vector_calculus.py deleted file mode 100644 index 9562fd741b31..000000000000 --- a/adv_math_engine/tests/test_vector_calculus.py +++ /dev/null @@ -1,56 +0,0 @@ -from __future__ import annotations - -import numpy as np - -from adv_math_engine.vector_calculus import ( - curl, - divergence, - gradient, - line_integral, - partial_derivative, - surface_integral, -) - - -def test_partial_derivative_central() -> None: - f = lambda x, y: x**2 + y**3 - value = partial_derivative(f, [2.0, 3.0], 0, method="central") - assert np.isclose(value, 4.0, atol=1e-4) - - -def test_gradient_matches_analytic() -> None: - f = lambda x, y, z: x**2 + y**2 + z**2 - grad = gradient(f, [1.0, -2.0, 0.5]) - assert np.allclose(grad, np.array([2.0, -4.0, 1.0]), atol=1e-4) - - -def test_divergence() -> None: - field = [lambda x, y, z: x**2, lambda x, y, z: y**2, lambda x, y, z: z**2] - div = divergence(field, [1.0, 2.0, 3.0]) - assert np.isclose(div, 12.0, atol=1e-3) - - -def test_curl() -> None: - field = [lambda x, y, z: -y, lambda x, y, z: x, lambda x, y, z: 0.0] - result = curl(field, [0.2, -0.3, 1.5]) - assert np.allclose(result, np.array([0.0, 0.0, 2.0]), atol=1e-3) - - -def test_line_integral_circle() -> None: - field = lambda r: np.column_stack((-r[:, 1], r[:, 0])) - curve = lambda t: np.column_stack((np.cos(t), np.sin(t))) - value = line_integral(field, curve, 0.0, 2.0 * np.pi, num_steps=2000) - assert np.isclose(value, 2.0 * np.pi, atol=2e-2) - - -def test_surface_integral_unit_sphere_area_like() -> None: - scalar_field = lambda points: np.ones(points.shape[:2]) - - def sphere(u: np.ndarray, v: np.ndarray) -> np.ndarray: - x = np.sin(u) * np.cos(v) - y = np.sin(u) * np.sin(v) - z = np.cos(u) - return np.stack((x, y, z), axis=2) - - value = surface_integral(scalar_field, sphere, (0.0, np.pi), (0.0, 2.0 * np.pi), num_u=150, num_v=150) - assert np.isclose(value, 4.0 * np.pi, atol=0.1) diff --git a/adv_math_engine/utils.py b/adv_math_engine/utils.py deleted file mode 100644 index d569f2c7caf5..000000000000 --- a/adv_math_engine/utils.py +++ /dev/null @@ -1,27 +0,0 @@ -"""Utility helpers for numerical computations.""" - -from __future__ import annotations - -from collections.abc import Callable -from time import perf_counter - -import numpy as np - -ArrayLike = float | list[float] | np.ndarray - - -def as_float_array(value: ArrayLike) -> np.ndarray: - """Convert any scalar/vector-like input to a NumPy float64 array.""" - return np.asarray(value, dtype=float) - - -def validate_tolerance(numerical: float | np.ndarray, actual: float | np.ndarray, epsilon: float = 1e-6) -> bool: - """Return True if |numerical - actual| < epsilon for all entries.""" - return bool(np.all(np.abs(as_float_array(numerical) - as_float_array(actual)) < epsilon)) - - -def timed_call(fn: Callable[..., object], *args: object, **kwargs: object) -> tuple[object, float]: - """Run a callable and return (result, elapsed_seconds).""" - start = perf_counter() - result = fn(*args, **kwargs) - return result, perf_counter() - start diff --git a/adv_math_engine/vector_algebra.py b/adv_math_engine/vector_algebra.py deleted file mode 100644 index 09316acaa74b..000000000000 --- a/adv_math_engine/vector_algebra.py +++ /dev/null @@ -1,88 +0,0 @@ -"""Vector algebra primitives. - -Formulas: -- Dot product: a·b = Σ a_i b_i -- Projection of a onto b: proj_b(a) = ((a·b)/||b||^2) b -- Angle: θ = arccos((a·b)/(||a|| ||b||)) -""" - -from __future__ import annotations - -import numpy as np - -from adv_math_engine.utils import as_float_array - - -def dot_product(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: - """Compute dot product, supporting vectorized leading dimensions.""" - a_arr = as_float_array(a) - b_arr = as_float_array(b) - if a_arr.shape[-1] != b_arr.shape[-1]: - msg = "Dot product requires matching trailing dimensions." - raise ValueError(msg) - return np.sum(a_arr * b_arr, axis=-1) - - -def cross_product(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: - """Compute 3D cross product for (..., 3) vectors.""" - a_arr = as_float_array(a) - b_arr = as_float_array(b) - if a_arr.shape[-1] != 3 or b_arr.shape[-1] != 3: - msg = "Cross product is defined only for 3D vectors." - raise ValueError(msg) - return np.cross(a_arr, b_arr) - - -def vector_projection(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: - """Project vector a onto vector b.""" - a_arr = as_float_array(a) - b_arr = as_float_array(b) - denom = np.sum(b_arr * b_arr, axis=-1, keepdims=True) - if np.any(np.isclose(denom, 0.0)): - msg = "Projection onto zero vector is undefined." - raise ValueError(msg) - scale = np.sum(a_arr * b_arr, axis=-1, keepdims=True) / denom - return scale * b_arr - - -def angle_between(a: np.ndarray | list[float], b: np.ndarray | list[float]) -> np.ndarray: - """Return the angle(s) between vectors in radians.""" - a_arr = as_float_array(a) - b_arr = as_float_array(b) - a_norm = np.linalg.norm(a_arr, axis=-1) - b_norm = np.linalg.norm(b_arr, axis=-1) - if np.any(np.isclose(a_norm, 0.0)) or np.any(np.isclose(b_norm, 0.0)): - msg = "Angle with zero vector is undefined." - raise ValueError(msg) - cosine = dot_product(a_arr, b_arr) / (a_norm * b_norm) - return np.arccos(np.clip(cosine, -1.0, 1.0)) - - -def gram_schmidt(vectors: np.ndarray | list[list[float]], tol: float = 1e-12) -> np.ndarray: - """Return an orthonormal basis via Gram-Schmidt.""" - matrix = as_float_array(vectors) - if matrix.ndim != 2: - msg = "Gram-Schmidt expects a 2D array (n_vectors, dimension)." - raise ValueError(msg) - orthonormal: list[np.ndarray] = [] - for vec in matrix: - work = vec.copy() - for basis in orthonormal: - work = work - np.dot(work, basis) * basis - norm = np.linalg.norm(work) - if norm > tol: - orthonormal.append(work / norm) - if not orthonormal: - msg = "No non-zero independent vectors were provided." - raise ValueError(msg) - return np.vstack(orthonormal) - - -def check_linear_independence(vectors: np.ndarray | list[list[float]], tol: float = 1e-10) -> bool: - """Check linear independence using numerical rank.""" - matrix = as_float_array(vectors) - if matrix.ndim != 2: - msg = "Input must be 2D with shape (n_vectors, dimension)." - raise ValueError(msg) - rank = np.linalg.matrix_rank(matrix, tol=tol) - return bool(rank == matrix.shape[0]) diff --git a/adv_math_engine/vector_calculus.py b/adv_math_engine/vector_calculus.py deleted file mode 100644 index 116049f6e8c5..000000000000 --- a/adv_math_engine/vector_calculus.py +++ /dev/null @@ -1,107 +0,0 @@ -"""Vector calculus operators and numerical integration routines.""" - -from __future__ import annotations - -from collections.abc import Callable - -import numpy as np - -from adv_math_engine.utils import as_float_array - - -DerivativeMethod = str - - -def partial_derivative( - f: Callable[..., float], - point: np.ndarray | list[float], - var_index: int, - h: float = 1e-5, - method: DerivativeMethod = "central", -) -> float: - """Estimate partial derivative using finite differences.""" - x0 = as_float_array(point).astype(float) - direction = np.zeros_like(x0) - direction[var_index] = 1.0 - - if method == "forward": - return float((f(*(x0 + h * direction)) - f(*x0)) / h) - if method == "backward": - return float((f(*x0) - f(*(x0 - h * direction))) / h) - if method == "central": - return float((f(*(x0 + h * direction)) - f(*(x0 - h * direction))) / (2.0 * h)) - msg = "method must be one of {'forward', 'backward', 'central'}" - raise ValueError(msg) - - -def gradient( - f: Callable[..., float], point: np.ndarray | list[float], h: float = 1e-5, method: DerivativeMethod = "central" -) -> np.ndarray: - """Compute gradient ∇f at a point.""" - x0 = as_float_array(point) - return np.array([partial_derivative(f, x0, i, h=h, method=method) for i in range(x0.size)]) - - -def divergence( - field: list[Callable[..., float]], point: np.ndarray | list[float], h: float = 1e-5, method: DerivativeMethod = "central" -) -> float: - """Compute divergence ∇·F for an n-dimensional vector field.""" - x0 = as_float_array(point) - if len(field) != x0.size: - msg = "Field dimension must match point dimension." - raise ValueError(msg) - return float(sum(partial_derivative(component, x0, i, h=h, method=method) for i, component in enumerate(field))) - - -def curl(field: list[Callable[..., float]], point: np.ndarray | list[float], h: float = 1e-5) -> np.ndarray: - """Compute curl ∇×F for a 3D vector field F = (P, Q, R).""" - x0 = as_float_array(point) - if len(field) != 3 or x0.size != 3: - msg = "Curl is defined for 3D vector fields only." - raise ValueError(msg) - p, q, r = field - d_r_dy = partial_derivative(r, x0, 1, h=h) - d_q_dz = partial_derivative(q, x0, 2, h=h) - d_p_dz = partial_derivative(p, x0, 2, h=h) - d_r_dx = partial_derivative(r, x0, 0, h=h) - d_q_dx = partial_derivative(q, x0, 0, h=h) - d_p_dy = partial_derivative(p, x0, 1, h=h) - return np.array([d_r_dy - d_q_dz, d_p_dz - d_r_dx, d_q_dx - d_p_dy]) - - -def line_integral( - field: Callable[[np.ndarray], np.ndarray], - curve: Callable[[np.ndarray], np.ndarray], - t_start: float, - t_end: float, - num_steps: int = 500, -) -> float: - """Numerically approximate ∫_C F·dr using trapezoidal rule.""" - t = np.linspace(t_start, t_end, num_steps) - r = curve(t) - dr_dt = np.gradient(r, t, axis=0) - f_values = field(r) - integrand = np.sum(f_values * dr_dt, axis=1) - return float(np.trapezoid(integrand, t)) - - -def surface_integral( - scalar_field: Callable[[np.ndarray], np.ndarray], - parametrization: Callable[[np.ndarray, np.ndarray], np.ndarray], - u_bounds: tuple[float, float], - v_bounds: tuple[float, float], - num_u: int = 120, - num_v: int = 120, -) -> float: - """Approximate ∬_S f dS using parameterization r(u,v).""" - u = np.linspace(u_bounds[0], u_bounds[1], num_u) - v = np.linspace(v_bounds[0], v_bounds[1], num_v) - uu, vv = np.meshgrid(u, v, indexing="ij") - surface_points = parametrization(uu, vv) - - r_u = np.gradient(surface_points, u, axis=0) - r_v = np.gradient(surface_points, v, axis=1) - normal_mag = np.linalg.norm(np.cross(r_u, r_v), axis=2) - field_vals = scalar_field(surface_points) - integrand = field_vals * normal_mag - return float(np.trapezoid(np.trapezoid(integrand, v, axis=1), u, axis=0)) diff --git a/adv_math_engine/visualizations/plot_gradient_field.py b/adv_math_engine/visualizations/plot_gradient_field.py deleted file mode 100644 index 9f3ad6759168..000000000000 --- a/adv_math_engine/visualizations/plot_gradient_field.py +++ /dev/null @@ -1,26 +0,0 @@ -from __future__ import annotations - -import matplotlib.pyplot as plt -import numpy as np - - -def main() -> None: - x = np.linspace(-2, 2, 25) - y = np.linspace(-2, 2, 25) - xx, yy = np.meshgrid(x, y) - - # f(x, y) = x^2 + y^2 -> ∇f = (2x, 2y) - u = 2 * xx - v = 2 * yy - - plt.figure(figsize=(6, 5)) - plt.quiver(xx, yy, u, v) - plt.title("Gradient Field of f(x,y)=x²+y²") - plt.xlabel("x") - plt.ylabel("y") - plt.tight_layout() - plt.savefig("adv_math_engine/visualizations/gradient_field.png", dpi=160) - - -if __name__ == "__main__": - main() diff --git a/adv_math_engine/visualizations/plot_series_approximation.py b/adv_math_engine/visualizations/plot_series_approximation.py deleted file mode 100644 index 5b0bd928fc37..000000000000 --- a/adv_math_engine/visualizations/plot_series_approximation.py +++ /dev/null @@ -1,28 +0,0 @@ -from __future__ import annotations - -import matplotlib.pyplot as plt -import numpy as np - -from adv_math_engine.series_expansion import maclaurin_series - - -def main() -> None: - x = np.linspace(-2, 2, 300) - y_actual = np.sin(x) - - plt.figure(figsize=(7, 5)) - plt.plot(x, y_actual, label="sin(x)", linewidth=2) - for order in (1, 3, 5, 7): - y_approx = maclaurin_series(x, order=order, function_name="sin") - plt.plot(x, y_approx, label=f"Order {order}") - - plt.title("Maclaurin Approximation of sin(x)") - plt.xlabel("x") - plt.ylabel("y") - plt.legend() - plt.tight_layout() - plt.savefig("adv_math_engine/visualizations/series_approximation.png", dpi=160) - - -if __name__ == "__main__": - main() diff --git a/main.py b/main.py deleted file mode 100644 index eb503d4ead0d..000000000000 --- a/main.py +++ /dev/null @@ -1,50 +0,0 @@ -"""CLI for adv_math_engine series approximations. - -Example: -python main.py --function "sin(x)" --series taylor --order 5 --x 0.5 --center 0.0 -""" - -from __future__ import annotations - -import argparse - -import numpy as np - -from adv_math_engine.series_expansion import maclaurin_series, taylor_series - -FUNCTION_MAP = { - "sin(x)": ("sin", np.sin), - "cos(x)": ("cos", np.cos), - "exp(x)": ("exp", np.exp), - "log(1+x)": ("log1p", np.log1p), -} - - -def main() -> None: - parser = argparse.ArgumentParser(description="Series approximation CLI") - parser.add_argument("--function", required=True, choices=tuple(FUNCTION_MAP.keys())) - parser.add_argument("--series", required=True, choices=("taylor", "maclaurin")) - parser.add_argument("--order", required=True, type=int) - parser.add_argument("--x", required=True, type=float) - parser.add_argument("--center", type=float, default=0.0) - args = parser.parse_args() - - function_name, actual_fn = FUNCTION_MAP[args.function] - if args.series == "maclaurin": - approximation = float(maclaurin_series(args.x, args.order, function_name=function_name)) - else: - approximation = float(taylor_series(args.x, args.order, center=args.center, function_name=function_name)) - - actual = float(actual_fn(args.x)) - abs_error = abs(actual - approximation) - - print(f"Function: {args.function}") - print(f"Series: {args.series} (order={args.order}, center={args.center})") - print(f"x = {args.x}") - print(f"Approximation = {approximation:.12f}") - print(f"Actual = {actual:.12f}") - print(f"Abs Error = {abs_error:.3e}") - - -if __name__ == "__main__": - main()