Dynamical integrators

src/StochasticDiffEq.jl

@@ -169,7 +169,7 @@ end
 
   export TauLeaping, CaoTauLeaping
 
-  export BAOAB
+  export BAOAB, ABOBA, OBABO
 
   export StochasticDiffEqRODEAlgorithm, StochasticDiffEqRODEAdaptiveAlgorithm,
          StochasticDiffEqRODECompositeAlgorithm

src/alg_utils.jl

@@ -122,6 +122,8 @@ alg_order(alg::TauLeaping) = 1 // 1
 alg_order(alg::CaoTauLeaping) = 1 // 1
 
 alg_order(alg::BAOAB) = 1 // 1
+alg_order(alg::ABOBA) = 2 // 1
+alg_order(alg::OBABO) = 2 // 1
 
 alg_order(alg::SKenCarp) = 2 // 1
 alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = maximum(alg_order.(alg.algs))
@@ -225,6 +227,8 @@ alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::IIF1M) = true
 alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::IIF2M) = true
 alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = max((alg_compatible(prob, a) for a in alg.algs)...)
 alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::BAOAB) = is_diagonal_noise(prob)
+alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::ABOBA) = is_diagonal_noise(prob)
+alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::OBABO) = is_diagonal_noise(prob)
 
 function alg_compatible(prob::JumpProblem, alg::Union{StochasticDiffEqJumpAdaptiveAlgorithm,StochasticDiffEqJumpAlgorithm})
     prob.prob isa DiscreteProblem
@@ -259,6 +263,7 @@ alg_needs_extra_process(alg::RS2) = true
 alg_needs_extra_process(alg::PL1WM) = true
 alg_needs_extra_process(alg::NON) = true
 alg_needs_extra_process(alg::NON2) = true
+alg_needs_extra_process(alg::OBABO) = true
 
 OrdinaryDiffEq._alg_autodiff(alg::StochasticDiffEqNewtonAlgorithm{CS,AD,FDT,ST,CJ,Controller}) where {CS,AD,FDT,ST,CJ,Controller} = Val{AD}()
 OrdinaryDiffEq._alg_autodiff(alg::StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,FDT,ST,CJ,Controller}) where {CS,AD,FDT,ST,CJ,Controller} = Val{AD}()

src/algorithms.jl

@@ -93,7 +93,7 @@ Defaults to solving the Ito problem, but RKMilCommute(interpretation=:Stratonovi
 Uses a 1.5/2.0 error estimate for adaptive time stepping.
 Default: ii_approx=IICommutative() does not approximate the Levy area.
 """
-struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm 
+struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm
   interpretation::Symbol
   ii_approx::T
 end
@@ -869,3 +869,35 @@ struct BAOAB{T} <: StochasticDiffEqAlgorithm
   scale_noise::Bool
 end
 BAOAB(;gamma=1.0, scale_noise=true) = BAOAB(gamma, scale_noise)
+
+
+@doc raw"""
+Leimkuhler B., Matthews C., Robust and efficient configurational molecular sampling via
+Langevin dynamics, J. Chem. Phys. 138, 174102 (2013)
+DOI:10.1063/1.4802990
+
+```math
+du = vdt \\
+dv = f(v,u) dt - \gamma v dt + g(u) \sqrt{2\gamma} dW
+```
+"""
+struct ABOBA{T} <: StochasticDiffEqAlgorithm
+  gamma::T
+  scale_noise::Bool
+end
+ABOBA(;gamma=1.0, scale_noise=true) = ABOBA(gamma, scale_noise)
+
+@doc raw"""
+Leimkuhler B., Matthews C., Molecular Dynamics With Deterministic and Stochastic Numerical Methods,Interdisciplinary Applied Mathematics; Springer International Publishing: Cham, 2015; Vol. 39
+DOI:10.1007/978-3-319-16375-8
+
+```math
+du = vdt \\
+dv = f(v,u) dt - \gamma v dt + g(u) \sqrt{2\gamma} dW
+```
+"""
+struct OBABO{T} <: StochasticDiffEqAlgorithm
+  gamma::T
+  scale_noise::Bool
+end
+OBABO(;gamma=1.0, scale_noise=true) = OBABO(gamma, scale_noise)

src/caches/dynamical_caches.jl

@@ -1,42 +1,188 @@
+abstract type StochasticDynamicalEqConstantCache <: StochasticDiffEqConstantCache end # Pourquoi faire ça, Si c'est pour avoir une seul function de check dans initialize!
+abstract type StochasticDynamicalEqMutableCache <: StochasticDiffEqMutableCache end
 
-mutable struct BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEqConstantCache
+
+mutable struct BAOABConstantCache{uType,uEltypeNoUnits,uCoeffType, uCoeffMType} <: StochasticDynamicalEqConstantCache
   k::uType
   half::uEltypeNoUnits
-  c1::uEltypeNoUnits
-  c2::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
 end
-@cache struct BAOABCache{uType,uEltypeNoUnits,rateNoiseType,uTypeCombined} <: StochasticDiffEqMutableCache
+@cache struct BAOABCache{uType,uEltypeNoUnits,rateNoiseType,uCoeffType, uCoeffMType,uTypeCombined} <: StochasticDynamicalEqMutableCache
   utmp::uType
+  dumid::uType
   dutmp::uType
+  dunoise::uType
   k::uType
-  gtmp::uType
-  noise::rateNoiseType
+  gtmp::rateNoiseType
+  noise::uType
   half::uEltypeNoUnits
-  c1::uEltypeNoUnits
-  c2::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
   tmp::uTypeCombined
 end
 
 function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
   k = zero(rate_prototype.x[1])
-  c1 = exp(-alg.gamma*dt)
-  c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
-  BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2))
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U# if scale_noise == false, c2 = 1
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)# if scale_noise == false, c2 = 1
+  end
+  BAOABConstantCache(k, uEltypeNoUnits(1//2),c2, c2)
 end
 
 function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  dumid = zero(u.x[1])
+  dutmp = zero(u.x[1])
+  dunoise = zero(u.x[1])
+  utmp = zero(u.x[2])
+  k = zero(rate_prototype.x[1])
+
+  gtmp = zero(noise_rate_prototype)
+  noise = zero(rate_prototype.x[1])
+
+  half = uEltypeNoUnits(1//2)
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U# if scale_noise == false, c2 = 1
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)# if scale_noise == false, c2 = 1
+  end
+
+  tmp = zero(u)
+
+  BAOABCache(utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ, tmp)
+end
+
+
+
+mutable struct ABOBAConstantCache{uType,uEltypeNoUnits, uCoeffType, uCoeffMType} <: StochasticDynamicalEqConstantCache
+  k::uType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+end
+@cache struct ABOBACache{uType,uEltypeNoUnits,rateNoiseType,uCoeffType, uCoeffMType,uTypeCombined} <: StochasticDynamicalEqMutableCache
+  utmp::uType
+  dumid::uType
+  dutmp::uType
+  dunoise::uType
+  k::uType
+  gtmp::rateNoiseType
+  noise::uType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+  tmp::uTypeCombined
+end
+
+function alg_cache(alg::ABOBA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  k = zero(rate_prototype.x[1])
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
+   # if scale_noise == false, c2 = 1
+  ABOBAConstantCache(k, uEltypeNoUnits(1//2), c2, σ)
+end
+
+function alg_cache(alg::ABOBA,prob,u,ΔW,ΔZ,p,rate_prototype, noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
   dutmp = zero(u.x[1])
+  dumid = zero(u.x[1])
+  dunoise = zero(u.x[1])
   utmp = zero(u.x[2])
   k = zero(rate_prototype.x[1])
 
-  gtmp = zero(rate_prototype.x[1])
+  gtmp = zero(noise_rate_prototype)
   noise = zero(rate_prototype.x[1])
 
   half = uEltypeNoUnits(1//2)
-  c1 = exp(-alg.gamma*dt)
-  c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
+
+  tmp = zero(u)
+
+  ABOBACache(utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ, tmp)
+end
+
+
+
+
+mutable struct OBABOConstantCache{uType,rateNoiseType, uEltypeNoUnits, uCoeffType, uCoeffMType} <: StochasticDynamicalEqConstantCache
+  k::uType
+  gt::rateNoiseType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+end
+
+@cache struct OBABOCache{uType,uEltypeNoUnits,rateNoiseType,uCoeffType, uCoeffMType,uTypeCombined} <: StochasticDynamicalEqMutableCache
+  utmp::uType
+  dumid::uType
+  dutmp::uType
+  dunoise::uType
+  k::uType
+  gtmp::rateNoiseType
+  noise::uType
+  half::uEltypeNoUnits
+  c2::uCoeffType
+  σ::uCoeffMType
+  tmp::uTypeCombined
+end
+
+function alg_cache(alg::OBABO,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  k = zero(rate_prototype.x[1])
+  gt = zero(noise_rate_prototype)
+  half=uEltypeNoUnits(1//2)
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*half*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*half*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
+   # if scale_noise == false, c2 = 1
+  OBABOConstantCache(k, gt, half, c2, σ)
+end
+
+function alg_cache(alg::OBABO,prob,u,ΔW,ΔZ,p,rate_prototype, noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  dutmp = zero(u.x[1])
+  dumid = zero(u.x[1])
+  dunoise = zero(u.x[1])
+  utmp = zero(u.x[2])
+  k = zero(rate_prototype.x[1])
+
+  gtmp = zero(noise_rate_prototype)
+  noise = zero(rate_prototype.x[1])
+
+
+  half = uEltypeNoUnits(1//2)
+
+  if typeof(alg.gamma) <: AbstractMatrix
+      c2 = exp(-alg.gamma*half*dt)
+      σ = cholesky(I - alg.scale_noise*c2*transpose(c2)).U
+  else
+      c2 = exp.(-alg.gamma*half*dt)
+      σ = sqrt.(1 .- alg.scale_noise*c2.^2)
+  end
 
   tmp = zero(u)
 
-  BAOABCache(utmp, dutmp, k, gtmp, noise, half, uEltypeNoUnits(c1), uEltypeNoUnits(c2), tmp)
+  OBABOCache(utmp, dumid, dutmp, dunoise, k, gtmp, noise, half, c2, σ, tmp)
 end