From fe43d9d9300c5929bbb6cdb6f58db40ce5cb9a10 Mon Sep 17 00:00:00 2001
From: baiyi <421586319@qq.com>
Date: Fri, 8 Dec 2023 22:15:24 +0800
Subject: [PATCH 1/2] fix dt scaling when scale_noise==false

---
 src/caches/dynamical_caches.jl | 4 ++--
 src/perform_step/dynamical.jl  | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/caches/dynamical_caches.jl b/src/caches/dynamical_caches.jl
index 22dae0cec..5980526fc 100644
--- a/src/caches/dynamical_caches.jl
+++ b/src/caches/dynamical_caches.jl
@@ -20,7 +20,7 @@ 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
+  c2 = alg.scale_noise ? sqrt((1 - c1^2)/abs(dt)) : 1 # if scale_noise == false, c2 = 1
   BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2))
 end
 
@@ -34,7 +34,7 @@ function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototy
 
   half = uEltypeNoUnits(1//2)
   c1 = exp(-alg.gamma*dt)
-  c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
+  c2 = alg.scale_noise ? sqrt((1 - c1^2)/abs(dt)) : 1 # if scale_noise == false, c2 = 1
 
   tmp = zero(u)
 
diff --git a/src/perform_step/dynamical.jl b/src/perform_step/dynamical.jl
index 074cbc7dd..68c2d1fa5 100644
--- a/src/perform_step/dynamical.jl
+++ b/src/perform_step/dynamical.jl
@@ -42,7 +42,7 @@ end
   u2 = u1 + half*dt*du2
 
   # O
-  noise = integrator.g(u2,p,t+dt*half).*W.dW / sqdt
+  noise = integrator.g(u2,p,t+dt*half).*W.dW
   du3 = c1*du2 + c2*noise
 
   # A
@@ -69,7 +69,7 @@ end
 
   # O
   integrator.g(gtmp,utmp,p,t+dt*half)
-  @.. noise = gtmp*W.dW / sqdt
+  @.. noise = gtmp*W.dW
   @.. dutmp = c1*dutmp + c2*noise
 
   # A

From 24026eb9bc86f03cc6fbe9d2b9131363fc76f930 Mon Sep 17 00:00:00 2001
From: baiyi <421586319@qq.com>
Date: Fri, 15 Dec 2023 10:01:41 +0800
Subject: [PATCH 2/2] add test

---
 test/sde/sde_dynamical.jl | 40 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/test/sde/sde_dynamical.jl b/test/sde/sde_dynamical.jl
index a0a259810..8e4475f04 100644
--- a/test/sde/sde_dynamical.jl
+++ b/test/sde/sde_dynamical.jl
@@ -45,3 +45,43 @@ end
     sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
     @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
 end
+
+@testset "Scalar u, scale_noise=false" begin
+    u0 = 0
+    v0 = 1
+
+    ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
+    prob1 = DynamicalSDEProblem(ff_harmonic,v0,u0,(0.0,5.0))
+
+    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=γ,scale_noise=false),(1/2)^10;trajectories=Int(2e2),use_noise_grid=false)
+    display(sim1.𝒪est)
+    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+end
+
+@testset "Vector u, scale_noise=false" begin
+
+    u0 = zeros(2)
+    v0 = ones(2)
+
+    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,v0,u0,(0.0,5.0))
+    sol1 = solve(prob1,BAOAB(gamma=γ,scale_noise=false);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=γ,scale_noise=false);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=γ,scale_noise=false),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
+    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+end