# This file includes the simulation of five node network simulation 

import Cairo, Fontconfig
using GraphPlot
using Graphs
using SimpleWeightedGraphs
using LinearAlgebra
using DynamicalNetworks
using PyPlot
pygui(true)
close("all")

# Helper functions 
getwave(x, i, d) = getindex.(x, (i - 1) * d + 1)

# Node settings 
nodetype = Lorenz           # Node type 
P = diagm([1.0, 0.0, 0])    # State coupling matrix 
ηt = 10.0                   # Threshold for maximum Lyapuno exponents

# Topology 
n = 8
graph = SimpleWeightedGraph(cycle_graph(n))
add_edge!(graph, 2, 7, 1)
add_edge!(graph, 3, 6, 1)
gplot(graph, nodelabel = 1:nv(graph), layout = circular_layout) |> display

# Coupling strength 
W = laplacian_matrix(graph)
Q = [
    0 0 0 0 0 0 0 1
    0 0 0 0 0 0 1 0
    0 0 0 0 0 1 0 0
    0 0 0 0 1 0 0 0
    0 0 0 1 0 0 0 0
    0 0 1 0 0 0 0 0
    0 1 0 0 0 0 0 0
    1 0 0 0 0 0 0 0
]
M = eigvecs(Q)
G = inv(M) * (W ./ diag(W)) * M
B = G[1:4, 1:4]
D = G[5:8, 5:8]
λt = eigvals(B)
λs = eigvals(D)
λt[1] > λs[2] || error("Synchronization is not possible")
ϵr = ηt ./ [λt[1], λs[2]]
ϵ = ϵr[1] + abs(ϵr[2] - ϵr[1]) * 0.5

# Network 
σ = 0.0  # Noise strength
H = getnoisemat(graph) * σ
net = SDENetwork(nodetype, ϵ, graph, H, P)
sol = solvenet!(net, (0.0, 100.0), solverkwargs = (saveat = 0.001,))
t = sol.t
x = sol.u

# Plots 
d = dimension(nodetype)
subplot(221)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 2, d)))
xlabel("Time")
ylabel("e12")
subplot(222)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 8, d)))
xlabel("Time")
ylabel("e18")
subplot(223)
plot(t, getwave(x, 1, d))
xlabel("Time")
ylabel("x1")
subplot(224)
plot(getindex.(x, 1), getindex.(x, 2))
xlabel("x1")
ylabel("x2")
tight_layout()

figure()
plot(λt, zeros(length(λt)), marker = "*", label = L"\lambda_t")
plot(λs, ones(length(λs)), marker = "o", label = L"\lambda_s")
axvline.([λt[1], λs[2]], linestyle = "dashed", color = "red")
axvline(ηt / ϵ, linestyle = "dashed", color = "green")
xlabel(L"\lambda")
title("Eigenvalue distribution")
legend()
