diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index f2c1332f..00000000 Binary files a/.DS_Store and /dev/null differ diff --git a/.gitignore b/.gitignore index d9282a28..d43b1e47 100644 --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,4 @@ lib *.so *.mexw64 *.asv +*.DS_Store diff --git a/CHANGELOG.md b/CHANGELOG.md index c73983ff..0b376df4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,25 @@ # Changelog TAPAS toolbox +## [3.2.0] 2019-09-dd +### Added +- HUGE: introduced object-oriented interface in addition to old interface +- HUGE: build-in unit tests +- HUGE: user manual +- PhysIO (details in tapas/PhysIO/CHANGELOG.md) + - more unit testing and integration testing for examples + - bandpass-filtering for cardiac data in preprocessing, user-defined max + heart rate for peak detection + +### Fixed +- PhysIO: Bugfix RVT/HRV convolution had erroneous half-width shift + - For details on the bug, its impact and fix, see our specific [Release + Info on RVT/HRV Bugfix](https://github.com/translationalneuromodeling/tapas/issues/65) + +### Changed +- HUGE: demo script reflect interface changes + + ## [3.1.0] 2019-03-26 ### Added diff --git a/PhysIO/CHANGELOG.md b/PhysIO/CHANGELOG.md index 6ec3c29f..fc8ca1e1 100644 --- a/PhysIO/CHANGELOG.md +++ b/PhysIO/CHANGELOG.md @@ -4,11 +4,46 @@ RELEASE INFORMATION Current Release --------------- -*Current version: PhysIO Toolbox Release R2019a, v7.1.0* +*Current version: PhysIO Toolbox Release R2019b, v7.2.0* -March 19, 2019 +August 20, 2019 +Minor Release Notes (R2019b, v7.2.0) +------------------------------------ + +### Added +- `requirements.txt` making dependencies on Matlab and specific toolboxes + explicit +- `max_heart_rate_bpm` now a user parameter to adjust prior on max allowed + heart rate for cardiac pulse detection (`method = 'auto_matched'`) +- bandpass-filtering of cardiac data during preprocessing now possible + (`preproc.cardiac.filter`) +- Added integration tests for all examples in `tests/integration` for both SPM +Batch Editor and Matlab-only configuration scripts. Reference data provided in +`examples/TestReferenceResults/examples` + +### Changed +- Toned down and replaced irrelevant peak height and missing cardiac pulse + warnings (github issue #51) +- Updated README to include external contributors and recent findings about + impact of physiological noise for serial correlations (Bollmann2018) +- Added unit tests for convolution and moved all to `tests/unit` + +### Fixed +- Corrected half-width shift of response functions for HRV and RVT regressors by + erroneous use of Matlab `conv` + - For details on the bug, its impact and fix, see our specific [Release + Info on RVT/HRV Bugfix](https://github.com/translationalneuromodeling/tapas/issues/65) + - other references: TNU gitlab issue #83, found and fixed by Sam Harrison, + TNU, see `tapas_physio_conv`) +- Bugfix `tapas_physio_init()` not working, because dependent on functions + in `utils` subfolder not in path; `utils` added to path +- `tapas_physio_review` for motion parameters (found and fixed by + Sam Harrison, TNU) +- visualization error for regressor orthogonalization (github issue #57), + when only `'RETROICOR'` set was chosen + Minor Release Notes (R2019a, v7.1.0) ------------------------------------ diff --git a/PhysIO/README.md b/PhysIO/README.md index 76502a41..108c0c1a 100644 --- a/PhysIO/README.md +++ b/PhysIO/README.md @@ -1,7 +1,7 @@ TAPAS PhysIO Toolbox ==================== -*Current version: Release 2019a, v7.1.0* +*Current version: Release 2019b, v7.2.0* > Copyright (C) 2012-2019 > Lars Kasper @@ -76,7 +76,8 @@ Getting Started ...following the installation, you can try out an example: -1. Download the TAPAS examples via running `tapas_download_example_data()` (found in `misc`-subfolder of TAPAS) +1. Download the TAPAS examples via running `tapas_download_example_data()` + (found in `misc`-subfolder of TAPAS) - The PhysIO Example files will be downloaded to `tapas/examples//PhysIO` 2. Run `philips_ecg3t_matlab_script.m` in subdirectory `Philips/ECG3T` 3. See subdirectory `physio/docs` and the next two section of this document for help. @@ -98,7 +99,9 @@ pointers and templates. Before you contact us, please try the following: [tapas@sympa.ethz.ch](https://sympa.ethz.ch/sympa/info/tapas), which has a searchable [archive](https://sympa.ethz.ch/sympa/arc/tapas). 3. For new requests, we would like to ask you to submit them as - [issues](https://github.com/translationalneuromodeling/tapas/issues) on our github release page for TAPAS, which is also an up-to-date resource to user-driven questions (since 2018). + [issues](https://github.com/translationalneuromodeling/tapas/issues) on our + github release page for TAPAS, which is also an up-to-date resource to + user-driven questions (since 2018). Documentation @@ -144,8 +147,14 @@ Facts about physiological noise in fMRI: become a particular concern at and above 3 Tesla (Kasper2009, Hutton2011). - In resting state fMRI, disregarding physiological noise leads to wrong connectivity results (Birn2006). +- Uncorrected physiological noise introduces serial correlations into the residual + voxel time series, that invalidate assumptions on noise correlations (e.g., AR(1)) + used in data prewhitening by all major analysis packages. This issue is particularly + aggravated at short TR (<1s), and most of its effects can be suitably addressed + by physiological noise correction (Bollmann2018) -Therefore, some kind of physiological noise correction is highly recommended for every statistical fMRI analysis. +Therefore, some kind of physiological noise correction is highly recommended for +every statistical fMRI analysis. Model-based correction of physiological noise: - Physiological noise can be decomposed into periodic time series following @@ -174,30 +183,39 @@ Features of this Toolbox - Flexible expansion orders to model different contributions of cardiac, respiratory and interaction terms (see Harvey2008, Hutton2011) - Data-driven noise regressors - - PCA extraction from nuisance ROIs (CSF, white matter), similar to aCompCor (Behzadi2007) + - PCA extraction from nuisance ROIs (CSF, white matter), similar to aCompCor + (Behzadi2007) ### Automatization and Performance Assessment - Automatic creation of nuisance regressors, full integration into standard GLMs, tested for SPM8/12 ("multiple_regressors.mat") -- Integration in SPM Batch Editor: GUI for parameter input, dependencies to integrate physiological noise correction in preprocessing pipeline -- Performance Assessment: Automatic F-contrast and tSNR Map creation and display for groups of physiological noise regressors, using SPM GLM tools +- Integration in SPM Batch Editor: GUI for parameter input, dependencies to + integrate physiological noise correction in preprocessing pipeline +- Performance Assessment: Automatic F-contrast and tSNR Map creation and display + for groups of physiological noise regressors, using SPM GLM tools via + `tapas_physio_report_contrasts()`. ### Flexible Read-in -The toolbox is dedicated to seamless integration into a clinical research s -etting and therefore offers correction methods to recover physiological -data from imperfect peripheral measures. +The toolbox is dedicated to seamless integration into a clinical research +setting and therefore offers correction methods to recover physiological +data from imperfect peripheral measures. Read-in of the following formats is +currently supported (alphabetic order): +- Biopac `.mat` and `.txt` -export +- Brain Imaging Data Structure (BIDS) +- Custom logfiles: should contain one amplitude value per line, one logfile per + device. Sampling interval(s) are provided as a separate parameter to the toolbox. - General Electric - Philips SCANPHYSLOG files (all versions from release 2.6 to 5.3) - Siemens VB (files `.ecg`, `.resp`, `.puls`) - Siemens VD (files `*_ECG.log`, `*_RESP.log`, `*_PULS.log`) - Siemens Human Connectome Project (preprocessed files `*Physio_log.txt`) -- Biopac .mat-export - - assuming the following variables (as columns): `data`, `isi`, `isi_units`, `labels`, `start_sample`, `units` - - See `tapas_physio_read_physlogfiles_biopac_mat.m` for details -- Custom logfiles: should contain one amplitude value per line, one logfile per device. Sampling interval(s) are provided as a separate parameter to the toolbox. + +See also the +[Wiki page on Read-In](https://gitlab.ethz.ch/physio/physio-doc/wikis/MANUAL_PART_READIN) +for a more detailed list and description of the supported formats. Compatibility @@ -217,7 +235,9 @@ Compatibility or as text file for export to any other package - raw and processed physiological logfile data - Graphical Batch Editor interface to SPM -- Part of the TAPAS Software Collection of the Translational Neuromodeling Unit (TNU) Zurich:long term support and ongoing development +- Part of the TAPAS Software Collection of the Translational Neuromodeling Unit + (TNU) Zurich + - ensures long term support and ongoing development Contributors @@ -229,12 +249,43 @@ Contributors - Project Team: - Steffen Bollmann, Centre for Advanced Imaging, University of Queensland, Australia - Saskia Bollmann, Centre for Advanced Imaging, University of Queensland, Australia -- Contributors: +- Contributors (Code): - Eduardo Aponte, TNU Zurich - Tobias U. Hauser, FIL London, UK + - Sam Harrison, TNU Zurich - Jakob Heinzle, TNU Zurich - Chloe Hutton, FIL London, UK (previously) - Miriam Sebold, Charite Berlin, Germany + - External TAPAS contributors listed in its [Contributor License Agreement] + (https://github.com/translationalneuromodeling/tapas/blob/master/Contributor-License-Agreement.md) +- Contributors (Examples): + - listed in [EXAMPLES.md](https://gitlab.ethz.ch/physio/physio-doc/wikis/EXAMPLES) + + +Requirements +------------ + +- All specific software requirements and their versions are in a separate file + in this folder, `requirements.txt`. +- In brief: + - PhysIO needs Matlab to run, and some of its toolboxes. + - Some functionality requires SPM (GUI, nuisance regression, contrast reporting, + writing residual and SNR images). + + +Acknowledgements +---------------- + +The PhysIO Toolbox ships with the following publicly available code from other +open source projects and gratefully acknowledges their use. + +- `utils\tapas_physio_propval.m` + - `propval` function from Princeton MVPA toolbox (GPL) + a nice wrapper function to create flexible propertyName/value optional + parameters +- `utils\tapas_physio_fieldnamesr.m` + - recursive parser for field names of a structure + - Matlab file exchange, adam.tudorjones@pharm.ox.ac.uk References @@ -242,55 +293,79 @@ References ### Main Toolbox Reference +Please cite the following paper in all of your publications that utilized the +PhysIO Toolbox. + 1. Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.-M., Pruessmann, K.P., Stephan, K.E., 2017. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. -Journal of Neuroscience Methods 276, 56–72. doi:10.1016/j.jneumeth.2016.10.019 +Journal of Neuroscience Methods 276, 56–72. https://doi.org/10.1016/j.jneumeth.2016.10.019 + +The [FAQ](https://gitlab.ethz.ch/physio/physio-doc/wikis/FAQ#3-how-do-i-cite-physio) +contains a complete suggestion for a snippet in your methods section. ### Related Papers (Implemented noise correction algorithms and optimal parameter choices) +The following sections list papers that +- first implemented specific noise correction algorithms +- determined optimal parameter choices for these algorithms, depending on the + targeted application +- demonstrate the impact of physiological noise and the importance of its correction + +It is loosely ordered by the dominant physiological noise model used in the +paper. The list is by no means complete, and we are happy to add any relevant papers +suggested to us. + #### RETROICOR 2. Glover, G.H., Li, T.Q. & Ress, D. Image‐based method for retrospective correction of PhysIOlogical motion effects in fMRI: RETROICOR. Magn Reson Med 44, 162-7 (2000). -3. Hutton, C. et al. The impact of PhysIOlogical noise correction on fMRI at 7 T. +3. Hutton, C. et al. The impact of Physiological noise correction on fMRI at 7 T. NeuroImage 57, 101‐112 (2011). 4. Harvey, A.K. et al. Brainstem functional magnetic resonance imaging: Disentangling signal from PhysIOlogical noise. Journal of Magnetic Resonance Imaging 28, 1337‐1344 (2008). +5. Bollmann, S., Puckett, A.M., Cunnington, R., Barth, M., 2018. +Serial correlations in single-subject fMRI with sub-second TR. +NeuroImage 166, 152–166. https://doi.org/10.1016/j.neuroimage.2017.10.043 + + #### aCompCor / Noise ROIs -5. Behzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. A component based noise +6. Behzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, -90–101. doi:10.1016/j.neuroimage.2007.04.042 +90–101. https://doi.org/10.1016/j.neuroimage.2007.04.042 #### RVT -6. Birn, R.M., Smith, M.A., Jones, T.B., Bandettini, P.A., 2008. The respiration response +7. Birn, R.M., Smith, M.A., Jones, T.B., Bandettini, P.A., 2008. The respiration response function: The temporal dynamics of fMRI s ignal fluctuations related to changes in respiration. NeuroImage 40, 644–654. doi:10.1016/j.neuroimage.2007.11.059 -7. Jo, H.J., Saad, Z.S., Simmons, W.K., Milbury, L.A., Cox, R.W., 2010. +8. Jo, H.J., Saad, Z.S., Simmons, W.K., Milbury, L.A., Cox, R.W., 2010. Mapping sources of correlation in resting state FMRI, with artifact detection and removal. NeuroImage 52, 571–582. https://doi.org/10.1016/j.neuroimage.2010.04.246 -*regressor delay suggestions* + - *regressor delay suggestions* #### HRV -8. Chang, C., Cunningham, J.P., Glover, G.H., 2009. Influence of heart rate on the +9. Chang, C., Cunningham, J.P., Glover, G.H., 2009. Influence of heart rate on the BOLD signal: The cardiac response function. NeuroImage 44, 857–869. doi:10.1016/j.neuroimage.2008.09.029 -9. Shmueli, K., van Gelderen, P., de Zwart, J.A., Horovitz, S.G., Fukunaga, M., +10. Shmueli, K., van Gelderen, P., de Zwart, J.A., Horovitz, S.G., Fukunaga, M., Jansma, J.M., Duyn, J.H., 2007. Low-frequency fluctuations in the cardiac rate as a source of variance in the resting-state fMRI BOLD signal. NeuroImage 38, 306–320. https://doi.org/10.1016/j.neuroimage.2007.07.037 -*regressor delay suggestions* + - *regressor delay suggestions* #### Motion (Censoring, Framewise Displacement) -10. Siegel, J.S., Power, J.D., Dubis, J.W., Vogel, A.C., Church, J.A., Schlaggar, B.L., +11. Siegel, J.S., Power, J.D., Dubis, J.W., Vogel, A.C., Church, J.A., Schlaggar, B.L., Petersen, S.E., 2014. Statistical improvements in functional magnetic resonance imaging analyses produced by censoring high-motion data points. Hum. Brain Mapp. -35, 1981–1996. doi:10.1002/hbm.22307 +35, 1981–1996. https://doi.org/10.1002/hbm.22307 -11. Power, J.D., Barnes, K.A., Snyder, A.Z., Schlaggar, B.L., Petersen, S.E., 2012. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. NeuroImage 59, 2142–2154. https://doi.org/10.1016/j.neuroimage.2011.10.018 +12. Power, J.D., Barnes, K.A., Snyder, A.Z., Schlaggar, B.L., Petersen, S.E., 2012. +Spurious but systematic correlations in functional connectivity MRI networks +arise from subject motion. NeuroImage 59, 2142–2154. +https://doi.org/10.1016/j.neuroimage.2011.10.018 Copying/License @@ -308,4 +383,4 @@ General Public License for more details. You should have received a copy of the GNU General Public License along with this program (see the file [LICENSE](LICENSE)). If not, see -. \ No newline at end of file +. diff --git a/PhysIO/code/assess/tapas_physio_review.m b/PhysIO/code/assess/tapas_physio_review.m index 759c4fcf..0c5066bc 100644 --- a/PhysIO/code/assess/tapas_physio_review.m +++ b/PhysIO/code/assess/tapas_physio_review.m @@ -131,17 +131,18 @@ hasRespData); end - if model.movement.include - censoring = model.movement.censoring; + rp = model.movement.rp; quality_measures = model.movement.quality_measures; + censoring = model.movement.censoring; + censoring_threshold = model.movement.censoring_threshold; switch lower(model.movement.censoring_method) case 'fd' - verbose.fig_handles(end+1) = tapas_physio_plot_movement_outliers_fd(rp, ... - quality_measures, censoring, movement.censoring_threshold); + verbose.fig_handles(end+1) = tapas_physio_plot_movement_outliers_fd( ... + rp, quality_measures, censoring, censoring_threshold); case 'maxval' - verbose.fig_handles(end+1) = tapas_physio_plot_movement_outliers_maxval(rp, ... - quality_measures, censoring, movement.censoring_threshold); + verbose.fig_handles(end+1) = tapas_physio_plot_movement_outliers_maxval( ... + rp, quality_measures, censoring, censoring_threshold); end end diff --git a/PhysIO/code/gui/tapas_physio_gui_filter.m b/PhysIO/code/gui/tapas_physio_gui_filter.m new file mode 100644 index 00000000..6004f35e --- /dev/null +++ b/PhysIO/code/gui/tapas_physio_gui_filter.m @@ -0,0 +1,108 @@ +function filter = tapas_physio_gui_filter() +% Creates sub-part for preproc filter options in Batch Editor GUI +% +% tapas_physio_gui_filter +% +% +% See also + +% Author: Lars Kasper +% Created: 2019-07-05 +% Copyright (C) 2019 Institute for Biomedical Engineering +% University of Zurich and ETH Zurich +% +% This file is part of the TAPAS PhysIO Toolbox, which is released under +% the terms of the GNU General Public License (GPL), version 3. You can +% redistribute it and/or modify it under the terms of the GPL (either +% version 3 or, at your option, any later version). For further details, +% see the file COPYING or . + +%-------------------------------------------------------------------------- +% filter_type +%-------------------------------------------------------------------------- +filter_type = cfg_menu; +filter_type.tag = 'type'; +filter_type.name = 'Filter Type'; +filter_type.help = {'Which infinite impulse response filter shall be used?' + ' ''Chebychev Type II (cheby2)'' Chebychev Type II filter, use for steep transition from' + ' start to stop band' + ' ''Butterworth (butter)'' Butterworth filter, standard filter with maximally flat' + ' passband (Infinite impulse response), but stronger' + ' ripples in transition band' +}; +filter_type.labels = {'Chebychev Type II', 'Butterworth'}; +filter_type.values = {'cheby2', 'butter'}; +filter_type.val = {'cheby2'}; + +%-------------------------------------------------------------------------- +% filter_stopband +%-------------------------------------------------------------------------- +filter_stopband = cfg_entry; +filter_stopband.tag = 'stopband'; +filter_stopband.name = 'Stopband'; +filter_stopband.help = { + '[f_min, f_max] frequency interval in Hz of all frequencies, s.th. frequencies' + ' outside this band should definitely *NOT* pass the filter' + ' Default: [] ' + ' NOTE: only relevant for ''cheby2'' filter type' + ' if empty, and passband is empty, no filtering is performed' + ' if empty, but passband exists, stopband interval is' + ' 10% increased passband interval' + }; +filter_stopband.strtype = 'e'; +filter_stopband.num = [0 Inf]; +filter_stopband.val = {[]}; + +%-------------------------------------------------------------------------- +% filter_passband +%-------------------------------------------------------------------------- +filter_passband = cfg_entry; +filter_passband.tag = 'passband'; +filter_passband.name = 'Passband'; +filter_passband.help = { + '[f_min, f_max] frequency interval in Hz of all frequency that should' + ' pass the passband filter' + ' default: [0.3 9] (to remove slow drifts, breathing' + ' and slice sampling artifacts)' + ' if empty, no filtering is performed' + }; +filter_passband.strtype = 'e'; +filter_passband.num = [0 Inf]; +filter_passband.val = {[0.3 9]}; + +%-------------------------------------------------------------------------- +%% filter_no +%-------------------------------------------------------------------------- + +filter_no = cfg_branch; +filter_no.tag = 'no'; +filter_no.name = 'No'; +filter_no.val = {}; +filter_no.help = {'Cardiac data remains unfiltered before pulse detection.'}; + +%-------------------------------------------------------------------------- +%% filter_yes +%-------------------------------------------------------------------------- + +filter_yes = cfg_branch; +filter_yes.tag = 'yes'; +filter_yes.name = 'Yes'; +filter_yes.val = {filter_type, filter_passband, ... +filter_stopband}; +filter_yes.help = {''}; + + +%-------------------------------------------------------------------------- +%% filter +%-------------------------------------------------------------------------- + +filter = cfg_choice; +filter.tag = 'filter'; +filter.name = 'Filter Raw Cardiac Time Series'; +filter.val = {filter_no}; +filter.values = {filter_no, filter_yes}; +filter.help = { + 'Filter properties for bandpass-filtering of cardiac signal before peak' + 'detection, phase extraction, and other physiological traces' + }; + diff --git a/PhysIO/code/model/tapas_physio_create_hrv_regressors.m b/PhysIO/code/model/tapas_physio_create_hrv_regressors.m index fb5d777e..f89f7ddb 100644 --- a/PhysIO/code/model/tapas_physio_create_hrv_regressors.m +++ b/PhysIO/code/model/tapas_physio_create_hrv_regressors.m @@ -40,6 +40,7 @@ model_hrv = physio.model.hrv; end + delays = model_hrv.delays; if nargin < 4 @@ -49,6 +50,8 @@ slicenum = 1:sqpar.Nslices; + +% Calculate HR sample_points = tapas_physio_get_sample_points(ons_secs, sqpar, slicenum); hr = tapas_physio_hr(ons_secs.cpulse, sample_points); @@ -61,24 +64,23 @@ ylabel('beats per min (bpm)'); end -% create convolution for whole time series first... + +% Generate CRF dt = sqpar.TR/sqpar.Nslices; -t = 0:dt:32; % 32 seconds regressor +t = 0:dt:30; % seconds crf = tapas_physio_crf(t); -crf = crf/max(abs(crf)); -% crf = spm_hrf(dt); +crf = crf / max(abs(crf)); + if verbose.level>=2 subplot(2,2,2) plot(t, crf,'r');xlabel('time (seconds)'); title('Cardiac response function'); end -% NOTE: the removal of the mean was implemented to avoid over/undershoots -% at the 1st and last scans of the session due to convolution -convHRV = conv(hr-mean(hr), crf, 'same'); -% rescaling to -1...1 for display purposes -convHRV = convHRV./max(abs(convHRV)); +% Convolve and rescale for display purposes +convHRV = tapas_physio_conv(hr, crf, 'causal'); +convHRV = convHRV / max(abs(convHRV)); if verbose.level>=2 subplot(2,2,3) @@ -87,21 +89,20 @@ end -% create shifted regressors convolved time series, which is equivalent to -% delayed response functions according to Wikipedia (convoution) +% Create shifted regressors convolved time series, which is equivalent to +% delayed response functions according to Wikipedia (convolution) % % "Translation invariance[edit] % The convolution commutes with translations, meaning that % % \tau_x ({f}*g) = (\tau_x f)*g = {f}*(\tau_x g)\, -% where \tau_x fis the translation of the function f by x defined by -% (\tau_x f)(y) = f(y-x).\. +% where \tau_x is the translation of the function f by x defined by +% (\tau_x f)(y) = f(y-x). % remove mean and linear trend to fulfill periodicity condition for % shifting convHRV = detrend(convHRV); - % TODO: what happens at the end/beginning of shifted convolutions? nDelays = numel(delays); nShiftSamples = ceil(delays/dt); @@ -133,3 +134,5 @@ xlabel('time (seconds)');ylabel('regessor'); legend([hp{1}(1), hp{2}(1)], 'heart rate (bpm)', 'cardiac response regressor'); end + +end \ No newline at end of file diff --git a/PhysIO/code/model/tapas_physio_create_rvt_regressors.m b/PhysIO/code/model/tapas_physio_create_rvt_regressors.m index cc33150b..d97d4621 100644 --- a/PhysIO/code/model/tapas_physio_create_rvt_regressors.m +++ b/PhysIO/code/model/tapas_physio_create_rvt_regressors.m @@ -50,9 +50,11 @@ slicenum = 1:sqpar.Nslices; + +% Calculate RVT sample_points = tapas_physio_get_sample_points(ons_secs, sqpar, slicenum); rvt = tapas_physio_rvt(ons_secs.fr, ons_secs.t, sample_points, verbose); -rvt = rvt/max(rvt); % normalize for reasonable range of regressor +rvt = rvt / max(abs(rvt)); % normalize for reasonable range of regressor if verbose.level >=2 verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); @@ -63,26 +65,23 @@ ylabel('a.u.'); end -% create convolution for whole time series first... -dt = sqpar.TR/sqpar.Nslices; -t = 0:dt:50; % 50 seconds regressor + +% Generate RRF +dt = sqpar.TR / sqpar.Nslices; +t = 0:dt:60; % seconds rrf = tapas_physio_rrf(t); -rrf = rrf/max(abs(rrf)); -% crf = spm_hrf(dt); +rrf = rrf / max(abs(rrf)); if verbose.level >= 2 subplot(2,2,2) - plot(t, rrf,'g');xlabel('time (seconds)'); + plot(t, rrf,'g'); xlabel('time (seconds)'); title('Respiratory response function'); end -% NOTE: the removal of the mean was implemented to avoid over/undershoots -% at the 1st and last scans of the session due to convolution -convRVT = conv(rvt-mean(rvt), rrf, 'same'); - -% rescaling to -1...1 for display purposes -convRVT = convRVT./max(abs(convRVT)); +% Convolve and rescale for display purposes +convRVT = tapas_physio_conv(rvt, rrf, 'causal'); +convRVT = convRVT / max(abs(convRVT)); if verbose.level >= 2 subplot(2,2,3) @@ -90,21 +89,21 @@ title('Resp vol time X resp response function'); end -% create shifted regressors convolved time series, which is equivalent to -% delayed response functions according to Wikipedia (convoution) + +% Create shifted regressors convolved time series, which is equivalent to +% delayed response functions according to Wikipedia (convolution) % % "Translation invariance[edit] % The convolution commutes with translations, meaning that % % \tau_x ({f}*g) = (\tau_x f)*g = {f}*(\tau_x g)\, -% where \tau_x fis the translation of the function f by x defined by -% (\tau_x f)(y) = f(y-x).\. +% where \tau_x is the translation of the function f by x defined by +% (\tau_x f)(y) = f(y-x). % remove mean and linear trend to fulfill periodicity condition for % shifting convRVT = detrend(convRVT); - % TODO: what happens at the end/beginning of shifted convolutions? nDelays = numel(delays); nShiftSamples = ceil(delays/dt); @@ -137,3 +136,5 @@ legend([hp{1}(1), hp{2}(1)], 'respiratory volume / time (a. u.)', ... 'respiratory response regressor'); end + +end \ No newline at end of file diff --git a/PhysIO/code/model/tapas_physio_orthogonalise_physiological_regressors.m b/PhysIO/code/model/tapas_physio_orthogonalise_physiological_regressors.m index d0ff0b03..6d907073 100644 --- a/PhysIO/code/model/tapas_physio_orthogonalise_physiological_regressors.m +++ b/PhysIO/code/model/tapas_physio_orthogonalise_physiological_regressors.m @@ -25,8 +25,9 @@ % 'c' or 'cardiac' - only cardiac regressors are orthogonalised % 'r' or 'resp' - only respiration regressors are orthogonalised % 'mult' - only multiplicative regressors are orthogonalised -% 'all' - all physiological regressors are -% orthogonalised to each other +% 'RETROICOR' - cardiac, resp and interaction (mult) +% regressors are orthogonalised +% 'all' - all regressors are orthogonalised to each other % verbose.level 0 = no output; % 1 or other = plot design matrix before and after % orthogonalisation of physiological regressors and difference @@ -60,10 +61,16 @@ R = [cardiac_sess, tapas_physio_scaleorthmean_regressors(respire_sess) mult_sess input_R]; case {'mult'} R = [cardiac_sess, respire_sess, tapas_physio_scaleorthmean_regressors(mult_sess) input_R]; - case 'all' + case 'retroicor' R = [tapas_physio_scaleorthmean_regressors([cardiac_sess, respire_sess, mult_sess]), input_R]; + case 'all' + R = tapas_physio_scaleorthmean_regressors([cardiac_sess, respire_sess, mult_sess, input_R]); case {'n', 'none'} R = [cardiac_sess, respire_sess, mult_sess input_R]; + otherwise + verbose = tapas_physio_log(... + sprintf('Orthogonalisation of regressor set %s is not supported yet', ... + orthogonalise), verbose, 2) end diff --git a/PhysIO/code/model/tapas_physio_rvt.m b/PhysIO/code/model/tapas_physio_rvt.m index 1dd61759..0a3d9904 100644 --- a/PhysIO/code/model/tapas_physio_rvt.m +++ b/PhysIO/code/model/tapas_physio_rvt.m @@ -47,15 +47,12 @@ % COPYING or . - -dt = t(2)-t(1); -dtBreath = round(2/dt); %in seconds, minimum distance between two breaths - % compute breathing "pulses" (occurence times "rpulse" of max inhalation % times) -thresh_cardiac = []; -thresh_cardiac.min = .1; -thresh_cardiac.method = 'auto_matched'; +pulse_detect_options = []; +pulse_detect_options.min = .1; +pulse_detect_options.method = 'auto_matched'; +pulse_detect_options.max_heart_rate_bpm = 30;% actually the breathing rate breaths/per minute if nargin < 4 verbose.level = 0; @@ -65,9 +62,9 @@ verbose_no = verbose; verbose_no.level = 0; timeRpulseMax = tapas_physio_get_cardiac_pulses(t, fr, ... - thresh_cardiac,'OXY', dtBreath, verbose); + pulse_detect_options, 'OXY', verbose); timeRpulseMin = tapas_physio_get_cardiac_pulses(t, -fr, ... - thresh_cardiac,'OXY', dtBreath, verbose); + pulse_detect_options, 'OXY', verbose); nMax = numel(timeRpulseMax); nMin = numel(timeRpulseMin); maxFr = max(abs(fr)); @@ -79,9 +76,19 @@ ampRpulseMax = interp1(timeRpulseMax, fr(iTimeRpulseMax), t, 'linear', 'extrap'); ampRpulseMin = interp1(timeRpulseMin, fr(iTimeRpulseMin), t, 'linear', 'extrap'); +% Interpolate breath duration, but don't extrapolate durationBreath = diff(timeRpulseMax); -interpDurationBreath = interp1(timeRpulseMax(2:end), durationBreath,t, ... - 'linear', 'extrap'); +interpDurationBreath = interp1( ... + timeRpulseMax(2:end), durationBreath, t, ... + 'linear'); +% Nearest-neighbour interpolation for before/after last breath +% Be more careful here as can't let breath duration go negative +if sum(isnan(interpDurationBreath)) > 0 + nan_inds = isnan(interpDurationBreath); + interpDurationBreath(nan_inds) = interp1( ... + timeRpulseMax(2:end), durationBreath, t(nan_inds), ... + 'nearest', 'extrap'); +end if verbose.level>=2 verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); diff --git a/PhysIO/code/plot/tapas_physio_plot_spectrum.m b/PhysIO/code/plot/tapas_physio_plot_spectrum.m new file mode 100644 index 00000000..2301a998 --- /dev/null +++ b/PhysIO/code/plot/tapas_physio_plot_spectrum.m @@ -0,0 +1,62 @@ +function [fh, f, P1] = tapas_physio_plot_spectrum(t, y, hAx) +% Plots one-sided amplitude spectrum of a time series, following Matlab +% help examples +% +% [fh, f, P1] = tapas_physio_plot_spectrum(t, y) +% +% IN +% t [nSamples,1] time vector +% y [nSamples,nTimeseries] different signal time series, one per +% columns +% hAx axis handle for plot (optional, otherwise new figure is created) +% +% OUT +% fh figure handle +% f [nSamples/2+1,1] frequency vector (x axis) +% P1 [nSamples/2+1,1] one-sided amplitude spectrum +% +% EXAMPLE +% tapas_physio_plot_spectrum +% +% See also + +% Author: Lars Kasper +% Created: 2019-07-03 +% Copyright (C) 2019 TNU, Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of the TAPAS PhysIO Toolbox, which is released under +% the terms of the GNU General Public License (GPL), version 3. You can +% redistribute it and/or modify it under the terms of the GPL (either +% version 3 or, at your option, any later version). For further details, +% see the file COPYING or . + +if nargin < 3 + fh = tapas_physio_get_default_fig_params(); +else + axes(hAx); + fh = gcf; +end + +% freq vector for time +L=size(y,1); +Fs = 1/(t(2)-t(1)); +f = Fs*(0:(L/2))/L; + +for iCols = 1:size(y,2) + Y=fft(y(:,iCols)); + + % single sided spectrum + P2 = abs(Y/L); + P1 = P2(1:L/2+1); + P1(2:end-1) = 2*P1(2:end-1); + plot(f,P1); hold all +end +stringTitle = 'Preproc: Single-Sided Amplitude Spectrum'; +title(stringTitle); +xlabel('f (Hz)') +ylabel('|P1(f)|') + +if nargin < 3 % update name of figure only for newly created one + set(gcf, 'Name', stringTitle); +end \ No newline at end of file diff --git a/PhysIO/code/plot/tapas_physio_print_figs_to_file.m b/PhysIO/code/plot/tapas_physio_print_figs_to_file.m index b0b0db88..6fc60e94 100644 --- a/PhysIO/code/plot/tapas_physio_print_figs_to_file.m +++ b/PhysIO/code/plot/tapas_physio_print_figs_to_file.m @@ -1,7 +1,7 @@ function verbose = tapas_physio_print_figs_to_file(verbose, save_dir) % prints all figure handles in verbose-struct to specified filename there % -% physio_print_figs_to_ps(verbose) +% verbose = tapas_physio_print_figs_to_file(verbose, save_dir) % % IN % verbose.fig_handles @@ -10,7 +10,7 @@ % OUT % % EXAMPLE -% physio_print_figs_to_ps +% verbose = tapas_physio_print_figs_to_file(verbose, save_dir) % % See also @@ -20,10 +20,12 @@ % Created: 2013-04-23 % Copyright (C) 2013 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. % -% This file is part of the TNU CheckPhysRETROICOR toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . +% +% This file is part of the TAPAS PhysIO Toolbox, which is released under +% the terms of the GNU General Public Licence (GPL), version 3. You can +% redistribute it and/or modify it under the terms of the GPL (either +% version 3 or, at your option, any later version). For further details, +% see the file COPYING or . if nargin > 1 verbose.fig_output_file = fullfile(save_dir, verbose.fig_output_file); diff --git a/PhysIO/code/preproc/tapas_physio_filter_cardiac.m b/PhysIO/code/preproc/tapas_physio_filter_cardiac.m new file mode 100644 index 00000000..1b38b1f9 --- /dev/null +++ b/PhysIO/code/preproc/tapas_physio_filter_cardiac.m @@ -0,0 +1,111 @@ +function [fc, verbose] = tapas_physio_filter_cardiac(t,c, options, verbose) +% Bandpass-filters cardiac data using butterworth/chebyshev filters +% +% [fc, verbose] = tapas_physio_filter_cardiac(t,c,options, verbose) +% +% IN +% t [nSamples,1] time vector of cardiac time series +% c [nSamples, 1] cardiac time series +% options filter options +% type 'butter' Butterworth, flat passband) or +% 'cheby2' Chebychev Type II, for steep transitions +% between pass/stop band +% passband [f_min, f_max] in Hz +% stopband [f_min, f_max] in Hz +% verbose if .level >=2, output plot comparing filtered/unfiltered +% time series +% +% OUT +% fc filtered time series +% verbose if plots were created, figure handles added here +% +% EXAMPLE +% tapas_physio_filter_cardiac +% +% See also tapas_physio_new + +% Author: Lars Kasper, based on code snippets at +% https://ch.mathworks.com/matlabcentral/answers/327475-band-pass-butterworth-filter +% Created: 2019-07-02 +% Copyright (C) 2019 TNU, Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of the TAPAS PhysIO Toolbox, which is released under +% the terms of the GNU General Public License (GPL), version 3. You can +% redistribute it and/or modify it under the terms of the GPL (either +% version 3 or, at your option, any later version). For further details, +% see the file COPYING or . + +doFilter = options.include && ~isempty(options.passband); + +if ~doFilter + fc = c; +else + % set default stopband 10% bigger than passband + if isempty(options.stopband) + widthPassBand = diff(options.passband); + marginTransitionBand = 0.1*widthPassBand; + options.stopband = [ options.passband(1)-marginTransitionBand, ... + options.passband(2)+marginTransitionBand ]; + end + + Fsp = 1/(t(2) - t(1)); % sampling rate + deltaF = 1/(t(end)-t(1)); + + % filter boundaries check! + options.passband(1) = max(2*deltaF, options.passband(1)); + options.passband(2) = min(Fsp-deltaF, options.passband(2)); + options.stopband(1) = max(deltaF, options.stopband(1)); + options.stopband(2) = min(Fsp, options.stopband(2)); + + + Fn = Fsp/2; + Wp = options.passband/Fn; + Ws = options.stopband/Fn; + + switch lower(options.type) + case {'butter', 'butterworth'} + + [z,p,k]=butter(8,Wp,'bandpass'); + [sos,g]=zp2sos(z,p,k); + + % Rp=1; + % Rs=25; + % [n,Wn] = buttord(Wp,Ws,Rp,Rs); + % [b,a]=butter(n,Wn); + % [sos,g]=tf2sos(b,a); + case {'cheby2', 'chebychev'} + Rp=10; + Rs=30; + [n,Ws] = cheb2ord(Wp,Ws,Rp,Rs); + [z,p,k] = cheby2(n,Rs,Ws); + [sos,g] = zp2sos(z,p,k); + end + fc=filtfilt(sos,g,c); + + + %% plot filter response + if verbose.level >=3 + fvtool(sos,'Analysis','freq') + end + + %% plot filtering results + if verbose.level >=2 + stringTitle = 'Preproc: Bandpass-filtered Filtered Cardiac time series'; + verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params(); + set(gcf, 'Name', stringTitle); + + % plot raw/filtered time series + subplot(2,1,1); + plot(t, c); hold all; + plot(t, fc); + xlabel('t(s)'); + ylabel('Cardiac Wave Amplitude (a.u.)'); + legend('raw', sprintf('filtered (%s)', options.type)); + title(stringTitle); + + subplot(2,1,2); + tapas_physio_plot_spectrum(t,[c fc], gca); + legend('raw', sprintf('filtered (%s)', options.type)); + end +end diff --git a/PhysIO/code/preproc/tapas_physio_find_ecg_r_peaks.m b/PhysIO/code/preproc/tapas_physio_find_ecg_r_peaks.m index c9c3c8a6..9a66c970 100644 --- a/PhysIO/code/preproc/tapas_physio_find_ecg_r_peaks.m +++ b/PhysIO/code/preproc/tapas_physio_find_ecg_r_peaks.m @@ -86,7 +86,16 @@ %% smooth ECG curve with R-wave kernel and plot autocorrelation -sy = conv(y./sqrt(sum(kRpeak.^2)),kRpeak/sqrt(sum(kRpeak.^2)),'same'); +% https://en.wikipedia.org/wiki/Matched_filter +filter = kRpeak / sqrt(sum(kRpeak.^2)); +if mod(length(filter), 2) == 0 + filter = [filter(:); 0.0]; % tapas_physio_conv needs odd length +end +sy = tapas_physio_conv(y / sqrt(sum(kRpeak.^2)), flip(filter), 'symmetric'); +% Note we don't necessarily know the phase here! We assume a symmetric +% filter (i.e. the central point is `t=0`), but that isn't necessarily the +% case (e.g. for a manual template). However, that gives us the detected +% peaks at the centre of the template, which is typically what we want. peaks_found = false; thresh_changed = false; diff --git a/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m b/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m index 463e7afe..6d71ee61 100644 --- a/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m +++ b/PhysIO/code/preproc/tapas_physio_get_cardiac_phase.m @@ -57,14 +57,14 @@ % add pulses in the end, if last onset time of scan after last recorded % pulses - use guess based on average heart rate if scannert(end) > pulset(end) - verbose = tapas_physio_log(... - 'Guessed additional cardiac pulse at time series end for phase estimation', ... - verbose, 1); % add more pulses before first one, using same average heartbeat % duration as in first nAverage cycles meanCycleDur = mean(diff(pulset((end-nAverage+1):end))); nAddPulses = ceil((scannert(end) - pulset(end))/meanCycleDur); pulset = [pulset; pulset(end) + meanCycleDur*(1:nAddPulses)']; + verbose = tapas_physio_log(... + sprintf('Note: Guessed %d additional cardiac pulse(s) at time series end for phase estimation', nAddPulses), ... + verbose, 0); end scannertpriorpulse = zeros(1,length(scannert)); diff --git a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m index 93ebd619..5cff7b5a 100644 --- a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m +++ b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulse_template.m @@ -118,7 +118,7 @@ plot(t, c, 'k'); title('Preproc: Finding first peak of cycle, backwards') end verbose = tapas_physio_log(['No peaks found in raw cardiac time series. Check raw ' ... - 'physiological recordings figure whether there is any non-constant' ... + 'physiological recordings figure whether there is any non-constant ' ... 'cardiac data'], verbose, 2); % error! end diff --git a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses.m b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses.m index 8475dc78..e5920ea0 100644 --- a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses.m +++ b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses.m @@ -1,13 +1,14 @@ function [cpulse, verbose] = tapas_physio_get_cardiac_pulses(t, c, ... - thresh_cardiac, cardiac_modality, dt120, verbose) + cpulse_detect_options, cardiac_modality, verbose) % extract heartbeat events from ECG or pulse oximetry time course % -% cpulse = tapas_physio_get_cardiac_pulses(t, c, thresh_cardiac, cardiac_modality, verbose); +% cpulse = tapas_physio_get_cardiac_pulses(t, c, cpulse_detect_options, cardiac_modality, verbose); % % IN % t vector of time series of log file (in seconds, corresponds to c) % c raw time series of ECG or pulse oximeter -% thresh_cardiac is a structure with the following elements +% cpulse_detect_options +% is a structure with the following elements % .method - 'auto_matched', 'manual_template', 'load_from_logfile', % 'load_template' % Specifies how to determine QRS-wave from noisy input @@ -28,14 +29,15 @@ % This file is saved after picking the QRS-wave % manually (i.e. if .ECG_min is set), so that % results are reproducible +% .max_heart_rate_bpm +% maximum allowed physiological heart rate (in beats +% per minute) for subject; default: 90 bpm % .krPeak [false] or true; if true, a user input is -% required to specify a characteristic R-peak interval in the ECG -% or pulse oximetry time series +% required to specify a characteristic R-peak interval in the ECG +% or pulse oximetry time series % cardiac_modality 'ECG', 'ECG_WiFi' electrocardiogram (with/without % wireless transmisssion for Philips data) % 'OXY'/'PPU' pulse oximetry unit -% dt120 - minimum distance between heart beats; default 120 -% bpm, i.e. 0.5 s % % verbose Substructure of Physio, holding verbose.level and % verbose.fig_handles with plotted figure handles @@ -64,16 +66,20 @@ %% detection of cardiac R-peaks -dt = t(2)-t(1); -if nargin < 5 || isempty(dt120) - dt120 = round(0.5/dt); % heart rate < 120 bpm +dt = t(2) - t(1); +minPulseDistanceSamples = ... + floor((1 / (cpulse_detect_options.max_heart_rate_bpm / 60)) / dt); + +if isempty(minPulseDistanceSamples) + minPulseDistanceSamples = round(0.5 / dt); % heart rate < 120 bpm end + switch lower(cardiac_modality) case 'oxy_old' [cpulse, verbose] = tapas_physio_get_oxy_pulses_filtered(c, t, ... - dt120, verbose); + minPulseDistanceSamples, verbose); otherwise % {'oxy','ppu', 'oxy_wifi', 'ppu_wifi','ecg', 'ecg_wifi'} etc., including ecg_raw o - switch thresh_cardiac.method + switch cpulse_detect_options.method case 'load_from_logfile' verbose = tapas_physio_log('How did you end up here? I better do nothing.', ... verbose, 1); @@ -81,10 +87,10 @@ case {'manual', 'manual_template', 'load', 'load_template'} % load/determine manual template [cpulse, verbose] = ... tapas_physio_get_cardiac_pulses_manual_template(... - c, t, thresh_cardiac, verbose); + c, t, cpulse_detect_options, verbose); case {'auto', 'auto_template', 'auto_matched'} [cpulse, verbose] = ... tapas_physio_get_cardiac_pulses_auto_matched( ... - c, t, thresh_cardiac.min, dt120, verbose); - end % switch thresh_cardiac.method + c, t, cpulse_detect_options.min, minPulseDistanceSamples, verbose); + end % switch cpulse_detect_options.method end diff --git a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_manual_template.m b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_manual_template.m index 38f59092..3b2e8bbe 100644 --- a/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_manual_template.m +++ b/PhysIO/code/preproc/tapas_physio_get_cardiac_pulses_manual_template.m @@ -1,5 +1,5 @@ function [cpulse, verbose] = tapas_physio_get_cardiac_pulses_manual_template(... - c, t, thresh_cardiac_initial_cpulse_select, verbose) + c, t, pulse_detect_options, verbose) % Detects R-peaks via matched-filter smoothing & peak detection using a % manually defined QRS-wave (or R-peak environment) % @@ -9,7 +9,7 @@ % IN % c [nSamples, 1] raw pulse oximeter samples % t [nSamples, 1] time vector corresponding to samples (un seconds) -% thresh_cardiac_initial_cpulse_select +% pulse_detect_options % physio.thresh.cardiac.initial_cpulse_select-substructure % with elements % .method 'manual' or 'load'/'load_template' @@ -45,30 +45,30 @@ % manual peak selection, if no file selected and loading is % specified -hasKrpeakLogfile = exist(thresh_cardiac_initial_cpulse_select.file,'file') || ... - exist([thresh_cardiac_initial_cpulse_select.file '.mat'],'file'); +hasKrpeakLogfile = exist(pulse_detect_options.file,'file') || ... + exist([pulse_detect_options.file '.mat'],'file'); % if no file exists, also do manual peak-find doSelectTemplateManually = any(strcmpi(... - thresh_cardiac_initial_cpulse_select.method, ... + pulse_detect_options.method, ... {'manual', 'manual_template'})) || ~hasKrpeakLogfile; if doSelectTemplateManually - thresh_cardiac_initial_cpulse_select.kRpeak = []; - hasECGMin = isfield(thresh_cardiac_initial_cpulse_select, 'min') && ~isempty(thresh_cardiac_initial_cpulse_select.min); + pulse_detect_options.kRpeak = []; + hasECGMin = isfield(pulse_detect_options, 'min') && ~isempty(pulse_detect_options.min); if ~hasECGMin - thresh_cardiac_initial_cpulse_select.min = 0.5; + pulse_detect_options.min = 0.5; end else - fprintf('Loading %s\n', thresh_cardiac_initial_cpulse_select.file); - ECGfile = load(thresh_cardiac_initial_cpulse_select.file); - thresh_cardiac_initial_cpulse_select.min = ECGfile.ECG_min; - thresh_cardiac_initial_cpulse_select.kRpeak = ECGfile.kRpeak; + fprintf('Loading %s\n', pulse_detect_options.file); + ECGfile = load(pulse_detect_options.file); + pulse_detect_options.min = ECGfile.ECG_min; + pulse_detect_options.kRpeak = ECGfile.kRpeak; end inp_events = []; -ECG_min = thresh_cardiac_initial_cpulse_select.min; -kRpeak = thresh_cardiac_initial_cpulse_select.kRpeak; +ECG_min = pulse_detect_options.min; +kRpeak = pulse_detect_options.kRpeak; if doSelectTemplateManually while ECG_min [cpulse, ECG_min_new, kRpeak] = tapas_physio_find_ecg_r_peaks(t,c, ECG_min, [], inp_events); @@ -83,5 +83,5 @@ % save manually found peak parameters to file if doSelectTemplateManually - save(thresh_cardiac_initial_cpulse_select.file, 'ECG_min', 'kRpeak'); + save(pulse_detect_options.file, 'ECG_min', 'kRpeak'); end \ No newline at end of file diff --git a/PhysIO/code/preproc/tapas_physio_get_oxy_pulses_filtered.m b/PhysIO/code/preproc/tapas_physio_get_oxy_pulses_filtered.m index 336c13c3..0c207dc0 100644 --- a/PhysIO/code/preproc/tapas_physio_get_oxy_pulses_filtered.m +++ b/PhysIO/code/preproc/tapas_physio_get_oxy_pulses_filtered.m @@ -40,8 +40,8 @@ c = c-mean(c); c = c./max(c); % normalize time series % smooth noisy pulse oximetry data to detect peaks -w = tapas_physio_gausswin(dt120,1); -sc = conv(c, w, 'same'); +w = tapas_physio_gausswin(2*floor(dt120/2)+1, 1); % Odd number of samples +sc = tapas_physio_conv(c, w, 'symmetric'); sc = sc-mean(sc); sc = sc./max(sc); % normalize time series % Highpass filter to remove drifts diff --git a/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m b/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m index 75fed38b..c0019dd1 100644 --- a/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m +++ b/PhysIO/code/preproc/tapas_physio_get_respiratory_phase.m @@ -75,9 +75,12 @@ % Calculate derivative of normalised pulse wrt time % over 1 sec of data as described in Glover et al. ksize = round(0.5 * (1/rsampint)); -kernel = [ones(1, ksize)*-1 0 ones(1, ksize)]; -dpulse = -conv(pulset, kernel); -dpulse = dpulse(ksize+1:end-ksize); +kernel = ksize:-1:-ksize; kernel = kernel ./ sum(kernel.^2); +dpulse = tapas_physio_conv(pulset, kernel, 'symmetric'); +% This uses a quadratic Savitzky-Golay filter, for which the coefficients +% have a simple linear form. See e.g. +% https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter +% `scipy.signal.savgol_coeffs(5, polyorder=2, deriv=1, use='conv')` % Tolerance to the derivative dpulse(abs(dpulse) < 1e-4) = 0; diff --git a/PhysIO/code/preproc/tapas_physio_simulate_pulse_samples.m b/PhysIO/code/preproc/tapas_physio_simulate_pulse_samples.m index 32f5e9da..90d7d127 100644 --- a/PhysIO/code/preproc/tapas_physio_simulate_pulse_samples.m +++ b/PhysIO/code/preproc/tapas_physio_simulate_pulse_samples.m @@ -76,6 +76,7 @@ end % ...and convolve with pulse template +% **TODO** tapas_physio_conv simulatedSamples = conv(simulatedPulses, pulseTemplate, 'same'); if doDebug diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m index be15fa9a..051214a7 100644 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m @@ -79,6 +79,7 @@ hasRespirationFile = ~isempty(log_files.respiration); hasCardiacFile = ~isempty(log_files.cardiac); +hasExplicitJsonFile = ~isempty(log_files.scan_timing); if hasCardiacFile fileName = log_files.cardiac; @@ -100,7 +101,9 @@ fileJson = regexprep(fileName, '\.tsv', '\.json'); end - +if hasExplicitJsonFile + fileJson = log_files.scan_timing; +end hasJsonFile = isfile(fileJson); @@ -123,14 +126,23 @@ end end -tRelStartScan = log_files.relative_start_acquisition; -if isempty(tRelStartScan) +% sum implicit (.json) and explicit relative shifts of log/scan acquisition +if isempty(log_files.relative_start_acquisition) if hasJsonFile % in BIDS, start of the phys logging is stated relative to the first volume scan start. % PhysIO defines the scan acquisiton relative to the phys log start tRelStartScan = -val.StartTime; else - tRelStartScan = 0; + verbose = tapas_physio_log(... + ['No .json file found and empty log_files.relative_start_acquisition. ' ... + 'Please specify explicitly.'], verbose, 2); + end +else + if hasJsonFile + % add both delays + tRelStartScan = log_files.relative_start_acquisition - val.StartTime; + else + tRelStartScan = log_files.relative_start_acquisition; end end diff --git a/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m b/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m index 0843e141..c62661a6 100755 --- a/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m +++ b/PhysIO/code/readin/tapas_physio_read_physlogfiles_siemens.m @@ -76,6 +76,12 @@ dt = log_files.sampling_interval; +explicit_relative_start_acquisition = log_files.relative_start_acquisition; + +if isempty(explicit_relative_start_acquisition) + explicit_relative_start_acquisition = 0; +end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Determine relative start of acquisition from dicom headers and @@ -133,7 +139,7 @@ % add arbitrary offset specified by user relative_start_acquisition = relative_start_acquisition + ... - log_files.relative_start_acquisition; + explicit_relative_start_acquisition; data_table = tapas_physio_siemens_line2table(lineData, cardiac_modality); @@ -204,7 +210,7 @@ % add arbitrary offset specified by user relative_start_acquisition = relative_start_acquisition + ... - log_files.relative_start_acquisition; + explicit_relative_start_acquisition; data_table = tapas_physio_siemens_line2table(lineData, 'RESP'); diff --git a/PhysIO/code/sync/tapas_physio_create_scan_timing_from_gradients_auto_philips.m b/PhysIO/code/sync/tapas_physio_create_scan_timing_from_gradients_auto_philips.m index 2e71d748..6c7f583c 100644 --- a/PhysIO/code/sync/tapas_physio_create_scan_timing_from_gradients_auto_philips.m +++ b/PhysIO/code/sync/tapas_physio_create_scan_timing_from_gradients_auto_philips.m @@ -366,6 +366,7 @@ set(gcf, 'Name', stringTitle); ampl = max(abs(G)); +% **TODO** tapas_physio_conv sG = conv(G, templateGradientVolume/sum(abs(templateGradientVolume)), ... 'same'); @@ -396,6 +397,7 @@ set(gcf, 'Name', stringTitle); ampl = max(abs(G)); +% **TODO** tapas_physio_conv sG = conv(G, templateGradientSlice/sum(abs(templateGradientSlice)), ... 'same'); diff --git a/PhysIO/code/sync/tapas_physio_rescale_gradient_gain_fluctuations.m b/PhysIO/code/sync/tapas_physio_rescale_gradient_gain_fluctuations.m index 51efbd13..239bc496 100644 --- a/PhysIO/code/sync/tapas_physio_rescale_gradient_gain_fluctuations.m +++ b/PhysIO/code/sync/tapas_physio_rescale_gradient_gain_fluctuations.m @@ -59,6 +59,7 @@ % determine positive and negative steps warning off tapas_physio_findpeaks:largeMinPeakHeight +warning off signal:findpeaks:largeMinPeakHeight [~, idxGainPlus] = tapas_physio_findpeaks((dmG), 'minpeakDistance', n, ... 'minpeakheight', minPeakHeight); @@ -66,6 +67,7 @@ 'minpeakheight', minPeakHeight); warning on tapas_physio_findpeaks:largeMinPeakHeight +warning on signal:findpeaks:largeMinPeakHeight % plus gains refer to max-changes in the future idxGainPlus = idxGainPlus + n; diff --git a/PhysIO/code/tapas_physio_cfg_matlabbatch.m b/PhysIO/code/tapas_physio_cfg_matlabbatch.m index 90bb87de..e24ff2d2 100644 --- a/PhysIO/code/tapas_physio_cfg_matlabbatch.m +++ b/PhysIO/code/tapas_physio_cfg_matlabbatch.m @@ -167,7 +167,8 @@ relative_start_acquisition.help = { ' Time (in seconds) when the 1st scan (or, if existing, dummy) started,' ' relative to the start of the logfile recording;' - ' e.g. 0 if simultaneous start' + ' [] (empty) to read from explicit acquisition timing info (s.b.)' + ' 0 if simultaneous start' ' 10, if 1st scan starts 10' ' seconds AFTER physiological' ' recording' @@ -176,12 +177,12 @@ ' NOTE: ' ' 1. For Philips SCANPHYSLOG, this parameter is ignored, if' ' scan_timing.sync is set.' - ' 2. If you specify an acquisition_info file, leave this parameter' - ' at 0 (e.g., for Siemens_Tics) since physiological recordings' - ' and acquisition timing are already synchronized by this' - ' information, and you would introduce another shift.' - ' 3. For BIDS, relative_start_acquisition is read as -StartTime from' - ' accompanying json-file, if existing' + ' 2. If you specify an acquisition_info file, leave this' + ' parameter empty or 0 (e.g., for Siemens_Tics, BIDS) since' + ' physiological recordings and acquisition timing are already' + ' synchronized by this information, and you would introduce an' + ' additional shift.' + }; relative_start_acquisition.strtype = 'e'; relative_start_acquisition.num = [Inf Inf]; @@ -555,6 +556,27 @@ initial_cpulse_select_file.num = [0 Inf]; initial_cpulse_select_file.val = {'initial_cpulse_kRpeakfile.mat'}; +%-------------------------------------------------------------------------- +% max_heart_rate_bpm +%-------------------------------------------------------------------------- +max_heart_rate_bpm = cfg_entry; +max_heart_rate_bpm.tag = 'max_heart_rate_bpm'; +max_heart_rate_bpm.name = 'Maximum heart rate (BPM)'; +max_heart_rate_bpm.help = { + 'Maximum expected heart rate in beats per minute. (default: 90)' + 'This only needs to be a rough guess and should be changed for specific' + 'subject populations.' + ' - If set too low, the auto_matched pulse detection might miss genuine' + ' cardiac pulses' + ' - If set too high, it might introduce artifactual pulse events, i.e.' + ' interpreting local maxima within a pulse as new pulse events' + ' You may need to increase this value if you have a subject with a very' + ' high heart rate, or decrease it if you have very pronounced local maxima' + ' in your wave form.' + }; +max_heart_rate_bpm.strtype = 'e'; +max_heart_rate_bpm.num = [0 Inf]; +max_heart_rate_bpm.val = {90}; %-------------------------------------------------------------------------- @@ -578,7 +600,8 @@ initial_cpulse_select_method_auto_template.name = 'auto_template'; initial_cpulse_select_method_auto_template.val = { min - initial_cpulse_select_file + initial_cpulse_select_file + max_heart_rate_bpm }; initial_cpulse_select_method_auto_template.help = { ... ' Auto generation of representative QRS-wave; detection via' @@ -595,6 +618,7 @@ initial_cpulse_select_method_auto_matched.val = { min initial_cpulse_select_file + max_heart_rate_bpm }; initial_cpulse_select_method_auto_matched.help = { ... 'Auto generation of template QRS wave, ' @@ -794,6 +818,10 @@ }; +%-------------------------------------------------------------------------- +% filter for cardiac time series +%-------------------------------------------------------------------------- +filter = tapas_physio_gui_filter(); %-------------------------------------------------------------------------- % cardiac @@ -801,7 +829,7 @@ cardiac = cfg_branch; cardiac.tag = 'cardiac'; cardiac.name = 'cardiac'; -cardiac.val = {modality initial_cpulse_select posthoc_cpulse_select}; +cardiac.val = {modality filter initial_cpulse_select posthoc_cpulse_select}; cardiac.help = {'...'}; diff --git a/PhysIO/code/tapas_physio_init.m b/PhysIO/code/tapas_physio_init.m index 0e0a990b..2544ea20 100644 --- a/PhysIO/code/tapas_physio_init.m +++ b/PhysIO/code/tapas_physio_init.m @@ -23,6 +23,12 @@ % COPYING or . % +% add path for utils, if physio not in path +if ~exist('tapas_physio_logo') + pathPhysio = fileparts(mfilename('fullpath')); + addpath(fullfile(pathPhysio, 'utils')); % needed for further path checks +end + tapas_physio_logo(); % print logo disp('Checking Matlab PhysIO paths, SPM paths and Batch Editor integration now...'); @@ -32,8 +38,8 @@ [isPhysioOnPath, pathPhysIO] = tapas_physio_check_path(); if ~isPhysioOnPath - addpath(pathPhysIO); - fprintf('added PhysIO path: %s\n', pathPhysIO); + addpath(genpath(pathPhysIO)); + fprintf('added PhysIO path recursively: %s\n', pathPhysIO); else fprintf('OK.\n'); end diff --git a/PhysIO/code/tapas_physio_main_create_regressors.m b/PhysIO/code/tapas_physio_main_create_regressors.m index 35750e11..bd235844 100644 --- a/PhysIO/code/tapas_physio_main_create_regressors.m +++ b/PhysIO/code/tapas_physio_main_create_regressors.m @@ -39,12 +39,14 @@ % Author: Lars Kasper % Created: 2011-08-01 -% Copyright (C) 2011-2018 TNU, Institute for Biomedical Engineering, University of Zurich and ETH Zurich. +% Copyright (C) 2011-2019 TNU, Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. % -% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public -% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL -% (either version 3 or, at your option, any later version). For further details, see the file -% COPYING or . +% This file is part of the TAPAS PhysIO Toolbox, which is released under +% the terms of the GNU General Public Licence (GPL), version 3. You can +% redistribute it and/or modify it under the terms of the GPL (either +% version 3 or, at your option, any later version). For further details, +% see the file COPYING or . %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -57,7 +59,6 @@ % These parameters could become toolbox inputs... minConstantIntervalAlertSeconds = 0.2; -maxHeartRateBpm = 90; if ~nargin error('Please specify a PhysIO-object as input to this function. See tapas_physio_new'); @@ -158,17 +159,19 @@ % preproc.cardiac.modality = 'OXY'; % 'ECG' or 'OXY' (for pulse oximetry) %% initial pulse select via load from logfile or autocorrelation with 1 %% cardiac pulse + + [ons_secs.c, verbose] = tapas_physio_filter_cardiac(... + ons_secs.t, ons_secs.c, preproc.cardiac.filter, verbose); + switch preproc.cardiac.initial_cpulse_select.method case {'load_from_logfile', ''} % do nothing otherwise - minCardiacCycleSamples = floor((1/(maxHeartRateBpm/60)/ons_secs.dt)); - % run one of the various cardiac pulse detection algorithms [ons_secs.cpulse, verbose] = ... tapas_physio_get_cardiac_pulses(ons_secs.t, ons_secs.c, ... preproc.cardiac.initial_cpulse_select, ... - preproc.cardiac.modality, minCardiacCycleSamples, verbose); + preproc.cardiac.modality, verbose); end diff --git a/PhysIO/code/tapas_physio_new.m b/PhysIO/code/tapas_physio_new.m index e1756f19..c6ee6036 100644 --- a/PhysIO/code/tapas_physio_new.m +++ b/PhysIO/code/tapas_physio_new.m @@ -166,7 +166,9 @@ % Time (in seconds) when the 1st scan (or, if existing, dummy) started, % relative to the start of the logfile recording; - % e.g. 0 if simultaneous start + % e.g. + % [] (empty) to read from explicit acquisition timing info (s.b.) + % 0 if simultaneous start % 10, if 1st scan starts 10 % seconds AFTER physiological % recording @@ -175,10 +177,11 @@ % NOTE: % 1. For Philips SCANPHYSLOG, this parameter is ignored, if % scan_timing.sync is set. - % 2. If you specify an acquisition_info file, leave this parameter - % at 0 (e.g., for Siemens_Tics) since physiological recordings - % and acquisition timing are already synchronized by this - % information, and you would introduce another shift. + % 2. If you specify an acquisition_info file, leave this + % parameter empty or 0 (e.g., for Siemens_Tics, BIDS) since + % physiological recordings and acquisition timing are already + % synchronized by this information, and you would introduce an + % additional shift. % log_files.relative_start_acquisition = 0; @@ -291,11 +294,45 @@ preproc = []; preproc.cardiac = []; - preproc.cardiac.modality = 'ecg_wifi'; % 'ECG','ECG_raw', or 'OXY'/'PPU' (for pulse oximetry), 'OXY_OLD', [deprecated] + + % Measurement modality of input cardiac signal + % 'ECG','ECG_raw', or 'OXY'/'PPU' (for pulse oximetry), 'OXY_OLD', [deprecated] + preproc.cardiac.modality = 'ecg_wifi'; + + % Filter properties for bandpass-filtering of cardiac signal before peak + % detection, phase extraction, and other physiological traces + preproc.cardiac.filter = []; + + preproc.cardiac.filter.include = 0; % 1 = filter executed; 0 = not used + + % filter type default: 'cheby2' + % 'cheby2' Chebychev Type II filter, use for steep transition from + % start to stop band + % 'butter' butterworth filter, standard filter with maximally flat + % passband (Infinite impulse response), but stronger + % ripples in transition band + preproc.cardiac.filter.type = 'butter'; + + % + % [f_min, f_max] frequency interval in Hz of all frequency that should + % pass the passband filter + % default: [0.3 9] (to remove slow drifts, breathing + % and slice sampling artifacts) + % if empty, no filtering is performed + preproc.cardiac.filter.passband = [0.3 9]; + + % [f_min, f_max] frequency interval in Hz of all frequencies, s.th. frequencies + % outside this band should definitely *NOT* pass the filter + % Default: [] + % NOTE: only relevant for 'cheby2' filter type + % if empty, and passband is empty, no filtering is performed + % if empty, but passband exists, stopband interval is + % 10% increased passband interval + preproc.cardiac.filter.stopband = []; % The initial cardiac pulse selection structure: Determines how the % majority of cardiac pulses is detected - % default: auto + % default: 'auto_matched' % % 'auto_matched' % - auto generation of representative QRS-wave; detection via @@ -306,6 +343,17 @@ % 'load' - from previous manual/auto run preproc.cardiac.initial_cpulse_select.method = 'auto_matched'; + % maximum allowed physiological heart rate (in beats per minute) + % for subject; default: 90 bpm + % - If set too low, the auto_mathed pulse detection might miss genuine + % cardiac pulses + % - If set too high, it might introduce artifactual pulse events, i.e. + % interpreting local maxima within a pulse as new pulse events + % Adjust this value, if you have a subject with very high heart rate + % (increase!), or if you have very pronounced local maxima in your wave form + % (decrease!). + preproc.cardiac.initial_cpulse_select.max_heart_rate_bpm = 90; + % file containing reference ECG-peak (variable kRpeak) % used for method 'manual' or 'load' [default: not set] % if method == 'manual', this file is saved after picking the QRS-wave diff --git a/PhysIO/code/tapas_physio_version.m b/PhysIO/code/tapas_physio_version.m index 9b01b8cf..f9b49671 100644 --- a/PhysIO/code/tapas_physio_version.m +++ b/PhysIO/code/tapas_physio_version.m @@ -1,7 +1,7 @@ function versionPhysio = tapas_physio_version() % returns version number of this copy of PhysIO % -% output = tapas_physio_version(input) +% versionPhysio = tapas_physio_version() % % IN % @@ -23,4 +23,4 @@ % version 3 or, at your option, any later version). For further details, % see the file COPYING or . % -versionPhysio = 'R2019a-v7.1.0'; \ No newline at end of file +versionPhysio = 'R2019b-v7.2.0'; \ No newline at end of file diff --git a/PhysIO/code/utils/tapas_physio_conv.m b/PhysIO/code/utils/tapas_physio_conv.m new file mode 100644 index 00000000..b6a65c41 --- /dev/null +++ b/PhysIO/code/utils/tapas_physio_conv.m @@ -0,0 +1,71 @@ +function [w] = tapas_physio_conv(u, v, filter_type, padding) +% Wrapper around `conv()` for convolution +% Deals with the padding and time offsets +% +% [w] = tapas_physio_filter_cardiac(u, v, filter_type) +% +% IN +% u Data time series [1-D] +% v Convolutional filter time series [1-D] +% filter_type ['causal', 'symmetric']. +% If 'causal', the filter is taken to be defined for +% `t >= 0`, with the first element corresponding to `t=0`. +% If 'symmetric', the filter is taken to be defined for both +% positive and negative time, with the central element +% corresponding to `t=0`. +% N.B. 'symmetric' implies an *odd* filter length. +% padding ['mean', 'zero']. Whether to pad with the mean of `u`, or +% with zeros. +% +% OUT +% w Convolved time series [size(u)] + +% Author: Sam Harrison +% Created: 2019-07-11 +% Copyright (C) 2019 TNU, Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of the TAPAS PhysIO Toolbox, which is released under +% the terms of the GNU General Public License (GPL), version 3. You can +% redistribute it and/or modify it under the terms of the GPL (either +% version 3 or, at your option, any later version). For further details, +% see the file COPYING or . + +if ~isvector(u) || ~isvector(v) + error('tapas_physio_conv: Both inputs must be vectors') +end +% Padding: use data mean to reduce transients by default +if nargin < 4 + padding = 'mean'; +end + +% Pad value +switch lower(padding) + case 'mean' + pad_val = mean(u); + case 'zero' + pad_val = 0.0; + otherwise + error('Unrecognised padding argument (%s)', padding) +end + +% Pad shape +switch lower(filter_type) + case 'causal' + u_pad = [pad_val * ones(length(v)-1, 1); u(:)]; + case 'symmetric' + if mod(length(v), 2) == 1 + pad_elem = (length(v) - 1) / 2; + u_pad = [pad_val * ones(pad_elem, 1); u(:); pad_val * ones(pad_elem, 1)]; + else + error('Symmetric filter lengths must be odd!') + end + otherwise + error('Unrecognised padding argument (%s)', filter_type) +end + +% Apply convolution and select portion of interest +w = conv(u_pad, v, 'valid'); +w = reshape(w, size(u)); + +end \ No newline at end of file diff --git a/PhysIO/code/utils/tapas_physio_job2physio.m b/PhysIO/code/utils/tapas_physio_job2physio.m index 0a696e86..bdc6a17f 100644 --- a/PhysIO/code/utils/tapas_physio_job2physio.m +++ b/PhysIO/code/utils/tapas_physio_job2physio.m @@ -6,11 +6,13 @@ % IN % % OUT +% physio physio input structure, as use by +% tapas_physio_main_create_regressors % % EXAMPLE % physio = tapas_physio_job2physio(job) % -% See also spm_physio_cfg_matlabbatch +% See also tapas_physio_cfg_matlabbatch tapas_physio_main_create_regressors % Author: Lars Kasper % Created: 2015-01-05 @@ -26,7 +28,7 @@ physio = tapas_physio_new(); -%% Use existing properties that are cfg_choices in job to overwrite +%% Use existing properties that are cfg_choices in job to overwrite % properties of physio and set corresponding method physio = tapas_physio_update_from_job(physio, job, ... @@ -43,8 +45,8 @@ 'noise_rois', 'other'}; physio = tapas_physio_update_from_job(physio, job, ... - strcat('model.', modelArray), strcat('model.', modelArray), ... - true, 'include'); + strcat('model.', modelArray), strcat('model.', modelArray), ... + true, 'include'); %% Convert yes => true (=1) and no => false (=0) nModels = numel(modelArray); @@ -53,6 +55,26 @@ physio.model.(modelArray{iModel}).include, 'yes'); end +%% Take over yes/no substructs as is, yes/no will become 'include' property +yesNoArray = ... + {'preproc.cardiac.filter'}; + +physio = tapas_physio_update_from_job(physio, job, ... + yesNoArray, yesNoArray, true, 'include'); + +%% Convert yes => true (=1) and no => false (=0) +nChoices = numel(yesNoArray); +for iChoice = 1:nChoices + try + eval(sprintf(['physio.%s.include = strcmpi(' ... + 'physio.%s.include, ''yes'');'], yesNoArray{iChoice}, ... + yesNoArray{iChoice})); + catch err + tapas_physio_log(sprintf('No property %s defined in job (error: %s)', ... + yesNoArray{iChoice}, err.message), [], 1); + end +end + %% Use existing properties in job to overwrite properties of physio physio = tapas_physio_update_from_job(physio, job, ... {'preproc.cardiac.modality', 'scan_timing.sqpar', ... diff --git a/PhysIO/code/utils/tapas_physio_save_batch_mat_script.m b/PhysIO/code/utils/tapas_physio_save_batch_mat_script.m index 8471c72a..52f2a93d 100644 --- a/PhysIO/code/utils/tapas_physio_save_batch_mat_script.m +++ b/PhysIO/code/utils/tapas_physio_save_batch_mat_script.m @@ -48,6 +48,7 @@ if ~ischar(fileBatchM) matlabbatch = fileBatchM; fileBatchM = fullfile(pwd, 'physio_job.m'); + fileBatchMat = fullfile(pwd, 'physio_job.mat'); else [fp,fn,ext] = fileparts(fileBatchM); switch ext diff --git a/PhysIO/code/utils/tapas_physio_update_from_job.m b/PhysIO/code/utils/tapas_physio_update_from_job.m index e6460dc5..aab6a9c7 100644 --- a/PhysIO/code/utils/tapas_physio_update_from_job.m +++ b/PhysIO/code/utils/tapas_physio_update_from_job.m @@ -83,34 +83,43 @@ for p = 1:nProperties - currentProperty = eval(sprintf('job.%s', jobPropertyArray{p})); - - % overwrite properties of physio with sub-properties of job, also - % set value of branchProperty to name of current property - if isBranch{p} - valueBranchArray = fields(currentProperty); - valueBranch = valueBranchArray{1}; - eval(sprintf('physio.%s.%s = valueBranch;', ... - physioPropertyArray{p}, branchProperty{p})); - - currentProperty = currentProperty.(valueBranch); + try + currentProperty = eval(sprintf('job.%s', jobPropertyArray{p})); + catch err + currentProperty = []; + tapas_physio_log(sprintf('No property %s defined in job (error: %s)', ... + jobPropertyArray{p}, err.message), [], 1); end - - % update property itself, if it has no sub-properties - if ~isstruct(currentProperty) - eval(sprintf('physio.%s = currentProperty;', ... - physioPropertyArray{p})); - else + if ~isempty(currentProperty) % no error in retrieval => continue parsing! - % update all existing sub-properties in job to physio + % overwrite properties of physio with sub-properties of job, also + % set value of branchProperty to name of current property + if isBranch{p} + valueBranchArray = fields(currentProperty); + valueBranch = valueBranchArray{1}; + eval(sprintf('physio.%s.%s = valueBranch;', ... + physioPropertyArray{p}, branchProperty{p})); + + currentProperty = currentProperty.(valueBranch); + end - subPropArray = fields(currentProperty); - nFields = numel(subPropArray); - for f = 1:nFields - eval(sprintf('physio.%s.%s = currentProperty.%s;', ... - physioPropertyArray{p},subPropArray{f}, subPropArray{f})); + % update property itself, if it has no sub-properties + if ~isstruct(currentProperty) + eval(sprintf('physio.%s = currentProperty;', ... + physioPropertyArray{p})); + else + + % update all existing sub-properties in job to physio + + subPropArray = fields(currentProperty); + nFields = numel(subPropArray); + + for f = 1:nFields + eval(sprintf('physio.%s.%s = currentProperty.%s;', ... + physioPropertyArray{p},subPropArray{f}, subPropArray{f})); + end end end end \ No newline at end of file diff --git a/PhysIO/docs/documentation.html b/PhysIO/docs/documentation.html index fde9bd5f..84d7714c 100644 --- a/PhysIO/docs/documentation.html +++ b/PhysIO/docs/documentation.html @@ -6,10 +6,10 @@ PhysIO Toolbox Documentation

TAPAS HUGE demo

Contents

Introduction

This toolbox implements Hierarchical Unsupervised Generative Embedding (HUGE) with variational Bayesian inversion. This demo shows how to call HUGE on two synthetic datasets. For more information, consult REF [1].

Author: Yu Yao (yao@biomed.ee.ethz.ch) Copyright (C) 2018 Translational Neuromodeling Unit Institute for Biomedical Engineering, University of Zurich and ETH Zurich.

This toolbox is part of TAPAS, which is released under the terms of the GNU General Public Licence (GPL), version 3. For further details, see http://www.gnu.org/licenses/

This software is intended for research only. Do not use for clinical purpose. Please note that this toolbox is in an early stage of development. Considerable changes are planned for future releases. For support please refer to: https://github.com/translationalneuromodeling/tapas/issues

Generate synthetic DCM fMRI datasets

We generate two synthetic datasets:

  1. A dataset based on a three-region bilinear DCM with 80 subjects divided into three groups (40, 20, 20).
  2. A dataset based on a two-region linear DCM with 20 subjects divided into two groups of 10 subjects each.
rng(8032,'twister')
-tapas_huge_generate_examples( 'example_data.mat' )
-load( 'example_data.mat' );
-

Analyzing first dataset

To perform analysis with HUGE use the following interface:

[DcmResults] = tapas_huge_invert(DCM, K, priors, verbose, randomize, seed)

DCM: cell array of DCM in SPM format

disp(DCMr3b{1})
-
          U: [1×1 struct]
-          Y: [1×1 struct]
-          n: 3
-          v: 256
-         TE: 0.0400
-    options: [1×1 struct]
-          a: [3×3 logical]
-          c: [3×2 logical]
-          b: [3×3×2 logical]
-          d: [3×3×0 double]
-
-

K: number of clusters in the HUGE model

K = 3;
-

priors: model priors stored in a struct containing the following fields:

-
-
    -
  • alpha: parameter of Dirichlet prior (α0 - - in Fig.1 of REF [1])
    • -
  • clustersMean: prior mean of clusters (m0 - - in Fig.1 of REF [1])
    • -
  • clustersTau: - τ0 - - in Fig.1 of REF [1]
    • -
  • clustersDeg: degrees of freedom of inverse-Wishart prior (ν0 - - in Fig.1 of REF [1])
  • clustersSigma: scale matrix of inverse-Wishart prior (S0 - - in Fig.1 of REF [1])
  • hemMean: prior mean of hemodynamic parameters (μ0 - - in Fig.1 of REF [1])
  • hemSigma: prior covariance of hemodynamic parameters(Σ0 - - in Fig.1 of REF [1])
    • -
  • noiseInvScale: prior inverse scale of observation noise (b0 - - in Fig.1 of REF [1])
    • -
  • noiseShape: prior shape parameter of observation noise (a0 - - in Fig.1 of REF [1])
    • -
-
-

- You may use tapas_huge_build_prior(DCM) to generate this struct. -

-
priors = tapas_huge_build_prior(DCMr3b);
-disp(priors)
-
-
-            alpha: 1
-     clustersMean: [-0.5000 0.0052 0.0052 0.0052 -0.5000 0.0052 -0.5000 0 0]
-      clustersTau: 0.1000
-      clustersDeg: 100
-    clustersSigma: [9×9 single]
-          hemMean: [0 0 0 0 0 0 0]
-         hemSigma: [7×7 double]
-    noiseInvScale: 0.0250
-       noiseShape: 1.2800
-
-
-

verbose: activate command line output (optional).

-
-  verbose = true;
-
-

randomize: randomize starting values (optional). If randomize is false (default), VB inversion will start from the prior values.

-
-  randomize = true;
-
-

seed: seed for random number generator (optional). Use this input to reproduce an earlier result.

-
-  seed = rng;
-
-

Starting the inversion:

-
-currentTimer = tic;
-[DcmResults] = tapas_huge_invert(DCMr3b, K, priors, verbose, randomize, seed);
-toc(currentTimer)
-
iteration 1, dF: 206273.7158
-iteration 2, dF: 5089.5705
-iteration 3, dF: 935.7616
-iteration 4, dF: 141.2354
-iteration 5, dF: 13.6022
-iteration 6, dF: 1044.9978
-iteration 7, dF: 168.1947
-iteration 8, dF: 49.7
-iteration 9, dF: 33.17
-iteration 10, dF: 15.7876
-iteration 11, dF: 7.9238
-iteration 12, dF: 33.4518
-iteration 13, dF: 21.0559
-iteration 14, dF: 14.291
-iteration 15, dF: 3.2596
-iteration 16, dF: 2.4821
-iteration 17, dF: 1.4834
-iteration 18, dF: 2.4299
-iteration 19, dF: 1.7049
-iteration 20, dF: 0.17127
-iteration 21, dF: 1.6542
-iteration 22, dF: 1.6965
-iteration 23, dF: 0.58714
-iteration 24, dF: 0.57489
-iteration 25, dF: 0.45
-iteration 26, dF: 0.34384
-iteration 27, dF: 0.29236
-iteration 28, dF: 0.14673
-iteration 29, dF: 0.13925
-iteration 30, dF: 0.14582
-iteration 31, dF: 0.078373
-iteration 32, dF: 0.099591
-iteration 33, dF: 0.075957
-iteration 34, dF: 0.017069
-iteration 35, dF: 0.054854
-iteration 36, dF: 0.053601
-iteration 37, dF: -0.02852
-Elapsed time is 213.671655 seconds.
-
-

The inference result is stored in DcmResults.posterior, which is a struct containing the following fields:

-
  • alpha: parameter of posterior over cluster weights (αk - - in Eq.(15) of REF [1])
    • -
  • softAssign: posterior assignment probability of subjects to clusters (qnk - - in Eq.(18) in REF [1])
    • -
  • clustersMean: posterior mean of clusters (mk - - in Eq.(16) of REF [1])
  • clustersTau: - τk - -in Eq.(16) of REF [1]
  • clustersDeg: posterior degrees of freedom (νk - - in Eq.(16) of REF [1])
  • clustersSigma: posterior scale matrix (Sk - -in Eq.(16) of REF [1])
    • -
  • logDetClustersSigma: log-determinant of Sk - -
    • -
  • dcmMean: posterior mean of DCM parameters (μn - - in Eq.(19) of REF [1])
    • -
  • dcmSigma: posterior covariance of hemodynamic parameters (Σn - - in Eq.(19) of REF [1])
    • -
  • logDetPostDcmSigma: log-determinant of - Σn - -
    • -
  • noiseInvScale: posterior inverse scale of observation noise (bn,r - - in Eq.(21) of REF [1])
    • -
  • noiseShape: posterior shape parameter of observation noise (an,r - - in Eq.(21) of REF [1])
    • -
  • meanNoisePrecision: posterior mean of precision of observation noise (λn,r - - in Eq.(23) of REF [1])
    • -
  • modifiedSumSqrErr: b′n,r - -in Eq.(22) of REF [1]

The negative free energy after convergence is stored in DcmResults.freeEnergy

disp(['Negative free energy is ' num2str(DcmResults.freeEnergy) ...
-    ' after ' num2str(DcmResults.nIterationsActual) ' iterations.'])
-
Negative free energy is 166233.7053 after 37 iterations.
-

We can plot the results using the command

tapas_huge_plot(DCMr3b,DcmResults);
-
-DCMr3b -

- The assignment probabilities shown in the upper left panel of the above figure indicates that all but three subjects could be correctly classified.

The panel in the top right shows the posterior cluster mean and 95% marginal credible intervals.

The lower panel shows the simulated BOLD measurement (black) and 25 samples from the posterior over (noise free) BOLD signals. -

-

- Analyzing the second dataset -

-

- Performing HUGE analysis on the second dataset proceeds analogously to the first dataset. Differences in the DCM's structure are handled automatically. -

- 
disp(DCMr2l{1})
-
-% choose number of clusters for the HUGE model
-K = 2;
-
-% generate priors
-priors = tapas_huge_build_prior(DCMr2l);
-disp(priors)
-
-% suppress command line output
-verbose = false;
-
-% randomize starting values
-randomize = true;
-
-% analysis
-currentTimer = tic;
-[DcmResults] = tapas_huge_invert(DCMr2l, K, priors, verbose, randomize);
-toc(currentTimer)
-
-disp(['Negative free energy is ' num2str(DcmResults.freeEnergy) ...
-    ' after ' num2str(DcmResults.nIterationsActual) ' iterations.'])
-
-% plot the results
-tapas_huge_plot(DCMr2l,DcmResults);
-
          U: [1×1 struct]
-          Y: [1×1 struct]
-          n: 2
-          v: 300
-         TE: 0.0400
-    options: [1×1 struct]
-          a: [2×2 logical]
-          c: [2×2 logical]
-          b: [2×2×2 logical]
-          d: [2×2×0 double]
-
-            alpha: 1
-     clustersMean: [-0.5000 0.0078 -0.5000 0 0]
-      clustersTau: 0.1000
-      clustersDeg: 100
-    clustersSigma: [5×5 single]
-          hemMean: [0 0 0 0 0]
-         hemSigma: [5×5 double]
-    noiseInvScale: 0.0250
-       noiseShape: 1.2800
-
-Elapsed time is 49.033623 seconds.
-Negative free energy is -10117.8605 after 56 iterations.
-
-IMG: DCMr2l -

Empirical Bayes analysis

To perform empirical Bayes analysis, set K to 1. Here, we perform empirical Bayesian analysis on the second dataset. For more information on empirical Bayes, see Fig.5 in REF [1].

% set K to 1 for empirical Bayes
-K = 1;
-
-% generate priors
-priors = tapas_huge_build_prior(DCMr2l);
-
-% suppress command line output
-verbose = false;
-
-% randomize starting values
-randomize = true;
-
-% analysis
-currentTimer = tic;
-[DcmResults] = tapas_huge_invert(DCMr2l, K, priors, verbose, randomize);
-toc(currentTimer)
-
-disp(['Negative free energy is ' num2str(DcmResults.freeEnergy) ...
-    ' after ' num2str(DcmResults.nIterationsActual) ' iterations.'])
-
-% plot the results
-tapas_huge_plot(DCMr2l,DcmResults);
-
Elapsed time is 31.027961 seconds.
-Negative free energy is -9913.7702 after 38 iterations.
-
-EB -

The above figure summarizes the result of empirical Bayes. The top panel shows the posterior group-level mean with 95% marginal credible intervals.

The lower panel shows a boxplot of first-level (i.e.: subject-level) MAP estimates of the DCM parameters.

References:

[1] Yao Y, Raman SS, Schiek M, Leff A, Frässle S, Stephan KE (2018). Variational Bayesian Inversion for Hierarchical Unsupervised Generative Embedding (HUGE). NeuroImage, 179: 604-619 https://doi.org/10.1016/j.neuroimage.2018.06.073


Published with MATLAB® R2017a

diff --git a/huge/tapas_huge_bold.m b/huge/tapas_huge_bold.m index 3b045e56..82323fad 100644 --- a/huge/tapas_huge_bold.m +++ b/huge/tapas_huge_bold.m @@ -1,23 +1,39 @@ -%% [response,x,s,f,v,q] = tapas_huge_bold( packedParameters, DcmInfo, iSubject ) +function [response, x, s, f, v, q] = tapas_huge_bold( A, B, C, D, tau, ... + kappa, epsilon, R, u, L, E_0, r_0, V_0, vartheta_0, alpha, gamma, TR, ... + TE, dt) +% Integrates the DCM forward equations to generate the predicted fMRI bold +% time series. % -% Compute the predicted response for the DCM-FMRI model using numerical -% integration of differential equations. -% -% INPUT: -% packedParameters - current value of parameters -% DcmInfo - struct containing DCM model specification and -% BOLD time series. -% iSubject - subject index +% INPUTS: +% A, B, C, D - DCM connectivity matrices. +% tau - Venous transit time. +% kappa - Decay of vasodilatory signal. +% epsilon - Ratio of intra- and extravascular signal. +% R - Number of regions. +% u - Experimental stimuli. +% L - Number of experimental stimuli. +% E_0 - Resting oxygen extraction fraction. +% r_0 - Slope of intravascular relaxation rate. +% V_0 - Resting venous volume. +% vartheta_0 - Frequency offset at the outer surface of magnetized +% vessels (Hz). +% alpha - Grubb's exponent. +% gamma - rate constant of feedback regulation. +% TR - Repetition time. +% TE - Echo time. +% dt - Sampling interval of inputs. % -% OUTPUT: -% response - matrix of predicted response for each region +% OUTPUTS: +% response - matrix of predicted response for each region % (column-wise) -% x - time series of neuronal states -% s - time series of vasodilatory signal -% f1 - time series of flow -% v1 - time series of blood volume -% q1 - time series of deoxyhemoglobin content. -% +% x - time series of neuronal states +% s - time series of vasodilatory signal +% f1 - time series of flow +% v1 - time series of blood volume +% q1 - time series of deoxyhemoglobin content. +% + +% % REFERENCE: % Klaas Enno Stephan, Nikolaus Weiskopf, Peter M. Drysdale, Peter A. % Robinson, Karl J. Friston (2007). Comparing hemodynamic models with @@ -26,111 +42,76 @@ % https://doi.org/10.1016/j.neuroimage.2007.07.040 % -% -% Author: Sudhir Shankar Raman -% Copyright (C) 2018 Translational Neuromodeling Unit +% Author: Yu Yao (yao@biomed.ee.ethz.ch), Sudhir Shankar Raman +% Copyright (C) 2019 Translational Neuromodeling Unit % Institute for Biomedical Engineering, % University of Zurich and ETH Zurich. % % This file is part of TAPAS, which is released under the terms of the GNU % General Public Licence (GPL), version 3. For further details, see -% . +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. % % This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -%% -function [response,x,s,f,v,q] = tapas_huge_bold( packedParameters, DcmInfo, iSubject ) +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% -paramList = [DcmInfo.timeStep(iSubject) DcmInfo.nTime(iSubject) ... - DcmInfo.nStates DcmInfo.nInputs iSubject DcmInfo.dcmTypeB ... - DcmInfo.dcmTypeD]; -[packedParameters] = tapas_huge_pack_params(packedParameters,paramList); +nt = size(u, 1); +rSmp = TR/dt; +C = C'/rSmp; -%---------------------------------------------------------------- -% unpack parameters -HEM_LIST = DcmInfo.hemParam.listHem; % kappa, tau, epsilon -% transpose of the original A -A = (packedParameters{1}.*DcmInfo.adjacencyA)'; -% transpose of the original C -C = (packedParameters{2}.*DcmInfo.adjacencyC)'; -% For Lorenz Dataset Analysis -C = C./DcmInfo.hemParam.scaleC; -B = (packedParameters{3}); % original B -D = packedParameters{4}; -tau = HEM_LIST(2)*exp(packedParameters{5}(1,:)); -kappa = HEM_LIST(1)*exp(packedParameters{5}(2,:)); -epsilon = HEM_LIST(3)*exp(packedParameters{5}(3,:)); -epsilon = epsilon(1)*ones(1, DcmInfo.nStates); -% estimated region-specific ratios of intra- to extra-vascular signal -%-------------------------------------------------------------------------- +% transform self-connections +A = A - .5*eye(R); +if isempty(D) + D = zeros(R, R, R); +end -%------------------------------------------------------------------ -% Hemodynamic constants -% TE -echoTime = DcmInfo.hemParam.echoTime; +tau = 2.*exp(tau); +kappa = .64.*exp(kappa); +epsilon = exp(epsilon); -% resting venous volume (%) restingVenousVolume -%-------------------------------------------------------------------------- -restingVenousVolume = DcmInfo.hemParam.restingVenousVolume; - -% slope relaxationRateSlope of intravascular relaxation rate R_iv as a -% function of oxygen saturation S: -% R_iv = relaxationRateSlope*[(1 - S)-(1 - S0)] (Hz) -%-------------------------------------------------------------------------- -relaxationRateSlope = DcmInfo.hemParam.relaxationRateSlope; - -% frequency offset at the outer surface of magnetized vessels (Hz) - nu0 -%-------------------------------------------------------------------------- -frequencyOffset = DcmInfo.hemParam.frequencyOffset; - -% resting oxygen extraction fraction - rho -%-------------------------------------------------------------------------- -oxygenExtractionFraction = DcmInfo.hemParam.oxygenExtractionFraction*... - ones(1,paramList(3)); +% resting oxygen extraction fraction +E_0 = repmat(E_0, 1, R); -%-Coefficients in BOLD signal model - -%========================================================================== -coefficientK1 = DcmInfo.hemParam.rho*frequencyOffset*echoTime*... - oxygenExtractionFraction; -coefficientK2 = epsilon.*(relaxationRateSlope*... - oxygenExtractionFraction*echoTime); -coefficientK3 = 1 - epsilon; +k1 = 4.3*vartheta_0*TE*E_0; +k2 = epsilon.*(r_0*E_0*TE); +k3 = 1 - epsilon; -mArrayB = B.*DcmInfo.adjacencyB; -mArrayB = permute(mArrayB,[2 1 3]); -mArrayD = D.*DcmInfo.adjacencyD; -mArrayD = permute(mArrayD,[2 1 3]); -% resting oxygen extraction fraction - change value to match what SPM does -%-------------------------------------------------------------------------- -% For Lorenz Dataset Analysis -oxygenExtractionFraction = DcmInfo.hemParam.oxygenExtractionFraction2*... - ones(1,paramList(3)); -%oxygenExtractionFraction = 0.32*ones(1,paramList(3)); +% Integrate the dynamical system +[x,s,f,v,q] = tapas_huge_int_euler(... + A',... + full(u*C),... + full(u),... + permute(B,[2 1 3]),... + permute(D,[2 1 3]),... + E_0,... + 1/alpha,... + tau,... + gamma,... + kappa,... + [dt,nt,R,L,0,any(B(:)),any(D(:))]); -% pre-calculation -C = DcmInfo.listU{paramList(5)}*C; -% Integrate the dynamical system -[x,s,f,v,q] = tapas_huge_int_euler(A,C,DcmInfo.listU{paramList(5)},... - mArrayB,mArrayD,oxygenExtractionFraction,DcmInfo.hemParam.alphainv,... - tau,DcmInfo.hemParam.gamma,kappa,paramList); - -% generate the responses per time point -response = restingVenousVolume*( ... - bsxfun(@times,coefficientK1,... - (1 - (q(DcmInfo.listResponseTimeIndices{paramList(5)},:)))) +... - bsxfun(@times,coefficientK2,... - (1 - (q(DcmInfo.listResponseTimeIndices{paramList(5)},:)./... - v(DcmInfo.listResponseTimeIndices{paramList(5)},:)))) +... - bsxfun(@times,coefficientK3,... - (1-v(DcmInfo.listResponseTimeIndices{paramList(5)},:)))); - +% generate the BOLD response +response = V_0*( ... + bsxfun(@times,k1,... + (1 - (q(rSmp:rSmp:end,:)))) +... + bsxfun(@times,k2,... + (1 - (q(rSmp:rSmp:end,:)./... + v(rSmp:rSmp:end,:)))) +... + bsxfun(@times,k3,... + (1-v(rSmp:rSmp:end,:)))); +% demean response +response = bsxfun(@minus, response, mean(response, 1)); +end \ No newline at end of file diff --git a/huge/tapas_huge_boxcar.m b/huge/tapas_huge_boxcar.m new file mode 100644 index 00000000..0124f867 --- /dev/null +++ b/huge/tapas_huge_boxcar.m @@ -0,0 +1,106 @@ +function [ u ] = tapas_huge_boxcar( dt, nBoxes, period, onRatio, padding ) +% Generate a boxcar function for use as experimental stimulus. All +% timing-related arguments must be specified in seconds. +% +% INPUTS: +% dt - Numeric scalar indicating sampling time interval. +% nBoxes - Vector indicating number of blocks. +% period - Vector containing time interval between block onsets. +% onRatio - Vector containing ratio between block length and 'period'. +% Must be between 0 and 1. +% +% OPTIONAL INPUTS: +% padding - Length of padding at the beginning and end. +% +% OUTPUTS: +% u - A cell array containing the boxcar functions. +% +% EXAMPLES: +% u = TAPAS_HUGE_BOXCAR(.01, 10, 3, 2/3, [4 0]) Generate boxcar +% function with 10 blocks, each 2 seconds long with 1 second inter +% block interval and onset of first block at 4 seconds. +% +% See also tapas_Huge.SIMULATE +% + +% Author: Yu Yao (yao@biomed.ee.ethz.ch) +% Copyright (C) 2019 Translational Neuromodeling Unit +% Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of TAPAS, which is released under the terms of the GNU +% General Public Licence (GPL), version 3. For further details, see +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. +% +% This software is intended for research only. Do not use for clinical +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% + + +%% check input +nSignals = numel(nBoxes); + +if isscalar(period) && nSignals > 1 + period = repmat(period, nSignals, 1); +else + assert(numel(period) == nSignals, 'TAPAS:HUGE:Boxcar:InputSize', ... + 'Size of period must match that of nBoxes') +end + +if isscalar(onRatio) && nSignals > 1 + onRatio = repmat(onRatio, nSignals, 1); +else + assert(numel(period) == nSignals, 'TAPAS:HUGE:Boxcar:InputSize', ... + 'Size of onRatio must match that of nBoxes') +end +assert(all(onRatio < 1) && all(onRatio > 0), 'TAPAS:HUGE:Boxcar:InputRange', ... + 'onRatio must be in range: 0 < onRatio < 1.') + +if nargin < 5 + padding = zeros(nSignals, 2); +else + if size(padding, 1) == 1 + padding = repmat(padding, nSignals, 1); + end + assert(size(padding, 1) == nSignals, 'TAPAS:HUGE:Boxcar:InputSize', ... + 'Number of rows in padding must match length of nBoxes') +end + + +%% generate signal +u = cell(1, nSignals); + +for iSignal = 1:nSignals + tMax = nBoxes(iSignal)*period(iSignal) + padding(iSignal, end); + + % amplitudes + amp = repmat([1;0], 1, nBoxes(iSignal)); + + % sample time points + boxStarts = (0:nBoxes(iSignal)-1)*period(iSignal); + boxDuration = [0; onRatio(iSignal)]*period(iSignal); + grid = bsxfun(@plus, boxStarts, boxDuration); + + % query time points + query = (0:dt:tMax)'; + + % padding + amp = [amp(:); 0]; + grid = [grid(:); tMax]; + if padding(iSignal, 1) + query = [-fliplr(dt:dt:padding(iSignal, 1))'; query(:)]; + amp = [0; amp(:)]; + grid = [-padding(iSignal, 1); grid(:)]; + end + % generate boxcar + u{iSignal} = interp1(grid, amp, query, 'previous'); + +end + diff --git a/huge/tapas_huge_bpurity.m b/huge/tapas_huge_bpurity.m new file mode 100644 index 00000000..9bcea07b --- /dev/null +++ b/huge/tapas_huge_bpurity.m @@ -0,0 +1,80 @@ +function [ balancedPurity ] = tapas_huge_bpurity( labels, estimates ) +% Calculate balanced purity (see Brodersen2014 Eq. 13 and 14) for a set of +% ground truth labels and a set of estimated labels +% +% INPUTS: +% labels - Vector of ground truth class labels. +% estimates - Clustering result as array of assignment probabilities or +% vector of cluster indices. +% +% OUTPUTS: +% balancedPurity - Balanced purity score according to Brodersen (2014) +% +% EXAMPLES: +% bp = TAPAS_HUGE_BPURITY(labels,estimates) +% + +% Author: Yu Yao (yao@biomed.ee.ethz.ch) +% Copyright (C) 2019 Translational Neuromodeling Unit +% Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of TAPAS, which is released under the terms of the GNU +% General Public Licence (GPL), version 3. For further details, see +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. +% +% This software is intended for research only. Do not use for clinical +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% + + +labels = labels(:); +if numel(estimates) == numel(labels) + estimates = estimates(:); +end +[N, K] = size(estimates); +assert(N == numel(labels), 'TAPAS:UTIL:InputSizeMismatch', ... + 'Number of rows in estimates must match number of elements in labels.'); +if K > 1 + [~, estimates] = max(estimates, [], 2); +else + assert(all(mod(estimates, 1) == 0) && all(estimates > 0), ... + 'TAPAS:UTIL:NonIntegerIndices', ... + 'Cluster indices must be positive integer.'); + K = max(estimates); +end + +% number of classes +C = max(max(labels), K); + +% degree of imbalance +xi = zeros(C, 1); +for c = 1:C + xi(c) = nnz(labels == c)/N; +end +xi = max(xi); + +% calculate purity (Brodersen2014 Eq. 13) +counts = zeros(C,C); +for k = 1:C + % class labels of samples grouped into cluster k + currentLabels = labels(estimates == k); + for c = 1:C + % number of samples belonging to class c + counts(c,k) = nnz(currentLabels == c); + end +end +purity = sum(max(counts))/N; + +% Brodersen2014 Eq. 14 (n->k) +balancedPurity = (1-1/C)*(purity - xi)/(1 - xi) + 1/C; + +end + diff --git a/huge/tapas_huge_build_prior.m b/huge/tapas_huge_build_prior.m index 62b191f4..3f643089 100644 --- a/huge/tapas_huge_build_prior.m +++ b/huge/tapas_huge_build_prior.m @@ -1,84 +1,92 @@ -%% [ priors, DcmInfo ] = tapas_huge_build_prior( DcmInfo ) -% -% Generates values for prior parameters for HUGE. Prior mean of cluster +function [ priors, DCM ] = tapas_huge_build_prior( DCM ) +% WARNING: This function is deprecated and will be removed in a future +% version of this toolbox. Please use the new object-oriented interface +% provided by the tapas_Huge class. +% +% Generate values for prior parameters for HUGE. Prior mean of cluster % centers and prior mean and covariance of hemodynamic parameters follow % SPM convention (SPM8 r6313). % % INPUT: -% DcmInfo - EITHER: struct containing DCM model specification and -% BOLD time series in DcmInfo format (see -% tapas_huge_simulate.m for an example) -% OR: cell array of DCM in SPM format +% DcmInfo - cell array of DCM in SPM format % % OUTPUT: % priors - struct containing priors % DcmInfo - struct containing DCM model specification and BOLD time % series in DcmInfo format % -% REFERENCE: -% [1] Yao Y, Raman SS, Schiek M, Leff A, Frssle S, Stephan KE (2018). -% Variational Bayesian Inversion for Hierarchical Unsupervised -% Generative Embedding (HUGE). NeuroImage, 179: 604-619 +% See also tapas_Huge, tapas_Huge.estimate, tapas_huge_demo % -% https://doi.org/10.1016/j.neuroimage.2018.06.073 -% % Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit +% Copyright (C) 2019 Translational Neuromodeling Unit % Institute for Biomedical Engineering, % University of Zurich and ETH Zurich. % % This file is part of TAPAS, which is released under the terms of the GNU % General Public Licence (GPL), version 3. For further details, see -% . +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. % % This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: % https://github.com/translationalneuromodeling/tapas/issues -function [ priors, DcmInfo ] = tapas_huge_build_prior( DcmInfo ) +% + +wnMsg = ['This function is deprecated and will be removed in a future ' ... + 'version of this toolbox. Please use the new object-oriented ' ... + 'interface provided by the tapas_Huge class.']; +warning('tapas:huge:deprecated',wnMsg) + %% check input format -if ~isfield(DcmInfo,'listBoldResponse') +if isvector(DCM)&&isstruct(DCM) try - DcmInfo = tapas_huge_import_spm(DcmInfo); - catch err - disp(['TAPAS:HUGE: Unsupported format. '... - 'Use cell array of DCM in SPM format as first input.']) - rethrow(err); + DCM = {DCM(:).DCM}'; + catch + DCM = num2cell(DCM); end +else + assert(iscell(DCM),'TAPAS:HUGE:inputFormat',... + 'DCM must be cell array of DCMs in SPM format'); end +dcm = DCM{1}; + %% set priors priors = struct(); - % parameter of Dirichlet prior (alpha_0 in Figure 1 of REF [1]) priors.alpha = 1; -tmp = DcmInfo.adjacencyA/64/DcmInfo.nStates; -tmp = tmp - diag(diag(tmp)) - .5*eye(DcmInfo.nStates); -tmp = [tmp(:); DcmInfo.adjacencyC(:)*0; DcmInfo.adjacencyB(:)*0; ... - DcmInfo.adjacencyD(:)*0]; +tmp = dcm.a/64/dcm.n; +tmp = tmp - diag(diag(tmp)) - .5*eye(dcm.n); +tmp = [tmp(:); dcm.b(:)*0; dcm.c(:)*0; ... + dcm.d(:)*0]; +connectionIndicator = find([dcm.a(:);dcm.b(:);dcm.c(:);dcm.d(:)]); % prior mean of clusters (m_0 in Figure 1 of REF [1]) -priors.clustersMean = tmp(DcmInfo.connectionIndicator)'; +priors.clustersMean = tmp(connectionIndicator)'; % tau_0 in Figure 1 of REF [1] priors.clustersTau = 0.1; % degrees of freedom of inverse-Wishart prior (nu_0 in Figure 1 of REF [1]) -priors.clustersDeg = max(100,1.5^DcmInfo.nConnections); +priors.clustersDeg = max(100,1.5^length(connectionIndicator)); priors.clustersDeg = min(priors.clustersDeg,double(realmax('single'))); % scale matrix of inverse-Wishart prior (S_0 in Figure 1 of REF [1]) -priors.clustersSigma = 0.01*eye(DcmInfo.nConnections)*... - (priors.clustersDeg - DcmInfo.nConnections - 1); +priors.clustersSigma = 0.01*eye(length(connectionIndicator))*... + (priors.clustersDeg - length(connectionIndicator) - 1); % prior mean of hemodynamic parameters (mu_h in Figure 1 of REF [1]) -priors.hemMean = zeros(1,DcmInfo.nStates*2 + 1); +priors.hemMean = zeros(1,dcm.n*2 + 1); % prior Covariance of hemodynamic parameters(Sigma_h in Figure 1 of % REF [1]) -priors.hemSigma = diag(zeros(1,DcmInfo.nStates*2 + 1)+exp(-6)); +priors.hemSigma = diag(zeros(1,dcm.n*2 + 1)+exp(-6)); % prior inverse scale of observation noise (b_0 in Figure 1 of REF [1]) priors.noiseInvScale = .025; diff --git a/huge/tapas_huge_compile.m b/huge/tapas_huge_compile.m index fbb9a875..7988dc60 100644 --- a/huge/tapas_huge_compile.m +++ b/huge/tapas_huge_compile.m @@ -1,27 +1,33 @@ -%% [ ] = tapas_huge_compile( ) -% -% Compiles integrator for DCM. Assumes that source code file -% 'tapas_huge_int_euler.c' is in the same folder as this function and -% places the compiled mex-file in the same folder. +function [ ] = tapas_huge_compile( ) +% Compile integrator for DCM fMRI forward model. This function assumes that +% source code file 'tapas_huge_int_euler.c' is in the same folder as this +% function and places the compiled mex-file in the same folder. Running +% this function requires a C compiler to be installed. For more information +% see: +% +% https://www.mathworks.com/support/requirements/supported-compilers.html % - % Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit +% Copyright (C) 2019 Translational Neuromodeling Unit % Institute for Biomedical Engineering, % University of Zurich and ETH Zurich. % % This file is part of TAPAS, which is released under the terms of the GNU % General Public Licence (GPL), version 3. For further details, see -% . +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. % % This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: % https://github.com/translationalneuromodeling/tapas/issues -% -function [ ] = tapas_huge_compile( ) +% + %% % check if mex file already exists if exist('tapas_huge_int_euler', 'file') == 3 diff --git a/huge/tapas_huge_count_params.m b/huge/tapas_huge_count_params.m deleted file mode 100644 index 42372c02..00000000 --- a/huge/tapas_huge_count_params.m +++ /dev/null @@ -1,59 +0,0 @@ -%% [ nParameters, idxParamsInf, idxSelfCon ] = tapas_huge_count_params( DcmInfo ) -% -% Calculate number of parameters and parameter indices for DCM. -% -% INPUT: -% DcmInfo - struct containing DCM model specification and BOLD time -% series. -% -% OUTPUT: -% nParameters - 2-by-3 array containing parameter counts -% idxParamsInf - indices of parameters being inferred -% idxSelfCon - indices of self connections -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ nParameters, idxParamsInf, idxSelfCon ] = tapas_huge_count_params( DcmInfo ) -% length of full DCM parameter vector -nDcmParametersCon = DcmInfo.nStates*(DcmInfo.nStates + ... % A - DcmInfo.nStates*DcmInfo.nInputs + ... % B - DcmInfo.nInputs + ... % C - DcmInfo.nStates*DcmInfo.nStates); % D -nDcmParametersHem = DcmInfo.nStates*3; % hem -nDcmParametersAll = nDcmParametersCon + nDcmParametersHem; - -% number of parameters to be inferred %%% dim(theta_c) dim(theta_h) -idxParamsInf = (1:nDcmParametersAll).'; -idxParamsInf(DcmInfo.noConnectionIndicator) = 0; -idxParamsInf(end-DcmInfo.nStates+2:end) = 0; % only first epsilon is used -idxParamsInf = idxParamsInf(idxParamsInf~=0); - -nDcmParamsInfCon = nnz(idxParamsInf<=nDcmParametersCon); %%% d_c -% nDcmParamsInfCon = DcmInfo.nConnections; % alternative -nDcmParamsInfHem = nnz(idxParamsInf>nDcmParametersCon); %%% d_h -nDcmParamsInfAll = length(idxParamsInf); %%% d - -% indices of diagonal elements of A (self connections) -idxSelfCon = find(eye(size(DcmInfo.adjacencyA))); -idxSelfCon = intersect(idxParamsInf,idxSelfCon); - -nParameters = [nDcmParametersCon,nDcmParametersHem,nDcmParametersAll;... - nDcmParamsInfCon,nDcmParamsInfHem,nDcmParamsInfAll;]; - - -end - diff --git a/huge/tapas_huge_demo.m b/huge/tapas_huge_demo.m index 4994a4fd..580c5703 100644 --- a/huge/tapas_huge_demo.m +++ b/huge/tapas_huge_demo.m @@ -1,230 +1,270 @@ -%% TAPAS HUGE demo +%% Hierarchical Unsupervised Generative Embedding %% Introduction -% This toolbox implements Hierarchical Unsupervised Generative Embedding -% (HUGE) with variational Bayesian inversion. This demo shows how to call -% HUGE on two synthetic datasets. For more information, consult REF [1]. -% - -%% -% -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This toolbox is part of TAPAS, which is released under the terms of the -% GNU General Public Licence (GPL), version 3. For further details, see -% +% This toolbox implements Hierarchical Unsupervised Generative Embedding (HUGE) +% with variational Bayesian inversion. This demo script provides a quick-start +% guide to using the HUGE toolbox to stratify heterogeneous cohorts via generative +% embedding or perform empirical Bayesian analysis. For a more detailed documentation, +% see the HUGE manual. For more information on the theory behind HUGE, please +% consult _Yao et al. Variational Bayesian Inversion for Hierarchical Unsupervised +% Generative Embedding (HUGE). NeuroImage (2018)_, available at . +%% Generating Synthetic Data +% Before applying HUGE, we need some data. Here, we use a simulated dataset, +% which we generate using the toolbox itself. The dataset is based on a three-region +% bilinear DCM with 20 subjects divided into two groups of 10 subjects each. % -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% -% - - -%% Generate synthetic DCM fMRI datasets -% We generate two synthetic datasets: +% First, we define the network structure of the DCM and generate the experimental +% stimuli. +%% +% fix random number generator seed for reproducible results +rng(8032, 'twister') + +% define DCM network structure +dcm = struct( ); +dcm.n = 3; +dcm.a = logical([ ... + 1 0 0; ... + 1 1 1; ... + 1 1 1; +]); +dcm.c = false(dcm.n, 3); +dcm.c([1, 5]) = true; +dcm.b = false(dcm.n, dcm.n, 3); +dcm.b(:, :, 3) = logical([ ... + 0 0 0; ... + 1 0 1; ... + 1 1 0; ... +]); +dcm.d = false(dcm.n, dcm.n, 0); + + +% generate experimental stimuli +U = struct(); +U.dt = 1.84/16; +tmp = tapas_huge_boxcar(U.dt, [24*13 24], [2 26], [3/4 16/26], [0 0;0 0]); +nSmp = length(tmp{1}) + 160; +tmp{1}(nSmp) = 0; +tmp{2}(nSmp) = 0; +tmp{2} = tmp{1}.*tmp{2}; +tmp{3} = zeros(1, 24); +tmp{3}([2 3 1 4 3 1] + (0:5)*4) = 1; +tmp{3} = reshape(repmat(tmp{3}, round(26/U.dt), 1), [], 1); +tmp{3}(nSmp) = 0; +tmp{3} = tmp{3}.*tmp{2}; +tmp{4} = zeros(1, 24); +tmp{4}([4 2 3 1 1 4] + (0:5)*4) = 1; +tmp{4} = reshape(repmat(tmp{4}, round(26/U.dt), 1), [], 1); +tmp{4}(nSmp) = 0; +tmp{4} = tmp{4}.*tmp{2}; +tmp = tmp(2:4); + +U.u = circshift(cell2mat(tmp), 6*16, 1); +U.name = {'1st stim', '2nd stim', '3rd stim'}; +dcm.U = U; +dcm.Y.dt = 16*dcm.U.dt; +dcm.Y.name = {'1st reg'; '2nd reg'; '3rd reg'}; +dcm.TE = .03; + +% plot stimuli +figure; +for ip = 1:3 + subplot(3, 1, ip); + plot(U.u(:, ip)); + ylabel(U.name{ip}) +end +xlabel('sample') +%% +% Now, we specify two subgroups within our simulated cohort, which differ +% in their effective connectivity profile. The first group prefers a bottom-up +% configuration (1st region -> 2nd region -> 3rd region) while the second group +% prefers a top-down configuration (1st region -> 3rd region -> 2nd region). +%% +sigma = .141; % group standard deviation +listGroups = cell(2, 1); + +% group 1 (bottom-up) +dcm.Ep.A = [-.1 .0 .0; ... + .2 -.1 .0; ... + .0 .1 -.1;]; +dcm.Ep.B = zeros(dcm.n, dcm.n, 3); +dcm.Ep.B(2, 1, 3) = .2; +dcm.Ep.B(3, 2, 3) = .35; +dcm.Ep.C = double(dcm.c)*.5; +dcm.Ep.D = double(dcm.d); +dcm.Ep.transit = zeros(dcm.n,1); +dcm.Ep.decay = zeros(dcm.n,1); +dcm.Ep.epsilon = 0; +tmp = [dcm.a(:);dcm.b(:);dcm.c(:);dcm.d(:)]; +dcm.Cp = diag([double(tmp).*sigma.^2; ones(2*dcm.n+1, 1)*exp(-6)]); +listGroups{1} = dcm; + +% group 2 (top-down) +dcm.Ep.A = [-.1 .0 .0; ... + .0 -.1 .1; ... + .4 -.15 -.1;]; +dcm.Ep.B = zeros(dcm.n, dcm.n, 3); +dcm.Ep.B(3, 1, 3) = .35; +dcm.Ep.B(2, 3, 3) = .35; +listGroups{2} = dcm; + +%% +% If you are familiar with SPM, you may have noticed that the data format +% we use to specify our desired DCM network structure resembles the DCM structure +% used by SPM. The HUGE toolbox uses the DCM structure of SPM to import experimental +% data, export results or specify the network structure for simulating data. % -% # A dataset based on a three-region bilinear DCM with 80 subjects divided -% into three groups (40, 20, 20). -% # A dataset based on a two-region linear DCM with 20 subjects divided -% into two groups of 10 subjects each. -% -rng(8032,'twister') -tapas_huge_generate_examples( 'example_data.mat' ) -load( 'example_data.mat' ); - -%% Analyzing first dataset -% To perform analysis with HUGE use the following interface: -% -% [DcmResults] = tapas_huge_invert(DCM, K, priors, verbose, randomize, seed) -% - +% The HUGE toolbox itself is implemented as a Matlab class. Hence, the first +% step to using HUGE is to create an instance of the class using the _tapas_HUGE_ +% command: %% -% |DCM|: cell array of DCM in SPM format -disp(DCMr3b{1}) - +hugeSim = tapas_Huge('Tag', 'simulated 3-region DCM'); +%% +% The _tapas_HUGE_ command takes optional parameters in the form of name-value +% pair arguments. In the above example, we added a tag to the newly created HUGE +% object, containing a short description of the object. To get a list of all available +% name-value pairs, type: + +help tapas_huge_property_names +%% +% Now, we are ready to generate the dataset by simulating fMRI time series +% for each subject. After simulation, the ground truth model parameters are stored +% in the class property _model_. The fMRI time series (simulated or experimental) +% are stored in the class property called _data_ and can be exported to SPM's +% DCM format using the _export_ function. %% -% |K|: number of clusters in the HUGE model -K = 3; - +groupSizes = [10 10]; % simulate 10 subjects for each group +snr = .5; % using a signal-to-noise ratio of 0.5 (-3 dB) +hugeSim = hugeSim.simulate(listGroups, groupSizes, 'Snr', snr); + +% plot simulated fMRI time series for one subject +n = 1; +figure; +for ip = 1:3 + subplot(3, 1, ip); + plot(hugeSim.data(n).bold(:, ip)); + ylabel(hugeSim.labels.regions{ip}) +end +xlabel('scans') + +% export data to SPM's DCM format +[ listDcms ] = hugeSim.export(); + +%% Stratification of Heterogeneous Cohorts +% Using the simulated dataset from the previous section, we demonstrate how +% to stratify the cohort into the bottom-up and top-down subgroups. All we have +% to do is to call the method _estimate_. %% -% |priors|: model priors stored in a struct containing the following -% fields: -% -% * |alpha|: parameter of Dirichlet prior ($\alpha_0$ in Fig.1 of -% REF [1]) -% * |clustersMean|: prior mean of clusters ($m_0$ in Fig.1 of REF [1]) -% * |clustersTau|: $\tau_0$ in Fig.1 of REF [1] -% * |clustersDeg|: degrees of freedom of inverse-Wishart prior ($\nu_0$ -% in Fig.1 of REF [1]) -% * |clustersSigma|: scale matrix of inverse-Wishart prior ($S_0$ in Fig.1 -% of REF [1]) -% * |hemMean|: prior mean of hemodynamic parameters ($\mu_h$ in -% Fig.1 of REF [1]) -% * |hemSigma|: prior covariance of hemodynamic parameters($\Sigma_h$ -% in Fig.1 of REF [1]) -% * |noiseInvScale|: prior inverse scale of observation noise ($b_0$ in -% Fig.1 of REF [1]) -% * |noiseShape|: prior shape parameter of observation noise ($a_0$ in -% Fig.1 of REF [1]) -% +% invert the HUGE model for the current dataset +% (this may take some time) +hugeSim = hugeSim.estimate('K', 2, 'Verbose', true); +%% +% This method accepts the same name-value pair arguments as the _tapas_Huge_ +% command we used to create the HUGE object in the last section. In the example +% above, we invert the HUGE model with two clusters and also activate command +% line outputs. The result is saved in the property _posterior_ and can be plotted +% using the method _plot_. Since, we are using simulated data for which ground +% truth group assignments are available, the model will automatically calculate +% the balanced purity of the estimation result. %% -% You may use |tapas_huge_build_prior(DCM)| to generate this struct. -% -priors = tapas_huge_build_prior(DCMr3b); -disp(priors) +% plot result +plot(hugeSim) -%% -% |verbose|: activate command line output (optional). -verbose = true; +% negative free energy +fprintf('Negative free energy: %e.\n', hugeSim.posterior.nfe) -%% -% |randomize|: randomize starting values (optional). If randomize is false -% (default), VB inversion will start from the prior values. -randomize = true; +% balanced accuracy +fprintf('Balanced accuracy: %f.\n', hugeSim.posterior.bPurity) -%% -% |seed|: seed for random number generator (optional). Use this input to -% reproduce an earlier result. -seed = rng; +% export result to SPM's DCM format +resultK2 = hugeSim.export(); -%% -% Starting the inversion: -% -currentTimer = tic; -[DcmResults] = tapas_huge_invert(DCMr3b, K, priors, verbose, randomize, seed); -toc(currentTimer) -%% -% The inference result is stored in |DcmResults.posterior|, which is a -% struct containing the following fields: +%% +% Note how plotting the object generates two graphs. The first one shows +% the estimated group assignments for each subject. The second one shows the posterior +% estimates of the group-level DCM connectivity parameters for each group. % -% * |alpha|: parameter of posterior over cluster weights -% ($\alpha_k$ in Eq.(15) of REF [1]) -% * |softAssign|: posterior assignment probability of subjects to -% clusters ($q_{nk}$ in Eq.(18) in REF [1]) -% * |clustersMean|: posterior mean of clusters ($m_k$ in Eq.(16) of -% REF [1]) -% * |clustersTau|: $\tau_k$ in Eq.(16) of REF [1] -% * |clustersDeg|: posterior degrees of freedom ($\nu_k$ in Eq.(16) -% of REF [1]) -% * |clustersSigma|: posterior scale matrix ($S_k$ in Eq.(16) of -% REF [1]) -% * |logDetClustersSigma|: log-determinant of $S_k$ -% * |dcmMean|: posterior mean of DCM parameters ($\mu_n$ in -% Eq.(19) of REF [1]) -% * |dcmSigma|: posterior covariance of hemodynamic parameters -% ($\Sigma_n$ in Eq.(19) of REF [1]) -% * |logDetPostDcmSigma|: log-determinant of $\Sigma_n$ -% * |noiseInvScale|: posterior inverse scale of observation noise -% ($b_{n,r}$ in Eq.(21) of REF [1]) -% * |noiseShape|: posterior shape parameter of observation noise -% ($a_{n,r}$ in Eq.(21) of REF [1]) -% * |meanNoisePrecision|: posterior mean of precision of observation noise -% ($\lambda_{n,r}$ in Eq.(23) of REF [1]) -% * |modifiedSumSqrErr|: $b'_{n,r}$ in Eq.(22) of REF [1] -% - +% The negative free energy can be used to compare different models. For example, +% we can invert the HUGE model with three clusters, see how the negative free +% energy changes and calculate the Bayes factor between the two and three cluster +% models. %% -% The negative free energy after convergence is stored in -% |DcmResults.freeEnergy| -disp(['Negative free energy is ' num2str(DcmResults.freeEnergy) ... - ' after ' num2str(DcmResults.nIterationsActual) ' iterations.']) - +% invert the HUGE model with 3 clusters +% (this may take some time) +hugeSimK3 = hugeSim.estimate('K', 3); + +% compare negative free energy between 2 and 3 cluster models +fprintf('Difference in negative free energy: %f.\n', hugeSim.posterior.nfe - hugeSimK3.posterior.nfe) + +% calculate Bayes factor +BF_23 = exp(hugeSim.posterior.nfe - hugeSimK3.posterior.nfe); +fprintf('Bayes factor: %e.\n', BF_23) +%% +% Note that options set via name-value pair arguments are persistent across +% the lifetime of the object. For example, the object remembers that the option +% _verbose_ has been set to _true_ in the previous section, and keeps generating +% command line outputs. +% +% The difference in the negative free energy shows that the model prefers +% the two-cluster solution, which is what we would expect, given that the generating +% model contained two subgroups. +%% Importing Data +% In the previous example, we used the HUGE object itself to generate data. +% It is of course possible to import (experimental) fMRI time series into a HUGE +% object. The most convenient method is to use either the _tapas_Huge_ or the +% _estimate_ commands with the appropriate name-value pair argument. %% -% We can plot the results using the command -tapas_huge_plot(DCMr3b,DcmResults); +% select the first 10 subjects (bottom-up group) from the previous example +listDcmsG1 = listDcms(1:10); + +% create a new HUGE object and import the data from these subjects +hugeEb = tapas_Huge('Tag', 'Group 1 - bottom-up', 'Dcm', listDcmsG1); + +%% +% Note that the data to be imported has to be stored in a cell array where +% each cell holds a DCM struct in SPM's DCM format. This allows for convenient +% transfer of data from SPM to HUGE. +%% Empirical Bayes +% Previously, we demonstrated how to use the HUGE toolbox to stratify a heterogeneous +% cohort. However, HUGE can also be used to perform empirical Bayesian analysis; +% i.e., to invert the DCMs while learning the prior over DCM parameters from the +% data. This is accomplished by setting the number of clusters to one. %% -% The assignment probabilities shown in the upper left panel of the above -% figure indicates that all but three subjects could be correctly -% classified. +% empirical Bayes +% (this may take some time) +hugeEb = hugeEb.estimate('K', 1, 'PriorClusterVariance', 0.02, 'Verbose', true); + +% plot result +plot(hugeEb); +%% +% Here, we also used a name-value pair argument to change the prior cluster +% variance. Note also that the _plot_ method generates different graphs for empirical +% Bayesian analysis than for stratification. The first graph shows a boxplot of +% the posterior mean DCM parameter estimates of each individual subject. The second +% plot shows the population-level posterior mean and standard deviation. +% +% In this case the second graph reveals that for example, the connection +% from the first to the second region in the current cohort is stronger than the +% standard prior value of zero. % -% The panel in the top right shows the posterior cluster mean and 95% -% marginal credible intervals. -% -% The lower panel shows the simulated BOLD measurement (black) and 25 -% samples from the posterior over (noise free) BOLD signals. -% -% - - -%% Analyzing the second dataset -% Performing HUGE analysis on the second dataset proceeds analogously to -% the first dataset. Differences in the DCM's structure are handled -% automatically. % -disp(DCMr2l{1}) - -% choose number of clusters for the HUGE model -K = 2; - -% generate priors -priors = tapas_huge_build_prior(DCMr2l); -disp(priors) - -% suppress command line output -verbose = false; - -% randomize starting values -randomize = true; - -% analysis -currentTimer = tic; -[DcmResultsR2l] = tapas_huge_invert(DCMr2l, K, priors, verbose, randomize); -toc(currentTimer) - -disp(['Negative free energy is ' num2str(DcmResultsR2l.freeEnergy) ... - ' after ' num2str(DcmResultsR2l.nIterationsActual) ' iterations.']) - -% plot the results -tapas_huge_plot(DCMr2l,DcmResultsR2l); - - -%% Empirical Bayes analysis -% To perform empirical Bayes analysis, set |K| to |1|. Here, we perform -% empirical Bayesian analysis on the second dataset. For more information -% on empirical Bayes, see Fig.5 in REF [1]. -% - -% set K to 1 for empirical Bayes -K = 1; - -% generate priors -priors = tapas_huge_build_prior(DCMr2l); - -% suppress command line output -verbose = false; - -% randomize starting values -randomize = true; - -% analysis -currentTimer = tic; -[DcmResultsEb] = tapas_huge_invert(DCMr2l, K, priors, verbose, randomize); -toc(currentTimer) - -disp(['Negative free energy is ' num2str(DcmResultsEb.freeEnergy) ... - ' after ' num2str(DcmResultsEb.nIterationsActual) ' iterations.']) - -% plot the results -tapas_huge_plot(DCMr2l,DcmResultsEb); -%% -% The above figure summarizes the result of empirical Bayes. The top panel -% shows the posterior group-level mean with 95% marginal credible -% intervals. -% -% The lower panel shows a boxplot of first-level (i.e.: subject-level) -% MAP estimates of the DCM parameters. -% - - -%% References: -% -% [1] Yao Y, Raman SS, Schiek M, Leff A, Frssle S, Stephan KE (2018). -% Variational Bayesian Inversion for Hierarchical Unsupervised -% Generative Embedding (HUGE). NeuroImage, 179: 604-619 -% +% Author: Yu Yao (yao@biomed.ee.ethz.ch) +% Copyright (C) 2019 Translational Neuromodeling Unit +% Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of TAPAS, which is released under the terms of the GNU +% General Public Licence (GPL), version 3. For further details, see +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. % +% This software is intended for research only. Do not use for clinical +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% + diff --git a/huge/tapas_huge_demo.mlx b/huge/tapas_huge_demo.mlx new file mode 100644 index 00000000..25f84025 Binary files /dev/null and b/huge/tapas_huge_demo.mlx differ diff --git a/huge/tapas_huge_export_spm.m b/huge/tapas_huge_export_spm.m deleted file mode 100644 index e1a3d091..00000000 --- a/huge/tapas_huge_export_spm.m +++ /dev/null @@ -1,80 +0,0 @@ -%% [ DcmSpm ] = tapas_huge_export_spm( DcmInfo ) -% -% Convert DCM from DcmInfo format to SPM format. -% -% INPUT: -% DcmInfo - DCM in DcmInfo format -% -% OUTPUT: -% DcmSpm - cell array of DCM structs in SPM format -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ DcmSpm ] = tapas_huge_export_spm( DcmInfo ) -%% assemble DcmInfo -template = struct(); -template.U = struct(); % place holder -template.Y = struct(); % place holder - -N = DcmInfo.nSubjects; % #subjects -R = DcmInfo.nStates; % #regions -template.n = R; - -template.v = 0; % place holder - -template.TE = DcmInfo.hemParam.echoTime; - -template.options = struct(); -template.options.nonlinear = DcmInfo.dcmTypeD; -template.options.two_state = false; -template.options.stochastic = false; -template.options.nograph = true; - -% DCM connection indicators -template.a = DcmInfo.adjacencyA; -template.c = DcmInfo.adjacencyC; -template.b = DcmInfo.adjacencyB; -if(any(DcmInfo.dcmTypeD)) - template.d = DcmInfo.adjacencyD; -else - template.d = zeros(R,R,0); -end - -% fill in data for individial subjects -DcmSpm = repmat({template},N,1); -for n = 1:N - DcmSpm{n}.U.dt = DcmInfo.timeStep(n); - DcmSpm{n}.U.u = DcmInfo.listU{n}; - DcmSpm{n}.Y.dt = DcmInfo.trSeconds(n); - DcmSpm{n}.Y.y = reshape(DcmInfo.listBoldResponse{n}(:),[],R); - DcmSpm{n}.Y.name = repmat({''},R,1); - DcmSpm{n}.v = size(DcmSpm{n}.Y.y,1); - DcmSpm{n}.Y.X0 = ones(DcmSpm{n}.v,1); - DcmSpm{n}.Y.secs = DcmSpm{n}.v*DcmSpm{n}.Y.dt; - DcmSpm{n}.Y.Q = cell(1,R); - for r = 1:R - base = zeros(R); - base(r,r) = 1; - tile = speye(DcmSpm{n}.v); - DcmSpm{n}.Y.Q{r} = kron(base,tile); - end -end - - - - -end diff --git a/huge/tapas_huge_generate_examples.m b/huge/tapas_huge_generate_examples.m deleted file mode 100644 index f4c3c6e2..00000000 --- a/huge/tapas_huge_generate_examples.m +++ /dev/null @@ -1,126 +0,0 @@ -%% [ ] = tapas_huge_generate_examples( tname ) -% -% generate two synthetic example datasets based on a 2-region and a -% 3-region DCM and save them in SPM format. -% INPUT: -% tname - filename for saving datasets -% -% REFERENCE: -% [1] Yao Y, Raman SS, Schiek M, Leff A, Frssle S, Stephan KE (2018). -% Variational Bayesian Inversion for Hierarchical Unsupervised -% Generative Embedding (HUGE). NeuroImage, 179: 604-619 -% -% https://doi.org/10.1016/j.neuroimage.2018.06.073 -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ ] = tapas_huge_generate_examples( tname ) -%% generate data from a two-region linear DCM -optionsGen = struct(); -optionsGen.snr = 1; % signal-to-noise-ratio -optionsGen.N_k = [10 10]; % number of subjects per cluster -optionsGen.R = 2; % number of regions - -% cluster -optionsGen.mu_k.idx = [1,3,4,5,8]; -optionsGen.mu_k.value = [... - -0.4, 0.5,-0.6, 0.3 , 0.8;... - -0.6,-0.2,-0.4, 0.8 , 0.3]; -optionsGen.sigma_k = 0.1; - -% hemodynamics -optionsGen.mu_h = zeros(1,optionsGen.R*2+1); -optionsGen.sigma_h = 0.01; - -% hemodynamic parameters -optionsGen.hemParam.listHem = [0.6400 2 1]; -optionsGen.hemParam.scaleC = 16; -optionsGen.hemParam.echoTime = 0.0400; -optionsGen.hemParam.restingVenousVolume = 4; -optionsGen.hemParam.relaxationRateSlope = 25; -optionsGen.hemParam.frequencyOffset = 40.3000; -optionsGen.hemParam.oxygenExtractionFraction = 0.4000; -optionsGen.hemParam.rho = 4.3000; -optionsGen.hemParam.gamma = 0.3200; -optionsGen.hemParam.alphainv = 3.1250; -optionsGen.hemParam.oxygenExtractionFraction2 = 0.3200; - -% input -nU = 6000; -idx1 = randi( nU - 400, 1, 75); -idx1 = bsxfun(@plus, idx1 + 200, (0:3)'); -optionsGen.input.u = zeros(nU, 2); -optionsGen.input.u(idx1(:), 1) = 1; -idx2 = bsxfun(@plus, 300:600:nU-301, (1:300)'); -optionsGen.input.u(idx2(:), 2) = 1; - -optionsGen.input.trSteps = 20; -optionsGen.input.trSeconds = 2; - -DcmInfo = tapas_huge_simulate(optionsGen); -DCMr2l = tapas_huge_export_spm(DcmInfo); %#ok - - -%% generate data from a three-region bilinear DCM -optionsGen = struct(); -optionsGen.snr = 1; % signal-to-noise-ratio -optionsGen.N_k = [40 20 20]; % number of subjects per cluster -optionsGen.R = 3; % number of regions - -% cluster -optionsGen.mu_k.idx = [1,2,3,4,5,6,9,10,27]; -optionsGen.mu_k.value = [... - -0.7, 0.2,-0.1,-0.2 ,-0.6, 0.3,-0.4,0.3, 0.3;... - -0.7, 0.1, 0.3, 0.1 ,-0.4,-0.1,-0.6,0.6, 0.1;... - -0.7,-0.2, 0.3, 0.25,-0.4,-0.1,-0.6,0.6,-0.2]; -optionsGen.sigma_k = 0.1; - -% hemodynamics -optionsGen.mu_h = zeros(1,optionsGen.R*2+1); -optionsGen.sigma_h = 0.01; - -% hemodynamic parameters -optionsGen.hemParam.listHem = [0.6400 2 1]; -optionsGen.hemParam.scaleC = 16; -optionsGen.hemParam.echoTime = 0.0400; -optionsGen.hemParam.restingVenousVolume = 4; -optionsGen.hemParam.relaxationRateSlope = 25; -optionsGen.hemParam.frequencyOffset = 40.3000; -optionsGen.hemParam.oxygenExtractionFraction = 0.4000; -optionsGen.hemParam.rho = 4.3000; -optionsGen.hemParam.gamma = 0.3200; -optionsGen.hemParam.alphainv = 3.1250; -optionsGen.hemParam.oxygenExtractionFraction2 = 0.3200; - -% input -optionsGen.input.u = double([... - reshape(repmat((1:2^9<2^8)'&(mod(1:2^9,2^5)==0)',1,2^3),2^12,1),... - reshape(repmat((1:2^10<2^8)',1,2^2),2^12,1)]); -optionsGen.input.u = circshift(optionsGen.input.u,117,1); -optionsGen.input.trSteps = 16; -optionsGen.input.trSeconds = 2; - -DcmInfo = tapas_huge_simulate(optionsGen); -DCMr3b = tapas_huge_export_spm(DcmInfo); %#ok - - -%% save DCM in SPM format -save(tname,'DCMr2l','DCMr3b'); - - -end - diff --git a/huge/tapas_huge_import_spm.m b/huge/tapas_huge_import_spm.m deleted file mode 100644 index da0b5a86..00000000 --- a/huge/tapas_huge_import_spm.m +++ /dev/null @@ -1,159 +0,0 @@ -%% [ DcmInfo ] = tapas_huge_import_spm( DcmSpm ) -% -% Convert DCMs from SPM format to DcmInfo format. -% -% INPUT: -% DcmSpm - (cell) array of DCM structs in SPM format. All DCM in -% DcmSpm must have same structure (i.e.: same number of -% regions and same connections.) They may have different -% number of scans. -% -% OUTPUT: -% DcmInfo - DCM in DcmInfo format -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ DcmInfo ] = tapas_huge_import_spm( DcmSpm ) -%% check input -if isvector(DcmSpm)&&isstruct(DcmSpm) - % convert to cell array - try - DcmSpm = {DcmSpm(:).DCM}'; - catch - DcmSpm = num2cell(DcmSpm); - end -else - assert(iscell(DcmSpm),'TAPAS:HUGE:InputFormat',... - 'Input must be cell array or array of structs'); -end - - -%% assemble DcmInfo -DcmInfo = struct(); - -N = length(DcmSpm); % #subjects -DcmInfo.nSubjects = N; -R = DcmSpm{1}.n; % #regions -DcmInfo.nStates = R; -L = size(DcmSpm{1}.U.u,2); % #inputs -DcmInfo.nInputs = L; - -% DCM connection indicators -DcmInfo.adjacencyA = logical(DcmSpm{1}.a); -DcmInfo.adjacencyC = logical(DcmSpm{1}.c); -DcmInfo.adjacencyB = logical(DcmSpm{1}.b); -if(~isempty(DcmSpm{1}.d)) - DcmInfo.adjacencyD = DcmSpm{1}.d; -else - DcmInfo.adjacencyD = false(R,R,R); -end -DcmInfo.dcmTypeB = any(DcmInfo.adjacencyB(:)); -DcmInfo.dcmTypeD = any(DcmInfo.adjacencyD(:)); - -adjacencyList = [(DcmInfo.adjacencyA(:))' (DcmInfo.adjacencyC(:))']; -adjacencyList = [adjacencyList (DcmInfo.adjacencyB(:))']; -adjacencyList = [adjacencyList (DcmInfo.adjacencyD(:))']; -connectionIndicator = find(adjacencyList ~= 0); -noConnectionIndicator = find(adjacencyList == 0); -DcmInfo.noConnectionIndicator = noConnectionIndicator; -DcmInfo.connectionIndicator = connectionIndicator; -DcmInfo.nNoConnections = length(noConnectionIndicator); -DcmInfo.nConnections = length(connectionIndicator); - -nParameters = (R^2)*(1 + L + R) + L*R + 3*R; % A, B, D, C, hem -DcmInfo.nParameters = nParameters; - - - -% hemodynamic model -DcmInfo.hemParam = struct(); -% decay (kappa), transit time (tau), ratio intra/extra (epsilon) -DcmInfo.hemParam.listHem = [0.64,2,1]; % SPM default values -DcmInfo.hemParam.scaleC = DcmSpm{1}.Y.dt/DcmSpm{1}.U.dt; -if isfield(DcmSpm{1},'TE') - DcmInfo.hemParam.echoTime = DcmSpm{1}.TE; -else - DcmInfo.hemParam.echoTime = 0.04; % default value hard-coded into SPM -end -% default value hard-coded into SPM -DcmInfo.hemParam.restingVenousVolume = 4; -DcmInfo.hemParam.relaxationRateSlope = 25; -DcmInfo.hemParam.frequencyOffset = 40.3; -DcmInfo.hemParam.oxygenExtractionFraction = .4; -DcmInfo.hemParam.rho = 4.3; -DcmInfo.hemParam.gamma = .32; -DcmInfo.hemParam.alphainv = 1/.32; -DcmInfo.hemParam.oxygenExtractionFraction2 = .4; - -% inputs and data -DcmInfo.trSeconds = zeros(1,N); -DcmInfo.trSteps = zeros(1,N); -DcmInfo.timeStep = zeros(1,N); -DcmInfo.nTime = zeros(1,N); - -DcmInfo.listU = cell(1,N); -DcmInfo.listBoldResponse = cell(1,N); -DcmInfo.listResponseTimeIndices = cell(1,N); - - -DcmInfo.listParameters = cell(N,1); -DcmInfo.trueLabels = NaN(1,N); - -for n = 1:N - - assert(DcmSpm{n}.n==R,'TAPAS:HUGE:InconsistentParameters',... - 'Number of regions inconsistent between subjects'); - assert(size(DcmSpm{n}.U.u,2)==L,'TAPAS:HUGE:InconsistentParameters',... - 'Number of inputs inconsistent between subjects'); - - assert(~any(DcmSpm{n}.a(~DcmInfo.adjacencyA)),... - 'TAPAS:HUGE:InconsistentParameters',... - 'A matrix inconsistent between subjects'); - assert(~any(DcmSpm{n}.b(~DcmInfo.adjacencyB)),... - 'TAPAS:HUGE:InconsistentParameters',... - 'B matrix inconsistent between subjects'); - assert(~any(DcmSpm{n}.c(~DcmInfo.adjacencyC)),... - 'TAPAS:HUGE:InconsistentParameters',... - 'C matrix inconsistent between subjects'); - if(~isempty(DcmSpm{n}.d)) - adjacencyD = DcmSpm{n}.d; - else - adjacencyD = false(R,R,R); - end - assert(~any(adjacencyD(~DcmInfo.adjacencyD)),... - 'TAPAS:HUGE:InconsistentParameters',... - 'D matrix inconsistent between subjects'); - - % inputs - DcmInfo.listU{n} = full(DcmSpm{n}.U.u); - DcmInfo.nTime(n) = size(DcmInfo.listU{n},1); - DcmInfo.timeStep(n) = DcmSpm{n}.U.dt; - - % BOLD time series - [nScans,nRegions] = size(DcmSpm{n}.Y.y); - assert(nRegions==R,'TAPAS:HUGE:DataSize',... - 'Size of data vector inconsistent with number of regions'); - DcmInfo.listBoldResponse{n} = DcmSpm{n}.Y.y(:); - DcmInfo.trSeconds(n) = DcmSpm{n}.Y.dt; - DcmInfo.trSteps(n) = DcmInfo.trSeconds(n)/DcmInfo.timeStep(n); - DcmInfo.listResponseTimeIndices{n} = ... - DcmInfo.trSteps(n):DcmInfo.trSteps(n):DcmInfo.trSteps(n)*nScans; - -end - - -end diff --git a/huge/tapas_huge_int_euler.c b/huge/tapas_huge_int_euler.c index 162f6899..d4c1fcbd 100644 --- a/huge/tapas_huge_int_euler.c +++ b/huge/tapas_huge_int_euler.c @@ -41,14 +41,19 @@ % % This file is part of TAPAS, which is released under the terms of the GNU % General Public Licence (GPL), version 3. For further details, see -% . +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. % % This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: % https://github.com/translationalneuromodeling/tapas/issues -% +% + %%*/ #include diff --git a/huge/tapas_huge_inv_vb.m b/huge/tapas_huge_inv_vb.m deleted file mode 100644 index 945e9057..00000000 --- a/huge/tapas_huge_inv_vb.m +++ /dev/null @@ -1,752 +0,0 @@ -%% [ DcmResults ] = tapas_huge_inv_vb( DcmInfo, DcmResults ) -% -% Inversion of hierarchical unsupervised generative embedding (HUGE) for -% DCM for fMRI using variational Bayes. -% -% INPUT: -% DcmInfo - struct containing DCM model specification and BOLD -% time series in DcmInfo format -% (see tapas_huge_simulate.m for an example) -% DcmResults - struct containing prior specification and settings. -% -% OUTPUT: -% DcmResults - struct used for storing the results from VB -% -% REFERENCE: -% [1] Yao Y, Raman SS, Schiek M, Leff A, Frssle S, Stephan KE (2018). -% Variational Bayesian Inversion for Hierarchical Unsupervised -% Generative Embedding (HUGE). NeuroImage, 179: 604-619 -% -% https://doi.org/10.1016/j.neuroimage.2018.06.073 -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ DcmResults ] = tapas_huge_inv_vb( DcmInfo, DcmResults ) -% -KMEANS_REP = 10; - - -%% extract DCM information -nStates = DcmInfo.nStates; %%% R -nSubjects = DcmInfo.nSubjects; %%% N - -listBoldResponse = DcmInfo.listBoldResponse; %%% y - -nMeasurements = zeros(nSubjects,1); -nMeasurementsPerState = zeros(nSubjects,1); - -for iSubject = 1:nSubjects - % %%% dim(y) - nMeasurements(iSubject) = length(listBoldResponse{iSubject}); - nMeasurementsPerState(iSubject) = nMeasurements(iSubject)/nStates; %%% q_r - -end - -[nParameters,idxParamsInf,idxSelfCon] = tapas_huge_count_params(DcmInfo); -% % % % % nDcmParametersAll = nParameters(1,3); -nDcmParamsInfCon = nParameters(2,1); %%% d_c -% % % % % nDcmParamsInfHem = nParameters(2,2); %%% d_h -nDcmParamsInfAll = nParameters(2,3); %%% d - -% measurement function %%% g(theta) -fnGenerateResponse = @tapas_huge_predict; - -counts.nSubjects = nSubjects; -counts.nParameters = nParameters; -counts.nMeasurements = nMeasurements; -counts.nMeasurementsPerState = nMeasurementsPerState; - - -%% parameters for variational Bayes: -% maximum number of iterations -nIterations = DcmResults.nIterations; - -% free energy threshold -epsEnergy = DcmResults.epsEnergy; - -% parameters for jacobian calculation -paramsJacobian.dimOut = nMeasurements; - -% method for calculating the jacobian matrix -fnJacobian = DcmResults.fnJacobian; - -% diagonal constant for numerical stability -diagConst = DcmResults.diagConst; -diagConst = diagConst*eye(nDcmParamsInfAll); - -% keep history of important variables -bKeepTrace = DcmResults.bKeepTrace; - -% keep history of response related variables -% note: has no effect if bKeepTrace is false -bKeepResp = DcmResults.bKeepResp; - -% switch off command line output -bVerbose = DcmResults.bVerbose; - -% update schedule -schedule = DcmResults.schedule; - -% number of clusters -nClusters = DcmResults.maxClusters; %%% K -if nClusters > nSubjects - nClusters = nSubjects; - DcmResults.maxClusters = nClusters; - warning('TAPAS:HUGE:NumberOfCluster',... - ['The number of clusters exceeds the number of subjects and '... - 'has been decreased for efficiency.']); -end -counts.nClusters = nClusters; - -% set default parameters to zero -dcmParametersDefault = zeros(1,DcmInfo.nParameters); - - -%% priors -priors = DcmResults.priors; -% component weights -if isscalar(priors.alpha) - priors.alpha = repmat(priors.alpha,nClusters,1); -end - -priors.hemPrecision = inv(priors.hemSigma); - -% prior noise level on fmri data -if isscalar(priors.noiseShape) %%% a_r: shape - priors.noiseShape = repmat(priors.noiseShape,1,nStates); -end -if isscalar(priors.noiseInvScale) %%% b_r: inverse scale - priors.noiseInvScale = repmat(priors.noiseInvScale,1,nStates); -end - - -%% initial values -% initialize posterior covariance of DCM parameters -bInit = isfield(DcmResults,'init'); - -if bInit&&isfield(DcmResults.init,'covFactor') - priors.clustersSigma = DcmResults.init.covFactor*priors.clustersSigma; - bCovFactor = true; -else - bCovFactor = false; -end - - -if bInit && isfield(DcmResults.init,'dcmSigma') - posterior.dcmSigma = DcmResults.init.dcmSigma; -else - posterior.dcmSigma(1:nDcmParamsInfCon,1:nDcmParamsInfCon) = ... - priors.clustersSigma; - posterior.dcmSigma(nDcmParamsInfCon+1:nDcmParamsInfAll,... - nDcmParamsInfCon+1:nDcmParamsInfAll) = priors.hemSigma; - posterior.dcmSigma = repmat(posterior.dcmSigma,[1,1,nSubjects]); -end -posterior.logDetPostDcmSigma = zeros(nSubjects,1); -for iSubject = 1:nSubjects - posterior.logDetPostDcmSigma(iSubject) = ... - tapas_util_logdet(posterior.dcmSigma(:,:,iSubject)); -end - -% initialize posterior mean of DCM parameters -if bInit && isfield(DcmResults.init,'dcmMean') - posterior.dcmMean = DcmResults.init.dcmMean; -else - posterior.dcmMean = repmat([priors.clustersMean,priors.hemMean],... - nSubjects,1); -end -% precalculate Jacobian, error term and sum of squared error -respError = cell(nSubjects,1); %%% epsilon -respJacobian = cell(nSubjects,1); %%% G -posterior.modifiedSumSqrErr = zeros(nSubjects,nStates); %%% hat(b)_nr - -for iSubject = 1:nSubjects - respCurrent = fnGenerateResponse(posterior.dcmMean(iSubject,:).',... - dcmParametersDefault,... - idxParamsInf,... - idxSelfCon,... - DcmInfo, iSubject); - respError{iSubject} = listBoldResponse{iSubject} - respCurrent; - paramsJacobian.dimOut = nMeasurements(iSubject); - respJacobian{iSubject} = fnJacobian(... - fnGenerateResponse,... - posterior.dcmMean(iSubject,:).',... %%% mu_n and current m_n - paramsJacobian, respCurrent,... - dcmParametersDefault, idxParamsInf, ... - idxSelfCon, DcmInfo, iSubject); - - % check if initial values cause divergence in response generating - % function - if fnc_check_response(respError{iSubject},... - respJacobian{iSubject}) - error('TAPAS:HUGE:badInit',... - 'Initial values of DCM Parameters cause divergence.'); - end - %%% epsilon^T*Q_r*epsilon + tr(Q_r*G_n*Sigma_n*G_n^T) - posterior.modifiedSumSqrErr(iSubject,:) = ... - sum(reshape(... - respError{iSubject}.^2 + ... - sum(... - (respJacobian{iSubject}*... - posterior.dcmSigma(:,:,iSubject)).*... - respJacobian{iSubject},2),... - nMeasurementsPerState(iSubject),nStates )); -end - - -% initialize posterior mean and covariance of clusters -if bInit && isfield(DcmResults.init,'clustersMean') - posterior.clustersMean = DcmResults.init.clustersMean; -else - posterior.clustersMean = repmat(priors.clustersMean,nClusters,1); -end - -if bInit && isfield(DcmResults.init,'clustersSigma') - posterior.clustersSigma = DcmResults.init.clustersSigma; -else - posterior.clustersSigma = repmat(priors.clustersSigma,[1,1,nClusters]); -end -posterior.logDetClustersSigma = zeros(nClusters,1); %%% log(det(bar(Sigma)_k)) -invClustersSigma = zeros(nDcmParamsInfCon,nDcmParamsInfCon,nClusters); %%% bar(Sigma)_k^-1 -for iCluster = 1:nClusters - posterior.logDetClustersSigma(iCluster) = ... - tapas_util_logdet(posterior.clustersSigma(:,:,iCluster)); - invClustersSigma(:,:,iCluster) = ... - inv(posterior.clustersSigma(:,:,iCluster)); -end - -posterior.noiseInvScale = repmat(priors.noiseInvScale,nSubjects,1); -% note: with the current parameterization 'posterior.noiseShape' is precomputable -posterior.noiseShape = ... %%% a_nr - repmat(priors.noiseShape + nMeasurementsPerState(iSubject)/2,... - nSubjects,1); -% mean precision -posterior.meanNoisePrecision = ... - posterior.noiseShape./posterior.noiseInvScale; - - -%% prelocating memory and initializing auxiliary variables -weightedSumDcmMean = zeros(nClusters,nDcmParamsInfCon); %%% mu_ck - -% initialize soft assign to one cluster only (empirical Bayes) -if bInit && isfield(DcmResults.init,'logSoftAssign') - logSoftAssign = DcmResults.init.logSoftAssign; -else - logSoftAssign = ones(nSubjects,nClusters); %%% log(q_nk) - logSoftAssign(:,1) = 100; - DcmResults.init.logSoftAssign = logSoftAssign; -end -posterior.softAssign = fnc_exp_norm(logSoftAssign); -nAssign = sum(posterior.softAssign,1).' + realmin; -posterior.alpha = priors.alpha + nAssign; %%% alpha -posterior.clustersTau = priors.clustersTau + nAssign; -posterior.clustersDeg = priors.clustersDeg + nAssign; - -partialDcmPrec = zeros(nDcmParamsInfAll); %%% Lambda'_n -partialDcmPrec(nDcmParamsInfCon+1:end,nDcmParamsInfCon+1:end) = ... - priors.hemPrecision; % note: the lower submatrix is constant - -partialDcmMean = zeros(1,nDcmParamsInfAll); %%% mu'_n -partialDcmMean(nDcmParamsInfCon+1:end) = priors.hemMean/priors.hemSigma; - -[freeEnergy,feAux] = tapas_huge_nfe(counts,priors,posterior); -feCurrent = freeEnergy; -dF = -1; - - - -%% history - -nItSubject = zeros(nSubjects,1); - -% keep history of important variables -if bKeepTrace - histClustersMean = cell(nIterations,1); - histClustersSigma = cell(nIterations,1); - histClustersTau = cell(nIterations,1); - histClustersDeg = cell(nIterations,1); - - histDcmMean = cell(nIterations,1); - histDcmSigma = cell(nIterations,1); - - histNoiseInvScale = cell(nIterations,1); - - histposterior.softAssign = cell(nIterations,1); - histPartialDcmMean = cell(nIterations,1); - histPartialDcmPrec = cell(nIterations,1); - if bKeepResp - histRespError = cell(nIterations,1); - histRespJacobian = cell(nIterations,1); - end - - histFreeEnergy = cell(nIterations,1); - histFeParts = cell(nIterations,1); - - timeSinceStart = cell(nIterations,1); - tic; -end - -histFe = zeros(nIterations,1); -itSatDcm = nIterations + 1; -itSatClusters = nIterations + 1; - - -% ------------------------------------------- -%% Variational Updates - Main loop -% ------------------------------------------- - -for iIteration = 1:nIterations - - - - - if schedule.itKmeans&&(iIteration == itSatClusters + schedule.itKmeans) - - [kmeansIdx,posterior.clustersMean] = ... - kmeans(posterior.dcmMean(:,1:nDcmParamsInfCon),... - nClusters,'Replicates',KMEANS_REP); - posterior.softAssign = zeros(nSubjects,nClusters); - posterior.softAssign(sub2ind([nSubjects,nClusters],... - (1:nSubjects)',kmeansIdx)) = 1; - nAssign = sum(posterior.softAssign,1).' + realmin; - posterior = fnc_set_cluster_cov(priors,posterior,kmeansIdx,nAssign); - - - % debug info - DcmResults.kmeans.kmIdx = kmeansIdx; - DcmResults.kmeans.dcmMean = posterior.dcmMean; - DcmResults.kmeans.iRelease = iIteration; - DcmResults.kmeans.dF = dF; -% DcmResults.kmeans.kmProb = KMEANS_PROB; - DcmResults.kmeans.kmRep = KMEANS_REP; - end - - - - if iIteration >= itSatClusters + schedule.itAssignment - %%% q_nk - % update soft assignments - postAlphaDigamma = psi(0,posterior.alpha); % Psi(sum(alpha)) is a constant - for iCluster = 1:nClusters - diffMeanDcmCluster = bsxfun(@minus,... %%% mu_cn - mu_k - posterior.dcmMean(:,1:nDcmParamsInfCon),... - posterior.clustersMean(iCluster,:)); - % note: logSoftAssign is only correct up to a constant - logSoftAssign(:,iCluster) = ... %%% log(q_nk) - -posterior.logDetClustersSigma(iCluster)/2 ... - +sum(psi(0,... - (posterior.clustersDeg(iCluster)-... - (0:nDcmParamsInfCon-1))/2))/2 ... - -nDcmParamsInfCon/2/posterior.clustersTau(iCluster)... - -posterior.clustersDeg(iCluster)/2*... - reshape( sum( sum( ... - bsxfun( @times, ... - invClustersSigma(:,:,iCluster),... - posterior.dcmSigma(1:nDcmParamsInfCon,... - 1:nDcmParamsInfCon,:) ), ... - 1), 2), nSubjects,1) ... %%% -nu/2tr(bar(Sigma)_k^-1*Sigma_cn) - -posterior.clustersDeg(iCluster)/2*... - sum((diffMeanDcmCluster/... - posterior.clustersSigma(:,:,iCluster)).*... - diffMeanDcmCluster,2) ... - +postAlphaDigamma(iCluster); - end - posterior.softAssign = fnc_exp_norm(logSoftAssign); %%% q_nk - % effective number of samples assigned to each cluster - nAssign = sum(posterior.softAssign,1).' + realmin; %%% q_k - %%% pi: alpha - % update parameters for mixture weights - posterior.alpha = priors.alpha + nAssign; - - end - - - - if iIteration >= itSatDcm + schedule.itCluster - - - %%% bar(mu)_k, bar(Sigma)_k - % update parameters of inverse Wishart - posterior.clustersTau = priors.clustersTau + nAssign; - posterior.clustersDeg = priors.clustersDeg + nAssign; - - for iCluster = 1:nClusters - weightedSumDcmMean(iCluster,:) = ... - sum(bsxfun(@times,... - posterior.dcmMean(:,1:nDcmParamsInfCon),... - posterior.softAssign(:,iCluster)),... - 1); - posterior.clustersMean(iCluster,:) = ... %%% (q_k*mu_ck + tau_0*bar(mu)_0)/tau_k - (weightedSumDcmMean(iCluster,:) + ... - priors.clustersTau*priors.clustersMean)/... - posterior.clustersTau(iCluster); - weightedSumDcmMean(iCluster,:) = ... - weightedSumDcmMean(iCluster,:)/nAssign(iCluster); - diffMeanDcmWeighted = bsxfun(@minus,... %%% mu_cn - mu_ck - posterior.dcmMean(:,1:nDcmParamsInfCon),... - weightedSumDcmMean(iCluster,:)); - posterior.clustersSigma(:,:,iCluster) = ... %%% bar(Sigma)_k = .. - +priors.clustersSigma ... %%% bar(Sigma)_0 + .. - +nAssign(iCluster)*priors.clustersTau/... - posterior.clustersTau(iCluster)*... %%% q_k*tau_0/tau_k*(mu_ck - bar(mu)_0)*(mu_ck - bar(mu)_0)^T + .. - (weightedSumDcmMean(iCluster,:) - priors.clustersMean).'* ... - (weightedSumDcmMean(iCluster,:) - priors.clustersMean) ... - +diffMeanDcmWeighted.'*... - bsxfun(@times,diffMeanDcmWeighted,... - posterior.softAssign(:,iCluster)) ... %%% sum(q_nk*(mu_cn - mu_ck)*(mu_cn - mu_ck)^T) + ... - +sum(... %%% Sigma_ck - bsxfun(... - @times,... - posterior.dcmSigma(1:nDcmParamsInfCon,1:... - nDcmParamsInfCon,:),... - reshape(posterior.softAssign(:,iCluster),1,1,nSubjects)... - ),... - 3); %%% Sigma_ck - posterior.logDetClustersSigma(iCluster) = ... - tapas_util_logdet(posterior.clustersSigma(:,:,iCluster)); - invClustersSigma(:,:,iCluster) = ... - inv(posterior.clustersSigma(:,:,iCluster)); %%% inv(bar(Sigma)_k) %%% TODO check rcond - end - end - - - % cache DCM parameters - dcmMeanBkp = posterior.dcmMean; - dcmSigmaBkp = posterior.dcmSigma; - logDetDcmSigmaBkp = posterior.logDetPostDcmSigma; - - respErrorBkp = respError; - respJacBkp = respJacobian; - sumSqErrBkp = posterior.modifiedSumSqrErr; - - - % DCM update - for iSubject = 1:nSubjects - - %%% mu_n, Sigma_n - % update parameters for distribution over DCM parameters - weightedClustersPrecision = bsxfun(... - @times,... %%% q_nk*nu_k*bar(Sigma)_k^-1 - invClustersSigma,... - reshape(transpose(... - posterior.softAssign(iSubject,:)).*... - posterior.clustersDeg,... - 1,1,nClusters)); - partialDcmPrec(1:nDcmParamsInfCon,1:nDcmParamsInfCon) = ... - sum(weightedClustersPrecision,3); %%% sum(q_nk*nu_k*bar(Sigma)_k^-1) - partialDcmMean(1:nDcmParamsInfCon) = zeros(1,nDcmParamsInfCon); - for iCluster = 1:nClusters - partialDcmMean(1:nDcmParamsInfCon) = ... - partialDcmMean(1:nDcmParamsInfCon) + ... - posterior.clustersMean(iCluster,:)*... - weightedClustersPrecision(:,:,iCluster); - end - expandedMeanNoisePrecision = ... - kron(posterior.meanNoisePrecision(iSubject,:),... - ones(1,nMeasurementsPerState(iSubject))); - - posterior.dcmSigma(:,:,iSubject) = ... %%% Sigma_n = ... - inv(... %%% (G^T*bar(Lambda)_n*G_n + ... - transpose(respJacobian{iSubject})*... %%% TODO check rcond - bsxfun(... - @times,... - respJacobian{iSubject},... - expandedMeanNoisePrecision.') + ... - partialDcmPrec + ... %%% Lambda'_n)^-1 - diagConst); % add a small constant to the diagonal elements - - posterior.logDetPostDcmSigma(iSubject) = ... - tapas_util_logdet(posterior.dcmSigma(:,:,iSubject)); - - posterior.dcmMean(iSubject,:) = ... - (((respError{iSubject}.' + ... - posterior.dcmMean(iSubject,:)*respJacobian{iSubject}.').*... - expandedMeanNoisePrecision)*respJacobian{iSubject} + ... - partialDcmMean)*... - posterior.dcmSigma(:,:,iSubject); - - % calculate auxiliary variables for next iteration - respCurrent = fnGenerateResponse(posterior.dcmMean(iSubject,:).',... - dcmParametersDefault,... - idxParamsInf,idxSelfCon, ... - DcmInfo, iSubject); - respError{iSubject} = listBoldResponse{iSubject} - respCurrent; - paramsJacobian.dimOut = nMeasurements(iSubject); - respJacobian{iSubject} = fnJacobian(... - fnGenerateResponse,... - posterior.dcmMean(iSubject,:).',... %%% mu_n and current m_n - paramsJacobian,... - respCurrent,... % current function value - dcmParametersDefault, idxParamsInf, ... - idxSelfCon, DcmInfo, iSubject); - - posterior.modifiedSumSqrErr(iSubject,:) = ... %%% epsilon^T*Q_r*epsilon + tr(Q_r*G_n*Sigma_n*G_n^T) - sum(reshape(... - respError{iSubject}.^2 + ... - sum(... - (respJacobian{iSubject}*... - posterior.dcmSigma(:,:,iSubject)).*... - respJacobian{iSubject},2),... - nMeasurementsPerState(iSubject),nStates )); - - % check free energy - feBkp = feCurrent; - feAuxBkp = feAux; - [feCurrent,feAux] = tapas_huge_nfe(counts,priors,posterior,... - feAux,iSubject); - - % if free energy decreases or response/jacobian contains NaNs - % cancel update and restore old values - if feCurrent < feBkp || ... - fnc_check_response(respError{iSubject},... - respJacobian{iSubject}) - - posterior.dcmMean(iSubject,:) = dcmMeanBkp(iSubject,:); - posterior.dcmSigma(:,:,iSubject) = dcmSigmaBkp(:,:,iSubject); - posterior.logDetPostDcmSigma(iSubject) = ... - logDetDcmSigmaBkp(iSubject); - - respError{iSubject} = respErrorBkp{iSubject}; - respJacobian{iSubject} = respJacBkp{iSubject}; - posterior.modifiedSumSqrErr(iSubject,:) = ... - sumSqErrBkp(iSubject,:); - feCurrent = feBkp; - feAux = feAuxBkp; - - else - nItSubject(iSubject) = nItSubject(iSubject) + 1; - end - end - - % noise precision update - noiseInvScaleBkp = posterior.noiseInvScale; - meanNoisePrecisionBkp = posterior.meanNoisePrecision; - - posterior.noiseInvScale = bsxfun(@plus,... %%% b_nr - posterior.modifiedSumSqrErr/2,... - priors.noiseInvScale); - posterior.meanNoisePrecision = ... - posterior.noiseShape./posterior.noiseInvScale; %%% bar(lambda)_nr - - % check free energy - feBkp = feCurrent; - feCurrent = tapas_huge_nfe(counts,priors,posterior); - if feCurrent < feBkp - - posterior.noiseInvScale = noiseInvScaleBkp; - posterior.meanNoisePrecision = meanNoisePrecisionBkp; - feCurrent = feBkp; %#ok - - end - - - - - - - -%-------------------- - % save complete history of parameters and important auxiliary variables - if bKeepTrace - histClustersMean{iIteration} = posterior.clustersMean; - histClustersSigma{iIteration} = posterior.clustersSigma; - histClustersTau{iIteration} = posterior.clustersTau; - histClustersDeg{iIteration} = posterior.clustersDeg; - - histDcmMean{iIteration} = posterior.dcmMean; - histDcmSigma{iIteration} = posterior.dcmSigma; - - histNoiseInvScale{iIteration} = posterior.noiseInvScale; - - histposterior.softAssign{iIteration} = posterior.softAssign; - histPartialDcmMean{iIteration} = partialDcmMean; - histPartialDcmPrec{iIteration} = partialDcmPrec; - if bKeepResp - histRespError{iIteration} = respError; - histRespJacobian{iIteration} = respJacobian; - end - - histFreeEnergy{iIteration} = freeEnergy; - histFeParts{iIteration} = feParts; - timeSinceStart{iIteration} = toc; - end - - - -%------------------- - - [feCurrent,feAux] = tapas_huge_nfe(counts,priors,posterior); - % check stopping condition - dF = feCurrent - freeEnergy; - freeEnergy = feCurrent; - histFe(iIteration) = freeEnergy; - if bVerbose - display(['iteration ' num2str(iIteration) ... - ', dF: ' num2str(dF)]); - end - - if itSatDcm > nIterations - if dF < schedule.dfDcm - itSatDcm = iIteration; - end - elseif itSatClusters > nIterations - if dF < schedule.dfClusters - itSatClusters = iIteration; - if bCovFactor - [priors, posterior] = fnc_reset_cov(priors,posterior,... - DcmResults,DcmInfo,nDcmParamsInfCon,nDcmParamsInfAll); - end - end - else - if (iIteration>=itSatClusters+schedule.itReturn) && (dF < epsEnergy) - break; - end - end - - -end - - - -%------------------------ End: Main loop ------------------------------ -%% save results -DcmResults.posterior = posterior; - -% save other info: number of actual iterations etc -DcmResults.freeEnergy = freeEnergy; -DcmResults.nIterationsActual = iIteration; -DcmResults.nIterationsSubject = nItSubject; -DcmResults.inversionScheme = 'VB'; - -% save prediction and residual -predictedResponse = cell(nSubjects,1); -for iSubject = 1:nSubjects - predictedResponse{iSubject} = ... - fnGenerateResponse(posterior.dcmMean(iSubject,:).',... - dcmParametersDefault,... - idxParamsInf,... - idxSelfCon,... - DcmInfo, iSubject); -end -DcmResults.predictedResponse = predictedResponse; -DcmResults.residuals = respError; - -DcmResults.histFe = histFe(1:iIteration); -DcmResults.itSatDcm = itSatDcm; -DcmResults.itSatClusters = itSatClusters; - - - -% save history of important variables -if bKeepTrace - % note: 'DcmResults.debug.trace' will be an array of structs not a - % struct of arrays - DcmResults.debug.trace = struct(... - 'clustersMean',histClustersMean(1:iIteration), ... - 'clustersSigma',histClustersSigma(1:iIteration), ... - 'clustersTau',histClustersTau(1:iIteration), ... - 'clustersDeg',histClustersDeg(1:iIteration), ... - 'dcmMean',histDcmMean(1:iIteration), ... - 'dcmSigma',histDcmSigma(1:iIteration), ... - 'noiseInvScale',histNoiseInvScale(1:iIteration), ... - 'posterior.softAssign',histposterior.softAssign(1:iIteration), ... - 'partialDcmMean',histPartialDcmMean(1:iIteration), ... - 'partialDcmPrec',histPartialDcmPrec(1:iIteration), ... - 'freeEnergy',histFreeEnergy(1:iIteration), ... - 'feParts',histFeParts(1:iIteration), ... - 'timeSinceStart',timeSinceStart(1:iIteration)); - DcmResults.debug.idxDcmParamsInfer = idxParamsInf; - DcmResults.debug.idxSelfCon = idxSelfCon; - DcmResults.debug.nParameters = nParameters; - - if bKeepResp - DcmResults.debug.resp = struct(... - 'responseError',histRespError(1:iIteration), ... - 'responseJacobian',histRespJacobian(1:iIteration)); - end -end - -end - - - - -function [ bError ] = fnc_check_response( response, jacobian ) -% function [ bError ] = dcm_fmri_check_response( response, jacobian ) -% Checks for NaNs and Infs in DCM response and jacobian matrix. - -bError = any(isnan(response(:))) || any(isinf(response(:))) || ... - any(isnan(jacobian(:))) || any(isinf(jacobian(:))); - -end - - -function [ normWeights ] = fnc_exp_norm(logWeights) -% exponentiate and normalize log weights -% subtract max of each row for numerical stability - -normWeights = exp(bsxfun(@minus,logWeights,max(logWeights,[],2))); -normWeights = bsxfun(@rdivide,normWeights,sum(normWeights,2)); - -end - - -function [priors, posterior] = fnc_reset_cov(priors,posterior,DcmResults,DcmInfo,nDcmParamsInfCon,nDcmParamsInfAll) - - priors.clustersSigma = DcmResults.priors.clustersSigma; - posterior.clustersSigma = repmat(priors.clustersSigma,... - [1,1,DcmResults.maxClusters]); - posterior.dcmSigma = ... - priors.clustersSigma/(priors.clustersDeg-DcmInfo.nConnections-1); - posterior.dcmSigma(nDcmParamsInfCon+1:nDcmParamsInfAll,... - nDcmParamsInfCon+1:nDcmParamsInfAll) = priors.hemSigma; - posterior.dcmSigma = repmat(posterior.dcmSigma,[1,1,DcmInfo.nSubjects]); - for iSubject = 1:DcmInfo.nSubjects - posterior.logDetPostDcmSigma(iSubject) = ... - tapas_util_logdet(posterior.dcmSigma(:,:,iSubject)); - end - -end - - -function [posterior] = fnc_set_cluster_cov(priors,posterior,kmeansIdx,nAssign) - - [nClusters,dimClusters] = size(posterior.clustersMean); - nSubjects = size(posterior.dcmMean,1); - - posterior.alpha = priors.alpha + nAssign; - posterior.clustersTau = priors.clustersTau + nAssign; - posterior.clustersDeg = priors.clustersDeg + nAssign; - - for iCluster = 1:nClusters - - S = zeros(size(priors.clustersSigma)); - for iSubject = 1:nSubjects - if kmeansIdx(iSubject) == iCluster - n = posterior.dcmMean(iSubject,1:dimClusters) - ... - posterior.clustersMean(iCluster,:); - S = S + n'*n; - end - end - - m = posterior.clustersMean(iCluster,:) - priors.clustersMean; - posterior.clustersSigma(:,:,iCluster) = priors.clustersSigma + ... - S + nAssign(iCluster)*priors.clustersTau/... - (nAssign(iCluster) + priors.clustersTau)*(m'*m); - - end - -end \ No newline at end of file diff --git a/huge/tapas_huge_invert.m b/huge/tapas_huge_invert.m index 635156fc..b16f535a 100644 --- a/huge/tapas_huge_invert.m +++ b/huge/tapas_huge_invert.m @@ -1,5 +1,8 @@ -%% [DcmResults] = tapas_huge_invert(DCM, K, priors, verbose, randomize, seed) -% +function [DcmResults] = tapas_huge_invert(DCM, K, priors, verbose, randomize, seed) +% WARNING: This function is deprecated and will be removed in a future +% version of this toolbox. Please use the new object-oriented interface +% provided by the tapas_Huge class. +% % Invert hierarchical unsupervised generative embedding (HUGE) model. % % INPUT: @@ -19,14 +22,12 @@ % of REF [1]) % hemMean: prior mean of hemodynamic parameters (mu_h in Fig.1 % of REF [1]) -% hemSigma: prior covariance of hemodynamic parameters(Sigma_h +% hemSigma: prior covariance of hemodynamic parameters (Sigma_h % in Fig.1 of REF [1]) % noiseInvScale: prior inverse scale of observation noise (b_0 in % Fig.1 of REF [1]) % noiseShape: prior shape parameter of observation noise (a_0 in % Fig.1 of REF [1]) -% (you may use tapas_huge_build_prior(DCM) to generate this -% struct) % verbose - activates command line output (prints free energy % difference, default: false) % randomize - randomize starting values (default: false). WARNING: @@ -62,121 +63,101 @@ % noise (lambda_n,r in Eq.(23) of REF [1]) % modifiedSumSqrErr: b'_n,r in Eq.(22) of REF [1] % -% REFERENCE: -% [1] Yao Y, Raman SS, Schiek M, Leff A, Frssle S, Stephan KE (2018). -% Variational Bayesian Inversion for Hierarchical Unsupervised -% Generative Embedding (HUGE). NeuroImage, 179: 604-619 +% See also tapas_Huge, tapas_Huge.estimate, tapas_huge_demo % -% https://doi.org/10.1016/j.neuroimage.2018.06.073 -% % Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit +% Copyright (C) 2019 Translational Neuromodeling Unit % Institute for Biomedical Engineering, % University of Zurich and ETH Zurich. % % This file is part of TAPAS, which is released under the terms of the GNU % General Public Licence (GPL), version 3. For further details, see -% . +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. % % This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: % https://github.com/translationalneuromodeling/tapas/issues -% -function [DcmResults] = tapas_huge_invert(DCM, K, priors, verbose, randomize, seed) -%% check input -if nargin >= 6 - rng(seed); -else - seed = rng(); -end +% + + +wnMsg = ['This function is deprecated and will be removed in a future ' ... + 'version of this toolbox. Please use the new object-oriented ' ... + 'interface provided by the tapas_Huge class.']; +warning('tapas:huge:deprecated',wnMsg) -if ~isfield(DCM,'listBoldResponse') +%% check input +if isvector(DCM)&&isstruct(DCM) try - DcmInfo = tapas_huge_import_spm(DCM); - catch err - disp('tapas_huge_invert: Unsupported format.'); - disp('Use cell array of DCM in SPM format as first input.'); - rethrow(err); + DCM = {DCM(:).DCM}'; + catch + DCM = num2cell(DCM); end else - DcmInfo = DCM; + assert(iscell(DCM),'TAPAS:HUGE:inputFormat',... + 'DCM must be cell array of DCMs in SPM format'); +end + +%% settings +opts = {'K', K}; +opts = [opts, {'Dcm', DCM}]; + +if nargin >= 6 + opts = [opts, {'Seed', seed}]; end -if nargin < 5 - randomize = false; +if nargin >= 5 + opts = [opts, {'Randomize', randomize}]; end -if nargin < 4 - verbose = false; +if nargin >= 4 + opts = [opts, {'Verbose', verbose}]; end -if nargin < 3 - priors = tapas_huge_build_prior(DcmInfo); +if nargin >= 3 + dcm = DCM{1}; + nConnections = nnz([dcm.a(:);dcm.b(:);dcm.c(:);dcm.d(:)]); + priors.clustersSigma = priors.clustersSigma/... + (priors.clustersDeg - nConnections - 1); + + opts = [opts, {'PriorVarianceRatio', priors.clustersTau, ... + 'PriorDegree', priors.clustersDeg, ... + 'PriorClusterVariance', priors.clustersSigma, ... + 'PriorClusterMean', priors.clustersMean}]; end assert(K>0,'TAPAS:HUGE:clusterSize',... 'Cluster size K must to be positive integer'); -% compile integrator -tapas_huge_compile(); +%% invert model +obj = tapas_Huge(opts{:}); +obj = obj.estimate(); - -%% settings -DcmResults = struct(); +DcmResults.freeEnergy = obj.posterior.nfe; DcmResults.maxClusters = K; -DcmResults.priors = priors; - -% variational parameters -% stopping criterion: minimum increase in free energy -DcmResults.epsEnergy = 1e-5; -% stopping cirterion: maximum number of iterations -DcmResults.nIterations = 1e3; - -% computational and technical parameters -% method for calculating jacobian matrix -DcmResults.fnJacobian = @tapas_huge_jacobian; -% small constant to be added to the diagonal of inv(postDcmSigma) for -% numerical stability -DcmResults.diagConst = 1e-10; -% keep history of parameters and important auxiliary variables -DcmResults.bKeepTrace = false; -% keep history of response related auxiliary variables -% has no effect if bKeepTrace is false -DcmResults.bKeepResp = false; -DcmResults.bVerbose = verbose; - -% set update schedule -DcmResults.schedule.dfDcm = 50; -DcmResults.schedule.dfClusters = 10; -DcmResults.schedule.itAssignment = 1; -DcmResults.schedule.itCluster = 1; -DcmResults.schedule.itReturn = 5; -DcmResults.schedule.itKmeans = 1; - - -%% randomize starting values -if randomize - init = struct(); - init.dcmMean = repmat([DcmResults.priors.clustersMean,... - DcmResults.priors.hemMean], DcmInfo.nSubjects, 1); - init.dcmMean = init.dcmMean + randn(size(init.dcmMean))*.05; - - init.clustersMean = repmat(... - DcmResults.priors.clustersMean, DcmResults.maxClusters,1); - init.clustersMean = init.clustersMean + ... - randn(size(init.clustersMean))*.05; - DcmResults.init = init; -end - - -%% call VB inversion -DcmResults = tapas_huge_inv_vb(DcmInfo, DcmResults); - -DcmResults.seed = seed; -DcmResults.ver = '201903'; +DcmResults.rngSeed = obj.posterior.seed; +DcmResults.posterior = struct( ... + 'alpha', obj.posterior.alpha, ... + 'softAssign', obj.posterior.q_nk, ... + 'clustersMean' , obj.posterior.m , ... + 'clustersTau', obj.posterior.tau , ... + 'clustersDeg' , obj.posterior.nu , ... + 'clustersSigma' , obj.posterior.S , ... + 'logDetClustersSigma', [] , ... + 'dcmMean', obj.posterior.mu_n , ... + 'dcmSigma', obj.posterior.Sigma_n , ... + 'logDetPostDcmSigma', [] , ... + 'noiseShape', obj.posterior.a, ... + 'noiseInvScale', obj.posterior.b , ... + 'meanNoisePrecision', obj.posterior.a./obj.posterior.b , ... + 'modifiedSumSqrErr', []); +DcmResults.residuals = obj.trace.epsilon; end diff --git a/huge/tapas_huge_jacobian.m b/huge/tapas_huge_jacobian.m deleted file mode 100644 index 5a476225..00000000 --- a/huge/tapas_huge_jacobian.m +++ /dev/null @@ -1,73 +0,0 @@ -%% [ J ] = tapas_huge_jacobian( fnFunction, argument, params, ~, varargin ) -% -% Calculates the jacobian matrix of a function via central differences -% -% INPUT: -% fnFunction - handle to function -% argument - argument at which the jacobian is calculated -% params - struct containing parameters like step-size, etc. -% ~ - a dummy input -% -% Optional: -% varargin - additional arguments of function -% -% OUTPUT: -% J - the jacobian matrix -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ J ] = tapas_huge_jacobian( fnFunction, argument, params, ~, varargin ) - -[dimIn,~] = size(argument); -dimOut = params.dimOut; - -J = zeros(dimOut,dimIn); - -if isfield(params,'stepSize') - stepSize = params.stepSize; -else - stepSize = 1e-3; % default step size -end -if isscalar(stepSize) - stepSize = repmat(stepSize,dimIn,1); -end - -for idxIn = 1:dimIn - - currentArgument = argument; - currentArgument(idxIn) = currentArgument(idxIn) + stepSize(idxIn); - valuePlus = fnFunction(currentArgument,varargin{:}); - - currentArgument = argument; - currentArgument(idxIn) = currentArgument(idxIn) - stepSize(idxIn); - valueMinus = fnFunction(currentArgument,varargin{:}); - - J(:,idxIn) = (valuePlus - valueMinus)/2/stepSize(idxIn); - -end - -% if requested, switch to denominator layout (gradient is a column vector) -% default: numerator layout (gradient is a row vector) -if isfield(params,'denominatorLayout')&¶ms.denominatorLayout - J = J.'; -end - - - -end - - diff --git a/huge/tapas_huge_logdet.m b/huge/tapas_huge_logdet.m new file mode 100644 index 00000000..89667f4e --- /dev/null +++ b/huge/tapas_huge_logdet.m @@ -0,0 +1,39 @@ +function [ld] = tapas_huge_logdet(A) +% Numerical stable calculation of log-determinant for positive-definite +% matrix. +% +% INPUT: +% A - Positive definite matrix. +% +% OUTPUT: +% ld - log(det(A)) +% +% EXAMPLE: +% ld = TAPAS_HUGE_LOGDET(eye(3)) Calculate log-determinant of 3x3 +% identity matrix. +% + +% Author: Yu Yao (yao@biomed.ee.ethz.ch) +% Copyright (C) 2019 Translational Neuromodeling Unit +% Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of TAPAS, which is released under the terms of the GNU +% General Public Licence (GPL), version 3. For further details, see +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. +% +% This software is intended for research only. Do not use for clinical +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% + +U = chol(A); +ld = 2*sum(log(diag(U))); + +return; \ No newline at end of file diff --git a/huge/tapas_huge_logit.m b/huge/tapas_huge_logit.m new file mode 100644 index 00000000..6d7b6453 --- /dev/null +++ b/huge/tapas_huge_logit.m @@ -0,0 +1,37 @@ +function [ y ] = tapas_huge_logit( x ) +% Numerical stable calculation of logit function. +% +% INPUTS: +% x - Array of double. +% +% OUTPUTS: +% y - logit of x. +% + +% Author: Yu Yao (yao@biomed.ee.ethz.ch) +% Copyright (C) 2019 Translational Neuromodeling Unit +% Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of TAPAS, which is released under the terms of the GNU +% General Public Licence (GPL), version 3. For further details, see +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. +% +% This software is intended for research only. Do not use for clinical +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% + +y = log(x./(1-x)); +y(y == Inf) = realmax; +y(y == -Inf) = -realmin; + + +end + diff --git a/huge/tapas_huge_manual.bib b/huge/tapas_huge_manual.bib new file mode 100644 index 00000000..284c183d --- /dev/null +++ b/huge/tapas_huge_manual.bib @@ -0,0 +1,25 @@ +@article{yao2018, + author = {Yao, Yu and Raman, Sudhir S. and Schiek, Michael and Leff, Alex and Frässle, Stefan and Stephan, Klaas E.}, + title = {Variational Bayesian inversion for hierarchical unsupervised generative embedding (HUGE)}, + journal = {NeuroImage}, + volume = {179}, + pages = {604-619}, + ISSN = {1053-8119}, + doi = {10.1016/j.neuroimage.2018.06.073}, + year = {2018}, + type = {Journal Article} +} + +@article{stephan2007, + author = {Stephan, Klaas E. and Weiskopf, Nikolaus and Drysdale, Peter M. and Robinson, Peter A. and Friston, Karl J.}, + title = {Comparing hemodynamic models with DCM}, + journal = {NeuroImage}, + volume = {38}, + number = {3}, + pages = {387-401}, + ISSN = {1053-8119}, + doi = {10.1016/j.neuroimage.2007.07.040}, + year = {2007}, + type = {Journal Article} +} + diff --git a/huge/tapas_huge_manual.pdf b/huge/tapas_huge_manual.pdf new file mode 100644 index 00000000..62bfd6e7 Binary files /dev/null and b/huge/tapas_huge_manual.pdf differ diff --git a/huge/tapas_huge_manual.tex b/huge/tapas_huge_manual.tex new file mode 100644 index 00000000..3e85e7d9 --- /dev/null +++ b/huge/tapas_huge_manual.tex @@ -0,0 +1,1038 @@ + +% This LaTeX was auto-generated from MATLAB code. +% To make changes, update the MATLAB code and republish this document. + +\documentclass{article} +\usepackage{graphicx} +\usepackage[table]{xcolor} +\usepackage[utf8]{inputenc} +\usepackage{natbib} +\usepackage{lmodern} +\usepackage{enumitem} + +\sloppy +\definecolor{lgray}{gray}{0.65} +\definecolor{llgray}{gray}{0.85} +\setlength{\parindent}{0pt} +\setlength\arrayrulewidth{0.8pt} + + +\title{User Manual Hierarchical Unsupervised Generative Embedding} +\author{\\Yu Yao\\% +\\Translational Neuromodeling Unit\\% +Institute for Biomedical Engineering\\% +University of Zurich \& ETH Zurich} +\begin{document} +\maketitle + +\tableofcontents + + + +\section{Introduction} + +\begin{par} +This document contains the user guide for the HUGE toolbox, which is part of TAPAS (Translational Algorithms for Psychiatry-Advancing Science), a collection of tools developed at the Translational Neuromodeling Unit, Zurich. This guide provides a detailed documentation of the user interface of the HUGE toolbox. To get started quickly, you may run the examples in the demo script \texttt{tapas\_huge\_demo.mlx}. +\end{par} \vspace{1em} +\begin{par} +HUGE stands for Hierarchical Unsupervised Generative Embedding. It is a generative model for (task-based) fMRI data. It can be used to stratify heterogeneous cohorts into subgroups or learn the prior distribution of DCM parameters from data via empirical Bayes. For more details on the theory behind the HUGE model, please refer to \cite{yao2018}. +\end{par} \vspace{1em} +\begin{par} +To use this toolbox, change your current working directory to the directory that contains the \texttt{@tapas\_Huge} folder or add this directory to your Matlab path. +\end{par} \vspace{1em} + + +\section{The \texttt{tapas\_Huge} Class} + +\begin{par} +This toolbox implements the HUGE model as a Matlab class called \texttt{tapas\_Huge}. Data, options and results are stored in class properties, while the main functionalities of the toolbox are provided by class methods. Additional functionalities are provided by regular Matlab functions. +\end{par} \vspace{1em} +\begin{par} +Instances of the \texttt{tapas\_Huge} class are created by calling the class constructor \texttt{tapas\_Huge()}. The constructor is a function that can be called without arguments, which creates an empty \texttt{tapas\_Huge} object: +\end{par} +\begin{verbatim}obj = tapas_Huge()\end{verbatim} +\begin{par} +However, the constructor also accepts optional arguments in the form of name-value pairs. For example, one can add a short description using the \texttt{Tag} property +\end{par} +\begin{verbatim}obj = tapas_Huge('Tag','my model')\end{verbatim} +\begin{par} +or import data using the \texttt{Dcm} property +\end{par} +\begin{verbatim}obj = tapas_Huge('Dcm',dcms)\end{verbatim} +\begin{par} +For more examples, see the demo script \texttt{tapas\_huge\_demo.mlx}. +\end{par} \vspace{1em} + + +\section{Class Properties} + +\begin{par} +Each instance of the \texttt{tapas\_Huge} class stores data, options and results in class properties, which are documented below. Properties can be accessed using dot notation: +\end{par} +\begin{verbatim}value = obj.property\end{verbatim} +\begin{par} +where \texttt{obj} is a class instance and \texttt{property} is the name of the property. +\end{par} \vspace{1em} +\begin{par} +Note that all properties, except for \texttt{K} (the number of clusters) and \texttt{tag} are read only, i.e. they cannot be changed by the user directly. Instead, their value is set automatically when importing data using the \texttt{import} method or inverting the model using the \texttt{estimate} method. +\end{par} \vspace{1em} +\begin{par} +Below is a list of all properties of the \texttt{tapas\_Huge} class. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{K}} + +\begin{par} +A scalar integer containing the number of clusters of the HUGE model. +\end{par} \vspace{1em} + + +\subsection*{ tag} + +\begin{par} +A character array that can be used for storing a short description of the model. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{L}} + +\begin{par} +A scalar integer containing the number of experimental inputs (i.e., the dimensionality of \texttt{obj.inputs.u}). +\end{par} \vspace{1em} + + +\subsection*{ \texttt{M}} + +\begin{par} +A scalar integer containing the number of group-level confound variables (like age, sex, handedness, etc). +\end{par} \vspace{1em} + + +\subsection*{ \texttt{N}} + +\begin{par} +A scalar integer containing the number of subjects. Permissions: read only. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{R}} + +\begin{par} +A scalar integer containing the number of brain regions in the model. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{idx}} + +\begin{par} +A Matlab struct storing parameter indices with the following fields: +\end{par} \vspace{1em} +\begin{itemize}[align=parleft, labelsep=10ex, leftmargin=15ex] +\setlength{\itemsep}{-1ex} + \item [ \texttt{clustering} ] Vector of indices of DCM parameters used for clustering. + \item [ \texttt{homogenous} ] Vector of indices of homogenous parameters. I.e. DCM parameters which are estimated but not used for clustering. + \item [ \texttt{P\_c} ] Number of clustering parameters. + \item [ \texttt{P\_c} ] Number of homogeneous parameters. + \item [ \texttt{P\_f} ] Number of parameters in a fully connected DCM with the same number of regions as the current one. +\end{itemize} + + +\subsection*{ \texttt{dcm}} + +\begin{par} +A Matlab struct in SPM's DCM format containing the DCM network structure. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{inputs}} + +\begin{par} +A struct array containing the experimental stimuli for all subject, with two fields: +\end{par} \vspace{1em} +\begin{itemize}[align=parleft, labelsep=10ex, leftmargin=15ex] +\setlength{\itemsep}{-1ex} + \item [ \texttt{u} ] Array containing experimental stimuli. + \item [ \texttt{dt} ] Sampling time interval. +\end{itemize} + + +\subsection*{ \texttt{data}} + +\begin{par} +A struct array containing the fMRI data for all subjects, with fields: +\end{par} \vspace{1em} +\begin{itemize}[align=parleft, labelsep=10ex, leftmargin=15ex] +\setlength{\itemsep}{-1ex} + \item [ \texttt{bold} ] Array containing the BOLD time series. + \item [ \texttt{te} ] The echo time (TE). + \item [ \texttt{tr} ] The repetition time TR. + \item [ \texttt{X0} ] Vector containing subject-level confounds. + \item [ \texttt{res} ] The residual forming matrix of X0. + \item [ \texttt{confounds} ] Vector containing the group-level confounds (e.g., sex, age, etc). + \item [ \texttt{spm} ] A struct containing the posterior from SPM (The Ep field in SPM's DCM format). +\end{itemize} + + +\subsection*{ \texttt{labels}} + +\begin{par} +A struct containing names for regions and inputs. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{options}} + +\begin{par} +A struct storing options for the HUGE model. Note: Set options using name-value pair arguments with the \texttt{taps\_huge} and \texttt{estimate} methods. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{prior}} + +\begin{par} +A struct storing priors of the HUGE model. Note: Set priors using name-value pair arguments with the \texttt{taps\_huge} and \texttt{estimate} methods. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{posterior}} + +\begin{par} +A struct storing the estimation results. The fields correspond to the parameters of the posterior distribution (see\cite{yao2018}). +\end{par} \vspace{1em} +\begin{itemize}[align=parleft, labelsep=10ex, leftmargin=15ex] +\setlength{\itemsep}{-1ex} + \item [ \texttt{alpha} ] Vector of posterior cluster weights ($\alpha$). + \item [ \texttt{m} ] Array of posterior cluster means ($m$). + \item [ \texttt{tau} ] Vector containing the $\tau$ parameter of the inverse-Wishart posterior distribution. + \item [ \texttt{nu} ] Vector of posterior degree of freedom parameters ($\nu$). + \item [ \texttt{S} ] Array containing the scale matrices of the inverse-Wishart posterior distribution ($S$). + \item [ \texttt{q\_nk} ] Array of posterior assignment probabilities ($q_{nk}$). + \item [ \texttt{mu\_n} ] Array of posterior mean subject-level DCM parameters ($\mu_n$). + \item [ \texttt{Sigma\_n} ] Array of posterior covariances of subject-level DCM parameters ($\Sigma_n$). + \item [ \texttt{b} ] Vector containing the $b$ parameter of the posterior Gamma distribution over noise precision. + \item [ \texttt{a} ] Vector containing the $a$ parameter of the posterior Gamma distribution over noise precision. + \item [ \texttt{m\_beta} ] When using group-level confounds, array containing posterior mean of coefficients, ($m_{\beta}$). Otherwise, empty array. + \item [ \texttt{S\_beta} ] When using group-level confounds, array containing posterior covariance of coefficients ($S_{\beta}$). Otherwise, empty array. + \item [ \texttt{nfe} ] The negative free energy ($F$). + \item [ \texttt{mu\_r} ] When using group-level confounds, array containing posterior mean of residual DCM parameters, ($m_r$). Otherwise, empty array. + \item [ \texttt{Sigma\_r} ] When using group-level confounds, array containing posterior covariance of residual DCM parameters, ($S_r$). Otherwise, empty array. + \item [ \texttt{nrv} ] Normalized residual variance. An array containing the amount of non-explained variance (1 minus variance explained) per subject and per region. + \item [ \texttt{bPurity} ] When using synthetic data, balanced purity of estimation result. + \item [ \texttt{method} ] Character array indicating inversion method. + \item [ \texttt{version} ] Character array indicating toolbox version. + \item [ \texttt{seed} ] Random number generator seed. +\end{itemize} + + +\subsection*{ \texttt{trace}} + +\begin{par} +A struct containing convergence diagnostics, with fields: +\end{par} \vspace{1em} +\begin{itemize}[align=parleft, labelsep=10ex, leftmargin=15ex] +\setlength{\itemsep}{-1ex} + \item [ \texttt{nfe} ] Vector containing the history of the negative free energy during iterative VB inversion. + \item [ \texttt{nDcmUpdate} ] Vector containing the number of times the subject-level DCM parameter update was accepted. + \item [ \texttt{convergence} ] The number of iterations it took VB to converge. + \item [ \texttt{kmeans} ] Struct containing information about the K-means initialization step. + \item [ \texttt{epsilon} ] Cell array containing the residual (observed minus predicted BOLD time series) for each subject at convergence. +\end{itemize} + + +\subsection*{ \texttt{aux}} + +\begin{par} +A struct used for storing auxiliary variables during model estimation. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{model}} + +\begin{par} +When using the model to simulate data, this struct contains the model parameters used for generating the synthetic dataset. Fields correspond to HUGE model parameters: +\end{par} \vspace{1em} +\begin{itemize}[align=parleft, labelsep=10ex, leftmargin=15ex] +\setlength{\itemsep}{-1ex} + \item [ \texttt{pi} ] Vector of cluster weights ($\pi$). + \item [ \texttt{d} ] Array of cluster labels in one-hot encoding ($d$). + \item [ \texttt{mu\_k} ] Array of cluster means ($\mu_k$). + \item [ \texttt{Sigma\_k} ] Array of cluster covariances ($\Sigma_k$). + \item [ \texttt{mu\_h} ] Vector containing mean of homogenous DCM parameters ($\mu_h$). + \item [ \texttt{Sigma\_h} ] Array containing covariance of homogenous DCM parameters ($\Sigma_h$). + \item [ \texttt{theta\_c} ] Array of clustering DCM parameters ($\theta_c$). + \item [ \texttt{theta\_h} ] Array of homogenous DCM parameters ($\theta_h$). + \item [ \texttt{lambda} ] Array of measurement noise precisions ($\lambda$). + \item [ \texttt{x\_n} ] Array of group-level confounds (e.g. age, sex, etc. $x_n$). + \item [ \texttt{beta} ] Array of confound coefficients ($\beta$). + \item [ \texttt{options} ] Struct containing options for simulation. + \item [ \texttt{seed} ] Random number generator seed. +\end{itemize} + + +\subsection*{ \texttt{constants}} + +\begin{par} +Struct storing constants used by the model. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{version}} + +\begin{par} +Character array indicating toolbox version. +\end{par} \vspace{1em} + + +\section{Class Methods} + +\begin{par} +The main functionalities of the HUGE toolbox are implemented as methods of the \texttt{tapas\_Huge} class. These methods are documented below. There are two equivalent ways to call class methods: +\end{par} +\begin{verbatim}obj = obj.method( ... )\end{verbatim} +\begin{par} +or +\end{par} +\begin{verbatim}obj = method( obj, ... )\end{verbatim} +\begin{par} +where \texttt{obj} is an instance of the \texttt{tapas\_Huge} class and \texttt{method} is the class method you want to call. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{tapas\_Huge.estimate}} + + \begin{verbatim} Estimate parameters of the HUGE model. + + INPUTS: + obj - A tapas_Huge object containing fMRI time series. + + OPTIONAL INPUTS: + This function accepts optional name-value pair arguments. For a list of + valid name-value pairs, see the user manual or type 'help + tapas_huge_property_names'. + + OUTPUTS: + obj - A tapas_Huge object containing the estimation result in the + 'posterior' property. + + EXAMPLES: + [obj] = ESTIMATE(obj) Invert the HUGE model stored in obj. + + [obj] = ESTIMATE(obj, 'K', 2) Set the number of clusters to 2 and + invert the HUGE model stored in obj. + + [obj] = ESTIMATE(obj, 'Verbose', true) Print progress of estimation + to command line. + + [obj] = ESTIMATE(obj, 'Dcm', dcms, 'OmitFromClustering', 1) Import + data stored in 'dcms' and omit self-connections from clustering. + + See also TAPAS_HUGE_DEMO + + +\end{verbatim} \color{black} + \begin{par} +Note that the HUGE toolbox uses a parameterization such that the DCM networks are self-inhibiting by default. This is achieved by subtracting 0.5 from the self-connections. Hence, specifying a prior mean of zero for all DCM connections implies that the effective self-connectivity is -0.5. This parametrization has been chosen for the convenience of the user, since it eliminates the need to identify the position of the self-connections in the parameter vector. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{tapas\_Huge.simulate}} + + \begin{verbatim} Generate synthetic task-based fMRI time series data, using HUGE as a + generative model. + + INPUTS: + obj - A tapas_Huge object. + clusters - A cell array containing DCM structs in SPM's DCM format, + indicating the DCM network structure and cluster mean + parameters. + sizes - A vector containing the number of subjects for each cluster. + + OPTIONAL INPUTS: + This function accepts optional name-value pair arguments. For a list of + valid name-value pairs, see examples below. + + OUTPUTS: + obj - A tapas_Huge object containing the simulated fMRI time series in + its 'data' property and the ground truth parameters in its + 'model' property. + + EXAMPLES: + [obj] = SIMULATE(obj,clusters,sizes) Simulate fMRI time series with + cluster mean parameters given in 'clusters' and number of subjects + given in 'sizes'. + + [obj] = SIMULATE(obj,clusters,sizes,'Snr',1) Set signal-to-noise + ratio of fMRI data to 1. + + [obj] = SIMULATE(obj,clusters,sizes,'NoiseFloor',0.1) Set minimum + noise variance to 0.1. + + [obj] = SIMULATE(obj,clusters,sizes,'confounds',confounds) Introduce + group-level confounds (like sex or age). + + [obj] = SIMULATE(obj,clusters,sizes,'beta',beta) Set coefficients + for group-level confounds. + + [obj] = SIMULATE(obj,clusters,sizes,'variant',2) Set confound + variant to 2 (i.e., clusters do not share confound coefficients). + + [obj] = SIMULATE(obj,clusters,sizes,'Inputs',U) Introduce subject- + specific experimental stimuli. 'U' must be an array or cell array + of structs with fields 'dt and 'u'. + + [obj] = SIMULATE(obj,clusters,sizes,'OmitFromClustering',omit) + Designate DCM parameters to be excluded from clustering model. + Excluded parameters still exist in the DCM network structure, but + are drawn from the same distribution for all subjects. 'omit' + should be a struct with fields a, b, c, and/or d which are bool + arrays with sizes matching the corresponding fields of the DCMs. If + omit is an array, it is interpreted as the field a. If omit is 1, + it is expanded to an identity matrix of suitable size. + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_Huge.plot}} + + \begin{verbatim} Plot cluster and subject-level estimation result from HUGE model. + + INPUTS: + obj - A tapas_Huge object containing estimation results. + + OPTIONAL INPUTS: + subjects - A vector containing indices of subjects for which to plot + detailed results. + + OUTPUTS: + fHdl - Handle of first figure. + + See also tapas_Huge.ESTIMATE + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_Huge.save}} + + \begin{verbatim} Save properties of HUGE object to mat file. + + INPUTS: + filename - File name. + obj - A tapas_Huge object. + OPTIONAL INPUTS: + Names of property to be saved. See examples below. + + EXAMPLES: + SAVE(filename,obj) Save properties of 'obj' as individual variables + file specified in 'filename'. + + SAVE(filename,obj,'options','posterior','prior') Only save the + 'options', 'posterior' and 'prior' properties of 'obj'. + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_Huge.import}} + + \begin{verbatim} Import fMRI time series data into HUGE object. + + WARNING: Importing data into a HUGE object will delete any data and + results which are already stored in that object. + + INPUTS: + obj - A tapas_Huge object. + dcms - A cell array containing DCM structs in SPM's DCM format. + + OPTIONAL INPUTS: + confounds - Group-level confounds (e.g., age, sex, etc). 'confounds' + must be empty or an array with as many rows as there are + elements in 'dcm'. + omit - specifies DCM parameters which should be omitted from + clustering. Parameters omitted from clustering will still + be estimated, but under a static Gaussian prior. 'omit' + should be a struct with fields a, b, c, and/or d which are + bool arrays with sizes matching the corresponding fields of + the DCMs. If omit is an array, it is interpreted as the + field a. If omit is 1, it is expanded to an identity matrix + of suitable size. + append - bool, if true keep current fMRI time series and append new + data in 'dcms'. Note: estimation results will still be + deleted. + + OUTPUTS: + obj - A tapas_Huge object containing the imported data. + + EXAMPLES: + [obj] = IMPORT(obj,dcms) Import the fMRI time series and DCM + network structure stored in dcms into obj. + + [obj] = IMPORT(obj,dcms,confounds) Import group-level confounds + (like age or sex) in addition to fMRI data. + + [obj] = IMPORT(obj,dcms,[],1) Import fMRI data and network + structure. Exclude self-connections from clustering. + + See also tapas_Huge.REMOVE, tapas_Huge.EXPORT + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_Huge.export}} + + \begin{verbatim} Export results and data from HUGE object to SPM's DCM format. + + INPUTS: + obj - A tapas_Huge object containing data. + + OUTPUTS: + dcms - A cell array containing DCM structs in SPM's DCM format. + confounds - An array containing group-level confounds (like age or + sex) if available. 'confounds' is an array with one row + per subject. + + EXAMPLES: + [dcms] = EXPORT(obj) Export fMRI time series and estimation results + stored in obj. + + [dcms,confounds] = EXPORT(obj) Also export group-level confounds + (like age or sex). + + See also tapas_Huge.IMPORT + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_Huge.remove}} + + \begin{verbatim} Remove data (fMRI time series, confounds, DCM network structure, ... ) + and estimation results from HUGE object. + + INPUTS: + obj - A tapas_Huge object. + + OPTIONAL INPUTS: + idx - Only remove data of the subjects indicated in 'idx'. 'idx' must + be a vector containing numeric or logical array indices. + + OUTPUTS: + obj - A tapas_Huge object with data and results removed. + + EXAMPLES: + [obj] = REMOVE(obj) Remove results and data of all subjects. + + [obj] = REMOVE(obj,1:5) Remove results and data for first 5 + subjects. + + [obj] = REMOVE(obj,'all') is the same as [obj] = REMOVE(obj) + + [obj] = REMOVE(obj,0) is the same as [obj] = REMOVE(obj) + + See also tapas_Huge.IMPORT + + +\end{verbatim} \color{black} + + +\section{Name-Value Pair Arguments} + +\begin{par} +Both the class constructor \texttt{tapas\_Huge} and the method \texttt{estimate} accept optional arguments in the form of name-value pairs. Options set using name-value pair arguments are persistent across the lifetime of the object. Below is a list of valid name-value pairs that can be used with these two functions. +\end{par} \vspace{1em} +\begin{par} + +{\renewcommand{\arraystretch}{1.5} +\renewcommand{\tabcolsep}{0.2cm} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+Confounds+ \\ +\hline +\rowcolor{llgray} Value: & double array\\ +\hline + \rowcolor{llgray} Description: & Specify confounds for group-level + analysis (e.g. age or sex) as double array with one row per subject and + one column per confound. Note: This property can only be used in + combination with the \verb+Dcm+ property.\\ +\rowcolor{llgray} & WARNING: This feature is experimental and has not +been extensively tested. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+ConfoundsVariant+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'none'+ \ensuremath{|} \verb+'global'+ \ensuremath{|} \verb+'cluster'+ (default: \verb+'global'+ if confounds specified, \verb+'none'+ otherwise)\\ +\hline + \rowcolor{llgray} Description: & Determines how confounds enter model. \verb+'none'+: Confounds are not used. \verb+'global'+: Confounds enter global regression (variant 1). \verb+'cluster'+: Confounds enter cluster-specific regression (variant 2). +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+Dcm+ \\ +\hline +\rowcolor{llgray} Value: & cell array of DCM structs in SPM format\\ +\hline + \rowcolor{llgray} Description: & Specify DCM structure and BOLD time series for all subjects as cell array with one DCM struct in SPM format per subject. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+K+ \\ +\hline +\rowcolor{llgray} Value: & positive integer (default: 1)\\ +\hline + \rowcolor{llgray} Description: & Number of clusters. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+Method+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'VB'+\\ +\hline + \rowcolor{llgray} Description: & Name of inversion method specified as character array. VB: variational Bayes +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+NumberOfIterations+ \\ +\hline +\rowcolor{llgray} Value: & positive integer (default: 999)\\ +\hline + \rowcolor{llgray} Description: & Maximum number of iterations of VB scheme. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+OmitFromClustering+ \\ +\hline +\rowcolor{llgray} Value: & array of logical \ensuremath{|} struct with fields \verb+a+, \verb+b+, \verb+c+ and \verb+d+ (default: empty struct)\\ +\hline + \rowcolor{llgray} Description: & Select DCM parameters to exclude from clustering. Parameters excluded from clustering will still be estimated, but under a static Gaussian prior. If input is array, it will be treated as the \verb+a+ field of a struct. Missing fields will be treated the same as arrays of \verb+false+. Note: This property can only be used in combination with the \verb+Dcm+ property. +\\ +\hline + \rowcolor{llgray} Example: & Exclude the self-connections from clustering: \\ + \rowcolor{llgray} & \verb+obj = obj.estimate('OmitFromClustering', 1);+ \\ + \rowcolor{llgray} & Exclude input strength from clustering: \\ + \rowcolor{llgray} & \verb+omit = struct('c',obj.dcm.c));+ \\ + \rowcolor{llgray} & \verb+obj = obj.estimate('OmitFromClustering', omit);+ \\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+PriorClusterMean+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'default'+ \ensuremath{|} row vector of double\\ +\hline + \rowcolor{llgray} Description: & Prior cluster mean. Scalar input will be expanded into vector. Default: \verb+[0, ... ,0]\verb+. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+PriorClusterVariance+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'default'+ \ensuremath{|} symmetric, positive definite matrix of double\\ +\hline + \rowcolor{llgray} Description: & Prior mean of cluster covariances. Must be a symmetric, positive definite matrix. Scalar input will be expanded into diagonal matrix. Default: \verb+diag([0.01, ..., 0.01])+. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+PriorDegree+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'default'+ \ensuremath{|} positive double\\ +\hline + \rowcolor{llgray} Description: & $\nu_0$ determines the prior precision of the cluster covariance. For VB, this is the degrees of freedom of the inverse-Wishart. Default: \verb+100+. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+PriorVarianceRatio+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'default'+ \ensuremath{|} positive double\\ +\hline + \rowcolor{llgray} Description: & Ratio $\tau_0$ between prior mean cluster covariance and prior covariance over cluster mean. Prior covariance over cluster mean equals prior cluster covariance divided $\tau_0$. Default: \verb+0.1+. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+Randomize+ \\ +\hline +\rowcolor{llgray} Value: & bool (default: \verb+false+)\\ +\hline + \rowcolor{llgray} Description: & If \verb+true+, starting values for subject level DCM parameter estimates are randomized. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+SaveTo+ \\ +\hline +\rowcolor{llgray} Value: & character array\\ +\hline + \rowcolor{llgray} Description: & Location for saving results consisting of path name and optional file name. Path name must end on file separator and point to an existing directory. If file name is not specified, it is set to 'huge' followed by date and time. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+Seed+ \\ +\hline +\rowcolor{llgray} Value: & double \ensuremath{|} cell array of double and rng name \ensuremath{|} random number generator seed obtained with \verb+rng+ command\\ +\hline + \rowcolor{llgray} Description: & Seed for random number generator. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+StartingValueDcm+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'default'+ \ensuremath{|} \verb+'spm'+ \ensuremath{|} double array (default: \verb+'default'+)\\ +\hline + \rowcolor{llgray} Description: & Starting values for subject-level DCM parameter estimates. \verb+'default'+ uses prior cluster mean for all subjects. 'spm' uses values supplied in the 'Ep' field of the SPM DCM structs. Use a double array with number of rows equal to number of subjects to specify custom starting values. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+StartingValueGmm+ \\ +\hline +\rowcolor{llgray} Value: & \verb+'default'+ \ensuremath{|} double array (default: \verb+'default'+)\\ +\hline + \rowcolor{llgray} Description: & Starting values for cluster-level DCM parameter estimates. \verb+'default'+ uses prior cluster mean for all clusters. Use a double array with number of rows equal to number of clusters to specify custom starting values. +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+Tag+ \\ +\hline +\rowcolor{llgray} Value: & character array\\ +\hline + \rowcolor{llgray} Description: & Model description +\\ +\hline +\end{tabular} +\vspace{1em} +\begin{tabular}{|l|p{10cm}|} +\hline +\rowcolor{lgray} Name: & \verb+Verbose+ \\ +\hline +\rowcolor{llgray} Value: & bool (default: \verb+false+)\\ +\hline + \rowcolor{llgray} Description: & Activate/deactivate command line output. +\\ +\hline +\end{tabular} +\vspace{1em} +} + +\end{par} \vspace{1em} + + +\section{Other Functions} + +\begin{par} +The HUGE toolbox provides additional functionalities in the form of regular Matlab functions, which are documented below. These function are used internally by the \texttt{tapas\_Huge} class. However, some of them may be useful for the user. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{tapas\_huge\_boxcar}} + + \begin{verbatim} Generate a boxcar function for use as experimental stimulus. All + timing-related arguments must be specified in seconds. + + INPUTS: + dt - Numeric scalar indicating sampling time interval. + nBoxes - Vector indicating number of blocks. + period - Vector containing time interval between block onsets. + onRatio - Vector containing ratio between block length and 'period'. + Must be between 0 and 1. + + OPTIONAL INPUTS: + padding - Length of padding at the beginning and end. + + OUTPUTS: + u - A cell array containing the boxcar functions. + + EXAMPLES: + u = TAPAS_HUGE_BOXCAR(.01, 10, 3, 2/3, [4 0]) Generate boxcar + function with 10 blocks, each 2 seconds long with 1 second inter + block interval and onset of first block at 4 seconds. + + See also tapas_Huge.SIMULATE + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_huge\_bpurity}} + + \begin{verbatim} Calculate balanced purity (see Brodersen2014 Eq. 13 and 14) for a set of + ground truth labels and a set of estimated labels + + INPUTS: + labels - Vector of ground truth class labels. + estimates - Clustering result as array of assignment probabilities or + vector of cluster indices. + + OUTPUTS: + balancedPurity - Balanced purity score according to Brodersen (2014) + + EXAMPLES: + bp = TAPAS_HUGE_BPURITY(labels,estimates) + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_huge\_bold}} + + \begin{verbatim} Integrates the DCM forward equations to generate the predicted fMRI bold + time series. + + INPUTS: + A, B, C, D - DCM connectivity matrices. + tau - Venous transit time. + kappa - Decay of vasodilatory signal. + epsilon - Ratio of intra- and extravascular signal. + R - Number of regions. + u - Experimental stimuli. + L - Number of experimental stimuli. + E_0 - Resting oxygen extraction fraction. + r_0 - Slope of intravascular relaxation rate. + V_0 - Resting venous volume. + vartheta_0 - Frequency offset at the outer surface of magnetized + vessels (Hz). + alpha - Grubb's exponent. + gamma - rate constant of feedback regulation. + TR - Repetition time. + TE - Echo time. + dt - Sampling interval of inputs. + + OUTPUTS: + response - matrix of predicted response for each region + (column-wise) + x - time series of neuronal states + s - time series of vasodilatory signal + f1 - time series of flow + v1 - time series of blood volume + q1 - time series of deoxyhemoglobin content. + + +\end{verbatim} \color{black} + \begin{par} +For more information on the fMRI BOLD model, please refer to \cite{stephan2007}. +\end{par} \vspace{1em} + + +\subsection*{ \texttt{tapas\_huge\_compile}} + +\begin{par} +This function is used internally to compile the mex function the HUGE toolbox needs. This happens automatically the first time you use the HUGE toolbox. In general, there is no need to call this function manually. However, you need to make sure that a C compiler is installed on your system. More details on this topic, including links to freely available compilers, can be found at: \begin{verbatim}https://www.mathworks.com/support/requirements/supported-compilers.html\end{verbatim} +\end{par} \vspace{1em} + + +\subsection*{ \texttt{tapas\_huge\_logdet}} + +\begin{par} +This function is intended for internal use only. Do not call directly +\end{par} + \begin{verbatim} Numerical stable calculation of log-determinant for positive-definite + matrix. + + INPUT: + A - Positive definite matrix. + + OUTPUT: + ld - log(det(A)) + + EXAMPLE: + ld = TAPAS_HUGE_LOGDET(eye(3)) Calculate log-determinant of 3x3 + identity matrix. + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_huge\_parse\_inputs}} + +\begin{par} +This function is intended for internal use only. Do not call directly +\end{par} + \begin{verbatim} Parse name-value pair type arguments into a struct. + + INPUTS: + opts - Struct containing all valid names as field names and + corresponding default values as field values. + + OPTIONAL INPUTS: + Name-value pair arguments. + + OUTPUTS: + opts - Struct containing name-value pair input arguments as fields. + + EXAMPLE: + opts = TAPAS_HUGE_PARSE_INPUTS(struct('a',0,'b',1),'a',10) Specify + 'a' and 'b' as valid property names and receive non-default value + for 'a'. + + +\end{verbatim} \color{black} + \begin{par} +The following two functions provide backward compatibility to the interface of the previous version of the HUGE toolbox (version 201903). +\end{par} \vspace{1em} + + +\subsection*{ \texttt{tapas\_huge\_invert}} + + \begin{verbatim} WARNING: This function is deprecated and will be removed in a future + version of this toolbox. Please use the new object-oriented interface + provided by the tapas_Huge class. + + Invert hierarchical unsupervised generative embedding (HUGE) model. + + INPUT: + DCM - cell array of DCM in SPM format + K - number of clusters (set K to one for empirical Bayes) + + OPTIONAL INPUT: + priors - model priors stored in a struct containing the + following fields: + alpha: parameter of Dirichlet prior (alpha_0 in Fig.1 of + REF [1]) + clustersMean: prior mean of clusters (m_0 in Fig.1 of REF [1]) + clustersTau: tau_0 in Fig.1 of REF [1] + clustersDeg: degrees of freedom of inverse-Wishart prior (nu_0 in + Fig.1 of REF [1]) + clustersSigma: scale matrix of inverse-Wishart prior (S_0 in Fig.1 + of REF [1]) + hemMean: prior mean of hemodynamic parameters (mu_h in Fig.1 + of REF [1]) + hemSigma: prior covariance of hemodynamic parameters (Sigma_h + in Fig.1 of REF [1]) + noiseInvScale: prior inverse scale of observation noise (b_0 in + Fig.1 of REF [1]) + noiseShape: prior shape parameter of observation noise (a_0 in + Fig.1 of REF [1]) + verbose - activates command line output (prints free energy + difference, default: false) + randomize - randomize starting values (default: false). WARNING: + randomizing starting values can cause divergence of DCM. + seed - seed for random number generator + + OUTPUT: + DcmResults - struct used for storing the results from VB. Posterior + parameters are stored in DcmResults.posterior, which is a + struct containing the following fields: + alpha: parameter of posterior over cluster weights + (alpha_k in Eq.(15) of REF [1]) + softAssign: posterior assignment probability of subjects + to clusters (q_nk in Eq.(18) in REF [1]) + clustersMean: posterior mean of clusters (m_k in Eq.(16) of + REF [1]) + clustersTau: tau_k in Eq.(16) of REF [1] + clustersDeg: posterior degrees of freedom (nu_k in Eq.(16) + of REF [1]) + clustersSigma: posterior scale matrix (S_k in Eq.(16) of + REF [1]) + logDetClustersSigma: log-determinant of S_k + dcmMean: posterior mean of DCM parameters (mu_n in + Eq.(19) of REF [1]) + dcmSigma: posterior covariance of hemodynamic + parameters (Sigma_n in Eq.(19) of REF [1]) + logDetPostDcmSigma: log-determinant of Sigma_n + noiseInvScale: posterior inverse scale of observation noise + (b_n,r in Eq.(21) of REF [1]) + noiseShape: posterior shape parameter of observation noise + (a_n,r in Eq.(21) of REF [1]) + meanNoisePrecision: posterior mean of precision of observation + noise (lambda_n,r in Eq.(23) of REF [1]) + modifiedSumSqrErr: b'_n,r in Eq.(22) of REF [1] + + See also tapas_Huge, tapas_Huge.estimate, tapas_huge_demo + + +\end{verbatim} \color{black} + + +\subsection*{ \texttt{tapas\_huge\_build\_prior}} + + \begin{verbatim} WARNING: This function is deprecated and will be removed in a future + version of this toolbox. Please use the new object-oriented interface + provided by the tapas_Huge class. + + Generate values for prior parameters for HUGE. Prior mean of cluster + centers and prior mean and covariance of hemodynamic parameters follow + SPM convention (SPM8 r6313). + + INPUT: + DcmInfo - cell array of DCM in SPM format + + OUTPUT: + priors - struct containing priors + DcmInfo - struct containing DCM model specification and BOLD time + series in DcmInfo format + + See also tapas_Huge, tapas_Huge.estimate, tapas_huge_demo + + +\end{verbatim} \color{black} + + +\section{Licence and Support} + +\begin{par} +This toolbox is part of TAPAS, which is released under the terms of the GNU General Public Licence (GPL), version 3. For further details, see: \begin{verbatim}https://www.gnu.org/licenses/\end{verbatim} +\end{par} \vspace{1em} +\begin{par} +The software contained in this toolbox is provided "as is", without warranty of any kind, express or implied, including, but not limited to the warranties of merchantability, fitness for a particular purpose and non-infringement. +\end{par} \vspace{1em} +\begin{par} +This software is intended for research only. Do not use for clinical purpose. Please note that this toolbox is under active development. Considerable changes may occur in future releases. +\end{par} \vspace{1em} +\begin{par} +Support for this toolbox is provided via the GitHub page of TAPAS. For questions and bug reports, please visit: \begin{verbatim}https://github.com/translationalneuromodeling/tapas/issues\end{verbatim} +\end{par} \vspace{1em} +\begin{par} + +\bibliographystyle{plainnat} +\bibliography{tapas_huge_manual} + +\end{par} \vspace{1em} + + + +\end{document} + diff --git a/huge/tapas_huge_nfe.m b/huge/tapas_huge_nfe.m deleted file mode 100644 index 2405f98c..00000000 --- a/huge/tapas_huge_nfe.m +++ /dev/null @@ -1,428 +0,0 @@ -%% [freeEnergyFull, feAux] = tapas_huge_nfe( counts, priors, posterior, feAux, listSubjects ) -% -% Calculates or updates negative free energy (NFE) for HUGE. -% -% INPUT: -% counts - number of parameters and subjects -% priors - parameters of prior distribution -% posterior - parameters of approximate posterior distribution. -% feAux - struct containing intermediate values of the NFE. -% listSubjects - indices of subjects for which to update the NFE. -% -% OUTPUT: -% freeEnergyFull - value of the negative free energy. -% feAux - struct containing intermediate values of the NFE. -% -% REFERENCE: -% -% Yao Y, Raman SS, Schiek M, Leff A, Frssle S, Stephan KE (2018). -% Variational Bayesian Inversion for Hierarchical Unsupervised Generative -% Embedding (HUGE). NeuroImage, 179: 604-619 -% -% https://doi.org/10.1016/j.neuroimage.2018.06.073 -% - -% -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [freeEnergyFull, feAux] = tapas_huge_nfe(counts,priors,posterior,feAux,listSubjects) -% -fepIdx = 0; -freeEnergyFullPart = zeros(100,1); - - -nParameters = counts.nParameters; -nDcmParamsInfCon = nParameters(2,1); %%% d_c -nDcmParamsInfHem = nParameters(2,2); %%% d_h -nDcmParamsInfAll = nParameters(2,3); %%% d - -if nargin < 4 - feAux.line2 = zeros(counts.nSubjects,counts.nClusters,2); - feAux.line4 = zeros(counts.nSubjects,2); -end -if nargin < 5 - listSubjects = 1:counts.nSubjects; -end - - - -%% line 1 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -sum(sum(bsxfun(@times,posterior.softAssign,... - posterior.logDetClustersSigma.')/2)); - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -sum(posterior.softAssign(:))*nDcmParamsInfCon*log(pi)/2; - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(sum(... - bsxfun(@times,... - posterior.softAssign,... - sum(psi(bsxfun(@plus,... - posterior.clustersDeg,... - 1-(1:nDcmParamsInfCon))/2),2).'... - )... - ))/2; - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -sum(sum(bsxfun(@times,... - posterior.softAssign,... - nDcmParamsInfCon./posterior.clustersTau.')... - ))/2; - - - -%% line 2 -fepIdx = fepIdx + 1; -for iSub = listSubjects - for iCluster = 1:counts.nClusters - feAux.line2(iSub,iCluster,1) = ... - -posterior.clustersDeg(iCluster)*... - posterior.softAssign(iSub,iCluster)*(... - trace(... - posterior.clustersSigma(:,:,iCluster)\... - posterior.dcmSigma(1:nDcmParamsInfCon,... - 1:nDcmParamsInfCon,... - iSub... - )... - )... - )/2; - end -end -freeEnergyFullPart(fepIdx) = sum(sum(feAux.line2(:,:,1))); -%%% - -fepIdx = fepIdx + 1; -for iSub = listSubjects - for iCluster = 1:counts.nClusters - feAux.line2(iSub,iCluster,2) = ... - -posterior.clustersDeg(iCluster)*... - posterior.softAssign(iSub,iCluster)*(... - (posterior.dcmMean(iSub,1:nDcmParamsInfCon)... - -posterior.clustersMean(iCluster,:))*... - (posterior.clustersSigma(:,:,iCluster)\... - (posterior.dcmMean(iSub,1:nDcmParamsInfCon)... - -posterior.clustersMean(iCluster,:)).' ... - )... - )/2; - end -end -freeEnergyFullPart(fepIdx) = sum(sum(feAux.line2(:,:,2))); -%%% - - -%% line 3 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -sum(counts.nMeasurements)*log(2*pi)/2; - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(sum(... - bsxfun(@times,... - counts.nMeasurementsPerState,... - psi(posterior.noiseShape) - log(posterior.noiseInvScale))... - ))/2; - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -posterior.meanNoisePrecision(:).'*posterior.modifiedSumSqrErr(:)/2; - - -%% line 4 -fepIdx = fepIdx + 1; -for iSub = listSubjects - feAux.line4(iSub,1) = ... - -trace(priors.hemSigma\posterior.dcmSigma(nDcmParamsInfCon+1:end,... - nDcmParamsInfCon+1:end,... - iSub... - ) ... - )/2; -end -freeEnergyFullPart(fepIdx) = sum(feAux.line4(:,1)); - - -fepIdx = fepIdx + 1; -for iSub = listSubjects - feAux.line4(iSub,2) = ... - -(posterior.dcmMean(iSub,nDcmParamsInfCon+1:end)-priors.hemMean)*... - (priors.hemSigma\... - (posterior.dcmMean(iSub,nDcmParamsInfCon+1:end)... - -priors.hemMean).'... - )/2; -end -freeEnergyFullPart(fepIdx) = sum(feAux.line4(:,2)); - - -%% line 5 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -counts.nSubjects*tapas_util_logdet(priors.hemSigma)/2; - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -counts.nSubjects*nDcmParamsInfHem*log(2*pi)/2; - - -%% line 6 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -counts.nSubjects*sum(gammaln(priors.noiseShape)); - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nSubjects*sum(priors.noiseShape.*log(priors.noiseInvScale)); - - -%% line 7 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(sum(bsxfun(@times, ... - priors.noiseShape-1, ... - (psi(posterior.noiseShape)-log(posterior.noiseInvScale)) ... - ))); - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -sum(sum(bsxfun(@times, ... - priors.noiseInvScale, ... - posterior.noiseShape./posterior.noiseInvScale ... - ))); - - -%% line 8 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -nDcmParamsInfCon*sum(priors.clustersTau./posterior.clustersTau)/2; - - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = 0; -for iCluster = 1:counts.nClusters - freeEnergyFullPart(fepIdx) = freeEnergyFullPart(fepIdx) + ... - -priors.clustersTau*posterior.clustersDeg(iCluster)/2 ... - *(posterior.clustersMean(iCluster,:)-priors.clustersMean) ... - *(posterior.clustersSigma(:,:,iCluster)\... - (posterior.clustersMean(iCluster,:)-priors.clustersMean).'); -end -%%% - - -%% line 9 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = 0; -for iCluster = 1:counts.nClusters - freeEnergyFullPart(fepIdx) = freeEnergyFullPart(fepIdx) + ... - -posterior.clustersDeg(iCluster)/2 ... - *trace(posterior.clustersSigma(:,:,iCluster)\... - priors.clustersSigma); -end -%%% - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -(priors.clustersDeg+nDcmParamsInfCon+2)/2*(... - sum(posterior.logDetClustersSigma) ... - -sum(sum(psi(bsxfun(@plus,... - posterior.clustersDeg,... - 1-(1:nDcmParamsInfCon)... - )/2))) ... - ); - - - -%% line 10 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -counts.nClusters*nDcmParamsInfCon*(nDcmParamsInfCon-1)/4*log(pi); -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -counts.nClusters*sum(gammaln((priors.clustersDeg+1-... - (1:nDcmParamsInfCon))/2)); - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nClusters*priors.clustersDeg*tapas_util_logdet(... - priors.clustersSigma)/2; - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nClusters*nDcmParamsInfCon*(nDcmParamsInfCon+2)/2*log(2); - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -counts.nClusters*nDcmParamsInfCon*log(2*pi)/2; -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nClusters*nDcmParamsInfCon*log(priors.clustersTau)/2; - - - -%% line 11 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(sum(bsxfun(@times,... - posterior.softAssign,... - psi(posterior.alpha.')-psi(sum(posterior.alpha))... - ))); - - -%% line 12 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +gammaln(sum(priors.alpha))-sum(gammaln(priors.alpha)); - - -%% line 13 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +(priors.alpha.'-1)*(psi(posterior.alpha)-psi(sum(posterior.alpha))); - - -%% line 14 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -sum(posterior.softAssign(:).*log(posterior.softAssign(:)+realmin)); - - -%% line 15 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(posterior.logDetPostDcmSigma)/2; -%%% neg dF - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nSubjects*nDcmParamsInfAll/2; - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nSubjects*nDcmParamsInfAll*log(2*pi)/2; - - - -%% line 16 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(gammaln(posterior.noiseShape(:))); - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -posterior.noiseShape(:).'*log(posterior.noiseInvScale(:)); - - -%% line 17 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -(posterior.noiseShape(:).'-1)*(psi(posterior.noiseShape(:))-... - log(posterior.noiseInvScale(:))); - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(posterior.noiseShape(:)); - - -%% line 18 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nClusters*nDcmParamsInfCon/2; - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +nDcmParamsInfCon/2*sum(posterior.clustersDeg); - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -counts.nClusters*nDcmParamsInfCon*(nDcmParamsInfCon+2)/2*log(2); - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nClusters*nDcmParamsInfCon/2*log(2*pi); -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +nDcmParamsInfCon/2*sum(log(1./posterior.clustersTau)); - - -%% line 19 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +(posterior.clustersDeg.'+nDcmParamsInfCon+2)/2*(... - posterior.logDetClustersSigma... - -sum(psi(bsxfun(@plus,... - posterior.clustersDeg,... - 1-(1:nDcmParamsInfCon)... - )/2),2)... - ); - - -%% line 20 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +counts.nClusters*nDcmParamsInfCon*(nDcmParamsInfCon-1)/4*log(pi); -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(sum(gammaln(bsxfun(@plus,... - posterior.clustersDeg,... - 1-(1:nDcmParamsInfCon)... - )/2))); - - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -posterior.clustersDeg.'*posterior.logDetClustersSigma/2; - - - -%% line 21 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - +sum(gammaln(posterior.alpha)); - -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -gammaln(sum(posterior.alpha)); - -%% line 22 -fepIdx = fepIdx + 1; -freeEnergyFullPart(fepIdx) = ... - -(posterior.alpha.'-1)*(psi(posterior.alpha)-psi(sum(posterior.alpha))); - - - -%% sum -freeEnergyFull = sum(freeEnergyFullPart); - - -%% aux -freeEnergyFullPart(61) = freeEnergyFullPart(2)... - +freeEnergyFullPart(13)... - +freeEnergyFullPart(34); - - - diff --git a/huge/tapas_huge_pack_params.m b/huge/tapas_huge_pack_params.m deleted file mode 100644 index a0a023b6..00000000 --- a/huge/tapas_huge_pack_params.m +++ /dev/null @@ -1,50 +0,0 @@ -%% [ packedParameters ] = tapas_huge_pack_params( unpackedParameters, paramList ) -% -% Transform vectorized parameters into a cell array format. -% -% INPUT: -% unpackedParameters - current value of parameters arranged as a -% vector -% paramList - supporting paramters -% -% OUTPUT: -% packedParameters - current value of parameters arranged as a cell -% array -% - -% -% Author: Sudhir Shankar Raman -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ packedParameters ] = tapas_huge_pack_params( unpackedParameters, paramList ) - -packedParameters{1} = reshape(unpackedParameters(... - 1:(paramList(3)*paramList(3))),paramList(3),paramList(3)); -nextIndex = paramList(3)*paramList(3)+1; -packedParameters{2} = reshape(unpackedParameters(nextIndex:... - (nextIndex+paramList(3)*paramList(4)-1)),paramList(3),paramList(4)); -nextIndex = nextIndex+paramList(3)*paramList(4); -packedParameters{3} = reshape(unpackedParameters(nextIndex:... - (nextIndex+paramList(4)*paramList(3)*paramList(3)-1)),paramList(3),... - paramList(3),paramList(4)); -nextIndex = nextIndex + paramList(3)*paramList(3)*paramList(4); -packedParameters{4} = reshape(unpackedParameters(nextIndex:... - (nextIndex+paramList(3)*paramList(3)*paramList(3)-1)),paramList(3),... - paramList(3),paramList(3)); -nextIndex = nextIndex + paramList(3)*paramList(3)*paramList(3); -packedParameters{5} = reshape(unpackedParameters(nextIndex:... - (nextIndex+paramList(3)*3-1)),paramList(3),3)'; - -return; \ No newline at end of file diff --git a/huge/tapas_huge_parse_inputs.m b/huge/tapas_huge_parse_inputs.m new file mode 100644 index 00000000..0bb6129e --- /dev/null +++ b/huge/tapas_huge_parse_inputs.m @@ -0,0 +1,62 @@ +function [ opts ] = tapas_huge_parse_inputs( opts, varargin ) +% Parse name-value pair type arguments into a struct. +% +% INPUTS: +% opts - Struct containing all valid names as field names and +% corresponding default values as field values. +% +% OPTIONAL INPUTS: +% Name-value pair arguments. +% +% OUTPUTS: +% opts - Struct containing name-value pair input arguments as fields. +% +% EXAMPLE: +% opts = TAPAS_HUGE_PARSE_INPUTS(struct('a',0,'b',1),'a',10) Specify +% 'a' and 'b' as valid property names and receive non-default value +% for 'a'. +% + +% Author: Yu Yao (yao@biomed.ee.ethz.ch) +% Copyright (C) 2019 Translational Neuromodeling Unit +% Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of TAPAS, which is released under the terms of the GNU +% General Public Licence (GPL), version 3. For further details, see +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. +% +% This software is intended for research only. Do not use for clinical +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% + +if isscalar(varargin) && iscell(varargin{1}) + varargin = varargin{1}; +end + +nIn = numel(varargin); +validNames = fieldnames(opts); + +for iIn = 1:2:nIn - 1 + + name = varargin{iIn}; + value = varargin{iIn + 1}; + + assert(ischar(name), 'TAPAPS:HUGE:Nvp:NoneChar', ... + 'Property name must be a character array.'); + assert(any(strcmpi(name, validNames)), 'TAPAPS:HUGE:Nvp:InvalidName', ... + '"%s" is not a valid property name.', name); + + opts.(lower(name)) = value; + +end + +end + diff --git a/huge/tapas_huge_plot.m b/huge/tapas_huge_plot.m deleted file mode 100644 index 26c1229b..00000000 --- a/huge/tapas_huge_plot.m +++ /dev/null @@ -1,130 +0,0 @@ -%% [ ] = tapas_huge_plot( DCM, DcmResults ) -% -% Generate simple graphical summary of inversion result. -% -% INPUT: -% DCM - cell array of DCM in SPM format -% DcmResults - struct used for storing the results from VB. Output of -% tapas_huge_invert.m -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ ] = tapas_huge_plot( DCM, DcmResults ) -%% clustering -if ~isfield(DCM,'listBoldResponse') - try - DcmInfo = tapas_huge_import_spm(DCM); - catch err - disp('tapas_huge_plot: Unsupported format.'); - disp('Use cell array of DCM in SPM format as first input.'); - rethrow(err); - end -else - DcmInfo = DCM; -end - - -figure -if DcmResults.maxClusters > 1 - % subject assignment - subplot(2,2,1) - bar(DcmResults.posterior.softAssign,'stacked'); - axis([0 DcmInfo.nSubjects+1 0 1]) - title('assignments') - ylabel('q_{nk}') - xlabel('subject index') - - % cluster - subplot(2,2,2) - hold on - line([0 DcmInfo.nConnections+1],[0 0],'color','k') - % plot posterior cluster mean and 95% marginal credible intervals - for k = 1:DcmResults.maxClusters - xOffset = ((k-1)/max(1,(DcmResults.maxClusters-1)) - .5)/2; - clMean = DcmResults.posterior.clustersMean(k,:); - clStd = sqrt(diag(DcmResults.posterior.clustersSigma(:,:,k))'/... - (DcmResults.posterior.clustersTau(k)*... - (DcmResults.posterior.clustersDeg(k) - DcmInfo.nConnections + 1))); - s = tinv(1-0.025,DcmResults.posterior.clustersDeg(k)); - errorbar((1:DcmInfo.nConnections) + xOffset,clMean,s*clStd,'d') - end - xlim([0 DcmInfo.nConnections+1]) - title('clusters') - ylabel('\mu_k') - xlabel('connection index'); - - % DCM - % plot 25 samples from posterior over (noise free) BOLD response - subplot(2,1,2) - hold on - n = 1; % for first subject - postMean = DcmResults.posterior.dcmMean(n,:); - postStd = chol(DcmResults.posterior.dcmSigma(:,:,n)); - for iSmp = 1:25 - % draw a sample from posterior over DCM parameters - postSmp = postMean + randn(size(postMean))*postStd; - tmp = zeros(1,DcmInfo.nParameters); - tmp(DcmInfo.connectionIndicator) = postSmp(1:DcmInfo.nConnections); - tmp(end-3*DcmInfo.nStates+1:end-DcmInfo.nStates+1) = ... - postSmp(DcmInfo.nConnections+1:end); - % predict BOLD response for current sample - pred = tapas_huge_bold(tmp,DcmInfo,n); - plot(pred(:),'b') - end - plot(DcmInfo.listBoldResponse{n}(:),'k') - % plot ground truth if available - try - tmp = DcmInfo.listParameters{n}; - pred = tapas_huge_bold([tmp{1}(:);tmp{2}(:);tmp{3}(:);tmp{4}(:);... - tmp{5}(:);]',DcmInfo,n); - plot(pred(:),'r') - catch - % omit ground truth - end - title('black: measurement - blue: posterior samples'); - ylabel('BOLD') - xlabel('sample index') - -else - %% empirical Bayes - subplot(2,1,1) - hold on - line([0 DcmInfo.nConnections+1],[0 0],'color','k') - k = 1; -% plot posterior cluster mean and 95% marginal credible intervals - clMean = DcmResults.posterior.clustersMean(k,:); - clStd = sqrt(diag(DcmResults.posterior.clustersSigma(:,:,k))'/... - (DcmResults.posterior.clustersTau(k)*... - (DcmResults.posterior.clustersDeg(k) - DcmInfo.nConnections + 1))); - s = tinv(1-0.025,DcmResults.posterior.clustersDeg(k)); - errorbar(1:DcmInfo.nConnections,clMean,s*clStd,'d') - ylabel('\mu_k') - xlabel('connection index'); - title('Empirical Bayes') - - % boxplot of MAP estimates of DCM parameters - subplot(2,1,2) - hold on - line([0 DcmInfo.nConnections+1],[0 0],'color','k') - boxplot(DcmResults.posterior.dcmMean) - ylabel('\mu_n') - xlabel('connection index'); - -end - -end - diff --git a/huge/tapas_huge_predict.m b/huge/tapas_huge_predict.m deleted file mode 100644 index 2f65e5ac..00000000 --- a/huge/tapas_huge_predict.m +++ /dev/null @@ -1,54 +0,0 @@ -%% function [response,x,s,f,v,q] = tapas_huge_predict( dcmParamsInf, dcmParamsDefault, idxDcmParamsInf, ~, DcmInfo, iSubject) -% -% Wrapper for the function for tapas_huge_bold.m -% -%-------------------------------------------------------------------------------------- -% INPUT: -% dcmParamsInf - values of DCM parameters being inferred -% dcmParamsDefault - full DCM parameters vector containing default -% values for remaining parameters -% idxDcmParamsInf - indices of DCM parameters being inferred -% idxDiagA - indices of self connections (dummy) -% DcmInfo - struct containing DCM model specification and -% BOLD time series. -% iSubject - subject index -% -% -%--------------------------------------------------------------------------------------- -% OUTPUT: -% response - matrix of predicted response for each region -% (column-wise) -% x - time series of neuronal states -% s - time series of vasodilatory signal -% f1 - time series of flow -% v1 - time series of blood volume -% q1 - time series of deoxyhemoglobin content. -% -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [response,x,s,f,v,q] = tapas_huge_predict( dcmParamsInf, dcmParamsDefault, idxDcmParamsInf, ~, DcmInfo, iSubject) - -dcmParamsDefault(idxDcmParamsInf) = dcmParamsInf; - - -[response,x,s,f,v,q] = tapas_huge_bold(dcmParamsDefault,... - DcmInfo,... - iSubject); -response = response(:); - - diff --git a/huge/tapas_huge_property_names.m b/huge/tapas_huge_property_names.m new file mode 100644 index 00000000..8cf9d7f4 --- /dev/null +++ b/huge/tapas_huge_property_names.m @@ -0,0 +1,131 @@ +%% List of name-value pair arguments accepted by tapas_Huge() and estimate() +% +% NAME: Confounds +% VALUE: double array +% DESCRIPTION: Specify confounds for group-level analysis (e.g. age or sex) +% as double array with one row per subject and one column per +% confound. Note: This property can only be used in +% combination with the Dcm property. +% WARNING: This feature is experimental. +% +% NAME: ConfoundsVariant +% VALUE: 'none' | 'global' | 'cluster' (default: 'global' if +% confounds specified, 'none' otherwise) +% DESCRIPTION: Determines how confounds enter model. 'none': Confounds are +% not used. 'global': Confounds enter global regression +% (variant 1). 'cluster': Confounds enter cluster-specific +% regression (variant 2). +% +% NAME: Dcm +% VALUE: cell array of DCM structs in SPM format +% DESCRIPTION: Specify DCM structure and BOLD time series for all subjects +% as cell array with one DCM struct in SPM format per subject. +% +% NAME: K +% VALUE: positive integer (default: 1) +% DESCRIPTION: Number of clusters. +% +% NAME: Method +% VALUE: 'VB' +% DESCRIPTION: Name of inversion method specified as character array. VB: +% variational Bayes +% +% NAME: NumberOfIterations +% VALUE: positive integer (default: 999) +% DESCRIPTION: Maximum number of iterations of VB scheme. +% +% NAME: OmitFromClustering +% VALUE: array of logical | struct with fields a, b, c and d +% (default: empty struct) +% DESCRIPTION: Select DCM parameters to exclude from clustering. Parameters +% excluded from clustering will still be estimated, but under +% a static Gaussian prior. If input is array, it will be +% treated as the a field in a struct. Missing fields will be +% treated the same as arrays of false. Note: This property can +% only be used in combination with the Dcm property. +% +% NAME: PriorClusterMean +% VALUE: 'default' | row vector of double +% DESCRIPTION: Prior cluster mean. Scalar input will be expanded into +% vector. Default: [0, ... ,0]. +% +% NAME: PriorClusterVariance +% VALUE: 'default' | symmetric, positive definite matrix of double +% DESCRIPTION: Prior mean of cluster covariances. Must be a symmetric, +% positive definite matrix. Scalar input will be expanded into +% diagonal matrix. Default: diag(..., 0.01, ...). +% +% NAME: PriorDegree +% VALUE: 'default' | positive double +% DESCRIPTION: nu_0 determines the prior precision of the cluster +% covariance. For VB, this is the degrees of freedom of the +% inverse-Wishart. Default: 100. +% +% NAME: PriorVarianceRatio +% VALUE: 'default' | positive double +% DESCRIPTION: Ratio tau_0 between prior mean cluster covariance and prior +% covariance over cluster mean. Prior covariance over cluster +% mean equals prior cluster covariance divided tau_0. Default: +% 0.1. +% +% NAME: Randomize +% VALUE: bool (default: false) +% DESCRIPTION: If true, starting values for subject level DCM parameter +% estimates are randomized. +% +% NAME: SaveTo +% VALUE: character array +% DESCRIPTION: Location for saving results consisting of path name and +% optional file name. Path name must end on file separator and +% point to an existing directory. If file name is not +% specified, it is set to 'huge' followed by date and time. +% +% NAME: Seed +% VALUE: double | cell array of double and rng name | random number +% generator seed obtained with rng() command +% DESCRIPTION: Seed for random number generator. +% +% NAME: StartingValueDcm +% VALUE: 'prior' | 'spm' | double array (default: 'prior') +% DESCRIPTION: Starting values for subject-level DCM parameter estimates. +% 'prior' uses prior cluster mean for all subjects. 'spm' uses +% values supplied in the 'Ep' field of the SPM DCM structs. +% Use a double array with number of rows equal to number of +% subjects to specify custom starting values. +% +% NAME: StartingValueGmm +% VALUE: 'prior' | double array (default: 'prior') +% DESCRIPTION: Starting values for cluster-level DCM parameter estimates. +% 'prior' uses prior cluster mean for all clusters. Use a +% double array with number of rows equal to number of clusters +% to specify custom starting values. +% +% NAME: Tag +% VALUE: character array +% DESCRIPTION: Model description +% +% NAME: Verbose +% VALUE: bool (default: false) +% DESCRIPTION: Activate/deactivate command line output. +% +% + +% Author: Yu Yao (yao@biomed.ee.ethz.ch) +% Copyright (C) 2019 Translational Neuromodeling Unit +% Institute for Biomedical Engineering, +% University of Zurich and ETH Zurich. +% +% This file is part of TAPAS, which is released under the terms of the GNU +% General Public Licence (GPL), version 3. For further details, see +% . +% +% This software is provided "as is", without warranty of any kind, express +% or implied, including, but not limited to the warranties of +% merchantability, fitness for a particular purpose and non-infringement. +% +% This software is intended for research only. Do not use for clinical +% purpose. Please note that this toolbox is under active development. +% Considerable changes may occur in future releases. For support please +% refer to: +% https://github.com/translationalneuromodeling/tapas/issues +% \ No newline at end of file diff --git a/huge/tapas_huge_simulate.m b/huge/tapas_huge_simulate.m deleted file mode 100644 index acf5f9ca..00000000 --- a/huge/tapas_huge_simulate.m +++ /dev/null @@ -1,214 +0,0 @@ -%% [ DcmInfo ] = tapas_huge_simulate( options ) -% -% Generates a synthetic HUGE dataset for DCM for fMRI. -% -% INPUT: -% options - struct containing information on parameters of the -% generative model to be used in the simulation. -% (see tapas_huge_generate_examples.m for an example) -% -% OUTPUT: -% DcmInfo - struct containing DCM model specification and BOLD time -% series. -% -% REFERENCE: -% -% Yao Y, Raman SS, Schiek M, Leff A, Frssle S, Stephan KE (2018). -% Variational Bayesian Inversion for Hierarchical Unsupervised Generative -% Embedding (HUGE). NeuroImage, 179: 604-619 -% -% https://doi.org/10.1016/j.neuroimage.2018.06.073 -% - -% Author: Yu Yao (yao@biomed.ee.ethz.ch) -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ DcmInfo ] = tapas_huge_simulate( options ) -%% process input -rngSeed = rng(); -% compile integrator -tapas_huge_compile(); - -K = numel(options.N_k); % number of clusters -N = sum(options.N_k); % number of subjects -R = options.R; % number of regions - -% inputs -u = options.input.u; -if iscell(u) - u = u(:)'; - assert(length(u)==N,... - 'TAPAS:HUGE:NumberOfInputs',... - 'Number of input arrays inconsistent with number of subjects'); - L = unique(cellfun(@size,u,repmat({2},size(u)))); - assert(isscalar(L(:)),... - 'TAPAS:HUGE:InputSize',... - 'Size of inputs inconsistent'); -else - L = size(u,2); - u = repmat({u},1,N); -end - -% number of parameters for fully connected DCM -nParametersConnect = R^2 + R^2*L + R*L + R^3; -nParameters = nParametersConnect + 3*R; - -% clusters -mu_k = zeros(K,nParametersConnect); -mu_k(:,options.mu_k.idx) = options.mu_k.value; -p_c = numel(options.mu_k.idx); -assert(size(mu_k,1)==K,... - 'TAPAS:HUGE:NumberOfClusters',... - 'Number of mean vectors inconsistent with number of clusters'); -sigma_k = options.sigma_k; -if isscalar(sigma_k) - sigma_k = repmat(sigma_k,size(options.mu_k.value)); -elseif iscolumn(sigma_k) - sigma_k = repmat(sigma_k,1,p_c); -end -assert(all(size(sigma_k)==size(options.mu_k.value)),... - 'TAPAS:HUGE:NumberOfClusters',... - 'Size of sigma_k inconsistent with mu_k'); - -% hemodynamics -p_h = 2*R + 1; -mu_h = options.mu_h; -assert(all(size(mu_h)==[1,p_h]),... - 'TAPAS:HUGE:HemodynamicsSize',... - 'Size of mu_h must be 1x2*R+1'); -sigma_h = options.sigma_h; -if isscalar(sigma_h) - sigma_h = repmat(sigma_h,1,p_h); -end -assert(all(size(sigma_h)==size(mu_h)),... - 'TAPAS:HUGE:HemodynamicsSize',... - 'mu_h and sigma_h must have same size'); - -% observation noise -snrBold = options.snr; - -%% assemble DcmInfo struct -DcmInfo = struct(); - -DcmInfo.nSubjects = N; -DcmInfo.nStates = R; -DcmInfo.nParameters = nParameters; -DcmInfo.connectionIndicator = options.mu_k.idx; -DcmInfo.noConnectionIndicator = 1:nParametersConnect; -DcmInfo.noConnectionIndicator(DcmInfo.connectionIndicator) = []; -DcmInfo.nNoConnections = nParametersConnect - p_c; -DcmInfo.nConnections = p_c; - -dummy = zeros(1,nParameters); -dummy(:,options.mu_k.idx) = 1; -dummy = tapas_huge_pack_params(dummy,[NaN NaN R L]); -DcmInfo.adjacencyA = logical(dummy{1}); -DcmInfo.adjacencyB = logical(dummy{3}); -DcmInfo.adjacencyC = logical(dummy{2}); -DcmInfo.adjacencyD = logical(dummy{4}); -DcmInfo.dcmTypeB = any(DcmInfo.adjacencyB(:)); -DcmInfo.dcmTypeD = any(DcmInfo.adjacencyD(:)); - -DcmInfo.hemParam = options.hemParam; -DcmInfo.nInputs = L; -DcmInfo.nTime = cellfun(@size,u,repmat({1},1,length(u))); -DcmInfo.listU = u; -% ratio of sampling rate between input and BOLD -if isscalar(options.input.trSteps) - DcmInfo.trSteps = repmat(options.input.trSteps,1,N); -else - assert(numel(options.input.trSteps)==N,... - 'TAPAS:HUGE:NumberOfInputs',... - 'Length of trSteps inconsistent with number of subjects'); - DcmInfo.trSteps = options.input.trSteps; -end -% repetition time -if isscalar(options.input.trSeconds) - DcmInfo.trSeconds = repmat(options.input.trSeconds,1,N); -else - assert(numel(options.input.trSeconds)==N,... - 'TAPAS:HUGE:NumberOfInputs',... - 'Length of trSeconds inconsistent with number of subjects'); - DcmInfo.trSeconds = options.input.trSeconds; -end -% sampling time step of input -DcmInfo.timeStep = DcmInfo.trSeconds./DcmInfo.trSteps; - -DcmInfo.ClusterParameters.clusterMeans = mu_k; -DcmInfo.ClusterParameters.clusterStds = sigma_k; -DcmInfo.ClusterParameters.clusterCovariances = zeros(p_c,p_c,K); -for k = 1:K - DcmInfo.ClusterParameters.clusterCovariances(:,:,k) = ... - diag(sigma_k(k,:).^2); -end - - - - -%% generate subject-specific BOLD responses - -listN = cumsum(options.N_k); -DcmInfo.listBoldResponse = cell(1,N); -DcmInfo.listParameters = cell(N,2); -DcmInfo.listResponseTimeIndices = cell(1,N); -DcmInfo.trueLabels = zeros(1,N); - -for n = 1:N - - k = nnz(n>listN) + 1; % current cluster index - DcmInfo.trueLabels(n) = k; - DcmInfo.listResponseTimeIndices{n} = ... - DcmInfo.trSteps(n):DcmInfo.trSteps(n):DcmInfo.nTime(n); - - bStable = false; - while ~bStable - - % draw subject-specific parameter vector - mu_n = zeros(1,nParameters); - mu_n(DcmInfo.connectionIndicator) = randn(1,p_c).*sigma_k(k,:); - mu_n(nParametersConnect+1:nParametersConnect+p_h) = ... - randn(1,p_h).*sigma_h; - mu_n = mu_n + [mu_k(k,:),mu_h,zeros(1,R-1)]; - - % simulate response BOLD - response = tapas_huge_bold(mu_n,DcmInfo,n); - - % check if DCM stable - dcmParameters = tapas_huge_pack_params(mu_n,[NaN NaN R L]); - responseSum = sum(response(:)); - bStable = ~isnan(responseSum) && ~isinf(abs(responseSum)) && ... - max(svd(dcmParameters{1}))<1; - - end - - % noise model - lambda_nr = snrBold./var(response,0,1); - DcmInfo.listBoldResponse{n} = response + ... - bsxfun(@rdivide,randn(size(response)),sqrt(lambda_nr)); - DcmInfo.listBoldResponse{n} = DcmInfo.listBoldResponse{n}(:); - - DcmInfo.listParameters{n,1} = dcmParameters; - DcmInfo.listParameters{n,2} = lambda_nr; - -end - - -% save rng seed -DcmInfo.rngSeed = rngSeed; -DcmInfo.tag = 'TAPAS:HUGE:v0'; - - -end - diff --git a/huge/tapas_util_logdet.m b/huge/tapas_util_logdet.m deleted file mode 100644 index 306ecaec..00000000 --- a/huge/tapas_util_logdet.m +++ /dev/null @@ -1,33 +0,0 @@ -%% [ ld ] = tapas_util_logdet( A ) -% -% Numerical stable calculation of log-determinant for positive-definite -% matrix. -% -% INPUT: -% A - positive definite matrix. -% -% OUTPUT: -% ld - log(det(A)) -% - -% Author: Sudhir Shankar Raman -% Copyright (C) 2018 Translational Neuromodeling Unit -% Institute for Biomedical Engineering, -% University of Zurich and ETH Zurich. -% -% This file is part of TAPAS, which is released under the terms of the GNU -% General Public Licence (GPL), version 3. For further details, see -% . -% -% This software is intended for research only. Do not use for clinical -% purpose. Please note that this toolbox is in an early stage of -% development. Considerable changes are planned for future releases. For -% support please refer to: -% https://github.com/translationalneuromodeling/tapas/issues -% -function [ld] = tapas_util_logdet(A) - -U = chol(A); -ld = 2*sum(log(diag(U))); - -return; \ No newline at end of file diff --git a/misc/log_tapas.txt b/misc/log_tapas.txt index 5ce21a68..da54dd5d 100644 --- a/misc/log_tapas.txt +++ b/misc/log_tapas.txt @@ -1,3 +1,4 @@ +3.2.0 https://www.tapas.tnu-zurich.com/examples_v3.2.0.zip 00137f7a332be031fdb2a748a8b1592c 3.1.0 https://www.tapas.tnu-zurich.com/examples_v3.1.0.zip 05cbb1df4eaeab59f88d2c6ad7a87946 3.0.1 https://www.tapas.tnu-zurich.com/examples_v3.0.0.zip 81ba0d8cb1267880e0bab5c45a47b18a 3.0.0 https://www.tapas.tnu-zurich.com/examples_v3.0.0.zip 81ba0d8cb1267880e0bab5c45a47b18a diff --git a/misc/tapas_version.m b/misc/tapas_version.m index 840ff8c2..f2827ed0 100644 --- a/misc/tapas_version.m +++ b/misc/tapas_version.m @@ -26,12 +26,15 @@ tapas_print_logo(); fprintf(1, '\n\nVersion %s.%s.%s\n', version{:}); fprintf(1, 'In your citation please include the current version.\n'); - fprintf(1, 'Please cite the corresponding publications according to the toolboxes used:\n') - fprintf(1, 'PhysIO: https://www.ncbi.nlm.nih.gov/pubmed/27832957\n') - fprintf(1, 'HGF: https://www.ncbi.nlm.nih.gov/pubmed/21629826\n') + fprintf(1, 'Please cite the corresponding publications according to the toolboxes used:\n'); + fprintf(1, 'PhysIO: https://www.ncbi.nlm.nih.gov/pubmed/27832957\n'); + fprintf(1, 'HGF: https://www.ncbi.nlm.nih.gov/pubmed/21629826\n'); fprintf(1, ' https://www.ncbi.nlm.nih.gov/pubmed/25477800\n'); fprintf(1, 'MPDCM: https://www.ncbi.nlm.nih.gov/pubmed/26384541\n'); fprintf(1, 'SERIA: https://www.ncbi.nlm.nih.gov/pubmed/28767650\n'); + fprintf(1, 'HUGE: https://www.ncbi.nlm.nih.gov/pubmed/29964187\n'); + fprintf(1, 'rDCM: https://www.ncbi.nlm.nih.gov/pubmed/29807151\n'); + fprintf(1, ' https://www.ncbi.nlm.nih.gov/pubmed/28259780\n'); end end