WIP: cleaning and revisiting spDCM tutorial

[david-hofmann]

Setting everything up so as to reproduce the nice results from the Friston papers in the spDCM tutorial

[david-hofmann]

@harisorgn I introduced this equation because otherwise it would complain that the system is unbalanced. However, I am a bit concerned that this way we will produce impulse like perturbations of the system at each t_stim point after which x will again be set to 0? Not sure but once we can simulate it we will know, perhaps you know already the expected behavior.

[harisorgn]

I think the solution is to remove the initial guess =0.0 from your x(t) state, but in general you could also define x as a parameter, e.g. this is what we have done in ImageStimulus as it is only affected by callbacks, it doesn't have dynamics.

[harisorgn]

I think the solution is to remove the initial guess =0.0 from your x(t) state

Sorry, that's not true here. It would work if x was just some placeholder state like jcn that is simplified away. But here you want to keep it irreducible=true so making it a parameter is my recommended solution. Otherwise you will end up with a DAE system.

[david-hofmann]

@harisorgn I've almost got it but I still can't simulate the system. If you have time to give it a look and see if you understand why maybe you are quicker than I am.

This is the error I get when calling prob = SDEProblem(simmodel, [], tspan)
ERROR: MethodError: no method matching SDEProblem(::ODESystem, ::Vector{Any}, ::Tuple{Float64, Float64})

[harisorgn]

This is the error I get when calling prob = SDEProblem(simmodel, [], tspan)
ERROR: MethodError: no method matching SDEProblem(::ODESystem, ::Vector{Any}, ::Tuple{Float64, Float64})

Looks like you have defined an ODESystem but you are trying to solve an SDEProblem. I guess you need system and problem to match.

[david-hofmann]

oh, whoops! thank you! Now the affect throws an error:
ERROR: MethodError: no method matching compile_affect(::Equation, ::ModelingToolkit.SymbolicDiscreteCallback, ::ODESystem, ::Vector{…}, ::Vector{…}; expression::DataType, postprocess_affect_expr!::Nothing, eval_expression::Bool, eval_module::Module)

Okay, need to figure how to do this right it seems, maybe need a different approach as compared to the list of pairs

[david-hofmann]

ignore the wrong keyword solver here ;)

[david-hofmann]

also if you want to run this you also need to run a few lines below these:

vars = matread(joinpath(@__DIR__, "../../../test/spm25_demo.mat"));
tspan = (0.0, 1124.0)
t_burnin = 102   # ignore the first 102 seconds
dt = 2   # 2 seconds as measurement interval for fMRI
times = t_burnin:dt:tspan[2]

[david-hofmann]

@harisorgn @MasonProtter the problem with this is that it will try to call the Dict for every time, so most keywords won't exist. What am I missing? I think I should do it the same way as Haris did, that is I define some interval like (t >= first(t_stims[idx])) && (t <= last(t_stims[idx])). Shouldn't I?

[MasonProtter]

What I was picturing was something like

stim_dict = Dict(t1 => (i, v1), t2 => (j, v2), t3 => (k, v3), ...)
stim_condition(t) = keys(stim_dict)
function stim_affect!(integ, u, p, ctx)
    i, v = stim_dict[integ.t]
    integ.u[u[i]] = v
end

inside of the ARBlox struct, but I'm not totally sure.

[harisorgn]

Yes David, to treat u as a state and write a symbolic equation for it you would need the ifelse I think. Because even in Mason's suggestion stim_dict[integ.t] this will be problematic, right? We don't know a priori all integ.t values to prefill a Dict.

[harisorgn]

or you keep the Dict and u state and just check for haskey(stim_dict, integ.t), which also has a control flow I guess similar to the ifelse solution.

[david-hofmann]

I've solved it with the ifelse solution or rather the ternary operator which seemed to me a bit more suited since I can avoid the + 1 as in this case:

Neuroblox.jl/src/blox/sources.jl

Line 104 in 36714a3

idx = floor(Int, t / t_trial) + 1

One question: wouldn't @MasonProtter 's solution also create 512 callbacks? In my understanding I'd have as many callbacks and each time I'd call stim_affect! or am I still not understanding how discrete events work? Which is quite likely. But [x ~ 0] => (affect!, [v, x], [p, q], [discretes...], ctx) from the docs is what we do and I'd have 512 such time points on the left hand side instead of one as is here I think.

@harisorgn , I do have all the times given and I do saveat with a time interval to match those times, so I suppose we do?

However, the general problem persists with the solution I have. The noise was not the issue after all, it would have been strange. Anyhow, digging deeper.

[harisorgn]

wouldn't @MasonProtter 's solution also create 512 callbacks?

No, you would have a single callback triggered at keys(stim_dict), each time it triggers it calls stim_affect! and inside this function you decide which value to set u to.

[harisorgn]

Suggested change

	return idx > 0 ? data[idx] : 0.0
	return ifelse(idx > 0, data[idx], 0.0)

use this instead of the ternary operator in symbolic expressions

[david-hofmann]

Oh my, I did not want to commit but to comment.... anyhow.

[david-hofmann]

the thing with ifelse is that it evaluates every expression whereas the ternary operator only evaluates the relevant one based on the condition. Thus I spare myself the little index "hack" and it works with the symbolic expression/function it seems. Or at least I didn't get any error?

[harisorgn]

This is my source
https://docs.sciml.ai/ModelingToolkit/stable/basics/FAQ/#How-do-I-handle-if-statements-in-my-symbolic-forms?

ifelse is that it evaluates every expression

that is a good thing for AD IIRC. In any case I think you can restructure the function to use ifelse with two branches that won't break on evaluation.