Merge pull request #186 from devmotion/dde_updates

Add support for AbstractContinuousCallback and arbitrary discontinuities

Author: devmotion

REQUIRE

@@ -1,5 +1,5 @@
 julia 0.6.0
-DiffEqBase 1.20.0
+DiffEqBase 1.23.0
 Parameters 0.5.0
 ForwardDiff 0.5.0
 GenericSVD 0.0.2

src/callbacks.jl

@@ -1,12 +1,12 @@
 # Use Recursion to find the first callback for type-stability
 
 # Base Case: Only one callback
-function find_first_continuous_callback(integrator,callback::ContinuousCallback)
+function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback)
   (find_callback_time(integrator,callback)...,1,1)
 end
 
 # Starting Case: Compute on the first callback
-function find_first_continuous_callback(integrator,callback::ContinuousCallback,args...)
+function find_first_continuous_callback(integrator, callback::AbstractContinuousCallback, args...)
   find_first_continuous_callback(integrator,find_callback_time(integrator,callback)...,1,1,args...)
 end
 

src/integrators/integrator_utils.jl

@@ -72,7 +72,7 @@ function modify_dt_for_tstops!(integrator)
   end
 end
 
-function savevalues!(integrator::ODEIntegrator,force_save=false)
+function savevalues!(integrator::ODEIntegrator,force_save=false,reduce_size=true)
   while !isempty(integrator.opts.saveat) && integrator.tdir*top(integrator.opts.saveat) <= integrator.tdir*integrator.t # Perform saveat
     integrator.saveiter += 1
     curt = pop!(integrator.opts.saveat)
@@ -132,7 +132,7 @@ function savevalues!(integrator::ODEIntegrator,force_save=false)
       copyat_or_push!(integrator.sol.alg_choice,integrator.saveiter,integrator.cache.current)
     end
   end
-  resize!(integrator.k,integrator.kshortsize)
+  reduce_size && resize!(integrator.k,integrator.kshortsize)
 end
 
 function postamble!(integrator)
@@ -302,7 +302,7 @@ function handle_callbacks!(integrator)
       continuous_modified,saved_in_cb = apply_callback!(integrator,continuous_callbacks[idx],time,upcrossing)
     end
   end
-  if !(typeof(discrete_callbacks)<:Tuple{})
+  if !integrator.force_stepfail && !(typeof(discrete_callbacks)<:Tuple{})
     discrete_modified,saved_in_cb = apply_discrete_callback!(integrator,discrete_callbacks...)
   end
   if !saved_in_cb
@@ -346,12 +346,14 @@ function apply_step!(integrator)
   end
 
   # Update fsal if needed
-  if isfsal(integrator.alg)
-    if !isempty(integrator.opts.d_discontinuities) && top(integrator.opts.d_discontinuities) == integrator.t
-      pop!(integrator.opts.d_discontinuities)
-      reset_fsal!(integrator)
-    elseif integrator.reeval_fsal || integrator.u_modified || (typeof(integrator.alg)<:DP8 && !integrator.opts.calck) || (typeof(integrator.alg)<:Union{Rosenbrock23,Rosenbrock32} && !integrator.opts.adaptive)
-      reset_fsal!(integrator)
+  if !isempty(integrator.opts.d_discontinuities) &&
+      top(integrator.opts.d_discontinuities) == integrator.t
+
+      handle_discontinuities!(integrator)
+      isfsal(integrator.alg) && reset_fsal!(integrator)
+  elseif isfsal(integrator.alg)
+    if integrator.reeval_fsal || integrator.u_modified || (typeof(integrator.alg)<:DP8 && !integrator.opts.calck) || (typeof(integrator.alg)<:Union{Rosenbrock23,Rosenbrock32} && !integrator.opts.adaptive)
+        reset_fsal!(integrator)
     else # Do not reeval_fsal, instead copy! over
       if isinplace(integrator.sol.prob)
         recursivecopy!(integrator.fsalfirst,integrator.fsallast)
@@ -362,6 +364,10 @@ function apply_step!(integrator)
   end
 end
 
+function handle_discontinuities!(integrator)
+    pop!(integrator.opts.d_discontinuities)
+end
+
 function calc_dt_propose!(integrator,dtnew)
   dtpropose = integrator.tdir*min(abs(integrator.opts.dtmax),abs(dtnew))
   dtpropose = integrator.tdir*max(abs(dtpropose),abs(integrator.opts.dtmin))

src/integrators/type.jl

@@ -1,10 +1,10 @@
-mutable struct DEOptions{uEltype,uEltypeNoUnits,tTypeNoUnits,tType,F2,F3,F4,F5,F6,tstopsType,ECType,SType,MI}
+mutable struct DEOptions{absType,relType,tTypeNoUnits,tType,F2,F3,F4,F5,F6,tstopsType,discType,ECType,SType,MI}
   maxiters::MI
   timeseries_steps::Int
   save_everystep::Bool
   adaptive::Bool
-  abstol::uEltype
-  reltol::uEltypeNoUnits
+  abstol::absType
+  reltol::relType
   gamma::tTypeNoUnits
   qmax::tTypeNoUnits
   qmin::tTypeNoUnits
@@ -17,7 +17,7 @@ mutable struct DEOptions{uEltype,uEltypeNoUnits,tTypeNoUnits,tType,F2,F3,F4,F5,F
   save_idxs::SType
   tstops::tstopsType
   saveat::tstopsType
-  d_discontinuities::tstopsType
+  d_discontinuities::discType
   userdata::ECType
   progress::Bool
   progress_steps::Int

src/solve.jl

@@ -71,7 +71,7 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}(
     error("Timespan is trivial")
   end
 
-  tstops_vec = collect(tType,Iterators.filter(x->tdir*tspan[1]<tdir*x≤tdir*tspan[end],Iterators.flatten((tstops,d_discontinuities,tspan[end]))))
+  tstops_vec = vec(collect(tType,Iterators.filter(x->tdir*tspan[1]<tdir*x≤tdir*tspan[end],Iterators.flatten((tstops,d_discontinuities,tspan[end])))))
 
   if tdir>0
     tstops_internal = binary_minheap(tstops_vec)
@@ -154,7 +154,7 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}(
     saveat_vec = collect(tType,tspan[1]+saveat:saveat:(tspan[end]-saveat))
     # Exclude the endpoint because of floating point issues
   else
-    saveat_vec = collect(tType,Iterators.filter(x->tdir*tspan[1]<tdir*x<tdir*tspan[end],saveat))
+    saveat_vec = vec(collect(tType,Iterators.filter(x->tdir*tspan[1]<tdir*x<tdir*tspan[end],saveat)))
   end
 
   if tdir>0
@@ -163,7 +163,7 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}(
     saveat_internal = binary_maxheap(saveat_vec)
   end
 
-  d_discontinuities_vec =  collect(tType,d_discontinuities)
+  d_discontinuities_vec = vec(collect(d_discontinuities))
 
   if tdir>0
     d_discontinuities_internal = binary_minheap(d_discontinuities_vec)
@@ -209,18 +209,21 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}(
     saveiter_dense = 0
   end
 
-  opts = DEOptions(maxiters,timeseries_steps,save_everystep,adaptive,abstol_internal,
-    reltol_internal,tTypeNoUnits(gamma),tTypeNoUnits(qmax),tTypeNoUnits(qmin),
-    tTypeNoUnits(qsteady_max),tTypeNoUnits(qsteady_min),
-    tTypeNoUnits(failfactor),tType(dtmax),tType(dtmin),internalnorm,save_idxs,
-    tstops_internal,saveat_internal,d_discontinuities_internal,
-    userdata,
-    progress,progress_steps,
-    progress_name,progress_message,
-    timeseries_errors,dense_errors,
-    tTypeNoUnits(beta1),tTypeNoUnits(beta2),tTypeNoUnits(qoldinit),dense,save_start,
-    callbacks_internal,isoutofdomain,unstable_check,verbose,calck,force_dtmin,
-    advance_to_tstop,stop_at_next_tstop)
+  opts = DEOptions{typeof(abstol_internal),typeof(reltol_internal),tTypeNoUnits,tType,
+                   typeof(internalnorm),typeof(callbacks_internal),typeof(isoutofdomain),
+                   typeof(progress_message),typeof(unstable_check),typeof(tstops_internal),
+                   typeof(d_discontinuities_internal),typeof(userdata),typeof(save_idxs),
+                   typeof(maxiters)}(
+                       maxiters,timeseries_steps,save_everystep,adaptive,abstol_internal,
+                       reltol_internal,tTypeNoUnits(gamma),tTypeNoUnits(qmax),
+                       tTypeNoUnits(qmin),tTypeNoUnits(qsteady_max),
+                       tTypeNoUnits(qsteady_min),tTypeNoUnits(failfactor),tType(dtmax),
+                       tType(dtmin),internalnorm,save_idxs,tstops_internal,saveat_internal,
+                       d_discontinuities_internal,userdata,progress,progress_steps,
+                       progress_name,progress_message,timeseries_errors,dense_errors,
+                       tTypeNoUnits(beta1),tTypeNoUnits(beta2),tTypeNoUnits(qoldinit),dense,
+                       save_start,callbacks_internal,isoutofdomain,unstable_check,verbose,
+                       calck,force_dtmin,advance_to_tstop,stop_at_next_tstop)
 
   progress ? (prog = Juno.ProgressBar(name=progress_name)) : prog = nothing