Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
176 changes: 77 additions & 99 deletions scijo/differentiate/deriv.mojo

Large diffs are not rendered by default.

34 changes: 16 additions & 18 deletions scijo/differentiate/jacob.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -18,25 +18,23 @@ 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
(f(x + h*e_j) - f(x - h*e_j)) / (2h), where e_j is the j-th unit vector.
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.
Expand All @@ -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:
Expand All @@ -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):
Expand All @@ -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:
Expand Down
61 changes: 23 additions & 38 deletions scijo/differentiate/utility.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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
Expand All @@ -166,31 +159,28 @@ 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,
1.0 / 3.0,
]

# 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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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
Expand All @@ -235,31 +223,28 @@ 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,
11.0 / 6.0,
]

# 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,
Expand All @@ -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,
Expand All @@ -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,
Expand Down
2 changes: 2 additions & 0 deletions scijo/integrate/fixed_sample.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -193,6 +193,7 @@ fn trapezoid[
# Simpson
# ===----------------------------------------------------------------------=== #


fn simpson[
dtype: DType
](
Expand Down Expand Up @@ -359,6 +360,7 @@ fn simpson[
# Romberg
# ===----------------------------------------------------------------------=== #


# TODO: fix the loop implementation.
fn romb[
dtype: DType
Expand Down
1 change: 1 addition & 0 deletions scijo/integrate/quadrature.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ from .utility import (
# Quad
# ===----------------------------------------------------------------------=== #


fn quad[
dtype: DType,
func: fn[dtype: DType](
Expand Down
4 changes: 4 additions & 0 deletions scijo/integrate/utility.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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:
Expand All @@ -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.

Expand Down Expand Up @@ -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.

Expand Down
Loading
Loading