# Problems fitting a model using LsqFit.jl

Hi everyone,

It’s a long shot as the problem is quite niche, but I would still appreciate any help/comments from the community.

### Problem Description

To put it briefly, I am interested in fitting a nonlinear model S(\lambda; C^v; \boldsymbol{\theta}) to data pairs (S_i,\, \lambda_i). The other variable C^v is the solution of a nonlinear ODE

\dot{C}^v = G(C^v, \lambda; \boldsymbol{\theta}),

where \boldsymbol{\theta} are the parameters I am eventually interested in. It is basically a constitutive model from nonlinear solid mechanics. The most natural way is to integrate the ODE numerically, get C^v and then solve a nonlinear least-squares problem. I have hard-coded my own time-integrator (a 5th order RK) and it seems to work fine (I know \boldsymbol{\theta} beforehand and just ran the forward model to see if it was working fine). Using DifferentialEquations.jl is next on my list. For the fitting, I am trying to use LsqFit.jl and specifically it’s curve_fit function along with the RK5 integrator.

### Data

1.000000000000000000e+00,0.000000000000000000e+00
1.080078125000000000e+00,1.007812500000000000e+01
1.157226562500000000e+00,1.660937500000000000e+01
1.255859375000000000e+00,2.221875000000000000e+01
1.336914062500000000e+00,2.584375000000000000e+01
1.452148437500000000e+00,2.982812500000000000e+01
1.618164062500000000e+00,3.306250000000000000e+01
1.771484375000000000e+00,3.612500000000000000e+01
1.945312500000000000e+00,3.881250000000000000e+01
2.111328125000000000e+00,4.153125000000000000e+01
2.302734375000000000e+00,4.421875000000000000e+01
2.515625000000000000e+00,4.743750000000000000e+01
2.716796875000000000e+00,5.050000000000000000e+01
2.996093750000000000e+00,5.446875000000000000e+01
2.923828125000000000e+00,5.156250000000000000e+01
2.800781250000000000e+00,4.721875000000000000e+01
2.664062500000000000e+00,4.343750000000000000e+01
2.535156250000000000e+00,3.962500000000000000e+01
2.400390625000000000e+00,3.584375000000000000e+01
2.259765625000000000e+00,3.168750000000000000e+01
2.156250000000000000e+00,2.879687500000000000e+01
2.054687500000000000e+00,2.553125000000000000e+01
1.956054687500000000e+00,2.118750000000000000e+01
1.845703125000000000e+00,1.757812500000000000e+01
1.743164062500000000e+00,1.304687500000000000e+01
1.584960937500000000e+00,6.894531250000000000e+00
1.469726562500000000e+00,7.329101562500000000e-01



### Code (S11 is the function)

using LsqFit;
using Plots;
using DelimitedFiles;
using Optim;

function S11(λ, Cv11, params)
mu1, mu2, alph1, alph2, m1, m2, a1, a2, K1, K2, bta1, bta2, eta0, etaInf = params
l2 = λ^(-0.5)
I1 = λ^2. + 2/λ
Cv22 = Cv11^(-0.5)
Ie1 = λ^2/Cv11 + 2*l2^2/Cv22
p = ((mu1 * (I1/3.)^(alph1 - 1.) + mu2 * (I1/3.)^(alph2 - 1.)) * l2 + (m1 * (Ie1/3.)^(a1 - 1.) + m2 * (Ie1/3.)^(a2 - 1.)) * l2/Cv22) * l2
s1 = ((mu1 * (I1/3.)^(alph1 - 1.) + mu2 * (I1/3.)^(alph2 - 1.)) * λ + (m1 * (Ie1/3.)^(a1 - 1.) + m2 * (Ie1/3.)^(a2 - 1.)) * λ/Cv11) - p/λ
return s1;
end

function Cvdot11(λ, Cv11, params)
mu1, mu2, alph1, alph2, m1, m2, a1, a2, K1, K2, bta1, bta2, eta0, etaInf = params
@assert λ > 0 "Stretches can't be negative"
l2 = λ^(-0.5)
# @assert Cv11 > 0 "Cv can't be negative" # this is where it fails
Cv22 = Cv11^(-0.5)

Ie1 = λ ^ 2. /Cv11 + 2.0 * l2^2. /Cv22
A2 = (m1 * (Ie1/3.)^(a1 - 1.) + m2 * (Ie1/3.)^(a2 - 1.))
Iv1 = Cv11 + 2. * Cv22
Ie1 = λ^2/Cv11 + 2*l2^2/Cv22
Cvinv_C_sq = λ^4/Cv11^2 + 2. * l2 ^4/Cv22^2
Ie2 = 1/2. * (Ie1 ^ 2. - Cvinv_C_sq)
J2NEq = (Ie1^2/3. - Ie2) * A2^2.

etaK = etaInf + (eta0 - etaInf + K1 * (Iv1 ^ (bta1) - 3^(bta1)))/(1. + (K2 * J2NEq)^(bta2))
cvd11 = A2/etaK * (λ^2. - 1/3. * Ie1 * Cv11)
return cvd11;
end

function kterms(dt, Cv_n, λ, λ_n,  λ_half, λ_quarter, λ_three_quarter, params)

G1 = Cvdot11(λ_n, Cv_n, params)
# @assert G1 > 0.
Cv_G2 = Cv_n + G1 * dt / 2.0
G2 = Cvdot11(
λ_half,
Cv_G2, params)

# @assert G2 > 1.
Cv_G3 = Cv_n + (3.0 * G1 + G2) * dt / 16.0
G3 = Cvdot11(
λ_quarter,
Cv_G3,
params)

# @assert G3 > 1.
Cv_G4 = Cv_n + G3 * dt / 2.0
G4 = Cvdot11(
λ_half,
Cv_G4,
params)

# @assert G4 > 1.
Cv_G5 = Cv_n + 3.0 * dt / 16.0 * (-G2 + 2.0 * G3 + 3.0 * G4)
G5 = Cvdot11(
λ_three_quarter,
Cv_G5,
params)

# @assert G5 > 1.
Cv_G6 = Cv_n + dt / 7.0 * (G1 + 4.0 * G2 + 6.0 * G3 - 12.0 * G4 + 8.0 * G5)
G6 = Cvdot11(
λ,
Cv_G6,
params)

# @assert G6 > 1.
Cv = Cv_n + dt / 90.0 * (7.0 * G1 + 32.0 * G3 + 12.0 * G4 + 32.0 * G5 + 7.0 * G6)
return Cv;
end

function correct_params()
params = [
(13.54 * 10^3),
(1.08 * 10^3),
(1.),
(-2.474),
(5.42 * 10^3),
(20.78 * 10^3),
(-10.),
(1.948),
(3507 * 10^3),
(10^(-6)),
(1.852),
(0.26),
(7014 * 10^3),
(0.1 * 10^3)];
return params;
end

function fit(λvals, params)
λdot = 0.01;
S11v = zeros(size(λvals, 1));
Cvvals = ones(size(λvals, 1));
for i ∈ range(1, size(λvals, 1) - 1)
Δt = abs(λvals[i+1] - λvals[i])/λdot;
λn = λvals[i];
λnp1 = λvals[i+1];
l1half = λn + 0.5 * (λnp1 - λn);
l1quarter = λn + 0.25 * (λnp1 - λn)
l1three_quarter = λn + 0.75 * (λnp1 - λn)
Cvvals[i+1] = kterms(Δt, Cvvals[i], λnp1, λn, l1half, l1quarter, l1three_quarter, params)
S11v[i+1] = S11(λnp1, Cvvals[i+1], params)
end
return S11v;
end

function plot_correct_plot(data)
λdot = 0.01;
timevals = [0.0];
for i ∈ 1:length(data[:, 1])-1
Δt = abs(data[i+1, 1] - data[i, 1])/λdot;
push!(timevals, timevals[i] + Δt);
end
length(timevals)
size(data[:, 1])
plot(timevals, data[:, 1]);
savefig("time.png");
size(data[:, 1])

# trial figure
params = correct_params();
Δtf = 0.01;
λdotf = 0.01;
λfinal = 3.0;
Tf = (λfinal - 1.0)/λdotf;
nsteps = Tf/Δtf + 1;
λvals = LinRange(1.0, λfinal, Int(nsteps));
λvals = vcat(λvals, LinRange(λfinal-λdotf * Δtf, 1.0, Int(nsteps)));
λvals
plot(λvals)
savefig("λvals.png");

plot(λvals, fit(λvals, params)./1E3, label="Correct Parameters");
scatter!(data[:, 1], data[:, 2], label="Data", legend=:bottomright);
savefig("correct_params.png");
end

plot_correct_plot(data);

p0 = [13E3, 1E3, 0.9, -3., 5E3, 20E3, -12., -10., 5E6, 1.5E-6, 2., 1., 1E7, 5E1];
ib = [1E2, 1E2, 0.5, -3.5, 1E2, 1E2, -15., -15., 1E6, 1E-7, 0.5, 0.1, 1E6, 1E1];
ub = [2E4, 2E4, 1.5, 2.5, 1E5, 1E5, 3., 3., 1E8, 3E-6, 2.5, 1.0, 1E7, 2E2];
myfit = curve_fit(fit, data[:, 1], data[:, 2] .* 1E3, p0, lower=ib, upper=ub);


### Result

If I run plot_correct_data above, nothing goes wrong (which led me to believe that the time-integrator was OK). Because I am calling the same function (which does the time-integration) but with the right set of parameters. ### Error

But If I run curve_fit on the data, it seems to throw an error specifically at the assertion statement I had (I was getting a DomainError before and I traced it back to the second @assert statement inside the Cvdot11 function). Here is the earlier stacktrace

ERROR: LoadError: DomainError with -0.30365595683307434:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
 throw_exp_domainerror(x::Float64)
@ Base.Math .\math.jl:37
 ^
@ .\math.jl:909 [inlined]
 Cvdot11(λ::Float64, Cv11::Float64, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:22
 kterms(dt::Float64, Cv_n::Float64, λ::Float64, λ_n::Float64, λ_half::Float64, λ_quarter::Float64, λ_three_quarter::Float64, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:70
 fit(λvals::Vector{Float64}, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:110



Now, if I try to do something like (Cv + 0im)^(-0.5) I get an InexactError, namely

ERROR: LoadError: InexactError: Float64(0.5141891571782113 - 1.6389735916551083e-6im)
Stacktrace:
 Real
@ .\complex.jl:44 [inlined]
 convert
@ .\number.jl:7 [inlined]
 setindex!
@ .\array.jl:903 [inlined]
 fit(λvals::Vector{Float64}, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:112


which seems really weird since line 112 is

Cvvals[i+1] = kterms(Δt, Cvvals[i], λnp1, λn, l1half, l1quarter, l1three_quarter, params).

It would be great if you could point me to what could be going wrong and where. Thanks in advance!

The DomainError is caused because you are taking the square root of a negative real number, and the workaround is to cast it to complex:

julia> (-1)^0.5
ERROR: DomainError with -1.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
...

julia> (Complex(-1))^0.5  # workaroud
0.0 + 1.0im


The InexactError is caused because you are trying to store the complex number calculated above in a vector of real numbers. Your Cvvals is a vector of Float64, but kterms() is generating a complex number because it calculates the output using the complex number obtained by square-rooting a negative number as above.

In the error message, the output Float64(0.5141891571782113 - 1.6389735916551083e-6im) generated by kterms() seems to be nearly real, and the imaginary part seems to be just a rounding error. If you are sure that kterms() should return always real numbers, discard its imaginary part and take only the real part when storing in Cvvals.

Thanks @wsshin, that was really helpful.

Not only that kterm should always return a positive real number. I did throw in abs(real(Cv)) when returning from Cvdot11 and similarly for kterms and so on. In principle I shouldn’t need to do this, right (?) I have done a similar fitting process

And for the correct set of parameters, the functionfit is already returning what it should.

But after doing abs(real(.)) on the output, the resulting fit is better (still not up to the mark; but definitely a good starting point to fiddle around with the algorithms/constraints etc.)

### Previously obtained parameters ### Obtained from LsqFit 