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

@@ -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")

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
 

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 (<[email protected]>), 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)
+<br>Institute for Biomedical Engineering
+<br>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_`.

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

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; 

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))));

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

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

rDCM/code/tapas_rdcm_spm_logdet.m

@@ -16,15 +16,14 @@
 % assume diagonal form
 %--------------------------------------------------------------------------
 TOL   = 1e-16;
-n     = length(C);
 s     = diag(C);
 i     = find(s > TOL & s < 1/TOL);
 C     = C(i,i);
 H     = sum(log(diag(C)));
 
 % invoke det if non-diagonal
 %--------------------------------------------------------------------------
-[i j] = find(C);
+[i, j] = find(C);
 if any(i ~= j)
       n = length(C);
       a = exp(H/n);

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)

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)])