diff --git a/scijo/differentiate/__init__.mojo b/scijo/differentiate/__init__.mojo index fbe2eb1..4398962 100644 --- a/scijo/differentiate/__init__.mojo +++ b/scijo/differentiate/__init__.mojo @@ -8,7 +8,8 @@ """Differentiate Module (scijo.differentiate) The `differentiate` module provides tools for numerical differentiation and gradient computation. -It includes functions for calculating derivatives, Jacobians, and other related operations essential for scientific computing tasks such as optimization, sensitivity analysis, and solving differential equations. +It includes functions for calculating derivatives, Jacobians, with much more to come in the future. """ -from .derivative import derivative -from .jacobian import jacobian + +from .deriv import derivative +from .jacob import jacobian diff --git a/scijo/differentiate/derivative.mojo b/scijo/differentiate/deriv.mojo similarity index 96% rename from scijo/differentiate/derivative.mojo rename to scijo/differentiate/deriv.mojo index 921642a..ac55a01 100644 --- a/scijo/differentiate/derivative.mojo +++ b/scijo/differentiate/deriv.mojo @@ -17,7 +17,6 @@ References: https://en.wikipedia.org/wiki/Finite_difference_coefficient """ -from numojo.prelude import * from .utility import ( DiffResult, @@ -27,6 +26,11 @@ from .utility import ( ) +# ===----------------------------------------------------------------------=== # +# Derivative +# ===----------------------------------------------------------------------=== # + + fn derivative[ dtype: DType, func: fn[dtype: DType]( @@ -72,6 +76,17 @@ fn derivative[ Raises: Error: If step_direction is not in {-1, 0, 1}. Error: If the specified order is not supported for the chosen method. + + Examples: + ```mojo + from scijo.differentiate import derivative + from scijo.prelude import * + + fn f[dtype: DType](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype]: + return x * x + + var res = derivative[f64, f](1.0) + ``` """ @parameter @@ -128,6 +143,11 @@ fn derivative[ ) +# ===----------------------------------------------------------------------=== # +# Internal methods +# ===----------------------------------------------------------------------=== # + + fn _derivative_central_difference[ dtype: DType, func: fn[dtype: DType]( diff --git a/scijo/differentiate/jacobian.mojo b/scijo/differentiate/jacob.mojo similarity index 89% rename from scijo/differentiate/jacobian.mojo rename to scijo/differentiate/jacob.mojo index 4331c28..c080e99 100644 --- a/scijo/differentiate/jacobian.mojo +++ b/scijo/differentiate/jacob.mojo @@ -44,12 +44,25 @@ fn jacobian[ tolerances: Tolerance dictionary with "abs" and "rel" keys (reserved for future use). maxiter: Maximum iterations (reserved for future use). + Raises: + Error: If function evaluation fails for any perturbation. + Returns: NDArray[dtype] of shape (m, n) representing the Jacobian matrix, where m is the output dimension and n is the input dimension. - Raises: - Error: If function evaluation fails for any perturbation. + Examples: + ```mojo + import numojo as nm + from scijo.differentiate import jacobian + from scijo.prelude import * + + fn f[dtype: DType](x: NDArray[dtype], args: Optional[List[Scalar[dtype]]]) raises -> NDArray[dtype]: + return x * x + + var x = nm.array[f64]([1.0, 2.0]) + var J = jacobian[f64, f](x) + ``` """ var n: Int = len(x) var f0: NDArray[dtype] = f(x, args) diff --git a/scijo/differentiate/utility.mojo b/scijo/differentiate/utility.mojo index 53181b4..9c2d362 100644 --- a/scijo/differentiate/utility.mojo +++ b/scijo/differentiate/utility.mojo @@ -18,6 +18,10 @@ References: Arbitrarily Spaced Grids. Mathematics of Computation, 51(184), 699-706. """ +# ===----------------------------------------------------------------------=== # +# Result types +# ===----------------------------------------------------------------------=== # + struct DiffResult[dtype: DType](ImplicitlyCopyable, Writable): """Result structure for numerical differentiation operations. @@ -88,6 +92,11 @@ struct DiffResult[dtype: DType](ImplicitlyCopyable, Writable): writer.write("Error displaying Result: " + String(e) + "\n") +# ===----------------------------------------------------------------------=== # +# Finite difference tables +# ===----------------------------------------------------------------------=== # + + @parameter fn generate_central_finite_difference_table[ dtype: DType diff --git a/scijo/integrate/fixed_sample.mojo b/scijo/integrate/fixed_sample.mojo index ecc2a84..133bf50 100644 --- a/scijo/integrate/fixed_sample.mojo +++ b/scijo/integrate/fixed_sample.mojo @@ -19,6 +19,7 @@ import numojo as nm # Trapezoid # ===----------------------------------------------------------------------=== # + fn trapezoid[ dtype: DType ]( diff --git a/scijo/optimize/min_scalar.mojo b/scijo/optimize/min_scalar.mojo index ad5e382..03359f8 100644 --- a/scijo/optimize/min_scalar.mojo +++ b/scijo/optimize/min_scalar.mojo @@ -103,7 +103,6 @@ fn _brent_minimize[ var tmpx = a var tmpf = fa a = b - fa = fb b = tmpx fb = tmpf @@ -113,7 +112,6 @@ fn _brent_minimize[ var bracket_iter = 0 while fb > fc and bracket_iter < 50: a = b - fa = fb b = c fb = fc c = b + (b - a) * Scalar[dtype](_phi) @@ -151,13 +149,10 @@ fn _brent_minimize[ nfev=nfev, ) - var p: Scalar[dtype] = 0 - var q: Scalar[dtype] = 0 - var r: Scalar[dtype] = 0 if abs(e) > tol1: - r = (x - w) * (fx - fv) - q = (x - v) * (fx - fw) - p = (x - v) * q - (x - w) * r + var r = (x - w) * (fx - fv) + var q = (x - v) * (fx - fw) + var p = (x - v) * q - (x - w) * r q = (q - r) * 2 if q > 0: p = -p @@ -167,7 +162,7 @@ fn _brent_minimize[ and p > q * (a - x) and p < q * (c - x) ): - d = p / q + var d = p / q var u1 = x + d if (u1 - a) < tol2 or (c - u1) < tol2: d = tol1 if x < m else -tol1 @@ -336,13 +331,10 @@ fn _bounded_minimize[ nfev=nfev, ) - var p: Scalar[dtype] = 0 - var q: Scalar[dtype] = 0 - var r: Scalar[dtype] = 0 if abs(e) > tol1: - r = (x - w) * (fx - fv) - q = (x - v) * (fx - fw) - p = (x - v) * q - (x - w) * r + var r = (x - w) * (fx - fv) + var q = (x - v) * (fx - fw) + var p = (x - v) * q - (x - w) * r q = (q - r) * 2 if q > 0: p = -p diff --git a/scijo/optimize/root_scalar.mojo b/scijo/optimize/root_scalar.mojo index eda4365..74ce5a3 100644 --- a/scijo/optimize/root_scalar.mojo +++ b/scijo/optimize/root_scalar.mojo @@ -18,7 +18,6 @@ methods (secant). # Root scalar # ===----------------------------------------------------------------------=== # - fn root_scalar[ dtype: DType, f: fn[dtype: DType]( @@ -70,37 +69,29 @@ fn root_scalar[ if method == "newton": if not fprime: raise Error( - "Scijo [root_scalar]: Derivative fprime must be provided for" - " Newton's method." + "Scijo [root_scalar]: Derivative fprime must be provided for Newton's method." ) return newton[dtype, f, fprime.value()](args, x0, xtol, rtol, maxiter) elif method == "bisect": if not bracket: - raise Error( - "Scijo [root_scalar]: Bracket must be provided for bisection" - " method." - ) + raise Error("Scijo [root_scalar]: Bracket must be provided for bisection method.") return bisect[dtype, f](args, bracket.value(), xtol, rtol, maxiter) elif method == "secant": if not (x0 and x1): raise Error( - "Scijo [root_scalar]: Initial guesses x0 and x1 must be" - " provided for secant method." + "Scijo [root_scalar]: Initial guesses x0 and x1 must be provided for secant method." ) return secant[dtype, f]( args, x0.value(), x1.value(), xtol, rtol, maxiter ) else: - raise Error( - "Scijo [root_scalar]: Unsupported method: " + String(method) - ) + raise Error("Scijo [root_scalar]: Unsupported method: " + String(method)) # ===----------------------------------------------------------------------=== # # Root scalar methods # ===----------------------------------------------------------------------=== # - fn newton[ dtype: DType, f: fn[dtype: DType]( @@ -145,10 +136,7 @@ fn newton[ if x0: xn = x0.value() else: - raise Error( - "Scijo [newton]: Initial guess x0 must be provided for Newton's" - " method." - ) + raise Error("Scijo [newton]: Initial guess x0 must be provided for Newton's method.") for _ in range(maxiter): var fx = f(xn, args) @@ -156,8 +144,7 @@ fn newton[ if fpx == 0: raise Error( - "Scijo [newton]: Derivative is zero. Newton-Raphson step would" - " divide by zero." + "Scijo [newton]: Derivative is zero. Newton-Raphson step would divide by zero." ) var delta = fx / fpx @@ -224,8 +211,8 @@ fn bisect[ if fa * fb > 0: raise Error( - "Scijo [newton]: f(a) and f(b) must have opposite signs (bracket" - " does not enclose a root)." + "Scijo [newton]: f(a) and f(b) must have opposite signs (bracket does not enclose a" + " root)." ) for _ in range(maxiter): @@ -294,10 +281,7 @@ fn secant[ var denom = f1 - f0 if denom == 0: - raise Error( - "Scijo [newton]: Secant method encountered zero slope (f1 - f0" - " == 0)." - ) + raise Error("Scijo [newton]: Secant method encountered zero slope (f1 - f0 == 0).") var xn = b - (f1 * (b - a)) / denom diff --git a/scijo/optimize/utility.mojo b/scijo/optimize/utility.mojo index 48d043a..8cd7eb3 100644 --- a/scijo/optimize/utility.mojo +++ b/scijo/optimize/utility.mojo @@ -14,7 +14,6 @@ Data structures for returning results from optimization and root-finding routine # RootResults # ===----------------------------------------------------------------------=== # - struct RootResults[dtype: DType = DType.float64](): """Result structure for scalar root-finding operations. diff --git a/tests/test_differentiate.mojo b/tests/test_differentiate.mojo index 91b115a..c9b5d7a 100644 --- a/tests/test_differentiate.mojo +++ b/tests/test_differentiate.mojo @@ -1,4 +1,4 @@ -from scijo.differentiate.derivative import derivative +from scijo.differentiate import derivative from testing import assert_almost_equal, assert_equal, assert_true, assert_false from testing import TestSuite diff --git a/tests/test_quad.mojo b/tests/test_quad.mojo index da13f41..166f5a3 100644 --- a/tests/test_quad.mojo +++ b/tests/test_quad.mojo @@ -11,7 +11,7 @@ from testing import assert_almost_equal, assert_equal, assert_true, assert_false from testing import TestSuite from math import sin, cos, exp, log, pi, sqrt -from scijo.integrate.quad import quad +from scijo.integrate import quad import scijo as sj