From e169929b3a96f5af16da1d55884bac776614d912 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 25 Sep 2023 14:44:59 -0400 Subject: [PATCH 01/28] Four half cosine parameterization Added for convolution option --- src/measurementmodels/hemodynamic_extras.jl | 39 +++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 src/measurementmodels/hemodynamic_extras.jl diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl new file mode 100644 index 00000000..79e3da5a --- /dev/null +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -0,0 +1,39 @@ +using Random, Plots + + +""" +Woolrich, Behrens and Smith. 2003. + +Specify widths (h₁, h₂, h₃, h₄) in seconds. +Amplitudes (f₁, f₂) relative to HRF height. +Output is HRF kernel in ms. + +""" + +function HRFFourHalfCosine(; h₁=nothing, h₂=nothing, h₃=nothing, h₄=nothing, f₁=nothing, f₂=nothing) + + h₁ = isnothing(h₁) ? 2*rand() : h₁ + h₂ = isnothing(h₂) ? 4*rand()+2 : h₂ + h₃ = isnothing(h₃) ? 4*rand()+2 : h₃ + h₄ = isnothing(h₄) ? 6*rand()+2 : h₄ + f₁ = isnothing(f₁) ? 0 : f₁ + f₂ = isnothing(f₂) ? 0.5*rand() : f₂ + + t₁ = 0:0.001:h₁ + t₂ = 0:0.001:h₂ + t₃ = 0:0.001:h₃ + t₄ = 0:0.001:h₄ + + cos₁ = 0.5.*f₁.*(cos.(π.*(t₁./h₁)) .- 1) + cos₂ = ((0.5*(f₁+1)).*(-cos.(π.*(t₂./h₂)) .+ 1)) .- f₁ + cos₃ = ((0.5*(f₂+1)).*(cos.(π.*(t₃./h₃)))) .+ ((0.5*(1-f₂))) + cos₄ = (0.5.*f₂.*(-cos.(π.*(t₄./h₄)) .+ 1)) .- f₂ + + hrf = vcat(cos₁, cos₂, cos₃, cos₄) + +end + + + +hmm = HRFFourHalfCosine(f₁=0) +plot(hmm) \ No newline at end of file From 519f1da9d21b040351ae71607c86689717cd5463 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 25 Sep 2023 15:10:14 -0400 Subject: [PATCH 02/28] Double gamma added for convolution Same format as the four half cosine model --- Project.toml | 1 + src/measurementmodels/hemodynamic_extras.jl | 23 +++++++++++++++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 51168ddc..e3c3bc69 100644 --- a/Project.toml +++ b/Project.toml @@ -52,6 +52,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" SignalAnalysis = "df1fea92-c066-49dd-8b36-eace3378ea47" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index 79e3da5a..fc8bc709 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -1,4 +1,4 @@ -using Random, Plots +using Random, SpecialFunctions, Plots """ @@ -33,7 +33,22 @@ function HRFFourHalfCosine(; h₁=nothing, h₂=nothing, h₃=nothing, h₄=noth end - +function HRFDoubleGamma(; A=nothing, α₁=nothing, α₂=nothing, β₁=nothing, β₂=nothing, c=nothing, tlen=nothing) -hmm = HRFFourHalfCosine(f₁=0) -plot(hmm) \ No newline at end of file + if (!isnothing(A) || !isnothing(α₁) || !isnothing(α₂) || !isnothing(β₁) || !isnothing(β₂) || !isnothing(c)) && isnothing(tlen) + throw(DomainError(nothing, "If you specify a different parameter, you need to specify tlen as well to ensure the kernel is wide enough.")) + end + + A = isnothing(A) ? 1 : A + α₁ = isnothing(α₁) ? 6 : α₁ + α₂ = isnothing(α₂) ? 16 : α₂ + β₁ = isnothing(β₁) ? 1 : β₁ + β₂ = isnothing(β₂) ? 1 : β₂ + c = isnothing(c) ? 1.0/6.0 : c + tlen = isnothing(tlen) ? 30 : tlen + + t = 0:0.001:tlen + + hrf = A.*((((t.^(α₁-1)).*(β₁^α₁).*exp.(-β₁.*t))./gamma(α₁)) .- (c.*((t.^(α₂-1)).*(β₂^α₂).*exp.(-β₂.*t)./ gamma(α₂)))) + +end From a8a2cbe0898009b66bfa29cdc55c1db363917907 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 25 Sep 2023 17:25:46 -0400 Subject: [PATCH 03/28] Citation details --- src/measurementmodels/hemodynamic_extras.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index fc8bc709..81fa9ad1 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -33,6 +33,16 @@ function HRFFourHalfCosine(; h₁=nothing, h₂=nothing, h₃=nothing, h₄=noth end +""" +Lindquist et al. 2012. + +All parameters are in seconds. +Output is in ms. +If you specify a change in any of the parameters, also change the tlen to ensure the kernel is wide enough. +It will throw an error if you don't. + +""" + function HRFDoubleGamma(; A=nothing, α₁=nothing, α₂=nothing, β₁=nothing, β₂=nothing, c=nothing, tlen=nothing) if (!isnothing(A) || !isnothing(α₁) || !isnothing(α₂) || !isnothing(β₁) || !isnothing(β₂) || !isnothing(c)) && isnothing(tlen) From 0c1bb5f554504a90ae8c364f28ba409dcf48ca76 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Tue, 26 Sep 2023 11:14:50 -0400 Subject: [PATCH 04/28] Minor testing updates + exports Cleaning up to prepare for more fMRI simulations --- src/Neuroblox.jl | 2 ++ test/jansen_rit_component_tests.jl | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 204483c2..a22c21db 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -95,6 +95,7 @@ include("gui/GUI.jl") include("blox/blox_utilities.jl") include("blox/connections.jl") include("Neurographs.jl") +include("measurementmodels/hemodynamic_extras.jl") function simulate(sys::ODESystem, u0, timespan, p, solver = AutoVern7(Rodas4()); kwargs...) prob = ODEProblem(sys, u0, timespan, p) @@ -169,4 +170,5 @@ export simulate, random_initials export system_from_graph, graph_delays export create_adjacency_edges! export get_namespaced_sys, namespace_expr, nameof +export HRFFourHalfCosine, HRFDoubleGamma end \ No newline at end of file diff --git a/test/jansen_rit_component_tests.jl b/test/jansen_rit_component_tests.jl index 3554cfdb..5ffdbddd 100644 --- a/test/jansen_rit_component_tests.jl +++ b/test/jansen_rit_component_tests.jl @@ -56,7 +56,7 @@ blox = [Str, GPe, STN, GPi, Th, EI, PY, II] # test graphs g = MetaDiGraph() -add_blox!(Ref(g), blox) +add_blox!.(Ref(g), blox) # Store parameters to be passed later on params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5 @@ -112,7 +112,7 @@ create_adjacency_edges!(g2, adj_matrix_lin) @named final_system = system_from_graph(g2, params) final_delays = graph_delays(g2) -sim_dur = 10.0 # Simulate for 10 Seconds +sim_dur = 600.0 # Simulate for 10 Seconds final_system_sys = structural_simplify(final_system) prob = DDEProblem(final_system_sys, [], From f781918d524997bb178e6cfe452c3eacec76a8c0 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Wed, 27 Sep 2023 19:41:19 -0400 Subject: [PATCH 05/28] Updates to run DCM on the new system_from_graph DO NOT PR THIS BRANCH! It's very hacky at the moment to get things running quickly. It shouldn't be included until we have some in-depth conversations about how to combine hemodynamic observer blocks with neural mass blocks. --- src/Neuroblox.jl | 5 +- src/blox/connections.jl | 87 +++++++--- src/blox/neural_mass.jl | 2 +- src/measurementmodels/hemodynamic_extras.jl | 109 ++++++++++++ test/data_fitting_v2.jl | 176 ++++++++++++++++++++ test/jansen_rit_hemodynamic_tests.jl | 43 +++++ 6 files changed, 394 insertions(+), 28 deletions(-) create mode 100644 test/data_fitting_v2.jl create mode 100644 test/jansen_rit_hemodynamic_tests.jl diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index a22c21db..88b6d3ad 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -48,6 +48,8 @@ abstract type Merger end abstract type AbstractNeuronBlox <: AbstractBlox end abstract type NeuralMassBlox <: AbstractBlox end abstract type SuperBlox <: AbstractBlox end +abstract type ObserverBlox <: AbstractBlox end +abstract type CompoundNOBlox <: AbstractBlox end #I know this is bad - adding for speed of implementing -AGC # abstract type SourceBlox <: Blox end will be added later # we define these in neural_mass.jl @@ -170,5 +172,6 @@ export simulate, random_initials export system_from_graph, graph_delays export create_adjacency_edges! export get_namespaced_sys, namespace_expr, nameof -export HRFFourHalfCosine, HRFDoubleGamma +export HRFFourHalfCosine, HRFDoubleGamma, BalloonModel, CompoundHemo, AlternativeBalloonModel, LinHemoCombo +export input_equations, BloxConnector, accumulate_equation!, get_sys #remove this line end \ No newline at end of file diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 508a1a83..d4fadd38 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -33,36 +33,47 @@ function (bc::BloxConnector)( accumulate_equation!(bc, eq) end -# function (bc::BloxConnector)( -# jc::JansenRit, -# bloxin; -# weight = 1, -# delay = 0 -# ) -# # Need t for the delay term -# @variables t - -# sys_out = get_namespaced_sys(jc) -# sys_in = get_namespaced_sys(bloxin) - -# # Define & accumulate delay parameter -# τ_name = Symbol("τ_$(nameof(sys_out))_$(nameof(sys_in))") -# τ = only(@parameters $(τ_name)=delay) -# push!(bc.delays, τ) - -# w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") -# w = only(@parameters $(w_name)=weight) -# push!(bc.weights, w) - -# x = namespace_expr(jc.connector, sys_out, nameof(sys_out)) -# eq = sys_in.jcn ~ x(t-τ)*w +function (bc::BloxConnector)( + bloxout::NeuralMassBlox, + bloxin::NeuralMassBlox; + weight=1, + delay=0, + density=0.1 +) + # Need t for the delay term + @variables t + + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + if typeof(bloxout.output) == Num + w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") + w = only(@parameters $(w_name)=weight) + push!(bc.weights, w) + x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) + eq = sys_in.jcn ~ x*w + else + # Define & accumulate delay parameter + # Don't accumulate if zero + τ_name = Symbol("τ_$(nameof(sys_out))_$(nameof(sys_in))") + τ = only(@parameters $(τ_name)=delay) + push!(bc.delays, τ) + + w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") + w = only(@parameters $(w_name)=weight) + push!(bc.weights, w) + + x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) + eq = sys_in.jcn ~ x(t-τ)*w + end -# accumulate_equation!(bc, eq) -# end + accumulate_equation!(bc, eq) +end +# additional dispatch to connect to hemodynamic observer blox function (bc::BloxConnector)( bloxout::NeuralMassBlox, - bloxin::NeuralMassBlox; + bloxin::ObserverBlox; weight=1, delay=0, density=0.1 @@ -97,6 +108,30 @@ function (bc::BloxConnector)( accumulate_equation!(bc, eq) end +# Ok yes this is a bad dispatch but the whole compound blocks implementation is hacky and needs fixing @@ +# Opening an issue to loop back to this during clean up week +function (bc::BloxConnector)( + bloxout::CompoundNOBlox, + bloxin::CompoundNOBlox; + weight=1, + delay=0, + density=0.1 +) + # Need t for the delay term + @variables t + + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") + w = only(@parameters $(w_name)=weight) + push!(bc.weights, w) + x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) + eq = sys_in.nmm₊jcn ~ x*w + + accumulate_equation!(bc, eq) +end + function (bc::BloxConnector)( wta_out::WinnerTakeAllBlox, wta_in::WinnerTakeAllBlox; diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index dc1887e9..a3459099 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -242,7 +242,7 @@ struct LinearNeuralMass <: NeuralMassBlox odesystem namespace function LinearNeuralMass(;name) - sts = @variables x(t) [output=true] jcn(t) [input=true] + sts = @variables x(t)=0.0 [output=true] jcn(t)=0.0 [input=true] eqs = [D(x) ~ jcn] sys = System(eqs, name=name) new(sts[1], sts[2], sys, nothing) diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index 81fa9ad1..e3b68f21 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -1,5 +1,7 @@ using Random, SpecialFunctions, Plots +@parameters t +D = Differential(t) """ Woolrich, Behrens and Smith. 2003. @@ -62,3 +64,110 @@ function HRFDoubleGamma(; A=nothing, α₁=nothing, α₂=nothing, β₁=nothing hrf = A.*((((t.^(α₁-1)).*(β₁^α₁).*exp.(-β₁.*t))./gamma(α₁)) .- (c.*((t.^(α₂-1)).*(β₂^α₂).*exp.(-β₂.*t)./ gamma(α₂)))) end + + +""" +### Input variables ### +adaptation of the Hemodynamics blox in fmri.jl +""" +struct BalloonModel <: ObserverBlox + params + output + jcn + odesystem + namespace + function BalloonModel(;name, lnκ=0.0, lnτ=0.0) + #= hemodynamic parameters + H(1) - signal decay d(ds/dt)/ds) + H(2) - autoregulation d(ds/dt)/df) + H(3) - transit time (t0) + H(4) - exponent for Fout(v) (alpha) + H(5) - resting state oxygen extraction (E0) + =# + + H = [0.64, 0.32, 2.00, 0.32, 0.4] + p = progress_scope(@parameters lnκ=lnκ lnτ=lnτ) # progress scope if needed + #p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters + lnκ, lnτ = p # assign the modified parameters + + sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 lnq(t)=1.0 jcn(t)=0.0 [input=true] + + eqs = [ + D(s) ~ jcn/1e2 - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1), + D(lnf) ~ s / exp(lnf), + D(lnν) ~ (exp(lnf) - exp(lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(lnν)), + D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) + ] + sys = System(eqs, name=name) + new(p, Num(0), sts[5], sys, nothing) + end +end + + +# This doesn't work for some good reasons. You can't call system from graph and expect to call it again later +struct CompoundHemo <: CompoundNOBlox + params + output + jcn + odesystem + namespace + function CompoundHemo(massChoice; name, lnκ=0.0, lnτ=0.0) #ONLY WORKS WITH DEFAULT PARAMETERS FOR NOW + p = progress_scope(@parameters lnκ lnτ) + @named hemo = BalloonModel(;lnκ=p[1], lnτ=p[2]) + @named nmm = massChoice() + @variables jcn(t) + g = MetaDiGraph() + add_blox!(g, nmm) + add_blox!(g, hemo) + add_edge!(g, 1, 2, Dict(:weight => 1.0)) + linhemo = system_from_graph(g; name=name) + new(p, states(linhemo)[1], states(linhemo)[2], linhemo, nothing) + end +end + +struct LinHemoCombo <: CompoundNOBlox + params + output + jcn + odesystem + namespace + function LinHemoCombo(;name, lnκ=0.0, lnτ=0.0) + p = progress_scope(@parameters lnκ=lnκ lnτ=lnτ) + lnκ, lnτ = p # assign the modified parameters + H = [0.64, 0.32, 2.00, 0.32, 0.4] + + sts = @variables nmm₊x(t)=0.0 [output=true] nmm₊jcn(t)=0.0 [input=true] hemo₊s(t)=1.0 hemo₊lnf(t)=1.0 hemo₊lnν(t)=1.0 hemo₊lnq(t)=1.0 hemo₊jcn(t)=0.0 + eqs = [D(nmm₊x) ~ nmm₊jcn, + D(hemo₊s) ~ nmm₊x - H[1]*exp(lnκ)*hemo₊s - H[2]*(exp(hemo₊lnf) - 1), + D(hemo₊lnf) ~ hemo₊s / exp(hemo₊lnf), + D(hemo₊lnν) ~ (exp(hemo₊lnf) - exp(hemo₊lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(hemo₊lnν)), + D(hemo₊lnq) ~ (exp(hemo₊lnf)/exp(hemo₊lnq)*((1 - (1 - H[5])^(exp(hemo₊lnf)^-1))/H[5]) - exp(hemo₊lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) + ] + sys = System(eqs, name=name) + + new(p, sts[1], sts[2], sys, nothing) + end +end + + +struct AlternativeBalloonModel <: ObserverBlox + params + output + jcn + odesystem + namespace + function AlternativeBalloonModel(;name, κ=0.65, α=0.32, τ=0.98, ρ=0.34, V₀=0.02, γ=0.41) + p = progress_scope(@parameters κ=κ α=α τ=τ ρ=ρ V₀=V₀ γ=γ) # progress scope if needed + κ, α, τ, ρ, V₀, γ = p # assign the modified parameters + sts = @variables s(t)=1.0 f(t)=1.0 v(t)=1.0 q(t)=1.0 b(t)=0.0 [output=true] jcn(t)=0.0 [input=true] + eqs = [ + D(s) ~ jcn/1e2 - κ*s - γ*(f - 1), + D(f) ~ s, + D(v) ~ (f - v^(1/α))/τ, + D(q) ~ (((f*(1 - (1 - ρ^(1/f))/ρ)) - (v^(1/α)*q)/v))/τ, + b ~ V₀ * ((7*ρ*(1 - q)) + (2*(1-q/v)) + ((2*ρ - 0.2)*(1 - v))) + ] + sys = System(eqs, name=name) + new(p, Num(0), sts[5], sys, nothing) + end +end \ No newline at end of file diff --git a/test/data_fitting_v2.jl b/test/data_fitting_v2.jl new file mode 100644 index 00000000..d6c3f48f --- /dev/null +++ b/test/data_fitting_v2.jl @@ -0,0 +1,176 @@ +using Neuroblox, Test, Graphs, MetaGraphs, OrderedCollections, LinearAlgebra, DataFrames +using MAT + +### Load data ### +vars = matread(joinpath(@__DIR__, "spectralDCM_toydata.mat")); +data = DataFrame(vars["data"], :auto) # turn data into DataFrame +x = vars["x"] # initial conditions +nd = ncol(data) # number of dimensions + +########## assemble the model ########## + +@parameters κ=0.0 # define brain-wide decay parameter for hemodynamics +g = MetaDiGraph() +for ii = 1:nd + region = LinHemo(;name=Symbol("r$ii"), lnκ=κ) + add_blox!(g, region) +end + +# add symbolic weights +@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) +for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) + if idx[1] == idx[2] + add_edge!(g, idx[1], idx[2], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 + else + add_edge!(g, idx[1], idx[2], :weight, A[i]) + end +end +# compose model +@named neuronmodel = ODEfromGraph(g) +neuronmodel = structural_simplify(neuronmodel) + +# measurement model +@named bold = boldsignal() + +# attribute initial conditions to states +all_s = states(neuronmodel) +initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() +rnames = [] +map(x->push!(rnames, split(string(x), "₊")[1]), all_s); +rnames = unique(rnames); +for (i, r) in enumerate(rnames) + for (j, s) in enumerate(all_s[r .== map(x -> x[1], split.(string.(all_s), "₊"))]) # TODO: fix this solution, it is not robust!! + initcond[s] = x[i, j] + end +end + +modelparam = OrderedDict() +for par in parameters(neuronmodel) + # while Symbolics.getdefaultval(par) isa Num + # par = Symbolics.getdefaultval(par) + # end + modelparam[par] = Symbolics.getdefaultval(par) +end +# Noise parameter mean +modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 +modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above +modelparam[:lnγ] = zeros(Float64, nd); # region specific observation noise +modelparam[:C] = ones(Float64, nd); # C as in equation 3. NB: whatever C is defined to be here, it will be replaced in csd_approx. Another strange thing of SPM12... + +for par in parameters(bold) + modelparam[par] = Symbolics.getdefaultval(par) +end + +# define prior variances +paramvariance = copy(modelparam) +paramvariance[:C] = zeros(Float64, nd); +paramvariance[:lnγ] = ones(Float64, nd)./64.0; +paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; +paramvariance[:lnβ] = ones(Float64, length(modelparam[:lnβ]))./64.0; +for (k, v) in paramvariance + if occursin("A[", string(k)) + paramvariance[k] = vars["pC"][1, 1] + elseif occursin("κ", string(k)) + paramvariance[k] = ones(length(v))./256.0; + elseif occursin("ϵ", string(k)) + paramvariance[k] = 1/256.0; + elseif occursin("τ", string(k)) + paramvariance[k] = 1/256.0; + end +end + +params = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) +hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix + :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector + ) + +csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) +results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) + +### COMPARE RESULTS WITH MATLAB RESULTS ### +@show results.F, vars["F"] +@test results.F < vars["F"]*0.99 +@test results.F > vars["F"]*1.01 + +########## assemble the model ########## + +@parameters κ=0.0 # define brain-wide decay parameter for hemodynamics +g2 = MetaDiGraph() +for ii = 1:nd + region = LinHemoCombo(;name=Symbol("r$ii"), lnκ=κ) + add_blox!(g2, region) +end + +# add symbolic weights +@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) +for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) + if idx[1] == idx[2] + add_edge!(g2, idx[1], idx[2], Dict(:weight => -exp(A[i])/2)) # treatement of diagonal elements in SPM12 + else + add_edge!(g2, idx[1], idx[2], Dict(:weight => A[i])) + end +end + +# compose model +@named neuronmodel2 = system_from_graph(g2) +neuronmodel2 = structural_simplify(neuronmodel2) + +# attribute initial conditions to states +all_s = states(neuronmodel2) +initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() +rnames = [] +map(x->push!(rnames, split(string(x), "₊")[1]), all_s); +rnames = unique(rnames); +for (i, r) in enumerate(rnames) + for (j, s) in enumerate(all_s[r .== map(x -> x[1], split.(string.(all_s), "₊"))]) # TODO: fix this solution, it is not robust!! + initcond[s] = x[i, j] + end +end + +modelparam = OrderedDict() +for par in parameters(neuronmodel2) + # while Symbolics.getdefaultval(par) isa Num + # par = Symbolics.getdefaultval(par) + # end + modelparam[par] = Symbolics.getdefaultval(par) +end +# Noise parameter mean +modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 +modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above +modelparam[:lnγ] = zeros(Float64, nd); # region specific observation noise +modelparam[:C] = ones(Float64, nd); # C as in equation 3. NB: whatever C is defined to be here, it will be replaced in csd_approx. Another strange thing of SPM12... + +for par in parameters(bold) + modelparam[par] = Symbolics.getdefaultval(par) +end + +# define prior variances +paramvariance = copy(modelparam) +paramvariance[:C] = zeros(Float64, nd); +paramvariance[:lnγ] = ones(Float64, nd)./64.0; +paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; +paramvariance[:lnβ] = ones(Float64, length(modelparam[:lnβ]))./64.0; +for (k, v) in paramvariance + if occursin("A[", string(k)) + paramvariance[k] = vars["pC"][1, 1] + elseif occursin("κ", string(k)) + paramvariance[k] = ones(length(v))./256.0; + elseif occursin("ϵ", string(k)) + paramvariance[k] = 1/256.0; + elseif occursin("τ", string(k)) + paramvariance[k] = 1/256.0; + end +end + +params2 = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) +hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix + :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector + ) + +csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) +results = spectralVI(data, neuronmodel2, bold, initcond, csdsetup, params2, hyperparams) + +### COMPARE RESULTS WITH MATLAB RESULTS ### +@show results.F, vars["F"] +@test results.F < vars["F"]*0.99 +@test results.F > vars["F"]*1.01 \ No newline at end of file diff --git a/test/jansen_rit_hemodynamic_tests.jl b/test/jansen_rit_hemodynamic_tests.jl new file mode 100644 index 00000000..42c796e3 --- /dev/null +++ b/test/jansen_rit_hemodynamic_tests.jl @@ -0,0 +1,43 @@ +using Neuroblox, DifferentialEquations, DataFrames, Test, Distributions, Statistics, LinearAlgebra, Graphs, MetaGraphs, Random + +# Store parameters to be passed later on +params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5 + +adj_matrix_lin = [0 0 0 0 0 0 0 0 1; + -0.5*C_BG_Th -0.5*C_BG_Th C_BG_Th 0 0 0 0 0 0; + 0 -0.5*C_BG_Th 0 0 0 0 C_Cor_BG_Th 0 0; + 0 -0.5*C_BG_Th C_BG_Th 0 0 0 0 0 0; + 0 0 0 -0.5*C_BG_Th 0 0 0 0 0; + 0 0 0 0 C_BG_Th_Cor 0 6*C_Cor 0 0; + 0 0 0 0 0 4.8*C_Cor 0 -1.5*C_Cor 0; + 0 0 0 0 0 0 1.5*C_Cor 3.3*C_Cor 0; + 0 0 0 0 0 0 0 0 0] + +# test new Jansen-Rit blox +@named Str = JansenRit(τ=0.0022, H=20, λ=300, r=0.3) +@named GPe = JansenRit(τ=0.04, cortical=false) # all default subcortical except τ +@named STN = JansenRit(τ=0.01, H=20, λ=500, r=0.1) +@named GPi = JansenRit(cortical=false) # default parameters subcortical Jansen Rit blox +@named Th = JansenRit(τ=0.002, H=10, λ=20, r=5) +@named EI = JansenRit(τ=0.01, H=20, λ=5, r=5) +@named PY = JansenRit(cortical=true) # default parameters cortical Jansen Rit blox +@named II = JansenRit(τ=2.0, H=60, λ=5, r=5) +@named hemo = AlternativeBalloonModel() +blox = [Str, GPe, STN, GPi, Th, EI, PY, II, hemo] + +# Alternative version using adjacency matrix +g = MetaDiGraph() +add_blox!.(Ref(g), blox) +create_adjacency_edges!(g, adj_matrix_lin) + +@named final_system = system_from_graph(g, params) +final_delays = graph_delays(g) +sim_dur = 600.0 # Simulate for 10 Seconds +final_system_sys = structural_simplify(final_system) +prob = DDEProblem(final_system_sys, + [], + (0.0, sim_dur), + constant_lags = final_delays) +alg = MethodOfSteps(Vern7()) +sol_dde_no_delays = solve(prob, alg, saveat=0.001) +sol3 = DataFrame(sol_dde_no_delays) \ No newline at end of file From 07eff35d35c3087889cc6fc721ebeb9e4ac620d1 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 25 Sep 2023 14:44:59 -0400 Subject: [PATCH 06/28] Four half cosine parameterization Added for convolution option --- src/measurementmodels/hemodynamic_extras.jl | 39 +++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 src/measurementmodels/hemodynamic_extras.jl diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl new file mode 100644 index 00000000..79e3da5a --- /dev/null +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -0,0 +1,39 @@ +using Random, Plots + + +""" +Woolrich, Behrens and Smith. 2003. + +Specify widths (h₁, h₂, h₃, h₄) in seconds. +Amplitudes (f₁, f₂) relative to HRF height. +Output is HRF kernel in ms. + +""" + +function HRFFourHalfCosine(; h₁=nothing, h₂=nothing, h₃=nothing, h₄=nothing, f₁=nothing, f₂=nothing) + + h₁ = isnothing(h₁) ? 2*rand() : h₁ + h₂ = isnothing(h₂) ? 4*rand()+2 : h₂ + h₃ = isnothing(h₃) ? 4*rand()+2 : h₃ + h₄ = isnothing(h₄) ? 6*rand()+2 : h₄ + f₁ = isnothing(f₁) ? 0 : f₁ + f₂ = isnothing(f₂) ? 0.5*rand() : f₂ + + t₁ = 0:0.001:h₁ + t₂ = 0:0.001:h₂ + t₃ = 0:0.001:h₃ + t₄ = 0:0.001:h₄ + + cos₁ = 0.5.*f₁.*(cos.(π.*(t₁./h₁)) .- 1) + cos₂ = ((0.5*(f₁+1)).*(-cos.(π.*(t₂./h₂)) .+ 1)) .- f₁ + cos₃ = ((0.5*(f₂+1)).*(cos.(π.*(t₃./h₃)))) .+ ((0.5*(1-f₂))) + cos₄ = (0.5.*f₂.*(-cos.(π.*(t₄./h₄)) .+ 1)) .- f₂ + + hrf = vcat(cos₁, cos₂, cos₃, cos₄) + +end + + + +hmm = HRFFourHalfCosine(f₁=0) +plot(hmm) \ No newline at end of file From 439c14e473650dc5c5f3302a742c583f9deac509 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 25 Sep 2023 15:10:14 -0400 Subject: [PATCH 07/28] Double gamma added for convolution Same format as the four half cosine model --- Project.toml | 1 + src/measurementmodels/hemodynamic_extras.jl | 23 +++++++++++++++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index dfe5b72f..87166e71 100644 --- a/Project.toml +++ b/Project.toml @@ -46,6 +46,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SignalAnalysis = "df1fea92-c066-49dd-8b36-eace3378ea47" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index 79e3da5a..fc8bc709 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -1,4 +1,4 @@ -using Random, Plots +using Random, SpecialFunctions, Plots """ @@ -33,7 +33,22 @@ function HRFFourHalfCosine(; h₁=nothing, h₂=nothing, h₃=nothing, h₄=noth end - +function HRFDoubleGamma(; A=nothing, α₁=nothing, α₂=nothing, β₁=nothing, β₂=nothing, c=nothing, tlen=nothing) -hmm = HRFFourHalfCosine(f₁=0) -plot(hmm) \ No newline at end of file + if (!isnothing(A) || !isnothing(α₁) || !isnothing(α₂) || !isnothing(β₁) || !isnothing(β₂) || !isnothing(c)) && isnothing(tlen) + throw(DomainError(nothing, "If you specify a different parameter, you need to specify tlen as well to ensure the kernel is wide enough.")) + end + + A = isnothing(A) ? 1 : A + α₁ = isnothing(α₁) ? 6 : α₁ + α₂ = isnothing(α₂) ? 16 : α₂ + β₁ = isnothing(β₁) ? 1 : β₁ + β₂ = isnothing(β₂) ? 1 : β₂ + c = isnothing(c) ? 1.0/6.0 : c + tlen = isnothing(tlen) ? 30 : tlen + + t = 0:0.001:tlen + + hrf = A.*((((t.^(α₁-1)).*(β₁^α₁).*exp.(-β₁.*t))./gamma(α₁)) .- (c.*((t.^(α₂-1)).*(β₂^α₂).*exp.(-β₂.*t)./ gamma(α₂)))) + +end From d382fa0dc7495943ac78e9a8c9dd17f7fc6050b4 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 25 Sep 2023 17:25:46 -0400 Subject: [PATCH 08/28] Citation details --- src/measurementmodels/hemodynamic_extras.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index fc8bc709..81fa9ad1 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -33,6 +33,16 @@ function HRFFourHalfCosine(; h₁=nothing, h₂=nothing, h₃=nothing, h₄=noth end +""" +Lindquist et al. 2012. + +All parameters are in seconds. +Output is in ms. +If you specify a change in any of the parameters, also change the tlen to ensure the kernel is wide enough. +It will throw an error if you don't. + +""" + function HRFDoubleGamma(; A=nothing, α₁=nothing, α₂=nothing, β₁=nothing, β₂=nothing, c=nothing, tlen=nothing) if (!isnothing(A) || !isnothing(α₁) || !isnothing(α₂) || !isnothing(β₁) || !isnothing(β₂) || !isnothing(c)) && isnothing(tlen) From a8d30c0625d4de37effce145390a17fbd7a1206a Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Tue, 26 Sep 2023 11:14:50 -0400 Subject: [PATCH 09/28] Minor testing updates + exports Cleaning up to prepare for more fMRI simulations --- src/Neuroblox.jl | 2 ++ test/jansen_rit_component_tests.jl | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 016db35e..42a8d26e 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -101,6 +101,7 @@ include("gui/GUI.jl") include("blox/connections.jl") include("blox/blox_utilities.jl") include("Neurographs.jl") +include("measurementmodels/hemodynamic_extras.jl") function simulate(sys::ODESystem, u0, timespan, p, solver = AutoVern7(Rodas4()); kwargs...) prob = ODEProblem(sys, u0, timespan, p) @@ -180,4 +181,5 @@ export create_adjacency_edges! export get_namespaced_sys, nameof export run_experiment! +export HRFFourHalfCosine, HRFDoubleGamma end \ No newline at end of file diff --git a/test/jansen_rit_component_tests.jl b/test/jansen_rit_component_tests.jl index 0f0d599e..8f77e5c0 100644 --- a/test/jansen_rit_component_tests.jl +++ b/test/jansen_rit_component_tests.jl @@ -118,7 +118,7 @@ create_adjacency_edges!(g2, adj_matrix_lin) @named final_system = system_from_graph(g2, params) final_delays = graph_delays(g2) -sim_dur = 10.0 # Simulate for 10 Seconds +sim_dur = 600.0 # Simulate for 10 Seconds final_system_sys = structural_simplify(final_system) prob = DDEProblem(final_system_sys, [], From 8fe5ec47f6496b243c48380d6b7944949cd598a8 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Wed, 27 Sep 2023 19:41:19 -0400 Subject: [PATCH 10/28] Updates to run DCM on the new system_from_graph DO NOT PR THIS BRANCH! It's very hacky at the moment to get things running quickly. It shouldn't be included until we have some in-depth conversations about how to combine hemodynamic observer blocks with neural mass blocks. --- src/Neuroblox.jl | 6 +- src/blox/connections.jl | 62 +++++++ src/blox/neural_mass.jl | 4 +- src/measurementmodels/hemodynamic_extras.jl | 109 ++++++++++++ test/data_fitting_v2.jl | 176 ++++++++++++++++++++ test/jansen_rit_hemodynamic_tests.jl | 43 +++++ 6 files changed, 395 insertions(+), 5 deletions(-) create mode 100644 test/data_fitting_v2.jl create mode 100644 test/jansen_rit_hemodynamic_tests.jl diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 42a8d26e..4a9954a6 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -51,6 +51,8 @@ abstract type Merger end abstract type AbstractNeuronBlox <: AbstractBlox end abstract type NeuralMassBlox <: AbstractBlox end abstract type SuperBlox <: AbstractBlox end +abstract type ObserverBlox <: AbstractBlox end +abstract type CompoundNOBlox <: AbstractBlox end #I know this is bad - adding for speed of implementing -AGC abstract type StimulusBlox <: AbstractBlox end # we define these in neural_mass.jl @@ -178,8 +180,6 @@ export vecparam, unvecparam, csd_Q, spectralVI export simulate, random_initials export system_from_graph, graph_delays export create_adjacency_edges! -export get_namespaced_sys, nameof -export run_experiment! - +export get_namespaced_sys, namespace_expr, nameof export HRFFourHalfCosine, HRFDoubleGamma end \ No newline at end of file diff --git a/src/blox/connections.jl b/src/blox/connections.jl index d6836cbc..ef0101d0 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -155,6 +155,68 @@ function (bc::BloxConnector)( accumulate_equation!(bc, eq) end +# additional dispatch to connect to hemodynamic observer blox +function (bc::BloxConnector)( + bloxout::NeuralMassBlox, + bloxin::ObserverBlox; + weight=1, + delay=0, + density=0.1 +) + # Need t for the delay term + @variables t + + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + if typeof(bloxout.output) == Num + w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") + w = only(@parameters $(w_name)=weight) + push!(bc.weights, w) + x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) + eq = sys_in.jcn ~ x*w + else + # Define & accumulate delay parameter + # Don't accumulate if zero + τ_name = Symbol("τ_$(nameof(sys_out))_$(nameof(sys_in))") + τ = only(@parameters $(τ_name)=delay) + push!(bc.delays, τ) + + w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") + w = only(@parameters $(w_name)=weight) + push!(bc.weights, w) + + x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) + eq = sys_in.jcn ~ x(t-τ)*w + end + + accumulate_equation!(bc, eq) +end + +# Ok yes this is a bad dispatch but the whole compound blocks implementation is hacky and needs fixing @@ +# Opening an issue to loop back to this during clean up week +function (bc::BloxConnector)( + bloxout::CompoundNOBlox, + bloxin::CompoundNOBlox; + weight=1, + delay=0, + density=0.1 +) + # Need t for the delay term + @variables t + + sys_out = get_namespaced_sys(bloxout) + sys_in = get_namespaced_sys(bloxin) + + w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") + w = only(@parameters $(w_name)=weight) + push!(bc.weights, w) + x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) + eq = sys_in.nmm₊jcn ~ x*w + + accumulate_equation!(bc, eq) +end + function (bc::BloxConnector)( wta_out::WinnerTakeAllBlox, wta_in::WinnerTakeAllBlox; diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 834367fb..3f8d9b82 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -300,8 +300,8 @@ struct LinearNeuralMass <: NeuralMassBlox jcn odesystem namespace - function LinearNeuralMass(;name, namespace=nothing) - sts = @variables x(t) [output=true] jcn(t) [input=true] + function LinearNeuralMass(;name) + sts = @variables x(t)=0.0 [output=true] jcn(t)=0.0 [input=true] eqs = [D(x) ~ jcn] sys = System(eqs, name=name) new(sts[1], sts[2], sys, namespace) diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index 81fa9ad1..e3b68f21 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -1,5 +1,7 @@ using Random, SpecialFunctions, Plots +@parameters t +D = Differential(t) """ Woolrich, Behrens and Smith. 2003. @@ -62,3 +64,110 @@ function HRFDoubleGamma(; A=nothing, α₁=nothing, α₂=nothing, β₁=nothing hrf = A.*((((t.^(α₁-1)).*(β₁^α₁).*exp.(-β₁.*t))./gamma(α₁)) .- (c.*((t.^(α₂-1)).*(β₂^α₂).*exp.(-β₂.*t)./ gamma(α₂)))) end + + +""" +### Input variables ### +adaptation of the Hemodynamics blox in fmri.jl +""" +struct BalloonModel <: ObserverBlox + params + output + jcn + odesystem + namespace + function BalloonModel(;name, lnκ=0.0, lnτ=0.0) + #= hemodynamic parameters + H(1) - signal decay d(ds/dt)/ds) + H(2) - autoregulation d(ds/dt)/df) + H(3) - transit time (t0) + H(4) - exponent for Fout(v) (alpha) + H(5) - resting state oxygen extraction (E0) + =# + + H = [0.64, 0.32, 2.00, 0.32, 0.4] + p = progress_scope(@parameters lnκ=lnκ lnτ=lnτ) # progress scope if needed + #p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters + lnκ, lnτ = p # assign the modified parameters + + sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 lnq(t)=1.0 jcn(t)=0.0 [input=true] + + eqs = [ + D(s) ~ jcn/1e2 - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1), + D(lnf) ~ s / exp(lnf), + D(lnν) ~ (exp(lnf) - exp(lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(lnν)), + D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) + ] + sys = System(eqs, name=name) + new(p, Num(0), sts[5], sys, nothing) + end +end + + +# This doesn't work for some good reasons. You can't call system from graph and expect to call it again later +struct CompoundHemo <: CompoundNOBlox + params + output + jcn + odesystem + namespace + function CompoundHemo(massChoice; name, lnκ=0.0, lnτ=0.0) #ONLY WORKS WITH DEFAULT PARAMETERS FOR NOW + p = progress_scope(@parameters lnκ lnτ) + @named hemo = BalloonModel(;lnκ=p[1], lnτ=p[2]) + @named nmm = massChoice() + @variables jcn(t) + g = MetaDiGraph() + add_blox!(g, nmm) + add_blox!(g, hemo) + add_edge!(g, 1, 2, Dict(:weight => 1.0)) + linhemo = system_from_graph(g; name=name) + new(p, states(linhemo)[1], states(linhemo)[2], linhemo, nothing) + end +end + +struct LinHemoCombo <: CompoundNOBlox + params + output + jcn + odesystem + namespace + function LinHemoCombo(;name, lnκ=0.0, lnτ=0.0) + p = progress_scope(@parameters lnκ=lnκ lnτ=lnτ) + lnκ, lnτ = p # assign the modified parameters + H = [0.64, 0.32, 2.00, 0.32, 0.4] + + sts = @variables nmm₊x(t)=0.0 [output=true] nmm₊jcn(t)=0.0 [input=true] hemo₊s(t)=1.0 hemo₊lnf(t)=1.0 hemo₊lnν(t)=1.0 hemo₊lnq(t)=1.0 hemo₊jcn(t)=0.0 + eqs = [D(nmm₊x) ~ nmm₊jcn, + D(hemo₊s) ~ nmm₊x - H[1]*exp(lnκ)*hemo₊s - H[2]*(exp(hemo₊lnf) - 1), + D(hemo₊lnf) ~ hemo₊s / exp(hemo₊lnf), + D(hemo₊lnν) ~ (exp(hemo₊lnf) - exp(hemo₊lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(hemo₊lnν)), + D(hemo₊lnq) ~ (exp(hemo₊lnf)/exp(hemo₊lnq)*((1 - (1 - H[5])^(exp(hemo₊lnf)^-1))/H[5]) - exp(hemo₊lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) + ] + sys = System(eqs, name=name) + + new(p, sts[1], sts[2], sys, nothing) + end +end + + +struct AlternativeBalloonModel <: ObserverBlox + params + output + jcn + odesystem + namespace + function AlternativeBalloonModel(;name, κ=0.65, α=0.32, τ=0.98, ρ=0.34, V₀=0.02, γ=0.41) + p = progress_scope(@parameters κ=κ α=α τ=τ ρ=ρ V₀=V₀ γ=γ) # progress scope if needed + κ, α, τ, ρ, V₀, γ = p # assign the modified parameters + sts = @variables s(t)=1.0 f(t)=1.0 v(t)=1.0 q(t)=1.0 b(t)=0.0 [output=true] jcn(t)=0.0 [input=true] + eqs = [ + D(s) ~ jcn/1e2 - κ*s - γ*(f - 1), + D(f) ~ s, + D(v) ~ (f - v^(1/α))/τ, + D(q) ~ (((f*(1 - (1 - ρ^(1/f))/ρ)) - (v^(1/α)*q)/v))/τ, + b ~ V₀ * ((7*ρ*(1 - q)) + (2*(1-q/v)) + ((2*ρ - 0.2)*(1 - v))) + ] + sys = System(eqs, name=name) + new(p, Num(0), sts[5], sys, nothing) + end +end \ No newline at end of file diff --git a/test/data_fitting_v2.jl b/test/data_fitting_v2.jl new file mode 100644 index 00000000..d6c3f48f --- /dev/null +++ b/test/data_fitting_v2.jl @@ -0,0 +1,176 @@ +using Neuroblox, Test, Graphs, MetaGraphs, OrderedCollections, LinearAlgebra, DataFrames +using MAT + +### Load data ### +vars = matread(joinpath(@__DIR__, "spectralDCM_toydata.mat")); +data = DataFrame(vars["data"], :auto) # turn data into DataFrame +x = vars["x"] # initial conditions +nd = ncol(data) # number of dimensions + +########## assemble the model ########## + +@parameters κ=0.0 # define brain-wide decay parameter for hemodynamics +g = MetaDiGraph() +for ii = 1:nd + region = LinHemo(;name=Symbol("r$ii"), lnκ=κ) + add_blox!(g, region) +end + +# add symbolic weights +@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) +for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) + if idx[1] == idx[2] + add_edge!(g, idx[1], idx[2], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 + else + add_edge!(g, idx[1], idx[2], :weight, A[i]) + end +end +# compose model +@named neuronmodel = ODEfromGraph(g) +neuronmodel = structural_simplify(neuronmodel) + +# measurement model +@named bold = boldsignal() + +# attribute initial conditions to states +all_s = states(neuronmodel) +initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() +rnames = [] +map(x->push!(rnames, split(string(x), "₊")[1]), all_s); +rnames = unique(rnames); +for (i, r) in enumerate(rnames) + for (j, s) in enumerate(all_s[r .== map(x -> x[1], split.(string.(all_s), "₊"))]) # TODO: fix this solution, it is not robust!! + initcond[s] = x[i, j] + end +end + +modelparam = OrderedDict() +for par in parameters(neuronmodel) + # while Symbolics.getdefaultval(par) isa Num + # par = Symbolics.getdefaultval(par) + # end + modelparam[par] = Symbolics.getdefaultval(par) +end +# Noise parameter mean +modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 +modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above +modelparam[:lnγ] = zeros(Float64, nd); # region specific observation noise +modelparam[:C] = ones(Float64, nd); # C as in equation 3. NB: whatever C is defined to be here, it will be replaced in csd_approx. Another strange thing of SPM12... + +for par in parameters(bold) + modelparam[par] = Symbolics.getdefaultval(par) +end + +# define prior variances +paramvariance = copy(modelparam) +paramvariance[:C] = zeros(Float64, nd); +paramvariance[:lnγ] = ones(Float64, nd)./64.0; +paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; +paramvariance[:lnβ] = ones(Float64, length(modelparam[:lnβ]))./64.0; +for (k, v) in paramvariance + if occursin("A[", string(k)) + paramvariance[k] = vars["pC"][1, 1] + elseif occursin("κ", string(k)) + paramvariance[k] = ones(length(v))./256.0; + elseif occursin("ϵ", string(k)) + paramvariance[k] = 1/256.0; + elseif occursin("τ", string(k)) + paramvariance[k] = 1/256.0; + end +end + +params = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) +hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix + :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector + ) + +csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) +results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) + +### COMPARE RESULTS WITH MATLAB RESULTS ### +@show results.F, vars["F"] +@test results.F < vars["F"]*0.99 +@test results.F > vars["F"]*1.01 + +########## assemble the model ########## + +@parameters κ=0.0 # define brain-wide decay parameter for hemodynamics +g2 = MetaDiGraph() +for ii = 1:nd + region = LinHemoCombo(;name=Symbol("r$ii"), lnκ=κ) + add_blox!(g2, region) +end + +# add symbolic weights +@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) +for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) + if idx[1] == idx[2] + add_edge!(g2, idx[1], idx[2], Dict(:weight => -exp(A[i])/2)) # treatement of diagonal elements in SPM12 + else + add_edge!(g2, idx[1], idx[2], Dict(:weight => A[i])) + end +end + +# compose model +@named neuronmodel2 = system_from_graph(g2) +neuronmodel2 = structural_simplify(neuronmodel2) + +# attribute initial conditions to states +all_s = states(neuronmodel2) +initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() +rnames = [] +map(x->push!(rnames, split(string(x), "₊")[1]), all_s); +rnames = unique(rnames); +for (i, r) in enumerate(rnames) + for (j, s) in enumerate(all_s[r .== map(x -> x[1], split.(string.(all_s), "₊"))]) # TODO: fix this solution, it is not robust!! + initcond[s] = x[i, j] + end +end + +modelparam = OrderedDict() +for par in parameters(neuronmodel2) + # while Symbolics.getdefaultval(par) isa Num + # par = Symbolics.getdefaultval(par) + # end + modelparam[par] = Symbolics.getdefaultval(par) +end +# Noise parameter mean +modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 +modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above +modelparam[:lnγ] = zeros(Float64, nd); # region specific observation noise +modelparam[:C] = ones(Float64, nd); # C as in equation 3. NB: whatever C is defined to be here, it will be replaced in csd_approx. Another strange thing of SPM12... + +for par in parameters(bold) + modelparam[par] = Symbolics.getdefaultval(par) +end + +# define prior variances +paramvariance = copy(modelparam) +paramvariance[:C] = zeros(Float64, nd); +paramvariance[:lnγ] = ones(Float64, nd)./64.0; +paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; +paramvariance[:lnβ] = ones(Float64, length(modelparam[:lnβ]))./64.0; +for (k, v) in paramvariance + if occursin("A[", string(k)) + paramvariance[k] = vars["pC"][1, 1] + elseif occursin("κ", string(k)) + paramvariance[k] = ones(length(v))./256.0; + elseif occursin("ϵ", string(k)) + paramvariance[k] = 1/256.0; + elseif occursin("τ", string(k)) + paramvariance[k] = 1/256.0; + end +end + +params2 = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) +hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix + :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector + ) + +csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) +results = spectralVI(data, neuronmodel2, bold, initcond, csdsetup, params2, hyperparams) + +### COMPARE RESULTS WITH MATLAB RESULTS ### +@show results.F, vars["F"] +@test results.F < vars["F"]*0.99 +@test results.F > vars["F"]*1.01 \ No newline at end of file diff --git a/test/jansen_rit_hemodynamic_tests.jl b/test/jansen_rit_hemodynamic_tests.jl new file mode 100644 index 00000000..42c796e3 --- /dev/null +++ b/test/jansen_rit_hemodynamic_tests.jl @@ -0,0 +1,43 @@ +using Neuroblox, DifferentialEquations, DataFrames, Test, Distributions, Statistics, LinearAlgebra, Graphs, MetaGraphs, Random + +# Store parameters to be passed later on +params = @parameters C_Cor=60 C_BG_Th=60 C_Cor_BG_Th=5 C_BG_Th_Cor=5 + +adj_matrix_lin = [0 0 0 0 0 0 0 0 1; + -0.5*C_BG_Th -0.5*C_BG_Th C_BG_Th 0 0 0 0 0 0; + 0 -0.5*C_BG_Th 0 0 0 0 C_Cor_BG_Th 0 0; + 0 -0.5*C_BG_Th C_BG_Th 0 0 0 0 0 0; + 0 0 0 -0.5*C_BG_Th 0 0 0 0 0; + 0 0 0 0 C_BG_Th_Cor 0 6*C_Cor 0 0; + 0 0 0 0 0 4.8*C_Cor 0 -1.5*C_Cor 0; + 0 0 0 0 0 0 1.5*C_Cor 3.3*C_Cor 0; + 0 0 0 0 0 0 0 0 0] + +# test new Jansen-Rit blox +@named Str = JansenRit(τ=0.0022, H=20, λ=300, r=0.3) +@named GPe = JansenRit(τ=0.04, cortical=false) # all default subcortical except τ +@named STN = JansenRit(τ=0.01, H=20, λ=500, r=0.1) +@named GPi = JansenRit(cortical=false) # default parameters subcortical Jansen Rit blox +@named Th = JansenRit(τ=0.002, H=10, λ=20, r=5) +@named EI = JansenRit(τ=0.01, H=20, λ=5, r=5) +@named PY = JansenRit(cortical=true) # default parameters cortical Jansen Rit blox +@named II = JansenRit(τ=2.0, H=60, λ=5, r=5) +@named hemo = AlternativeBalloonModel() +blox = [Str, GPe, STN, GPi, Th, EI, PY, II, hemo] + +# Alternative version using adjacency matrix +g = MetaDiGraph() +add_blox!.(Ref(g), blox) +create_adjacency_edges!(g, adj_matrix_lin) + +@named final_system = system_from_graph(g, params) +final_delays = graph_delays(g) +sim_dur = 600.0 # Simulate for 10 Seconds +final_system_sys = structural_simplify(final_system) +prob = DDEProblem(final_system_sys, + [], + (0.0, sim_dur), + constant_lags = final_delays) +alg = MethodOfSteps(Vern7()) +sol_dde_no_delays = solve(prob, alg, saveat=0.001) +sol3 = DataFrame(sol_dde_no_delays) \ No newline at end of file From e5225ef16e6fe0c2042bd4748ec669b7f1f6d174 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Mon, 2 Oct 2023 16:33:03 -0400 Subject: [PATCH 11/28] Again DO NOT PR THIS BRANCH Just committing so that I have a backup as I start messing with some other things - this is still NOT ready for merging so needs to remain quarantined. For reference: This includes some modular implementations to extend the fitting procedures to blocks beyond LinearNeuralMass. As it is, there's hard-coded state number issues with the old implementation, and hierarchy issues with the `system_from_graph` implementation. So we'll need to talk through some design choices after the deadline is met for other stuff. --- src/Neuroblox.jl | 3 +- src/blox/connections.jl | 8 +- src/datafitting/spectralDCM.jl | 54 ++++++++++ src/measurementmodels/hemodynamic_extras.jl | 19 ++++ test/data_fitting_v2.jl | 108 +++++++++++++++++--- 5 files changed, 176 insertions(+), 16 deletions(-) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 4a9954a6..95d387fc 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -181,5 +181,6 @@ export simulate, random_initials export system_from_graph, graph_delays export create_adjacency_edges! export get_namespaced_sys, namespace_expr, nameof -export HRFFourHalfCosine, HRFDoubleGamma +export HRFFourHalfCosine, HRFDoubleGamma, BalloonModel, CompoundHemo, AlternativeBalloonModel, LinHemoCombo, JRHemo, spectralVI2 +export input_equations, BloxConnector, accumulate_equation!, get_sys #remove this line end \ No newline at end of file diff --git a/src/blox/connections.jl b/src/blox/connections.jl index ef0101d0..8c2d0275 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -202,14 +202,16 @@ function (bc::BloxConnector)( delay=0, density=0.1 ) - # Need t for the delay term - @variables t sys_out = get_namespaced_sys(bloxout) sys_in = get_namespaced_sys(bloxin) w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") - w = only(@parameters $(w_name)=weight) + if typeof(weight) == Symbol + w = weight + else + w = only(@parameters $(w_name)=weight) + end push!(bc.weights, w) x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) eq = sys_in.nmm₊jcn ~ x*w diff --git a/src/datafitting/spectralDCM.jl b/src/datafitting/spectralDCM.jl index 1f4ca755..e6e8d147 100644 --- a/src/datafitting/spectralDCM.jl +++ b/src/datafitting/spectralDCM.jl @@ -511,6 +511,60 @@ function spectralVI(data, neuraldynmodel, observationmodel, initcond, csdsetup, idx_A = findall(occursin.("A[", string.(calculate_jacobian(neuraldynmodel)))) + ### Compute the variational Bayes with Laplace approximation ### + return variationalbayes(idx_A, y_csd, derivatives, freqs, V, p, priors, 128) +end + +function spectralVI2(data, neuraldynmodel, observationmodel, initcond, csdsetup, params, hyperparams) + # compute cross-spectral density + y = Matrix(data); + nd = ncol(data); # dimension of the data, number of regions + dt = csdsetup[:dt]; # order of MAR. Hard-coded in SPM12 with this value. We will use the same for now. + freqs = csdsetup[:freq]; # frequencies at which the CSD is evaluated + p = csdsetup[:p]; # order of MAR + mar = mar_ml(y, p); # compute MAR from time series y and model order p + y_csd = mar2csd(mar, freqs, dt^-1); # compute cross spectral densities from MAR parameters at specific frequencies freqs, dt^-1 is sampling rate of data + + grad_full = function(p, grad, sts, nd) + tmp = zeros(typeof(p), nd, length(sts)) + for i in 1:nd + tmp[i, (i-1)*6 .+ (5:6)] = grad(vcat([1], sts[(i-1)*6 .+ (6:-1:5)]), p, t)[4:-1:3] + end + return tmp + end + jac_f = generate_jacobian(neuraldynmodel, expression = Val{false})[1] + grad_g = generate_jacobian(observationmodel, expression = Val{false})[1] + + statevals = [v for v in values(initcond)] + derivatives = Dict(:∂f => par -> jac_f(statevals, par, t), + :∂g => par -> grad_full(par, grad_g, statevals, nd)) + + θΣ = diagm(vecparam(OrderedDict(params.name .=> params.variance))) + # depending on the definition of the priors (note that we take it from the SPM12 code), some dimensions are set to 0 and thus are not changed. + # Extract these dimensions and remove them from the remaining computation. I find this a bit odd and further thoughts would be necessary to understand + # to what extend this is a the most reasonable approach. + idx = findall(x -> x != 0, θΣ); + V = zeros(size(θΣ, 1), length(idx)); + order = sortperm(θΣ[idx], rev=true); + idx = idx[order]; + for i = 1:length(idx) + V[idx[i][1], i] = 1.0 + end + θΣ = V'*θΣ*V; # reduce dimension by removing columns and rows that are all 0 + + ### Collect prior means and covariances ### + Q = csd_Q(y_csd); # compute prior of Q, the precision (of the data) components. See Friston etal. 2007 Appendix A + priors = Dict(:μ => OrderedDict(params.name .=> params.mean), + :Σ => Dict( + :Πθ_pr => inv(θΣ), # prior model parameter precision + :Πλ_pr => hyperparams[:Πλ_pr], # prior metaparameter precision + :μλ_pr => hyperparams[:μλ_pr], # prior metaparameter mean + :Q => Q # decomposition of model parameter covariance + ) + ); + + idx_A = findall(occursin.("A[", string.(calculate_jacobian(neuraldynmodel)))) + ### Compute the variational Bayes with Laplace approximation ### return variationalbayes(idx_A, y_csd, derivatives, freqs, V, p, priors, 128) end \ No newline at end of file diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index e3b68f21..38c39881 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -170,4 +170,23 @@ struct AlternativeBalloonModel <: ObserverBlox sys = System(eqs, name=name) new(p, Num(0), sts[5], sys, nothing) end +end + +mutable struct JRHemo <: AbstractComponent + connector::Num + bloxinput::Num + odesystem::ODESystem + function JRHemo(;name, lnκ=0.0, lnτ=0.0) + params = progress_scope(lnκ, lnτ) + @named hemo = Hemodynamics(;lnκ=params[1], lnτ=params[2]) + @named nmm = JansenRitCBlox() + @variables jcn(t) + + g = MetaDiGraph() + add_vertex!(g, Dict(:blox => nmm, :jcn => jcn)) + add_vertex!(g, :blox, hemo) + add_edge!(g, 1, 2, :weight, 1.0) + linhemo = ODEfromGraph(g; name=name) + new(linhemo.nmm₊x, linhemo.jcn, linhemo) + end end \ No newline at end of file diff --git a/test/data_fitting_v2.jl b/test/data_fitting_v2.jl index d6c3f48f..4a1f7390 100644 --- a/test/data_fitting_v2.jl +++ b/test/data_fitting_v2.jl @@ -94,31 +94,115 @@ results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperp ########## assemble the model ########## +# @parameters κ=0.0 # define brain-wide decay parameter for hemodynamics +# g2 = MetaDiGraph() +# for ii = 1:nd +# region = LinHemoCombo(;name=Symbol("r$ii"), lnκ=κ) +# add_blox!(g2, region) +# end + +# # add symbolic weights +# @parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) +# for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) +# if idx[1] == idx[2] +# add_edge!(g2, idx[1], idx[2], Dict(:weight => -exp(A[i])/2)) # treatement of diagonal elements in SPM12 +# else +# add_edge!(g2, idx[1], idx[2], Dict(:weight => A[i])) +# end +# end + +# # compose model +# @named neuronmodel2 = system_from_graph(g2) +# neuronmodel2 = structural_simplify(neuronmodel2) + +# # attribute initial conditions to states +# all_s = states(neuronmodel2) +# initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() +# rnames = [] +# map(x->push!(rnames, split(string(x), "₊")[1]), all_s); +# rnames = unique(rnames); +# for (i, r) in enumerate(rnames) +# for (j, s) in enumerate(all_s[r .== map(x -> x[1], split.(string.(all_s), "₊"))]) # TODO: fix this solution, it is not robust!! +# initcond[s] = x[i, j] +# end +# end + +# modelparam = OrderedDict() +# for par in parameters(neuronmodel2) +# while Symbolics.getdefaultval(par) isa Num +# par = Symbolics.getdefaultval(par) +# end +# modelparam[par] = Symbolics.getdefaultval(par) +# end +# # Noise parameter mean +# modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 +# modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above +# modelparam[:lnγ] = zeros(Float64, nd); # region specific observation noise +# modelparam[:C] = ones(Float64, nd); # C as in equation 3. NB: whatever C is defined to be here, it will be replaced in csd_approx. Another strange thing of SPM12... + +# for par in parameters(bold) +# modelparam[par] = Symbolics.getdefaultval(par) +# end + +# # define prior variances +# paramvariance = copy(modelparam) +# paramvariance[:C] = zeros(Float64, nd); +# paramvariance[:lnγ] = ones(Float64, nd)./64.0; +# paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; +# paramvariance[:lnβ] = ones(Float64, length(modelparam[:lnβ]))./64.0; +# for (k, v) in paramvariance +# if occursin("A[", string(k)) +# paramvariance[k] = vars["pC"][1, 1] +# elseif occursin("κ", string(k)) +# paramvariance[k] = ones(length(v))./256.0; +# elseif occursin("ϵ", string(k)) +# paramvariance[k] = 1/256.0; +# elseif occursin("τ", string(k)) +# paramvariance[k] = 1/256.0; +# end +# end + +# params2 = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) +# hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix +# :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector +# ) + +# csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) +# results = spectralVI(data, neuronmodel2, bold, initcond, csdsetup, params2, hyperparams) + +# ### COMPARE RESULTS WITH MATLAB RESULTS ### +# @show results.F, vars["F"] +# @test results.F < vars["F"]*0.99 +# @test results.F > vars["F"]*1.01 + @parameters κ=0.0 # define brain-wide decay parameter for hemodynamics -g2 = MetaDiGraph() +g = MetaDiGraph() for ii = 1:nd - region = LinHemoCombo(;name=Symbol("r$ii"), lnκ=κ) - add_blox!(g2, region) + region = JRHemo(;name=Symbol("r$ii"), lnκ=κ) + add_blox!(g, region) end # add symbolic weights @parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) if idx[1] == idx[2] - add_edge!(g2, idx[1], idx[2], Dict(:weight => -exp(A[i])/2)) # treatement of diagonal elements in SPM12 + add_edge!(g, idx[1], idx[2], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 else - add_edge!(g2, idx[1], idx[2], Dict(:weight => A[i])) + add_edge!(g, idx[1], idx[2], :weight, A[i]) end end - # compose model -@named neuronmodel2 = system_from_graph(g2) -neuronmodel2 = structural_simplify(neuronmodel2) +@named neuronmodel = ODEfromGraph(g) +neuronmodel = structural_simplify(neuronmodel) + +# measurement model +@named bold = boldsignal() # attribute initial conditions to states -all_s = states(neuronmodel2) +all_s = states(neuronmodel) initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() rnames = [] +x = rand(3, 6) #WHY IS x 3x5? map(x->push!(rnames, split(string(x), "₊")[1]), all_s); rnames = unique(rnames); for (i, r) in enumerate(rnames) @@ -128,7 +212,7 @@ for (i, r) in enumerate(rnames) end modelparam = OrderedDict() -for par in parameters(neuronmodel2) +for par in parameters(neuronmodel) # while Symbolics.getdefaultval(par) isa Num # par = Symbolics.getdefaultval(par) # end @@ -162,13 +246,13 @@ for (k, v) in paramvariance end end -params2 = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) +params = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector ) csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) -results = spectralVI(data, neuronmodel2, bold, initcond, csdsetup, params2, hyperparams) +results = spectralVI2(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) ### COMPARE RESULTS WITH MATLAB RESULTS ### @show results.F, vars["F"] From 2e5f581957b95eaf5b36c001572963e2f2fff935 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Wed, 8 Nov 2023 17:05:47 +0100 Subject: [PATCH 12/28] WIP: David's integration into the new system_from_graph structure. --- src/blox/blox_utilities.jl | 2 +- src/blox/connections.jl | 12 +++++++++--- src/datafitting/spectralDCM.jl | 2 +- src/measurementmodels/hemodynamic_extras.jl | 8 ++++---- test/data_fitting_v2.jl | 20 ++++++++++++++++---- 5 files changed, 31 insertions(+), 13 deletions(-) diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index 048e2dc2..efb9ec3d 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -48,7 +48,7 @@ end function compileparameterlist(;kwargs...) paramlist = [] for (kw, v) in kwargs - if v isa Float64 + if v isa Float64 # note that Num is also subtype of Real. If we want to be more inclusive we need to create a union of types. paramlist = vcat(paramlist, @parameters $kw = v) else paramlist = vcat(paramlist, v) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index d4fadd38..d46ba251 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -124,11 +124,17 @@ function (bc::BloxConnector)( sys_in = get_namespaced_sys(bloxin) w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") - w = only(@parameters $(w_name)=weight) + @show "3", weight, w_name + # Main.foo[] = weight, w_name + if weight isa Num # note that Num is also subtype of Real. If we want to be more inclusive we need to create a union of types. + w = weight + else + w = only(@parameters $(w_name)=weight) + end push!(bc.weights, w) x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) eq = sys_in.nmm₊jcn ~ x*w - + accumulate_equation!(bc, eq) end @@ -188,7 +194,7 @@ function (bc::BloxConnector)( weight=1, delay=0, density=0.1 -) +) sys_out = get_namespaced_sys(stim) sys_in = get_namespaced_sys(neuron) diff --git a/src/datafitting/spectralDCM.jl b/src/datafitting/spectralDCM.jl index 02ae6efe..ea23d7f5 100644 --- a/src/datafitting/spectralDCM.jl +++ b/src/datafitting/spectralDCM.jl @@ -509,7 +509,7 @@ function spectralVI(data, neuraldynmodel, observationmodel, initcond, csdsetup, :Q => Q # decomposition of model parameter covariance ) ); - + Main.foo[] = idx_A = findall(occursin.("A[", string.(calculate_jacobian(neuraldynmodel)))) ### Compute the variational Bayes with Laplace approximation ### diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index e3b68f21..c2879d0b 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -113,15 +113,15 @@ struct CompoundHemo <: CompoundNOBlox namespace function CompoundHemo(massChoice; name, lnκ=0.0, lnτ=0.0) #ONLY WORKS WITH DEFAULT PARAMETERS FOR NOW p = progress_scope(@parameters lnκ lnτ) - @named hemo = BalloonModel(;lnκ=p[1], lnτ=p[2]) - @named nmm = massChoice() + @named hemo = BalloonModel(;lnκ=p[1], lnτ=p[2], namespace=name) + @named nmm = massChoice(;namespace=name) @variables jcn(t) g = MetaDiGraph() add_blox!(g, nmm) add_blox!(g, hemo) add_edge!(g, 1, 2, Dict(:weight => 1.0)) linhemo = system_from_graph(g; name=name) - new(p, states(linhemo)[1], states(linhemo)[2], linhemo, nothing) + new(p, states(linhemo)[1], states(linhemo)[2], linhemo, name) end end @@ -132,7 +132,7 @@ struct LinHemoCombo <: CompoundNOBlox odesystem namespace function LinHemoCombo(;name, lnκ=0.0, lnτ=0.0) - p = progress_scope(@parameters lnκ=lnκ lnτ=lnτ) + p = progress_scope(lnκ, lnτ) lnκ, lnτ = p # assign the modified parameters H = [0.64, 0.32, 2.00, 0.32, 0.4] diff --git a/test/data_fitting_v2.jl b/test/data_fitting_v2.jl index d6c3f48f..fee8b621 100644 --- a/test/data_fitting_v2.jl +++ b/test/data_fitting_v2.jl @@ -111,6 +111,11 @@ for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) end end +g, bc, blox_syss = foo[] + +ODESystem(bc.eqs, t, [], vcat(bc.weights, bc.delays); name=:name) + +foo = Ref{Any}() # compose model @named neuronmodel2 = system_from_graph(g2) neuronmodel2 = structural_simplify(neuronmodel2) @@ -129,10 +134,17 @@ end modelparam = OrderedDict() for par in parameters(neuronmodel2) - # while Symbolics.getdefaultval(par) isa Num - # par = Symbolics.getdefaultval(par) - # end - modelparam[par] = Symbolics.getdefaultval(par) + if Symbolics.getdefaultval(par) isa Num + ex = Symbolics.getdefaultval(par) + @show ex + p = only(Symbolics.get_variables(ex)) + @show p + # Symbolics.value(Symbolics.substitute(Symbolics.getdefaultval(parameters(neuronmodel2)[1]), Dict(p => Symbolics.getdefaultval(p)))) + par = Symbolics.substitute(ex, Dict(p => Symbolics.getdefaultval(p))) + modelparam[p] = Symbolics.value(par) + else + modelparam[par] = Symbolics.getdefaultval(par) + end end # Noise parameter mean modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 From 8162d2611dd7fd6a6e92b68ae62b0a88c3e54074 Mon Sep 17 00:00:00 2001 From: agchesebro <76024790+agchesebro@users.noreply.github.com> Date: Thu, 9 Nov 2023 11:32:01 -0500 Subject: [PATCH 13/28] Indexing and finding hemodynamic observer states --- src/Neuroblox.jl | 4 ++-- src/blox/blox_utilities.jl | 12 ++++++++++ src/blox/neural_mass.jl | 2 +- src/measurementmodels/hemodynamic_extras.jl | 4 ++-- test/data_fitting_extra_examples.jl | 25 +++++++++++++++++++++ 5 files changed, 42 insertions(+), 5 deletions(-) create mode 100644 test/data_fitting_extra_examples.jl diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 95d387fc..72f267a0 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -27,7 +27,7 @@ using Distributions using ModelingToolkit: get_namespace, get_systems, renamespace, namespace_equation, namespace_variables, namespace_parameters, namespace_expr, AbstractODESystem -import ModelingToolkit: inputs, nameof +import ModelingToolkit: inputs, nameof, outputs, getdescription using Symbolics: @register_symbolic using IfElse @@ -181,6 +181,6 @@ export simulate, random_initials export system_from_graph, graph_delays export create_adjacency_edges! export get_namespaced_sys, namespace_expr, nameof -export HRFFourHalfCosine, HRFDoubleGamma, BalloonModel, CompoundHemo, AlternativeBalloonModel, LinHemoCombo, JRHemo, spectralVI2 +export HRFFourHalfCosine, HRFDoubleGamma, BalloonModel, CompoundHemo, AlternativeBalloonModel, LinHemoCombo, JRHemo, spectralVI2, get_hemodynamic_observers export input_equations, BloxConnector, accumulate_equation!, get_sys #remove this line end \ No newline at end of file diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index c4cd42f7..adf68a9d 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -227,3 +227,15 @@ function count_spikes(x::AbstractVector{T}; minprom=zero(T), maxprom=nothing, mi return length(spikes) end + +function get_hemodynamic_observers(sys_from_graph) + obs_idx = [] + obs_states = [] + for (i, s) in enumerate(outputs(sys_from_graph)) + if isequal(getdescription(s), "hemodynamic_observer") + push!(obs_idx, i) + push!(obs_states, s) + end + end + return (obs_idx, obs_states) +end \ No newline at end of file diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 3f8d9b82..580a6d91 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -304,7 +304,7 @@ struct LinearNeuralMass <: NeuralMassBlox sts = @variables x(t)=0.0 [output=true] jcn(t)=0.0 [input=true] eqs = [D(x) ~ jcn] sys = System(eqs, name=name) - new(sts[1], sts[2], sys, namespace) + new(sts[1], sts[2], sys, nothing) end end diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index 38c39881..52fc3cc8 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -1,4 +1,4 @@ -using Random, SpecialFunctions, Plots +using Random, SpecialFunctions @parameters t D = Differential(t) @@ -90,7 +90,7 @@ struct BalloonModel <: ObserverBlox #p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters lnκ, lnτ = p # assign the modified parameters - sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 lnq(t)=1.0 jcn(t)=0.0 [input=true] + sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 [output=true, description="hemodynamic_observer"] lnq(t)=1.0 [output=true, description="hemodynamic_observer"] jcn(t)=0.0 [input=true] eqs = [ D(s) ~ jcn/1e2 - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1), diff --git a/test/data_fitting_extra_examples.jl b/test/data_fitting_extra_examples.jl new file mode 100644 index 00000000..44842aed --- /dev/null +++ b/test/data_fitting_extra_examples.jl @@ -0,0 +1,25 @@ +using Neuroblox, MetaGraphs +import ModelingToolkit: outputs + +@named ob1 = BalloonModel() +@named ob2 = BalloonModel() +@named nmm1 = LinearNeuralMass() +@named nmm2 = LinearNeuralMass() + +blox = [nmm1, ob1, nmm2, ob2] + +g = MetaDiGraph() +add_blox!.(Ref(g), blox) + +add_edge!(g, 1, 2, Dict(:weight => 1.0)) # Connect hemodynamic observer #1 +add_edge!(g, 3, 4, Dict(:weight => 1.0)) # Connect hemodynamic observer #2 +add_edge!(g, 1, 3, Dict(:weight => 1.0)) # Connect neural mass +add_edge!(g, 3, 1, Dict(:weight => 1.0)) # Bidirectional connection because that's the assumption in spectral spectralDCM + +@named final_system = system_from_graph(g) + +obs_idx, obs_state_names = get_hemodynamic_observers(final_system) + +# To use obs_idx +all_outputs = outputs(final_system) +obs_outputs = all_outputs[obs_idx] #equivalent to obs_state_names \ No newline at end of file From 9f8be598752ffd2bee71c0a30a0e260935ca2d32 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Mon, 27 Nov 2023 15:24:52 +0100 Subject: [PATCH 14/28] changed connections such that new connection parameters are only created if a float is passed. If a graph weight is a symbolic expression use that directly instead. --- src/blox/connections.jl | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index b7b0d872..7399058f 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -30,7 +30,11 @@ function generate_weight_param(blox_out, blox_in; kwargs...) weight = get_weight(kwargs, name_out, name_in) w_name = Symbol("w_$(name_out)_$(name_in)") - w = only(@parameters $(w_name)=weight) + if typeof(weight) == Num # Symbol + w = weight + else + w = only(@parameters $(w_name)=weight) + end return w end @@ -59,12 +63,13 @@ end function params(bc::BloxConnector) weights = [] for w in bc.weights - if Symbolics.getdefaultval(w) isa Num - p = Symbolics.get_variables(Symbolics.getdefaultval(w)) - append!(weights, p) - else - append!(weights, w) - end + append!(weights, Symbolics.get_variables(w)) + # if Symbolics.getdefaultval(w) isa Num + # p = Symbolics.get_variables(Symbolics.getdefaultval(w)) + # append!(weights, p) + # else + # append!(weights, w) + # end end return vcat(reduce(vcat, weights), bc.delays) # return vcat(bc.weights, bc.delays) @@ -126,7 +131,6 @@ function (bc::BloxConnector)( neurons_in = get_inh_neurons(cb_in) bc(asc_out, neurons_in[end]; kwargs...) - end @@ -181,7 +185,11 @@ function (bc::BloxConnector)( if typeof(bloxout.output) == Num w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") - w = only(@parameters $(w_name)=weight) + if typeof(weight) == Num # Symbol + w = weight + else + w = only(@parameters $(w_name)=weight) + end push!(bc.weights, w) x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) eq = sys_in.jcn ~ x*w @@ -217,7 +225,7 @@ function (bc::BloxConnector)( sys_in = get_namespaced_sys(bloxin) w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") - if typeof(weight) == Symbol + if typeof(weight) == Num # Symbol w = weight else w = only(@parameters $(w_name)=weight) From d37e923eb58499d7eac9522e76d1b34117cbddaa Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Mon, 27 Nov 2023 15:26:00 +0100 Subject: [PATCH 15/28] just removed some white spaces in cortical_blox.jl --- src/blox/cortical_blox.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/blox/cortical_blox.jl b/src/blox/cortical_blox.jl index 3b03669f..525f62cb 100644 --- a/src/blox/cortical_blox.jl +++ b/src/blox/cortical_blox.jl @@ -38,16 +38,14 @@ struct CorticalBlox <: AbstractComponent τ_inhib ) end - + n_ff_inh = HHNeuronInhibBlox( name = "ff_inh", - namespace = namespaced_name(namespace, name), - E_syn = E_syn_inhib, - G_syn = G_syn_ff_inhib, + namespace = namespaced_name(namespace, name), + E_syn = E_syn_inhib, + G_syn = G_syn_ff_inhib, τ = τ_inhib - ) - - + ) g = MetaDiGraph() add_blox!.(Ref(g), vcat(wtas, n_ff_inh)) From 57e0270134096acf3df877fefba16479d9e4f712 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Mon, 27 Nov 2023 15:27:29 +0100 Subject: [PATCH 16/28] add the namespace and more importantly use compileparameterlist instead of adding an @parameters in progress_scope, otherwise scoping is not working properly --- src/measurementmodels/hemodynamic_extras.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index 893f22f1..f4b9b235 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -76,7 +76,7 @@ struct BalloonModel <: ObserverBlox jcn odesystem namespace - function BalloonModel(;name, lnκ=0.0, lnτ=0.0) + function BalloonModel(;name, namespace=nothing, lnκ=0.0, lnτ=0.0) #= hemodynamic parameters H(1) - signal decay d(ds/dt)/ds) H(2) - autoregulation d(ds/dt)/df) @@ -86,8 +86,8 @@ struct BalloonModel <: ObserverBlox =# H = [0.64, 0.32, 2.00, 0.32, 0.4] - p = progress_scope(@parameters lnκ=lnκ lnτ=lnτ) # progress scope if needed - #p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters + p = progress_scope(lnκ, lnτ) # progress scope if needed + p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters lnκ, lnτ = p # assign the modified parameters sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 [output=true, description="hemodynamic_observer"] lnq(t)=1.0 [output=true, description="hemodynamic_observer"] jcn(t)=0.0 [input=true] @@ -99,7 +99,7 @@ struct BalloonModel <: ObserverBlox D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) ] sys = System(eqs, name=name) - new(p, Num(0), sts[5], sys, nothing) + new(p, Num(0), sts[5], sys, namespace) end end @@ -133,6 +133,7 @@ struct LinHemoCombo <: CompoundNOBlox namespace function LinHemoCombo(;name, lnκ=0.0, lnτ=0.0) p = progress_scope(lnκ, lnτ) + p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters lnκ, lnτ = p # assign the modified parameters H = [0.64, 0.32, 2.00, 0.32, 0.4] From 5b614cb0acd1493976ca4dae2f365cfe8a72ef4c Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Mon, 27 Nov 2023 15:28:19 +0100 Subject: [PATCH 17/28] add some namespace parameters --- src/blox/neural_mass.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index 580a6d91..becad3d2 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -300,11 +300,11 @@ struct LinearNeuralMass <: NeuralMassBlox jcn odesystem namespace - function LinearNeuralMass(;name) + function LinearNeuralMass(;name, namespace=nothing) sts = @variables x(t)=0.0 [output=true] jcn(t)=0.0 [input=true] eqs = [D(x) ~ jcn] sys = System(eqs, name=name) - new(sts[1], sts[2], sys, nothing) + new(sts[1], sts[2], sys, namespace) end end From fe874bdcbf46f5bcd9dc097ee4b905ff5206dec3 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Mon, 27 Nov 2023 15:29:47 +0100 Subject: [PATCH 18/28] add tunable to compileparameterlist, this way every parameter that is created within blox is automatically tunable. Note that weight parameters in comparison are currently not tunable by default --- src/blox/blox_utilities.jl | 46 ++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index 843394a4..1b2dc0e9 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -1,31 +1,30 @@ -function progress_scope(params; lvl=0) - para_list = [] - for p in params - pp = ModelingToolkit.unwrap(p) - if ModelingToolkit.hasdefault(pp) - d = ModelingToolkit.getdefault(pp) - if typeof(d)==SymbolicUtils.BasicSymbolic{Real} - if lvl==0 - pp = ParentScope(pp) - else - pp = DelayParentScope(pp,lvl) - end - end - end - push!(para_list,ModelingToolkit.wrap(pp)) - end - return para_list -end +# function progress_scope(params; lvl=0) +# para_list = [] +# for p in params +# pp = ModelingToolkit.unwrap(p) +# if ModelingToolkit.hasdefault(pp) +# d = ModelingToolkit.getdefault(pp) +# if typeof(d)==SymbolicUtils.BasicSymbolic{Real} +# if lvl==0 +# pp = ParentScope(pp) +# else +# pp = DelayParentScope(pp,lvl) +# end +# end +# end +# push!(para_list,ModelingToolkit.wrap(pp)) +# end +# return para_list +# end """ This function progresses the scope of parameters and leaves floating point values untouched """ function progress_scope(args...) paramlist = [] + @show args for p in args - if p isa Float64 - push!(paramlist, p) - else + if p isa Num p = ParentScope(p) # pp = ModelingToolkit.unwrap(p) # if ModelingToolkit.hasdefault(pp) @@ -36,8 +35,11 @@ function progress_scope(args...) # end # push!(para_list,ModelingToolkit.wrap(pp)) push!(paramlist, p) + else + push!(paramlist, p) end end + @show paramlist return paramlist end @@ -49,7 +51,7 @@ function compileparameterlist(;kwargs...) paramlist = [] for (kw, v) in kwargs if v isa Float64 # note that Num is also subtype of Real. If we want to be more inclusive we need to create a union of types. - paramlist = vcat(paramlist, @parameters $kw = v) + paramlist = vcat(paramlist, @parameters $kw = v [tunable=true]) else paramlist = vcat(paramlist, v) end From 6cd83075b944e0f94cd386c5aaeaee02a1607709 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Mon, 27 Nov 2023 15:57:26 +0100 Subject: [PATCH 19/28] augmented get_hemodynamic_observers such that indices and states are now grouped by region --- src/blox/blox_utilities.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index 1b2dc0e9..d609b33e 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -230,13 +230,14 @@ function count_spikes(x::AbstractVector{T}; minprom=zero(T), maxprom=nothing, mi return length(spikes) end -function get_hemodynamic_observers(sys_from_graph) - obs_idx = [] - obs_states = [] +function get_hemodynamic_observers(sys_from_graph, nr) + obs_idx = Dict([k => [] for k in 1:nr]) + obs_states = Dict([k => [] for k in 1:nr]) for (i, s) in enumerate(outputs(sys_from_graph)) if isequal(getdescription(s), "hemodynamic_observer") - push!(obs_idx, i) - push!(obs_states, s) + regionidx = parse(Int64, split(string(s), "₊")[1][end]) + push!(obs_idx[regionidx], i) + push!(obs_states[regionidx], s) end end return (obs_idx, obs_states) From 7c949ff501e66241438d78ee8a3ea98a8bf35949 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Mon, 27 Nov 2023 22:33:37 +0100 Subject: [PATCH 20/28] added function that completes the list of tunable parameters with those that are non-tunable which are also needed when computing the jacobian --- src/blox/blox_utilities.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index d609b33e..3ef40632 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -241,4 +241,18 @@ function get_hemodynamic_observers(sys_from_graph, nr) end end return (obs_idx, obs_states) +end + +function addnontunableparams(param, model) + newparam = [] + k = 0 + for p in parameters(model) + if istunable(p) + k += 1 + push!(newparam, param[k]) + else + push!(newparam, Symbolics.getdefaultval(p)) + end + end + return newparam end \ No newline at end of file From 41aca52c61b7d74a7ab2507cb9633eca5d9da57d Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Tue, 28 Nov 2023 01:35:14 +0100 Subject: [PATCH 21/28] had to switch ordering of variables in bold function to ensure a match with the ordering in the balloon model, otherwise we run into issues when computing the gradient --- src/measurementmodels/fmri.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/measurementmodels/fmri.jl b/src/measurementmodels/fmri.jl index 3125f1a4..54cdaedf 100644 --- a/src/measurementmodels/fmri.jl +++ b/src/measurementmodels/fmri.jl @@ -39,7 +39,7 @@ mutable struct Hemodynamics <: AbstractComponent params = progress_scope(lnκ, lnτ) # progress scope if needed params = compileparameterlist(lnκ=params[1], lnτ=params[2]) # finally compile all parameters lnκ, lnτ = params # assign the modified parameters - + states = @variables s(t) lnf(t) lnν(t) lnq(t) jcn(t) eqs = [ @@ -106,7 +106,7 @@ function boldsignal(;name, lnϵ=0.0) k1 = 4.3*nu0*E0*TE params = @parameters lnϵ=lnϵ - vars = @variables bold(t) lnq(t) lnν(t) + 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! 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ν)) From d9decb7d8818f5cad48cfed30a2326ea0a01f5aa Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Tue, 28 Nov 2023 01:36:16 +0100 Subject: [PATCH 22/28] removed the division by 100 of the neuronal signal in the balloon model --- src/measurementmodels/hemodynamic_extras.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl index f4b9b235..2d62f613 100644 --- a/src/measurementmodels/hemodynamic_extras.jl +++ b/src/measurementmodels/hemodynamic_extras.jl @@ -93,7 +93,7 @@ struct BalloonModel <: ObserverBlox sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 [output=true, description="hemodynamic_observer"] lnq(t)=1.0 [output=true, description="hemodynamic_observer"] jcn(t)=0.0 [input=true] eqs = [ - D(s) ~ jcn/1e2 - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1), + D(s) ~ jcn - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1), D(lnf) ~ s / exp(lnf), D(lnν) ~ (exp(lnf) - exp(lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(lnν)), D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) From 813b7926fa3e03c77b645b291f1f16354eb57f93 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Tue, 28 Nov 2023 01:37:59 +0100 Subject: [PATCH 23/28] WIP: integration almost finished. However, fitting values still wrong. Gradient of bold signal is correct, jacobian of hemodynamic response is erroneous --- src/Neuroblox.jl | 2 +- src/blox/blox_utilities.jl | 1 + src/datafitting/spectralDCM.jl | 36 ++++++++++------- test/data_fitting_v2.jl | 70 ++++++++++++++++++++++++---------- 4 files changed, 75 insertions(+), 34 deletions(-) diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 72f267a0..4feb792f 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -181,6 +181,6 @@ export simulate, random_initials export system_from_graph, graph_delays export create_adjacency_edges! export get_namespaced_sys, namespace_expr, nameof -export HRFFourHalfCosine, HRFDoubleGamma, BalloonModel, CompoundHemo, AlternativeBalloonModel, LinHemoCombo, JRHemo, spectralVI2, get_hemodynamic_observers +export HRFFourHalfCosine, HRFDoubleGamma, BalloonModel, CompoundHemo, AlternativeBalloonModel, LinHemoCombo, JRHemo, addnontunableparams, get_hemodynamic_observers export input_equations, BloxConnector, accumulate_equation!, get_sys #remove this line end \ No newline at end of file diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index 3ef40632..ba0480e7 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -254,5 +254,6 @@ function addnontunableparams(param, model) push!(newparam, Symbolics.getdefaultval(p)) end end + append!(newparam, param[k+1:end]) return newparam end \ No newline at end of file diff --git a/src/datafitting/spectralDCM.jl b/src/datafitting/spectralDCM.jl index 9df2aa21..8ca491e2 100644 --- a/src/datafitting/spectralDCM.jl +++ b/src/datafitting/spectralDCM.jl @@ -78,6 +78,7 @@ function transferfunction_fmri(w, idx_A, derivatives, params) C = params[(6+2nd+nd^2):(5+3nd+nd^2)] C /= 16.0 # TODO: unclear why C is devided by 16 but see spm_fx_fmri.m:49 ∂f = derivatives[:∂f](params[1:(nd^2+nd+1)]) #convert(Array{Real}, substitute(derivatives[:∂f], params)) +Main.bar[] = ∂f if ∂f isa Vector ∂f = reshape(∂f, nd*5, nd*5) # TODO: generalize this to arbitrary number of states! This is specific for LinHemo! end @@ -108,6 +109,8 @@ function transferfunction_fmri(w, idx_A, derivatives, params) end end end + Main.bar2[] = S, ∂g, params + return S end @@ -290,6 +293,7 @@ end for k = 1:niter state.iter = k dfdp = ForwardDiff.jacobian(f_prep, μθ_po) * V +Main.foo[] = dfdp, μθ_po, V norm_dfdp = matlab_norm(dfdp, Inf); revert = isnan(norm_dfdp) || norm_dfdp > exp(32); @@ -464,26 +468,30 @@ Input: function spectralVI(data, neuraldynmodel, observationmodel, initcond, csdsetup, params, hyperparams) # compute cross-spectral density y = Matrix(data); - nd = ncol(data); # dimension of the data, number of regions - dt = csdsetup[:dt]; # order of MAR. Hard-coded in SPM12 with this value. We will use the same for now. - freqs = csdsetup[:freq]; # frequencies at which the CSD is evaluated - p = csdsetup[:p]; # order of MAR + nr = ncol(data); # number of regions + ns = length(initcond) # number of states in total + dt = csdsetup[:dt]; # order of MAR. Hard-coded in SPM12 with this value. We will use the same for now. + freqs = csdsetup[:freq]; # frequencies at which the CSD is evaluated + p = csdsetup[:p]; # order of MAR mar = mar_ml(y, p); # compute MAR from time series y and model order p y_csd = mar2csd(mar, freqs, dt^-1); # compute cross spectral densities from MAR parameters at specific frequencies freqs, dt^-1 is sampling rate of data - grad_full = function(p, grad, sts, nd) - tmp = zeros(typeof(p), nd, length(sts)) - for i in 1:nd - tmp[i, (i-1)*5 .+ (4:5)] = grad(vcat([1], sts[(i-1)*5 .+ (5:-1:4)]), p, t)[3:-1:2] + grad_full = function(grad, obsstates, obsidx, params, nr, ns) + tmp = zeros(typeof(params), nr, ns) + for i in 1:nr + tmp[i, obsidx[i]] = grad(vcat([1], obsstates[i]), params, t)[2:end] # [(length(obsstates[i])+1):-1:2] end return tmp end + jac_f = generate_jacobian(neuraldynmodel, expression = Val{false})[1] - grad_g = generate_jacobian(observationmodel, expression = Val{false})[1] - + grad_g = generate_jacobian(observationmodel, expression = Val{false})[1] # computes gradient since output is a scalar + statevals = [v for v in values(initcond)] - derivatives = Dict(:∂f => par -> jac_f(statevals, par, t), - :∂g => par -> grad_full(par, grad_g, statevals, nd)) + obs = get_hemodynamic_observers(neuraldynmodel, nr) + obsstates = map(obs -> [initcond[s] for s in obs], values(obs[2])) + derivatives = Dict(:∂f => par -> jac_f(statevals, addnontunableparams(par, neuraldynmodel), t), + :∂g => par -> grad_full(grad_g, obsstates, obs[1], par, nr, ns)) θΣ = diagm(vecparam(OrderedDict(params.name .=> params.variance))) # depending on the definition of the priors (note that we take it from the SPM12 code), some dimensions are set to 0 and thus are not changed. @@ -508,7 +516,9 @@ function spectralVI(data, neuraldynmodel, observationmodel, initcond, csdsetup, :Q => Q # decomposition of model parameter covariance ) ); - Main.foo[] = + # the following is needed to properly deal with C. + # Fixing this likely means departing from the SPM12 version, which should happen at some point since it is an odd part, + # see transferfunction idx_A = findall(occursin.("A[", string.(calculate_jacobian(neuraldynmodel)))) ### Compute the variational Bayes with Laplace approximation ### diff --git a/test/data_fitting_v2.jl b/test/data_fitting_v2.jl index e51557ab..a4eda84a 100644 --- a/test/data_fitting_v2.jl +++ b/test/data_fitting_v2.jl @@ -94,24 +94,36 @@ nd = ncol(data) # number of dimensions ########## assemble the model ########## -@parameters κ=0.0 # define brain-wide decay parameter for hemodynamics -g2 = MetaDiGraph() +g = MetaDiGraph() +# shouldn't the following be somehow dealt with a system_from_parts() call? But that does not allow for connections between subsystems it seems +# I would like to do so to create a region that is composed of different blox that are interconnected. +# can't just build systems from graphs first and then connect them with another graph...? +regions = Dict() +@parameters κ=0.0 [tunable = true] # define brain-wide decay parameter for hemodynamics for ii = 1:nd - region = LinHemoCombo(;name=Symbol("r$ii"), lnκ=κ) - add_blox!(g2, region) + region = LinearNeuralMass(;name=Symbol("r$(ii)₊lm")) + add_blox!(g, region) + regions[ii] = 2ii - 1 # store index of neural mass model + # add hemodynamic observer + observer = BalloonModel(;name=Symbol("r$(ii)₊bm"), lnκ=κ) + add_blox!(g, observer) + # connect observer with neuronal signal + add_edge!(g, 2ii - 1, 2ii, Dict(:weight => 1.0)) + # region = LinHemoCombo(;name=Symbol("r$ii"), lnκ=κ) end # add symbolic weights -@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) +@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) [tunable = true] for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) if idx[1] == idx[2] - add_edge!(g, idx[1], idx[2], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 + add_edge!(g, regions[idx[1]], regions[idx[2]], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 else - add_edge!(g, idx[1], idx[2], :weight, A[i]) + add_edge!(g, regions[idx[2]], regions[idx[1]], :weight, A[i]) end end + # compose model -@named neuronmodel = ODEfromGraph(g) +@named neuronmodel = system_from_graph(g) neuronmodel = structural_simplify(neuronmodel) # measurement model @@ -121,7 +133,13 @@ neuronmodel = structural_simplify(neuronmodel) all_s = states(neuronmodel) initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() rnames = [] -x = rand(3, 6) #WHY IS x 3x5? +# ns = Int(length(all_s)/nd) +# for ii = 1:nd +# for jj = 1:ns +# initcond[all_s[(ii-1)*ns+jj]] = 0.0 +# end +# end +# x = rand(3, 6) #WHY IS x 3x5? map(x->push!(rnames, split(string(x), "₊")[1]), all_s); rnames = unique(rnames); for (i, r) in enumerate(rnames) @@ -132,15 +150,16 @@ end modelparam = OrderedDict() for par in parameters(neuronmodel) - if Symbolics.getdefaultval(par) isa Num - ex = Symbolics.getdefaultval(par) - @show ex - p = only(Symbolics.get_variables(ex)) - @show p - # Symbolics.value(Symbolics.substitute(Symbolics.getdefaultval(parameters(neuronmodel2)[1]), Dict(p => Symbolics.getdefaultval(p)))) - par = Symbolics.substitute(ex, Dict(p => Symbolics.getdefaultval(p))) - modelparam[p] = Symbolics.value(par) - else + # if Symbolics.getdefaultval(par) isa Num + # ex = Symbolics.getdefaultval(par) + # @show ex + # p = only(Symbolics.get_variables(ex)) + # @show p + # # Symbolics.value(Symbolics.substitute(Symbolics.getdefaultval(parameters(neuronmodel2)[1]), Dict(p => Symbolics.getdefaultval(p)))) + # par = Symbolics.substitute(ex, Dict(p => Symbolics.getdefaultval(p))) + # modelparam[p] = Symbolics.value(par) + # else + if istunable(par) modelparam[par] = Symbolics.getdefaultval(par) end end @@ -156,6 +175,7 @@ end # define prior variances paramvariance = copy(modelparam) +# [paramvariance[k] = 0.0 for k in keys(paramvariance) if occursin("w_lm", string(k)) && occursin("bm", string(k))] paramvariance[:C] = zeros(Float64, nd); paramvariance[:lnγ] = ones(Float64, nd)./64.0; paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; @@ -176,9 +196,19 @@ params = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(m hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector ) - +foo = Ref{Any}() +bar = Ref{Any}() +bar2 = Ref{Any}() csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) -results = spectralVI2(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) +results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) + +J, iΣ, Πθ_pr = foo[] +S, ∂g, ps = bar2[] +∂f = bar[] + +obs = get_hemodynamic_observers(neuronmodel, nd) +obsstates = map(obs -> [initcond[s] for s in obs], values(obs[2])) + ### COMPARE RESULTS WITH MATLAB RESULTS ### @show results.F, vars["F"] From 7658a3b93972937d3127bf6c203ca1faaf62a2f7 Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Wed, 29 Nov 2023 15:25:19 +0100 Subject: [PATCH 24/28] fixed issue. Needed to use states rather than outputs in get_hemodynamic_observers to have the correct indices. However, there is still a very subtle deviation of the results with previous version. Free Energy is 709 instead of 704. Need to figure where that comes from. Should be the last issue! --- src/blox/blox_utilities.jl | 2 +- src/datafitting/spectralDCM.jl | 61 ++-------------------------------- test/data_fitting_v2.jl | 11 +----- 3 files changed, 4 insertions(+), 70 deletions(-) diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index ba0480e7..4f5212fa 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -233,7 +233,7 @@ end function get_hemodynamic_observers(sys_from_graph, nr) obs_idx = Dict([k => [] for k in 1:nr]) obs_states = Dict([k => [] for k in 1:nr]) - for (i, s) in enumerate(outputs(sys_from_graph)) + for (i, s) in enumerate(states(sys_from_graph)) if isequal(getdescription(s), "hemodynamic_observer") regionidx = parse(Int64, split(string(s), "₊")[1][end]) push!(obs_idx[regionidx], i) diff --git a/src/datafitting/spectralDCM.jl b/src/datafitting/spectralDCM.jl index 8ca491e2..0d68d622 100644 --- a/src/datafitting/spectralDCM.jl +++ b/src/datafitting/spectralDCM.jl @@ -78,9 +78,8 @@ function transferfunction_fmri(w, idx_A, derivatives, params) C = params[(6+2nd+nd^2):(5+3nd+nd^2)] C /= 16.0 # TODO: unclear why C is devided by 16 but see spm_fx_fmri.m:49 ∂f = derivatives[:∂f](params[1:(nd^2+nd+1)]) #convert(Array{Real}, substitute(derivatives[:∂f], params)) -Main.bar[] = ∂f if ∂f isa Vector - ∂f = reshape(∂f, nd*5, nd*5) # TODO: generalize this to arbitrary number of states! This is specific for LinHemo! + ∂f = reshape(∂f, sqrt(length(∂f)), sqrt(length(∂f))) end dfdu = zeros(eltype(C), size(∂f, 1), length(C)) @@ -109,7 +108,6 @@ Main.bar[] = ∂f end end end - Main.bar2[] = S, ∂g, params return S end @@ -293,7 +291,6 @@ end for k = 1:niter state.iter = k dfdp = ForwardDiff.jacobian(f_prep, μθ_po) * V -Main.foo[] = dfdp, μθ_po, V norm_dfdp = matlab_norm(dfdp, Inf); revert = isnan(norm_dfdp) || norm_dfdp > exp(32); @@ -479,7 +476,7 @@ function spectralVI(data, neuraldynmodel, observationmodel, initcond, csdsetup, grad_full = function(grad, obsstates, obsidx, params, nr, ns) tmp = zeros(typeof(params), nr, ns) for i in 1:nr - tmp[i, obsidx[i]] = grad(vcat([1], obsstates[i]), params, t)[2:end] # [(length(obsstates[i])+1):-1:2] + tmp[i, obsidx[i]] = grad(vcat([1], obsstates[i]), params, t)[2:end] # [(length(obsstates[i])+1):-1:2] TODO: using the reverse will improve results but is wrong. end return tmp end @@ -521,60 +518,6 @@ function spectralVI(data, neuraldynmodel, observationmodel, initcond, csdsetup, # see transferfunction idx_A = findall(occursin.("A[", string.(calculate_jacobian(neuraldynmodel)))) - ### Compute the variational Bayes with Laplace approximation ### - return variationalbayes(idx_A, y_csd, derivatives, freqs, V, p, priors, 128) -end - -function spectralVI2(data, neuraldynmodel, observationmodel, initcond, csdsetup, params, hyperparams) - # compute cross-spectral density - y = Matrix(data); - nd = ncol(data); # dimension of the data, number of regions - dt = csdsetup[:dt]; # order of MAR. Hard-coded in SPM12 with this value. We will use the same for now. - freqs = csdsetup[:freq]; # frequencies at which the CSD is evaluated - p = csdsetup[:p]; # order of MAR - mar = mar_ml(y, p); # compute MAR from time series y and model order p - y_csd = mar2csd(mar, freqs, dt^-1); # compute cross spectral densities from MAR parameters at specific frequencies freqs, dt^-1 is sampling rate of data - - grad_full = function(p, grad, sts, nd) - tmp = zeros(typeof(p), nd, length(sts)) - for i in 1:nd - tmp[i, (i-1)*6 .+ (5:6)] = grad(vcat([1], sts[(i-1)*6 .+ (6:-1:5)]), p, t)[4:-1:3] - end - return tmp - end - jac_f = generate_jacobian(neuraldynmodel, expression = Val{false})[1] - grad_g = generate_jacobian(observationmodel, expression = Val{false})[1] - - statevals = [v for v in values(initcond)] - derivatives = Dict(:∂f => par -> jac_f(statevals, par, t), - :∂g => par -> grad_full(par, grad_g, statevals, nd)) - - θΣ = diagm(vecparam(OrderedDict(params.name .=> params.variance))) - # depending on the definition of the priors (note that we take it from the SPM12 code), some dimensions are set to 0 and thus are not changed. - # Extract these dimensions and remove them from the remaining computation. I find this a bit odd and further thoughts would be necessary to understand - # to what extend this is a the most reasonable approach. - idx = findall(x -> x != 0, θΣ); - V = zeros(size(θΣ, 1), length(idx)); - order = sortperm(θΣ[idx], rev=true); - idx = idx[order]; - for i = 1:length(idx) - V[idx[i][1], i] = 1.0 - end - θΣ = V'*θΣ*V; # reduce dimension by removing columns and rows that are all 0 - - ### Collect prior means and covariances ### - Q = csd_Q(y_csd); # compute prior of Q, the precision (of the data) components. See Friston etal. 2007 Appendix A - priors = Dict(:μ => OrderedDict(params.name .=> params.mean), - :Σ => Dict( - :Πθ_pr => inv(θΣ), # prior model parameter precision - :Πλ_pr => hyperparams[:Πλ_pr], # prior metaparameter precision - :μλ_pr => hyperparams[:μλ_pr], # prior metaparameter mean - :Q => Q # decomposition of model parameter covariance - ) - ); - - idx_A = findall(occursin.("A[", string.(calculate_jacobian(neuraldynmodel)))) - ### Compute the variational Bayes with Laplace approximation ### return variationalbayes(idx_A, y_csd, derivatives, freqs, V, p, priors, 128) end \ No newline at end of file diff --git a/test/data_fitting_v2.jl b/test/data_fitting_v2.jl index a4eda84a..9c4411a9 100644 --- a/test/data_fitting_v2.jl +++ b/test/data_fitting_v2.jl @@ -196,19 +196,10 @@ params = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(m hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector ) -foo = Ref{Any}() -bar = Ref{Any}() -bar2 = Ref{Any}() + csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) -J, iΣ, Πθ_pr = foo[] -S, ∂g, ps = bar2[] -∂f = bar[] - -obs = get_hemodynamic_observers(neuronmodel, nd) -obsstates = map(obs -> [initcond[s] for s in obs], values(obs[2])) - ### COMPARE RESULTS WITH MATLAB RESULTS ### @show results.F, vars["F"] From ea55407c1eb0b1bf0737467ad1b6ffe6078c3dbe Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Fri, 1 Dec 2023 03:21:57 +0100 Subject: [PATCH 25/28] fixed issue of reduce over empty set when there are no weights as for instance in Striatum and other subcortical blox. --- src/blox/blox_utilities.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/blox/blox_utilities.jl b/src/blox/blox_utilities.jl index b2cb5cb6..5c368571 100644 --- a/src/blox/blox_utilities.jl +++ b/src/blox/blox_utilities.jl @@ -22,7 +22,6 @@ This function progresses the scope of parameters and leaves floating point value """ function progress_scope(args...) paramlist = [] - @show args for p in args if p isa Num p = ParentScope(p) @@ -39,7 +38,6 @@ function progress_scope(args...) push!(paramlist, p) end end - @show paramlist return paramlist end @@ -50,7 +48,7 @@ end function compileparameterlist(;kwargs...) paramlist = [] for (kw, v) in kwargs - if v isa Float64 # note that Num is also subtype of Real. If we want to be more inclusive we need to create a union of types. + if v isa Union{Float64, Int} # note that Num is also subtype of Real. Thus union of types seems to be the solution. paramlist = vcat(paramlist, @parameters $kw = v [tunable=true]) else paramlist = vcat(paramlist, v) From cb71273ae3f256cc7e885601de41cd3a756a95de Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Fri, 1 Dec 2023 03:29:23 +0100 Subject: [PATCH 26/28] fixed issue of reduce over empty set when there are no weights as for instance in Striatum and other subcortical blox. And in compileparameterlist I added an Int alongside the Float64 since some blox pass integers --- src/blox/connections.jl | 61 ++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 32 deletions(-) diff --git a/src/blox/connections.jl b/src/blox/connections.jl index 025a5e7e..47eed8df 100644 --- a/src/blox/connections.jl +++ b/src/blox/connections.jl @@ -85,15 +85,12 @@ function params(bc::BloxConnector) weights = [] for w in bc.weights append!(weights, Symbolics.get_variables(w)) - # if Symbolics.getdefaultval(w) isa Num - # p = Symbolics.get_variables(Symbolics.getdefaultval(w)) - # append!(weights, p) - # else - # append!(weights, w) - # end end - return vcat(reduce(vcat, weights), bc.delays) - # return vcat(bc.weights, bc.delays) + if isempty(weights) + return vcat(weights, bc.delays) + else + return vcat(reduce(vcat, weights), bc.delays) + end end function (bc::BloxConnector)( @@ -229,31 +226,31 @@ function (bc::BloxConnector)( accumulate_equation!(bc, eq) end -# Ok yes this is a bad dispatch but the whole compound blocks implementation is hacky and needs fixing @@ -# Opening an issue to loop back to this during clean up week -function (bc::BloxConnector)( - bloxout::CompoundNOBlox, - bloxin::CompoundNOBlox; - weight=1, - delay=0, - density=0.1 -) - - sys_out = get_namespaced_sys(bloxout) - sys_in = get_namespaced_sys(bloxin) - - w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") - if typeof(weight) == Num # Symbol - w = weight - else - w = only(@parameters $(w_name)=weight) - end - push!(bc.weights, w) - x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) - eq = sys_in.nmm₊jcn ~ x*w +# # Ok yes this is a bad dispatch but the whole compound blocks implementation is hacky and needs fixing @@ +# # Opening an issue to loop back to this during clean up week +# function (bc::BloxConnector)( +# bloxout::CompoundNOBlox, +# bloxin::CompoundNOBlox; +# weight=1, +# delay=0, +# density=0.1 +# ) + +# sys_out = get_namespaced_sys(bloxout) +# sys_in = get_namespaced_sys(bloxin) + +# w_name = Symbol("w_$(nameof(sys_out))_$(nameof(sys_in))") +# if typeof(weight) == Num # Symbol +# w = weight +# else +# w = only(@parameters $(w_name)=weight) +# end +# push!(bc.weights, w) +# x = namespace_expr(bloxout.output, sys_out, nameof(sys_out)) +# eq = sys_in.nmm₊jcn ~ x*w - accumulate_equation!(bc, eq) -end +# accumulate_equation!(bc, eq) +# end function (bc::BloxConnector)( wta_out::WinnerTakeAllBlox, From c3b23bbd3aa47a969c024bedb72f4c54bbebb94f Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Fri, 1 Dec 2023 03:32:22 +0100 Subject: [PATCH 27/28] it is important to note that the use of progress_scope with the @parameters is wrong, the correct use is to call compileparameterlist which will then produce symbolics if numbers are passed or directly use the symbol if that is passed instead. --- src/blox/neural_mass.jl | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/blox/neural_mass.jl b/src/blox/neural_mass.jl index becad3d2..243f6e2f 100644 --- a/src/blox/neural_mass.jl +++ b/src/blox/neural_mass.jl @@ -20,7 +20,8 @@ mutable struct HarmonicOscillatorBlox <: NeuralMassBlox initial::Dict{Num, Tuple{Float64, Float64}} odesystem::ODESystem function HarmonicOscillatorBlox(;name, ω=25*(2*pi), ζ=1.0, k=625*(2*pi), h=35.0) - params = progress_scope(@parameters ω=ω ζ=ζ k=k h=h) + params = progress_scope(ω, ζ, k, h) + params = compileparameterlist(ω=params[1], ζ=params[2], k=params[3], h=params[4]) sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0 ω, ζ, k, h = params eqs = [D(x) ~ y-(2*ω*ζ*x)+ k*(2/π)*(atan((jcn)/h)) @@ -45,7 +46,8 @@ mutable struct JansenRitCBlox <: NeuralMassBlox initial::Dict{Num, Tuple{Float64, Float64}} odesystem::ODESystem function JansenRitCBlox(;name, τ=0.001, H=20.0, λ=5.0, r=0.15) - params = progress_scope(@parameters τ=τ H=H λ=λ r=r) + params = progress_scope(τ, H, λ, r) + params = compileparameterlist(τ=params[1], H=params[2], λ=params[3], r=params[4]) sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0 τ, H, λ, r = params eqs = [D(x) ~ y - ((2/τ)*x), @@ -66,7 +68,8 @@ mutable struct JansenRitSCBlox <: NeuralMassBlox initial::Dict{Num, Tuple{Float64, Float64}} odesystem::ODESystem function JansenRitSCBlox(;name, τ=0.014, H=20.0, λ=400.0, r=0.1) - params = progress_scope(@parameters τ=τ H=H λ=λ r=r) + params = progress_scope(τ, H, λ, r) + params = compileparameterlist(τ=params[1], H=params[2], λ=params[3], r=params[4]) sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0 τ, H, λ, r = params eqs = [D(x) ~ y - ((2/τ)*x), @@ -319,7 +322,8 @@ struct HarmonicOscillator <: NeuralMassBlox odesystem namespace function HarmonicOscillator(;name, namespace=nothing, ω=25*(2*pi)*0.001, ζ=1.0, k=625*(2*pi), h=35.0) - p = progress_scope(@parameters ω=ω ζ=ζ k=k h=h) + p = progress_scope(ω, ζ, k, h) + p = compileparameterlist(ω=p[1], ζ=p[2], k=p[3], h=p[4]) sts = @variables x(t)=1.0 [output=true] y(t)=1.0 jcn(t)=0.0 [input=true] ω, ζ, k, h = p eqs = [D(x) ~ y-(2*ω*ζ*x)+ k*(2/π)*(atan((jcn)/h)) @@ -353,7 +357,8 @@ struct JansenRit <: NeuralMassBlox λ = isnothing(λ) ? (cortical ? 5.0 : 400.0) : λ r = isnothing(r) ? (cortical ? 0.15 : 0.1) : r - p = progress_scope(@parameters τ=τ H=H λ=λ r=r) + p = progress_scope(τ, H, λ, r) + p = compileparameterlist(τ=p[1], H=p[2], λ=p[3], r=p[4]) τ, H, λ, r = p sts = @variables x(..)=1.0 [output=true] y(t)=1.0 jcn(t)=0.0 [input=true] eqs = [D(x(t)) ~ y - ((2/τ)*x(t)), @@ -389,8 +394,8 @@ struct WilsonCowan <: NeuralMassBlox θ_I=3.5, η=1.0 ) - p = progress_scope(@parameters τ_E=τ_E τ_I=τ_I a_E=a_E a_I=a_I c_EE=c_EE c_IE=c_IE c_EI=c_EI c_II=c_II θ_E=θ_E θ_I=θ_I η=η) - + p = progress_scope(τ_E, τ_I, a_E, a_I, c_EE, c_IE, c_EI, c_II, θ_E, θ_I, η) + p = compileparameterlist(τ_E=p[1], τ_I=p[2], a_E=p[3], a_I=p[4], c_EE=p[5], c_IE=p[6], c_EI=p[7], c_II=p[8], θ_E=p[9], θ_I=p[10], η=p[11]) τ_E, τ_I, a_E, a_I, c_EE, c_IE, c_EI, c_II, θ_E, θ_I, η = p sts = @variables E(t)=1.0 [output=true] I(t)=1.0 jcn(t)=0.0 [input=true] #P(t)=0.0 eqs = [D(E) ~ -E/τ_E + 1/(1 + exp(-a_E*(c_EE*E - c_IE*I - θ_E + η*(jcn)))), #old form: D(E) ~ -E/τ_E + 1/(1 + exp(-a_E*(c_EE*E - c_IE*I - θ_E + P + η*(jcn)))), @@ -443,7 +448,8 @@ struct LarterBreakspear <: NeuralMassBlox r_NMDA=0.25, C=0.35 ) - p = progress_scope(@parameters C=C δ_VZ=δ_VZ T_Ca=T_Ca δ_Ca=δ_Ca g_Ca=g_Ca V_Ca=V_Ca T_K=T_K δ_K=δ_K g_K=g_K V_K=V_K T_Na=T_Na δ_Na=δ_Na g_Na=g_Na V_Na=V_Na V_L=V_L g_L=g_L V_T=V_T Z_T=Z_T Q_Vmax=Q_Vmax Q_Zmax=Q_Zmax IS=IS a_ee=a_ee a_ei=a_ei a_ie=a_ie a_ne=a_ne a_ni=a_ni b=b τ_K=τ_K ϕ=ϕ r_NMDA=r_NMDA) + p = progress_scope(C, δ_VZ, T_Ca, δ_Ca, g_Ca, V_Ca, T_K, δ_K, g_K, V_K, T_Na, δ_Na, g_Na, V_Na, V_L, g_L, V_T, Z_T, Q_Vmax, Q_Zmax, IS, a_ee, a_ei, a_ie, a_ne, a_ni, b, τ_K, ϕ,r_NMDA) + p = compileparameterlist(C=p[1], δ_VZ=p[2], T_Ca=p[3], δ_Ca=p[4], g_Ca=p[5], V_Ca=p[6], T_K=p[7], δ_K=p[8], g_K=p[9], V_K=p[10], T_Na=p[11], δ_Na=p[12], g_Na=p[13],V_Na=p[14], V_L=p[15], g_L=p[16], V_T=p[17], Z_T=p[18], Q_Vmax=p[19], Q_Zmax=p[20], IS=p[21], a_ee=p[22], a_ei=p[23], a_ie=p[24], a_ne=p[25], a_ni=p[26], b=p[27], τ_K=p[28], ϕ=p[29], r_NMDA=p[30]) C, δ_VZ, T_Ca, δ_Ca, g_Ca, V_Ca, T_K, δ_K, g_K, V_K, T_Na, δ_Na, g_Na,V_Na, V_L, g_L, V_T, Z_T, Q_Vmax, Q_Zmax, IS, a_ee, a_ei, a_ie, a_ne, a_ni, b, τ_K, ϕ, r_NMDA = p sts = @variables V(t)=0.5 Z(t)=0.5 W(t)=0.5 jcn(t)=0.0 [input=true] Q_V(t) [output=true] Q_Z(t) m_Ca(t) m_Na(t) m_K(t) From 1f06d6f3547893b626c335a9dc1b43e849deb68d Mon Sep 17 00:00:00 2001 From: david-hofmann Date: Fri, 1 Dec 2023 03:32:59 +0100 Subject: [PATCH 28/28] some cleanup to wrap up --- src/Neuroblox.jl | 6 +- src/measurementmodels/fmri.jl | 47 ++--- src/measurementmodels/hemodynamic_extras.jl | 193 ------------------ test/data_fitting_extra_examples.jl | 25 --- test/data_fitting_v2.jl | 207 -------------------- test/datafitting.jl | 30 +-- 6 files changed, 38 insertions(+), 470 deletions(-) delete mode 100644 src/measurementmodels/hemodynamic_extras.jl delete mode 100644 test/data_fitting_extra_examples.jl delete mode 100644 test/data_fitting_v2.jl diff --git a/src/Neuroblox.jl b/src/Neuroblox.jl index 41aa6860..c79062d0 100644 --- a/src/Neuroblox.jl +++ b/src/Neuroblox.jl @@ -52,6 +52,7 @@ abstract type AbstractNeuronBlox <: AbstractBlox end abstract type NeuralMassBlox <: AbstractBlox end abstract type CompositeBlox <: AbstractBlox end abstract type StimulusBlox <: AbstractBlox end +abstract type ObserverBlox <: AbstractBlox end # we define these in neural_mass.jl # abstract type HarmonicOscillatorBlox <: NeuralMassBlox end @@ -101,7 +102,6 @@ include("gui/GUI.jl") include("blox/connections.jl") include("blox/blox_utilities.jl") include("Neurographs.jl") -include("measurementmodels/hemodynamic_extras.jl") function simulate(sys::ODESystem, u0, timespan, p, solver = AutoVern7(Rodas4()); kwargs...) prob = ODEProblem(sys, u0, timespan, p) @@ -173,12 +173,12 @@ export LinearConnections, SynapticConnections, ODEfromGraph, ODEfromGraphNeuron, export add_blox! export powerspectrum, complexwavelet, bandpassfilter, hilberttransform, phaseangle, mar2csd, csd2mar, mar_ml export learningrate, ControlError -export Hemodynamics, LinHemo, boldsignal +export boldsignal, BalloonModel export vecparam, unvecparam, csd_Q, spectralVI export simulate, random_initials export system_from_graph, graph_delays export create_adjacency_edges! export get_namespaced_sys, namespace_expr, nameof -export HRFFourHalfCosine, HRFDoubleGamma, BalloonModel, CompoundHemo, AlternativeBalloonModel, LinHemoCombo, JRHemo, addnontunableparams, get_hemodynamic_observers +export addnontunableparams, get_hemodynamic_observers export input_equations, BloxConnector, accumulate_equation!, get_sys #remove this line end \ No newline at end of file diff --git a/src/measurementmodels/fmri.jl b/src/measurementmodels/fmri.jl index 54cdaedf..50dd2328 100644 --- a/src/measurementmodels/fmri.jl +++ b/src/measurementmodels/fmri.jl @@ -3,8 +3,8 @@ fmri.jl Models of fMRI signals. -hemodynamics : computes hemodynamic responses and its Jacobian -boldsignal : computes BOLD signal and gradient +BalloonModel : computes hemodynamic responses +boldsignal : computes BOLD signal """ """ @@ -23,10 +23,13 @@ lnτ : logarithmic prefactor to transit time H[3], set to 0 for standard parame ### Return variables ### returns an ODESystem of the biophysical model for the hemodynamics """ -mutable struct Hemodynamics <: AbstractComponent - connector::Num - odesystem::ODESystem - function Hemodynamics(;name, lnκ=0.0, lnτ=0.0) +struct BalloonModel <: ObserverBlox + params + output + jcn + odesystem + namespace + function BalloonModel(;name, namespace=nothing, lnκ=0.0, lnτ=0.0) #= hemodynamic parameters H(1) - signal decay d(ds/dt)/ds) H(2) - autoregulation d(ds/dt)/df) @@ -36,11 +39,11 @@ mutable struct Hemodynamics <: AbstractComponent =# H = [0.64, 0.32, 2.00, 0.32, 0.4] - params = progress_scope(lnκ, lnτ) # progress scope if needed - params = compileparameterlist(lnκ=params[1], lnτ=params[2]) # finally compile all parameters - lnκ, lnτ = params # assign the modified parameters - - states = @variables s(t) lnf(t) lnν(t) lnq(t) jcn(t) + p = progress_scope(lnκ, lnτ) # progress scope if needed + p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters + lnκ, lnτ = p # assign the modified parameters + + sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 [output=true, description="hemodynamic_observer"] lnq(t)=1.0 [output=true, description="hemodynamic_observer"] jcn(t)=0.0 [input=true] eqs = [ D(s) ~ jcn - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1), @@ -48,30 +51,12 @@ mutable struct Hemodynamics <: AbstractComponent D(lnν) ~ (exp(lnf) - exp(lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(lnν)), D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) ] - odesys = ODESystem(eqs, t, states, params; name=name) - new(Num(0), odesys) + sys = System(eqs, name=name) + new(p, Num(0), sts[5], sys, namespace) end end -mutable struct LinHemo <: AbstractComponent - connector::Num - bloxinput::Num - odesystem::ODESystem - function LinHemo(;name, lnκ=0.0, lnτ=0.0) - params = progress_scope(lnκ, lnτ) - @named hemo = Hemodynamics(;lnκ=params[1], lnτ=params[2]) - @named nmm = LinearNeuralMassBlox() - @variables jcn(t) - - g = MetaDiGraph() - add_vertex!(g, Dict(:blox => nmm, :jcn => jcn)) - add_vertex!(g, :blox, hemo) - add_edge!(g, 1, 2, :weight, 1.0) - linhemo = ODEfromGraph(g; name=name) - new(linhemo.nmm₊x, linhemo.jcn, linhemo) - end -end """ BOLD signal model as described in: diff --git a/src/measurementmodels/hemodynamic_extras.jl b/src/measurementmodels/hemodynamic_extras.jl deleted file mode 100644 index 2d62f613..00000000 --- a/src/measurementmodels/hemodynamic_extras.jl +++ /dev/null @@ -1,193 +0,0 @@ -using Random, SpecialFunctions - -@parameters t -D = Differential(t) - -""" -Woolrich, Behrens and Smith. 2003. - -Specify widths (h₁, h₂, h₃, h₄) in seconds. -Amplitudes (f₁, f₂) relative to HRF height. -Output is HRF kernel in ms. - -""" - -function HRFFourHalfCosine(; h₁=nothing, h₂=nothing, h₃=nothing, h₄=nothing, f₁=nothing, f₂=nothing) - - h₁ = isnothing(h₁) ? 2*rand() : h₁ - h₂ = isnothing(h₂) ? 4*rand()+2 : h₂ - h₃ = isnothing(h₃) ? 4*rand()+2 : h₃ - h₄ = isnothing(h₄) ? 6*rand()+2 : h₄ - f₁ = isnothing(f₁) ? 0 : f₁ - f₂ = isnothing(f₂) ? 0.5*rand() : f₂ - - t₁ = 0:0.001:h₁ - t₂ = 0:0.001:h₂ - t₃ = 0:0.001:h₃ - t₄ = 0:0.001:h₄ - - cos₁ = 0.5.*f₁.*(cos.(π.*(t₁./h₁)) .- 1) - cos₂ = ((0.5*(f₁+1)).*(-cos.(π.*(t₂./h₂)) .+ 1)) .- f₁ - cos₃ = ((0.5*(f₂+1)).*(cos.(π.*(t₃./h₃)))) .+ ((0.5*(1-f₂))) - cos₄ = (0.5.*f₂.*(-cos.(π.*(t₄./h₄)) .+ 1)) .- f₂ - - hrf = vcat(cos₁, cos₂, cos₃, cos₄) - -end - -""" -Lindquist et al. 2012. - -All parameters are in seconds. -Output is in ms. -If you specify a change in any of the parameters, also change the tlen to ensure the kernel is wide enough. -It will throw an error if you don't. - -""" - -function HRFDoubleGamma(; A=nothing, α₁=nothing, α₂=nothing, β₁=nothing, β₂=nothing, c=nothing, tlen=nothing) - - if (!isnothing(A) || !isnothing(α₁) || !isnothing(α₂) || !isnothing(β₁) || !isnothing(β₂) || !isnothing(c)) && isnothing(tlen) - throw(DomainError(nothing, "If you specify a different parameter, you need to specify tlen as well to ensure the kernel is wide enough.")) - end - - A = isnothing(A) ? 1 : A - α₁ = isnothing(α₁) ? 6 : α₁ - α₂ = isnothing(α₂) ? 16 : α₂ - β₁ = isnothing(β₁) ? 1 : β₁ - β₂ = isnothing(β₂) ? 1 : β₂ - c = isnothing(c) ? 1.0/6.0 : c - tlen = isnothing(tlen) ? 30 : tlen - - t = 0:0.001:tlen - - hrf = A.*((((t.^(α₁-1)).*(β₁^α₁).*exp.(-β₁.*t))./gamma(α₁)) .- (c.*((t.^(α₂-1)).*(β₂^α₂).*exp.(-β₂.*t)./ gamma(α₂)))) - -end - - -""" -### Input variables ### -adaptation of the Hemodynamics blox in fmri.jl -""" -struct BalloonModel <: ObserverBlox - params - output - jcn - odesystem - namespace - function BalloonModel(;name, namespace=nothing, lnκ=0.0, lnτ=0.0) - #= hemodynamic parameters - H(1) - signal decay d(ds/dt)/ds) - H(2) - autoregulation d(ds/dt)/df) - H(3) - transit time (t0) - H(4) - exponent for Fout(v) (alpha) - H(5) - resting state oxygen extraction (E0) - =# - - H = [0.64, 0.32, 2.00, 0.32, 0.4] - p = progress_scope(lnκ, lnτ) # progress scope if needed - p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters - lnκ, lnτ = p # assign the modified parameters - - sts = @variables s(t)=1.0 lnf(t)=1.0 lnν(t)=1.0 [output=true, description="hemodynamic_observer"] lnq(t)=1.0 [output=true, description="hemodynamic_observer"] jcn(t)=0.0 [input=true] - - eqs = [ - D(s) ~ jcn - H[1]*exp(lnκ)*s - H[2]*(exp(lnf) - 1), - D(lnf) ~ s / exp(lnf), - D(lnν) ~ (exp(lnf) - exp(lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(lnν)), - D(lnq) ~ (exp(lnf)/exp(lnq)*((1 - (1 - H[5])^(exp(lnf)^-1))/H[5]) - exp(lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) - ] - sys = System(eqs, name=name) - new(p, Num(0), sts[5], sys, namespace) - end -end - - -# This doesn't work for some good reasons. You can't call system from graph and expect to call it again later -struct CompoundHemo <: CompoundNOBlox - params - output - jcn - odesystem - namespace - function CompoundHemo(massChoice; name, lnκ=0.0, lnτ=0.0) #ONLY WORKS WITH DEFAULT PARAMETERS FOR NOW - p = progress_scope(lnκ, lnτ) - @named hemo = BalloonModel(;lnκ=p[1], lnτ=p[2]) - @named nmm = massChoice() - @variables jcn(t) - g = MetaDiGraph() - add_blox!(g, nmm) - add_blox!(g, hemo) - add_edge!(g, 1, 2, Dict(:weight => 1.0)) - linhemo = system_from_graph(g; name=name) - new(p, states(linhemo)[1], states(linhemo)[2], linhemo, nothing) - end -end - -struct LinHemoCombo <: CompoundNOBlox - params - output - jcn - odesystem - namespace - function LinHemoCombo(;name, lnκ=0.0, lnτ=0.0) - p = progress_scope(lnκ, lnτ) - p = compileparameterlist(lnκ=p[1], lnτ=p[2]) # finally compile all parameters - lnκ, lnτ = p # assign the modified parameters - H = [0.64, 0.32, 2.00, 0.32, 0.4] - - sts = @variables nmm₊x(t)=0.0 [output=true] nmm₊jcn(t)=0.0 [input=true] hemo₊s(t)=1.0 hemo₊lnf(t)=1.0 hemo₊lnν(t)=1.0 hemo₊lnq(t)=1.0 hemo₊jcn(t)=0.0 - eqs = [D(nmm₊x) ~ nmm₊jcn, - D(hemo₊s) ~ nmm₊x - H[1]*exp(lnκ)*hemo₊s - H[2]*(exp(hemo₊lnf) - 1), - D(hemo₊lnf) ~ hemo₊s / exp(hemo₊lnf), - D(hemo₊lnν) ~ (exp(hemo₊lnf) - exp(hemo₊lnν)^(H[4]^-1)) / (H[3]*exp(lnτ)*exp(hemo₊lnν)), - D(hemo₊lnq) ~ (exp(hemo₊lnf)/exp(hemo₊lnq)*((1 - (1 - H[5])^(exp(hemo₊lnf)^-1))/H[5]) - exp(hemo₊lnν)^(H[4]^-1 - 1))/(H[3]*exp(lnτ)) - ] - sys = System(eqs, name=name) - - new(p, sts[1], sts[2], sys, nothing) - end -end - - -struct AlternativeBalloonModel <: ObserverBlox - params - output - jcn - odesystem - namespace - function AlternativeBalloonModel(;name, κ=0.65, α=0.32, τ=0.98, ρ=0.34, V₀=0.02, γ=0.41) - p = progress_scope(@parameters κ=κ α=α τ=τ ρ=ρ V₀=V₀ γ=γ) # progress scope if needed - κ, α, τ, ρ, V₀, γ = p # assign the modified parameters - sts = @variables s(t)=1.0 f(t)=1.0 v(t)=1.0 q(t)=1.0 b(t)=0.0 [output=true] jcn(t)=0.0 [input=true] - eqs = [ - D(s) ~ jcn/1e2 - κ*s - γ*(f - 1), - D(f) ~ s, - D(v) ~ (f - v^(1/α))/τ, - D(q) ~ (((f*(1 - (1 - ρ^(1/f))/ρ)) - (v^(1/α)*q)/v))/τ, - b ~ V₀ * ((7*ρ*(1 - q)) + (2*(1-q/v)) + ((2*ρ - 0.2)*(1 - v))) - ] - sys = System(eqs, name=name) - new(p, Num(0), sts[5], sys, nothing) - end -end - -mutable struct JRHemo <: AbstractComponent - connector::Num - bloxinput::Num - odesystem::ODESystem - function JRHemo(;name, lnκ=0.0, lnτ=0.0) - params = progress_scope(lnκ, lnτ) - @named hemo = Hemodynamics(;lnκ=params[1], lnτ=params[2]) - @named nmm = JansenRitCBlox() - @variables jcn(t) - - g = MetaDiGraph() - add_vertex!(g, Dict(:blox => nmm, :jcn => jcn)) - add_vertex!(g, :blox, hemo) - add_edge!(g, 1, 2, :weight, 1.0) - linhemo = ODEfromGraph(g; name=name) - new(linhemo.nmm₊x, linhemo.jcn, linhemo) - end -end \ No newline at end of file diff --git a/test/data_fitting_extra_examples.jl b/test/data_fitting_extra_examples.jl deleted file mode 100644 index 44842aed..00000000 --- a/test/data_fitting_extra_examples.jl +++ /dev/null @@ -1,25 +0,0 @@ -using Neuroblox, MetaGraphs -import ModelingToolkit: outputs - -@named ob1 = BalloonModel() -@named ob2 = BalloonModel() -@named nmm1 = LinearNeuralMass() -@named nmm2 = LinearNeuralMass() - -blox = [nmm1, ob1, nmm2, ob2] - -g = MetaDiGraph() -add_blox!.(Ref(g), blox) - -add_edge!(g, 1, 2, Dict(:weight => 1.0)) # Connect hemodynamic observer #1 -add_edge!(g, 3, 4, Dict(:weight => 1.0)) # Connect hemodynamic observer #2 -add_edge!(g, 1, 3, Dict(:weight => 1.0)) # Connect neural mass -add_edge!(g, 3, 1, Dict(:weight => 1.0)) # Bidirectional connection because that's the assumption in spectral spectralDCM - -@named final_system = system_from_graph(g) - -obs_idx, obs_state_names = get_hemodynamic_observers(final_system) - -# To use obs_idx -all_outputs = outputs(final_system) -obs_outputs = all_outputs[obs_idx] #equivalent to obs_state_names \ No newline at end of file diff --git a/test/data_fitting_v2.jl b/test/data_fitting_v2.jl deleted file mode 100644 index 9c4411a9..00000000 --- a/test/data_fitting_v2.jl +++ /dev/null @@ -1,207 +0,0 @@ -using Neuroblox, Test, Graphs, MetaGraphs, OrderedCollections, LinearAlgebra, DataFrames -using MAT - -### Load data ### -vars = matread(joinpath(@__DIR__, "spectralDCM_toydata.mat")); -data = DataFrame(vars["data"], :auto) # turn data into DataFrame -x = vars["x"] # initial conditions -nd = ncol(data) # number of dimensions - -########## assemble the model ########## - -# @parameters κ=0.0 # define brain-wide decay parameter for hemodynamics -# g = MetaDiGraph() -# for ii = 1:nd -# region = LinHemo(;name=Symbol("r$ii"), lnκ=κ) -# add_blox!(g, region) -# end - -# # add symbolic weights -# @parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) -# for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) -# if idx[1] == idx[2] -# add_edge!(g, idx[1], idx[2], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 -# else -# add_edge!(g, idx[1], idx[2], :weight, A[i]) -# end -# end -# # compose model -# @named neuronmodel = ODEfromGraph(g) -# neuronmodel = structural_simplify(neuronmodel) - -# # measurement model -# @named bold = boldsignal() - -# # attribute initial conditions to states -# all_s = states(neuronmodel) -# initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() -# rnames = [] -# map(x->push!(rnames, split(string(x), "₊")[1]), all_s); -# rnames = unique(rnames); -# for (i, r) in enumerate(rnames) -# for (j, s) in enumerate(all_s[r .== map(x -> x[1], split.(string.(all_s), "₊"))]) # TODO: fix this solution, it is not robust!! -# initcond[s] = x[i, j] -# end -# end - -# modelparam = OrderedDict() -# for par in parameters(neuronmodel) -# # while Symbolics.getdefaultval(par) isa Num -# # par = Symbolics.getdefaultval(par) -# # end -# modelparam[par] = Symbolics.getdefaultval(par) -# end -# # Noise parameter mean -# modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 -# modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above -# modelparam[:lnγ] = zeros(Float64, nd); # region specific observation noise -# modelparam[:C] = ones(Float64, nd); # C as in equation 3. NB: whatever C is defined to be here, it will be replaced in csd_approx. Another strange thing of SPM12... - -# for par in parameters(bold) -# modelparam[par] = Symbolics.getdefaultval(par) -# end - -# # define prior variances -# paramvariance = copy(modelparam) -# paramvariance[:C] = zeros(Float64, nd); -# paramvariance[:lnγ] = ones(Float64, nd)./64.0; -# paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; -# paramvariance[:lnβ] = ones(Float64, length(modelparam[:lnβ]))./64.0; -# for (k, v) in paramvariance -# if occursin("A[", string(k)) -# paramvariance[k] = vars["pC"][1, 1] -# elseif occursin("κ", string(k)) -# paramvariance[k] = ones(length(v))./256.0; -# elseif occursin("ϵ", string(k)) -# paramvariance[k] = 1/256.0; -# elseif occursin("τ", string(k)) -# paramvariance[k] = 1/256.0; -# end -# end - -# params = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) -# hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix -# :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector -# ) - -# csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) -# results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) - -# ### COMPARE RESULTS WITH MATLAB RESULTS ### -# @show results.F, vars["F"] -# @test results.F < vars["F"]*0.99 -# @test results.F > vars["F"]*1.01 - -########## assemble the model ########## - -g = MetaDiGraph() -# shouldn't the following be somehow dealt with a system_from_parts() call? But that does not allow for connections between subsystems it seems -# I would like to do so to create a region that is composed of different blox that are interconnected. -# can't just build systems from graphs first and then connect them with another graph...? -regions = Dict() -@parameters κ=0.0 [tunable = true] # define brain-wide decay parameter for hemodynamics -for ii = 1:nd - region = LinearNeuralMass(;name=Symbol("r$(ii)₊lm")) - add_blox!(g, region) - regions[ii] = 2ii - 1 # store index of neural mass model - # add hemodynamic observer - observer = BalloonModel(;name=Symbol("r$(ii)₊bm"), lnκ=κ) - add_blox!(g, observer) - # connect observer with neuronal signal - add_edge!(g, 2ii - 1, 2ii, Dict(:weight => 1.0)) - # region = LinHemoCombo(;name=Symbol("r$ii"), lnκ=κ) -end - -# add symbolic weights -@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) [tunable = true] -for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) - if idx[1] == idx[2] - add_edge!(g, regions[idx[1]], regions[idx[2]], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 - else - add_edge!(g, regions[idx[2]], regions[idx[1]], :weight, A[i]) - end -end - -# compose model -@named neuronmodel = system_from_graph(g) -neuronmodel = structural_simplify(neuronmodel) - -# measurement model -@named bold = boldsignal() - -# attribute initial conditions to states -all_s = states(neuronmodel) -initcond = OrderedDict{typeof(all_s[1]), eltype(x)}() -rnames = [] -# ns = Int(length(all_s)/nd) -# for ii = 1:nd -# for jj = 1:ns -# initcond[all_s[(ii-1)*ns+jj]] = 0.0 -# end -# end -# x = rand(3, 6) #WHY IS x 3x5? -map(x->push!(rnames, split(string(x), "₊")[1]), all_s); -rnames = unique(rnames); -for (i, r) in enumerate(rnames) - for (j, s) in enumerate(all_s[r .== map(x -> x[1], split.(string.(all_s), "₊"))]) # TODO: fix this solution, it is not robust!! - initcond[s] = x[i, j] - end -end - -modelparam = OrderedDict() -for par in parameters(neuronmodel) - # if Symbolics.getdefaultval(par) isa Num - # ex = Symbolics.getdefaultval(par) - # @show ex - # p = only(Symbolics.get_variables(ex)) - # @show p - # # Symbolics.value(Symbolics.substitute(Symbolics.getdefaultval(parameters(neuronmodel2)[1]), Dict(p => Symbolics.getdefaultval(p)))) - # par = Symbolics.substitute(ex, Dict(p => Symbolics.getdefaultval(p))) - # modelparam[p] = Symbolics.value(par) - # else - if istunable(par) - modelparam[par] = Symbolics.getdefaultval(par) - end -end -# Noise parameter mean -modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 -modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above -modelparam[:lnγ] = zeros(Float64, nd); # region specific observation noise -modelparam[:C] = ones(Float64, nd); # C as in equation 3. NB: whatever C is defined to be here, it will be replaced in csd_approx. Another strange thing of SPM12... - -for par in parameters(bold) - modelparam[par] = Symbolics.getdefaultval(par) -end - -# define prior variances -paramvariance = copy(modelparam) -# [paramvariance[k] = 0.0 for k in keys(paramvariance) if occursin("w_lm", string(k)) && occursin("bm", string(k))] -paramvariance[:C] = zeros(Float64, nd); -paramvariance[:lnγ] = ones(Float64, nd)./64.0; -paramvariance[:lnα] = ones(Float64, length(modelparam[:lnα]))./64.0; -paramvariance[:lnβ] = ones(Float64, length(modelparam[:lnβ]))./64.0; -for (k, v) in paramvariance - if occursin("A[", string(k)) - paramvariance[k] = vars["pC"][1, 1] - elseif occursin("κ", string(k)) - paramvariance[k] = ones(length(v))./256.0; - elseif occursin("ϵ", string(k)) - paramvariance[k] = 1/256.0; - elseif occursin("τ", string(k)) - paramvariance[k] = 1/256.0; - end -end - -params = DataFrame(name=[k for k in keys(modelparam)], mean=[m for m in values(modelparam)], variance=[v for v in values(paramvariance)]) -hyperparams = Dict(:Πλ_pr => vars["ihC"]*ones(1,1), # prior metaparameter precision, needs to be a matrix - :μλ_pr => [vars["hE"]] # prior metaparameter mean, needs to be a vector - ) - -csdsetup = Dict(:p => 8, :freq => vec(vars["Hz"]), :dt => vars["dt"]) -results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperparams) - - -### COMPARE RESULTS WITH MATLAB RESULTS ### -@show results.F, vars["F"] -@test results.F < vars["F"]*0.99 -@test results.F > vars["F"]*1.01 \ No newline at end of file diff --git a/test/datafitting.jl b/test/datafitting.jl index c6ae1cc9..f3a7de5c 100644 --- a/test/datafitting.jl +++ b/test/datafitting.jl @@ -9,24 +9,32 @@ nd = ncol(data) # number of dimensions ########## assemble the model ########## -@parameters κ=0.0 # define brain-wide decay parameter for hemodynamics g = MetaDiGraph() +regions = Dict() +@parameters κ=0.0 [tunable = true] # define brain-wide decay parameter for hemodynamics for ii = 1:nd - region = LinHemo(;name=Symbol("r$ii"), lnκ=κ) + region = LinearNeuralMass(;name=Symbol("r$(ii)₊lm")) add_blox!(g, region) + regions[ii] = 2ii - 1 # store index of neural mass model + # add hemodynamic observer + observer = BalloonModel(;name=Symbol("r$(ii)₊bm"), lnκ=κ) + add_blox!(g, observer) + # connect observer with neuronal signal + add_edge!(g, 2ii - 1, 2ii, Dict(:weight => 1.0)) end # add symbolic weights -@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) +@parameters A[1:length(vars["pE"]["A"])] = vec(vars["pE"]["A"]) [tunable = true] for (i, idx) in enumerate(CartesianIndices(vars["pE"]["A"])) if idx[1] == idx[2] - add_edge!(g, idx[1], idx[2], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 + add_edge!(g, regions[idx[1]], regions[idx[2]], :weight, -exp(A[i])/2) # treatement of diagonal elements in SPM12 else - add_edge!(g, idx[1], idx[2], :weight, A[i]) + add_edge!(g, regions[idx[2]], regions[idx[1]], :weight, A[i]) end end + # compose model -@named neuronmodel = ODEfromGraph(g) +@named neuronmodel = system_from_graph(g) neuronmodel = structural_simplify(neuronmodel) # measurement model @@ -46,11 +54,11 @@ end modelparam = OrderedDict() for par in parameters(neuronmodel) - # while Symbolics.getdefaultval(par) isa Num - # par = Symbolics.getdefaultval(par) - # end - modelparam[par] = Symbolics.getdefaultval(par) + if istunable(par) + modelparam[par] = Symbolics.getdefaultval(par) + end end + # Noise parameter mean modelparam[:lnα] = [0.0, 0.0]; # intrinsic fluctuations, ln(α) as in equation 2 of Friston et al. 2014 modelparam[:lnβ] = [0.0, 0.0]; # global observation noise, ln(β) as above @@ -90,4 +98,4 @@ results = spectralVI(data, neuronmodel, bold, initcond, csdsetup, params, hyperp ### COMPARE RESULTS WITH MATLAB RESULTS ### @show results.F, vars["F"] @test results.F < vars["F"]*0.99 -@test results.F > vars["F"]*1.01 +@test results.F > vars["F"]*1.01 \ No newline at end of file