Merge pull request #102 from JuliaDynamics/add_tests

Project.toml

@@ -72,11 +72,12 @@ Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "ModelingToolkit", "Documenter", "Attractors", "ChaosTools", "Statistics", "ExplicitImports", "Aqua", "JET"]
+test = ["Test", "Random", "ModelingToolkit", "Documenter", "Attractors", "ChaosTools", "Statistics", "ExplicitImports", "Aqua", "JET", "IntervalArithmetic"]

src/CoupledSDEs_utils.jl

@@ -0,0 +1,62 @@
+diagonal_cov(l) = Diagonal(ones(l))
+
+"""
+$(TYPEDSIGNATURES)
+
+Translates the stochastic process specified in `sys` into the language required by the
+`SDEProblem` of `DynamicalSystems.jl`.
+"""
+function noise_process(sys::CoupledSDEs)
+    return sys.integ.W
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Gives the covariance matrix specified in `sys`.
+"""
+function covariance_matrix(sys::CoupledSDEs)
+    noise = noise_process(sys)
+    covariance = isnothing(noise) ? nothing : noise.covariance
+    if isnothing(noise) || isnothing(covariance)
+        return diagonal_cov(dimension(sys))
+    else
+        return covariance
+    end
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Gives the noise strength specified in `sys`.
+"""
+function noise_strength(sys::CoupledSDEs)
+    return sys.noise_strength
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Returns the drift field ``b(x)`` of the CoupledSDEs `sys` at the state vector `x`.
+"""
+function drift(sys::CoupledSDEs{IIP}, x) where {IIP}
+    # assumes the drift is time independent
+    f = dynamic_rule(sys)
+    if IIP
+        dx = similar(x)
+        f(dx, x, sys.p0, 0)
+        return dx
+    else
+        return f(x, sys.p0, 0)
+    end
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Computes the divergence of the drift field `sys.f` at the given point `x`.
+"""
+function div_drift(sys::CoupledSDEs, x)
+    b(x) = drift(sys, x)
+    return tr(ForwardDiff.jacobian(b, x))
+end;

src/CriticalTransitions.jl

@@ -27,7 +27,7 @@ using StochasticDiffEq:
     u_modified!
 
 using ForwardDiff: ForwardDiff
-using IntervalArithmetic: interval, IntervalBox
+using IntervalArithmetic: IntervalArithmetic, interval
 using Dierckx: Dierckx, ParametricSpline
 using Optim: Optim, LBFGS
 using Symbolics: Symbolics
@@ -52,6 +52,7 @@ using Reexport: @reexport
 include("extention_functions.jl")
 include("utils.jl")
 include("CoupledSDEs.jl")
+include("CoupledSDEs_utils.jl")
 include("io.jl")
 include("trajectories/simulation.jl")
 include("trajectories/transition.jl")

src/noiseprocesses/stochprocess.jl

@@ -1,62 +1 @@
-diagonal_cov(l) = Diagonal(ones(l))
 
-"""
-$(TYPEDSIGNATURES)
-
-Translates the stochastic process specified in `sys` into the language required by the
-`SDEProblem` of `DynamicalSystems.jl`.
-"""
-function noise_process(sys::CoupledSDEs)
-    return sys.integ.W
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Gives the covariance matrix specified in `sys`.
-"""
-function covariance_matrix(sys::CoupledSDEs)
-    noise = noise_process(sys)
-    covariance = isnothing(noise) ? nothing : noise.covariance
-    if isnothing(noise) || isnothing(covariance)
-        return diagonal_cov(dimension(sys))
-    else
-        return covariance
-    end
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Gives the noise strength specified in `sys`.
-"""
-function noise_strength(sys::CoupledSDEs)
-    return sys.noise_strength
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Returns the drift field ``b(x)`` of the CoupledSDEs `sys` at the state vector `x`.
-"""
-function drift(sys::CoupledSDEs{IIP}, x) where {IIP}
-    # assumes the drift is time independent
-    f = dynamic_rule(sys)
-    if IIP
-        dx = similar(x)
-        f(dx, x, sys.p0, 0)
-        return dx
-    else
-        return f(x, sys.p0, 0)
-    end
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Computes the divergence of the drift field `sys.f` at the given point `x`.
-"""
-function div_drift(sys::CoupledSDEs, x)
-    b(x) = drift(sys, x)
-    return tr(ForwardDiff.jacobian(b, x))
-end;

src/utils.jl

@@ -60,7 +60,7 @@ function intervals_to_box(bmin::Vector, bmax::Vector)
     end
     box = intervals[1]
     for i in 2:dim
-        box = IntervalBox(box, intervals[i])
+        box = IntervalArithmetic.cross(box, intervals[i])
     end
     return box
 end;

test/CoupledSDEs.jl

@@ -42,9 +42,11 @@
     @testset "Scalar noise Wiener" begin
         # Scalar noise Wiener
         # a single random variable is applied to all dependent variables
+        f!(du, u, p, t) = du .= 1.01u # deterministic part
+        σ = 0.25 # noise strength
         W = WienerProcess(0.0, 0.0, 0.0)
-        prob = SDEProblem(meier_stein, diagonal_noise(σ), zeros(2), (0.0, Inf), (); noise=W)
-        sde = CoupledSDEs(meier_stein, idfunc, zeros(2), (), σ; noise=W)
+        prob = SDEProblem(f!, diagonal_noise!(σ), zeros(2), (0.0, Inf), (); noise=W)
+        sde = CoupledSDEs(f!, idfunc!, zeros(2), (), σ; noise=W)
 
         @test sde.integ.sol.prob.f isa SDEFunction
         @test sde.integ.sol.prob.f.f.f == prob.f.f
@@ -60,6 +62,10 @@
         @test noise == W
         @test length(noise.dW) == 1
         @test W.covariance == nothing
+
+        sol = simulate(sde, 1.0; dt=0.01, alg=SOSRA())
+        # every vatiable has the same noise and dynamic_rule
+        @test all(sol[1, :] .≈ sol[2, :])
     end
 
     @testset "multiplicative noise Wiener" begin

test/largedeviations/gMAM.jl

@@ -5,38 +5,38 @@
     x_i = SA[sqrt(2 / 3), sqrt(2 / 27)]
     x_f = SA[0.001, 0.0]
     N = 100
-    res = geometric_min_action_method(fhn, x_i, x_f; N=75, maxiter=200)
+    res = geometric_min_action_method(fhn, x_i, x_f; N=75, maxiter=200, verbose=false)
     S = geometric_action(fhn, res[1][end])
     @test isapprox(S, 0.18, atol=0.01)
 end
 
-"""
 @testset "gMAM Meier Stein" begin
     function meier_stein(u, p, t) # out-of-place
         x, y = u
         dx = x - x^3 - 10 * x * y^2
         dy = -(1 + x^2) * y
-        [dx, dy]
+        return [dx, dy]
     end
     σ = 0.25
-    sys = StochSystem(meier_stein, [], zeros(2), σ, idfunc, nothing, I(2), "WhiteGauss")
+    # sys = CoupledSDEs(fitzhugh_nagumo, idfunc, zeros(2), p, σ)
+    sys = CoupledSDEs(meier_stein, idfunc, zeros(2), (), σ)
 
     # initial path: parabola
-    xx = range(-1.0, 1.0, length = 30)
+    xx = range(-1.0, 1.0; length=30)
     yy = 0.3 .* (-xx .^ 2 .+ 1)
     init = Matrix([xx yy]')
 
     x_i = init[:, 1]
     x_f = init[:, end]
 
     @testset "LBFGS" begin
-        gm = geometric_min_action_method(sys, x_i, x_f, maxiter = 10, verbose = false)# runtest
-        gm = geometric_min_action_method(sys, init, maxiter = 100, verbose=false)
+        gm = geometric_min_action_method(sys, x_i, x_f; maxiter=10, verbose=false)
+        gm = geometric_min_action_method(sys, init; maxiter=100, verbose=false)
 
         path = gm[1][end]
         action_val = gm[2][end]
-        @test all(isapprox.(path[2, :][(end - 5):end], 0, atol = 1e-3))
-        @test all(isapprox.(action_val, 0.3375, atol = 1e-3))
+        @test all(isapprox.(path[2, :][(end - 5):end], 0, atol=1e-3))
+        @test all(isapprox.(action_val, 0.3375, atol=1e-3))
     end
 
     @testset "HeymannVandenEijnden" begin # broken
@@ -49,4 +49,3 @@ end
         # @test all(isapprox.(action_val, 0.3375, atol = 1e-3))
     end # HeymannVandenEijnden
 end # gMAM Meier Stein
-"""

test/runtests.jl

@@ -58,13 +58,9 @@ end
     include("largedeviations/gMAM.jl")
 end
 
-# @testset "utilities" begin
-#     include("utils.jl")
-# end
-
-# @testset "basin" begin
-#     include("basin/basin_boundary.jl")
-# end
+@testset "utilities" begin
+    include("utils.jl")
+end
 
 @testset "Trajactories" begin
     include("trajactories/simulate.jl")

test/trajactories/transition.jl

@@ -20,6 +20,12 @@
     ensemble = transitions(sys, fp1, fp2, 10)
     @test ensemble.success_rate ≈ 1.0
     @test ensemble.t_trans ≈ 4.493941793363376 atol = 1e-2
-    @test ensemble.t_res ≈ 5449.261866107592 skip = true # SEED is not working on github
-    @test length(ensemble.times) ≈ 11 skip = true # SEED is not working on github
+    # SEED is different on github
+    if ensemble.t_res ≈ 5299.98
+        @test length(ensemble.times) == 11
+        @test ensemble.t_res ≈ 5299.98 atol = 1e-1
+    else
+        @test length(ensemble.times) == 10
+        @test ensemble.t_res ≈ 4953.88 atol = 1e-1
+    end
 end