From af934b43618ad446c251811cc08b5221d26e7aff Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Thu, 7 Apr 2022 13:12:16 +0200 Subject: [PATCH 01/11] Add common tools for MATLAB and Python testing Common base for the matlab and python testing branch to minimize merge conflicts down the line --- .gitignore | 11 ++-- .pre-commit-config.yaml | 7 +++ CONTRIBUTING.md | 56 +++++++++++++++++++ Makefile | 116 ++++++++++++++++++++++++++++++++++++++++ README.md | 35 ++++++------ requirements_dev.txt | 11 ++++ 6 files changed, 217 insertions(+), 19 deletions(-) create mode 100644 .pre-commit-config.yaml create mode 100644 CONTRIBUTING.md create mode 100644 Makefile create mode 100644 requirements_dev.txt diff --git a/.gitignore b/.gitignore index 7c82424..7bd97a6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,10 @@ -# +# vscode files .vscode/ +# test related +tests/outputs/ -# data folder +# data folders data/ # example output folders @@ -13,6 +15,9 @@ figures/ report.html ##*.png +# pyenv +Pipfile + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] @@ -54,6 +59,7 @@ pip-delete-this-directory.txt htmlcov/ .tox/ .coverage +coverage_html .coverage.* .cache nosetests.xml @@ -153,4 +159,3 @@ codegen/ # Octave session info octave-workspace - diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..e9f04ad --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,7 @@ +# See https://pre-commit.com for more information +# See https://pre-commit.com/hooks.html for more hooks +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.1.0 + hooks: + - id: check-yaml diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..baaa794 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,56 @@ +# CONTRIBUTING + +Information for anyone who would like to contribute to this repository. + +## Repository map + +```bash +├── .git +├── .github +│ └── workflows # Github continuous integration set up +├── examples +│ ├── data +│ ├── example1outputs +│ ├── example2outputs +├── glmsingle # Python implementation +│ ├── cod +│ ├── design +│ ├── gmm +│ ├── hrf +│ ├── ols +│ ├── ssq +│ └── utils +├── matlab # Matlab implementation +│ ├── examples +│ ├── fracridge +│ └── utilities +└── tests # Python and Matlab tests + └── data + +``` + +## Generic + +### Makefile + +### pre-commit + +## Matlab + +### Style guide + +### Tests + +#### Demos + +### Continuous integration + +## Python + +### Style guide + +### Tests + +#### Demos + +### Continuous integration \ No newline at end of file diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..75e39ce --- /dev/null +++ b/Makefile @@ -0,0 +1,116 @@ +.DEFAULT_GOAL := help + +BROWSER := python -c "$$BROWSER_PYSCRIPT" + +# TODO make more general to use the local matlab version +MATLAB = /usr/local/MATLAB/R2017a/bin/matlab +MATLAB_ARG = -nodisplay -nosplash -nodesktop + +define BROWSER_PYSCRIPT +import os, webbrowser, sys + +from urllib.request import pathname2url + +webbrowser.open("file://" + pathname2url(os.path.abspath(sys.argv[1]))) +endef +export BROWSER_PYSCRIPT + +# determines what "make help" will show +define PRINT_HELP_PYSCRIPT +import re, sys + +for line in sys.stdin: + match = re.match(r'^([a-zA-Z_-]+):.*?## (.*)$$', line) + if match: + target, help = match.groups() + print("%-20s %s" % (target, help)) +endef +export PRINT_HELP_PYSCRIPT + +################################################################################ +# GENERIC +.PHONY: help clean clean-test lint install_dev + +help: + @python -c "$$PRINT_HELP_PYSCRIPT" < $(MAKEFILE_LIST) +clean: clean-build clean-pyc clean-test ## remove all build, test, coverage artifacts + +clean-test: ## remove test and coverage artifacts + rm -rf .tox/ + rm -rf .coverage + rm -rf htmlcov/ + rm -rf .pytest_cache + rm -rf tests/data + rm -rf test/outputs + +install_dev: ## install for both matlab and python developpers + pip install -e . + pip install -r requirements_dev.txt + +lint: lint/black lint/flake8 lint/miss_hit ## check style + +test: test-matlab test-python + +tests/data/nsdcoreexampledataset.mat: + mkdir tests/data + curl -fsSL --retry 5 -o "tests/data/nsdcoreexampledataset.mat" https://osf.io/k89b2/download + +################################################################################ + +################################################################################ +# MATLAB + +.PHONY: lint/miss_hit + +lint/miss_hit: ## lint and checks matlab code + mh_style --fix tests && mh_metric --ci tests && mh_lint tests + +test-matlab: tests/data/nsdcoreexampledataset.mat + $(MATLAB) $(MATLAB_ARG) -r "run_tests; exit()" + +coverage-matlab: test-matlab + $(BROWSER) coverage_html/index.html + +################################################################################ + +################################################################################ +# PYTHON + +.PHONY: clean-build clean-pyc coverage-python install lint/flake8 lint/black + +clean-build: ## remove build artifacts + rm -fr build/ + rm -fr dist/ + rm -fr .eggs/ + find . -name '*.egg-info' -exec rm -fr {} + + find . -name '*.egg' -exec rm -f {} + + +clean-pyc: ## remove Python file artifacts + find . -name '*.pyc' -exec rm -f {} + + find . -name '*.pyo' -exec rm -f {} + + find . -name '*~' -exec rm -f {} + + find . -name '__pycache__' -exec rm -fr {} + + +lint/flake8: ## check style with flake8 + flake8 tests +lint/black: ## check style with black + black tests + +test-python: tests/data/nsdcoreexampledataset.mat ## run tests quickly with the default Python + pytest + +test-notebooks: + pytest --nbmake --nbmake-timeout=3000 "./examples" +test-all: ## run tests on every Python version with tox + tox + +coverage-python: ## check code coverage quickly with the default Python + coverage run --source glmsingle -m pytest + coverage report -m + coverage html + $(BROWSER) htmlcov/index.html + +install: clean ## install the package to the active Python's site-packages + python setup.py install + +################################################################################ \ No newline at end of file diff --git a/README.md b/README.md index 1cc319e..336ed01 100644 --- a/README.md +++ b/README.md @@ -33,20 +33,6 @@ This will also clone [`fracridge`](https://github.com/nrdg/fracridge) as a submo To use the GLMsingle toolbox, add it and `fracridge` to your MATLAB path by running the `setup.m` script. -## Example scripts - -We provide a number of example scripts that demonstrate usage of GLMsingle. You can browse these example scripts here: - -(Python Example 1 - event-related design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/examples/example1.html - -(Python Example 2 - block design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/examples/example2.html - -(MATLAB Example 1 - event-related design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/matlab/examples/example1preview/example1.html - -(MATLAB Example 2 - block design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/matlab/examples/example2preview/example2.html - -If you would like to run these example scripts, the Python versions are available in `/GLMsingle/examples`, and the MATLAB versions are available in `/GLMsingle/matlab/examples`. Each notebook contains a full walkthrough of the process of loading an example dataset and design matrix, estimating neural responses using GLMsingle, estimating the reliability of responses at each voxel, and comparing those achieved via GLMsingle to those achieved using a baseline GLM. - ## Python To install: @@ -63,12 +49,25 @@ Running the demos requires: pip install jupyterlab ``` -Code dependencies: see requirements.txt +Code dependencies: see [requirements.txt](./requirements.txt) Notes: -* Please note that GLMsingle is not (yet) compatible with Python 3.9 (due to an incompatibility between scikit-learn and Python 3.9). Please use Python 3.8 or earlier. * Currently, numpy has a 4GB limit for the pickle files it writes; thus, GLMsingle will crash if the file outputs exceed that size. One workaround is to turn off "disk saving" and instead get the outputs of GLMsingle in your workspace and save the outputs yourself to HDF5 format. +## Example scripts + +We provide a number of example scripts that demonstrate usage of GLMsingle. You can browse these example scripts here: + +(Python Example 1 - event-related design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/examples/example1.html + +(Python Example 2 - block design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/examples/example2.html + +(MATLAB Example 1 - event-related design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/matlab/examples/example1preview/example1.html + +(MATLAB Example 2 - block design) https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/matlab/examples/example2preview/example2.html + +If you would like to run these example scripts, the Python versions are available in `/GLMsingle/examples`, and the MATLAB versions are available in `/GLMsingle/matlab/examples`. Each notebook contains a full walkthrough of the process of loading an example dataset and design matrix, estimating neural responses using GLMsingle, estimating the reliability of responses at each voxel, and comparing those achieved via GLMsingle to those achieved using a baseline GLM. + ## Additional information For additional information, please visit the Wiki page associated with this @@ -80,6 +79,10 @@ If you use GLMsingle in your research, please cite the following paper: * [Prince, J.S., Charest, I., Kurzawski, J.W., Pyles, J.A., Tarr, M.J., Kay, K.N. GLMsingle: a toolbox for improving single-trial fMRI response estimates. bioRxiv (2022).](https://www.biorxiv.org/content/10.1101/2022.01.31.478431v1) +## Contributing + +If you want to contribute to GLMsingle see the [contributing](./CONTRIBUTING.md) documentation to help you know what is where and how to set things up. + ## Change history * 2021/10/12 - Version 1.0 of GLMsingle is now released. A git tag has been added to the repo. diff --git a/requirements_dev.txt b/requirements_dev.txt new file mode 100644 index 0000000..8bf205e --- /dev/null +++ b/requirements_dev.txt @@ -0,0 +1,11 @@ +# Matlab dev +miss_hit + +# Python dev +flake8 +tox +coverage +pytest +black +nbmake +pytest-cov From ada50fd400dee78bd539946bbddd6b7b77dc9d8b Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 9 Apr 2022 11:04:19 +0200 Subject: [PATCH 02/11] initial sphinx set up with sphinx-quickstart --- docs/Makefile | 20 ++++++++++++++++ docs/make.bat | 35 +++++++++++++++++++++++++++ docs/source/conf.py | 55 +++++++++++++++++++++++++++++++++++++++++++ docs/source/index.rst | 20 ++++++++++++++++ requirements_dev.txt | 4 ++++ 5 files changed, 134 insertions(+) create mode 100644 docs/Makefile create mode 100644 docs/make.bat create mode 100644 docs/source/conf.py create mode 100644 docs/source/index.rst diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d0c3cbf --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..dc1312a --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 0000000..9a6a912 --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,55 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +# import os +# import sys +# sys.path.insert(0, os.path.abspath('.')) + + +# -- Project information ----------------------------------------------------- + +project = 'GLMsingle' +copyright = '2022, Prince, J.S., Charest, I., Kurzawski, J.W., Pyles, J.A., Tarr, M.J., Kay, K.N.' +author = 'Prince, J.S., Charest, I., Kurzawski, J.W., Pyles, J.A., Tarr, M.J., Kay, K.N.' + +# The full version, including alpha/beta/rc tags +release = '0.0.1' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = [] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'alabaster' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 0000000..a804c6d --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,20 @@ +.. GLMsingle documentation master file, created by + sphinx-quickstart on Sat Apr 9 11:02:56 2022. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to GLMsingle's documentation! +===================================== + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/requirements_dev.txt b/requirements_dev.txt index 8bf205e..1e19a65 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -9,3 +9,7 @@ pytest black nbmake pytest-cov + +# Documentation +sphinx +sphinxcontrib-matlabdomain \ No newline at end of file From 9b67d55e5e5489c28284221a94d6602e79706a74 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 9 Apr 2022 11:44:45 +0200 Subject: [PATCH 03/11] get basic set up for matlab doc --- docs/source/conf.py | 19 +++++++++++++------ docs/source/index.rst | 3 ++- docs/source/matlab.rst | 6 ++++++ docs/source/python.rst | 9 +++++++++ requirements_dev.txt | 3 ++- 5 files changed, 32 insertions(+), 8 deletions(-) create mode 100644 docs/source/matlab.rst create mode 100644 docs/source/python.rst diff --git a/docs/source/conf.py b/docs/source/conf.py index 9a6a912..564364a 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -10,9 +10,11 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # -# import os -# import sys -# sys.path.insert(0, os.path.abspath('.')) +import os +import sys +root_dir = os.path.abspath('../..') +sys.path.insert(0, root_dir) +matlab_src_dir = root_dir # -- Project information ----------------------------------------------------- @@ -30,8 +32,7 @@ # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. -extensions = [ -] +extensions = ['sphinxcontrib.matlab', 'sphinx.ext.autodoc'] # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -41,13 +42,19 @@ # This pattern also affects html_static_path and html_extra_path. exclude_patterns = [] +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# The master toctree document. +master_doc = 'index' + # -- Options for HTML output ------------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -html_theme = 'alabaster' +html_theme = 'sphinx_rtd_theme' # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, diff --git a/docs/source/index.rst b/docs/source/index.rst index a804c6d..49263ed 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -10,7 +10,8 @@ Welcome to GLMsingle's documentation! :maxdepth: 2 :caption: Contents: - + matlab + python Indices and tables ================== diff --git a/docs/source/matlab.rst b/docs/source/matlab.rst new file mode 100644 index 0000000..121eb92 --- /dev/null +++ b/docs/source/matlab.rst @@ -0,0 +1,6 @@ +MATLAB doc +********** + +.. mat:automodule:: matlab + +.. mat:autofunction:: GLMestimatesingletrial \ No newline at end of file diff --git a/docs/source/python.rst b/docs/source/python.rst new file mode 100644 index 0000000..c4ee38a --- /dev/null +++ b/docs/source/python.rst @@ -0,0 +1,9 @@ +Python doc +********** + +.. module:: glmsingle.glmsingle + +.. autoclass:: GLM_single +.. module:: glmsingle.glmsingle + +.. autoclass:: GLM_single \ No newline at end of file diff --git a/requirements_dev.txt b/requirements_dev.txt index 1e19a65..ea0a774 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -12,4 +12,5 @@ pytest-cov # Documentation sphinx -sphinxcontrib-matlabdomain \ No newline at end of file +sphinxcontrib-matlabdomain +sphinx_rtd_theme \ No newline at end of file From 1669dcbfa1a2897e08be618ebbace7c7507178f9 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 9 Apr 2022 11:54:42 +0200 Subject: [PATCH 04/11] make sure both matlab and python doc work --- docs/source/matlab.rst | 6 +++++- docs/source/python.rst | 8 ++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/docs/source/matlab.rst b/docs/source/matlab.rst index 121eb92..942d756 100644 --- a/docs/source/matlab.rst +++ b/docs/source/matlab.rst @@ -3,4 +3,8 @@ MATLAB doc .. mat:automodule:: matlab -.. mat:autofunction:: GLMestimatesingletrial \ No newline at end of file +.. mat:autofunction:: GLMestimatesingletrial + +.. mat:automodule:: matlab.utilities + +.. mat:autofunction:: calccod \ No newline at end of file diff --git a/docs/source/python.rst b/docs/source/python.rst index c4ee38a..de77628 100644 --- a/docs/source/python.rst +++ b/docs/source/python.rst @@ -4,6 +4,10 @@ Python doc .. module:: glmsingle.glmsingle .. autoclass:: GLM_single -.. module:: glmsingle.glmsingle + :members: + :undoc-members: + :show-inheritance: + +.. module:: glmsingle.cod.calc_cod -.. autoclass:: GLM_single \ No newline at end of file +.. autofunction:: calc_cod \ No newline at end of file From f9bd6081e4c38474a04e49dd983bc3d29b7f5dc3 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 9 Apr 2022 11:58:52 +0200 Subject: [PATCH 05/11] set up for read the docs --- .readthedocs.yml | 26 ++++++++++++++++++++++++++ requirements_dev.txt | 3 +++ 2 files changed, 29 insertions(+) create mode 100644 .readthedocs.yml diff --git a/.readthedocs.yml b/.readthedocs.yml new file mode 100644 index 0000000..c4e3841 --- /dev/null +++ b/.readthedocs.yml @@ -0,0 +1,26 @@ +# .readthedocs.yml +# Read the Docs configuration file +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Build documentation in the docs/ directory with Sphinx +sphinx: + configuration: docs/source/conf.py + builder: html + fail_on_warning: true + +# Build documentation with MkDocs +#mkdocs: +# configuration: mkdocs.yml + +# Optionally build your docs in additional formats such as PDF +formats: + - pdf + +# Optionally set the version of Python and requirements required to build your docs +python: + version: 3.6 + install: + - requirements: requirements_dev.txt diff --git a/requirements_dev.txt b/requirements_dev.txt index ea0a774..d80402b 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -1,3 +1,6 @@ +# Base requirements +requirements.txt + # Matlab dev miss_hit From c0b154da8dbd9f542ab5c5c53b344aca281af921 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 9 Apr 2022 12:03:47 +0200 Subject: [PATCH 06/11] fix typo --- requirements_dev.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements_dev.txt b/requirements_dev.txt index d80402b..d3e0966 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -1,5 +1,5 @@ # Base requirements -requirements.txt +-r requirements.txt # Matlab dev miss_hit From 604cba504a64c2372de484e87e029eee45707c0a Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 9 Apr 2022 12:06:49 +0200 Subject: [PATCH 07/11] do not fail build on warning in RTD --- .readthedocs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index c4e3841..583a43b 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -9,7 +9,7 @@ version: 2 sphinx: configuration: docs/source/conf.py builder: html - fail_on_warning: true + fail_on_warning: false # Build documentation with MkDocs #mkdocs: From fba98c21ff96f25070b1148e486d647cb496bf7b Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Thu, 28 Apr 2022 17:08:01 +0200 Subject: [PATCH 08/11] add Myst parsing to support markd own in the documentation --- docs/source/conf.py | 27 ++++++++++++++++----------- requirements_dev.txt | 3 ++- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 564364a..c0bd06d 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -12,19 +12,24 @@ # import os import sys -root_dir = os.path.abspath('../..') + +root_dir = os.path.abspath("../..") sys.path.insert(0, root_dir) matlab_src_dir = root_dir # -- Project information ----------------------------------------------------- -project = 'GLMsingle' -copyright = '2022, Prince, J.S., Charest, I., Kurzawski, J.W., Pyles, J.A., Tarr, M.J., Kay, K.N.' -author = 'Prince, J.S., Charest, I., Kurzawski, J.W., Pyles, J.A., Tarr, M.J., Kay, K.N.' +project = "GLMsingle" +copyright = """ +2022, Prince, J.S., Charest, I., Kurzawski, +J.W., Pyles, J.A., Tarr, M.J., Kay, K.N.""" +author = """ +Prince, J.S., Charest, I., Kurzawski, +J.W., Pyles, J.A., Tarr, M.J., Kay, K.N.""" # The full version, including alpha/beta/rc tags -release = '0.0.1' +release = "0.0.1" # -- General configuration --------------------------------------------------- @@ -32,10 +37,10 @@ # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. -extensions = ['sphinxcontrib.matlab', 'sphinx.ext.autodoc'] +extensions = ["sphinxcontrib.matlab", "sphinx.ext.autodoc", "myst_parser"] # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. @@ -43,10 +48,10 @@ exclude_patterns = [] # The name of the Pygments (syntax highlighting) style to use. -pygments_style = 'sphinx' +pygments_style = "sphinx" # The master toctree document. -master_doc = 'index' +master_doc = "index" # -- Options for HTML output ------------------------------------------------- @@ -54,9 +59,9 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -html_theme = 'sphinx_rtd_theme' +html_theme = "sphinx_rtd_theme" # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] \ No newline at end of file +# html_static_path = ["_static"] diff --git a/requirements_dev.txt b/requirements_dev.txt index d3e0966..4817aff 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -16,4 +16,5 @@ pytest-cov # Documentation sphinx sphinxcontrib-matlabdomain -sphinx_rtd_theme \ No newline at end of file +sphinx_rtd_theme +myst-parser From bc77081d99c090cbe210592b459733dd7e94cdcf Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Thu, 28 Apr 2022 17:08:23 +0200 Subject: [PATCH 09/11] add wiki to the doc and reformat it --- docs/source/index.rst | 3 +- docs/source/wiki.md | 293 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 295 insertions(+), 1 deletion(-) create mode 100644 docs/source/wiki.md diff --git a/docs/source/index.rst b/docs/source/index.rst index 49263ed..be19329 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -7,11 +7,12 @@ Welcome to GLMsingle's documentation! ===================================== .. toctree:: - :maxdepth: 2 + :maxdepth: 1 :caption: Contents: matlab python + wiki Indices and tables ================== diff --git a/docs/source/wiki.md b/docs/source/wiki.md new file mode 100644 index 0000000..10fdc4f --- /dev/null +++ b/docs/source/wiki.md @@ -0,0 +1,293 @@ +# WIKI + +## Basic information + +GLMsingle is introduced and described in the following pre-print: + +[Prince, J.S., Charest, I., Kurzawski, J.W., Pyles, J.A., Tarr, M.J., Kay, K.N. GLMsingle: a toolbox for improving single-trial fMRI response estimates. bioRxiv (2022).](https://doi.org/10.1101/2022.01.31.478431) + +GLMsingle is an analysis technique (an algorithm) for obtaining accurate +estimates of single-trial beta weights in fMRI data. + +If you have questions or discussion points, please use the Discussions feature +of this github repository, or alternatively, e-mail Kendrick (kay@umn.edu). If +you find a bug, please let us know by raising a Github issue. + +## Example scripts + +We provide a number of example scripts that demonstrate usage of GLMsingle. + +You can browse these example scripts here: + +- [Python - example 1](https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/examples/example1.html) +- [Python - example 2](https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/examples/example2.html) +- [MATLAB - example 1](https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/matlab/examples/example1preview/example1.html) +- [MATLAB - example 2](https://htmlpreview.github.io/?https://github.com/kendrickkay/GLMsingle/blob/main/matlab/examples/example2preview/example2.html) + +## Tips on usage + +### What's the deal with `sessionindicator`? + +If your experiment wants to analyze/combine data from multiple distinct scan +sessions (multiple scanner visits), there is the possibility that after basic +pre-processing, there may be substantial abrupt differences in percent BOLD +signal change across responses observed on different days. The +`sessionindicator` option allows you to tell the code how runs are grouped, and +it will internally use this information to z-score responses within sessions in +order to better estimate cross-validation performance. This normalization is +just done internally (under the hood) and does not propagate to the final +outputs. + +### Can I exert control over the noise pool voxels? + +The default behavior is to automatically select noise pool voxels based on +passing a simple signal intensity threshold (i.e. a simple mask that ignores +out-of-brain voxels) and on having very low amounts of BOLD variance related to +the experiment (i.e. a "R2" threshold). If you want to specifically +control the voxels that are used for the noise pool, you can achieve this by +setting `brainthresh` to [99 0] (which allows all voxels to pass the intensity +threshold), `brainR2` to 100 (which allows potentially all voxels in) and then +setting `brainexclude` to be all voxels that you do NOT want in the noise pool. + +## FAQ + +### What are the main things that GLMsingle does? + +The main components of GLMsingle include: + +1. a "library of HRFs" technique where an empirically derived set of HRF + timecourses (from the subjects in the NSD dataset) are used as potential HRFs + for each voxel in your dataset, +2. the GLMdenoise technique where data-driven nuisance regressors are obtained + and added into the GLM (using cross-validation for determining how many + nuisance regressors to add), and +3. ridge regression as a way to improve robustness of single-trial beta + estimates. The technique relies on heavy amounts of computation, but is + implemented in a relatively efficient manner, and could therefore serve as a + go-to tool for deriving beta estimates from many types of fMRI experimental + data, and especially for condition-rich designs with few repeats per + condition. + +### How does GLMsingle achieve denoising? + +We can consider each of the three components. + +1. Mismodeling the timecourse of a voxel can lead to suboptimal results (i.e. + results that are, in a sense, "noisy"); GLMsingle attempts to remedy this by + using a well regularized HRF selection strategy where a simple "index" is + learned for each voxel. +2. fMRI data suffer from very large amounts of spatially correlated noise (e.g. + due to head motion, physiological noise, etc.) --- by deriving these noise + sources from the data themselves, GLMsingle is able to provide some amount of + "modeling" out the noise. +3. Finally, fMRI designs often involve substantial overlap of the response + across successive trials. From a statistical perspective, this is going to + hurt estimation efficiency. Ridge regression is a method that induces some + amount of shrinkage bias to help improve out-of-sample generalization. + +Note that these three components are heterogenous with regards to the idea of +"fitting" the data. For the HRF component, the presumption is that we probably +have at least enough data to estimate an HRF index for each voxel; hence the +philosophy is to indeed fit the data. For the GLMdenoise component, the goal is +to try to add regressors to the model to better fit the data, while +acknowledging that at some point, overfitting is going to occur (and the +algorithm attempts to determine that point). For the ridge regression component, +things are a bit different. The point of ridge regression is to allow for the +possibility that due to measurement noise, single-trial beta estimates are +inaccurate and we want to actually limit the extent to which we fit the data. +Thus, somewhat counter-intuitively, ridge regression _increases_ the residuals +of the model fit. + +### In GLMsingle, the GLMdenoise and ridge regression (RR) components of the method require experimental conditions to repeat across runs. How should I think about whether this is appropriate for my experiment? + +In order to determine the hyperparameter settings, it is indeed necessary for +some amount of repeated trials to exist in the data that are given to GLMsingle. +It is true that, in a sense, any component of the fMRI data that does not repeat +for the repeated trials associated with a condition is thought of a "noise" by +the cross-validation process. Hence, there is some potential risk that +components of interest might be removed. However, for the most part, we believe +that the risk of this is actually minimal. This can be seen by thinking about +the nature of the effects of GLMdenoise and ridge regression. GLMdenoise is +simply trying to capture variance that appears to be spatially/globally +correlated across voxels in the dataset; hence, its flexibility is actually +quite limited. Ridge regression can only dampen the instabilities of beta +estimates that persists across temporally nearby trials and in a uniform manner +for all trials in the experiment; hence, its flexibility is also quite limited. + +### What are the metrics of quality that are relevant here? + +The philosophy behind GLMsingle lies primarily in cross-validation. Internal to +the algorithm is the use of cross-validation to determine the two main +"hyperparameters" that are used --- one is the number of nuisance regressors +(which are derived via PCA on a noise pool of voxels) to add into the GLM model, +and the other is the amount of ridge regularization to be applied to each voxel. +The idea is that the hyperparameters that best generalize to the unregularized +single-trial betas in some left-out data are the appropriate parameters to use. + +More generally, the quality metric being invoked here is "reliability" or +"reproducibility". Typically, in an experiment you have experimental conditions +that are presented multiple times. The assumption is that the signal induced by +a given condition should reproduce across different trials. + +The concept of noise ceiling (described at length in the NSD data paper) is +essentially a sophisticated method to quantify reliability. A simpler metric +that can assess reliability is simply taking the responses obtained for a voxel +and checking its reliability across different trials. For example, you could +split the data into two halves and correlate the signal estimated on one half +with the signal estimated on the other half. + +### I noticed that GLMsingle involves some internal cross-validation. Is this a problem for decoding-style analyses where we want to divide the data into a training set and a test set? + +It is true that GLMsingle makes use of all of the data that it is presented +with: for example, if you give it 10 runs of time-series data, its outputs are +dependent (in complex ways) on all of the data. However, we don't think that +this poses any major "circularity"-type problems. All that the algorithm knows +about are the onsets of your experimental conditions and some specification of +when you think that conditions repeat, insofar that you expect some component of +the response to be reproducible across trials. The repeats are used (in +leave-one-run-out cross-validation) to determine the setting of the +hyperparameters. GLMsingle has no access to your scientific hypotheses, and it +is hard to see how it could bias the single-trial beta estimates in favor of one +hypothesis over another. There is one exception, however. If the primary outcome +you are trying to demonstrate is that responses to the same condition are +reliable across runs, then that is, in a sense, exactly what the algorithm is +trying to achieve --- so in that scenario, you might not want to give all of the +data to GLMsingle. Instead, you could divide your data into two halves (for +example), independently give each half to GLMsingle, and then test your +hypothesis that indeed responses are reliable across the two halves. + +### I noticed that GLMsingle needs some conditions to repeat across runs. But my experiment is not fully balanced in this regard (or I have just a few repeats). Is this a problem? + +In general, we don't think that perfect balancing is that critical. The reason +lies in thinking about the nature of the technique --- the repeats are simply +used to set the hyperparameters (number of nuisance regressors; fractional +regularization amount for each voxel). In our experience, even just a handful of +repeats appears to be sufficient to robustly guide the estimation of the +hyperparameters. Thus, even if you can code only a subset of the trials in your +experiment as condition repeats, this may be sufficient to guide the +hyperparameter estimation, and the overall results from GLMsingle might be quite +good. (In fact, there are very few repeats in the NSD and BOLD5000 datasets, +which are the main two datasets demonstrated in the GLMsingle pre-print.) + +### How does R2 start becoming meaningful in the full FITHRF_GLMDENOISE_RR version? + +For a single-trial design matrix, note that in a sense the predictors are +extremely flexible and are happy to capture essentially all or almost all of the +variance in the time-series data from a voxel, even if that voxel contains no +actual signal. Thus, for `ASSUMEHRF` or `FITHRF` or `FITHRF_GLMDENOISE` type +models, the R2 values from these models are more or less meaningless +(everything looks "good"). However, the RR technique (as a direct consequence of +the fact that it will shrink betas to 0 in accordance to the cross-validated +generalizability of the single-trial beta estimates) will essentially leave +unperturbed the good voxels that have good SNR and aggressively shrink the bad +voxels with little or no SNR. As a consequence, the variance explained by the +beta estimates that are produced by ridge regression will be directly related to +the "goodness" of the voxel as determined by the cross-validation procedure. +Thus, the ridge regression results will have high R2 values for the +voxels that seem to have reproducible beta estimates across runs. + +### Why does ridge regression improve beta estimates? + +Ridge regression imposes a shrinkage prior on regression weights, with the exact +amount of regularization controlled by the hyperparameter. In regression +problems where there are correlations across the predictors, ordinary +least-squares estimates of regression weights are unbiased but can suffer from +high variance (i.e. affected by noise). By imposing some amount of shrinkage, +the ability of the estimated regression weights to be more generalizable to +unseen data can be improved. However, these are general statistical concepts. In +the specific case of fMRI designs involving closely spaced trials, there are +high levels of positive correlations between temporally adjacent trials. In a +sense, you can think of this situation as there being limited amounts of data +(statistical power) to disentangle the responses due to adjacent trials. +Ordinary least squares tries to estimate these responses but will generally lead +to noisy beta estimates with large magnitudes (either positive and/or negative). +Ridge regression allows the estimation to impose a prior that effectively +downweights the noisy and limited amounts of data that inform the distinction +between adjacent trials; we can view it as imposing some amount of "temporal +smoothing" (pushing the beta weights from adjacent trials to be more similar in +magnitude) with the goal of attaining an overall solution that has better +out-of-sample performance. Also, in GLMsingle, note that a different shrinkage +hyperparameter is selected for each voxel: this is important since voxels vary +substantially in their signal-to-noise ratios. + +### What are the units of the single-trial beta weights, and are there any interpretation issues? + +The default behavior of GLMsingle is to return beta weights in units of percent +signal change (by dividing by the mean signal intensity observed at each voxel +and multiplying by 100). However, there can be tricky issues having to do with +the selection of the HRF. If the HRF used in a GLM is incorrect or different +from the underlying HRF, there can easily be gross "gain" or "scale" differences +in the obtained percent signal change units of the betas. One of the things that +GLMsingle attempts to do for you is to estimate a more proper HRF timecourse (as +opposed to assuming a canonical HRF). It is the case that depending on how you +use GLMsingle, you may encounter some scale issues. For example, if you use +GLMsingle to analyze data from one scan session and then separately analyze data +from a second scan session, there will likely be some differences across the two +sets of results. First, there could be real (physiological) changes in the +percent signal change across different days. Second, each set of results may +have a different HRF identified and used for a given voxel; thus, the percent +signal change units in the betas for a given voxel will have a different +ultimate meaning for the two sets of results. It is an open question how to +"best" normalize or handle the betas from different scan sessions to more +accurately get at the underlying neural activity. However, a quick and dirty +heuristic is to normalize the betas, e.g., by z-scoring, or scaling, depending +on your overall scientific goals. + +### If the HRF changes from voxel to voxel (or area to area), doesn't that pose some interpretation difficulties or confounding issues? + +This is an important issue. In essence, what's at stake is the interpretation of +the absolute magnitudes of beta weights that are derived from a GLM analysis. +Typically, the amplitude of an evoked hemodynamic response is expressed in terms +of percent signal change (PSC) (e.g. by dividing the amplitude increase by some +baseline signal level and then multiplying by 100). Now, if all +voxels/areas/subjects shared exactly the same hemodynamic timecourse shape (i.e. +HRF) and we used this HRF in the analysis, then things would be easy. However, +there are real differences in timecourse shape across voxels/areas/subjects. +(And timecourses can actually also change across scan sessions for the same +voxel(s), due to physiology changes.) Thus, an approach that just assumes a +single fixed HRF has the advantage of being easy to think about, but already +will produce biases in PSC magnitudes. (For example, if you just change the +single fixed HRF used in an analysis, you will likely see that some voxel's +magnitudes go up and other voxels' magnitudes go down.) Now, consider the +approach of tailoring the HRF used in the analysis for individual voxels. +Certainly, the interpretation has to follow carefully -- the set of betas we +observe from a given voxel needs to be interpreted as the amplitude of responses +that appear to be present in the data _assuming_ the HRF selected for that +voxel. Note that in GLMsingle, the library of HRFs are normalized such that the +peak HRF value is 1 for each of the HRFs, so you can conveniently think of the +betas that result for a given HRF as literally the amplitude. Nonetheless, it +remains an open question how much we can actually interpret differences in BOLD +amplitudes across voxels/areas. For example, if one region seems to have higher +PSC than another region (even within the same subject), this doesn't necessarily +mean that the underlying local neural activity is actually higher in the first +region. Finally, we note that one approach to "get rid" of these types of +amplitude issues is to normalize (e.g. z-score) the responses you observe from a +voxel (or region) over the course of a scan session. Of course, one should think +carefully about the implications of such an approach for downstream analyses... + +### What is the interpretation of the scale and offset? + +By default, after the application of ridge regression, GLMsingle applies a +post-hoc scale and offset to the single-trial betas obtained for a given voxel +to best match what is obtained in the unregularized case. The reason for this is +that due to shrinkage, the betas for a voxel are generally "shrunken" down close +to 0, and therefore have some bias to be small in size. The theory is that we +can undo some of the bias by applying a simple scale and offset to the beta +weights. For many post-hoc uses of the beta estimates, a scale and offset is +probably not a big deal, but certainly one should think hard about whether or +not this matters to the analyses you are trying to do. One can always omit the +scale and offset (via appropriate setting of input options) and/or avoid ridge +regression altogether. + +### In the NSD data paper, in the calculation of the noise ceiling, the beta weights are z-scored. What's going on there? + +The NSD experiment involves aggregating responses across many distinct scan +sessions. Z-scoring is proposed as a potential (simple) method to reduce +instabilities that exist across different scan sessions. Obviously, z-scoring +can throw away relevant information, so one should be careful in that regard. + +## Things to watch out for + +If you use the input option to use a canonical HRF (instead of the library of +HRFs), the filename that is written is still "`FITHRF`" (even though the results +reflect what the user specified). From a820914da8c1cf569606d82c846dcec404bf3c20 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Thu, 28 Apr 2022 17:22:55 +0200 Subject: [PATCH 10/11] reformat one help section as example of minimalist change for a good render --- matlab/utilities/calccod.m | 47 ++++++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/matlab/utilities/calccod.m b/matlab/utilities/calccod.m index f696273..bfef701 100644 --- a/matlab/utilities/calccod.m +++ b/matlab/utilities/calccod.m @@ -1,29 +1,37 @@ function f = calccod(x,y,dim,wantgain,wantmeansub,wantsafe) - -% function f = calccod(x,y,dim,wantgain,wantmeansub,wantsafe) +% +% USAGE:: +% +% f = calccod(x,y,dim,wantgain,wantmeansub,wantsafe) % % , are matrices with the same dimensions +% % (optional) is the dimension of interest. -% default to 2 if is a (horizontal) vector and to 1 if not. -% special case is 0 which means to calculate globally. +% default to 2 if is a (horizontal) vector and to 1 if not. +% special case is 0 which means to calculate globally. +% % (optional) is -% 0 means normal -% 1 means allow a gain to be applied to each case of +% +% - 0 means normal +% - 1 means allow a gain to be applied to each case of % to minimize the squared error with respect to . % in this case, there cannot be any NaNs in or . -% 2 is like 1 except that gains are restricted to be non-negative. +% - 2 is like 1 except that gains are restricted to be non-negative. % so, if the gain that minimizes the squared error is negative, % we simply set the gain to be applied to be 0. -% default: 0. +% default: 0. +% % (optional) is -% 0 means do not subtract any mean. this makes it such that +% +% - 0 means do not subtract any mean. this makes it such that % the variance quantification is relative to 0. -% 1 means subtract the mean of each case of from both +% - 1 means subtract the mean of each case of from both % and before performing the calculation. this makes % it such that the variance quantification % is relative to the mean of each case of . % note that occurs before . -% default: 1. +% default: 1. +% % (optional) is whether to protect against NaNs. Default: 1. % % calculate the coefficient of determination (R^2) indicating @@ -45,18 +53,23 @@ % if there are no valid data points (i.e. all data points are % ignored because of NaNs), we return NaN for that case. % -% note some weird cases: -% calccod([],[]) is [] +% Note some weird cases: +% +% calccod([],[]) is [] +% +% Example:: +% +% x = randn(1,100); +% calccod(x,x+0.1*randn(1,100)) % + + % history: % 2013/08/18 - fix pernicious case where is all zeros and is 1 or 2. % 2010/11/28 - add ==2 case % 2010/11/23 - changed the output range to percentages. thus, the range is (-Inf,100]. % also, we removed the input since it was dumb. -% -% example: -% x = randn(1,100); -% calccod(x,x+0.1*randn(1,100)) + % input if ~exist('dim','var') || isempty(dim) From 5da4ecdde0933d2300ddd45daa3e27e7eb0a6c96 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Thu, 28 Apr 2022 17:40:04 +0200 Subject: [PATCH 11/11] do minimal reformating of help of GLMestimatesingletrial to silence warnings and ensure decent rendering --- matlab/GLMestimatesingletrial.m | 306 +++++++++++++++++++------------- 1 file changed, 180 insertions(+), 126 deletions(-) diff --git a/matlab/GLMestimatesingletrial.m b/matlab/GLMestimatesingletrial.m index 15ae55c..789aeca 100644 --- a/matlab/GLMestimatesingletrial.m +++ b/matlab/GLMestimatesingletrial.m @@ -1,123 +1,145 @@ function results = GLMestimatesingletrial(design,data,stimdur,tr,outputdir,opt) - -% function results = GLMestimatesingletrial(design,data,stimdur,tr,outputdir,opt) +% +% USAGE:: +% +% results = GLMestimatesingletrial(design,data,stimdur,tr,outputdir,opt) % % is the experimental design. There are two possible cases: +% % 1. A where A is a matrix with dimensions time x conditions. % Each column should be zeros except for ones indicating condition onsets. % 2. {A1 A2 A3 ... An} where each of the A's are like the previous case. % The different A's correspond to different runs, and different runs % can have different numbers of time points. However, all A's must have % the same number of conditions. -% Note that we ultimately compute single-trial response estimates (one estimate for -% each condition onset), and these will be provided in chronological order. However, -% by specifying that a given condition occurs more than one time over the course -% of the experiment, this information can and will be used for cross-validation purposes. +% +% Note that we ultimately compute single-trial response estimates (one estimate for +% each condition onset), and these will be provided in chronological order. However, +% by specifying that a given condition occurs more than one time over the course +% of the experiment, this information can and will be used for cross-validation purposes. +% % is the time-series data with dimensions X x Y x Z x time or a cell vector of -% elements that are each X x Y x Z x time. XYZ can be collapsed such that the data -% are given as a 2D matrix (units x time), which is useful for surface-format data. -% The dimensions of should mirror that of . For example, and -% should have the same number of runs, the same number of time points, etc. -% should not contain any NaNs. We automatically convert to single -% format if not already in single format. +% elements that are each X x Y x Z x time. XYZ can be collapsed such that the data +% are given as a 2D matrix (units x time), which is useful for surface-format data. +% The dimensions of should mirror that of . For example, and +% should have the same number of runs, the same number of time points, etc. +% should not contain any NaNs. We automatically convert to single +% format if not already in single format. +% % is the duration of a trial in seconds. For example, 3.5 means that you -% expect the neural activity from a given trial to last for 3.5 s. +% expect the neural activity from a given trial to last for 3.5 s. +% % is the sampling rate in seconds. For example, 1 means that we get a new -% time point every 1 s. Note that applies to both and . +% time point every 1 s. Note that applies to both and . +% % (optional) is a directory to which files will be written. (If the -% directory does not exist, we create it; if the directory already exists, -% we delete its contents so we can start fresh.) If you set -% to NaN, we will not create a directory and no files will be written. -% If you provide {outputdir figuredir}, we will save the large output files -% to and the small figure files to (either or both can be NaN). -% Default is 'GLMestimatesingletrialoutputs' (created in the current working directory). +% directory does not exist, we create it; if the directory already exists, +% we delete its contents so we can start fresh.) If you set +% to NaN, we will not create a directory and no files will be written. +% If you provide {outputdir figuredir}, we will save the large output files +% to and the small figure files to (either or both can be NaN). +% Default is 'GLMestimatesingletrialoutputs' (created in the current working directory). +% % (optional) is a struct with the following optional fields: % % *** MAJOR, HIGH-LEVEL FLAGS *** % -% (optional) is -% 0 means use an assumed HRF -% 1 means determine the best HRF for each voxel using the library-of-HRFs approach -% Default: 1. +% (optional) is: +% +% - 0 means use an assumed HRF +% - 1 means determine the best HRF for each voxel using the library-of-HRFs approach +% - Default: 1. +% % (optional) is -% 0 means do not perform GLMdenoise -% 1 means perform GLMdenoise -% Default: 1. +% +% - 0 means do not perform GLMdenoise +% - 1 means perform GLMdenoise +% - Default: 1. +% % (optional) is -% 0 means do not perform ridge regression -% 1 means perform ridge regression -% Default: 1. +% +% - 0 means do not perform ridge regression +% - 1 means perform ridge regression +% - Default: 1. +% % (optional) is the number of voxels that we will process at the same time. -% This number should be large in order to speed computation, but should not be so -% large that you run out of RAM. Default: 50000. +% This number should be large in order to speed computation, but should not be so +% large that you run out of RAM. Default: 50000. +% % (optional) is a cell vector of vectors of run indices, indicating the -% cross-validation scheme. For example, if we have 8 runs, we could use -% {[1 2] [3 4] [5 6] [7 8]} which indicates to do 4 folds of cross-validation, -% first holding out the 1st and 2nd runs, then the 3rd and 4th runs, etc. -% Default: {[1] [2] [3] ... [n]} where n is the number of runs. +% cross-validation scheme. For example, if we have 8 runs, we could use +% {[1 2] [3 4] [5 6] [7 8]} which indicates to do 4 folds of cross-validation, +% first holding out the 1st and 2nd runs, then the 3rd and 4th runs, etc. +% Default: {[1] [2] [3] ... [n]} where n is the number of runs. +% % (optional) is 1 x n (where n is the number of runs) with -% positive integers indicating the run groupings that are interpreted as -% "sessions". The purpose of this input is to allow for session-wise z-scoring -% of single-trial beta weights for the purposes of hyperparameter evaluation. -% For example, if you are analyzing data aggregated from multiple scan sessions, -% you may want beta weights to be z-scored per voxel within each session in order -% to compensate for any potential gross changes in betas across scan sessions. -% Note that the z-scoring has effect only INTERNALLY: it is used merely to -% calculate the cross-validation performance and the associated hyperparameter -% selection; the outputs of this function do not reflect z-scoring, and the user -% may wish to post-hoc apply z-scoring. Default: 1*ones(1,n) which means to interpret -% all runs as coming from the same session. +% positive integers indicating the run groupings that are interpreted as +% "sessions". The purpose of this input is to allow for session-wise z-scoring +% of single-trial beta weights for the purposes of hyperparameter evaluation. +% For example, if you are analyzing data aggregated from multiple scan sessions, +% you may want beta weights to be z-scored per voxel within each session in order +% to compensate for any potential gross changes in betas across scan sessions. +% Note that the z-scoring has effect only INTERNALLY: it is used merely to +% calculate the cross-validation performance and the associated hyperparameter +% selection; the outputs of this function do not reflect z-scoring, and the user +% may wish to post-hoc apply z-scoring. Default: 1*ones(1,n) which means to interpret +% all runs as coming from the same session. % % *** I/O FLAGS *** % % (optional) is a logical vector [A B C D] indicating which of the -% four model types to save to disk (assuming that they are computed). -% A = 0/1 for saving the results of the ONOFF model -% B = 0/1 for saving the results of the FITHRF model -% C = 0/1 for saving the results of the FITHRF_GLMDENOISE model -% D = 0/1 for saving the results of the FITHRF_GLMDENOISE_RR model -% Default: [1 1 1 1] which means save all computed results to disk. +% four model types to save to disk (assuming that they are computed). +% +% - A = 0/1 for saving the results of the ONOFF model +% - B = 0/1 for saving the results of the FITHRF model +% - C = 0/1 for saving the results of the FITHRF_GLMDENOISE model +% - D = 0/1 for saving the results of the FITHRF_GLMDENOISE_RR model +% - Default: [1 1 1 1] which means save all computed results to disk. +% % (optional) is a logical vector [A B C D] indicating which of the -% four model types to return in the output . The user must be careful with this, -% as large datasets can require a lot of RAM. If you do not request the various model types, -% they will be cleared from memory (but still potentially saved to disk). -% Default: [0 0 0 1] which means return only the final type-D model. +% four model types to return in the output . The user must be careful with this, +% as large datasets can require a lot of RAM. If you do not request the various model types, +% they will be cleared from memory (but still potentially saved to disk). +% Default: [0 0 0 1] which means return only the final type-D model. % % *** GLM FLAGS *** % % (optional) is time x regressors or a cell vector -% of elements that are each time x regressors. The dimensions of -% should mirror that of (i.e. same number of -% runs, same number of time points). The number of extra regressors -% does not have to be the same across runs, and each run can have zero -% or more extra regressors. If [] or not supplied, we do -% not use extra regressors in the model. +% of elements that are each time x regressors. The dimensions of +% should mirror that of (i.e. same number of +% runs, same number of time points). The number of extra regressors +% does not have to be the same across runs, and each run can have zero +% or more extra regressors. If [] or not supplied, we do +% not use extra regressors in the model. +% % (optional) is a non-negative integer with the maximum -% polynomial degree to use for polynomial nuisance functions, which -% are used to capture low-frequency noise fluctuations in each run. -% Can be a vector with length equal to the number of runs (this -% allows you to specify different degrees for different runs). -% Default is to use round(L/2) for each run where L is the -% duration in minutes of a given run. +% polynomial degree to use for polynomial nuisance functions, which +% are used to capture low-frequency noise fluctuations in each run. +% Can be a vector with length equal to the number of runs (this +% allows you to specify different degrees for different runs). +% Default is to use round(L/2) for each run where L is the +% duration in minutes of a given run. +% % (optional) is whether to convert amplitude estimates -% to percent BOLD change. This is done as the very last step, and is -% accomplished by dividing by the absolute value of 'meanvol' and -% multiplying by 100. (The absolute value prevents negative values in -% 'meanvol' from flipping the sign.) Default: 1. +% to percent BOLD change. This is done as the very last step, and is +% accomplished by dividing by the absolute value of 'meanvol' and +% multiplying by 100. (The absolute value prevents negative values in +% 'meanvol' from flipping the sign.) Default: 1. % % *** HRF FLAGS *** % % (optional) is time x 1 with an assumed HRF that characterizes the evoked -% response to each trial. We automatically divide by the maximum value so that the -% peak is equal to 1. Default is to generate a canonical HRF (see getcanonicalhrf.m). -% Note that the HRF supplied in is used in only two instances: -% (1) it is used for the simple ONOFF type-A model, and (2) if the user sets -% to 0, it is also used for the type-B, type-C, and type-D models. +% response to each trial. We automatically divide by the maximum value so that the +% peak is equal to 1. Default is to generate a canonical HRF (see getcanonicalhrf.m). +% Note that the HRF supplied in is used in only two instances: +% (1) it is used for the simple ONOFF type-A model, and (2) if the user sets +% to 0, it is also used for the type-B, type-C, and type-D models. +% % (optional) is time x H with H different HRFs to choose from for the -% library-of-HRFs approach. We automatically normalize each HRF to peak at 1. -% Default is to generate a library of 20 HRFs (see getcanonicalhrflibrary.m). -% Note that if is 0, is clobbered with the contents -% of , which in effect causes a single assumed HRF to be used. +% library-of-HRFs approach. We automatically normalize each HRF to peak at 1. +% Default is to generate a library of 20 HRFs (see getcanonicalhrflibrary.m). +% Note that if is 0, is clobbered with the contents +% of , which in effect causes a single assumed HRF to be used. % % *** MODEL TYPE A (ONOFF) FLAGS *** % @@ -126,58 +148,66 @@ % *** MODEL TYPE B (FITHRF) FLAGS *** % % (optional) is 0/1 indicating whether "least-squares-separate" estimates -% are desired. If 1, then the type-B model will be estimated using the least-squares- -% separate method (as opposed to ordinary least squares). Default: 0. +% are desired. If 1, then the type-B model will be estimated using the least-squares- +% separate method (as opposed to ordinary least squares). Default: 0. % % *** MODEL TYPE C (FITHRF_GLMDENOISE) FLAGS *** % % (optional) is a non-negative integer indicating the maximum -% number of PCs to enter into the model. Default: 10. +% number of PCs to enter into the model. Default: 10. +% % (optional) is [A B] where A is a percentile for voxel intensity -% values and B is a fraction to apply to the percentile. These parameters -% are used in the selection of the noise pool. Default: [99 0.1]. +% values and B is a fraction to apply to the percentile. These parameters +% are used in the selection of the noise pool. Default: [99 0.1]. +% % (optional) is an R^2 value (percentage). After fitting the type-A model, -% voxels whose R^2 is below this value are allowed to enter the noise pool. -% Default is [] which means to automatically determine a good value. +% voxels whose R^2 is below this value are allowed to enter the noise pool. +% Default is [] which means to automatically determine a good value. +% % (optional) is X x Y x Z (or XYZ x 1) with 1s indicating voxels to -% specifically exclude when selecting the noise pool. 0 means all voxels can be -% potentially chosen. Default: 0. +% specifically exclude when selecting the noise pool. 0 means all voxels can be +% potentially chosen. Default: 0. +% % (optional) is an R^2 value (percentage). To decide the number -% of PCs to include, we examine a subset of the available voxels. Specifically, -% we examine voxels whose type-A model R^2 is above . Default is [] -% which means to automatically determine a good value. +% of PCs to include, we examine a subset of the available voxels. Specifically, +% we examine voxels whose type-A model R^2 is above . Default is [] +% which means to automatically determine a good value. +% % (optional) is X x Y x Z (or XYZ x 1) with 1s indicating all possible -% voxels to consider when selecting the subset of voxels. 1 means all voxels -% can be potentially selected. Default: 1. +% voxels to consider when selecting the subset of voxels. 1 means all voxels +% can be potentially selected. Default: 1. +% % (optional) is -% A: a number greater than or equal to 1 indicating when to stop adding PCs -% into the model. For example, 1.05 means that if the cross-validation -% performance with the current number of PCs is within 5% of the maximum -% observed, then use that number of PCs. (Performance is measured -% relative to the case of 0 PCs.) When is 1, the selection -% strategy reduces to simply choosing the PC number that achieves -% the maximum. The advantage of stopping early is to achieve a selection -% strategy that is robust to noise and shallow performance curves and -% that avoids overfitting. -% -B: where B is the number of PCs to use for the final model. B can be any -% integer between 0 and opt.numpcstotry. Note that if -B case is used, -% cross-validation is NOT performed for the type-C model, and instead we -% blindly use B PCs. -% Default: 1.05. +% +% - A: a number greater than or equal to 1 indicating when to stop adding PCs +% into the model. For example, 1.05 means that if the cross-validation +% performance with the current number of PCs is within 5% of the maximum +% observed, then use that number of PCs. (Performance is measured +% relative to the case of 0 PCs.) When is 1, the selection +% strategy reduces to simply choosing the PC number that achieves +% the maximum. The advantage of stopping early is to achieve a selection +% strategy that is robust to noise and shallow performance curves and +% that avoids overfitting. +% - B: where B is the number of PCs to use for the final model. B can be any +% integer between 0 and opt.numpcstotry. Note that if -B case is used, +% cross-validation is NOT performed for the type-C model, and instead we +% blindly use B PCs. +% - Default: 1.05. % % *** MODEL TYPE D (FITHRF_GLMDENOISE_RR) FLAGS *** % % (optional) is a vector of fractions that are greater than 0 -% and less than or equal to 1. We automatically sort in descending order and -% ensure the fractions are unique. These fractions indicate the regularization -% levels to evaluate using fractional ridge regression (fracridge) and -% cross-validation. Default: fliplr(.05:.05:1). A special case is when -% is specified as a single scalar value. In this case, cross-validation -% is NOT performed for the type-D model, and we instead blindly use -% the supplied fractional value for the type-D model. +% and less than or equal to 1. We automatically sort in descending order and +% ensure the fractions are unique. These fractions indicate the regularization +% levels to evaluate using fractional ridge regression (fracridge) and +% cross-validation. Default: fliplr(.05:.05:1). A special case is when +% is specified as a single scalar value. In this case, cross-validation +% is NOT performed for the type-D model, and we instead blindly use +% the supplied fractional value for the type-D model. +% % (optional) is whether to automatically scale and offset -% the model estimates from the type-D model to best match the unregularized -% estimates. Default: 1. +% the model estimates from the type-D model to best match the unregularized +% estimates. Default: 1. % % This function computes up to four model outputs (called type-A (ONOFF), % type-B (FITHRF), type-C (FITHRF_GLMDENOISE), and type-D (FITHRF_GLMDENOISE_RR)), @@ -187,24 +217,30 @@ % There are a variety of cases that you can achieve. Here are some examples: % % - wantlibrary=1, wantglmdenoise=1, wantfracridge=1 [Default] -% A = simple ONOFF model -% B = single-trial estimates using a tailored HRF for every voxel -% C = like B but with GLMdenoise regressors added into the model -% D = like C but with ridge regression regularization (tailored to each voxel) +% +% - A = simple ONOFF model +% - B = single-trial estimates using a tailored HRF for every voxel +% - C = like B but with GLMdenoise regressors added into the model +% - D = like C but with ridge regression regularization (tailored to each voxel) % % - wantlibrary=0 +% % A fixed assumed HRF is used in all model types. % % - wantglmdenoise=0, wantfracridge=0 +% % Model types C and D are not computed. % % - wantglmdenoise=0, wantfracridge=1 +% % Model type C is not computed; model type D is computed using 0 GLMdenoise regressors. % % - wantglmdenoise=1, wantfracridge=0 +% % Model type C is computed; model type D is not computed. % % - wantlss=1 +% % Model type B is computed, but using least-squares-separate instead of OLS. % Other model types, if computed, use OLS. % @@ -218,30 +254,48 @@ % a single scalar opt.fracs. In other words, you can perform wantfracridge % without any cross-validation, but you need to provide opt.fracs as a scalar. % -% OUTPUTS: +% *** OUTPUTS: *** % % There are various outputs for each of the four model types: % % is either -% (1) the HRF (time x 1) and ON-OFF beta weights (X x Y x Z) -% (2) the full set of single-trial beta weights (X x Y x Z x TRIALS) +% +% 1. the HRF (time x 1) and ON-OFF beta weights (X x Y x Z) +% 2. the full set of single-trial beta weights (X x Y x Z x TRIALS) +% % is model accuracy expressed in terms of R^2 (percentage). +% % is R2 separated by run +% % is the mean of all volumes +% % is the R2 for each of the different HRFs in the library +% % is separated by run +% % is the 1-index of the best HRF +% % is HRFindex separated by run +% % indicates voxels selected for the noise pool +% % indicates the full set of candidate GLMdenoise regressors that were found +% % is the cross-validation results for GLMdenoise +% % is the set of voxels used to summarize GLMdenoise cross-validation results +% % is the summary GLMdenoise cross-validation result on which pcnum selection is done +% % is the number of PCs that were selected for the final model +% % is the fractional ridge regression regularization level chosen for each voxel +% % is the cross-validation results for the ridge regression +% % is the scale and offset applied to RR estimates to best match the unregularized result % + % History (keep in line with Python version): % - 2021/06/03 - add the option for use to specify two directories for % - 2020/08/22 - Implement opt.sessionindicator. Also, now the cross-validation units now reflect