Parameter estimation of an ODE in Julia is slower than in R

General Usage Performance
, , ,
19

I finally just opened it up. This is the easy optimizations, just eliminating the big stuff from the profile:

using DifferentialEquations
using Plots
using Optim
using DataFrames
using Distributions
using BenchmarkTools
using LinearAlgebra
using StatsBase
using CSV

data_s = DataFrame(CSV.File("data_sim.csv"));

const nplot = length(unique(data_s.ID))
const AGE = data_s.AGE
const HDOM = data_s.HDOM
const UNIQUE_AGES = unique(AGE)

function f(du, u, p, t)
    du .= p[1] .* (u) .* (1 / t .^ (1 + p[2]))
    nothing
end

function DEsol_l(Age, parms)

    time_int = (minimum(Age), maximum(Age))
    u0 = parms[1:nplot]
    prob = ODEProblem{true}(f, u0, time_int, @view parms[(nplot+1):(nplot+2)])
    sol = solve(prob, Tsit5(), saveat=UNIQUE_AGES, save_everystep=false)
    sol_v = vec(transpose(Array(sol)))
    return sol_v
end

function loglik1(time, var, θ)

    sigma = exp(θ[nplot+3])
    hpred = DEsol_l(time, θ)
    residual = hpred .- var
    result = sum(logpdf.(Normal(0, sigma), residual))

    return -result

end

function get_inits()
    inits = combine(groupby(data_s, :ID), :HDOM => minimum)
    inits = inits[!, 2]
    inits = [inits; 5; 0.5; 1]
end

const inits1 = get_inits();

@profview result1 = optimize(b -> loglik1(AGE, HDOM, b), inits1, NelderMead(), Optim.Options(iterations=1000000))

@time result1 = optimize(b -> loglik1(AGE, HDOM, b), inits1, NelderMead(), Optim.Options(iterations=1000000))
6 Likes