Fix broken default constructors (#29)

* Add tests for default constructors

* Add default value and fix docs

* Add default value and fix broken tests

* Add missing parameter value and fix typo

* Update Project.toml

---------

Co-authored-by: George Datseris <[email protected]>

Author: csimal

Project.toml

@@ -1,7 +1,7 @@
 name = "PredefinedDynamicalSystems"
 uuid = "31e2f376-db9e-427a-b76e-a14f56347a14"
 repo = "https://github.com/JuliaDynamics/PredefinedDynamicalSystems.jl.git"
-version = "1.2"
+version = "1.3"
 
 [deps]
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"

src/continuous_famous_systems.jl

@@ -85,7 +85,7 @@ end
     return SVector{3}(du1, du2, du3)
 end
 function chua_jacob(u, p, t)
-    return SMatrix{3,3}(-p[1]*(1 + chua_element_derivative(u[1], p[3], p[4])),1.0,0.0,p[1],-1.0,-p[2],0.0,1.0,0.0)    
+    return SMatrix{3,3}(-p[1]*(1 + chua_element_derivative(u[1], p[3], p[4])),1.0,0.0,p[1],-1.0,-p[2],0.0,1.0,0.0)
 end
 # Helper functions for Chua's circuit.
 function chua_element(x, m0, m1)
@@ -130,7 +130,7 @@ end
     return SVector{3}(du1, du2, du3)
 end
 function roessler_jacob(u, p, t)
-    return SMatrix{3,3}(0.0,1.0,u[3],-1.0,p[1],0.0,-1.0,0.0,u[1]-p[3])          
+    return SMatrix{3,3}(0.0,1.0,u[3],-1.0,p[1],0.0,-1.0,0.0,u[1]-p[3])
 end
 
 """
@@ -170,7 +170,7 @@ end
     sin_ϕ = sin(φ); cos_ϕ = cos(φ)
     sin_θ₂ = sin(u[3]); sin_θ₁ = sin(u[1])
     Δ = (p[4] + p[5]) - p[5] * cos_ϕ^2
-    
+
     du1 = u[2]
     du2 = (p[5] * (cos_ϕ * (p[2] * u[2]^2 * sin_ϕ + p[1] * sin_θ₂) +
             p[3] * u[4]^2 * sin_ϕ) -
@@ -295,7 +295,7 @@ end
 `N` is the chain length, `F` the forcing. Jacobian is created automatically.
 (parameter container only contains `F`)
 """
-function lorenz96(N::Int, u0 = range(0; length = N, step = 0.1); F=0.01)
+function lorenz96(N::Int = 4, u0 = range(0; length = N, step = 0.1); F=0.01)
     @assert N ≥ 4 "`N` must be at least 4"
     lor96 = Lorenz96{N}() # create struct
     return CoupledODEs(lor96, u0, [F])
@@ -398,7 +398,7 @@ function gissinger_rule(u, p, t)
     return SVector{3}(du1, du2, du3)
 end
 function gissinger_jacob(u, p, t)
-    μ, ν, Γ = p            
+    μ, ν, Γ = p
    return SMatrix{3,3}(μ,u[3],u[2],-u[3],-ν,u[1],-u[2],u[1],-1)
 end
 
@@ -433,7 +433,7 @@ function rikitake_jacob(u, p, t)
     μ, α = p
     x,y,z = u
 
-    return SMatrix{3,3}(-μ,z-α,-y,z,-μ,-x,y,x,0)                   
+    return SMatrix{3,3}(-μ,z-α,-y,z,-μ,-x,y,x,0)
 end
 
 """
@@ -472,8 +472,8 @@ function nosehoover_rule(u, p, t)
 end
 function nosehoover_jacob(u, p, t)
     x,y,z = u
-                     
-	return SMatrix{3,3}(0,-1,0,1,z,-2y,0,y,0)                     
+
+	return SMatrix{3,3}(0,-1,0,1,z,-2y,0,y,0)
 end
 
 """
@@ -613,7 +613,7 @@ end
 function ueda_jacob(u, p, t)
     x,y = u
     k, B = p
-    return SMatrix{2,2}(0,-3*x^2,1,-k)               
+    return SMatrix{2,2}(0,-3*x^2,1,-k)
 end
 
 
@@ -778,9 +778,9 @@ function stommel_thermohaline_jacob(x, p, t)
     η1, η2, η3 = p
     q = abs(T-S)
     if T ≥ S
-		return SMatrix{2,2}((-1 - 2T + S), -S,T,(-η3 - T + 2S))                     
+		return SMatrix{2,2}((-1 - 2T + S), -S,T,(-η3 - T + 2S))
     else
-        return SMatrix{2,2}((-1 + 2T - S), S,-T,(-η3 + T - 2S))                  
+        return SMatrix{2,2}((-1 + 2T - S), S,-T,(-η3 + T - 2S))
     end
 end
 
@@ -823,7 +823,7 @@ end
 end
 function lorenz84_rule_jacob(u, p, t)
     F, G, a, b = p
-	x, y, z = u                     
+	x, y, z = u
 	return SMatrix{3,3}(-a,y-b*z,b*y+z,-2y,x-1,b*x,-2z,-b*x,x-1)
 end
 
@@ -867,7 +867,7 @@ end
     return SVector{3}(dx, dy, dz)
 end
 function lorenzdl_rule_jacob(u, p, t)
-    x, y, z = u                      
+    x, y, z = u
     return SMatrix{3,3}(-1,-z,y,1,0,x,0,-x,0)
 end
 
@@ -909,7 +909,7 @@ end
 
 
 """
-    kuramoto(D = 20, u0 = range(0, 2π; length = D);
+    kuramoto(D = 25, u0 = range(0, 2π; length = D);
         K = 0.3, ω = range(-1, 1; length = D)
     )
 The Kuramoto model[^Kuramoto1975] of `D` coupled oscillators with equation
@@ -977,7 +977,7 @@ end
 function sprott_dissipative_conservative_jacob(u, p, t)
     a, b, c = p
     x, y, z = u
-                    
+
     return SMatrix{3,3}(a*y + z,-4x,c - 2x,1 + a*x,b*z,-2y,x,b*y,0)
 end
 
@@ -1188,15 +1188,15 @@ function hindmarshrose_rule(u, p, t)
 end
 function hindmarshrose_jacob(u, p, t)
     @inbounds begin
-        a,b,c,d,r,s, xr, I = p        
+        a,b,c,d,r,s, xr, I = p
         return SMatrix{3,3}(-3*a*u[1]^2 + 2*b*u[1],-2*d*u[1],r*s,1,-1,0,-1,0,-r)
     end
 end
 
 """
 ```julia
 hindmarshrose_two_coupled(u0=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6];
-a = 1.0, b = 3.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0,
+a = 1.0, b = 3.0, c=1.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0,
 k1 = -0.17, k2 = -0.17, k_el = 0.0, xv = 2.0)
 ```
 ```math
@@ -1214,7 +1214,7 @@ The default parameter values are taken from article "Dragon-king-like extreme ev
 coupled bursting neurons", DOI:https://doi.org/10.1103/PhysRevE.97.062311.
 """
 function hindmarshrose_two_coupled(u0=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6];
-			a = 1.0, b = 3.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0,
+			a = 1.0, b = 3.0, c=1.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0,
 			k1 = -0.17, k2 = -0.17, k_el = 0.0, xv = 2.0)
 	return CoupledODEs(hindmarshrose_coupled_rule, u0, [a, b, c, d, r, s, xr, I, k1, k2, k_el, xv])
 end
@@ -1225,11 +1225,11 @@ function hindmarshrose_coupled_rule(u, p, t)
     a, b, c, d, r, s, xr, I, k1, k2, k_el, xv = p
     x1, y1, z1, x2, y2, z2 = u
 
-    du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
+    du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - xv ) * sigma(x2) + k_el * ( x2 - x1 )
     du2 = c - d * x1 ^2 - y1
     du3 = r * ( s * ( x1 - xr ) - z1 )
 
-    du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
+    du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - xv ) * sigma(x1) + k_el * ( x1 - x2 )
     du5 = c - d * x2 ^2 - y2
     du6 = r * ( s * ( x2 - xr ) - z2 )
     return SVector(du1, du2, du3, du4, du5, du6)
@@ -1288,7 +1288,7 @@ end
 function stuartlandau_jacob(u, p, t)
     @inbounds begin
         μ, ω, b = p
-            
+
         return SMatrix{2,2}(μ - 3*u[1]^2 -u[2]^2 -2*b*u[1]*u[2],-2*u[1]*u[2] +ω +b*u[2]^2 +3*b*u[1]^2,
         -2*u[1]*u[2] -ω -b*u[1]^2 -3*b*u[2]^2,μ -u[1]^2 -3*u[2]^2 +2*b*u[1]*u[2])
     end

src/discrete_famous_systems.jl

@@ -111,7 +111,7 @@ The first `M` entries of the state are the angles, the last `M` are the momenta.
 """
 function coupledstandardmaps end
 using SparseArrays
-function coupledstandardmaps(M::Int, u0 = 0.001rand(2M);
+function coupledstandardmaps(M::Int = 2, u0 = 0.001rand(2M);
     ks = ones(M), Γ = 1.0)
 
     SV = SVector{M, Int}
@@ -130,7 +130,7 @@ function coupledstandardmaps(M::Int, u0 = 0.001rand(2M);
     sparseJ = sparse(J)
     p = vcat(ks, Γ)
     csm(sparseJ, u0, p, 0)
-    return DeterministicIteratedMap(csm, u0, p, csm, sparseJ)
+    return DeterministicIteratedMap(csm, u0, p)
 end
 struct CoupledStandardMaps{N}
     idxs::SVector{N, Int}
@@ -148,7 +148,7 @@ function (f::CoupledStandardMaps{N})(xnew::AbstractVector, x, p, n) where {N}
 
         xnew[i] = mod2pi(x[i] + xnew[i+N])
     end
-    return xnew
+    return nothing
 end
 function (f::CoupledStandardMaps{M})(
     J::AbstractMatrix, x, p, n) where {M}
@@ -250,23 +250,23 @@ x_n(1 + |2x_n|^{z-1}), & \\quad |x_n| \\le 0.5 \\\\
 [2] : Meyer et al., New. J. Phys **20** (2019)
 """
 function pomeau_manneville(u0 = 0.2, z = 2.5)
-    return DeterministicIteratedMap(pm_rule, u0, [z], pm_jac)
+    return DeterministicIteratedMap(pm_rule, SVector{1}(u0), [z])
 end
 function pm_rule(x, p, n)
-    if x < -0.5
-        -4x - 3
-    elseif -0.5 ≤ x ≤ 0.5
-        @inbounds x*(1 + abs(2x)^(p[1]-1))
+    if x[1] < -0.5
+        SVector{1}(-4x[1] - 3)
+    elseif -0.5 ≤ x[1] ≤ 0.5
+        @inbounds SVector{1}(x[1]*(1 + abs(2x[1])^(p[1]-1)))
     else
-        -4x + 3
+        SVector{1}(-4x[1] + 3)
     end
 end
 function pomeau_manneville_jacob(x, p, n)
-    if x < -0.5
+    if x[1] < -0.5
         -4.0
-    elseif -0.5 ≤ x ≤ 0.5
+    elseif -0.5 ≤ x[1] ≤ 0.5
         @inbounds z = p[1]
-        0.5(x^2 * 2^z * (z-1)*abs(x)^(z-3) + 2^z * abs(x)^(z-1) + 2)
+        0.5(x[1]^2 * 2^z * (z-1)*abs(x[1])^(z-3) + 2^z * abs(x[1])^(z-1) + 2)
     else
         -4.0
     end
@@ -288,13 +288,13 @@ function's documentation string.
 [^Manneville1980]: Manneville, P. (1980). Intermittency, self-similarity and 1/f spectrum in dissipative dynamical systems. [Journal de Physique, 41(11), 1235–1243](https://doi.org/10.1051/jphys:0198000410110123500)
 """
 function manneville_simple(x0=0.4; ε = 0.1)
-    return DeterministicIteratedMap(manneville_f, x0, [ε], manneville_j)
+    return DeterministicIteratedMap(manneville_f, SVector{1}(x0), [ε])
 end
 
 function manneville_f(x, p, t)
     e = p[1]
-    y = (1+e)*x + (1-e)*x*x
-    return y%1
+    y = (1+e)*x[1] + (1-e)*x[1]*x[1]
+    return SVector{1}(y%1)
 end
 manneville_simple_jacob(x, p, n) = (1+p[1]) + (1-p[1])*2x
 
@@ -411,15 +411,15 @@ function tentmap(u0 = 0.25, μ = 2.0)
 end
 function tentmap_rule(x, p, n)
     μ = p[1]
-    if x < 0.5
-        μ*x
+    if x[1] < 0.5
+        SVector{1}(μ*x[1])
     else
-        μ*(1 - x)
+        SVector{1}(μ*(1 - x[1]))
     end
 end
 function tentmap_jacob(x, p, n)
     μ = p[1]
-    if x < -0.5
+    if x[1] < -0.5
         μ
     else
         -μ
@@ -440,14 +440,14 @@ The parameter β controls the dynamics of the map. Its Lyapunov exponent can be
 At β=2, it becomes the dyadic transformation, also known as the bit shift map, the 2x mod 1 map, the Bernoulli map or the sawtooth map. The typical trajectory for this case is chaotic, though there are countably infinite periodic orbits [^Ott2002].
 """
 function betatransformationmap(u0 = 0.25; β=2.0)
-    return DeterministicIteratedMap(betatransformation_rule, u0, [β])
+    return DeterministicIteratedMap(betatransformation_rule, SVector{1}(u0), [β])
 end
 function betatransformation_rule(x, p, n)
     @inbounds β = p[1]
-    if 0 ≤ x < 1/β
-        β*x
+    if 0 ≤ x[1] < 1/β
+        SVector{1}(β*x[1])
     else
-        β*x - 1
+        SVector{1}(β*x - 1)
     end
 end
 function betatransformation_jacob(x, p, n)
@@ -520,9 +520,9 @@ end
     a,b,c,d = p
     x,y = u
     t = c - d/(1 + x^2 + y^2)
-    aux = 2*d/(1+x^2+y^2)      
-    return SMatrix{2,2}(b*(cos(t)-x^2*sin(t)*aux -x*y*cos(t)*aux), b*(sin(t) +x^2*cos(t)*aux)-x*y*sin(t)*aux, 
-    b*(-sin(t) -x*y*sin(t)*aux -y^2*cos(t)*aux), b*(cos(t) -x*y*cos(t)*aux -y^2*sin(t)*aux))                
+    aux = 2*d/(1+x^2+y^2)
+    return SMatrix{2,2}(b*(cos(t)-x^2*sin(t)*aux -x*y*cos(t)*aux), b*(sin(t) +x^2*cos(t)*aux)-x*y*sin(t)*aux,
+    b*(-sin(t) -x*y*sin(t)*aux -y^2*cos(t)*aux), b*(cos(t) -x*y*cos(t)*aux -y^2*sin(t)*aux))
 end
 
 

test/constructors.jl

@@ -0,0 +1,81 @@
+
+
+@testset "Discrete Systems Default Constructors" begin
+    systems = [
+        :towel,
+        :standardmap,
+        :coupledstandardmaps,
+        :henon,
+        :logistic,
+        :pomeau_manneville,
+        :manneville_simple,
+        :arnoldcat,
+        :grebogi_map,
+        :nld_coupled_logistic_maps,
+        :tentmap,
+        :betatransformationmap,
+        :rulkovmap,
+        :ikedamap,
+        :ulam,
+        ]
+        for system in systems
+            @test @eval PredefinedDynamicalSystems.$system() isa DeterministicIteratedMap
+        end
+end
+
+@testset "Continous Systems Default Constructors" begin
+    systems = [
+        :lorenz,
+        :chua,
+        :roessler,
+        :double_pendulum,
+        :henonheiles,
+        :qbh,
+        :lorenz96,
+        :duffing,
+        :shinriki,
+        :gissinger,
+        :rikitake,
+        :nosehoover,
+        :antidots,
+        :ueda,
+        :magnetic_pendulum,
+        :fitzhugh_nagumo,
+        :more_chaos_example,
+        :thomas_cyclical,
+        :stommel_thermohaline,
+        :lorenz84,
+        :lorenzdl,
+        :coupled_roessler,
+        :kuramoto,
+        :sprott_dissipative_conservative,
+        :hodgkinhuxley,
+        :vanderpol,
+        :lotkavolterra,
+        :hindmarshrose,
+        :hindmarshrose_two_coupled,
+        :stuartlandau_oscillator,
+        :forced_pendulum,
+        :riddled_basins,
+        :morris_lecar,
+        :sakarya,
+        :lorenz_bounded,
+        :swinging_atwood,
+        :guckenheimer_holmes,
+        :halvorsen,
+        :multispecies_competition,
+        :hyper_roessler,
+        :hyper_lorenz,
+        :hyper_qi,
+        :hyper_jha,
+        :hyper_wang,
+        :hyper_xu,
+        :hyper_bao,
+        :hyper_cai,
+        :hyper_lu,
+        :hyper_pang,
+    ]
+    for system in systems
+        @test @eval PredefinedDynamicalSystems.$system() isa CoupledODEs
+    end
+end

test/runtests.jl

@@ -2,3 +2,5 @@ using PredefinedDynamicalSystems
 using Test
 
 @test PredefinedDynamicalSystems.henon() isa DeterministicIteratedMap
+
+include("constructors.jl")