Add BAOAB algorithm #397
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,39 @@ | ||
|
|
||
| struct BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEqConstantCache | ||
| k::uType | ||
| half::uEltypeNoUnits | ||
| c1::uEltypeNoUnits | ||
| c2::uEltypeNoUnits | ||
| end | ||
| @cache struct BAOABCache{uType,uEltypeNoUnits,rateNoiseType} <: StochasticDiffEqMutableCache | ||
| utmp::uType | ||
| dutmp::uType | ||
| k::uType | ||
| gtmp::uType | ||
| noise::rateNoiseType | ||
| half::uEltypeNoUnits | ||
| c1::uEltypeNoUnits | ||
| c2::uEltypeNoUnits | ||
| end | ||
|
|
||
| function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) | ||
| k = zero(rate_prototype.x[1]) | ||
| c1 = exp(-alg.gamma*dt) | ||
| c2 = sqrt(1 - c1^2) | ||
| BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2)) | ||
| end | ||
|
|
||
| function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}}) | ||
| dutmp = zero(u.x[1]) | ||
| utmp = zero(u.x[2]) | ||
| k = zero(rate_prototype.x[1]) | ||
|
|
||
| gtmp = zero(rate_prototype.x[1]) | ||
| noise = zero(rate_prototype.x[1]) | ||
|
|
||
| half = uEltypeNoUnits(1//2) | ||
| c1 = exp(-alg.gamma*dt) | ||
| c2 = sqrt(1 - c1^2) | ||
|
|
||
| BAOABCache(utmp, dutmp, k, gtmp, noise, half, uEltypeNoUnits(c1), uEltypeNoUnits(c2)) | ||
| end |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,81 @@ | ||
| function verify_f2(f, p, q, pa, t, integrator, ::BAOABConstantCache) | ||
| res = f(p, q, pa, t) | ||
| res != p && throwex(integrator) | ||
| end | ||
| function verify_f2(f, res, p, q, pa, t, integrator, ::BAOABCache) | ||
| f(res, p, q, pa, t) | ||
| res != p && throwex(integrator) | ||
| end | ||
| function throwex(integrator) | ||
| algn = typeof(integrator.alg) | ||
| throw(ArgumentError("Algorithm $algn is not applicable if f2(p, q, t) != p")) | ||
| end | ||
|
|
||
| function initialize!(integrator, cache::BAOABConstantCache) | ||
| @unpack t,dt,uprev,u,p,W = integrator | ||
| du1 = integrator.uprev.x[1] | ||
| u1 = integrator.uprev.x[2] | ||
|
|
||
| verify_f2(integrator.f.f2, du1, u1, p, t, integrator, cache) | ||
| cache.k .= integrator.f.f1(du1,u1,p,t) | ||
| end | ||
|
|
||
| function initialize!(integrator, cache::BAOABCache) | ||
| @unpack t,dt,uprev,u,p,W = integrator | ||
| du1 = integrator.uprev.x[1] | ||
| u1 = integrator.uprev.x[2] | ||
|
|
||
| verify_f2(integrator.f.f2, cache.k, du1, u1, p, t, integrator, cache) | ||
| integrator.f.f1(cache.k,du1,u1,p,t) | ||
| end | ||
|
|
||
| @muladd function perform_step!(integrator,cache::BAOABConstantCache,f=integrator.f) | ||
| @unpack t,dt,sqdt,uprev,u,p,W = integrator | ||
| @unpack k, half, c1, c2 = cache | ||
| du1 = uprev.x[1] | ||
| u1 = uprev.x[2] | ||
|
|
||
| # B | ||
| du2 = du1 + half*dt*k | ||
|
|
||
| # A | ||
| u2 = u1 + half*dt*du2 | ||
|
|
||
| # O | ||
| noise = integrator.g(u2,p,t+dt*half).*W.dW / sqdt | ||
| du3 = c1*du2 + c2*noise | ||
|
|
||
| # A | ||
| u = u2 + half*dt*du3 | ||
|
|
||
| # B | ||
| k .= f.f1(du3,u,p,t+dt) | ||
| du = du3 + half*dt*k | ||
|
|
||
| integrator.u = ArrayPartition((du, u)) | ||
| end | ||
|
|
||
| @muladd function perform_step!(integrator,cache::BAOABCache,f=integrator.f) | ||
| @unpack t,dt,sqdt,uprev,u,p,W = integrator | ||
| @unpack utmp, dutmp, k, gtmp, noise, half, c1, c2 = cache | ||
| du1 = uprev.x[1] | ||
| u1 = uprev.x[2] | ||
|
|
||
| # B | ||
| @.. dutmp = du1 + half*dt*k | ||
|
|
||
| # A | ||
| @.. utmp = u1 + half*dt*dutmp | ||
|
|
||
| # O | ||
| integrator.g(gtmp,utmp,p,t+dt*half) | ||
| @.. noise = gtmp*W.dW / sqdt | ||
| @.. dutmp = c1*dutmp + c2*noise | ||
|
|
||
| # A | ||
| @.. u.x[2] = utmp + half*dt*dutmp | ||
|
|
||
| # B | ||
| f.f1(k,dutmp,u.x[2],p,t+dt) | ||
| @.. u.x[1] = dutmp + half*dt*k | ||
| end |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,29 @@ | ||
| using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools, Random | ||
| Random.seed!(1) | ||
|
|
||
| u0 = zeros(2) | ||
| v0 = ones(2) | ||
| γ = 1 | ||
|
|
||
| f1_harmonic(v,u,p,t) = -u | ||
| f2_harmonic(v,u,p,t) = v | ||
| g(u,p,t) = 1 | ||
|
|
||
| f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t) | ||
| f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t) | ||
| g_iip(du,u,p,t) = du .= g(u,p,t) | ||
|
|
||
| ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g) | ||
| prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0)) | ||
| sol1 = solve(prob1,BAOAB(gamma=γ);dt=1/10,save_noise=true) | ||
|
|
||
| prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W)) | ||
| sol2 = solve(prob2,BAOAB(gamma=γ);dt=1/10) | ||
|
|
||
| @test sol1[:] ≈ sol2[:] | ||
|
|
||
| dts = (1/2) .^ (8:-1:4) | ||
|
|
||
| # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v. | ||
| sim1 = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false) | ||
|
There was a problem hiding this comment. That's definitely not enough trajectories to produce a weak convergence estimate. Do a local plot with like 5e5 trajectories and show me the There was a problem hiding this comment.
Here's the plot with 5e4 trajectories with tspan=(0.0,0.5). I seemed to be running into memory issues if I try more trajectories or a longer tspan. However, it does seem to be converged now. I also tried with 1e4 trajectories and the plot is only marginally different. Here's the estimate I obtain: julia> sim1.𝒪est
Dict{Any,Any} with 2 entries:
:weak_final => 1.45594
:final => 0.802602There was a problem hiding this comment. What if you make the dts a bit smaller? the strong convergence estimate seems to not entire its asymtopic regime until half way through. There was a problem hiding this comment.
This is with only 1e4 trajectories as it takes a very long time with smaller timesteps. The test dt here is 0.5^12. julia> sim1.𝒪est
Dict{Any,Any} with 2 entries:
:weak_final => 1.19034
:final => 0.777283This doesn't quite match up with the previous plot however so I guess it's not properly converged. I'm not sure I can run it for many more trajectories though. There was a problem hiding this comment. More trajectories won't help. You'd see that because it would flatline in a noisy fashion on the left end due to not having enough values for the sampling error. This doesn't show sampling eror. It shows 1st order convergence. There was a problem hiding this comment. Lower There was a problem hiding this comment. Oh ok, does this mean the algorithm is 1st order then? I don't know much about SDEs, I was just going off the plots in the paper, but that's converging observables calculated from the distributions obtained from the integration. Is that different? There was a problem hiding this comment. A lot of the convergence diagrams aren't well-behaved in their paper either. Can you try and recreate one and see if it's around the same? There was a problem hiding this comment. It turns out I had made a mistake in the algorithm which I've now fixed. The convergence is a lot more well behaved now and I was able to roughly reproduce one of their plots. I did 2000 trajectories with 1e6 steps and they did 2000 with 1e8 so I think it's safe to assume I'd converge to the same result given enough timesteps. |
||
| @test abs(sim1.𝒪est[:weak_final]-1) < 0.3 | ||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
this should get convergence tests