Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

LAPACKException when using two_stage_method #108

Open
atrophiedbrain opened this issue Jul 21, 2019 · 1 comment · May be fixed by #114
Open

LAPACKException when using two_stage_method #108

atrophiedbrain opened this issue Jul 21, 2019 · 1 comment · May be fixed by #114

Comments

@atrophiedbrain
Copy link

atrophiedbrain commented Jul 21, 2019

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 atrophiedbrain changed the title LAPackException when using two_stage_method LAPACKException when using two_stage_method Jul 21, 2019
@atrophiedbrain
Copy link
Author

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.

@ChrisRackauckas ChrisRackauckas linked a pull request Jul 28, 2019 that will close this issue
# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant