This post is meant to provide some summary about the things that have helped so far.
Indeed, the error of the forward differentiation I got was fixed by transforming things to abstract type. As @ChrisRackauckas wrote:
this helped me in changing the type of the initial conditions of the ODE solver from Float64 to abstract type.
Then I had to do more changes in the code to make sure everything gets abstract types.
Now, when I run optimize I don’t get the error, but the differentiation takes a long time (the code is running for more than 15 hours, without any output so far).
In the meantime I am trying to solve the problem using Optim.jl and BlackBoxOptim.jl. See this post for more information.
There is still a chance that NLOptControl.jl will be helpful, but I am not sure yet about it.
Update - automatic differentiation (ForwardDiff.jl) returns Nan
My optimization function needs to calculate elliptic integrals. After converting the code in Elliptic.jl to be able to work with automatic differentiation, I get that the derivative of the complete elliptic integral of the 1st kind is NaN. The derivative is actually given by other elliptic integrals (as can be seen in the last link at the end of the section).
Here’s the minimal code (note that the credit to the code of calculating the elliptic integral belongs to nolta’s package Elliptic.jl):
julia> function DRF(X::Real, Y::Real, Z::Real)
D1MACH1 = typeof(X)(1e-308)
D1MACH2 = typeof(X)(1e308)
D1MACH3 = typeof(X)(5e-16)/typeof(X)(2.0)
D1MACH4 = typeof(X)(5e-16)
D1MACH5 = log10(typeof(X)(2.))
ERRTOL = (typeof(X)(4.0)*D1MACH3)^typeof(X)(1.0/6.0)
LOLIM = typeof(X)(5.0) * D1MACH1
UPLIM = D1MACH2/typeof(X)(5.0)
C1 = typeof(X)(1.0/24.0)
C2 = typeof(X)(3.0/44.0)
C3 = typeof(X)(1.0/14.0)
ans = zero(X)
if min(X,Y,Z) < zero(X)
return ans, one(X)
end
if max(X,Y,Z) > UPLIM
return ans, typeof(X)(3)
end
if min(X+Y,X+Z,Y+Z) < LOLIM
return ans, typeof(X)(2.0)
end
XN = X
YN = Y
ZN = Z
MU = zero(X)
XNDEV = zero(X)
YNDEV = zero(X)
ZNDEV = zero(X)
while true
MU = (XN+YN+ZN)/typeof(X)(3.0)
XNDEV = typeof(X)(2.0) - (MU+XN)/MU
YNDEV = typeof(X)(2.0) - (MU+YN)/MU
ZNDEV = typeof(X)(2.0) - (MU+ZN)/MU
EPSLON = max(abs(XNDEV),abs(YNDEV),abs(ZNDEV))
if (EPSLON < ERRTOL) break end
XNROOT = sqrt(XN)
YNROOT = sqrt(YN)
ZNROOT = sqrt(ZN)
LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
XN = (XN+LAMDA)*typeof(X)(0.250)
YN = (YN+LAMDA)*typeof(X)(0.250)
ZN = (ZN+LAMDA)*typeof(X)(0.250)
end
E2 = XNDEV*YNDEV - ZNDEV*ZNDEV
E3 = XNDEV*YNDEV*ZNDEV
S = typeof(X)(1.0) + (C1*E2-typeof(X)(0.10)-C2*E3)*E2 + C3*E3
ans = S/sqrt(MU)
return ans, typeof(X)(0)
end
DRF (generic function with 1 method)
julia> function Kdiff(m::Real)
if m < one(m)
drf,ierr = DRF(zero(m), one(m) - m, one(m))
@assert ierr == zero(m)
return drf
elseif m == one(m)
return Inf
elseif isnan(m)
return NaN
else
throw(DomainError(m, "argument m not <= 1"))
end
end
Kdiff (generic function with 1 method)
julia> Kdiff(0.5)
1.854074677301372
julia> using ForwardDiff
julia> ForwardDiff.derivative(Kdiff, 0.5)
NaN
Is there a way to fix this Nan issue? or to provide somehow the derivative only to this part of the code? (I mean that my fitness function needs to calculate the elliptic integral at some point, but I do not know the derivatives of other parts of the fitness function).
In the meantime, I gave this fitness function to JuMP and it runs. I get some results which I try to compare to other algorithm, and they seem for now to make sense.
But, shouldn’t it be a problem that the derivative is Nan? In what way does it affect the solver (I am using Ipopt).