SciML / DiffEqParamEstim.jl

Easy scientific machine learning (SciML) parameter estimation with pre-built loss functions
https://docs.sciml.ai/DiffEqParamEstim/stable/
Other
61 stars 34 forks source link

LAPACKException when using two_stage_method #108

Open atrophiedbrain opened 5 years ago

atrophiedbrain commented 5 years ago

I am receiving a LAPACKException when using two_stage_method. A MVE is included in this post.

Exception:

chklapackerror at lapack.jl:38 [inlined] trtrs!(::Char, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at lapack.jl:3348 inv at triangular.jl:583 [inlined] inv(::Array{Float64,2}) at dense.jl:728 construct_estimated_solution_and_derivative!(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float64,2}, ::Function, ::Array{Float64,1}, ::Float64, ::Int64) at two_stage_method.jl:68

two_stage_method#43(::Symbol, ::Type, ::Bool, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,typeof(f2f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Array{Float64,1}, ::Array{Float64,2}) at two_stage_method.jl:86

two_stage_method(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,typeof(f2f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Array{Float64,1}, ::Array{Float64,2}) at two_stage_method.jl:78 top-level scope at none:0

Both calls to two_stage_method in the MVE below, one using parameters and one not, both return this LAPACKException.


using DifferentialEquations, DiffEqParamEstim
const time_step = 0.5
const times = collect(7.0:time_step:85.0)
const rc = 62
const connf = randn(rc,rc,length(times))
const dataf = randn(rc,length(times))

function f2f(du, u, p, t)
    a, b = p
    t_index = floor(Int,(t - times[1]) / time_step + 1)
    if t_index >= size(connf)[3]
        t_index = size(connf)[3] - 1
    end

    @inbounds for i in 1:rc
        t2 = 0.0
        @inbounds for j in 1:rc
            if j != i
                t2 += u[j].*connf[j, i, t_index]
            end
        end

        du[i] = u[i].*b[i] .+ t2.*a[i]
    end
end

function f2f2(du, u, t)
    du .= u .* t
end

tspan = (minimum(times), maximum(times))
const u0 = dataf[:, 1]

bs = repeat([-0.005], rc)
as = repeat([0.0002], rc)

ps = [as, bs]
prob = ODEProblem(f2f, u0, tspan, ps)
prob2 = ODEProblem(f2f2, u0, tspan)

res = two_stage_method(prob,times,dataf)
res = two_stage_method(prob2,times,dataf)
atrophiedbrain commented 5 years ago

The default Epanechnikov_kernel returns 0 if abs(t) > 0. All my time points are positive. When I replaced my time vector with const times = collect(7.0:time_step:85.0)/85.0 the two_stage_method function runs sorry it returns 0 if abs(t) > 1 so this line in two_stage_method was returning 0 for all of my time points: W[i] = kernel_function((tpoints[i]-t)/h)/h h = (n^(-1/5))(n^(-3/35))((log(n))^(-1/16)) where n is the number of time points this is a decaying function that is less than 1, so (tpoints[i]-t)/h will be greater than 1 for all of my times that are greater than one Switching to a kernel that doesn't return 0 for time values > 1 works, for example Logistic_Kernel.

@ChrisRackauckas commented that perhaps we should change the default kernel so that the default kernel supports time series > 1.