Fix Python 3.14 stack overflow in NonlinearSolveSciPy tests

Project.toml

@@ -108,7 +108,7 @@ Reexport = "1.2.2"
 ReverseDiff = "1.15"
 SciMLLogging = "1.3"
 SIAMFANLEquations = "1.0.1"
-SciMLBase = "2.116"
+SciMLBase = "2.127"
 SimpleNonlinearSolve = "2.11"
 SparseArrays = "1.10"
 SparseConnectivityTracer = "1"

lib/NonlinearSolveBase/Project.toml

@@ -82,7 +82,7 @@ Preferences = "1.4"
 Printf = "1.10"
 RecursiveArrayTools = "3"
 ReverseDiff = "1.15"
-SciMLBase = "2.116"
+SciMLBase = "2.127"
 SciMLJacobianOperators = "0.1.1"
 SciMLLogging = "1.3.1"
 SciMLOperators = "1.7"

lib/NonlinearSolveBase/src/solve.jl

@@ -136,6 +136,15 @@ function solve_call(_prob, args...; merge_callbacks = true, kwargshandle = nothi
     end
 
     checkkwargs(kwargshandle; kwargs...)
+
+    # Check bounds support for problems with bounds
+    if (_prob isa SciMLBase.NonlinearProblem || _prob isa SciMLBase.NonlinearLeastSquaresProblem) &&
+       (hasfield(typeof(_prob), :lb) && hasfield(typeof(_prob), :ub)) &&
+       (_prob.lb !== nothing || _prob.ub !== nothing) &&
+       length(args) > 0 && !SciMLBase.allowsbounds(args[1])
+        error("Algorithm $(args[1]) does not support bounds. Use an algorithm with allowsbounds(alg) == true.")
+    end
+
     if isdefined(_prob, :u0)
         if _prob.u0 isa Array
             if !isconcretetype(RecursiveArrayTools.recursive_unitless_eltype(_prob.u0))

lib/NonlinearSolveSciPy/CondaPkg.toml

@@ -0,0 +1,8 @@
+# Pin Python to < 3.14 due to stack overflow detection bug in Python 3.14
+# See: https://github.com/JuliaPy/PythonCall.jl/issues/694
+# Python 3.14 introduced a buggy stack overflow detection mechanism that causes
+# false positives with negative memory values when interacting with Julia's task system.
+
+[deps]
+python = ">=3.9,<3.14"
+scipy = ""

lib/NonlinearSolveSciPy/Project.toml

@@ -4,6 +4,7 @@ authors = ["SciML"]
 version = "1.2.0"
 
 [deps]
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -15,6 +16,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 path = "../NonlinearSolveBase"
 
 [compat]
+CommonSolve = "0.2"
 ConcreteStructs = "0.2.3"
 Hwloc = "3"
 InteractiveUtils = "<0.0.1, 1"
@@ -23,7 +25,7 @@ PrecompileTools = "1.2"
 PythonCall = "0.9"
 ReTestItems = "1.24"
 Reexport = "1.2.2"
-SciMLBase = "2.116"
+SciMLBase = "2.127"
 Test = "1.10"
 julia = "1.10"
 

lib/NonlinearSolveSciPy/src/NonlinearSolveSciPy.jl

@@ -3,7 +3,7 @@ module NonlinearSolveSciPy
 using ConcreteStructs: @concrete
 using Reexport: @reexport
 
-using PythonCall: pyimport, pyfunc, Py
+using PythonCall: pyimport, pyfunc, pyconvert, Py
 
 const scipy_optimize = Ref{Union{Py, Nothing}}(nothing)
 const PY_NONE = Ref{Union{Py, Nothing}}(nothing)
@@ -19,7 +19,9 @@ function __init__()
     end
 end
 
+using CommonSolve
 using SciMLBase
+using SciMLBase: allowsbounds
 using NonlinearSolveBase: AbstractNonlinearSolveAlgorithm,
                           construct_extension_function_wrapper
 
@@ -99,7 +101,7 @@ Internal: wrap a Julia residual function into a Python callable
 """
 function _make_py_residual(f::F, p) where F
     return pyfunc(x_py -> begin
-        x = Vector{Float64}(x_py)
+        x = pyconvert(Vector{Float64}, x_py)
         r = f(x, p)
         return r
     end)
@@ -110,7 +112,7 @@ Internal: wrap a Julia scalar function into a Python callable
 """
 function _make_py_scalar(f::F, p) where F
     return pyfunc(x_py -> begin
-        x = Float64(x_py)
+        x = pyconvert(Float64, x_py)
         return f(x, p)
     end)
 end
@@ -122,36 +124,51 @@ function SciMLBase.__solve(
     # Construct Python residual
     py_f = _make_py_residual(prob.f, prob.p)
 
-    # Bounds handling (lb/ub may be missing)
-    has_lb = hasproperty(prob, :lb)
-    has_ub = hasproperty(prob, :ub)
-    if has_lb || has_ub
-        lb = has_lb ? getproperty(prob, :lb) : fill(-Inf, length(prob.u0))
-        ub = has_ub ? getproperty(prob, :ub) : fill(Inf, length(prob.u0))
+    # Bounds handling from problem fields
+    if prob.lb !== nothing || prob.ub !== nothing
+        lb = prob.lb !== nothing ? prob.lb : fill(-Inf, length(prob.u0))
+        ub = prob.ub !== nothing ? prob.ub : fill(Inf, length(prob.u0))
         bounds = (lb, ub)
     else
         bounds = nothing
     end
 
-    res = scipy_optimize[].least_squares(py_f, collect(prob.u0);
-        method = alg.method,
-        loss = alg.loss,
-        max_nfev = maxiters,
-        bounds = bounds === nothing ? PY_NONE[] : bounds,
-        kwargs...)
+    # Filter out Julia-specific kwargs that scipy doesn't understand
+    scipy_kwargs = Tuple(k => v for (k, v) in pairs(kwargs) if k ∉ (:alias, :verbose))
 
-    u_vec = Vector{Float64}(res.x)
-    resid = Vector{Float64}(res.fun)
+    # Call scipy with conditional bounds argument
+    if bounds === nothing
+        res = scipy_optimize[].least_squares(py_f, collect(prob.u0);
+            method = alg.method,
+            loss = alg.loss,
+            max_nfev = maxiters,
+            scipy_kwargs...)
+    else
+        res = scipy_optimize[].least_squares(py_f, collect(prob.u0);
+            method = alg.method,
+            loss = alg.loss,
+            max_nfev = maxiters,
+            bounds = bounds,
+            scipy_kwargs...)
+    end
+
+    u_vec = pyconvert(Vector{Float64}, res.x)
+    resid = pyconvert(Vector{Float64}, res.fun)
 
     u = prob.u0 isa Number ? u_vec[1] : reshape(u_vec, size(prob.u0))
 
-    ret = res.success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure
+    ret = pyconvert(Bool, res.success) ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure
+    nfev = try
+        pyconvert(Int, res.nfev)
+    catch
+        0
+    end
     njev = try
-        Int(res.njev)
+        pyconvert(Int, res.njev)
     catch
         0
     end
-    stats = SciMLBase.NLStats(res.nfev, njev, 0, 0, res.nfev)
+    stats = SciMLBase.NLStats(nfev, njev, 0, 0, nfev)
 
     return SciMLBase.build_solution(prob, alg, u, resid; retcode = ret,
         original = res, stats = stats)
@@ -163,35 +180,36 @@ function SciMLBase.__solve(prob::SciMLBase.NonlinearProblem, alg::SciPyRoot;
     f!, u0, resid = construct_extension_function_wrapper(prob; alias_u0)
 
     py_f = pyfunc(x_py -> begin
-        x = Vector{Float64}(x_py)
+        x = pyconvert(Vector{Float64}, x_py)
         f!(resid, x)
         return resid
     end)
 
     tol = abstol === nothing ? nothing : abstol
 
+    # Filter out Julia-specific kwargs that scipy doesn't understand
+    scipy_kwargs = Tuple(k => v for (k, v) in pairs(kwargs) if k ∉ (:alias, :verbose))
+
     res = scipy_optimize[].root(py_f, collect(u0);
         method = alg.method,
         tol = tol,
         options = Dict("maxiter" => maxiters),
-        kwargs...)
+        scipy_kwargs...)
 
-    u_vec = Vector{Float64}(res.x)
+    u_vec = pyconvert(Vector{Float64}, res.x)
     f!(resid, u_vec)
 
     u_out = prob.u0 isa Number ? u_vec[1] : reshape(u_vec, size(prob.u0))
 
-    ret = res.success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure
+    ret = pyconvert(Bool, res.success) ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure
     nfev = try
-        Int(res.nfev)
+        pyconvert(Int, res.nfev)
     catch
-        ;
         0
     end
     niter = try
-        Int(res.nit)
+        pyconvert(Int, res.nit)
     catch
-        ;
         0
     end
     stats = SciMLBase.NLStats(nfev, 0, 0, 0, niter)
@@ -200,35 +218,36 @@ function SciMLBase.__solve(prob::SciMLBase.NonlinearProblem, alg::SciPyRoot;
         original = res, stats = stats)
 end
 
-function SciMLBase.__solve(prob::SciMLBase.IntervalNonlinearProblem, alg::SciPyRootScalar;
+function CommonSolve.solve(prob::SciMLBase.IntervalNonlinearProblem, alg::SciPyRootScalar, args...;
         abstol = nothing, maxiters = 10_000, kwargs...)
     f = prob.f
     p = prob.p
     py_f = _make_py_scalar(f, p)
 
     a, b = prob.tspan
 
+    # Filter out Julia-specific kwargs that scipy doesn't understand
+    scipy_kwargs = Tuple(k => v for (k, v) in pairs(kwargs) if k ∉ (:alias, :verbose))
+
     res = scipy_optimize[].root_scalar(py_f;
         method = alg.method,
         bracket = (a, b),
         maxiter = maxiters,
         xtol = abstol,
-        kwargs...)
+        scipy_kwargs...)
 
-    u_root = res.root
+    u_root = pyconvert(Float64, res.root)
     resid = f(u_root, p)
 
-    ret = res.converged ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure
+    ret = pyconvert(Bool, res.converged) ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure
     nfev = try
-        Int(res.function_calls)
+        pyconvert(Int, res.function_calls)
     catch
-        ;
         0
     end
     niter = try
-        Int(res.iterations)
+        pyconvert(Int, res.iterations)
     catch
-        ;
         0
     end
     stats = SciMLBase.NLStats(nfev, 0, 0, 0, niter)
@@ -237,6 +256,11 @@ function SciMLBase.__solve(prob::SciMLBase.IntervalNonlinearProblem, alg::SciPyR
         original = res, stats = stats)
 end
 
+# Trait declarations
+SciMLBase.allowsbounds(::SciPyLeastSquares) = true
+SciMLBase.allowsbounds(::SciPyRoot) = false
+SciMLBase.allowsbounds(::SciPyRootScalar) = false
+
 @reexport using SciMLBase, NonlinearSolveBase
 
 export SciPyLeastSquares, SciPyLeastSquaresTRF, SciPyLeastSquaresDogbox,