Merge pull request #99 from JuliaDynamics/invalidities

Project.toml

@@ -9,12 +9,11 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
@@ -38,40 +37,46 @@ ChaosToolsExt = ["ChaosTools"]
 CoupledSDEsBaisin = ["ChaosTools", "Attractors"]
 
 [compat]
-Attractors = "1.15.0"
-ChaosTools = "3.1.2"
+Aqua = "0.8.7"
+Attractors = "1.18"
+ChaosTools = "3.1"
+Dates = ">=1.9.0"
 Dierckx = "0.5.3"
-DiffEqNoiseProcess = "5.22.0"
+DiffEqNoiseProcess = "5.22"
 DocStringExtensions = "0.9.3"
-Documenter = "1.4.1"
-DrWatson = "2.14.1"
+Documenter = "^1.4.1"
 DynamicalSystemsBase = "3.7.1"
+ExplicitImports = "1.9"
 Format = "1"
 ForwardDiff = "0.10.36"
 HDF5 = "0.17.1"
-IntervalRootFinding = "0.5.11"
+IntervalArithmetic = "0.20"
+JET = "0.9"
 JLD2 = "0.4.46"
 Markdown = ">=1.9.0"
+ModelingToolkit = "9.17"
 Optim = "1.9.3"
-OrdinaryDiffEq = "6.74.1"
+OrdinaryDiffEq = "6.82"
 ProgressBars = "1.5.1"
 ProgressMeter = "1.10.0"
 Reexport = "1.2.2"
 StaticArrays = "1.9.3"
 Statistics = ">=1.9"
 StochasticDiffEq = "6.65.1"
-Symbolics = "5.26.0"
+Symbolics = "5.30"
 julia = "1.9"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "IntervalRootFinding", "ModelingToolkit", "Documenter", "Attractors", "ChaosTools", "StatsBase"]
+test = ["Test", "Random", "ModelingToolkit", "Documenter", "Attractors", "ChaosTools", "Statistics", "ExplicitImports", "Aqua", "JET"]

ext/CoupledSDEsBaisin.jl

@@ -8,6 +8,7 @@ using ProgressMeter
 using ProgressBars
 using DynamicalSystemsBase: ParallelDynamicalSystem
 using DocStringExtensions
+using LinearAlgebra
 
 include("basin/planeofbox.jl")
 include("basin/basinsofattraction.jl")

ext/basin/edgetrack.jl

@@ -233,6 +233,8 @@ vector `attractors` containing as elements the state vectors of the stable fixed
 and a parameter `eps` to control the mapping algorithm. For more info, see the
 [`docs`](https://juliadynamics.github.io/Attractors.jl/dev/attractors/#Attractors.AttractorsViaProximity).
 """
-function attractor_mapper(sys::CoupledSDEs, attractors, eps=0.01; kwargs...)
+function Attractors.AttractorsViaProximity(
+    sys::CoupledSDEs, attractors, eps=0.01; kwargs...
+)
     return Attractors.AttractorsViaProximity(CoupledODEs(sys), attrs, eps; kwargs...)
 end

src/CoupledSDEs.jl

@@ -1,7 +1,13 @@
-using DynamicalSystemsBase: CoupledODEs, isinplace, __init, SciMLBase, correct_state
-import DynamicalSystemsBase: successful_step, _set_parameter!, current_state
-using StochasticDiffEq: SDEProblem, SDEIntegrator
-using StochasticDiffEq: EM, SDEProblem
+using DynamicalSystemsBase.SciMLBase: __init
+using DynamicalSystemsBase:
+    DynamicalSystemsBase,
+    CoupledODEs,
+    isinplace,
+    SciMLBase,
+    correct_state,
+    _set_parameter!,
+    current_state
+using StochasticDiffEq: SDEProblem
 
 ###########################################################################################
 # DiffEq options
@@ -53,7 +59,7 @@ If you want to specify a solver, do so by using the keyword `alg`, e.g.:
 `using OrdinaryDiffEq` to access the solvers. The default `diffeq` is:
 
 ```julia
-$(CriticalTransitions.DEFAULT_DIFFEQ)
+$(DEFAULT_DIFFEQ)
 ```
 
 `diffeq` keywords can also include `callback` for [event handling
@@ -160,7 +166,7 @@ function CoupledSDEs(
 )
     return CoupledSDEs(
         dynamic_rule(ds),
-        prob.g,
+        g,
         current_state(ds),
         p,
         noise_strength;
@@ -177,7 +183,9 @@ $(TYPEDSIGNATURES)
 Converts a [`CoupledSDEs`](@ref) into [`CoupledODEs`](https://juliadynamics.github.io/DynamicalSystems.jl/stable/tutorial/#DynamicalSystemsBase.CoupledODEs)
 from DynamicalSystems.jl.
 """
-function CoupledODEs(sys::CoupledSDEs; diffeq=DynamicalSystemsBase.DEFAULT_DIFFEQ, t0=0.0)
+function DynamicalSystemsBase.CoupledODEs(
+    sys::CoupledSDEs; diffeq=DynamicalSystemsBase.DEFAULT_DIFFEQ, t0=0.0
+)
     return DynamicalSystemsBase.CoupledODEs(
         sys.integ.f, SVector{length(sys.integ.u)}(sys.integ.u), sys.p0; diffeq=diffeq, t0=t0
     )
@@ -197,7 +205,7 @@ SciMLBase.isinplace(::CoupledSDEs{IIP}) where {IIP} = IIP
 StateSpaceSets.dimension(::CoupledSDEs{IIP,D}) where {IIP,D} = D
 DynamicalSystemsBase.current_state(ds::CoupledSDEs) = current_state(ds.integ)
 
-function set_parameter!(ds::CoupledSDEs, args...)
+function DynamicalSystemsBase.set_parameter!(ds::CoupledSDEs, args...)
     _set_parameter!(ds, args...)
     u_modified!(ds.integ, true)
     return nothing
@@ -216,7 +224,7 @@ SciMLBase.step!(ds::CoupledSDEs, args...) = (step!(ds.integ, args...); ds)
 # Besides, within DynamicalSystems.jl the integration is never expected to terminate.
 # Nevertheless here we extend explicitly only for ODE stuff because it may be that for
 # other type of DEIntegrators a different step interruption is possible.
-function successful_step(integ::SciMLBase.AbstractSDEIntegrator)
+function DynamicalSystemsBase.successful_step(integ::SciMLBase.AbstractSDEIntegrator)
     rcode = integ.sol.retcode
     return rcode == SciMLBase.ReturnCode.Default || rcode == SciMLBase.ReturnCode.Success
 end

src/CriticalTransitions.jl

@@ -1,18 +1,53 @@
 module CriticalTransitions
 
-using Reexport
-@reexport using DynamicalSystemsBase
+# Base
+using Statistics: Statistics, mean
+using LinearAlgebra: LinearAlgebra, Diagonal, I, norm, tr, dot
+using StaticArrays: StaticArrays, SVector
+
+# core
+using DynamicalSystemsBase:
+    DynamicalSystemsBase,
+    ContinuousTimeDynamicalSystem,
+    StateSpaceSets,
+    dimension,
+    dynamic_rule
+using DiffEqNoiseProcess: DiffEqNoiseProcess
+using OrdinaryDiffEq: OrdinaryDiffEq, Tsit5
+using StochasticDiffEq:
+    StochasticDiffEq,
+    DiscreteCallback,
+    ODEProblem,
+    SDEFunction,
+    SOSRA,
+    remake,
+    solve,
+    step!,
+    terminate!,
+    u_modified!
+
+using ForwardDiff: ForwardDiff
+using IntervalArithmetic: interval, IntervalBox
+using Dierckx: Dierckx, ParametricSpline
+using Optim: Optim, LBFGS
+using Symbolics: Symbolics
+
+# io and documentation
+using Format: Format
+using Dates: Dates
+using Printf: Printf
+using Markdown: Markdown
+using DocStringExtensions: TYPEDSIGNATURES
+using HDF5: HDF5, h5open, push!
+using JLD2: JLD2, jldopen
+using ProgressBars: ProgressBars, tqdm
+using ProgressMeter: ProgressMeter
+
+# reexport
+using Reexport: @reexport
 @reexport using StaticArrays
-@reexport using OrdinaryDiffEq
 @reexport using StochasticDiffEq
 @reexport using DiffEqNoiseProcess
-@reexport using LinearAlgebra
-using Format, Dates, JLD2, HDF5, ProgressBars, ProgressMeter, DocStringExtensions
-using IntervalRootFinding
-using ForwardDiff
-using Symbolics
-using Optim, Dierckx
-using Printf, DrWatson, Dates, Statistics, Markdown
 
 include("extention_functions.jl")
 include("utils.jl")
@@ -31,16 +66,23 @@ using .CTLibrary
 
 # Core types
 export CoupledSDEs,
-    idfunc!, idfunc, add_noise_strength, noise_process, covariance_matrix, noise_strength
+    idfunc!,
+    idfunc,
+    add_noise_strength,
+    noise_process,
+    covariance_matrix,
+    noise_strength,
+    CoupledODEs
 
 # Methods
-export equilib, fixedpoints, basins, basinboundary, basboundary
+export equilib, basins, basinboundary, basboundary
 export simulate, relax
 export transition, transitions
 export fw_action, om_action, action, geometric_action
 export min_action_method, geometric_min_action_method
-export basins, basinboundary
-export edgetracking, bisect_to_edge, attractor_mapper
 export make_jld2, make_h5, intervals_to_box
-
+# export basins, basinboundary
+# export edgetracking, bisect_to_edge, AttractorsViaProximity
+# export fixedpoints
+# ^ extention tests needed
 end # module CriticalTransitions

src/largedeviations/action.jl

@@ -75,11 +75,11 @@ Computes the action functional specified by `functional` for a given CoupledSDEs
 * `functional = "FW"`: Returns the Freidlin-Wentzell action ([`fw_action`](@ref))
 * `functional = "OM"`: Returns the Onsager-Machlup action ([`om_action`](@ref))
 """
-function action(sys::CoupledSDEs, path::Matrix, time, functional; kwargs...)
+function action(sys::CoupledSDEs, path::Matrix, time, functional)
     if functional == "FW"
-        action = fw_action(sys, path, time; kwargs...)
+        action = fw_action(sys, path, time)
     elseif functional == "OM"
-        action = om_action(sys, path, time; kwargs...)
+        action = om_action(sys, path, time)
     end
     return action
 end;

src/largedeviations/geometric_min_action_method.jl

@@ -141,6 +141,8 @@ function heymann_vandeneijnden_step(sys::CoupledSDEs, path, N, L; tau=0.1, diff_
     lambdas, lambdas_prime = zeros(N), zeros(N)
     x_prime = path_velocity(path, 0:dx:L; order=diff_order)
 
+    A = inv(covariance_matrix(sys))
+
     for i in 2:(N - 1)
         lambdas[i] = anorm(drift(sys, path[:, i]), A) / anorm(path[:, i], A)
     end

src/trajectories/simulation.jl

@@ -12,8 +12,8 @@ This function uses the [`SDEProblem`](https://diffeq.sciml.ai/stable/types/sde_t
 For more info, see [`SDEProblem`](https://diffeq.sciml.ai/stable/types/sde_types/#SciMLBase.SDEProblem).
 
 """
-function simulate(sys::CoupledSDEs, T, init=current_state(sys);
-    alg=sys.integ.alg, kwargs...
+function simulate(
+    sys::CoupledSDEs, T, init=current_state(sys); alg=sys.integ.alg, kwargs...
 )
     prob = remake(referrenced_sciml_prob(sys); u0=init, tspan=(0, T))
     return solve(prob, alg; kwargs...)
@@ -34,10 +34,8 @@ For more info, see [`ODEProblem`](https://diffeq.sciml.ai/stable/types/ode_types
 For stochastic integration, see [`simulate`](@ref).
 
 """
-function relax(sys::CoupledSDEs, T, init=current_state(sys);
-    alg=Tsit5(), kwargs...
-)
+function relax(sys::CoupledSDEs, T, init=current_state(sys); alg=Tsit5(), kwargs...)
     sde_prob = referrenced_sciml_prob(sys)
     prob = ODEProblem{isinplace(sde_prob)}(dynamic_rule(sys), init, (0, T), sys.p0)
     return solve(prob, alg; kwargs...)
-end;
+end;

src/utils.jl

@@ -56,11 +56,11 @@ function intervals_to_box(bmin::Vector, bmax::Vector)
     intervals = []
     dim = length(bmin)
     for i in 1:dim
-        push!(intervals, bmin[i] .. bmax[i])
+        push!(intervals, interval(bmin[i], bmax[i]))
     end
     box = intervals[1]
     for i in 2:dim
-        box = box × intervals[i]
+        box = IntervalBox(box, intervals[i])
     end
     return box
 end;

systems/CTLibrary.jl

@@ -1,6 +1,8 @@
 module CTLibrary
 
-using CriticalTransitions
+using CriticalTransitions: CriticalTransitions, smoothabs
+using IntervalArithmetic: interval
+using StaticArrays: SA, SVector
 
 include("fitzhughnagumo.jl")
 include("truscottbrindley_mod.jl")
@@ -10,7 +12,7 @@ include("rooth.jl")
 include("stommel.jl")
 include("rivals.jl")
 
-export fitzhughnagumo, fitzhughnagumo!, stommel, rivals!, rival, cessi, rooth_smooth
+export fitzhugh_nagumo!, fitzhugh_nagumo, stommel, rivals!, rivals, cessi, rooth_smooth
 modifiedtruscottbrindleywithdimensions!, modifiedtruscottbrindleywithdimensions
 originaltruscottbrindley!, originaltruscottbrindley
 rampedoriginaltruscottbrindley!, rampedoriginaltruscottbrindley

systems/fitzhughnagumo.jl

@@ -34,10 +34,6 @@ function fitzhugh_nagumo(u, p, t)
     return SA[dx, dy]
 end
 
-# For backwards compatibility
-FitzHughNagumo(u, p, t) = fitzhugh_nagumo(u, p, t)
-FitzHughNagumo!(u, p, t) = fitzhugh_nagumo!(u, p, t)
-
 # """
 #     fhn_ϵσ(ϵ,σ)
 # A shortcut command for returning a CoupledSDEs of the FitzHugh Nagumo system in a default setup with additive isotropic noise.

systems/truscottbrindley_mod.jl

@@ -82,7 +82,7 @@ function rampedmodifiedtruscottbrindley!(du, u, p, t)
         P₁ *
         (α * (P / P₁) * (1 - β * (P / P₁)) - γ * (Z / Z₁) * (P / P₁)^2 / (1 + (P / P₁)^2))
     du[2] = ξ * Z₁ * ((Z / Z₁) * (P / P₁)^2 / (1 + (P / P₁)^2) - (Z / Z₁)^2)
-    return du[3] = t ∈ Ttrans .. (Ttrans + Tramp) ? v : 0
+    return du[3] = t ∈ interval(Ttrans, (Ttrans + Tramp)) ? v : 0
 end
 
 """
@@ -99,7 +99,7 @@ function rampedmodifiedtruscottbrindley(u, p, t)
         P₁ *
         (α * (P / P₁) * (1 - β * (P / P₁)) - γ * (Z / Z₁) * (P / P₁)^2 / (1 + (P / P₁)^2))
     dZ = ξ * Z₁ * ((Z / Z₁) * (P / P₁)^2 / (1 + (P / P₁)^2) - (Z / Z₁)^2)
-    dα = t ∈ Ttrans .. (Ttrans + Tramp) ? v : 0
+    dα = t ∈ interval(Ttrans, (Ttrans + Tramp)) ? v : 0
 
     return SA[dP, dZ, dα]
 end

systems/truscottbrindley_orig.jl

@@ -46,7 +46,7 @@ function rampedoriginaltruscottbrindley!(du, u, p, t)
 
     du[1] = r * P * (1 - P / K) - Rₘ * Z * P^2 / (α^2 + P^2)
     du[2] = γ * Rₘ * Z * P^2 / (α^2 + P^2) - μ * Z
-    return du[3] = t ∈ Ttrans .. (Ttrans + Tramp) ? v : 0
+    return du[3] = t ∈ interval(Ttrans, (Ttrans + Tramp)) ? v : 0
 end
 
 """
@@ -61,7 +61,7 @@ function rampedoriginaltruscottbrindley(u, p, t)
 
     dP = r * P * (1 - P / K) - Rₘ * Z * P^2 / (α^2 + P^2)
     dZ = γ * Rₘ * Z * P^2 / (α^2 + P^2) - μ * Z
-    dr = t ∈ Ttrans .. (Ttrans + Tramp) ? v : 0
+    dr = t ∈ interval(Ttrans, (Ttrans + Tramp)) ? v : 0
 
     return SA[dP, dZ, dr]
 end

systems/truscottbrindley_orig1.jl

@@ -46,7 +46,7 @@ function rampedoriginaltruscottbrindley1!(du, u, p, t)
 
     du[1] = (1 / γ) * (r * P * (1 - P / K) - Rₘ * Z * P^2 / (α^2 + P^2))
     du[2] = Rₘ * Z * P^2 / (α^2 + P^2) - μ * Z
-    return du[3] = t ∈ Ttrans .. (Ttrans + Tramp) ? v : 0
+    return du[3] = t ∈ interval(Ttrans, (Ttrans + Tramp)) ? v : 0
 end
 
 """
@@ -61,7 +61,7 @@ function rampedoriginaltruscottbrindley1(u, p, t)
 
     dP = (1 / γ) * (r * P * (1 - P / K) - Rₘ * Z * P^2 / (α^2 + P^2))
     dZ = Rₘ * Z * P^2 / (α^2 + P^2) - μ * Z
-    dr = t ∈ Ttrans .. (Ttrans + Tramp) ? v : 0
+    dr = t ∈ interval(Ttrans, (Ttrans + Tramp)) ? v : 0
 
     return SA[dP, dZ, dr]
 end