Skip to content
Merged
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
4 changes: 3 additions & 1 deletion benchmark/bench_derivative.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,9 @@ fn poly_func[

fn sin_func[
dtype: DType
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype] where dtype.is_floating_point():
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[
dtype
] where dtype.is_floating_point():
return sin(x)


Expand Down
6 changes: 3 additions & 3 deletions pixi.toml
Original file line number Diff line number Diff line change
Expand Up @@ -42,12 +42,12 @@ backend = { name = "pixi-build-mojo", version = "0.*"}
name = "scijo"

[package.host-dependencies]
modular = ">=26.2.0.dev2026022717,<27"
modular = ">=26.2.0.dev2026030605,<27"

[package.build-dependencies]
modular = ">=26.2.0.dev2026022717,<27"
modular = ">=26.2.0.dev2026030605,<27"
numojo = { git = "https://github.com/shivasankarka/NuMojo.git", branch = "pre-0.9"}

[package.run-dependencies]
modular = ">=26.2.0.dev2026022717,<27"
modular = ">=26.2.0.dev2026030605,<27"
numojo = { git = "https://github.com/shivasankarka/NuMojo.git", branch = "pre-0.9"}
5 changes: 1 addition & 4 deletions scijo/fft/fastfourier.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -144,10 +144,7 @@ fn _ifft_unnormalized[

for k in range(half_size):
var angle = (
2.0
* Constants.pi
* Scalar[dtype.dtype](k)
/ Scalar[dtype.dtype](n)
2.0 * Constants.pi * Scalar[dtype.dtype](k) / Scalar[dtype.dtype](n)
)
var twiddle = ComplexSIMD[dtype](
cos(angle).cast[dtype.dtype](), sin(angle).cast[dtype.dtype]()
Expand Down
21 changes: 4 additions & 17 deletions scijo/integrate/fixed_sample.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ Includes the composite trapezoidal rule, Simpson's rule, and Romberg integration
"""

from numojo.core.ndarray import NDArray, NDArrayShape
from numojo.core.error import NumojoError
import numojo as nm


Expand Down Expand Up @@ -48,26 +47,14 @@ fn trapezoid[

if y.ndim != 1:
raise Error(
NumojoError(
category="shape",
message=String(
"Expected y to be 1-D, received ndim={}. Pass a 1-D NDArray"
" for y (e.g. shape (N,))."
).format(y.ndim),
location="trapezoid(y, dx=1.0)",
)
t"Scijo [trapezoid]: Expected y to be 1-D array, received"
t" ndim={{y.ndim}}."
)

if y.size == 0:
raise Error(
NumojoError(
category="value",
message=(
"Cannot integrate over an empty array. Provide a non-empty"
" array for y."
),
location="trapezoid(y, dx=1.0)",
)
t"Scijo [trapezoid]: y.size = 0, Cannot interage over an empty"
t" array."
)

if y.size == 1:
Expand Down
89 changes: 61 additions & 28 deletions scijo/optimize/root_scalar.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,10 @@ methods (secant).

# TODO: check if we are using the right tolerance conditions in all methods.

# ===----------------------------------------------------------------------=== #
# Root scalar
# ===----------------------------------------------------------------------=== #


fn root_scalar[
dtype: DType,
Expand All @@ -27,6 +31,7 @@ fn root_scalar[
dtype
]
] = None,
*,
method: String = "bisect",
](
args: Optional[List[Scalar[dtype]]] = None,
Expand All @@ -36,7 +41,6 @@ fn root_scalar[
xtol: Scalar[dtype] = 1e-8,
rtol: Scalar[dtype] = 1e-8,
maxiter: Int = 100,
# options: SolverOptions
) raises -> Scalar[dtype]:
"""Finds a root of a scalar function using the specified method.

Expand All @@ -55,30 +59,46 @@ fn root_scalar[
rtol: Relative tolerance for convergence.
maxiter: Maximum number of iterations.

Returns:
The approximate root as a Scalar[dtype].

Raises:
Error: If required inputs for the chosen method are missing or invalid.

Returns:
The approximate root of the function.
"""

@parameter
if method == "newton" and fprime:
if method == "newton":
if not fprime:
raise Error(
"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("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(
"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("Unsupported method: " + String(method))
raise Error(
"Scijo [root_scalar]: Unsupported method: " + String(method)
)


# ===----------------------------------------------------------------------=== #
# Root scalar methods
# ===----------------------------------------------------------------------=== #


fn newton[
Expand Down Expand Up @@ -115,25 +135,29 @@ fn newton[
rtol: Relative tolerance for convergence.
maxiter: Maximum number of iterations.

Returns:
The approximate root as a Scalar[dtype].

Raises:
Error: If x0 is not provided or the derivative is zero at any step.

Returns:
The approximate root as a Scalar[dtype].
"""
var xn: Scalar[dtype]
if x0:
xn = x0.value()
else:
raise Error("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)
var fpx = fprime(xn, args)

if fpx == 0:
raise Error(
"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
Expand Down Expand Up @@ -181,11 +205,11 @@ fn bisect[
rtol: Relative tolerance for convergence.
maxiter: Maximum number of iterations.

Returns:
The approximate root as a Scalar[dtype].

Raises:
Error: If f(a) and f(b) do not have opposite signs.

Returns:
The approximate root as a Scalar[dtype].
"""
var a: Scalar[dtype] = bracket[0]
var b: Scalar[dtype] = bracket[1]
Expand All @@ -200,8 +224,8 @@ fn bisect[

if fa * fb > 0:
raise Error(
"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):
Expand All @@ -212,16 +236,15 @@ fn bisect[
return c

var tol_x = max(xtol, rtol * abs(c))
var half_width = (b - a) / 2
var half_width = abs(b - a) / 2
if half_width <= tol_x:
return c

if fa * fc < 0:
# Root is in [a, c]
b = c
else:
# Root is in [c, b]
a = c
fa = fc

return (a + b) / 2

Expand All @@ -238,7 +261,7 @@ fn secant[
xtol: Scalar[dtype] = 1e-8,
rtol: Scalar[dtype] = 1e-8,
maxiter: Int = 100,
) -> Scalar[dtype]:
) raises -> Scalar[dtype]:
"""Finds a root using the secant method.

A derivative-free method that approximates the derivative using finite
Expand All @@ -256,19 +279,29 @@ fn secant[
rtol: Relative tolerance for convergence.
maxiter: Maximum number of iterations.

Raises:
Error: If zero slope is encountered.

Returns:
The approximate root as a Scalar[dtype].
"""
var a: Scalar[dtype] = x0
var b: Scalar[dtype] = x1

for _ in range(maxiter):
var f0 = f[dtype](a, args)
var f1 = f[dtype](b, args)
var f0 = f(a, args)
var f1 = f(b, args)

var denom = f1 - f0
if denom == 0:
raise Error(
"Scijo [newton]: Secant method encountered zero slope (f1 - f0"
" == 0)."
)

var xn = b - (f1 * (b - a)) / (f1 - f0)
var xn = b - (f1 * (b - a)) / denom

var fxn = f[dtype](xn, args)
var fxn = f(xn, args)
var tol_x = max(xtol, rtol * abs(xn))
var tol_f = max(xtol, rtol * abs(fxn))

Expand All @@ -278,7 +311,7 @@ fn secant[
if abs(xn - b) <= tol_x:
return xn

a = b
b = xn
a = x1

return x1
return b
4 changes: 4 additions & 0 deletions scijo/optimize/utility.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,10 @@
Data structures for returning results from optimization and root-finding routines.
"""

# ===----------------------------------------------------------------------=== #
# RootResults
# ===----------------------------------------------------------------------=== #


struct RootResults[dtype: DType = DType.float64]():
"""Result structure for scalar root-finding operations.
Expand Down
13 changes: 10 additions & 3 deletions tests/test_differentiate.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,9 @@ fn cubic_function[

fn sin_function[
dtype: DType
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype] where dtype.is_floating_point():
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[
dtype
] where dtype.is_floating_point():
"""
F(x) = sin(x), f'(x) = cos(x).
"""
Expand All @@ -52,7 +54,9 @@ fn sin_function[

fn cos_function[
dtype: DType
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype] where dtype.is_floating_point():
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[
dtype
] where dtype.is_floating_point():
"""
F(x) = cos(x), f'(x) = -sin(x).
"""
Expand All @@ -61,7 +65,9 @@ fn cos_function[

fn exp_function[
dtype: DType
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[dtype] where dtype.is_floating_point():
](x: Scalar[dtype], args: Optional[List[Scalar[dtype]]]) -> Scalar[
dtype
] where dtype.is_floating_point():
"""
F(x) = e^x, f'(x) = e^x.
"""
Expand Down Expand Up @@ -384,5 +390,6 @@ fn test_step_size_parameters() raises:
msg="Result should be consistent across step factors",
)


def main():
TestSuite.discover_tests[__functions_in_module()]().run()
1 change: 1 addition & 0 deletions tests/test_fft.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -263,5 +263,6 @@ fn test_error_conditions() raises:
except:
print("Non-power-of-2 error handling - PASSED")


def main():
TestSuite.discover_tests[__functions_in_module()]().run()
9 changes: 5 additions & 4 deletions tests/test_interpolate.mojo
Original file line number Diff line number Diff line change
Expand Up @@ -313,9 +313,7 @@ fn test_edge_cases() raises:
)

# Test with very small intervals
var x_small = nm.array[nm.f64](
[0.0, 1e-10, 2e-10], [3]
)
var x_small = nm.array[nm.f64]([0.0, 1e-10, 2e-10], [3])
var y_small = nm.array[nm.f64]([0.0, 1.0, 2.0], [3])
var interp_small = LinearInterpolator(x_small, y_small)

Expand All @@ -327,7 +325,9 @@ fn test_edge_cases() raises:

var small_test_point = 1.5e-10
var result_small = interp_small(small_test_point)
var scipy_small = Float64(py=py_interp_small(PythonObject(small_test_point)))
var scipy_small = Float64(
py=py_interp_small(PythonObject(small_test_point))
)
assert_almost_equal(
result_small,
scipy_small,
Expand Down Expand Up @@ -516,5 +516,6 @@ fn test_performance_comparison() raises:
msg="Large dataset interpolation should match SciPy",
)


def main():
TestSuite.discover_tests[__functions_in_module()]().run()
Loading
Loading