change from internal integrator to DiffEq.jl (only IZ)

This is only a first try and so far introduces breaking changes!

What works:
Integration of several IZ neurons over time (see examples/iz_neuron.jl)

Project.toml

@@ -10,11 +10,11 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-julia = "1"
 Parameters = "0.12"
 Reexport = "0.2"
 Requires = "0.5"
 Unitful = "0.17"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

examples/if_neuron.jl

@@ -1,8 +1,10 @@
-using Plots, SNN
+using Plots, SpikingNeuralNetworks
+
+const SNN = SpikingNeuralNetworks
 
 E = SNN.IF(;N = 1)
 E.I = [11]
 SNN.monitor(E, [:v, :fire])
 
-SNN.sim!([E], []; duration = 300ms)
+SNN.sim!([E], []; duration = 400SNN.ms)
 SNN.vecplot(E, :v) |> display

examples/iz_neuron.jl

@@ -1,24 +1,35 @@
-using Plots, SNN
+using Plots, SpikingNeuralNetworks
+const SNN = SpikingNeuralNetworks
 
-RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
-IB = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -55, d = 4))
-CH = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -50, d = 2))
-FS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
-TC1 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
-TC2 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
-RZ = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.26, c = -65, d = 2))
-LTS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.25, c = -65, d = 2))
-P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]
+v0 = fill(-65.0, 6)
+u0 = fill(0.2*-65, 6)
+vu0 = SNN.ArrayPartition(v0, u0)
+p = SNN.IZParameter()
+tspan = [0.0,100.0]
+cb = SNN.DiscreteCallback(SNN.fire,SNN.affect!)
+prob = SNN.ODEProblem(SNN.integrate!, vu0, tspan, p)
 
-SNN.monitor(P, [:v])
-T = 2second
-for t = 0:T
-    for p in [RS, IB, CH, FS, LTS]
-        p.I = [10]
-    end
-    TC1.I = [(t < 0.2T) ? 0mV : 2mV]
-    TC2.I = [(t < 0.2T) ? -30mV : 0mV]
-    RZ.I =  [(0.5T < t < 0.6T) ? 10mV : 0mV]
-    SNN.sim!(P, [], 0.1ms)
-end
-SNN.vecplot(P, :v) |> display
+sol = SNN.solve(prob,SNN.Tsit5();callback = cb)
+
+#RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
+#IB = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -55, d = 4))
+#CH = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -50, d = 2))
+#FS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
+#TC1 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
+#TC2 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
+#RZ = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.26, c = -65, d = 2))
+#LTS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.25, c = -65, d = 2))
+#P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]
+
+#SNN.monitor(P, [:v])
+#T = 2second
+#for t = 0:T
+#    for p in [RS, IB, CH, FS, LTS]
+#        p.I = [10]
+#    end
+#    TC1.I = [(t < 0.2T) ? 0mV : 2mV]
+#    TC2.I = [(t < 0.2T) ? -30mV : 0mV]
+#    RZ.I =  [(0.5T < t < 0.6T) ? 10mV : 0mV]
+#    SNN.sim!(P, [], 0.1ms)
+#end
+#SNN.vecplot(P, :v) |> display

src/SpikingNeuralNetworks.jl

@@ -1,11 +1,17 @@
 module SpikingNeuralNetworks
 
+export SNN
+
 const SNN = SpikingNeuralNetworks
 
 using SparseArrays
 using Reexport
 using Parameters
 using Requires
+using OrdinaryDiffEq
+using DifferentialEquations
+using Parameters
+using RecursiveArrayTools
 #using Unitful
 #using Unitful.DefaultSymbols
 #@reexport using Utils

src/neuron/iz.jl

@@ -1,31 +1,24 @@
 @with_kw struct IZParameter
-    a::SNNFloat = 0.01
-    b::SNNFloat = 0.2
-    c::SNNFloat = -65
-    d::SNNFloat = 2
+    a::Float64 = 0.01
+    b::Float64 = 0.2
+    c::Float64 = -65
+    d::Float64 = 2
+    I::Float64 = 11
 end
 
-@with_kw mutable struct IZ
-    param::IZParameter = IZParameter()
-    N::SNNInt = 100
-    v::Vector{SNNFloat} = fill(-65.0, N)
-    u::Vector{SNNFloat} = param.b * v
-    fire::Vector{Bool} = zeros(Bool, N)
-    I::Vector{SNNFloat} = zeros(N)
-    records::Dict = Dict()
+function integrate!(dv, v, p, t)
+    @unpack a, b, c, d, I = p
+    dv.x[1] .= 0.04f0 .* v.x[1] .^2 .+ 5f0 .* v.x[1] .+ 140f0 .- v.x[2] .+ I
+    dv.x[2] .= a .* (b .* v.x[1] - v.x[2])
 end
 
-function integrate!(p::IZ, param::IZParameter, dt::SNNFloat)
-    @unpack N, v, u, fire, I = p
-    @unpack a, b, c, d = param
-    @inbounds for i = 1:N
-        v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
-        v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
-        u[i] += dt * (a * (b * v[i] - u[i]))
-    end
-    @inbounds for i = 1:N
-        fire[i] = v[i] > 30f0
-        v[i] = ifelse(fire[i], c, v[i])
-        u[i] += ifelse(fire[i], d, 0f0)
-    end
+function fire(v,t,integrator)
+    any(30f0 .- v.x[1] .> 0)
 end
+
+function affect!(integrator)
+    F = integrator.u.x[1] .- 30f0 .> 0
+    integrator.u.x[1][F] .= integrator.p.c
+    integrator.u.x[2][F] .+= integrator.p.d
+end
+