Auto optimize when some functions depend on data?

[rjrosati]

This is probably a dumb question, but is something like this auto_optimize-able?

function odefunc(du,u,p,t)
    du[:] .= Eh(t).*u
end

prob = ODEProblem(odefunc,[0.1,0.2],(0.0,10.0))

I only know the function Eh(t) numerically (from a previous ODE solve).

[ChrisRackauckas]

I don't think it'll know what to do with that function in modelingtookitize, so it would skip over to other forms of analysis and probably find the sparsity at least.

[agoldin]

I found it refuses to optimize in my case (KeplerSolve is a numerical Float64 function, solved by iterating) :

prob_1opt,alg_ = auto_optimize(prob_1)
Try ModelingToolkitization
MethodError(KeplerSolve, (α₃, α₁ * t), 0x0000000000006dfa)
┌ Warning: ModelingToolkitization Approach Failed
└ @ AutoOptimize /Users/alexey/.julia/packages/AutoOptimize/4YFOU/src/AutoOptimize.jl:65
MethodError: no method matching KeplerSolve(::Operation, ::Operation)

Stacktrace:
[1] auto_optimize(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(rot_1d_expl!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Nothing; verbose::Bool, stiff::Bool, mtkify::Bool, sparsify::Bool, gpuify::Bool, static::Bool, gpup::Nothing) at /Users/alexey/.julia/packages/AutoOptimize/4YFOU/src/AutoOptimize.jl:67
[2] auto_optimize at /Users/alexey/.julia/packages/AutoOptimize/4YFOU/src/AutoOptimize.jl:37 [inlined] (repeats 2 times)
[3] top-level scope at In[110]:1
[4] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

[ChrisRackauckas]

I'd need to see code.

[rjrosati]

I think this reproduces it

using AutoOptimize
using DifferentialEquations

function odefunc1(du,u,p,t)
    du[:] .= u
end
prob = ODEProblem(odefunc1,[0.1,0.2],(0.0,10.0))
prob,_ = auto_optimize(prob)
sol = solve(prob)

numfunc(t) = sol(t)[end]

function odefunc2(du,u,p,t)
    du[:] .= numfunc(t).*u
end

prob = ODEProblem(odefunc2,[0.1,0.2],(0.0,10.0))
prob,_ = auto_optimize(prob)
sol = solve(prob)

The second auto_optimize needs mtkify=false, gpuify=false to work on my system.

[agoldin]

My code can be reduced, but so far:

# From here: https://github.com/agoldin/2006HY51
using DifferentialEquations;
using AutoOptimize;

T = 4.17*365.25

n = 2*π /T;   # mean motion of asteroid  rad/day, period of 4.17 years
e = 0.9684;
σ = 0.2;

M(t, n) = mod(n*t, 2pi);
Ta(E, e) = 2* atan(sqrt((1 + e)/(1 - e))* sin(E/2.),cos(E/2.))

# Read http://murison.alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf
# Efficiently solve M = y - esin(y)

function KeplerStart3(e::Float64, M::Float64)
    t34=e*e
    t35 = e*t34
    t33= cos(M)
    M+(-.5*t35+e+(t34+1.5*t33*t35)*t33)*sin(M);
end


function eps3(e::Float64,M::Float64,x::Float64)

    t1 = cos(x);
    t2 = -1+e*t1;
    t3 = sin(x);
    t4 = e*t3;
    t5 = -x+t4+M;
    t6 = t5/(.5*t5*t4/t2+t2);
    t5/((.5*t3 - 1.0/6*t1*t6)*e*t6+t2);
end;

function KeplerSolve( e::Float64, M::Float64)#; tol=1.0e-14 )
    tol=1.5e-14
    Mnorm = mod(M,2π)
    E0 = KeplerStart3(e,Mnorm);
    dE = tol + 1.0;
    count = 0;
    E = E0;
    while dE > tol 
        E = E0 - eps3(e,Mnorm,E0);
        dE = abs(E-E0);
        E0 = E;
        count = count + 1;
        if count >10
            return E
        end
  
    end
    return E
end 

function rot_1d_expl!(du, u, p,t)
    n,σ,e = p
    e_anom = KeplerSolve(e, n*t)

    du[1] = -3.0/2.0*n^2*σ/((1.0 - e*cos(e_anom))^3.0)*(sin(2*u[2] - 2*Ta(e_anom,e))) 
    du[2] = u[1] # u[1] - velocity, u[2] - angle
    nothing
end

tspan = [-T/2,T/2]
u0 = [800.0*n, π*1.5]
params = [n,σ, e]
save_at = range(tspan[1],stop=tspan[2],length=100000);
prob_1 = ODEProblem(rot_1d_expl!,u0, tspan,params )

prob_1opt,alg_ = auto_optimize(prob_1)

[agoldin]

I managed to express solution for eccentric anomaly as DAE. auto_optimize works if solution is expressed in form of mass matrix problem, but does not work for DAEProblem yet. Impressive piece of software though, it will hopefully get even better.

[ChrisRackauckas]

I managed to express solution for eccentric anomaly as DAE. auto_optimize works if solution is expressed in form of mass matrix problem, but does not work for DAEProblem yet.

Yes, right now it's just on ODEProblem, but when I get a second I'll extend it to some others.