diff --git a/Project.toml b/Project.toml
index cebc5c0..1e558f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 Cairo = "159f3aea-2a34-519c-b102-8c37f9878175"
 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"
@@ -17,6 +18,7 @@ MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377"
 Permutations = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
+
+
diff --git a/src/StateTransitionNetworks.jl b/src/StateTransitionNetworks.jl
index a4767d7..4ab24af 100644
--- a/src/StateTransitionNetworks.jl
+++ b/src/StateTransitionNetworks.jl
@@ -1,21 +1,24 @@
 module StateTransitionNetworks
 
-import Graphs: DiGraph,add_edge!, inneighbors,outneighbors, nv, ne, weights, is_strongly_connected,strongly_connected_components,edges,degree_histogram
+import Graphs: DiGraph,add_edge!, inneighbors,outneighbors, nv, ne, weights, is_strongly_connected,strongly_connected_components,edges,degree_histogram,vertices
 import StatsBase: sample, mean, var, Weights
 using MetaGraphsNext
 using GraphPlot,Cairo,Compose
 import SparseArrays: spzeros 
 using DynamicalSystemsBase: DynamicalSystemsBase,poincaresos,StateSpaceSet
-import LinearAlgebra: eigen, Diagonal
+import LinearAlgebra: eigen, Diagonal, I, nullspace
+import DataStructures: OrderedDict
 
 include("create_stn.jl")
 include("network_measures.jl")
 include("plot_stn.jl")
 include("timeseries_analysis.jl")
+include("operations.jl")
 
-export timeseries_to_grid, create_stn, check_stn!, prob_matrix, weight_matrix, calculate_weight_matrix, random_walk_on_stn, randomwalk_step,isnormalized
+export timeseries_to_grid, create_stn, check_stn!, prob_matrix, weight_matrix, state_distribution, calculate_weight_matrix, random_walk_on_stn, randomwalk_step,isnormalized
 export network_measures, sinai_kolmogorov_entropy, measure_convergence,lyapunov_measure
 export plot_stn
 export stn_analysis,read_bin,ndensity
+export add_timeseries,are_equal,timeseries_to_common_grid
 
 end
diff --git a/src/create_stn.jl b/src/create_stn.jl
index ef87ef4..90a3387 100644
--- a/src/create_stn.jl
+++ b/src/create_stn.jl
@@ -4,29 +4,30 @@ Discretizes a 2D timeseries/trajectory on a grid. Returns
 a discrete timeseries containing the the cell coordinates and the list
 of vertices with cell coordinates. 
 """
-function timeseries_to_grid(timeseries, grid)    
+function timeseries_to_grid(timeseries, grid; grid_edges = [])    
     M = zeros(grid,grid)
     T = length(timeseries[:,1])
-    x_min = minimum(timeseries[:, 1])
-    y_min = minimum(timeseries[:, 2])
-    x_max = maximum(timeseries[:, 1])
-    y_max = maximum(timeseries[:, 2])
     
-    #partitioning with little extra space on both ends
-    
-    dx = 0.5*(x_max-x_min)/grid
-    dy = 0.5*(y_max-y_min)/grid
+    if isempty(grid_edges)		
+		x_min = minimum(timeseries[:, 1])
+		y_min = minimum(timeseries[:, 2])
+		x_max = maximum(timeseries[:, 1])
+		y_max = maximum(timeseries[:, 2])
+    else
+    	x_min,y_min,x_max,y_max = grid_edges 
+    end
     
-    x_grid = range(x_min-dx, x_max+dx, grid);
+	#partitioning between [`x_min`, `x_max`] into `grid` number of cells
+    x_grid = range(x_min, nextfloat(x_max, 100), grid+1);
     x_min = x_grid[1]
-    y_grid = range(y_min-dy, y_max+dy, grid);
+    y_grid = range(y_min, nextfloat(y_max, 100), grid+1);
     y_min = y_grid[1]
     
     #arrays for space-discrete timeseries 
     x_n = Vector{Int64}(undef, T)
     y_n = Vector{Int64}(undef, T)
     num_vertex = 0
-    vertex_names = [];
+    vertex_names = OrderedDict{Tuple{Int64,Int64},Int64}()
     x_n, y_n = [], [];
 
     for row in timeseries
@@ -34,18 +35,16 @@ function timeseries_to_grid(timeseries, grid)
         x = floor(Int,(row[1]-x_min)/Float64(x_grid.step))+1
         if M[y,x] == 0
             num_vertex += 1 
-            push!(vertex_names, [num_vertex, x, y])
+            vertex_names[(x,y)] = num_vertex
             M[y,x] = 1
         end
         push!(x_n, x)
         push!(y_n, y)
     end
-    vertex_names = reduce(hcat, vertex_names)'
-    d_timeseries = [x_n y_n]
+    d_timeseries = collect(zip(x_n,y_n))
     return d_timeseries, vertex_names
 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 `DynamicalSystemsBase.poincaresos`. The returned `stn` is a directed metagraph object. \\
@@ -60,6 +59,7 @@ function create_stn(ts,grid::Int64,plane,idxs;make_ergodic=false, verbose=false,
 	
 	#psos
 	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)
 	
@@ -90,23 +90,29 @@ Creates a state transition network (STN) using the discrete timeseries and verte
 	stn[i,j] -> (:prob => P[i,j],:weight => Q[i,j])
 """
 function create_stn(discrete_timeseries,vertex_names;make_ergodic=false,verbose=false)
-	nr_vertices = length(vertex_names[:,1])
-	states = zeros(Int32,length(discrete_timeseries[:,1])) #discrete timeseries but only with vertex indices
-	
-	for v in 1:nr_vertices
-		timepoints_at_v = findall(x -> x == vertex_names[v,2:end],[eachrow(discrete_timeseries)...]) #find all times when state was v
-		states[timepoints_at_v] .= v  
-	end		
+	nr_vertices = length(vertex_names)
 	
 	#weight and transition probability matrices
 	Q = spzeros(nr_vertices, nr_vertices)
     P = spzeros(nr_vertices, nr_vertices)
     
+    #probability distribution of states
+    states_distrib = zeros(nr_vertices)
+    
     #count transitions
-    next_states = circshift(states,-1)
-    for i in eachindex(states[1:end-1])
-        Q[states[i],next_states[i]] += 1
+    for i in eachindex(discrete_timeseries[1:end-1])
+    	state = discrete_timeseries[i]
+    	next_state = discrete_timeseries[i+1]
+        Q[vertex_names[state],vertex_names[next_state]] += 1
+        states_distrib[vertex_names[state]] += 1
     end
+    
+    #update end (final) state 
+    end_state = discrete_timeseries[end]
+    states_distrib[vertex_names[end_state]] += 1
+    
+    #normalize state distribution
+	states_distrib = states_distrib ./ sum(states_distrib)
 
 	#normalize Q and fill P by normalizing rows
     Q = Q./sum(Q)
@@ -116,7 +122,7 @@ function create_stn(discrete_timeseries,vertex_names;make_ergodic=false,verbose=
 	stn = MetaGraph(
         DiGraph(),
         Int64,
-        Dict{Symbol, Int64},
+        Dict{Symbol, Union{Int64,Float64}},
         Dict{Symbol, Float64},
         nothing,
         edge_data -> 1.0,
@@ -126,13 +132,12 @@ function create_stn(discrete_timeseries,vertex_names;make_ergodic=false,verbose=
 	#add edges and properties
 	#Properties: vertices ->  Dict{Symbol,Int64}, edges -> Dict{Symbol,Float64}
 	
-	for v in 1:nr_vertices
-		x,y = vertex_names[v,2:end] 
-		stn[v] = Dict(:x => x,:y => y)
+	for state in keys(vertex_names)
+		stn[vertex_names[state]] = Dict(:x => state[1],:y => state[2],:prob => states_distrib[vertex_names[state]])
 	end
 	
-	for i in 1:nr_vertices
-        for j in 1:nr_vertices
+	for i in 1:length(vertex_names)
+        for j in 1:length(vertex_names)
             if P[i, j] != 0
 				stn[i,j] = Dict(:prob => P[i,j],:weight => Q[i,j])
             end
@@ -160,14 +165,16 @@ function create_stn(P::AbstractMatrix;make_ergodic=false,verbose=false)
 	stn = MetaGraph(
         DiGraph(),
         Int64,
-        Dict{Symbol, Int64},
+        Dict{Symbol, Union{Int64,Float64}},
         Dict{Symbol, Float64},
         nothing,
         edge_data -> 1.0,
         0.0)
 
+	# add dummy coordinates
 	for v in 1:nr_vertices
-		stn[v] = Dict{Symbol, Int64}()
+		#stn[v] = Dict{Symbol, Int64}()
+		stn[v] = Dict(:x => v, :y => 0, :prob => 0.)
 	end
 	
 	Q = calculate_weight_matrix(P)
@@ -185,6 +192,27 @@ function create_stn(P::AbstractMatrix;make_ergodic=false,verbose=false)
 	return stn,retcode	
 end
 
+
+"""
+	state_distribution(stn) -> prob_states,pos_states
+Returns the probability distribution of states and the positions of the vertices that correspond to these states.  
+"""
+function state_distribution(stn)
+	prob_states = zeros(nv(stn))	#probability of vertices(states)
+	pos_states = zeros(Int32,nv(stn),2)	#position of vertices
+	
+	for v in collect(vertices(stn))
+		# v is the code of the vertex
+		v_label = label_for(stn,v)
+		prob_states[v] = stn[v_label][:prob] 
+		pos_states[v,:] .= [stn[v_label][:x],stn[v_label][:y]] 
+	end
+	
+	return prob_states,pos_states
+end
+
+
+
 function prob_matrix(stn)
 	nr_vertices = nv(stn)
 	P = spzeros(Float64,(nr_vertices,nr_vertices))
diff --git a/src/network_measures.jl b/src/network_measures.jl
index 4277c2b..8a6fe37 100644
--- a/src/network_measures.jl
+++ b/src/network_measures.jl
@@ -46,9 +46,9 @@ end
 Calculates the Sinai-Kolmogorov Entropy and Lyapunov measure of a STN
 by using the analytical definitions of both quantities
 """
-function network_measures(P::AbstractMatrix)
-   	entropy, ret_code_entr = sinai_kolmogorov_entropy(P)
-    lyapunov, variance, covariance, ret_code_lyap = lyapunov_measure(P)
+function network_measures(P::AbstractMatrix;x=nothing)
+   	entropy, ret_code_entr = sinai_kolmogorov_entropy(P;x=x)
+    lyapunov, variance, covariance, ret_code_lyap = lyapunov_measure(P;x=x)
 	return entropy, lyapunov
 end
 
@@ -109,21 +109,19 @@ end
 Calculates analytically the Sinai-Kolmogorov entropy given the P transition probability matrix of the STN. 
 
 """
-function sinai_kolmogorov_entropy(P::AbstractMatrix)
-	λ, V = eigen(Matrix(P))
-	λ, X = eigen(transpose(Matrix(P)))
-	
-	if real(λ[end]) ≈ 1
-	   x = transpose(X[:,end]./sum(X[:,end]))
-	   v = V[:,end]./V[1,end]
-	else
-	   return -1, :StochasticMatrixError
+function sinai_kolmogorov_entropy(P::AbstractMatrix;x=nothing)
+		
+	if isnothing(x)
+		x = nullspace(Matrix(P' - I))
+		x = x ./sum(x)
 	end
+	
+	v = ones(length(x))
 
 	L = Matrix(-log.(P))
 	replace!(L, Inf=>0.0)
 	L = P.*L
-	entropy = x*L*v
+	entropy = (x'*L*v)[1]
 	return real(entropy), :Success
 end
 
@@ -133,28 +131,25 @@ end
 Calculates analytically the Lyapunov measure given the the P transition probability matrix of the STN. 
 
 """
-function lyapunov_measure(P::AbstractMatrix)
-	λ, V = eigen(Matrix(P))
-	λ, X = eigen(transpose(Matrix(P)))
+function lyapunov_measure(P::AbstractMatrix;x=nothing)
 	
-	if real(λ[end]) ≈ 1
-	   x = transpose(X[:,end]./sum(X[:,end]))
-	   v = V[:,end]./V[1,end]
-	else
-	   return -1, -1, -1, :StochasticMatrixError
+	if isnothing(x)
+		x = nullspace(Matrix(P' - I))
+		x = x ./sum(x)
 	end
- 
+	x = x'
+	v = ones(length(x))
+	
 	L = Matrix(-log.(P))
 	replace!(L, Inf=>0.0)
 	L2 = P.*L.^2
 	L = P.*L
-	I = Diagonal(ones(length(x)))
 	S = (I-v*x)+(P-v*x)*inv(I-P+v*x)
-	covariance = x*L*S*L*v
-	variance = x*L2*v-(x*L*v)^2
-	lyapunov =  variance + 2*real(covariance)
+	covariance = (x*L*S*L*v)[1]
+	variance = (x*L2*v)[1] -(x*L*v)[1]^2
+	lyapunov =  variance + 2*covariance
 	if imag(covariance) < 1.0e-3
-	   return real(lyapunov), real(variance), real(covariance), :Success
+	   return lyapunov, variance, covariance, :Success
 	else
 	   @show covariance
 	   return real(lyapunov), real(variance), real(covariance), :ComplexCovariancveWarning
diff --git a/src/operations.jl b/src/operations.jl
new file mode 100644
index 0000000..7f95720
--- /dev/null
+++ b/src/operations.jl
@@ -0,0 +1,87 @@
+function get_grid_edges(psections)
+    x_max, y_max = maximum(reduce(vcat, maximum.(Matrix.(psections); dims=1)); dims=1)
+    x_min, y_min = minimum(reduce(vcat, minimum.(Matrix.(psections); dims=1)); dims=1)
+    return [x_min, y_min, x_max, y_max]
+end
+
+function timeseries_to_common_grid(list_of_timeseries, grid)
+    grid_edges = get_grid_edges(list_of_timeseries)
+    
+    symbolic_ts_list = []
+    
+    for ts in list_of_timeseries
+    	sym_ts,vnames = timeseries_to_grid(ts, grid; grid_edges = grid_edges)
+    	push!(symbolic_ts_list,sym_ts)
+    end
+	return symbolic_ts_list 
+end
+
+
+
+function add_timeseries(symbolic_ts_list, grid; make_ergodic=false, verbose=false)
+    Q_null = zeros(Int32, grid*grid, grid*grid) # Null matrix with all possible transitions
+    vertex_names = OrderedDict{Tuple{Int64,Int64},Int64}() 
+    vertex_place = [];  # Rows, and column number in Q_null for a given vertex
+    M = zeros(grid,grid)
+    nr_vertices = 0;
+    for ts in symbolic_ts_list
+        O = zeros(Int32, grid*grid, grid*grid)
+        
+        
+        for (i,state) in enumerate(ts)
+            x, y = state
+            v = (x-1)*grid + y
+            
+            if M[y,x] == 0
+                nr_vertices += 1 
+                vertex_names[(x,y)] = nr_vertices
+                push!(vertex_place, v)
+                M[y,x] = 1
+            end
+            
+            if i < length(ts) 
+            	x_next,y_next = ts[i+1]
+            	v_next = (x_next-1)*grid + y_next
+            	O[v,v_next] +=1 
+            end
+            
+        end
+        
+        Q_null = Q_null + O
+    end
+    Q = Q_null[vertex_place, vertex_place]
+	#weight and transition probability matrices
+    P = spzeros(nr_vertices, nr_vertices)
+	#normalize Q and fill P by normalizing rows
+    Q = Q./sum(Q)
+
+	P = renormalize(Q)
+	#create directed metagraph with static label and metadata types and default weight 0
+	stn, ret_code = create_stn(P; make_ergodic=make_ergodic, verbose=verbose)
+	
+  	for state in keys(vertex_names)
+		stn[vertex_names[state]] = Dict(:x => state[1],:y => state[2])
+	end
+	
+    return stn, ret_code
+end
+
+function are_equal(stn1, stn2)
+    Q1 = weight_matrix(stn1)
+    Q2 = weight_matrix(stn2)
+    if size(Q1) != size(Q2)
+        return false
+    end
+    vertices_2 = Vector(vertices(stn2));
+    permutation = [];
+    for i in vertices(stn1)
+        for j in vertices_2
+            if stn1[i] == stn2[j]
+                push!(permutation, j)
+                deleteat!(vertices_2, findall(x->x==j,vertices_2))
+                break
+            end
+        end
+    end
+    return Q1 ≈ Q2[permutation, permutation]
+end
diff --git a/src/plot_stn.jl b/src/plot_stn.jl
index fe88fb2..aaaa6a9 100644
--- a/src/plot_stn.jl
+++ b/src/plot_stn.jl
@@ -1,15 +1,10 @@
 function plot_stn(stn;filename="stn.pdf",nodesize=1,nodefillc="orange",linetype="straight",max_edgelinewidth=1,weight_exponent=1,nodelabels=false,kwargs...)
-	nr_vertices = nv(stn)
-	x = []
-	y = []
-		
-	for i in 1:nr_vertices
-		ilabel = label_for(stn,i)
-		push!(x,stn[ilabel][:x])
-		push!(y,stn[ilabel][:y])
-	end
-	x = Int32.(x)
-	y = Int32.(y)
+	
+	prob_states,pos_states = state_distribution(stn)
+	#vertex alphas could be a function of their probability
+	
+	x = pos_states[:,1]
+	y = pos_states[:,2]
 	
 	w = weight_matrix(stn)
 	
diff --git a/src/timeseries_analysis.jl b/src/timeseries_analysis.jl
index c7a4027..5d8032b 100644
--- a/src/timeseries_analysis.jl
+++ b/src/timeseries_analysis.jl
@@ -1,4 +1,6 @@
-stn_analysis(timeseries::Matrix;grid,plane,save_idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false) = stn_analysis(DynamicalSystemsBase.StateSpaceSet(timeseries);grid=grid,plane=plane,save_idxs=save_idxs,ensemble=ensemble,N_steps=N_steps,make_ergodic=make_ergodic, verbose=verbose,return_stn=return_stn,use_analytic=use_analytic)
+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::DynamicalSystemsBase.StateSpaceSet;grid,plane,idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false) ->  S, Λ
@@ -13,8 +15,11 @@ 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::DynamicalSystemsBase.StateSpaceSet;grid,plane,save_idxs,ensemble=100,N_steps=1000,make_ergodic=false, verbose=false,return_stn=false,use_analytic=false)
+
+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)
 	
@@ -23,7 +28,7 @@ function stn_analysis(timeseries::DynamicalSystemsBase.StateSpaceSet;grid,plane,
 	if retcode ==:Success
 		if use_analytic
 			P = prob_matrix(stn)
-			entropy,lyapunov = network_measures(P)
+			entropy,lyapunov = use_stored_distribution ? network_measures(P;x=state_distribution(stn)[1]) : network_measures(P)
 		else
 			entropy,lyapunov = network_measures(stn,ensemble,N_steps)
 		end
diff --git a/tests/plot_measures_henon.jl b/tests/plot_measures_henon.jl
index 09b518c..97c9b8a 100644
--- a/tests/plot_measures_henon.jl
+++ b/tests/plot_measures_henon.jl
@@ -10,14 +10,93 @@ using LinearAlgebra
 using DelimitedFiles
 using LaTeXStrings
 
+###############################
+### Measures for a single value
+###############################
+
+grid_size = 200;
+b = 0.3;
+a = 1.4;
+ds = Systems.henon([0.1, 0.123]; a=a, b=b);
+λ = lyapunov(ds, 10000; d0 = 1e-7, threshold = 1e-4, Ttr = 500);
+timeseries = trajectory(ds, 100000, [0, 0]; Ttr=1000);
+discrete_timeseries, vertex_names = timeseries_to_grid(timeseries, grid_size);
+@time stn,retcode = create_stn(discrete_timeseries, vertex_names);
+P = prob_matrix(stn);
+@time S, Λ = network_measures(P)
+
+###############################
+### Measures as a function of the grid
+###############################
+
+data0 = collect(10:1:500);
+data1 = zeros(length(data0)); 
+data2 = zeros(length(data0)); 
+b = 0.3;
+a = 1.4;
+a = 1.2265;
+ens = 5
+ds = Systems.henon([0.1, 0.123]; a=a, b=b);
+λ = lyapunov(ds, 10000; d0 = 1e-7, threshold = 1e-4, Ttr = 500)
+#plot(timeseries[end-1000:end,1])
+for i in 1:ens
+    @show i 
+    u0 = 0.2*rand(2)
+    timeseries = trajectory(ds, 5000000, u0; Ttr=1000);
+    for (g,grid_size) in enumerate(data0);
+        @show grid_size
+        discrete_timeseries, vertex_names = timeseries_to_grid(timeseries, grid_size);
+        @time stn,retcode = create_stn(discrete_timeseries, vertex_names);
+        P = prob_matrix(stn);
+        @time S, Λ = network_measures(P)
+        data1[g] += S
+        data2[g] += Λ
+    end
+end
+data1 ./= ens
+data2 ./= ens
+
+# save data
+f_name1 = "./tests/henon_S-gr_a=1.4_b=0.3_gr=5-500_dgr=5_t=5x10^6_ttrans=1000.dat"
+f_name1 = "./tests/henon_S-gr_a=1.2265_b=0.3_gr=5-500_dgr=5_t=5x10^6_ttrans=1000.dat"
+f_name1 = "./tests/henon_S-gr_a=1.2265_b=0.3_gr=10-500_dgr=1_t=5x10^6_ttrans=1000_ens=5.dat"
+writedlm(f_name1,data1)
+f_name2 = "./tests/henon_Lyap-gr_a=1.4_b=0.3_gr=5-500_dgr=5_t=5x10^6_ttrans=1000.dat"
+f_name2 = "./tests/henon_Lyap-gr_a=1.2265_b=0.3_gr=5-500_dgr=5_t=5x10^6_ttrans=1000.dat"
+f_name2 = "./tests/henon_Lyap-gr_a=1.2265_b=0.3_gr=10-500_dgr=1_t=5x10^6_ttrans=1000_ens=5.dat"
+writedlm(f_name2,data2)
+
+# load data
+f_name1 = "./tests/henon_S-gr_a=1.4_b=0.3_gr=5-500_dgr=5_t=5x10^6_ttrans=1000.dat"
+f_name2 = "./tests/henon_Lyap-gr_a=1.4_b=0.3_gr=5-500_dgr=5_t=5x10^6_ttrans=1000.dat"
+f_name1cr = "./tests/henon_S-gr_a=1.2265_b=0.3_gr=10-500_dgr=1_t=5x10^6_ttrans=1000_ens=5.dat"
+f_name2cr = "./tests/henon_Lyap-gr_a=1.2265_b=0.3_gr=10-500_dgr=1_t=5x10^6_ttrans=1000_ens=5.dat"
+data0 = collect(5:5:500);
+data0cr = collect(10:1:500);
+data1 = readdlm(f_name1);
+data2 = readdlm(f_name2);
+data1cr = readdlm(f_name1cr);
+data2cr = readdlm(f_name2cr);
+
+plot(data0, data1, lw=2, color=:orange, label=L"S, a=1.4")
+plot!(data0, data2, lw=2, color=:orange, linestyle=:dash, label=L"\Lambda, a=1.4")
+plot!(data0cr, data1cr, lw=2, color=:red, label=L"S, a=1.2265")
+plot!(data0cr, data2cr, lw=2, color=:red, linestyle=:dash, label=L"\Lambda, a=1.2265")
+plot!(xlabel="# of grid cells", ylabel=L"S, \Lambda", ylim=[0,3], xguidefontsize=22, yguidefontsize=22, tickfontsize=14, lw=2, legendfontsize=16, dpi=300)
+savefig("./tests/henon_S-Lyap-gr_a=1.4_1.2265_gr=5-500_t=5x10^6_ttrans=1000_ens=5.pdf")
+
+
 #################
 ### orbit diagram
 #################
 
 b = 0.3;
-a = 1.0;
+a = 1.4;
 a_values = 1:0.001:1.4;
+
 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
@@ -64,7 +143,7 @@ sim_lyapunov_eponent = readdlm(f_name)
 pl = plot()
 plot!(pl, a_values, zeros(length(a_values)), linestyle=:dash, color="black", label=nothing)
 plot!(pl, a_values, sim_lyapunov_eponent, label=nothing, lw=2)
-plot!(pl, xlabel=L"a", ylabel=L"\lambda", xlim=[1,1.4], xticks=1:0.1:1.4, xguidefontsize=22, yguidefontsize=22, tickfontsize=14, lw=2, fontfamily="Serif", legendfontsize=16,
+plot!(pl, xlabel=L"a", ylabel=L"\lambda", xlim=[1,1.4], xticks=1:0.1:1.4, 
 dpi=300)
 
 ####################
@@ -73,10 +152,10 @@ dpi=300)
 
 b = 0.3;
 a_values = 1:0.001:1.4;
-a_values = 1.21:0.0001:1.23;
-traj_length = 30000;
+#a_values = 1.21:0.0001:1.23;
+traj_length = 1000000;
 trans = 1000;
-grid_size = 20;
+grid_size = 200;
 ensemble = 1000;
 N_steps = 10000;
 
@@ -101,12 +180,15 @@ end
 data = hcat(a_values, sim_entropy_measures, sim_lyapunov_measures, theor_entropy_measures, theor_lyapunov_measures)
 
 #save data
-f_name = "./tests/henon_measures_a=1.0-1.4_da=0.001_t=10^5_ens=1000_tmax=30000_ttrans=1000.dat"
-f_name = "./tests/henon_measures_a=1.21-1.23_da=0.0001_t=10^5_ens=1000_tmax=30000_ttrans=1000.dat"
+#f_name = "./tests/henon_measures_a=1.0-1.4_da=0.001_t=10^5_ens=1000_tmax=30000_ttrans=1000.dat"
+#f_name = "./tests/henon_measures_a=1.21-1.23_da=0.0001_t=10^5_ens=1000_tmax=30000_ttrans=1000.dat"
+#f_name = "./tests/henon_measures_a=1.0-1.4_da=0.001_ens=1000_tmax=300000_ttrans=1000_grid=100.dat"
+f_name = "./tests/henon_measures_a=1.0-1.4_da=0.001_ens=1000_tmax=1000000_ttrans=1000_grid=200.dat"
 writedlm(f_name,data)
 # load data
 f_name = "./tests/henon_measures_a=1.0-1.4_da=0.001_t=10^5_ens=1000_tmax=30000_ttrans=1000.dat"
 f_name = "./tests/henon_measures_a=1.21-1.23_da=0.0001_t=10^5_ens=1000_tmax=30000_ttrans=1000.dat"
+f_name = "./tests/henon_measures_a=1.0-1.4_da=0.001_ens=1000_tmax=300000_ttrans=1000.dat"
 data = readdlm(f_name)
 
 pl = scatter(data[:,1], data[:,2], label = "Random walk", lw=2, color="gray", alpha=0.5)
diff --git a/tests/plot_measures_lorenz.jl b/tests/plot_measures_lorenz.jl
index bd8bb95..693a369 100644
--- a/tests/plot_measures_lorenz.jl
+++ b/tests/plot_measures_lorenz.jl
@@ -10,6 +10,29 @@ using LinearAlgebra
 using DelimitedFiles
 using LaTeXStrings
 
+
+###############################
+### Measures for a single value
+###############################
+
+Δt = 0.001;
+plane = (1,15.0);
+grid_size = 1000;
+u0 = rand(3)
+ρ = 180.
+T = 100000;
+trans = 5000
+
+
+system = Systems.lorenz(u0; ρ=ρ);
+timeseries = trajectory(system, T; Δt=Δt, Ttr=trans);
+psection = ChaosTools.poincaresos(timeseries, plane; direction=+1, idxs=[2,3]);
+d_traj, v_names = timeseries_to_grid(psection, grid_size);
+stn, ret_code_stn = create_stn(d_traj, v_names);
+P = prob_matrix(stn);
+network_measures(P)
+
+
 #################
 ### orbit diagram
 #################
@@ -204,4 +227,7 @@ savefig(pl, "./tests/lorenz_corr_lyapexp-lyapmes.png")
 pl = plot()
 scatter!(pl, (data_zoom[:,1] .- ρ_c)/ρ_c, data_zoom[:,5] ./ data_zoom[:,4], label=L"\rho_c=%$(ρ_c)", lw=0, ms=6, markerstrokewidth=0, color="red", alpha=0.5)
 plot!(pl, xaxis=:log, yaxis=:log, xlabel=L"(\rho-\rho_c)/\rho_c", ylabel=L"\Lambda/S", xlim=[1e-6,1e-3], xguidefontsize=22, yguidefontsize=22, tickfontsize=14, lw=2, fontfamily="sans-serif", legendfontsize=16, dpi=300)
-savefig(pl, "./tests/lorenz_scaling_rhoc=181.669.png")
\ No newline at end of file
+savefig(pl, "./tests/lorenz_scaling_rhoc=181.669.png")
+
+
+
diff --git a/tests/test_logistic_map.jl b/tests/test_logistic_map.jl
new file mode 100644
index 0000000..b5c3c24
--- /dev/null
+++ b/tests/test_logistic_map.jl
@@ -0,0 +1,175 @@
+using StateTransitionNetworks
+using Plots
+using DynamicalSystems
+using Graphs
+using GraphPlot
+using Random
+using Statistics
+
+using LinearAlgebra
+using DelimitedFiles
+using LaTeXStrings
+
+function transition_matrix(f, grid, grid_edges)
+    P = zeros(grid,grid)
+    x_min,x_max = grid_edges 
+    dx = 0.5*(x_max-x_min)/grid
+    x = collect(range(grid_edges[1],grid_edges[2], grid))
+    for i in 1:grid
+        P[i,Int(f(i))]
+    end
+end
+
+function bit_number_measures(P::AbstractMatrix)
+	λ, X = eigen(transpose(Matrix(P)))
+	
+	if real(λ[end]) ≈ 1
+	   x = transpose(X[:,end]./sum(X[:,end]))
+	else
+	   return -1, :StochasticMatrixError
+	end
+
+	l = -log.(x)
+	replace!(l, Inf=>0.0)
+	entropy = sum((x .* l))
+    #variance = sum(x .* ((l .- entropy).^2))
+    variance =  sum(x .* l .* l) - entropy^2
+	return real(entropy), real(variance), :Success
+end
+
+###############################
+### Measures for a single value
+###############################
+
+grid_size = 20;
+r = 4.;
+ds = Systems.logistic(0.4; r = r)
+λ = lyapunov(ds, 10000; Ttr = 500)
+timeseries = trajectory(ds, 10^8, 0.4; Ttr=1000);
+ys = zeros(length(timeseries))
+timeseries = hcat(timeseries, ys)
+discrete_timeseries, vertex_names = timeseries_to_grid(timeseries, grid_size; grid_edges = [0., 0., 1., 1.]);
+@time stn,retcode = create_stn(discrete_timeseries, vertex_names);
+P = prob_matrix(stn);
+@time S, Λ = network_measures(P)
+@time C1, C2 = bit_number_measures(P) 
+
+####################
+### network measures
+####################
+
+r_values = 3.995:0.001:4.;
+traj_length = 10^8;
+trans = 1000;
+grid_size = 20;
+ensemble = 1000;
+N_steps = 10000;
+
+sim_entropy_measures = zeros(length(r_values))
+sim_lyapunov_measures = zeros(length(r_values))
+theor_entropy_measures = zeros(length(r_values))
+theor_lyapunov_measures = zeros(length(r_values))
+C1_measures = zeros(length(r_values))
+C2_measures = zeros(length(r_values))
+
+for (i,r) in enumerate(r_values)
+    system = Systems.logistic(0.4; r=r)
+    @show r
+    timeseries = trajectory(system, traj_length, 0.4; Ttr=trans)
+    ys = zeros(length(timeseries))
+    timeseries = hcat(timeseries, ys)
+    discrete_timeseries, vertex_names = timeseries_to_grid(timeseries, grid_size; grid_edges = [0., 0., 1., 1.]);
+    stn,retcode = create_stn(discrete_timeseries, vertex_names)
+    #sim_entropy_measures[i], sim_lyapunov_measures[i] = network_measures(stn, ensemble, N_steps)
+    P = prob_matrix(stn);
+    theor_entropy_measures[i], theor_lyapunov_measures[i] = network_measures(P)
+    #C1_measures[i], C2_measures[i] = bit_number_measures(P)
+end
+
+data = hcat(r_values, sim_entropy_measures, sim_lyapunov_measures, theor_entropy_measures, theor_lyapunov_measures, C1_measures, C2_measures)
+
+#save data
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.001_ens=1000_tmax=100000_ttrans=1000_grid=200.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=10^2_ttrans=1000_grid=200.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=10^3_ttrans=1000_grid=200.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=10^4_ttrans=1000_grid=200.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=10^5_ttrans=1000_grid=200.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=200.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=10^7_ttrans=1000_grid=200.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=1000.dat"
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=1000000_ttrans=1000_grid=2000.dat"
+writedlm(f_name,data)
+# load data
+f_name = "./tests/logistic_measures_r=0-4_da=0.01_ens=1000_tmax=100000_ttrans=1000_grid=200.dat"
+data = readdlm(f_name)
+
+# plot bit number measures C2
+pl = plot()
+plot!(pl, data[:,1], data[:,5], label = L"\Lambda", lw=2, alpha=0.5, color="red")
+plot!(pl, data[:,1], data[:,7], label = L"C_2", lw=2, alpha=0.5, color="orange")
+plot!(pl, xlabel=L"r", ylabel=L"\Lambda, C_2", xguidefontsize=22, yguidefontsize=22, xlim=[3.5,4], tickfontsize=14, lw=2, fontfamily="Serif", legendfontsize=16,
+dpi=300)
+plot!(ylim=[0,1], xlim=[3.6,3.7])
+savefig(pl, "./tests/logistic_var_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.svg")
+savefig(pl, "./tests/logistic_var_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.pdf")
+
+# plot bit number measures C1
+pl = plot()
+plot!(pl, data[:,1], data[:,4], label = L"S", lw=2, alpha=0.5, color="black")
+plot!(pl, data[:,1], data[:,6], label = L"C_1", lw=2, alpha=0.5, color="brown")
+plot!(pl, xlabel=L"r", ylabel=L"S, C_1", xguidefontsize=22, yguidefontsize=22, xlim=[3.5,4], tickfontsize=14, lw=2, fontfamily="Serif", legendfontsize=16,
+dpi=300)
+savefig(pl, "./tests/logistic_entr_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.svg")
+savefig(pl, "./tests/logistic_entr_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.pdf")
+
+pl = plot()
+#scatter!(pl, data[:,1], data[:,2], label = "Random walk", lw=2, color="gray", alpha=0.5)
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=1000.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,4], label = "grid=1000", lw=2, alpha=0.5, color="black")
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=200.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,4], label = "grid=200", lw=2, alpha=0.5, color="red")
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,4], label = "grid=20", lw=2, alpha=0.5, color="blue")
+plot!(pl, xlabel=L"r", ylabel=L"S", xguidefontsize=22, yguidefontsize=22, xlim=[3.5,4], tickfontsize=14, lw=2, fontfamily="Serif", legendfontsize=16,
+dpi=300)
+savefig(pl, "./tests/logistic_entr_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20-200-1000.svg")
+savefig(pl, "./tests/logistic_entr_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20-200-1000.pdf")
+
+
+# lyapunov measure
+pl = plot()
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=1000.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,5].+0.01, label = "grid=1000", lw=2, alpha=0.5, color="black")
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=200.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,5].+0.01, label = "grid=200", lw=2, alpha=0.5, color="red")
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,5].+0.01, label = "grid=20", lw=2, alpha=0.5, color="blue")
+plot!(pl, xlabel=L"r", ylabel=L"\Lambda", xguidefontsize=22, yguidefontsize=22, xlim=[3.5,4], tickfontsize=14, lw=2, fontfamily="Serif", legendfontsize=16,
+dpi=300, yscale=:linear)
+plot!(ylim=[0,2])
+savefig(pl, "./tests/logistic_lyap_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20-200-1000.svg")
+savefig(pl, "./tests/logistic_lyap_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20-200-1000.pdf")
+
+# compare gris in r\in [3.6,3.7]
+pl = plot()
+#scatter!(pl, data[:,1], data[:,3], label = "Random walk", lw=2, color="gray", alpha=0.5)
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=1000.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,5].+0.01, label = "grid=1000", lw=2, alpha=0.5, color="black")
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=200.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,5].+0.01, label = "grid=200", lw=2, alpha=0.5, color="red")
+f_name = "./tests/logistic_measures_r=3.5-4_da=0.0001_tmax=100000_ttrans=1000_grid=20.dat"
+data = readdlm(f_name)
+plot!(pl, data[:,1], data[:,5].+0.01, label = "grid=20", lw=2, alpha=0.5, color="blue")
+plot!(pl, xlabel=L"r", ylabel=L"\Lambda+0.01", xguidefontsize=22, yguidefontsize=22, xlim=[3.6,3.7], tickfontsize=14, lw=2, fontfamily="Serif", legendfontsize=16,
+dpi=300, yscale=:log)
+savefig(pl, "./tests/logistic_lyap_r=3.6-3.7_da=0.0001_tmax=100000_ttrans=1000_grid=20-200-1000.svg")
+savefig(pl, "./tests/logistic_lyap_r=3.6-3.7_da=0.0001_tmax=100000_ttrans=1000_grid=20-200-1000.pdf")
\ No newline at end of file
diff --git a/tests/test_random_stns.jl b/tests/test_random_stns.jl
new file mode 100644
index 0000000..b636117
--- /dev/null
+++ b/tests/test_random_stns.jl
@@ -0,0 +1,140 @@
+using StateTransitionNetworks
+using Plots
+using DynamicalSystems
+using Graphs
+using GraphPlot
+using Random
+using Statistics
+
+using LinearAlgebra
+using DelimitedFiles
+using LaTeXStrings
+
+N = 4
+ens = 1000
+
+
+function random_transition_matrix(N, A)
+    P = zeros(N,N)
+    for i in 1:N
+        for j in 1:N
+            if A[i,j]==1
+                P[i,j] = rand()
+            end
+        end
+    end
+    for i in 1:N
+        if sum(P[i,:]) !=0 
+            P[i,:] = P[i,:]/sum(P[i,:])
+        end
+    end
+    return P
+end
+
+a = 2^(N*N)-1
+A = Float64.(reshape(digits(a, base=2, pad=N*N), (N,N)))
+P = random_transition_matrix(N,A)
+network_measures(P)
+stn, retcode = create_stn(P; make_ergodic=true)
+P = prob_matrix(stn)
+network_measures(P)
+
+data0 = [];
+data1 = [];
+data2 = [];
+data3 = [];
+data4 = [];
+for a in 1:2^(N*N)-1
+    @show a
+    A = Float64.(reshape(digits(a, base=2, pad=N*N), (N,N)))
+    P = random_transition_matrix(N,A)
+    stn, retcode = create_stn(P; make_ergodic=true)
+    if retcode == :Success
+        av_entropy = 0.
+        av_lyapunov = 0.
+        var_entropy = 0.
+        var_lyapunov = 0.
+        for i in 1:ens
+            P = random_transition_matrix(N,A)
+            stn, retcode = create_stn(P; make_ergodic=true)
+            P = prob_matrix(stn)
+            entropy, lyapunov = network_measures(P)
+            av_entropy += entropy
+            av_lyapunov += lyapunov
+            var_entropy += entropy^2
+            var_lyapunov += lyapunov^2
+        end
+    else
+        av_entropy = -1. * ens
+        av_lyapunov = -1. * ens
+        var_entropy = 1. * ens
+        var_lyapunov = 1. * ens
+    end
+    push!(data0, a)
+    push!(data1, av_entropy/ens)
+    push!(data2, av_lyapunov/ens)
+    push!(data3, sqrt(var_entropy/ens-(av_entropy/ens)^2))
+    push!(data4, sqrt(var_lyapunov/ens-(av_lyapunov/ens)^2))
+end
+
+mask1 = data1 .!= -1
+plot(data0[mask1], data3[mask1], label=L"\mathrm{std}(S)")
+plot!(data0[mask1], data1[mask1], label=L"\langle S \rangle")
+plot!(xlabel=L"a", ylabel=L"S", xguidefontsize=22, yguidefontsize=22, tickfontsize=14, xtickfontsize=12, lw=2, legendfontsize=16)
+savefig("random_networks_entropy_N=4_uniform.pdf")
+mask2 = data2 .!= -1
+plot(data0[mask2], data4[mask2], label=L"\mathrm{std}(\Lambda)")
+plot!(data0[mask2], data2[mask2], label=L"\langle \Lambda \rangle")
+plot!(ylim=[0.01,2])
+plot!(xlabel=L"a", ylabel=L"\Lambda", xguidefontsize=22, yguidefontsize=22, tickfontsize=14, xtickfontsize=12, lw=2, legendfontsize=16)
+savefig("random_networks_lyapunov_N=4_uniform.pdf")
+
+
+idx = findall(x->x>1.15, data1)
+A = reshape(digits(idx[1], base=2, pad=N*N), (N,N))
+P = random_transition_matrix(N,A)
+stn, retcode = create_stn(P; make_ergodic=true)
+P = prob_matrix(stn)
+network_measures(P)
+plot_stn(stn)
+
+idx = findall(x->x>10, data4)
+
+# 0 lyapunov measure
+idx = findall(x->x≈0, data2)
+A = reshape(digits(idx[end], base=2, pad=N*N), (N,N))
+P = random_transition_matrix(N,A)
+stn, retcode = create_stn(P; make_ergodic=true)
+P = prob_matrix(stn)
+network_measures(P)
+
+idx = findmax(data1)[end]
+A = reshape(digits(idx[end], base=2, pad=N*N), (N,N))
+P = random_transition_matrix(N,A)
+stn, retcode = create_stn(P; make_ergodic=true)
+P = prob_matrix(stn)
+network_measures(P)
+
+
+# fully connected STN with random weights => S ~ N, Λ = 0.25
+data0 = [];
+data1 = [];
+data2 = [];
+for N in 1:500
+    @show N
+    A = ones(N,N)
+    P = random_transition_matrix(N,A)
+    S, Λ = network_measures(P)
+    push!(data0, N)
+    push!(data1, S)
+    push!(data2, Λ)
+    #lyapunov_measure(P)
+end
+mask1 = data1 .!= -1
+#mask1 = 1 ./ data1 .!= Inf
+plot(data0[mask1], data1[mask1], label=L"S", linewidth=4)
+plot!(data0[mask1], data2[mask1], label=L"\Lambda", linewidth=4)
+plot!(data0[mask1], log.(data0[mask1]), label=L"S=\log(N)", linewidth=2, linestyle= :dash,color=:black)
+plot!(data0[mask1], 0.25*ones(500), label=L"\Lambda=1/4", linewidth=2, linestyle= :dash,color=:gray)
+plot!(xlabel=L"N", ylabel=L"S,\Lambda", xguidefontsize=22, yguidefontsize=22, tickfontsize=14, xtickfontsize=12, lw=2, legendfontsize=16)
+#plot!(xaxis=:log, yaxis=:log)
\ No newline at end of file
diff --git a/unit_tests/operations.jl b/unit_tests/operations.jl
new file mode 100644
index 0000000..db621b0
--- /dev/null
+++ b/unit_tests/operations.jl
@@ -0,0 +1,26 @@
+using StateTransitionNetworks
+using DynamicalSystems
+using Test
+
+
+@testset "stn addition test: addition" begin
+	stn_added_test,retcode = add_timeseries([[(1,2),(2,3),(3,4),(1,2),(3,4)],[(3,4),(10,4),(1,2)]],10)
+	@test retcode == :Success
+	@test length(prob_matrix(stn_added_test).nzval) == 6
+end
+
+
+ds = Systems.henon()
+traj1 = trajectory(ds,10000;Ttr = 1000)
+traj2 = trajectory(ds,10000;Ttr = 11000)
+
+
+@testset "stn addition test: commutativity" begin
+	symts_list12 = timeseries_to_common_grid([traj1,traj2],20)
+	symts_list21 = timeseries_to_common_grid([traj2,traj1],20)
+	
+	stn_added12,retcode = add_timeseries(symts_list12,20)
+	stn_added21,retcode = add_timeseries(symts_list21,20)
+	@test are_equal(stn_added12,stn_added21)
+	
+end