Merge pull request #566 from ChrisRackauckas-Claude/fix-zero-type-any-error

Fix MethodError: no method matching zero(::Type{Any}) in SINDy solve

Author: ChrisRackauckas

lib/DataDrivenSparse/test/sparse_linear_solve.jl

@@ -70,3 +70,40 @@ end
         @test vec(rescoeff)≈[0.25; 0.0; -1.0; 0.11; 0.0; -0.5] atol=5e-2
     end
 end
+
+# Issue #564: Test that solve doesn't throw MethodError when coefficients are all zero
+# This can happen with very small data values or high regularization
+@testset "Zero coefficients handling (Issue #564)" begin
+    rng = StableRNG(1111)
+
+    # Test case 1: Very small data values that lead to zero coefficients after regularization
+    N = 3
+    X̂ = randn(rng, N, 50) * 1e-10
+    Ŷ = randn(rng, 1, 50) * 1e-10
+
+    @variables u[1:N]
+    b = polynomial_basis(u, 2)
+    basis = Basis(b, u)
+    problem = DirectDataDrivenProblem(X̂, Ŷ)
+
+    λ = 1e-1
+    opt = ADMM(λ)
+    options = DataDrivenCommonOptions()
+
+    # This should not throw MethodError: no method matching zero(::Type{Any})
+    result = @test_nowarn solve(problem, basis, opt, options = options)
+    @test result isa DataDrivenSolution
+    @test eltype(result.prob) == Float64
+
+    # Test case 2: High regularization that forces all coefficients to zero
+    X̂2 = randn(rng, N, 50)
+    Ŷ2 = randn(rng, 1, 50)
+    problem2 = DirectDataDrivenProblem(X̂2, Ŷ2)
+
+    λ_high = 1e10  # Very high regularization
+    opt_high = ADMM(λ_high)
+
+    result2 = @test_nowarn solve(problem2, basis, opt_high, options = options)
+    @test result2 isa DataDrivenSolution
+    @test eltype(result2.prob) == Float64
+end

src/basis/type.jl

@@ -532,7 +532,13 @@ This extends `getmetadata` in a way that all parameters have a numeric value.
 Values are unwrapped from symbolic wrappers to ensure compatibility with ODEProblem.
 """
 function get_parameter_values(x::Basis)
-    map(parameters(x)) do p
+    ps = parameters(x)
+    # Return a properly typed empty vector if there are no parameters
+    # to avoid creating Any[] which causes type instability issues
+    if isempty(ps)
+        return Float64[]
+    end
+    map(ps) do p
         val = if hasmetadata(p, Symbolics.VariableDefaultValue)
             Symbolics.getdefaultval(p)
         else

src/solution.jl

@@ -35,7 +35,23 @@ function DataDrivenSolution(b::AbstractBasis, p::AbstractDataDrivenProblem,
     ps = get_parameter_values(b)
     prob = remake_problem(p, p = ps)
 
-    rss = sum(abs2, get_implicit_data(prob) .- b(prob))
+    # Calculate residual sum of squares, handling potential type instability
+    # when basis has no equations (which can result in Any-typed matrices)
+    Y = get_implicit_data(prob)
+    Ŷ = b(prob)
+    T = eltype(p)
+    residuals = Y .- Ŷ
+    # Handle case where residuals may have element type Any
+    # (e.g., when basis has no equations due to all-zero coefficients)
+    if eltype(residuals) === Any
+        if isempty(residuals)
+            rss = zero(T)
+        else
+            rss = sum(abs2, T.(residuals))
+        end
+    else
+        rss = sum(abs2, residuals)
+    end
 
     return DataDrivenSolution{eltype(p)}(b,
         retcode,