LsqFit.jl error on ODE data fitting

question

#1

Hi all, my first post here so bare with me please.
I’m also new to the Julia language and enjoying it very much so far.

Now to the point.
I’ve been trying LsqFit and wanted to implement a simple curve fitting model using a single shooting method. My goal is to estimate the rate constants of a chemical kinetics model. My code is the following:

using Optim
using LsqFit
using DifferentialEquations

function floudas_one(dz_dt, z, phi, t)
    r_1 = phi[1]*z[1]
    r_2 = phi[2]*z[2]

    dz_dt[1] = -r_1
    dz_dt[2] = r_1 - r_2
end

ode_fun = floudas_one

data = Float64[1.0 0.57353 0.328937 0.188654 0.108198 0.0620545 0.0355906 0.0204132 0.011708 0.00671499; 
        0.0 0.401566 0.589647 0.659731 0.666112 0.639512 0.597179 0.54867 0.499168 0.451377]
t = Float64[0.0, 0.111111, 0.222222, 0.333333, 0.444444, 0.555556, 0.666667, 0.777778, 0.888889, 1.0]

function lsq_ss_estimator(time_array, phi)
    tspan = (t[1], t[end])
    ini_cond = data[:,1]
    oprob = ODEProblem(ode_fun, ini_cond, tspan, phi)
    osol  = solve(oprob, Tsit5(), saveat=t)
    estimated = reduce(hcat, osol.u)
    return estimated
end

p0 = [5., 5.]
fit = curve_fit(lsq_ss_estimator, t, data, p0)

However, I get the following error:

MethodError: no method matching mul!(::Array{Float64,1}, ::LinearAlgebra.Transpose{Float64,Array{Float64,2}}, ::Array{Float64,2})
Closest candidates are:
  mul!(::AbstractArray, !Matched::Number, ::AbstractArray) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:26
  
(...)

` in top-level scope at [base/none](#)

in curve_fit at [LsqFit/TBrzN/src/curve_fit.jl:87](#)

in #curve_fit#5 at [LsqFit/TBrzN/src/curve_fit.jl:89](#)

in lmfit at [LsqFit/TBrzN/src/curve_fit.jl:41](#)

in #lmfit#3 at [LsqFit/TBrzN/src/curve_fit.jl:44](#)

in lmfit at [LsqFit/TBrzN/src/curve_fit.jl:48](#)

in #lmfit#4 at [LsqFit/TBrzN/src/curve_fit.jl:48](#)

in levenberg_marquardt at [LsqFit/TBrzN/src/levenberg_marquardt.jl:36](#)

in #levenberg_marquardt#1 at [LsqFit/TBrzN/src/levenberg_marquardt.jl:107](#)
`

Since I’m new to Julia and LsqFit, I’m not so sure what I’m missing here.
Any help would be appreciated.

PS: I’m aware of the DiffEqParamEstim package and it does an awesome job, I just wanted to try a few things by myself and ran into this issue.


Optimization and arbitrary precision
#2

The function you pass to curve_fit must return a vector, but lsq_ss_estimator returns a 2 by 10 matrix.

You could

  1. replace return estimated by return vec(estimated) in lsq_ss_estimator and
  2. call fit = curve_fit(lsq_ss_estimator, t, vec(data), p0) instead of your original call to curve_fit.

I don’t understand if there is a fundamental reason curve_fit cannot fit a data set with arbitrary dimensions, not just a vector.


#3

Oh, I see. That did the trick!
Slowly but surely I’m finding my way in Julia.
Thanks @mzaffalon!