From d7b3ba2ece3ccddc862b0b5366cff8434a94b6d1 Mon Sep 17 00:00:00 2001 From: Imre Kertesz <36957615+ImreKertesz@users.noreply.github.com> Date: Tue, 3 Dec 2024 11:07:30 +0100 Subject: [PATCH] rDCM toolbox v1.5.0 (#290) * replace deprecated caxis with clim * clean up code * remove r_dt factor and fix bug in calculation of free energy for sparse rDCM --- Contributor-License-Agreement.md | 1 + rDCM/CHANGELOG.md | 9 +++++++++ rDCM/README.md | 17 +++++++++-------- rDCM/code/tapas_rdcm_get_convolution_bm.m | 1 - rDCM/code/tapas_rdcm_ridge.m | 5 ++--- rDCM/code/tapas_rdcm_sparse.m | 6 +++--- rDCM/code/tapas_rdcm_spm_dcm_fmri_priors.m | 6 +++--- rDCM/code/tapas_rdcm_spm_dcm_generate.m | 8 ++++---- rDCM/code/tapas_rdcm_spm_logdet.m | 3 +-- rDCM/code/tapas_rdcm_to_spm_struct.m | 2 ++ rDCM/code/tapas_rdcm_visualize.m | 16 ++++++++-------- 11 files changed, 42 insertions(+), 32 deletions(-) diff --git a/Contributor-License-Agreement.md b/Contributor-License-Agreement.md index 5ec43e1b..5e236044 100644 --- a/Contributor-License-Agreement.md +++ b/Contributor-License-Agreement.md @@ -25,6 +25,7 @@ Matthias Müller-Schrader | TNU, University of Zurich | Zurich | CH Saskia Bollmann | The University of Queensland | Brisbane | AUS | SaskiaBollmann | 1.1 Anna Yurova | LMU University Hospital | Munich | GER | a-yur | 1.1 Irene Sophia Plank | LMU University Hospital | Munich | GER | IreneSophia | 1.1 +Imre Kertesz | TNU, University of Zurich | Zurich | CH | ImreKertesz | 1.1 **- Add Entry here -** | **- Add Entry here -** | **Add** | **Add** | **Add** | 1.1 (hereinafter referred to as "Contributor") diff --git a/rDCM/CHANGELOG.md b/rDCM/CHANGELOG.md index a99e03bf..d284629c 100755 --- a/rDCM/CHANGELOG.md +++ b/rDCM/CHANGELOG.md @@ -1,6 +1,15 @@ # Changelog Regression Dynamic Causal Modeling (rDCM) toolbox +## [1.5] 2024-12-03 +### Added +none +### Changed +- Fix bug in calculation of negative free energy for sparse rDCM +- remove r_dt from update equations +- clean up code +### Removed +none ## [1.4] 2022-03-22 diff --git a/rDCM/README.md b/rDCM/README.md index 9c39d3f6..655eb624 100755 --- a/rDCM/README.md +++ b/rDCM/README.md @@ -8,18 +8,19 @@ rDCM - regression Dynamic Causal Modeling. This ReadMe file provides the relevant information on the regression dynamic causal modeling (rDCM) toolbox. +> **_NOTE:_** There is a new implementation of rDCM in Julia which can be found [here](https://github.com/ComputationalPsychiatry/RegressionDynamicCausalModeling.jl). ------------------- General information ------------------- - Authors: Stefan Frässle (), Ekaterina I. Lomakina -- Copyright (C) 2016-2021 -- Translational Neuromodeling Unit (TNU) -- Institute for Biomedical Engineering -- University of Zurich & ETH Zurich - +- Updated by [Imre Kertesz](https://github.com/ImreKertesz) +- Copyright (C) 2016-2024 +Translational Neuromodeling Unit (TNU) +
Institute for Biomedical Engineering +
University of Zurich & ETH Zurich -------- Download @@ -37,7 +38,7 @@ Download Purpose ------- -The regression dynamic causal modeling (rDCM) toobox implements a novel variant +The regression dynamic causal modeling (rDCM) toolbox implements a novel variant of DCM for fMRI that enables computationally efficient inference on effective (i.e., directed) connectivity parameters among brain regions. Due to its computational efficiency, inversion of large (whole-brain) networks becomes feasible. @@ -65,7 +66,7 @@ Important Notes --------------- Please note that rDCM is a method that is still in an early stage of development and thus subject to -various limiations. Due to these limitations, requirements of rDCM in terms of +various limitations. Due to these limitations, requirements of rDCM in terms of fMRI data quality (i.e., fast TR, high SNR) are high. For data that does not meet these conditions, the method might not give reliable results. It remains the responsibility of each user to ensure that his/her dataset fulfills these @@ -132,6 +133,6 @@ Acknowledgment -------------- We would like to highlight and acknowledge that the rDCM toolbox uses some -functions that were publised as part of the Statistical Parameteric Mapping +functions that were published as part of the Statistical Parametric Mapping ([SPM](https://www.fil.ion.ucl.ac.uk/spm/software/)) toolbox. The respective functions are marked with the prefix `tapas_rdcm_spm_`. diff --git a/rDCM/code/tapas_rdcm_get_convolution_bm.m b/rDCM/code/tapas_rdcm_get_convolution_bm.m index ca986d9a..81bb54f5 100755 --- a/rDCM/code/tapas_rdcm_get_convolution_bm.m +++ b/rDCM/code/tapas_rdcm_get_convolution_bm.m @@ -71,6 +71,5 @@ % sample the HRF at the sampling rate of the data h = y(1:r_dt:end,1); -h = h; end diff --git a/rDCM/code/tapas_rdcm_ridge.m b/rDCM/code/tapas_rdcm_ridge.m index f52d1608..38892c57 100755 --- a/rDCM/code/tapas_rdcm_ridge.m +++ b/rDCM/code/tapas_rdcm_ridge.m @@ -84,12 +84,11 @@ Y_r = Y(idx_y,k); % effective number of data points - N_eff = sum(idx_y)/(args.r_dt); + N_eff = sum(idx_y); % effective dimensionality D_r = sum(idx(k,:)); - %% read priors % prior precision matrix @@ -109,7 +108,7 @@ t = a0/b0; % estimate alpha_N - aN_r = a0 + N_eff/ (2*args.r_dt); + aN_r = a0 + N_eff/2; % cycle stops after 500 iterations count = 500; diff --git a/rDCM/code/tapas_rdcm_sparse.m b/rDCM/code/tapas_rdcm_sparse.m index e8f2e9ee..3828cafa 100755 --- a/rDCM/code/tapas_rdcm_sparse.m +++ b/rDCM/code/tapas_rdcm_sparse.m @@ -181,7 +181,7 @@ % intialize z, t & aN per regions z_r = p0; t = a0/b0; - aN_r = a0 + N_eff/(2*args.r_dt); + aN_r = a0 + N_eff/2; % cycle stops after 500 iterations; @@ -263,7 +263,7 @@ %% compute model evidence % compute the components of the model evidence - log_lik = N_eff*(psi(aN_r) - log(bN_r))/2 - N_eff*log(2*pi)/2 - t*QF/2; + log_lik = N_eff*(psi(aN_r) - log(bN_r))/2 - N_eff*log(2*pi)/2 - t*QF; log_p_weight = 1/2*tapas_rdcm_spm_logdet(l0_r) - D*log(2*pi)/2 - (mN_r-m0_r)'*l0_r*(mN_r-m0_r)/2 - trace(l0_r*sN_r)/2; log_p_prec = a0*log(b0) - gammaln(a0) + (a0-1)*(psi(aN_r) - log(bN_r)) - b0*t; log_p_z = sum(log(1-p0(z_idx)) + z_r(z_idx).*log(p0(z_idx)./(1-p0(z_idx)))); @@ -303,7 +303,7 @@ %% compute model evidence % compute the components of the model evidence - log_lik_iter(iter) = N_eff*(psi(aN_r) - log(bN_r))/2 - N_eff*log(2*pi)/2 - t*QF/2; + log_lik_iter(iter) = N_eff*(psi(aN_r) - log(bN_r))/2 - N_eff*log(2*pi)/2 - t*QF; log_p_weight_iter(iter) = 1/2*tapas_rdcm_spm_logdet(l0_r) - D*log(2*pi)/2 - (mN_r-m0_r)'*l0_r*(mN_r-m0_r)/2 - trace(l0_r*sN_r)/2; log_p_prec_iter(iter) = a0*log(b0) - gammaln(a0) + (a0-1)*(psi(aN_r) - log(bN_r)) - b0*t; log_p_z_iter(iter) = sum(log(1-p0(z_idx)) + z_r(z_idx).*log(p0(z_idx)./(1-p0(z_idx)))); diff --git a/rDCM/code/tapas_rdcm_spm_dcm_fmri_priors.m b/rDCM/code/tapas_rdcm_spm_dcm_fmri_priors.m index 8cf004af..6e922c73 100755 --- a/rDCM/code/tapas_rdcm_spm_dcm_fmri_priors.m +++ b/rDCM/code/tapas_rdcm_spm_dcm_fmri_priors.m @@ -35,9 +35,9 @@ % check options and D (for nonlinear coupling) %-------------------------------------------------------------------------- -try, options.two_state; catch, options.two_state = 0; end -try, options.endogenous; catch, options.endogenous = 0; end -try, D; catch, D = zeros(n,n,0); end +try options.two_state; catch, options.two_state = 0; end +try options.endogenous; catch, options.endogenous = 0; end +try D; catch, D = zeros(n,n,0); end % prior (initial) states and shrinkage priors on A for endogenous DCMs diff --git a/rDCM/code/tapas_rdcm_spm_dcm_generate.m b/rDCM/code/tapas_rdcm_spm_dcm_generate.m index 8552f1eb..deee59e4 100755 --- a/rDCM/code/tapas_rdcm_spm_dcm_generate.m +++ b/rDCM/code/tapas_rdcm_spm_dcm_generate.m @@ -23,7 +23,7 @@ if isstruct(syn_model) DCM = syn_model; else - load(syn_model) + load(syn_model,"DCM") end if nargin < 2 || isempty(SNR) SNR = 1; @@ -35,7 +35,7 @@ U = DCM.U; % inputs v = DCM.v; % number of scans n = DCM.n; % number of regions -m = size(U.u,2); % number of inputs +% m = size(U.u,2); % number of inputs % Check whether the model is stable by examining the eigenvalue @@ -81,8 +81,8 @@ % fMRI slice time sampling %-------------------------------------------------------------------------- -try, M.delays = DCM.delays; end -try, M.TE = DCM.TE; end +try M.delays = DCM.delays; catch, end +try M.TE = DCM.TE; catch, end % Integrate and compute hemodynamic response at v sample points diff --git a/rDCM/code/tapas_rdcm_spm_logdet.m b/rDCM/code/tapas_rdcm_spm_logdet.m index d92ded05..668bad03 100644 --- a/rDCM/code/tapas_rdcm_spm_logdet.m +++ b/rDCM/code/tapas_rdcm_spm_logdet.m @@ -16,7 +16,6 @@ % assume diagonal form %-------------------------------------------------------------------------- TOL = 1e-16; -n = length(C); s = diag(C); i = find(s > TOL & s < 1/TOL); C = C(i,i); @@ -24,7 +23,7 @@ % invoke det if non-diagonal %-------------------------------------------------------------------------- -[i j] = find(C); +[i, j] = find(C); if any(i ~= j) n = length(C); a = exp(H/n); diff --git a/rDCM/code/tapas_rdcm_to_spm_struct.m b/rDCM/code/tapas_rdcm_to_spm_struct.m index 855aac03..f0d3722b 100644 --- a/rDCM/code/tapas_rdcm_to_spm_struct.m +++ b/rDCM/code/tapas_rdcm_to_spm_struct.m @@ -200,6 +200,8 @@ try DCM.y = abs(reshape(output.signal.yd_pred_rdcm_fft,DCM.v,DCM.n)).^2; DCM.R = abs(reshape(output.signal.yd_source_fft,DCM.v,DCM.n)).^2 - DCM.y; +catch + end % predicted data (temporal derivative in Fourier space) diff --git a/rDCM/code/tapas_rdcm_visualize.m b/rDCM/code/tapas_rdcm_visualize.m index 9378f9ee..d789d6e8 100755 --- a/rDCM/code/tapas_rdcm_visualize.m +++ b/rDCM/code/tapas_rdcm_visualize.m @@ -87,7 +87,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(output.Ep.A) title('estimated','FontSize',14) axis square - caxis([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) + clim([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) @@ -100,7 +100,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(DCM.Tp.A) title('true','FontSize',14) axis square - caxis([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) + clim([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) @@ -151,7 +151,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(output.Ep.A) title('estimated','FontSize',16) axis square - caxis([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) + clim([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) @@ -164,7 +164,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(DCM.Tp.A) title('true','FontSize',16) axis square - caxis([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) + clim([-1*max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A))))) max(max(abs(DCM.Tp.A-diag(diag(DCM.Tp.A)))))]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) @@ -184,7 +184,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(output.Ep.A) title('estimated','FontSize',14) axis square - caxis([-1*max(max(abs(output.Ep.A-diag(diag(output.Ep.A))))) max(max(abs(output.Ep.A-diag(diag(output.Ep.A)))))]) + clim([-1*max(max(abs(output.Ep.A-diag(diag(output.Ep.A))))) max(max(abs(output.Ep.A-diag(diag(output.Ep.A)))))]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) @@ -197,7 +197,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(output.Ip.A) title('Pp binary indicator','FontSize',14) axis square - caxis([0 1]) + clim([0 1]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) @@ -249,7 +249,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(output.Ep.A) title('estimated','FontSize',16) axis square - caxis([-1*max(max(abs(output.Ep.A-diag(diag(output.Ep.A))))) max(max(abs(output.Ep.A-diag(diag(output.Ep.A)))))]) + clim([-1*max(max(abs(output.Ep.A-diag(diag(output.Ep.A))))) max(max(abs(output.Ep.A-diag(diag(output.Ep.A)))))]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) @@ -262,7 +262,7 @@ function tapas_rdcm_visualize(output, DCM, options, plot_regions, plot_mode) imagesc(output.Ip.A) title('Pp binary indicator','FontSize',16) axis square - caxis([0 1]) + clim([0 1]) colorbar set(gca,'xtick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)]) set(gca,'ytick',[1 round(size(output.Ep.A,1)/2) size(output.Ep.A,1)])