Merge pull request #2 from rusandris/ds-v3

Update functions to DynamicalSystems v3.0

Project.toml

@@ -5,23 +5,20 @@ version = "0.1.0"
 
 [deps]
 Cairo = "159f3aea-2a34-519c-b102-8c37f9878175"
-ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f"
 DiscreteMarkovChains = "8abcb7ef-b365-4f7b-ac38-56893fb62f9f"
+DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377"
 Permutations = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 
-[compat]
-ChaosTools = "2.9.0"
-DelayEmbeddings = "<2.6"
+

README.md

@@ -21,7 +21,7 @@ using StateTransitionNetworks
 Constructing state-transition network for the [Henon map](https://juliadynamics.github.io/DynamicalSystems.jl/dev/ds/predefined/#DynamicalSystemsBase.Systems.henon) (2D discrete dynamical system):
 ```julia
 using DynamicalSystems
-ds = Systems.henon()
+ds = PredefinedDynamicalSystems.henon()
 traj = trajectory(ds,10000;Ttr = 1000) #generate timeseries
 ```
 This can be fed into `timeseries_to_grid` which discretizes the timeseries and returns the name of the vertices:
@@ -102,9 +102,9 @@ entropy_measures = zeros(length(a_values))
 lyapunov_measures = zeros(length(a_values))
 
 for (i,a) in enumerate(a_values)
-    system = Systems.henon([0.0, 0.0]; a=a, b=b)
+    system = PredefinedDynamicalSystems.henon([0.0, 0.0]; a=a, b=b)
     @show a
-    timeseries = trajectory(system, traj_length, [0, 0]; Ttr=trans)
+    timeseries,  = trajectory(system, traj_length, [0, 0]; Ttr=trans)
     discrete_timeseries, vertex_names = timeseries_to_grid(timeseries, grid_size);
     stn,retcode = create_stn(discrete_timeseries, vertex_names)
     entropy_measures[i], lyapunov_measures[i] = network_measures(stn,ensemble, N_steps)
@@ -116,7 +116,7 @@ The workflow is similar for continuous systems. The difference is that one must
 
 Constructing state-transition network for the [Lorenz system](https://juliadynamics.github.io/DynamicalSystems.jl/dev/ds/predefined/#DynamicalSystemsBase.Systems.lorenz) (3D continuous dynamical system):
 ```julia
-ds = Systems.lorenz()
+ds = PredefinedDynamicalSystems.lorenz()
 plane = (1,15.0) #plane in 3D phase space with x = 15.0
 ```
 Calculate the `PSOS` with `poincaresos`
@@ -145,9 +145,9 @@ plane = (1,15.0) #plane in 3D phase space with x = 15.0
 lyap_measures = zeros(length(rho_values))
 entropy_measures = zeros(length(rho_values))
 for (i,ρ) in enumerate(rho_values)
-    ds = Systems.lorenz(ρ=ρ)
+    ds = PredefinedDynamicalSystems.lorenz(ρ=ρ)
     @show ρ
-    psection = ChaosTools.poincaresos(ds, plane, T; Ttr=500, idxs = [2,3],direction=+1,rootkw = (xrtol = 1e-8, atol = 1e-8))
+    psection = DynamicalSystemsBase.poincaresos(ds, plane, T; Ttr=500, idxs = [2,3],direction=+1,rootkw = (xrtol = 1e-8, atol = 1e-8))
     discrete_timeseries, vertex_names = timeseries_to_grid(psection, grid_size)
     stn,retcode = create_stn(discrete_timeseries, vertex_names)
     entropy_measures[i], lyap_measures[i] = network_measures(stn,ensemble, N_steps)

src/StateTransitionNetworks.jl

@@ -5,8 +5,7 @@ import StatsBase: sample, mean, var, Weights
 using MetaGraphsNext
 using GraphPlot,Cairo,Compose
 import SparseArrays: spzeros 
-import DelayEmbeddings: Dataset,AbstractDataset
-import ChaosTools: poincaresos
+using DynamicalSystemsBase: DynamicalSystemsBase,poincaresos,StateSpaceSet
 import LinearAlgebra: eigen, Diagonal, I, nullspace
 import DataStructures: OrderedDict
 

src/create_stn.jl

@@ -30,7 +30,7 @@ function timeseries_to_grid(timeseries, grid; grid_edges = [])
     vertex_names = OrderedDict{Tuple{Int64,Int64},Int64}()
     x_n, y_n = [], [];
 
-    for row in eachrow(timeseries)
+    for row in timeseries
         y = floor(Int,(row[2]-y_min)/Float64(y_grid.step))+1
         x = floor(Int,(row[1]-x_min)/Float64(x_grid.step))+1
         if M[y,x] == 0
@@ -47,18 +47,19 @@ end
 
 """
 	create_stn(ts,grid::Int64,plane,idxs;make_ergodic=false, verbose=false,kwargs...) -> stn
-Higher level function that creates a state transition network (STN) directly from the timeseries of a continuous system `ts` using `ChaosTools.poincaresos`. The returned `stn` is a directed metagraph object. \\
-`plane` is the same as in `ChaosTools.poincaresos`. With `idxs` you can choose the variables you want to save.
+Higher level function that creates a state transition network (STN) directly from the timeseries of a continuous system `ts` using `DynamicalSystemsBase.poincaresos`. The returned `stn` is a directed metagraph object. \\
+`plane` is the same as in `DynamicalSystemsBase.poincaresos`. With `idxs` you can choose the variables you want to save.
 `make_ergodic=true` returns an STN with only one strongly connected components. Defaults to `false`.   
 `verbose=true` logs info about the network checking process. Defaults to `false`.
-Other `kwargs` get propageted into `ChaosTools.poincaresos`.
+Other `kwargs` get propageted into `DynamicalSystemsBase.poincaresos`.
 
-For more info about PSOS  go to https://juliadynamics.github.io/DynamicalSystems.jl/v1.3/chaos/orbitdiagram/#ChaosTools.poincaresos
+For more info about PSOS  go to https://juliadynamics.github.io/DynamicalSystems.jl/v1.3/chaos/orbitdiagram/#DynamicalSystemsBase.poincaresos
 """
 function create_stn(ts,grid::Int64,plane,idxs;make_ergodic=false, verbose=false,kwargs...)
 
 	#psos
-	psection = poincaresos(Dataset(ts), plane; idxs=idxs,warning=true,kwargs...)
+	psection = DynamicalSystemsBase.poincaresos(DynamicalSystemsBase.StateSpaceSet(ts), plane; save_idxs=idxs,warning=true,kwargs...)
+
 	#method for time-discrete trajectory
 	create_stn(psection,grid; make_ergodic=make_ergodic, verbose=verbose)
 

src/timeseries_analysis.jl

@@ -1,7 +1,9 @@
-stn_analysis(timeseries::Matrix;grid,plane,idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false,use_stored_distribution=false) = stn_analysis(Dataset(timeseries);grid=grid,plane=plane,idxs=idxs,ensemble=ensemble,N_steps=N_steps,make_ergodic=make_ergodic, verbose=verbose,return_stn=return_stn,use_analytic=use_analytic,use_stored_distribution=use_stored_distribution)
+stn_analysis(timeseries::Matrix;grid,plane,idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false,use_stored_distribution=false) = stn_analysis(DynamicalSystemsBase.StateSpaceSet(timeseries);grid=grid,plane=plane,idxs=idxs,ensemble=ensemble,N_steps=N_steps,make_ergodic=make_ergodic, verbose=verbose,return_stn=return_stn,use_analytic=use_analytic,use_stored_distribution=use_stored_distribution)
+
+
 
 """
-	stn_analysis(timeseries::Dataset;grid,plane,idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false) ->  S, Λ
+	stn_analysis(timeseries::DynamicalSystemsBase.StateSpaceSet;grid,plane,idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false) ->  S, Λ
 	
 Calculates the Sinai-Kolmogorov Entropy and Lyapunov measure of a STN created from `timeseries`. If `retcode` is other than `:Success`, `NaN`s are returned.
 ## Keyword arguments
@@ -13,10 +15,13 @@ Calculates the Sinai-Kolmogorov Entropy and Lyapunov measure of a STN created fr
 * `return_stn` : returns the `stn` and the `retcode` as well
 * `make_ergodic` : returns an `stn` with only one strongly connected component. Defaults to `false`.  
 * `verbose` : logs the connectedness checking process
+* `use_analytic` : use analytic formula for calculating network measures
+* `use_stored_distribution` : use the already stored stationary distribution
 """
-function stn_analysis(timeseries::Dataset;grid,plane,idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false,use_stored_distribution=false)
 
-	stn,retcode = create_stn(timeseries,grid::Int64,plane,idxs;make_ergodic=make_ergodic, verbose=verbose)
+function stn_analysis(timeseries::DynamicalSystemsBase.StateSpaceSet;grid,plane,idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false,use_stored_distribution=false)
+
+	stn,retcode = create_stn(timeseries,grid::Int64,plane,save_idxs;make_ergodic=make_ergodic, verbose=verbose)
 
 	entropy,lyapunov = NaN,NaN #initialize variables
 

tests/adding_stn.jl

@@ -95,16 +95,16 @@ T = 500;
 ρ = 180.1;
 
 u0 = rand(Float64,3).*50 .-25';
-system = Systems.lorenz(u0; ρ=ρ);
+system = PredefinedDynamicalSystems.lorenz(u0; ρ=ρ);
 timeseries_1 = trajectory(system, T; Δt=Δt, Ttr=1000);
-psection_1 = ChaosTools.poincaresos(timeseries_1, plane; direction=+1, idxs=[2,3]);
+psection_1 = DynamicalSystemsBase.poincaresos(timeseries_1, plane; direction=+1, save_idxs=[2,3]);
 dtraj_1, v1 = timeseries_to_grid(psection_1, grid);
 stn1, ret_code_stn = create_stn(dtraj_1, v1; make_ergodic=true,verbose=false);
 
 u0 = rand(Float64,3).*50 .-25';
-system = Systems.lorenz(u0; ρ=ρ);
+system = PredefinedDynamicalSystems.lorenz(u0; ρ=ρ);
 timeseries_2 = trajectory(system, T; Δt=Δt, Ttr=1000);
-psection_2 = ChaosTools.poincaresos(timeseries_2, plane; direction=+1, idxs=[2,3]);
+psection_2 = DynamicalSystemsBase.poincaresos(timeseries_2, plane; direction=+1, save_idxs=[2,3]);
 dtraj_2, v2 = timeseries_to_grid(psection_2, grid);
 stn2, ret_code_stn = create_stn(dtraj_2, v2; make_ergodic=true,verbose=false);
 

tests/plot_convergence_henon.jl

@@ -3,7 +3,7 @@ using StateTransitionNetworks
 using Plots
 using LaTeXStrings
 
-ds = Systems.henon()
+ds = PredefinedDynamicalSystems.henon()
 params = [1.2,1.2265,1.24,1.27]
 ensemble = 10000
 N_max = 5000

tests/plot_convergence_lorenz.jl

@@ -4,7 +4,7 @@ using Plots
 using LaTeXStrings
 using OrdinaryDiffEq	
 
-ds = Systems.lorenz()
+ds = PredefinedDynamicalSystems.lorenz()
 tol = 1e-8
 params = [180.7,180.1,180.78,181.1]
 ensemble = 10000
@@ -23,7 +23,7 @@ pl = plot()
 for (i,ρ) in enumerate(params)
 	set_parameter!(ds,2,ρ)
 	@show ρ
-	psection = ChaosTools.poincaresos(ds, plane, T; Ttr=Ttr, direction=+1,rootkw = (xrtol = tol, atol = tol));
+	psection = DynamicalSystemsBase.poincaresos(ds, plane, T; Ttr=Ttr, direction=+1,rootkw = (xrtol = tol, atol = tol));
 	traj = psection[:,2:end] 
 	traj_grid, vertex_names = timeseries_to_grid(traj,grid_size) 
 	stn,retcode = create_stn(traj_grid,vertex_names)

tests/plot_functions.jl

@@ -39,7 +39,7 @@ end
 
 function plot_psection(timeseries::Matrix;plane,idxs,color)
 	dim = size(timeseries)[2]
-	psection = poincaresos(Dataset(timeseries), plane;idxs=idxs); 
+	psection = poincaresos(Dataset(timeseries), plane;save_idxs=idxs); 
 	i = idxs[1]
 	j = idxs[2]
 

tests/plot_measures_henon.jl

@@ -93,7 +93,10 @@ savefig("./tests/henon_S-Lyap-gr_a=1.4_1.2265_gr=5-500_t=5x10^6_ttrans=1000_ens=
 b = 0.3;
 a = 1.4;
 a_values = 1:0.001:1.4;
-ds = Systems.henon([0.1, 0.123]; a=a, b=b)
+
+ds = PredefinedDynamicalSystems.henon([0.0, 0.0]; a=a, b=b)
+#ds = PredefinedDynamicalSystems.henon([0.1, 0.123]; a=a, b=b)
+
 i = 1
 ics = [rand() for m in 1:10]
 n = 500
@@ -123,7 +126,7 @@ dpi=300)
 
 sim_lyapunov_eponent = zeros(length(a_values))
 for (i,a) in enumerate(a_values)
-    henon = Systems.henon([0.0, 0.0]; a=a, b=b)
+    henon = PredefinedDynamicalSystems.henon([0.0, 0.0]; a=a, b=b)
     @show a
     λ = lyapunov(henon, 10000; d0 = 1e-7, threshold = 1e-4, Ttr = 500)
     sim_lyapunov_eponent[i] = λ
@@ -164,9 +167,9 @@ theor_lyapunov_measures = zeros(length(a_values))
 #a = 1.4
 #i = 1
 for (i,a) in enumerate(a_values)
-    system = Systems.henon([0.0, 0.0]; a=a, b=b)
+    system = PredefinedDynamicalSystems.henon([0.0, 0.0]; a=a, b=b)
     @show a
-    timeseries = trajectory(system, traj_length, [0, 0]; Ttr=trans)
+    timeseries,  = trajectory(system, traj_length, [0, 0]; Ttr=trans)
     discrete_timeseries, vertex_names = timeseries_to_grid(timeseries, grid_size);
     stn,retcode = create_stn(discrete_timeseries, vertex_names)
     sim_entropy_measures[i], sim_lyapunov_measures[i] = network_measures(stn, ensemble, N_steps)
@@ -277,7 +280,7 @@ savefig(pl, "./tests/henon_scaling_ac=1.2266.png")
 #############
 ### intermittency plot
 #############
-ds = Systems.henon([0.01, 0.20]; a=1.2265, b=0.3)
+ds = PredefinedDynamicalSystems.henon([0.01, 0.20]; a=1.2265, b=0.3)
 traj = trajectory(ds,10000;Ttr = 1000) #generate timeseries
 pl = plot()
 plot!(pl, 1:10001, traj[:,1], label=nothing, lw=2, color="red")

tests/plot_measures_lorenz.jl

@@ -40,7 +40,7 @@ plane = (1,15.0);
 u0 = rand(3)
 rho_values = 180.:0.001:182.;
 T = 20000;
-ds = Systems.lorenz(u0)
+ds = PredefinedDynamicalSystems.lorenz(u0)
 Ttr = 19990
 
 # computation
@@ -72,7 +72,7 @@ savefig(pl, "./tests/lorenz_bif.png")
 
 sim_lyapunov_eponent = zeros(length(rho_values))
 for (i,ρ) in enumerate(rho_values)
-    ds = Systems.lorenz(; ρ=ρ)
+    ds = PredefinedDynamicalSystems.lorenz(; ρ=ρ)
     @show ρ
     λ = lyapunov(ds, 10000; d0=1e-7, Ttr=5000)
     sim_lyapunov_eponent[i] = λ
@@ -119,9 +119,9 @@ for (i,ρ) in enumerate(rho_values)
         @show ρ
         while true # STN must be 'healthy'
             u0 = rand(Float64,3).*50 .-25;
-            system = Systems.lorenz(u0; ρ=ρ);
-            timeseries = trajectory(system, T; Δt=Δt, Ttr=trans);
-            psection = ChaosTools.poincaresos(timeseries, plane; direction=+1, idxs=[2,3]);
+            system = PredefinedDynamicalSystems.lorenz(u0; ρ=ρ);
+            timeseries,  = trajectory(system, T; Δt=Δt, Ttr=trans);
+            psection = DynamicalSystemsBase.poincaresos(timeseries, plane; direction=+1, save_idxs=[2,3]);
             d_traj, v_names = timeseries_to_grid(psection, grid);
             stn, ret_code_stn = create_stn(d_traj, v_names);
             if ret_code_stn == :Success

tests/roessler_duffingmap/measures_od_duffingmap.jl

@@ -1,15 +1,13 @@
 using StateTransitionNetworks
-using Plots
 using DynamicalSystems
 using Graphs
-using LaTeXStrings
 using StatsBase
 using DelimitedFiles
 
 
 T = 30000
 Ttr = 1000
-a_vals = [2.505:0.000005:2.51;]
+a_vals = [2.505:0.0005:2.51;]
 
 #stn_colors = Dict(zip([180.10,180.70,180.78,181.10],[colorant"orange",colorant"red",colorant"green",colorant"blue"]))
 
@@ -18,7 +16,7 @@ duffingmap(u,p,n) = SVector(u[2],-p[2]*u[1] +p[1]*u[2] - u[2]^3)
 p = [2.3,0.1]
 u0 = [0.5,0.2]
 
-ds = DiscreteDynamicalSystem(duffingmap,u0,p);
+ds = DeterministicIteratedMap(duffingmap,u0,p);
 grid_size = 20
 rw_ensemble = 100 
 nr_steps = 10000
@@ -29,7 +27,7 @@ nr_steps = 10000
 
 println("Orbit diagram calculation starting...")
 od = orbitdiagram(ds, 1, 1, a_vals; n=1000, Ttr=Ttr)
-writedlm("orbit_diagram_duffingmap_a_saved_x_T_$T"*"_Ttr_$Ttr.txt",od)
+writedlm("data/orbit_diagram_duffingmap_a_saved_x_T_$T"*"_Ttr_$Ttr.txt",od)
 println("Done.")
 
 #---------------normal lyapunovs--------------------
@@ -46,7 +44,7 @@ for (i,a) in enumerate(a_vals)
 
 end
 
-writedlm("lyapunov_exponent_duffingmap_a_T_$T"*"_Ttr_$Ttr.txt",hcat(a_vals,lyapunov_exponents))
+writedlm("data/lyapunov_exponent_duffingmap_a_T_$T"*"_Ttr_$Ttr.txt",hcat(a_vals,lyapunov_exponents))
 println("Done.")
 
 #----------------calculate network measures-----------------
@@ -63,9 +61,8 @@ for a in a_vals
 
 	set_parameter!(ds,1,a)
 
-	timeseries = trajectory(ds, T;Ttr=Ttr);
-	d_traj, v_names = timeseries_to_grid(timeseries, grid_size);
-	stn, retcode = create_stn(d_traj, v_names;make_ergodic=true,verbose=true);
+	timeseries,  = trajectory(ds, T;Ttr=Ttr);
+	stn, retcode = create_stn(timeseries,grid_size;make_ergodic=true,verbose=true);
 
 	if retcode == :Success
 
@@ -95,7 +92,7 @@ for a in a_vals
 
 end
 
-writedlm("data_netmeasures_duffingmap_b_T_$T"*"_Ttr_$Ttr"*"_rwens_$rw_ensemble"*"_nsteps_$nr_steps"*"_grid_$grid_size"*".txt",hcat(a_vals,analytic_entropies,analytic_lyapunovs,num_entropies,num_lyapunovs))
+writedlm("data/data_netmeasures_duffingmap_b_T_$T"*"_Ttr_$Ttr"*"_rwens_$rw_ensemble"*"_nsteps_$nr_steps"*"_grid_$grid_size"*".txt",hcat(a_vals,analytic_entropies,analytic_lyapunovs,num_entropies,num_lyapunovs))
 
 
 

tests/roessler_duffingmap/measures_od_roessler.jl

@@ -1,22 +1,21 @@
 using StateTransitionNetworks
-using Plots
 using DynamicalSystems
 using Graphs
-using LaTeXStrings
 using StatsBase
 using DelimitedFiles
 
-
 T = 5000
+nr_points = 5000
 Ttr = 500
-b_vals = [0.2:0.001:1.0;]
+b_vals = [0.2:0.01:1.0;]
 
 #stn_colors = Dict(zip([180.10,180.70,180.78,181.10],[colorant"orange",colorant"red",colorant"green",colorant"blue"]))
 
 
 
 plane = (2,0.0)
-ds = Systems.roessler();
+ds = PredefinedDynamicalSystems.roessler();
+pmap = PoincareMap(ds,plane;rootkw=(xrtol=1e-10, atol=1e-10), direction=1)
 grid_size = 20
 Δt = 0.01
 rw_ensemble = 100 
@@ -26,8 +25,8 @@ nr_steps = 10000
 #---------------orbit diagram-----------------------
 
 println("Orbit diagram calculation starting...")
-od = produce_orbitdiagram(ds, plane, 1, 2, b_vals; tfinal=T, Ttr=Ttr, direction=1, printparams=true, rootkw=(xrtol=1e-10, atol=1e-10))
-writedlm("orbit_diagram_roessler_b_saved_z_T_$T"*"_Ttr_$Ttr.txt",od)
+od = orbitdiagram(pmap, 1, 2, b_vals; n=nr_points, Ttr=Ttr, show_progress=true)
+writedlm("data/orbit_diagram_roessler_b_saved_z_T_$T"*"_Ttr_$Ttr.txt",od)
 println("Done.")
 
 #---------------normal lyapunovs--------------------
@@ -44,7 +43,7 @@ for (i,b) in enumerate(b_vals)
 
 end
 
-writedlm("lyapunov_exponent_roessler_b_T_$T"*"_Ttr_$Ttr.txt",hcat(b_vals,lyapunov_exponents))
+writedlm("data/lyapunov_exponent_roessler_b_T_$T"*"_Ttr_$Ttr.txt",hcat(b_vals,lyapunov_exponents))
 println("Done.")
 
 
@@ -63,10 +62,8 @@ for b in b_vals
 
 	set_parameter!(ds,2,b)
 
-	timeseries = trajectory(ds, T; Δt=Δt, Ttr=Ttr);
-	psection = ChaosTools.poincaresos(timeseries, plane; direction=+1, idxs=[1,3]);
-	d_traj, v_names = timeseries_to_grid(psection, grid_size);
-	stn, retcode = create_stn(d_traj, v_names;make_ergodic=true,verbose=true);
+	timeseries,  = trajectory(ds, T; Δt=Δt, Ttr=Ttr);
+	stn, retcode = create_stn(timeseries, grid_size, plane, [1,3];make_ergodic=true,verbose=true,direction=+1);
 
 	if retcode == :Success