ERROR: LinearAlgebra.LAPACKException(5) when adding bounds to sinuisoidal fit

[Datseris]

Hi there, I'd appreciate some support on the following: I am trying to fit a function with arbitrary amount of sine modes into some data. The fit works if I don't include fit bounds to the parameters, but fails with bounds.

Here is what I have:

using Dates, LsqFit

function harmonic_fit(cyclic_time::Vector{<:Real}, y, cycles = (1, 2))
    @assert cycles[1] == 1
    function harmonic(x, p)
        r = fill(p[1], length(x))
        for (c, i) in enumerate(1:2:2length(cycles))
            @. r += p[i+1]*sin(cycles[c]*x + p[i+2])
        end
        return r
    end
    p0 = [mean(y), std(y), 0.0]
    ub = Float64[maximum(y), maximum(y)-minimum(y), 2π] # set bounds for the fit
    lb = Float64[minimum(y), 0, 0]
    for i in 2:length(cycles)
        push!(p0, 1.0, 0.0) # assume amplitude 1 and phase 0 for all others
        push!(ub, maximum(y)-minimum(y), 2π)
        push!(lb, 0, 0)
    end
    fit = LsqFit.curve_fit(harmonic, cyclic_time, Vector(y), p0; lower=lb, upper=ub)
    c = LsqFit.coef(fit)
    cσ = LsqFit.stderror(fit)
    f = harmonic(cyclic_time, c)
    return f, nrmse(y, f), c, cσ
end

and I fit it using e.g.:

cyclic_time(t) = @. (dayofyear(t) - 80)/365.25 * 2π
t = Date("2000-03-15"):Month(1):Date("2019-03-15")
ct = cyclic_time(t)
harmonic_fit(ct, cos.(ct))

and I get:

julia> harmonic_fit(ct, cos.(ct))
ERROR: LinearAlgebra.LAPACKException(5)
Stacktrace:
 [1] chklapackerror at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\lapack.jl:38 [inlined]
 [2] trtrs!(::Char, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\lapack.jl:3413
 [3] ldiv! at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\triangular.jl:656 [inlined]
 [4] inv at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\triangular.jl:703 [inlined]
 [5] inv(::Array{Float64,2}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\dense.jl:766
 [6] estimate_covar(::LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,1}}) at C:\Users\m300808\.julia\packages\LsqFit\NjkFI\src\curve_fit.jl:203
 [7] #stderror#49(::Float64, ::Int64, ::typeof(StatsBase.stderror), ::LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,1}})
at C:\Users\m300808\.julia\packages\LsqFit\NjkFI\src\curve_fit.jl:217
 [8] stderror at C:\Users\m300808\.julia\packages\LsqFit\NjkFI\src\curve_fit.jl:217 [inlined]

which is quite vauge. I also don't really understand why the error happens with parameter bounds but not without. The error happensi n the line cσ = LsqFit.stderror(fit).

[Magalame]

I'm getting the error even without bounds it seems.

It's particularly weird since J[1:end,1] is pretty much just a vector of 1, which is contradictory with the fact that the fit has converged

[pkofod]

I'm not sure. It seems to end up close to a bound...