Merge pull request #14 from wi11dey/master

SecondOrderDDEProblem (and underlying DynamicalDDEProblem)

Author: ChrisRackauckas

src/SciMLBase.jl

@@ -546,6 +546,8 @@ export SplitSDEProblem
 export RODEProblem, RODESolution, SDEProblem
 export DAEProblem, DAESolution
 export DDEProblem
+export DynamicalDDEFunction, DynamicalDDEProblem,
+       SecondOrderDDEProblem
 export SDDEProblem
 export PDEProblem
 export IncrementingODEProblem

src/problems/dde_problems.jl

@@ -1,7 +1,12 @@
 """
 $(TYPEDEF)
 """
-struct DDEProblem{uType,tType,lType,lType2,isinplace,P,F,H,K} <:
+struct StandardDDEProblem end
+
+"""
+$(TYPEDEF)
+"""
+struct DDEProblem{uType,tType,lType,lType2,isinplace,P,F,H,K,PT} <:
                           AbstractDDEProblem{uType,tType,lType,isinplace}
   f::F
   u0::uType
@@ -13,18 +18,20 @@ struct DDEProblem{uType,tType,lType,lType2,isinplace,P,F,H,K} <:
   kwargs::K
   neutral::Bool
   order_discontinuity_t0::Int
+  problem_type::PT
 
   @add_kwonly function DDEProblem{iip}(f::AbstractDDEFunction{iip}, u0, h, tspan, p=NullParameters();
                                        constant_lags = (),
                                        dependent_lags = (),
                                        neutral = f.mass_matrix !== I && det(f.mass_matrix) != 1,
                                        order_discontinuity_t0 = 0,
+                                       problem_type = StandardDDEProblem(),
                                        kwargs...) where {iip}
     _tspan = promote_tspan(tspan)
     new{typeof(u0),typeof(_tspan),typeof(constant_lags),typeof(dependent_lags),isinplace(f),
-        typeof(p),typeof(f),typeof(h),typeof(kwargs)}(
+        typeof(p),typeof(f),typeof(h),typeof(kwargs),typeof(problem_type)}(
           f, u0, h, _tspan, p, constant_lags, dependent_lags, kwargs, neutral,
-          order_discontinuity_t0)
+          order_discontinuity_t0, problem_type)
   end
 
   function DDEProblem{iip}(f::AbstractDDEFunction{iip}, h, tspan::Tuple, p=NullParameters();
@@ -43,3 +50,112 @@ DDEProblem(f, args...; kwargs...) =
 
 DDEProblem(f::AbstractDDEFunction, args...; kwargs...) =
   DDEProblem{isinplace(f)}(f, args...; kwargs...)
+
+"""
+$(TYPEDEF)
+"""
+abstract type AbstractDynamicalDDEProblem end
+
+"""
+$(TYPEDEF)
+"""
+struct DynamicalDDEProblem{iip} <: AbstractDynamicalDDEProblem end
+
+# u' = f1(v,h)
+# v' = f2(t,u,h)
+"""
+    DynamicalDDEProblem(f::DynamicalDDEFunction,v0,u0,tspan,p=NullParameters(),callback=CallbackSet())
+
+Define a dynamical DDE problem from a [`DynamicalDDEFunction`](@ref).
+"""
+function DynamicalDDEProblem(f::DynamicalDDEFunction,v0,u0,h,tspan,p=NullParameters();dependent_lags=(),kwargs...)
+  DDEProblem(f,ArrayPartition(v0,u0),h,tspan,p;
+             problem_type=DynamicalDDEProblem{isinplace(f)}(),
+             dependent_lags=ntuple(i->(u,p,t)->dependent_lags[i](u[1],u[2],p,t),length(dependent_lags)),
+             kwargs...)
+end
+function DynamicalDDEProblem(f::DynamicalDDEFunction,h,tspan,p=NullParameters();kwargs...)
+  DynamicalDDEProblem(f,h(p,first(tspan))...,h,tspan,p;kwargs...)
+end
+function DynamicalDDEProblem(f1::Function,f2::Function,args...;kwargs...)
+  DynamicalDDEProblem(DynamicalDDEFunction(f1,f2),args...;kwargs...)
+end
+
+"""
+    DynamicalDDEProblem{isinplace}(f1,f2,v0,u0,h,tspan,p=NullParameters(),callback=CallbackSet())
+
+Define a dynamical DDE problem from the two functions `f1` and `f2`.
+
+# Arguments
+* `f1` and `f2`: The functions in the DDE.
+* `v0` and `u0`: The initial conditions.
+* `h`: The initial history function.
+* `tspan`: The timespan for the problem.
+* `p`: Parameter values for `f1` and `f2`.
+* `callback`: A callback to be applied to every solver which uses the problem. Defaults to nothing.
+
+`isinplace` optionally sets whether the function is inplace or not.
+This is determined automatically, but not inferred.
+"""
+function DynamicalDDEProblem{iip}(f1::Function,f2::Function,args...;kwargs...) where iip
+  DynamicalDDEProblem(DynamicalDDEFunction{iip}(f1,f2),args...;kwargs...)
+end
+
+# u'' = f(du,u,h,p,t)
+"""
+$(TYPEDEF)
+"""
+struct SecondOrderDDEProblem{iip} <: AbstractDynamicalDDEProblem end
+function SecondOrderDDEProblem(f,args...;kwargs...)
+  iip = isinplace(f,6)
+  SecondOrderDDEProblem{iip}(f,args...;kwargs...)
+end
+
+"""
+    SecondOrderDDEProblem{isinplace}(f,du0,u0,h,tspan,p=NullParameters(),callback=CallbackSet())
+
+Define a second order DDE problem with the specified function.
+
+# Arguments
+* `f`: The function for the second derivative.
+* `du0`: The initial derivative.
+* `u0`: The initial condition.
+* `h`: The initial history function.
+* `tspan`: The timespan for the problem.
+* `p`: Parameter values for `f`.
+* `callback`: A callback to be applied to every solver which uses the problem. Defaults to nothing.
+
+`isinplace` optionally sets whether the function is inplace or not.
+This is determined automatically, but not inferred.
+"""
+function SecondOrderDDEProblem{iip}(f,args...;kwargs...) where iip
+  if iip
+    f2 = function (du,v,u,h,p,t)
+      du .= v
+    end
+  else
+    f2 = function (v,u,h,p,t)
+      v
+    end
+  end
+  DynamicalDDEProblem{iip}(f,f2,args...;problem_type=SecondOrderDDEProblem{iip}(),kwargs...)
+end
+function SecondOrderDDEProblem(f::DynamicalDDEFunction,args...;kwargs...)
+  iip = isinplace(f.f1, 6)
+  if f.f2.f === nothing
+    if iip
+      f2 = function (du,v,u,h,p,t)
+        du .= v
+      end
+    else
+      f2 = function (v,u,h,p,t)
+        v
+      end
+    end
+    return DynamicalDDEProblem(DynamicalDDEFunction{iip}(f.f1,f2;mass_matrix=f.mass_matrix,analytic=f.analytic),
+                               args...;problem_type=SecondOrderDDEProblem{iip}(),kwargs...)
+  else
+    return DynamicalDDEProblem(DynamicalDDEFunction{iip}(f.f1,f.f2;mass_matrix=f.mass_matrix,analytic=f.analytic),
+                               args...;problem_type=SecondOrderDDEProblem{iip}(),kwargs...)
+  end
+end

src/scimlfunctions.jl

@@ -109,6 +109,27 @@ struct DDEFunction{iip,F,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: Abstrac
   colorvec::TCV
 end
 
+"""
+$(TYPEDEF)
+"""
+struct DynamicalDDEFunction{iip,F1,F2,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} <: AbstractDDEFunction{iip}
+  f1::F1
+  f2::F2
+  mass_matrix::TMM
+  analytic::Ta
+  tgrad::Tt
+  jac::TJ
+  jvp::JVP
+  vjp::VJP
+  jac_prototype::JP
+  sparsity::SP
+  Wfact::TW
+  Wfact_t::TWt
+  paramjac::TPJ
+  syms::S
+  colorvec::TCV
+end
+
 """
 $(TYPEDEF)
 """
@@ -322,6 +343,21 @@ end
 (f::DDEFunction)(args...) = f.f(args...)
 (f::DDEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
 
+function (f::DynamicalDDEFunction)(u,h,p,t)
+  ArrayPartition(f.f1(u.x[1],u.x[2],h,p,t),f.f2(u.x[1],u.x[2],h,p,t))
+end
+function (f::DynamicalDDEFunction)(du,u,h,p,t)
+  f.f1(du.x[1],u.x[1],u.x[2],h,p,t)
+  f.f2(du.x[2],u.x[1],u.x[2],h,p,t)
+end
+function Base.getproperty(f::DynamicalDDEFunction, name::Symbol)
+  if name === :f
+    # Use the f property as an alias for calling the function itself, so DynamicalDDEFunction fits the same interface as DDEFunction as expected by the ODEFunctionWrapper in DelayDiffEq.jl.
+    return f
+  end
+  return getfield(f, name)
+end
+
 (f::SDEFunction)(args...) = f.f(args...)
 (f::SDEFunction)(::Type{Val{:analytic}},args...) = f.analytic(args...)
 (f::SDEFunction)(::Type{Val{:tgrad}},args...) = f.tgrad(args...)
@@ -933,6 +969,63 @@ DDEFunction{iip}(f::DDEFunction; kwargs...) where iip = f
 DDEFunction(f; kwargs...) = DDEFunction{isinplace(f, 5),RECOMPILE_BY_DEFAULT}(f; kwargs...)
 DDEFunction(f::DDEFunction; kwargs...) = f
 
+@add_kwonly function DynamicalDDEFunction{iip}(f1,f2,mass_matrix,analytic,tgrad,jac,jvp,vjp,
+                                               jac_prototype,sparsity,Wfact,Wfact_t,paramjac,
+                                               syms,colorvec) where iip
+  f1 = typeof(f1) <: AbstractDiffEqOperator ? f1 : DDEFunction(f1)
+  f2 = DDEFunction(f2)
+  DynamicalDDEFunction{isinplace(f2),typeof(f1),typeof(f2),typeof(mass_matrix),
+                       typeof(analytic),typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),
+                       typeof(jac_prototype),
+                       typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
+                       typeof(colorvec)}(f1,f2,mass_matrix,analytic,tgrad,jac,jvp,vjp,
+                                         jac_prototype,sparsity,Wfact,Wfact_t,paramjac,syms,colorvec)
+end
+function DynamicalDDEFunction{iip,true}(f1,f2;mass_matrix=I,
+                                        analytic=nothing,
+                                        tgrad=nothing,
+                                        jac=nothing,
+                                        jvp=nothing,
+                                        vjp=nothing,
+                                        jac_prototype=nothing,
+                                        sparsity=jac_prototype,
+                                        Wfact=nothing,
+                                        Wfact_t=nothing,
+                                        paramjac = nothing,
+                                        syms = nothing,
+                                        colorvec = nothing) where iip
+  DynamicalDDEFunction{iip,typeof(f1),typeof(f2),typeof(mass_matrix),
+                       typeof(analytic),
+                       typeof(tgrad),typeof(jac),typeof(jvp),typeof(vjp),typeof(jac_prototype),typeof(sparsity),
+                       typeof(Wfact),typeof(Wfact_t),typeof(paramjac),typeof(syms),
+                       typeof(colorvec)}(
+                         f1,f2,mass_matrix,analytic,tgrad,jac,jvp,vjp,jac_prototype,sparsity,
+                         Wfact,Wfact_t,paramjac,syms,colorvec)
+end
+function DynamicalDDEFunction{iip,false}(f1,f2;mass_matrix=I,
+                                         analytic=nothing,
+                                         tgrad=nothing,
+                                         jac=nothing,
+                                         jvp=nothing,
+                                         vjp=nothing,
+                                         jac_prototype=nothing,
+                                         sparsity=jac_prototype,
+                                         Wfact=nothing,
+                                         Wfact_t=nothing,
+                                         paramjac = nothing,
+                                         syms = nothing,
+                                         colorvec = nothing) where iip
+  DynamicalDDEFunction{iip,Any,Any,Any,Any,Any,Any,Any,
+                       Any,Any,Any,Any,Any}(
+                         f1,f2,mass_matrix,analytic,tgrad,
+                         jac,jvp,vjp,jac_prototype,sparsity,
+                         Wfact,Wfact_t,paramjac,syms,colorvec)
+end
+DynamicalDDEFunction(f1,f2=nothing; kwargs...) = DynamicalDDEFunction{isinplace(f1, 6)}(f1, f2; kwargs...)
+DynamicalDDEFunction{iip}(f1,f2; kwargs...) where iip =
+  DynamicalDDEFunction{iip,RECOMPILE_BY_DEFAULT}(DDEFunction{iip}(f1), DDEFunction{iip}(f2); kwargs...)
+DynamicalDDEFunction(f::DynamicalDDEFunction; kwargs...) = f
+
 function SDDEFunction{iip,true}(f,g;
                                 mass_matrix=I,
                                 analytic=nothing,
@@ -1135,14 +1228,14 @@ has_Wfact_t(f::Union{SplitFunction,SplitSDEFunction}) = has_Wfact_t(f.f1)
 has_paramjac(f::Union{SplitFunction,SplitSDEFunction}) = has_paramjac(f.f1)
 has_colorvec(f::Union{SplitFunction,SplitSDEFunction}) = has_colorvec(f.f1)
 
-has_jac(f::DynamicalODEFunction) = has_jac(f.f1)
-has_jvp(f::DynamicalODEFunction) = has_jvp(f.f1)
-has_vjp(f::DynamicalODEFunction) = has_vjp(f.f1)
-has_tgrad(f::DynamicalODEFunction) = has_tgrad(f.f1)
-has_Wfact(f::DynamicalODEFunction) = has_Wfact(f.f1)
-has_Wfact_t(f::DynamicalODEFunction) = has_Wfact_t(f.f1)
-has_paramjac(f::DynamicalODEFunction) = has_paramjac(f.f1)
-has_colorvec(f::DynamicalODEFunction) = has_colorvec(f.f1)
+has_jac(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_jac(f.f1)
+has_jvp(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_jvp(f.f1)
+has_vjp(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_vjp(f.f1)
+has_tgrad(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_tgrad(f.f1)
+has_Wfact(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_Wfact(f.f1)
+has_Wfact_t(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_Wfact_t(f.f1)
+has_paramjac(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_paramjac(f.f1)
+has_colorvec(f::Union{DynamicalODEFunction,DynamicalDDEFunction}) = has_colorvec(f.f1)
 
 has_jac(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_jac(f.f)
 has_jvp(f::Union{UDerivativeWrapper,UJacobianWrapper}) = has_jvp(f.f)