WIP: cleaning and revisiting spDCM tutorial

docs/Project.toml

@@ -8,13 +8,15 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Neuroblox = "769b91e5-4c60-41ee-bfae-153c84203cb2"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [compat]

docs/src/tutorials/spectralDCM.jl

@@ -31,6 +31,7 @@ using OrderedCollections
 using CairoMakie
 using ModelingToolkit
 using Random
+using StatsBase
 
 # # Model simulation
 # ## Define the model
@@ -54,28 +55,29 @@ for i = 1:nr
     push!(regions, region)          # store neural mass model in list. We need this list below. If you haven't seen the Julia command `push!` before [see here](http://jlhub.com/julia/manual/en/function/push-exclamation).
 
     ## add Ornstein-Uhlenbeck block as noisy input to the current region
-    input = OUBlox(;name=Symbol("r$(i)₊ou"), σ=0.1)
+    input = OUBlox(;name=Symbol("r$(i)₊ou"), σ=0.2, τ=2)
     add_edge!(g, input => region, weight=1/16)   # Note that 1/16 is taken from SPM12, this stabilizes the balloon model simulation. Alternatively the noise of the Ornstein-Uhlenbeck block or the weight of the edge connecting neuronal activity and balloon model could be reduced to guarantee numerical stability.
 
     ## simulate fMRI signal with BalloonModel which includes the BOLD signal on top of the balloon model dynamics
     measurement = BalloonModel(;name=Symbol("r$(i)₊bm"))
     add_edge!(g, region => measurement, weight=1.0)
 end
 # Next we define the between-region connectivity matrix and connect regions; we use the same matrix as is defined in [3]
-A_true = [[-0.5 -2 0]; [0.4 -0.5 -0.3]; [0 0.2 -0.5]]
+A_true = [[-0.5 -0.2 0]; [0.4 -0.5 -0.3]; [0 0.2 -0.5]]
+# Note that in SPM DCM connection matrices column variables denote output from and rows denote inputs to a particular region.
+# This is different from the usual Neuroblox definition of connection matrices. Thus we flip the indices in what follows:
 for idx in CartesianIndices(A_true)
-    add_edge!(g, regions[idx[1]] => regions[idx[2]], weight=A_true[idx[1], idx[2]])
+    add_edge!(g, regions[idx[1]] => regions[idx[2]], weight=A_true[idx[2], idx[1]])   # Note the definition of columns as outputs and rows as inputs. For consistency with SPM we keep this notation.
 end
 
 # finally we compose the simulation model
 @named simmodel = system_from_graph(g)
 
 # ## Run the simulation and plot the results
-
+tspan = (0, 1022)
+dt = 2   # 2 seconds as measurement interval for fMRI
 # setup simulation of the model, time in seconds
-tspan = (0.0, 512.0)
 prob = SDEProblem(simmodel, [], tspan)
-dt = 2   # 2 seconds (units are seconds) as measurement interval for fMRI
 sol = solve(prob, ImplicitRKMil(), saveat=dt);
 
 # we now want to extract all the variables in our model which carry the tag "measurement". For this purpose we can use the Neuroblox function `get_idx_tagged_vars`
@@ -89,30 +91,51 @@ ax = Axis(f[1, 1],
     xlabel = "Time [ms]",
     ylabel = "BOLD",
 )
-lines!(ax, sol, idxs=idx_m)
+lines!(ax, sol, idxs=idx_m);
 f
 
 # We note that the initial spike is not meaningful and a result of the equilibration of the stochastic process thus we remove it.
-dfsol = DataFrame(sol[ceil(Int, 101/dt):end]);
+dfsol = DataFrame(sol);
+
+# ## Add measurement noise and rescale data
+data = Matrix(dfsol[:, idx_m .+ 1]);    # +1 due to the additional time-dimension in the data frame.
+# add measurement noise
+data += randn(size(data))/4;
+# center and rescale data (as done in SPM):
+data .-= mean(data, dims=1);
+data *= 1/std(data[:])/4;
+dfsol = DataFrame(data, :auto);
+# Add correct names to columns of the data frame
+_, obsvars = get_eqidx_tagged_vars(simmodel, "measurement");  # get index of equation of bold state
+rename!(dfsol, Symbol.(obsvars))
 
 # ## Estimate and plot the cross-spectral densities
-data = Matrix(dfsol[:, idx_m]);
+
 # We compute the cross-spectral density by fitting a linear model of order `p` and then compute the csd analytically from the parameters of the multivariate autoregressive model
-p = 8
-mar = mar_ml(data, p)   # maximum likelihood estimation of the MAR coefficients and noise covariance matrix
-ns = size(data, 1)
-freq = range(min(128, ns*dt)^-1, max(8, 2*dt)^-1, 32)
+p = 8;
+mar = mar_ml(data, p);   # maximum likelihood estimation of the MAR coefficients and noise covariance matrix
+ns = size(data, 1);
+freq = range(min(128, ns*dt)^-1, max(8, 2*dt)^-1, 32);
 csd = mar2csd(mar, freq, dt^-1);
-# Now plot the cross-spectrum:
+# Now plot the real part of the cross-spectra. Most part of the signal is in the lower frequencies:
 fig = Figure(size=(1200, 800))
 grid = fig[1, 1] = GridLayout()
 for i = 1:nr
     for j = 1:nr
-        ax = Axis(grid[i, j])
+        if i == 1 && j == 1
+            ax = Axis(grid[i, j], xlabel="Frequency [Hz]", ylabel="real value of CSD")
+        else
+            ax = Axis(grid[i, j])
+        end
         lines!(ax, freq, real.(csd[:, i, j]))
     end
 end
+Label(grid[1, 1:3, Top()], "Cross-spectral densities", valign = :bottom,
+    font = :bold,
+    fontsize = 32,
+    padding = (0, 0, 5, 0))
 fig
+# These cross-spectral densities are the data we use in spectral DCM to fit our model to and perform the inference of connection strengths.
 
 # # Model Inference
 
@@ -143,6 +166,10 @@ end
 # Here we define the prior expectation values of the effective connectivity matrix we wish to infer:
 A_prior = 0.01*randn(nr, nr)
 A_prior -= diagm(diag(A_prior))    # remove the diagonal
+# These two parameters are not present in the ground truth thus set them to zero and set their tuning parameter to false: 
+A_prior[3] = 0.0
+A_prior[7] = 0.0
+
 # Since we want to optimize these weights we turn them into symbolic parameters:
 # Add the symbolic weights to the edges and connect regions.
 A = []
@@ -160,20 +187,18 @@ for (i, idx) in enumerate(CartesianIndices(A_prior))
 end
 # Avoid simplification of the model in order to be able to exclude some parameters from fitting
 @named fitmodel = system_from_graph(g, simplify=false)
-# With the function `changetune`` we can provide a dictionary of parameters whose tunable flag should be changed, for instance set to false to exclude them from the optimizatoin procedure.
-# For instance the the effective connections that are set to zero in the simulation:
-untune = Dict(A[3] => false, A[7] => false)
-fitmodel = changetune(fitmodel, untune)                 # 3 and 7 are not present in the simulation model
-fitmodel = structural_simplify(fitmodel)   # and now simplify the euqations; the `split` parameter is necessary for some ModelingToolkit peculiarities and will soon be removed. So don't lose time with it ;)
+# With the function `changetune`` we can provide a dictionary of parameters whose tunable flag should be changed, for instance set to false to exclude them from the optimization procedure.
+# For instance the effective connections that are set to zero in the simulation and the self-connections:
+untune = Dict(A[3] => false, A[7] => false, A[1] => false, A[5] => false, A[9] => false)
+fitmodel = changetune(fitmodel, untune)           # 3 and 7 are not present in the simulation model
+fitmodel = structural_simplify(fitmodel)          # and now simplify the euqations
 
 # ## Setup spectral DCM
 max_iter = 128;            # maximum number of iterations
 ## attribute initial conditions to states
 sts, _ = get_dynamic_states(fitmodel);
-sts = unknowns(fitmodel)
-rv = rand(length(sts))
 # the following step is needed if the model's Jacobian would give degenerate eigenvalues when expanded around the fixed point 0 (which is the default expansion). We simply add small random values to avoid this degeneracy:
-perturbedfp = OrderedDict(sts .=> abs.(0.001*rv[1:length(sts)]))     # slight noise to avoid issues with Automatic Differentiation.
+perturbedfp = Dict(sts .=> abs.(10^-10*rand(length(sts))))     # slight noise to avoid issues with Automatic Differentiation.
 # For convenience we can use the default prior function to use standardized prior values as given in SPM:
 pmean, pcovariance, indices = defaultprior(fitmodel, nr)
 
@@ -186,11 +211,8 @@ hyperpriors = (Πλ_pr = 128.0*ones(1, 1),   # prior metaparameter precision, ne
               );
 # To compute the cross spectral densities we need to provide the sampling interval of the time series, the frequency axis and the order of the multivariate autoregressive model:
 csdsetup = (mar_order = p, freq = freq, dt = dt);
-# earlier we used the function `get_idx_tagged_vars` to get the indices of tagged variables. Here we don't want to get the indices but rather the symbolic variable names themselves.
-# and in particular we need to get the measurement variables in the same ordering as the model equations are defined.
-_, s_bold = get_eqidx_tagged_vars(fitmodel, "measurement");    # get bold signal variables
 # Prepare the DCM. This function will setup the computation of the Dynamic Causal Model. The last parameter specifies that wer are using fMRI time series as opposed to LFPs.
-(state, setup) = setup_sDCM(dfsol[:, String.(Symbol.(s_bold))], fitmodel, perturbedfp, csdsetup, priors, hyperpriors, indices, pmean, "fMRI");
+(state, setup) = setup_sDCM(dfsol, fitmodel, perturbedfp, csdsetup, priors, hyperpriors, indices, pmean, "fMRI");
 
 # We are now ready to run the optimization procedure! :)
 # That is we loop over run_sDCM_iteration! which will alter `state` after each optimization iteration. It essentially computes the Variational Laplace estimation of expectation and variance of the tunable parameters. 
@@ -206,6 +228,7 @@ for iter in 1:max_iter
         end
     end
 end
+# Note that the output `F` is the free energy at each iteration step and `dF` is the predicted change of free energy at each step which approximates the actual free energy change and is used as stopping criterion by requiring that it does not excede the `tolerance` level for 4 consecutive times.
 
 # # Plot Results
 # Free energy is the objective function of the optimization scheme of spectral DCM. Note that in the machine learning literature this it is called Evidence Lower Bound (ELBO). 

src/Neuroblox.jl

@@ -236,7 +236,7 @@ export LearningBlox
 export CosineSource, CosineBlox, NoisyCosineBlox, PhaseBlox, ImageStimulus, ConstantInput, ExternalInput, SpikeSource, PoissonSpikeTrain, generate_spike_times
 export DBS, ProtocolDBS, detect_transitions, compute_transition_times, compute_transition_values, get_protocol_duration
 export BandPassFilterBlox
-export OUBlox, OUCouplingBlox
+export OUBlox, OUCouplingBlox, ARBlox
 export phase_inter, phase_sin_blox, phase_cos_blox
 export SynapticConnections, create_rl_loop
 export add_blox!, get_system

src/Neurographs.jl

@@ -116,7 +116,7 @@ function LIF_spike_affect!(integ, u, p, ctx)
     integ.p[p[4]] = 1
 
     SciMLBase.add_tstop!(integ, t_refract_end)
-    
+
     c = 1
     for i in eachindex(u)[2:end]
         integ.u[u[i]] += integ.p[p[c + 4]]
@@ -221,9 +221,9 @@ function generate_discrete_callbacks(blox::Union{LIFExciNeuron, LIFInhNeuron}, b
 
     states_affect = get_states_spikes_affect(sa, name_blox)
     params_affect = get_params_spikes_affect(sa, name_blox)
-   
+
     param_pairs = make_unique_param_pairs(params_affect)
-    
+
     ps = vcat([
         sys.V_reset => Symbol(sys.V_reset), 
         sys.t_refract_duration => Symbol(sys.t_refract_duration), 

src/blox/stochastic.jl

@@ -22,12 +22,35 @@ mutable struct OUBlox <: NeuralMassBlox
         sts = @variables x(t)=0.0 [output=true] jcn(t) [input=true]
         @brownian w
 
-        eqs = [D(x) ~ -(x-μ)/τ + jcn + sqrt(2/τ)*σ*w]
+        eqs = [D(x) ~ (-x + μ + jcn)/τ + sqrt(2/τ)*σ*w]
         sys = System(eqs, t; name=name)
         new(namespace, true, sys)
     end
 end
 
+function get_ts_data(t, dt::Real, data::Array{Float64})
+    idx = ceil(Int, t / dt)
+
+    return idx > 0 ? data[idx] : 0.0
+end
+
+@register_symbolic get_ts_data(t, dt::Real, data::Array{Float64})
+
+mutable struct ARBlox <: StimulusBlox
+    namespace
+    system
+    function ARBlox(;name, namespace=nothing, dt, data)
+        sts = @variables u(t) = 0.0 [output=true]
+
+        eq = [u ~ get_ts_data(t, dt, data)]
+
+        sys = System(eq, t; name)
+
+        new(namespace, sys)
+    end
+end
+
+
 # """
 # Ornstein-Uhlenbeck Coupling Blox
 # This blox takes an input and multiplies that input with

src/datafitting/spDCM_VL.jl

@@ -311,15 +311,16 @@ function integration_step(dfdx, f, v, solenoid=false)
         dfdx = dfdx - Q*dfdx;        
     end
 
-    # (expm(dfdx*t) - I)*inv(dfdx)*f ~~~ could also be done with expv (expv(t, dFdθθ, dFdθθ \ dFdθ) - dFdθθ \ dFdθ) but doesn't work with Dual.
-    # could also be done with exponential! but isn't numerically stable
+    # NB: (exp(dfdx*t) - I)*inv(dfdx)*f, could also be done with expv (expv(t, dFdθθ, dFdθθ \ dFdθ) - dFdθθ \ dFdθ) but doesn't work with Dual.
+    # Could also be done with `exponential!` but isn't numerically stable.
+    # Thus, just use `exp`.
     n = length(f)
     t = exp(v - spm_logdet(dfdx)/n)
 
     if t > exp(16)
         dx = - dfdx \ f   # -inv(dfdx)*f
     else
-        dx = (exp(t * dfdx) - I) * inv(dfdx)*f # (expm(dfdx*t) - I)*inv(dfdx)*f
+        dx = (exp(t * dfdx) - I) * inv(dfdx) * f # (expm(dfdx*t) - I)*inv(dfdx)*f
     end
 
     return dx

src/measurementmodels/fmri.jl

@@ -74,7 +74,7 @@ struct BalloonModel <: ObserverBlox
         p = paramscoping(lnκ=lnκ, lnτ=lnτ, lnϵ=lnϵ)     # finally compile all parameters
         lnκ, lnτ, lnϵ = p                               # assign the modified parameters
 
-        sts = @variables s(t)=1.0 lnu(t)=0.0 lnν(t)=0.0 lnq(t)=0.0 bold(t)=0.0 [irreducible=true, output=true, description="measurement"] jcn(t) [input=true]
+        sts = @variables s(t)=0.0 lnu(t)=0.0 lnν(t)=0.0 lnq(t)=0.0 bold(t)=0.0 [irreducible=true, output=true, description="measurement"] jcn(t) [input=true]
 
         eqs = [
             D(s)   ~ jcn - H[1]*exp(lnκ)*s - H[2]*(exp(lnu) - 1),
@@ -129,7 +129,7 @@ function boldsignal(;name, lnϵ=0.0)
     k1  = 4.3*nu0*E0*TE
 
     params = @parameters lnϵ=lnϵ
-    vars = @variables bold(t) lnν(t) lnq(t)   # TODO: got to be really careful with the current implementation! A simple permutation of this breaks the algorithm!
+    vars = @variables bold(t) lnν(t) lnq(t)
 
     eqs = [
         bold ~ V0*(k1 - k1*exp(lnq) + exp(lnϵ)*r0*E0*TE - exp(lnϵ)*r0*E0*TE*exp(lnq)/exp(lnν) + 1-exp(lnϵ) - (1-exp(lnϵ))*exp(lnν))