From 69902828d2370722c7b0199e28c8ca329c8a0d9b Mon Sep 17 00:00:00 2001 From: shivasankar Date: Fri, 13 Mar 2026 22:38:21 +0900 Subject: [PATCH 1/4] remove dtype param from derivative module --- scijo/differentiate/deriv.mojo | 184 ++++++++++++++----------------- scijo/differentiate/jacob.mojo | 34 +++--- scijo/differentiate/utility.mojo | 61 ++++------ 3 files changed, 124 insertions(+), 155 deletions(-) diff --git a/scijo/differentiate/deriv.mojo b/scijo/differentiate/deriv.mojo index ac55a01..8be2ad4 100644 --- a/scijo/differentiate/deriv.mojo +++ b/scijo/differentiate/deriv.mojo @@ -32,21 +32,20 @@ from .utility import ( fn derivative[ - dtype: DType, - func: fn[dtype: DType]( - x: Scalar[dtype], args: Optional[List[Scalar[dtype]]] - ) -> Scalar[dtype], + func: fn( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] + ) -> Scalar[f64], *, step_direction: Int = 0, ]( - x0: Scalar[dtype], - args: Optional[List[Scalar[dtype]]] = None, - tolerances: Dict[String, Scalar[dtype]] = {"atol": 1e-6, "rtol": 1e-6}, + x0: Scalar[f64], + args: Optional[List[Scalar[f64]]] = None, + tolerances: Dict[String, Scalar[f64]] = {"atol": 1e-6, "rtol": 1e-6}, max_iter: Int = 10, order: Int = 8, - initial_step: Scalar[dtype] = 0.5, - step_factor: Scalar[dtype] = 2.0, -) raises -> DiffResult[dtype] where dtype.is_floating_point(): + initial_step: Scalar[f64] = 0.5, + step_factor: Scalar[f64] = 2.0, +) raises -> DiffResult[f64]: """Computes the first derivative of a scalar function using finite differences. Provides a unified interface for computing first-order derivatives using @@ -54,8 +53,7 @@ fn derivative[ size reduction for improved accuracy. Parameters: - dtype: The floating-point data type. - func: Function to differentiate with signature fn(x, args) -> Scalar[dtype]. + func: Function to differentiate with signature fn(x, args) -> Scalar[f64]. step_direction: Direction of finite difference (central=0, forward=1, backward=-1). Keyword-only. @@ -70,7 +68,7 @@ fn derivative[ step_factor: Factor by which to reduce step size in each iteration (must be > 1). Returns: - DiffResult[dtype] containing the derivative, convergence status, + DiffResult[f64] containing the derivative, convergence status, estimated error, iteration count, and function evaluation count. Raises: @@ -82,19 +80,17 @@ fn derivative[ from scijo.differentiate import derivative from scijo.prelude import * - fn f[dtype: DType](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype]: + fn f(x: Scalar[f64], args: Optional[List[Scalar[f64]]]) -> Scalar[f64]: return x * x - var res = derivative[f64, f](1.0) + var res = derivative[f](1.0) ``` """ @parameter if step_direction == 0: - comptime first_order_coefficients = generate_central_finite_difference_table[ - dtype - ]() - return _derivative_central_difference[dtype, func]( + comptime first_order_coefficients = generate_central_finite_difference_table() + return _derivative_central_difference[func]( x0, args, tolerances=tolerances, @@ -104,10 +100,8 @@ fn derivative[ max_iter=max_iter, ) elif step_direction == 1: - comptime first_order_coefficients = generate_forward_finite_difference_table[ - dtype - ]() - return _derivative_forward_difference[dtype, func]( + comptime first_order_coefficients = generate_forward_finite_difference_table() + return _derivative_forward_difference[func]( x0, args, tolerances, @@ -117,10 +111,8 @@ fn derivative[ max_iter, ) elif step_direction == -1: - comptime first_order_coefficients = generate_backward_finite_difference_table[ - dtype - ]() - return _derivative_backward_difference[dtype, func]( + comptime first_order_coefficients = generate_backward_finite_difference_table() + return _derivative_backward_difference[func]( x0, args, tolerances, @@ -149,27 +141,25 @@ fn derivative[ fn _derivative_central_difference[ - dtype: DType, - func: fn[dtype: DType]( - x: Scalar[dtype], args: Optional[List[Scalar[dtype]]] - ) -> Scalar[dtype], + func: fn( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] + ) -> Scalar[f64], ]( - x0: Scalar[dtype], - args: Optional[List[Scalar[dtype]]], - tolerances: Dict[String, Scalar[dtype]] = {"atol": 1e-6, "rtol": 1e-6}, + x0: Scalar[f64], + args: Optional[List[Scalar[f64]]], + tolerances: Dict[String, Scalar[f64]] = {"atol": 1e-6, "rtol": 1e-6}, order: Int = 8, - initial_step: Scalar[dtype] = 0.5, - step_factor: Scalar[dtype] = 2.0, + initial_step: Scalar[f64] = 0.5, + step_factor: Scalar[f64] = 2.0, max_iter: Int = 10, -) raises -> DiffResult[dtype]: +) raises -> DiffResult[f64]: """Computes first derivative using central finite difference method. Uses symmetric stencils around the evaluation point with adaptive step size reduction until convergence within specified tolerances. Parameters: - dtype: The floating-point data type. - func: Function to differentiate with signature fn(x, args) -> Scalar[dtype]. + func: Function to differentiate with signature fn(x, args) -> Scalar[f64]. Args: x0: Point at which to evaluate the derivative. @@ -180,24 +170,24 @@ fn _derivative_central_difference[ step_factor: Factor by which to reduce step size in each iteration. max_iter: Maximum number of iterations. - Returns: - DiffResult[dtype] containing the derivative and convergence information. - Raises: Error: If the specified order is not in {2, 4, 6, 8}. Error: If tolerances, step size, step factor, or max_iter are invalid. + + Returns: + DiffResult[f64] containing the derivative and convergence information. """ comptime first_order_coefficients_compiletime: Dict[ - Int, List[Scalar[dtype]] - ] = generate_central_finite_difference_table[dtype]() + Int, List[Scalar[f64]] + ] = generate_central_finite_difference_table() var first_order_coefficients = materialize[ first_order_coefficients_compiletime ]() - var diff_estimate: Scalar[dtype] = 0.0 - var prev_diff: Scalar[dtype] = 0.0 - var atol: Scalar[dtype] = tolerances["atol"] - var rtol: Scalar[dtype] = tolerances["rtol"] + var diff_estimate: Scalar[f64] = 0.0 + var prev_diff: Scalar[f64] = 0.0 + var atol: Scalar[f64] = tolerances["atol"] + var rtol: Scalar[f64] = tolerances["rtol"] if atol < 0: raise Error( @@ -247,7 +237,7 @@ fn _derivative_central_difference[ " computation." ) - var coefficients: List[Scalar[dtype]] + var coefficients: List[Scalar[f64]] if order in (2, 4, 6, 8): coefficients = first_order_coefficients[order].copy() else: @@ -260,14 +250,14 @@ fn _derivative_central_difference[ " function evaluations.\n Available orders map to truncation" " errors: 2→O(h²), 4→O(h⁴), 6→O(h⁶), 8→O(h⁸)" ) - var step: Scalar[dtype] = initial_step + var step: Scalar[f64] = initial_step for i in range(max_iter): diff_estimate = 0.0 var j: Int = 0 for ref coeff in coefficients: diff_estimate += coeff * func( - x0 + step * Scalar[dtype](j - len(coefficients) // 2), args + x0 + step * Scalar[f64](j - len(coefficients) // 2), args ) j += 1 diff_estimate /= step @@ -275,7 +265,7 @@ fn _derivative_central_difference[ var diff_change = abs(diff_estimate - prev_diff) var tolerance_threshold = atol + rtol * abs(diff_estimate) if diff_change < tolerance_threshold: - return DiffResult[dtype]( + return DiffResult[f64]( success=True, df=diff_estimate, error=diff_change, @@ -287,7 +277,7 @@ fn _derivative_central_difference[ prev_diff = diff_estimate step /= step_factor - return DiffResult[dtype]( + return DiffResult[f64]( success=False, df=diff_estimate, error=0.0, @@ -298,19 +288,18 @@ fn _derivative_central_difference[ fn _derivative_forward_difference[ - dtype: DType, - func: fn[dtype: DType]( - x: Scalar[dtype], args: Optional[List[Scalar[dtype]]] - ) -> Scalar[dtype], + func: fn( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] + ) -> Scalar[f64], ]( - x0: Scalar[dtype], - args: Optional[List[Scalar[dtype]]], - tolerances: Dict[String, Scalar[dtype]] = {"atol": 1e-6, "rtol": 1e-6}, + x0: Scalar[f64], + args: Optional[List[Scalar[f64]]], + tolerances: Dict[String, Scalar[f64]] = {"atol": 1e-6, "rtol": 1e-6}, order: Int = 8, - initial_step: Scalar[dtype] = 0.5, - step_factor: Scalar[dtype] = 2.0, + initial_step: Scalar[f64] = 0.5, + step_factor: Scalar[f64] = 2.0, max_iter: Int = 10, -) raises -> DiffResult[dtype]: +) raises -> DiffResult[f64]: """Computes first derivative using forward finite difference method. Uses one-sided stencils in the forward direction with adaptive step size @@ -318,8 +307,7 @@ fn _derivative_forward_difference[ are not available. Parameters: - dtype: The floating-point data type. - func: Function to differentiate with signature fn(x, args) -> Scalar[dtype]. + func: Function to differentiate with signature fn(x, args) -> Scalar[f64]. Args: x0: Point at which to evaluate the derivative. @@ -331,23 +319,23 @@ fn _derivative_forward_difference[ max_iter: Maximum number of iterations. Returns: - DiffResult[dtype] containing the derivative and convergence information. + DiffResult[f64] containing the derivative and convergence information. Raises: Error: If the specified order is not in {1, 2, 3, 4, 5, 6}. Error: If tolerances, step size, step factor, or max_iter are invalid. """ comptime first_order_coefficients_compiletime: Dict[ - Int, List[Scalar[dtype]] - ] = generate_forward_finite_difference_table[dtype]() + Int, List[Scalar[f64]] + ] = generate_forward_finite_difference_table() var first_order_coefficients = materialize[ first_order_coefficients_compiletime ]() - var diff_estimate: Scalar[dtype] = 0.0 - var prev_diff: Scalar[dtype] = 0.0 - var atol: Scalar[dtype] = tolerances["atol"] - var rtol: Scalar[dtype] = tolerances["rtol"] + var diff_estimate: Scalar[f64] = 0.0 + var prev_diff: Scalar[f64] = 0.0 + var atol: Scalar[f64] = tolerances["atol"] + var rtol: Scalar[f64] = tolerances["rtol"] if atol < 0: raise Error( @@ -409,20 +397,20 @@ fn _derivative_forward_difference[ " numerical stability.\n Available orders map to truncation" " errors: 1→O(h), 2→O(h²), ..., 6→O(h⁶)" ) - var step: Scalar[dtype] = initial_step + var step: Scalar[f64] = initial_step for i in range(max_iter): diff_estimate = 0.0 var j: Int = 0 for ref coeff in coefficients: - diff_estimate += coeff * func(x0 + step * Scalar[dtype](j), args) + diff_estimate += coeff * func(x0 + step * Scalar[f64](j), args) j += 1 diff_estimate /= step if i > 0: var diff_change = abs(diff_estimate - prev_diff) var tolerance_threshold = atol + rtol * abs(diff_estimate) if diff_change < tolerance_threshold: - return DiffResult[dtype]( + return DiffResult[f64]( success=True, df=diff_estimate, error=diff_change, @@ -434,7 +422,7 @@ fn _derivative_forward_difference[ prev_diff = diff_estimate step /= step_factor - return DiffResult[dtype]( + return DiffResult[f64]( success=False, df=diff_estimate, error=0.0, @@ -445,19 +433,18 @@ fn _derivative_forward_difference[ fn _derivative_backward_difference[ - dtype: DType, - func: fn[dtype: DType]( - x: Scalar[dtype], args: Optional[List[Scalar[dtype]]] - ) -> Scalar[dtype], + func: fn( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] + ) -> Scalar[f64], ]( - x0: Scalar[dtype], - args: Optional[List[Scalar[dtype]]], - tolerances: Dict[String, Scalar[dtype]] = {"atol": 1e-6, "rtol": 1e-6}, + x0: Scalar[f64], + args: Optional[List[Scalar[f64]]], + tolerances: Dict[String, Scalar[f64]] = {"atol": 1e-6, "rtol": 1e-6}, order: Int = 8, - initial_step: Scalar[dtype] = 0.5, - step_factor: Scalar[dtype] = 2.0, + initial_step: Scalar[f64] = 0.5, + step_factor: Scalar[f64] = 2.0, max_iter: Int = 10, -) raises -> DiffResult[dtype]: +) raises -> DiffResult[f64]: """Computes first derivative using backward finite difference method. Uses one-sided stencils in the backward direction with adaptive step size @@ -465,8 +452,7 @@ fn _derivative_backward_difference[ are not available. Parameters: - dtype: The floating-point data type. - func: Function to differentiate with signature fn(x, args) -> Scalar[dtype]. + func: Function to differentiate with signature fn(x, args) -> Scalar[f64]. Args: x0: Point at which to evaluate the derivative. @@ -478,23 +464,23 @@ fn _derivative_backward_difference[ max_iter: Maximum number of iterations. Returns: - DiffResult[dtype] containing the derivative and convergence information. + DiffResult[f64] containing the derivative and convergence information. Raises: Error: If the specified order is not in {1, 2, 3, 4, 5, 6}. Error: If tolerances, step size, step factor, or max_iter are invalid. """ comptime first_order_coefficients_compiletime: Dict[ - Int, List[Scalar[dtype]] - ] = generate_backward_finite_difference_table[dtype]() + Int, List[Scalar[f64]] + ] = generate_backward_finite_difference_table() var first_order_coefficients = materialize[ first_order_coefficients_compiletime ]() - var diff_estimate: Scalar[dtype] = 0.0 - var prev_diff: Scalar[dtype] = 0.0 - var atol: Scalar[dtype] = tolerances["atol"] - var rtol: Scalar[dtype] = tolerances["rtol"] + var diff_estimate: Scalar[f64] = 0.0 + var prev_diff: Scalar[f64] = 0.0 + var atol: Scalar[f64] = tolerances["atol"] + var rtol: Scalar[f64] = tolerances["rtol"] if atol < 0: raise Error( @@ -556,11 +542,11 @@ fn _derivative_backward_difference[ " numerical stability.\n Available orders map to truncation" " errors: 1→O(h), 2→O(h²), ..., 6→O(h⁶)" ) - var step: Scalar[dtype] = initial_step + var step: Scalar[f64] = initial_step for i in range(max_iter): diff_estimate = 0.0 - var j: Scalar[dtype] = 0 + var j: Scalar[f64] = 0 for ref coeff in coefficients: diff_estimate += coeff * func(x0 + step * j, args) j += 1 @@ -569,7 +555,7 @@ fn _derivative_backward_difference[ var diff_change = abs(diff_estimate - prev_diff) var tolerance_threshold = atol + rtol * abs(diff_estimate) if diff_change < tolerance_threshold: - return DiffResult[dtype]( + return DiffResult[f64]( success=True, df=diff_estimate, error=diff_change, @@ -581,7 +567,7 @@ fn _derivative_backward_difference[ prev_diff = diff_estimate step /= step_factor - return DiffResult[dtype]( + return DiffResult[f64]( success=False, df=diff_estimate, error=0.0, diff --git a/scijo/differentiate/jacob.mojo b/scijo/differentiate/jacob.mojo index c080e99..1e26d0c 100644 --- a/scijo/differentiate/jacob.mojo +++ b/scijo/differentiate/jacob.mojo @@ -18,16 +18,15 @@ from algorithm.functional import parallelize fn jacobian[ - dtype: DType, - f: fn[dtype: DType]( - x: NDArray[dtype], args: Optional[List[Scalar[dtype]]] - ) raises -> NDArray[dtype], + f: fn( + x: NDArray[f64], args: Optional[List[Scalar[f64]]] + ) raises -> NDArray[f64], ]( - x: NDArray[dtype], - args: Optional[List[Scalar[dtype]]] = None, - tolerances: Dict[String, Scalar[dtype]] = {"abs": 1e-5, "rel": 1e-3}, + x: NDArray[f64], + args: Optional[List[Scalar[f64]]] = None, + tolerances: Dict[String, Scalar[f64]] = {"abs": 1e-5, "rel": 1e-3}, maxiter: Int = 10, -) raises -> NDArray[dtype]: +) raises -> NDArray[f64]: """Computes the Jacobian matrix of a vector-valued function using central finite differences. Evaluates J[i, j] = ∂f_i/∂x_j using the central difference formula @@ -35,8 +34,7 @@ fn jacobian[ Each column of the Jacobian is computed in parallel. Parameters: - dtype: The floating-point data type. - f: Vector-valued function with signature fn(x, args) -> NDArray[dtype]. + f: Vector-valued function with signature fn(x, args) -> NDArray[f64]. Args: x: Input vector of shape (n,) at which to evaluate the Jacobian. @@ -48,7 +46,7 @@ fn jacobian[ Error: If function evaluation fails for any perturbation. Returns: - NDArray[dtype] of shape (m, n) representing the Jacobian matrix, + NDArray[f64] of shape (m, n) representing the Jacobian matrix, where m is the output dimension and n is the input dimension. Examples: @@ -57,19 +55,19 @@ fn jacobian[ from scijo.differentiate import jacobian from scijo.prelude import * - fn f[dtype: DType](x: NDArray[dtype], args: Optional[List[Scalar[dtype]]]) raises -> NDArray[dtype]: + fn f(x: NDArray[f64], args: Optional[List[Scalar[f64]]]) raises -> NDArray[f64]: return x * x var x = nm.array[f64]([1.0, 2.0]) - var J = jacobian[f64, f](x) + var J = jacobian[f](x) ``` """ var n: Int = len(x) - var f0: NDArray[dtype] = f(x, args) + var f0: NDArray[f64] = f(x, args) var m: Int = len(f0) - var jacob: NDArray[dtype] = zeros[dtype](Shape(m, n)) - var step: NDArray[dtype] = full[dtype](Shape(n), fill_value=0.5) + var jacob: NDArray[f64] = zeros[f64](Shape(m, n)) + var step: NDArray[f64] = full[f64](Shape(n), fill_value=0.5) @parameter fn closure(j: Int): @@ -82,10 +80,10 @@ fn jacobian[ var f_plus = f(x_plus, args) var f_minus = f(x_minus, args) - var col: NDArray[dtype] = (f_plus - f_minus) / (2.0 * hj) + var col: NDArray[f64] = (f_plus - f_minus) / (2.0 * hj) for i in range(m): - var val: Scalar[dtype] = col.load(i) + var val: Scalar[f64] = col.load(i) var flat_idx: Int = i * n + j jacob.store(flat_idx, val=val) except: diff --git a/scijo/differentiate/utility.mojo b/scijo/differentiate/utility.mojo index 9c2d362..dc48956 100644 --- a/scijo/differentiate/utility.mojo +++ b/scijo/differentiate/utility.mojo @@ -23,7 +23,7 @@ References: # ===----------------------------------------------------------------------=== # -struct DiffResult[dtype: DType](ImplicitlyCopyable, Writable): +struct DiffResult[dtype: DType = DType.float64](ImplicitlyCopyable, Writable): """Result structure for numerical differentiation operations. Encapsulates the results of derivative computations, including the computed @@ -98,28 +98,23 @@ struct DiffResult[dtype: DType](ImplicitlyCopyable, Writable): @parameter -fn generate_central_finite_difference_table[ - dtype: DType -]() -> Dict[Int, List[Scalar[dtype]]]: +fn generate_central_finite_difference_table() -> Dict[Int, List[Scalar[f64]]]: """Generates central finite difference coefficients for first-order derivatives. Creates a lookup table of coefficients for central difference approximations using symmetric stencils: f'(x) ≈ Σ(c_i * f(x + i*h)) / h. - Parameters: - dtype: The floating-point data type for the coefficients. - Returns: Coefficient arrays indexed by accuracy order. Available orders: 2, 4, 6, 8 with truncation errors O(h²), O(h⁴), O(h⁶), O(h⁸). """ - var coefficients = Dict[Int, List[Scalar[dtype]]]() + var coefficients = Dict[Int, List[Scalar[f64]]]() # Order 2: points [-1, 0, 1] - coefficients[2]: List[Scalar[dtype]] = [-0.5, 0.0, 0.5] + coefficients[2]: List[Scalar[f64]] = [-0.5, 0.0, 0.5] # Order 4: points [-2, -1, 0, 1, 2] - coefficients[4]: List[Scalar[dtype]] = [ + coefficients[4]: List[Scalar[f64]] = [ 1.0 / 12.0, -2.0 / 3.0, 0.0, @@ -128,7 +123,7 @@ fn generate_central_finite_difference_table[ ] # Order 6: points [-3, -2, -1, 0, 1, 2, 3] - coefficients[6]: List[Scalar[dtype]] = [ + coefficients[6]: List[Scalar[f64]] = [ -1.0 / 60.0, 3.0 / 20.0, -3.0 / 4.0, @@ -139,7 +134,7 @@ fn generate_central_finite_difference_table[ ] # Order 8: points [-4, -3, -2, -1, 0, 1, 2, 3, 4] - coefficients[8]: List[Scalar[dtype]] = [ + coefficients[8]: List[Scalar[f64]] = [ 1.0 / 280.0, -4.0 / 105.0, 1.0 / 5.0, @@ -155,9 +150,7 @@ fn generate_central_finite_difference_table[ @parameter -fn generate_forward_finite_difference_table[ - dtype: DType -]() -> Dict[Int, List[Scalar[dtype]]]: +fn generate_forward_finite_difference_table() -> Dict[Int, List[Scalar[f64]]]: """Generates forward finite difference coefficients for first-order derivatives. Creates a lookup table of coefficients for forward difference approximations @@ -166,23 +159,20 @@ fn generate_forward_finite_difference_table[ Forward differences are necessary at left domain boundaries or when backward evaluations are not feasible. - Parameters: - dtype: The floating-point data type for the coefficients. - Returns: Coefficient arrays indexed by accuracy order. Available orders: 1, 2, 3, 4, 5, 6 with truncation errors O(h) through O(h⁶). """ - var coefficients = Dict[Int, List[Scalar[dtype]]]() + var coefficients = Dict[Int, List[Scalar[f64]]]() # Order 1: points [0, 1] - coefficients[1]: List[Scalar[dtype]] = [-1.0, 1.0] + coefficients[1]: List[Scalar[f64]] = [-1.0, 1.0] # Order 2: points [0, 1, 2] - coefficients[2]: List[Scalar[dtype]] = [-3.0 / 2.0, 2.0, -1.0 / 2.0] + coefficients[2]: List[Scalar[f64]] = [-3.0 / 2.0, 2.0, -1.0 / 2.0] # Order 3: points [0, 1, 2, 3] - coefficients[3]: List[Scalar[dtype]] = [ + coefficients[3]: List[Scalar[f64]] = [ -11.0 / 6.0, 3.0, -3.0 / 2.0, @@ -190,7 +180,7 @@ fn generate_forward_finite_difference_table[ ] # Order 4: points [0, 1, 2, 3, 4] - coefficients[4]: List[Scalar[dtype]] = [ + coefficients[4]: List[Scalar[f64]] = [ -25.0 / 12.0, 4.0, -3.0, @@ -199,7 +189,7 @@ fn generate_forward_finite_difference_table[ ] # Order 5: points [0, 1, 2, 3, 4, 5] - coefficients[5]: List[Scalar[dtype]] = [ + coefficients[5]: List[Scalar[f64]] = [ -137.0 / 60.0, 5.0, -5.0, @@ -209,7 +199,7 @@ fn generate_forward_finite_difference_table[ ] # Order 6: points [0, 1, 2, 3, 4, 5, 6] - coefficients[6]: List[Scalar[dtype]] = [ + coefficients[6]: List[Scalar[f64]] = [ -49.0 / 20.0, 6.0, -15.0 / 2.0, @@ -223,9 +213,7 @@ fn generate_forward_finite_difference_table[ @parameter -fn generate_backward_finite_difference_table[ - dtype: DType -]() -> Dict[Int, List[Scalar[dtype]]]: +fn generate_backward_finite_difference_table() -> Dict[Int, List[Scalar[f64]]]: """Generates backward finite difference coefficients for first-order derivatives. Creates a lookup table of coefficients for backward difference approximations @@ -235,23 +223,20 @@ fn generate_backward_finite_difference_table[ stencil and adjusting signs. They are necessary at right domain boundaries or when forward evaluations are not feasible. - Parameters: - dtype: The floating-point data type for the coefficients. - Returns: Coefficient arrays indexed by accuracy order. Available orders: 1, 2, 3, 4, 5, 6 with truncation errors O(h) through O(h⁶). """ - var coefficients = Dict[Int, List[Scalar[dtype]]]() + var coefficients = Dict[Int, List[Scalar[f64]]]() # Order 1: points [-1, 0] - coefficients[1]: List[Scalar[dtype]] = [-1.0, 1.0] + coefficients[1]: List[Scalar[f64]] = [-1.0, 1.0] # Order 2: points [-2, -1, 0] - coefficients[2]: List[Scalar[dtype]] = [1.0 / 2.0, -2.0, 3.0 / 2.0] + coefficients[2]: List[Scalar[f64]] = [1.0 / 2.0, -2.0, 3.0 / 2.0] # Order 3: points [-3, -2, -1, 0] - coefficients[3]: List[Scalar[dtype]] = [ + coefficients[3]: List[Scalar[f64]] = [ -1.0 / 3.0, 3.0 / 2.0, -3.0, @@ -259,7 +244,7 @@ fn generate_backward_finite_difference_table[ ] # Order 4: points [-4, -3, -2, -1, 0] - coefficients[4]: List[Scalar[dtype]] = [ + coefficients[4]: List[Scalar[f64]] = [ 1.0 / 4.0, -4.0 / 3.0, 3.0, @@ -268,7 +253,7 @@ fn generate_backward_finite_difference_table[ ] # Order 5: points [-5, -4, -3, -2, -1, 0] - coefficients[5]: List[Scalar[dtype]] = [ + coefficients[5]: List[Scalar[f64]] = [ -1.0 / 5.0, 5.0 / 4.0, -10.0 / 3.0, @@ -278,7 +263,7 @@ fn generate_backward_finite_difference_table[ ] # Order 6: points [-6, -5, -4, -3, -2, -1, 0] - coefficients[6]: List[Scalar[dtype]] = [ + coefficients[6]: List[Scalar[f64]] = [ 1.0 / 6.0, -6.0 / 5.0, 15.0 / 4.0, From ae6a5a77f2523391b79dcbfcb6352f3c66bc7fe6 Mon Sep 17 00:00:00 2001 From: shivasankar Date: Fri, 13 Mar 2026 22:41:45 +0900 Subject: [PATCH 2/4] format differentiate module files --- scijo/differentiate/deriv.mojo | 16 ++++------------ scijo/differentiate/jacob.mojo | 6 +++--- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/scijo/differentiate/deriv.mojo b/scijo/differentiate/deriv.mojo index 8be2ad4..1002677 100644 --- a/scijo/differentiate/deriv.mojo +++ b/scijo/differentiate/deriv.mojo @@ -32,9 +32,7 @@ from .utility import ( fn derivative[ - func: fn( - x: Scalar[f64], args: Optional[List[Scalar[f64]]] - ) -> Scalar[f64], + func: fn(x: Scalar[f64], args: Optional[List[Scalar[f64]]]) -> Scalar[f64], *, step_direction: Int = 0, ]( @@ -141,9 +139,7 @@ fn derivative[ fn _derivative_central_difference[ - func: fn( - x: Scalar[f64], args: Optional[List[Scalar[f64]]] - ) -> Scalar[f64], + func: fn(x: Scalar[f64], args: Optional[List[Scalar[f64]]]) -> Scalar[f64], ]( x0: Scalar[f64], args: Optional[List[Scalar[f64]]], @@ -288,9 +284,7 @@ fn _derivative_central_difference[ fn _derivative_forward_difference[ - func: fn( - x: Scalar[f64], args: Optional[List[Scalar[f64]]] - ) -> Scalar[f64], + func: fn(x: Scalar[f64], args: Optional[List[Scalar[f64]]]) -> Scalar[f64], ]( x0: Scalar[f64], args: Optional[List[Scalar[f64]]], @@ -433,9 +427,7 @@ fn _derivative_forward_difference[ fn _derivative_backward_difference[ - func: fn( - x: Scalar[f64], args: Optional[List[Scalar[f64]]] - ) -> Scalar[f64], + func: fn(x: Scalar[f64], args: Optional[List[Scalar[f64]]]) -> Scalar[f64], ]( x0: Scalar[f64], args: Optional[List[Scalar[f64]]], diff --git a/scijo/differentiate/jacob.mojo b/scijo/differentiate/jacob.mojo index 1e26d0c..7cdb81a 100644 --- a/scijo/differentiate/jacob.mojo +++ b/scijo/differentiate/jacob.mojo @@ -18,9 +18,9 @@ from algorithm.functional import parallelize fn jacobian[ - f: fn( - x: NDArray[f64], args: Optional[List[Scalar[f64]]] - ) raises -> NDArray[f64], + f: fn(x: NDArray[f64], args: Optional[List[Scalar[f64]]]) raises -> NDArray[ + f64 + ], ]( x: NDArray[f64], args: Optional[List[Scalar[f64]]] = None, From 33e200f339c6702892341ef478a31bb9c840f518 Mon Sep 17 00:00:00 2001 From: shivasankar Date: Fri, 13 Mar 2026 22:45:57 +0900 Subject: [PATCH 3/4] Update test_differentiate.mojo --- tests/test_differentiate.mojo | 147 ++++++++++++++++------------------ 1 file changed, 71 insertions(+), 76 deletions(-) diff --git a/tests/test_differentiate.mojo b/tests/test_differentiate.mojo index c9b5d7a..daee50a 100644 --- a/tests/test_differentiate.mojo +++ b/tests/test_differentiate.mojo @@ -1,82 +1,77 @@ from scijo.differentiate import derivative +from scijo.prelude import * from testing import assert_almost_equal, assert_equal, assert_true, assert_false from testing import TestSuite import math -fn constant_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype]: +fn constant_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = 5, f'(x) = 0. """ return 5.0 -fn linear_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype]: +fn linear_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = 3x + 2, f'(x) = 3. """ return 3.0 * x + 2.0 -fn quadratic_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype]: +fn quadratic_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = 2x^2 + 3x + 1, f'(x) = 4x + 3. """ return 2.0 * x * x + 3.0 * x + 1.0 -fn cubic_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype]: +fn cubic_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = x^3 - 2x^2 + x - 5, f'(x) = 3x^2 - 4x + 1. """ return x * x * x - 2.0 * x * x + x - 5.0 -fn sin_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[ - dtype -] where dtype.is_floating_point(): +fn sin_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = sin(x), f'(x) = cos(x). """ return math.sin(x) -fn cos_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[ - dtype -] where dtype.is_floating_point(): +fn cos_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = cos(x), f'(x) = -sin(x). """ return math.cos(x) -fn exp_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[ - dtype -] where dtype.is_floating_point(): +fn exp_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = e^x, f'(x) = e^x. """ return math.exp(x) -fn parameterized_function[ - dtype: DType -](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype]: +fn parameterized_function( + x: Scalar[f64], args: Optional[List[Scalar[f64]]] +) -> Scalar[f64]: """ F(x) = a*x^2 + b*x + c, f'(x) = 2*a*x + b. """ @@ -90,9 +85,9 @@ fn test_basic_derivatives() raises: """Test derivatives of basic polynomial functions.""" # Test constant function: F(x) = 5, f'(x) = 0 - var result_const = derivative[ - DType.float64, constant_function, step_direction=0 - ](x0=2.0, args=None) + var result_const = derivative[constant_function, step_direction=0]( + x0=2.0, args=None + ) assert_true( result_const.success, "Constant function derivative should converge" ) @@ -104,9 +99,9 @@ fn test_basic_derivatives() raises: ) # Test linear function: F(x) = 3x + 2, f'(x) = 3 - var result_linear = derivative[ - DType.float64, linear_function, step_direction=0 - ](x0=1.5, args=None) + var result_linear = derivative[linear_function, step_direction=0]( + x0=1.5, args=None + ) assert_true( result_linear.success, "Linear function derivative should converge" ) @@ -117,9 +112,9 @@ fn test_basic_derivatives() raises: # Test quadratic function: F(x) = 2x^2 + 3x + 1, f'(x) = 4x + 3 var x_quad = 2.0 var expected_quad = 4.0 * x_quad + 3.0 - var result_quad = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=x_quad, args=None) + var result_quad = derivative[quadratic_function, step_direction=0]( + x0=x_quad, args=None + ) assert_true( result_quad.success, "Quadratic function derivative should converge" ) @@ -140,9 +135,9 @@ fn test_cubic_derivatives() raises: for i in range(len(test_points)): var x = test_points[i] var expected = 3.0 * x * x - 4.0 * x + 1.0 - var result = derivative[ - DType.float64, cubic_function, step_direction=0 - ](x0=x, args=None) + var result = derivative[cubic_function, step_direction=0]( + x0=x, args=None + ) assert_true(result.success, "Cubic function derivative should converge") assert_almost_equal( result.df, @@ -158,7 +153,7 @@ fn test_trigonometric_derivatives() raises: # Test sin(x): f'(x) = cos(x) var x_sin = 0.5 var expected_sin = math.cos(x_sin) - var result_sin = derivative[DType.float64, sin_function, step_direction=0]( + var result_sin = derivative[sin_function, step_direction=0]( x0=x_sin, args=None, order=6 ) assert_true(result_sin.success, "Sin function derivative should converge") @@ -172,7 +167,7 @@ fn test_trigonometric_derivatives() raises: # Test cos(x): f'(x) = -sin(x) var x_cos = 1.0 var expected_cos = -math.sin(x_cos) - var result_cos = derivative[DType.float64, cos_function, step_direction=0]( + var result_cos = derivative[cos_function, step_direction=0]( x0=x_cos, args=None, order=6 ) assert_true(result_cos.success, "Cos function derivative should converge") @@ -190,7 +185,7 @@ fn test_exponential_derivative() raises: # Test exp(x): f'(x) = exp(x) var x_exp = 1.0 var expected_exp = math.exp(x_exp) - var result_exp = derivative[DType.float64, exp_function, step_direction=0]( + var result_exp = derivative[exp_function, step_direction=0]( x0=x_exp, args=None, order=6 ) assert_true(result_exp.success, "Exp function derivative should converge") @@ -210,9 +205,9 @@ fn test_parameterized_function() raises: var args: List[Scalar[DType.float64]] = [2.0, 5.0, 3.0] var x_param = 1.5 var expected_param = 4.0 * x_param + 5.0 # = 11.0 - var result_param = derivative[ - DType.float64, parameterized_function, step_direction=0 - ](x0=x_param, args=args^) + var result_param = derivative[parameterized_function, step_direction=0]( + x0=x_param, args=args^ + ) assert_true( result_param.success, "Parameterized function derivative should converge", @@ -232,9 +227,9 @@ fn test_different_step_directions() raises: var expected = 4.0 * x_test + 3.0 # Central differences - var result_central = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=x_test, args=None) + var result_central = derivative[quadratic_function, step_direction=0]( + x0=x_test, args=None + ) assert_true(result_central.success, "Central difference should converge") assert_almost_equal( result_central.df, @@ -244,9 +239,9 @@ fn test_different_step_directions() raises: ) # Forward differences - var result_forward = derivative[ - DType.float64, quadratic_function, step_direction=1 - ](x0=x_test, args=None, order=6, max_iter=50) + var result_forward = derivative[quadratic_function, step_direction=1]( + x0=x_test, args=None, order=6, max_iter=50 + ) assert_true(result_forward.success, "Forward difference should converge") assert_almost_equal( result_forward.df, @@ -256,9 +251,9 @@ fn test_different_step_directions() raises: ) # Backward differences - var result_backward = derivative[ - DType.float64, quadratic_function, step_direction= -1 - ](x0=x_test, args=None, order=6, max_iter=50) + var result_backward = derivative[quadratic_function, step_direction= -1]( + x0=x_test, args=None, order=6, max_iter=50 + ) assert_true(result_backward.success, "Backward difference should converge") assert_almost_equal( result_backward.df, @@ -278,9 +273,9 @@ fn test_different_orders() raises: for i in range(len(orders)): var order = orders[i] - var result = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=x_test, args=None, order=order) + var result = derivative[quadratic_function, step_direction=0]( + x0=x_test, args=None, order=order + ) assert_true(result.success, "Higher order should converge") assert_almost_equal( result.df, @@ -300,9 +295,9 @@ fn test_tolerance_settings() raises: strict_tolerance["atol"] = 1e-10 strict_tolerance["rtol"] = 1e-10 - var result_strict = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=x_test, args=None, tolerances=strict_tolerance) + var result_strict = derivative[quadratic_function, step_direction=0]( + x0=x_test, args=None, tolerances=strict_tolerance + ) assert_almost_equal( result_strict.df, expected, @@ -314,9 +309,9 @@ fn test_tolerance_settings() raises: loose_tolerance["atol"] = 1e-3 loose_tolerance["rtol"] = 1e-3 - var result_loose = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=x_test, args=None, tolerances=loose_tolerance) + var result_loose = derivative[quadratic_function, step_direction=0]( + x0=x_test, args=None, tolerances=loose_tolerance + ) assert_almost_equal( result_loose.df, expected, @@ -328,9 +323,9 @@ fn test_tolerance_settings() raises: fn test_convergence_properties() raises: """Test convergence properties and diagnostic information.""" - var result = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=1.0, args=None, max_iter=5) + var result = derivative[quadratic_function, step_direction=0]( + x0=1.0, args=None, max_iter=5 + ) assert_true(result.nit > 0, "Number of iterations should be positive") assert_true( @@ -346,7 +341,7 @@ fn test_error_conditions() raises: """Test error conditions and invalid parameters.""" try: - var _ = derivative[DType.float64, quadratic_function, step_direction=5]( + var _ = derivative[quadratic_function, step_direction=5]( x0=1.0, args=None ) assert_false(True, "Invalid step direction should raise an error") @@ -364,9 +359,9 @@ fn test_step_size_parameters() raises: for i in range(len(step_sizes)): var step = step_sizes[i] - var result = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=x_test, args=None, initial_step=step) + var result = derivative[quadratic_function, step_direction=0]( + x0=x_test, args=None, initial_step=step + ) assert_true(result.success, "Different initial step sizes should work") assert_almost_equal( result.df, @@ -379,9 +374,9 @@ fn test_step_size_parameters() raises: for i in range(len(factors)): var factor = factors[i] - var result = derivative[ - DType.float64, quadratic_function, step_direction=0 - ](x0=x_test, args=None, step_factor=factor) + var result = derivative[quadratic_function, step_direction=0]( + x0=x_test, args=None, step_factor=factor + ) assert_true(result.success, "Different step factors should work") assert_almost_equal( result.df, From cd5db6e7daa7d03640875948988542394a8c1ce7 Mon Sep 17 00:00:00 2001 From: shivasankar Date: Fri, 13 Mar 2026 22:47:08 +0900 Subject: [PATCH 4/4] update format. --- scijo/integrate/fixed_sample.mojo | 2 ++ scijo/integrate/quadrature.mojo | 1 + scijo/integrate/utility.mojo | 4 ++++ scijo/optimize/root_scalar.mojo | 34 +++++++++++++++++++++++-------- scijo/optimize/utility.mojo | 1 + 5 files changed, 33 insertions(+), 9 deletions(-) diff --git a/scijo/integrate/fixed_sample.mojo b/scijo/integrate/fixed_sample.mojo index 133bf50..00392ce 100644 --- a/scijo/integrate/fixed_sample.mojo +++ b/scijo/integrate/fixed_sample.mojo @@ -193,6 +193,7 @@ fn trapezoid[ # Simpson # ===----------------------------------------------------------------------=== # + fn simpson[ dtype: DType ]( @@ -359,6 +360,7 @@ fn simpson[ # Romberg # ===----------------------------------------------------------------------=== # + # TODO: fix the loop implementation. fn romb[ dtype: DType diff --git a/scijo/integrate/quadrature.mojo b/scijo/integrate/quadrature.mojo index 04a0324..203df46 100644 --- a/scijo/integrate/quadrature.mojo +++ b/scijo/integrate/quadrature.mojo @@ -52,6 +52,7 @@ from .utility import ( # Quad # ===----------------------------------------------------------------------=== # + fn quad[ dtype: DType, func: fn[dtype: DType]( diff --git a/scijo/integrate/utility.mojo b/scijo/integrate/utility.mojo index 11d412d..eeb674d 100644 --- a/scijo/integrate/utility.mojo +++ b/scijo/integrate/utility.mojo @@ -31,6 +31,7 @@ comptime largest_positive_dtype[dtype: DType] = max_finite[dtype]() # Helpers # ===----------------------------------------------------------------------=== # + fn machine_epsilon[dtype: DType]() -> Float64 where dtype.is_floating_point(): """Returns the machine epsilon for the given floating-point dtype. @@ -40,6 +41,7 @@ fn machine_epsilon[dtype: DType]() -> Float64 where dtype.is_floating_point(): Returns: The machine epsilon as a Float64 value. """ + # TODO: Check if these values are correct lol @parameter if dtype == DType.float16: @@ -54,6 +56,7 @@ fn machine_epsilon[dtype: DType]() -> Float64 where dtype.is_floating_point(): # Adaptive intervals # ===----------------------------------------------------------------------=== # + struct QAGSInterval[dtype: DType](ImplicitlyCopyable, Movable): """Represents an integration subinterval with error estimate for adaptive subdivision. @@ -243,6 +246,7 @@ fn get_quad_error_message(ier: Int) -> String: # Result types # ===----------------------------------------------------------------------=== # + struct IntegralResult[dtype: DType](Copyable, Movable, Writable): """Result structure for numerical integration operations. diff --git a/scijo/optimize/root_scalar.mojo b/scijo/optimize/root_scalar.mojo index 74ce5a3..eda4365 100644 --- a/scijo/optimize/root_scalar.mojo +++ b/scijo/optimize/root_scalar.mojo @@ -18,6 +18,7 @@ methods (secant). # Root scalar # ===----------------------------------------------------------------------=== # + fn root_scalar[ dtype: DType, f: fn[dtype: DType]( @@ -69,29 +70,37 @@ 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]( @@ -136,7 +145,10 @@ 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) @@ -144,7 +156,8 @@ 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 @@ -211,8 +224,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): @@ -281,7 +294,10 @@ 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 8cd7eb3..48d043a 100644 --- a/scijo/optimize/utility.mojo +++ b/scijo/optimize/utility.mojo @@ -14,6 +14,7 @@ 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.