Commit 7a5b2e07 authored by lindnemi's avatar lindnemi

some adjustments

parent 68b465a3
# For explanations check the notebook
cd(@__DIR__)
......@@ -9,17 +8,22 @@ using DelimitedFiles
using LightGraphs
using NetworkDynamics
using OrdinaryDiffEq
using Plots
# In this script we are not plotting anything
#using Plots
#G = readdlm(joinpath(@__DIR__, "PIAB_16T5Z.txt"), '\t', Float64, '\n')
G = readdlm("/Users/jess/Desktop/Different_Variations_DTI/Root6_DTI.txt", '\t', Float64, '\n')
G = readdlm(joinpath(@__DIR__, "Root3_DTI.txt"), '\t', Float64, '\n')
g = SimpleDiGraph(G)
N = nv(g)
#synaptic_weights = G + G'
synaptic_weights = G
#synaptic_weights ./= max(synaptic_weights...)
# Extract edge weights
edge_weights = Array{Float64}(undef, 1, ne(g))
for (i, edge) in zip(1:ne(g), collect(edges(g)))
edge_weights[i] = synaptic_weights[edge.src, edge.dst]
end
# FHN parameters
const ϵ = 0.1 # time separation
......@@ -32,19 +36,8 @@ const α = 1. / 3.
const = a * ϵ
const = b * ϵ
#we only need random x0 for the first step
x0 = rand(2N) * 4 .- 2
for sig in [10,5,2,1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0]
#for sig in [0,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,2,5,10]
global x0
σ = sig
# Weight parameters
edgep = Array{Float64}(undef, 1, ne(g))
for (i, edge) in zip(1:ne(g), collect(edges(g)))
edgep[i] = synaptic_weights[edge.src, edge.dst]
end
edgep .*= σ # multiply by coupling strength beforehand to save computations later
# Defining network functions
begin
@inline Base.@propagate_inbounds function fhn_stp!(, ξ, edges, p, t)
# u, v = ξ[1:2]
# FHN model
......@@ -62,110 +55,28 @@ for sig in [10,5,2,1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0]
nothing
end
electricaledge = StaticEdge(f! = electrical_edge!, dim = 1)
electricaledge = StaticEdge(f! = electrical_edge!, dim = 1, coupling = :directed)
odevertex = ODEVertex(f! = fhn_stp!, dim = 2, sym = [:u, :v])
fhn_network! = network_dynamics(odevertex, electricaledge, g);
end
#we only need random x0 for the first step
x0 = rand(2N) * 4 .- 2
for σ in [5.,2.,1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.]
#for σ in [0,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,2,5,10]
global x0
# Edge parameters are edge_weights multiplied with coupling strength
edgep = σ * edge_weights
p = (nothing, edgep) # network_dynamics uses tuple syntax for parameters
tspan = (0., 10000.)
tspan = (0., 10.)
prob = ODEProblem(fhn_network!, x0, tspan, p)
@time sol = solve(prob, Tsit5(), saveat=tspan[1]:1:tspan[end]);
@time sol = solve(prob, Tsit5(), reltol = 1e-6, saveat=tspan[1]:1:tspan[end])
writedlm( "/Volumes/004915257822556/sigma-outputs-model2/Root6/res1_rev_10t_u_v_$σ.csv", sol[:,end-1000:end], ',')
# writedlm( "/Volumes/004915257822556/sigma-outputs-model2/Root6/res1_rev_10t_u_v_$σ.csv", sol[:,end-1000:end], ',')
x0 = sol[:,end]
# plot(sol, tspan = tspan, vars = idx_containing(fhn_network!, :u), legend=false)
# uv_pl = plot!(sol, tspan = tspan, vars = idx_containing(fhn_network!, :v), legend=false, color = "lightblue")
# savefig(uv_pl,"/Volumes/004915257822556/sigma-outputs-model2/test_uv_pl_$σ.png")
# tstart = 1
# plot(sol[tstart:end], vars=(1,2))
# plot!(sol[tstart:end], vars=(5,6))
# plot!(sol[tstart:end], vars=(33,34))
# heatmap(Array(sol)[1:2:end,1:end])
# anim = @animate for t in 250:4:2000
# display_time = round(.2t, digits=3)
# scatter(sol[1:2:end, t], ylims = [-2, 2], title="t=$display_time", legend=false)
# end
# #gif(anim,"output/gif_osci.gif", fps=12)
# #############################
# using FFTW
# using Statistics
# spectrum= (abs.(fft(sol[1,250:1000] .- mean(sol[1,250:1000]))).^2)
# findall(x -> x == maximum(spectrum), spectrum)
# #vline([-0.813, 0.813]./2, label = "T = 2.46", linewidth = 2.5, color=:red)
# plot( collect(-75:0.2:75) ./ 150, abs.(fft(sol[1,250:1000] .- mean(sol[1,250:1000]))).^2, color = :darkblue, xlabel="Hz")
# #xticks!([-1, -0.813, -0.5,0.,0.5,0.813,1.] ./2)
# const affected_vars = collect(1:2:40)
# const input_strength = 100.
# """@docs
# Wrapper for the fhn_network!. Has the same calling signature and will add the inputs when
# p[1][1] = true. p is assumed to be a tuple ([trigger], (vertexp, edgep)). The trigger is
# wrapped in an array since tuples are immutable, but may contain Arrays with mutable
# elements.
# """
# function fhn_input_network!(dx, x, p, t)
# fhn_network!(dx, x, p[2], t)
# if p[1][1]
# dx[affected_vars] .+= input_strength
# end
# end;
# input_network! = ODEFunction(fhn_input_network!,
# syms = fhn_network!.syms,
# mass_matrix = fhn_network!.mass_matrix);
# function affect!(integrator)
# # a ? b : c is a fancy julia syntax for if a do b else c
# integrator.p[1][1] ? integrator.p[1][1] = false : integrator.p[1][1] = true
# end;
# using DiffEqCallbacks
# Δt = 200. # after Δt the callback is applied
# cb = PeriodicCallback(affect!, Δt, initial_affect=false)
# input_p = ([false], p)
# input_prob = ODEProblem(input_network!, x0, tspan, input_p)
# sol_input = solve(input_prob, AutoTsit5(TRBDF2()), saveat=tspan[1]:.5:tspan[end],
# callback= cb);
# plot(sol_input, tspan = tspan, vars = idx_containing(fhn_network!, :u), legend=false, color = "darkred")
# uv_pl_input = plot!(sol_input, tspan = tspan, vars = idx_containing(fhn_network!, :v), legend=false, color = "lightblue")
# savefig(uv_pl_input,"test_uv_pl_input_$σ.png")
#writedlm( "/Volumes/004915257822556/sigma-output-compactInfo/reltol_rev_100t_u_v_x_y_$σ.csv", sol[:,end-1000:end], ',')
# display(@benchmark solve(input_prob, AutoTsit5(TRBDF2()), saveat=tspan[1]:.5:tspan[end], callback= cb))
# using Statistics: cov, mean
# window_length = 60
# step_width = 20
# uvars = idx_containing(fhn_network!, :u)
# function get_covmat_list(sol, vars; window_length = 50, step_width = 10)
# covmat_list = Array{Float64,2}[]
# start_list = Float64[]
# for start = 1:step_width:(length(sol) - window_length)
# stop = start + window_length
# push!(covmat_list, cov(sol[vars, start:stop], sol[vars, start:stop]; dims=2))
# push!(start_list, start)
# end
# return zip(covmat_list, start_list)
# end
# covmat_list, start_list = get_covmat_list(sol, uvars)
# anim2 = @animate for (covmat, start) ∈ get_covmat_list(sol_input, uvars, window_length= window_length, step_width =step_width)
# covmat = covmat
# heatmap(log.(covmat ./ maximum(abs.(covmat)) .+1), clims=(-1,1), title="t=$((start-1)/2)", c= :vikO)
# end
# gif(anim2, "cov_window_60_fps3.gif", fps = 3)
# # include("speed_experiments_mtk.jl")
end
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment