tspan should probably have a default value or warning

[elbert5770]

This is related to the termination criteria potentially stopping too early with default values of abstol and reltol. If the abstol or reltol are set low to prevent early termination, then if they can't be met the solve runs to tspan = infinity. This crashes with "ERROR: ArgumentError: matrix contains Infs or NaNs", which doesn't really give the sense that time has run to 1.446951554681122e307.

Having a default value of say tspan = 1000000 or a warning when t got large could be helpful.

The following code has abstol and reltol too small for FiniteDiff gradient, but just right for ForwardDiff gradient. Thus it crashes with FiniteDiff gradient due to running to t = 1e307.

using FiniteDiff
using ForwardDiff
using Optimization
using OptimizationOptimJL
using NonlinearSolve
using DifferentialEquations
using SteadyStateDiffEq
using Plots

function sys!(du, u, (p,p_fix,kf,plateau,kd,start_time,end_time), t)
    #@show t,forcing_function(t,plateau,kd,start_time,end_time),p,p_fix,kf
    du[1] = kf[1]*(1.0 - forcing_function(t,plateau,kd,start_time,end_time)) - p[1]*u[1]
    du[2] = p[1]*u[1] - p_fix[1]*u[2]
    du[3] = p_fix[1]*u[2] - p[2]*u[3]
    
     #   @show u
      #  @show du
      @show t,dt
    
    return nothing
end

function forcing_function(t,plateau,kd,start_time,end_time)
    if t < start_time
        return plateau*t
    else
        if t <= end_time
            return plateau
        else
            if t > end_time && t < Inf
                return plateau*(exp(-kd*(t-end_time)))
            else
                return 0.0
            end
        end
    end
end



function set_kf_outer(probss,u03)
    function set_kf_inner(kf_guess,(u0,p_guess,p_fix,plateau,kd,start_time,end_time))
        # Solve for steady state with current p_guess and an initial guess of kf and u0
        plateau = 0.0
        sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol = 1e-16,reltol = 1e-16))
        # Compare kf at steady state to measured u0[3]
        # @show sol.u
        # @show sol.resid
        # @show sol.prob
        @show (sol[3]-u03)^2
        return (sol[3]-u03)^2
    end
end

function main()
    #### Parameters and state variables
    kf0= [10.0] #true value of zero order production rate
    p0 = [1.0,0.25,0.5] #true values of parameters
    p_fix = [p0[2]] # a parameter that is known, note that the system is otherwise unidentifiable
    p_var = [p0[1],p0[3]] # Unknown parameters
    u0_guess = [1000.0,750.0,20.0] #initial guess at steady state of state variables, except u0[3] is known to be 20.0
    u03 = 20.0 # measured value
    plateau = 0.5 # Zero order production rate changes at time = 0, reaches a fractional change given by plateau at start_time
    kd = 1.0 # Rate constant for exponential decay back to original value starting at end_time
    start_time = 1.0
    end_time = 10.0

    #### Set up ODE problem
    t_end = 36.0 # Length of simulation
    tspan = (0.0,t_end) # timespan of simulation
    probODE = ODEProblem(sys!,u0_guess,tspan,(p_var,p_fix,kf0,plateau,kd,start_time,end_time))
    
    #### Set up steady state problem
    probss = SteadyStateProblem(probODE)
    sol_ss = solve(probss)
    u0_ss = sol_ss.u
    @show "u0_ss true",sol_ss.u
    
    #### Initial guesses
    p_guess = [0.1,3.0] #inital guess at values in p_val
    kf_guess = [100.0] #initial guess at kf0
    

    #### Set up optimization problem to determine kf for current value of p_guess
    set_kf = set_kf_outer(probss,u03)
    optkf = OptimizationFunction(set_kf,Optimization.AutoForwardDiff())
    prob_set_kf = OptimizationProblem(optkf,kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
    
    
    @show FiniteDiff.finite_difference_gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess)
    @show ForwardDiff.gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess )
    # @show set_kf(kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
    # @show set_kf(kf_guess,(u0_guess,[p_guess[1]*1.01,p_guess[2]],p_fix,plateau,kd,start_time,end_time))
    # @show set_kf(kf_guess,(u0_guess,[p_guess[1],p_guess[2]*1.01],p_fix,plateau,kd,start_time,end_time))
    
   

    return nothing
end

main()