kanav99 on nlsolve
(wip) change to NonlinearSolve.… (compare)
ChrisRackauckas on master
Fix typo in ReadMe (compare)
[slack] <Ramiro Vignolo> Hi guys. I have a struct that contains a ODEProblem inside (among other things) and has a call method that remakes the problem by changing its tspan (this is an ODE problem with TERMINAL, not INITIAL, condition) and solves it.
```struct ZCB{P}
prob::P
end
function (zcb::ZCB)(T::Real)
prob = zcb.prob
new_prob = remake(prob, tspan = (T, zero(T)))
sol = sol(new_prob)
endSo far, so good. However, I would like to store the solution for a given `T` in a cache, because if the user performs a new call using the same `T` as before, I would like to avoid solving the ODEProblem. In this context, should I define the struct to be something like:
mutable struct ZCB{P,C}
prob::P
cache::C
end``
How can I construct an instance of ZCB? Should I perform a
solvein its constructor? On the other hand, should I use a dictionary as a
cachesuch as
Dict{Real,SolutionObject}so I can search the solution for a specific
T` as its keyword?
occursin(xs::Vector, expr)
which tells you if each element of xs
appears in expr would help you speed up stuff? I see that this could help for example in ismassaction
. What are some other primitives like this that maybe helpful?
Well I'll just leave this here: I'm trying to implement in DiffEq the update rule for a spiking neural network. For now I've implemented with plain Julia and works decently, but I'd like to speed it up a little:
The function simulate
just initialize the variables and then runs step! for 10000 times.
The function step!
takes as input:
V: a matrix(neurons x time), where values are potentials of individual neurons at a certain time
sd: a vector(neurons), where values represent how long ago a given neuron spiked
w: a matrix(neurons x neurons) of weights (connectivity matrix)
np, ep are parameters
tottime is the total time of simulation
function simulate(tottime, np, sp, ep)
w = genweights(np, sp)
V = zeros(np.N, tottime)
sd = zeros(np.N) .+ 0.2
V[:, 1:4] .= np.Vᵣ
V[126:132,1:4] .= 0.
for t in 3:tottime-1
step!(V, sd, t, w, np, ep, tottime)
end
V[:, 1:tottime], sd
end
function step!(V, sd, t, w, np, ep, tottime)
V[:, t+1] = dv(view(V, :, t), sd, w, np, ep) # This is the update
V[view(V,:, t) .> np.Vₜ, t+1] .= 0.0 # This makes the neurons over the threshold spike
V[view(V,:, t) .== 0.0, t+1] .= np.Vᵣ # This reset the neurons that spiked to a low potential
sd .+= ep.dt # Here we adjust the spike delays for the next timestep
sd[view(V,:, t-2) .== 0.0] .= 3*ep.dt # This sets the spike delays to 3 timesteps for the neurons who spiked 2 timesteps ago, emulating transmission delays or something
end
end
function dv(v, sd, w, np, ep)
k = sd .* exp.(-sd ./ np.τₛ) # This is just to save computations
v.+ ((Eₗₜ .- v ./ np.τₘ ).- (v.-np.Eₑ) .*( w[2] * k) .- (v.-np.Eᵢ) .*( w[1] * k)) .* ep.dt # This is the differential equation, evaluated with Euler
end
The function step!
runs 10000 times, simulating 10 seconds of activity, and it takes ~350ms
I'd like to implement this in DiffEq, but I don't understand how to write a function in the way DiffEq wants it and having it computing the spikes.
Specifically I'm confused by:
solve
do to my matrix: should I multiply by dt, or will solve
do it?I came up with this, I don't know if it makes any sense. I think there is a main problem in updating sd
# This function initializes all the needed arrays
function DEsimulate(tottime, np, sp, ep)
alg = FunctionMap{true}()
w = genweights(np, sp)
v = zeros(np.N) .+ np.Vᵣ
sd = zeros(np.N) .+ 0.2
v[1:5] .= 0.
p = (w, np, ep)
prob = DiscreteProblem(update!, [v, sd], (ep.dt, tottime*ep.dt), p)
sim = solve(prob, alg, dt=ep.dt)
return sim
end
function update!(dv, u, p, t)
w, np, ep = p
v, sd = u
v[v .== 0.0] .= np.Vᵣ
v[v .> np.Vₜ] .= 0.0
k = sd .* exp.(-sd ./ np.τₛ)
dv[1] += ((Eₗₜ .- v ./ np.τₘ ).- (v.-np.Eₑ) .*( w[2] * k) .- (v.-np.Eᵢ) .*( w[1] * k))
dv[2][v .== 0.0] .= 2*ep.dt
dv[2] .+= 1
end
Any kind of hint, also on style or everything else is very well appreciated!
solve
step of the very first example
findfirst
or similar as the index into an array, but it returned nothing
vGLCc
and so forth are algebraic equations written within the model
function. Anything stand out as to why this error could've occurred?
LabelledArrays.jl
, hence the du.
and p.
syntax.