Repository: MICA-MNI/BrainSpace (language: python)

Issue:
Subject list for group connectivity matrices
I was wondering if the list of HCP subjects used for the group functional connectivity matrices with `load_group_fc`/`load_group_mpc` was available. My understanding/assumption is that these subsets of HCP data are based on the selection procedures described in [Vos de Wael et al (2018)](https://www.pnas.org/content/115/40/10154), but I'm not 100% certain in that. 

I'm looking to compute similar connectivity matrices on the samples used in Vos de Wael et al (2018), but I also think it would be great if the documentation could clarify the source of the group connectivity matrices. 
diffusion time default value discrepancy 
Hi,

The following may be on purpose (I find it a bit inconsistent though).

In gradient/embedding.py , the class DiffusionMaps has a default diffusion time value set to 1, I believe it should be 0.
The documentation of the class says the default is 0, also the function doing the work "diffusion_mapping" has default at 0.

A GradientMaps would inherit this time diffusion of 1, when doing something like in the documentation ie:

```python
gm = GradientMaps(n_components=2, approach='dm', kernel='normalized_angle')
gm.fit(m) #this would have diffusion time of 1
```

--
Nicolas Ninin
[GIT] Changed Travis CI for GitHub actions
List of changes:
- Removed Travis CI
- Implemented Github actions to run:
  - pytest on all tests using windows-latest, macos-latest, and ubuntu-latest on Python 3.6, 3.7, and 3.8. 
  - pytest with code coverage on macos-latest. Code coverage runs only on this branch and the master branch.
 - Added code coverage and github actions badges to README.rst and ReadTheDocs main page.
 
 To do:
 - Fix (out of memory?) errors in ubuntu and windows. 
diffusion time default value discrepancy 
Hi,

The following may be on purpose (I find it a bit inconsistent though).

In gradient/embedding.py , the class DiffusionMaps has a default diffusion time value set to 1, I believe it should be 0.
The documentation of the class says the default is 0, also the function doing the work "diffusion_mapping" has default at 0.

A GradientMaps would inherit this time diffusion of 1, when doing something like in the documentation ie:

```python
gm = GradientMaps(n_components=2, approach='dm', kernel='normalized_angle')
gm.fit(m) #this would have diffusion time of 1
```

--
Nicolas Ninin


Discussion / hints:

Fix patch:
diff --git a/.DS_Store b/.DS_Store
deleted file mode 100644
index d4c2034..0000000
Binary files a/.DS_Store and /dev/null differ
diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
new file mode 100644
index 0000000..1495e17
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -0,0 +1,24 @@
+---
+name: Bug report
+about: Create a report to help us improve
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+* **What is the current behavior?**
+
+
+* **Please provide the steps to reproduce and if possible a minimal demo of the problem**.
+
+
+* **What is the expected behavior?**
+
+
+
+* **Please tell us about your computing environment:**
+ 
+
+
+* **Other information**
diff --git a/.github/ISSUE_TEMPLATE/development_todo_list.md b/.github/ISSUE_TEMPLATE/development_todo_list.md
new file mode 100644
index 0000000..1f9c3d3
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/development_todo_list.md
@@ -0,0 +1,14 @@
+---
+name: Development to-do list.
+about: Set lists of things to-do.
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+A short summary of what should be done. 
+
+- [ ] List item 1.
+- [ ] List item 2.
+- [ ] List item 3.
diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md
new file mode 100644
index 0000000..37c65e7
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/feature_request.md
@@ -0,0 +1,19 @@
+---
+name: Feature request
+about: Suggest an idea for BrainSpace
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+* **What is the current behavior?**
+
+
+* **What is the motivation / use case for changing the behavior?**
+
+
+* **If relevant, please tell us about your computing environment:** 
+  
+
+* **Other information**
diff --git a/.github/ISSUE_TEMPLATE/support-request.md b/.github/ISSUE_TEMPLATE/support-request.md
new file mode 100644
index 0000000..b47e252
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/support-request.md
@@ -0,0 +1,16 @@
+---
+name: Support request
+about: Ask for help running BrainSpace
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+* **What are you trying to do?**
+
+* **If you've tried running this with BrainSpace code already, how did you do so and what went wrong?**
+
+* **Please tell us about your computing environment:**
+  
+* **Other information**
diff --git a/.github/workflows/github_actions_ci.yml b/.github/workflows/github_actions_ci.yml
new file mode 100644
index 0000000..68a5811
--- /dev/null
+++ b/.github/workflows/github_actions_ci.yml
@@ -0,0 +1,57 @@
+name: Python CI
+
+on:
+  push:
+    branches:
+      - master
+  pull_request:
+
+jobs:
+  build:
+
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, windows-latest, macos-latest]
+        python-version: [3.6, 3.7, 3.8, 3.9]
+
+    steps:
+    - uses: actions/checkout@v2
+    
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: ${{ matrix.python-version }}
+    
+    - name: Test Windows
+      if: matrix.os == 'windows-latest'
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install pytest
+        python -m pip install -r requirements.txt
+        python -m pytest brainspace --ignore brainspace/tests/test_plotting.py
+
+
+    - name: Prepare Linux
+      if: matrix.os == 'ubuntu-latest'
+      run: |
+        python3 -m pip install --upgrade pip
+        python3 -m pip install pytest
+        python3 -m pip install -r requirements.txt
+          
+    - name: Test Linux
+      if: matrix.os == 'ubuntu-latest'
+      uses: GabrielBB/xvfb-action@v1.5
+      with:
+        options: -screen 0 1600x1200x24
+        run: python3 -m pytest brainspace
+
+
+    - name: Test macOS 
+      if: matrix.os == 'macos-latest'
+      run: |
+        python3 -m pip install --upgrade pip
+        python3 -m pip install pytest
+        python3 -m pip install -r requirements.txt
+        python3 -m pytest brainspace
diff --git a/.github/workflows/python_coverage.yml b/.github/workflows/python_coverage.yml
new file mode 100644
index 0000000..1488cea
--- /dev/null
+++ b/.github/workflows/python_coverage.yml
@@ -0,0 +1,29 @@
+name: Python Code Coverage
+on:
+  push:
+    branches:
+      - 'master'
+
+jobs:
+  codecov:
+    name: Codecov Workflow
+    runs-on: macos-latest
+    
+    steps:
+      - uses: actions/checkout@v1
+      - name: Set up Python
+        uses: actions/setup-python@master
+        with:
+          python-version: 3.8
+      - name: Generate coverage report
+        run: |
+          python3 -m pip install --upgrade pip
+          python3 -m pip install pytest pytest-cov
+          python3 -m pip install -r requirements.txt
+          python3 -m pytest --cov=./ --cov-report=xml 
+      - name: Upload coverage to Codecov
+        uses: codecov/codecov-action@v1
+        with:
+          token: ${{ secrets.CODECOV_TOKEN }}
+          file: ./coverage.xml
+          flags: unittests
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index 1ad885e..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,111 +0,0 @@
-language: python
-sudo: false
-dist: xenial
-
-
-matrix:
-  include:
-
-    - name: "Python 3.5.6 on Xenial Linux"
-      python: 3.5
-      services: xvfb
-
-    - name: "Python 3.6.7 on Xenial Linux"
-      python: 3.6
-      services: xvfb
-
-    - name: "Python 3.7.1 on Xenial Linux"
-      python: 3.7
-      services: xvfb
-      env: REPORT_COVERAGE=true
-
-    - name: "Python 3.7.1 setup.py install"
-      python: 3.7
-      services: xvfb
-      install:
-        - python setup.py install
-
-    - name: "Python 3.7.1 sdist"
-      python: 3.7
-      services: xvfb
-      install:
-        - python setup.py sdist;
-        - pip install dist/*.tar.gz;
-
-    - name: "Python 3.7.1 wheel"
-      python: 3.7
-      services: xvfb
-      install:
-        - python setup.py bdist_wheel;
-        - pip install dist/*.whl;
-
-    - name: "Python 3.7.1 conda"
-      python: 3.7
-      services: xvfb
-      before_install: []
-      install:
-        # Install Anaconda
-        - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then
-            wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh -O miniconda.sh;
-          else
-            wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
-          fi
-        - bash miniconda.sh -b -p $HOME/miniconda
-        - export PATH="$HOME/miniconda/bin:$PATH"
-        - hash -r
-        - conda config --set always_yes yes --set changeps1 no
-        - conda update -q conda
-        # Useful for debugging any issues with conda
-        - conda info -a
-
-        - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy>=1.11.0 scipy>=0.17.0 scikit-learn>=0.20.0 matplotlib>=2.0.0 vtk>=8.1.0 pillow pytest
-        - source activate test-environment
-        - pip install nibabel
-        - pip install --no-deps .
-
-    # ======= OSX ========
-    - name: "Python 3.6.5 on macOS 10.13"
-      os: osx
-      osx_image: xcode9.4
-      language: generic
-
-    - name: "Python 3.7.4 on macOS 10.14"
-      os: osx
-      osx_image: xcode11
-      language: generic
-
-
-before_install:
-  - if [[ "${TRAVIS_OS_NAME}" == "osx" ]]; then
-      python3 --version;
-      pip3 install -U pip;
-    else
-      python --version;
-      pip install -U pip;
-    fi
-
-
-# command to install dependencies
-install:
-  - if [[ "${TRAVIS_OS_NAME}" == "osx" ]]; then
-      pip3 install -e ".[test]";
-    else
-      pip install -e ".[test]";
-    fi
-
-
-# command to run tests
-script:
-  - if [[ "${TRAVIS_OS_NAME}" == "osx" ]]; then
-      python3 -m pytest;
-    else
-      py.test brainspace;
-    fi
-
-
-after_success:
-  # Upload coverage data to codecov
-  - if [[ "$REPORT_COVERAGE" == "true" ]]; then
-      echo "Uploading coverage to Codecov";
-      codecov;
-    fi
diff --git a/BrainSpace.mltbx b/BrainSpace.mltbx
new file mode 100644
index 0000000..c6b0379
Binary files /dev/null and b/BrainSpace.mltbx differ
diff --git a/BrainSpace.prj b/BrainSpace.prj
new file mode 100644
index 0000000..c681963
--- /dev/null
+++ b/BrainSpace.prj
@@ -0,0 +1,162 @@
+<deployment-project plugin="plugin.toolbox" plugin-version="1.0">
+  <configuration file="/Users/reinder/GitHub/BrainSpace/BrainSpace.prj" location="/Users/reinder/GitHub/BrainSpace" name="BrainSpace" target="target.toolbox" target-name="Package Toolbox">
+    <param.appname>BrainSpace</param.appname>
+    <param.authnamewatermark>MICA Lab</param.authnamewatermark>
+    <param.email>reinder.vosdewael@gmail.com</param.email>
+    <param.company>McGill University</param.company>
+    <param.summary>BrainSpace is a lightweight cross-platform toolbox primarily intended for macroscale gradient mapping and analysis.</param.summary>
+    <param.description>BrainSpace is a lightweight cross-platform toolbox primarily intended for macroscale gradient mapping and analysis of neuroimaging and connectome level data. The current version of BrainSpace is available in Python and MATLAB, programming languages widely used by the neuroimaging and network neuroscience communities. The toolbox also contains several maps that allow for exploratory analysis of gradient correspondence with other brain-derived features, together with tools to generate spatial null models.
+
+If you consider using BrainSpace, please cite our manuscript: Vos de Wael R, Benkarim O, Paquola C, Lariviere S, Royer J, Tavakol S, Xu T, Hong S, Langs G, Valk S, Misic B, Milham M, Margulies D, Smallwood J, Bernhardt B (2020). BrainSpace: a toolbox for the analysis of macroscale gradients in neuroimaging and connectomics datasets. Commun Biol 3, 103.</param.description>
+    <param.screenshot>/Users/reinder/Desktop/1612370158.png</param.screenshot>
+    <param.version>0.1.2</param.version>
+    <param.output>${PROJECT_ROOT}/BrainSpace.mltbx</param.output>
+    <param.products.name>
+      <item>MATLAB</item>
+    </param.products.name>
+    <param.products.id>
+      <item>1</item>
+    </param.products.id>
+    <param.products.version>
+      <item>9.10</item>
+    </param.products.version>
+    <param.platforms />
+    <param.guid>cc685c6b-b946-4d79-b3d0-762743da5f8e</param.guid>
+    <param.exclude.filters />
+    <param.exclude.pcodedmfiles>true</param.exclude.pcodedmfiles>
+    <param.examples />
+    <param.demosxml />
+    <param.apps />
+    <param.registered.apps />
+    <param.docs />
+    <param.getting.started.guide />
+    <param.matlabpath.excludes />
+    <param.javaclasspath.excludes />
+    <param.exported.on.package>false</param.exported.on.package>
+    <param.required.addons />
+    <param.matlab.project.id />
+    <param.matlab.project.name />
+    <param.release.start>R2018b</param.release.start>
+    <param.release.end>latest</param.release.end>
+    <param.release.current.only>false</param.release.current.only>
+    <param.compatiblity.windows>true</param.compatiblity.windows>
+    <param.compatiblity.macos>true</param.compatiblity.macos>
+    <param.compatiblity.linux>true</param.compatiblity.linux>
+    <param.compatiblity.matlabonline>true</param.compatiblity.matlabonline>
+    <param.installation.map />
+    <param.additional.sw.names />
+    <param.additional.sw.licenses />
+    <param.additional.sw.win.url />
+    <param.additional.sw.mac.url />
+    <param.additional.sw.linux.url />
+    <unset>
+      <param.output />
+      <param.platforms />
+      <param.exclude.filters />
+      <param.exclude.pcodedmfiles />
+      <param.examples />
+      <param.demosxml />
+      <param.apps />
+      <param.registered.apps />
+      <param.docs />
+      <param.getting.started.guide />
+      <param.matlabpath.excludes />
+      <param.javaclasspath.excludes />
+      <param.exported.on.package />
+      <param.required.addons />
+      <param.matlab.project.id />
+      <param.matlab.project.name />
+      <param.release.current.only />
+      <param.compatiblity.windows />
+      <param.compatiblity.macos />
+      <param.compatiblity.linux />
+      <param.compatiblity.matlabonline />
+      <param.installation.map />
+      <param.additional.sw.names />
+      <param.additional.sw.licenses />
+      <param.additional.sw.win.url />
+      <param.additional.sw.mac.url />
+      <param.additional.sw.linux.url />
+    </unset>
+    <fileset.rootdir>
+      <file>${PROJECT_ROOT}/matlab</file>
+    </fileset.rootdir>
+    <fileset.rootfiles>
+      <file>${PROJECT_ROOT}/matlab/@GradientMaps</file>
+      <file>${PROJECT_ROOT}/matlab/analysis_code</file>
+      <file>${PROJECT_ROOT}/matlab/example_data_loaders</file>
+      <file>${PROJECT_ROOT}/matlab/plot_data</file>
+      <file>${PROJECT_ROOT}/matlab/surface_manipulation</file>
+    </fileset.rootfiles>
+    <fileset.depfun.included />
+    <fileset.depfun.excluded>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/Contents.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/display.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/export.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/fieldnames.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/gifti.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/isfield.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/plot.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/base64decode.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/base64encode.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/getdict.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/isintent.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/mvtk_write.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/read_freesurfer_file.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/read_gifti_file_standalone.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/xml_parser.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/private/zstream.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/save.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/saveas.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/struct.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/subsasgn.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@gifti/subsref.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/@variogram/variogram.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/Contents.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/display.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/export.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/fieldnames.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/gifti.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/isfield.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/plot.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/base64decode.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/base64encode.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/getdict.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/isintent.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/mvtk_write.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/read_freesurfer_file.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/read_gifti_file_standalone.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/xml_parser.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/private/zstream.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/save.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/saveas.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/struct.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/subsasgn.m</file>
+      <file>/Users/reinder/OneDrive/Academics/MATLAB/packages/cifti-matlab/@gifti/subsref.m</file>
+    </fileset.depfun.excluded>
+    <fileset.package />
+    <build-deliverables>
+      <file location="${PROJECT_ROOT}" name="BrainSpace.mltbx" optional="false">/Users/reinder/GitHub/BrainSpace/BrainSpace.mltbx</file>
+    </build-deliverables>
+    <workflow />
+    <matlab>
+      <root>/Applications/MATLAB_R2021a.app</root>
+      <toolboxes />
+    </matlab>
+    <platform>
+      <unix>true</unix>
+      <mac>true</mac>
+      <windows>false</windows>
+      <win2k>false</win2k>
+      <winxp>false</winxp>
+      <vista>false</vista>
+      <linux>false</linux>
+      <solaris>false</solaris>
+      <osver>10.14.6</osver>
+      <os32>false</os32>
+      <os64>true</os64>
+      <arch>maci64</arch>
+      <matlab>true</matlab>
+    </platform>
+  </configuration>
+</deployment-project>
\ No newline at end of file
diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index 0ab606f..40b3d27 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,3 +1,32 @@
+v0.1.2, XXXXX 2021 -- Let there be null models update
+MATLAB
+- New functionality
+    - Added a verbose flag to GradientMaps (default value: false).
+    - plot_hemispheres has several new functionalities:
+        - Now returns an object rather than the graphics handles.
+          Old graphics handles are stored in object.handles. 
+        - Added the following methods to the plot_hemispheres object for 
+        modifying figure properties:
+            - object.colorlimits: modifies color limits.
+            - object.colormaps: modifies colormaps. 
+            - object.labels: modifies labeltext properties. 
+            - Added a "views" name-value pair that enables anterior, 
+              posterior, superior, and inferior views. 
+    - Added a MATLAB implementation of BrainSMASH (i.e. variogram 
+    matching null model)
+    - procrustes_alignment now outputs the rotation matrices as a second
+    output variable. 
+
+Python:
+  - New functionality
+    - brainspace.mesh.mesh_io.read_surface now supports gzipped files. 
+
+ReadTheDocs
+ - Updated MATLAB documentation of GradientMaps with the verbose flag.
+ - Updated MATLAB documentation of plot_hemispheres to reflect the new 
+   changes.  
+ - Added MATLAB documentation for the variogram function. 
+
 Hotfix, January 2021 -- MATLAB Joint Alignment
 MATLAB:
 - Fixed faulty dimensionality extraction in joint alignment of asymmetric
diff --git a/README.rst b/README.rst
index 5f59ee0..5b995a1 100644
--- a/README.rst
+++ b/README.rst
@@ -3,8 +3,11 @@ BrainSpace
 ====================
 
 
-.. image:: https://img.shields.io/travis/MICA-MNI/BrainSpace/master?style=flat&logo=travis
-   :target: https://travis-ci.org/MICA-MNI/BrainSpace
+.. image:: https://github.com/MICA-LAB/brainspace/workflows/Python%20package/badge.svg
+   :target: https://github.com/MICA-LAB/brainspace/actions
+
+.. image:: https://codecov.io/gh/mica-mni/brainspace/branch/master/graph/badge.svg
+   :target: https://codecov.io/gh/mica-mni/brainspace
 
 .. image:: https://img.shields.io/appveyor/build/OualidBenkarim/brainspace/master?logo=appveyor
    :alt: AppVeyor branch
diff --git a/brainspace/.DS_Store b/brainspace/.DS_Store
deleted file mode 100644
index d98352b..0000000
Binary files a/brainspace/.DS_Store and /dev/null differ
diff --git a/brainspace/_version.py b/brainspace/_version.py
index 68cfb48..2462202 100644
--- a/brainspace/_version.py
+++ b/brainspace/_version.py
@@ -1,3 +1,3 @@
 """BrainSpace version"""
 
-__version__ = '0.1.1'
+__version__ = '0.1.2'
diff --git a/brainspace/datasets/base.py b/brainspace/datasets/base.py
index 4b88d9d..bf73d5c 100644
--- a/brainspace/datasets/base.py
+++ b/brainspace/datasets/base.py
@@ -12,7 +12,8 @@
 def load_group_fc(parcellation, scale=400, group='main'):
     """ Load group level connectivity matrix for a given parcellation.
 
-    Connectivity is derived from a subset of HCP data.
+    Connectivity is derived from a the discovery (main) and validation (holdout) 
+    datasets of (Vos de Wael et al., 2018, PNAS).
 
     Parameters
     ----------
@@ -40,7 +41,7 @@ def load_group_fc(parcellation, scale=400, group='main'):
 def load_group_mpc(parcellation, scale=400):
     """ Load group level connectivity matrix for a given parcellation.
 
-    Connectivity is derived from a subset of HCP data.
+    MPC is derived from the dataset used by (Paquola et al., 2019, PLoS biology).
 
     Parameters
     ----------
@@ -266,6 +267,9 @@ def _load_feat(feat_name, parcellation=None, mask=None):
 def load_marker(name, join=False):
     """ Load cortical data for conte69.
 
+    Markers are derived from the discovery (main) datasets of (Vos de Wael et
+    al., 2018, PNAS).
+
     Parameters
     ----------
     name : {'curvature', 'thickness', 't1wt2w'}
diff --git a/brainspace/examples/plot_tutorial3.py b/brainspace/examples/plot_tutorial3.py
index 04df1f6..f09e8c9 100644
--- a/brainspace/examples/plot_tutorial3.py
+++ b/brainspace/examples/plot_tutorial3.py
@@ -5,9 +5,10 @@
 canonical gradient and data from other modalities (curvature, cortical
 thickness and T1w/T2w image intensity). A normal test of the significance of
 the correlation cannot be used, because the spatial auto-correlation in MRI
-data may bias the test statistic. In this tutorial we will show two approaches
-for null hypothesis testing: spin permutations and Moran spectral
-randomization.
+data may bias the test statistic. In this tutorial we will show three
+approaches for null hypothesis testing: spin permutations, Moran spectral
+randomization, and autocorrelation-preserving surrogates based on variogram
+matching.
 
 .. note::
     When using either approach to compare gradients to non-gradient markers,
@@ -117,14 +118,13 @@
     pv_spin = np.mean(np.abs(r_spin) >= np.abs(r_obs))
 
     # Plot null dist
-    axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))
+    axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
     axs[k].axvline(r_obs, lw=2, ls='--', color='k')
-    axs[k].set_xlabel('Correlation with {}'.format(fn))
+    axs[k].set_xlabel(f'Correlation with {fn}')
     if k == 0:
         axs[k].set_ylabel('Density')
 
-    print('{}:\n Obs : {:.5e}\n Spin: {:.5e}\n'.
-          format(fn.capitalize(), pv_obs, pv_spin))
+    print(f'{fn.capitalize()}:\n Obs : {pv_obs:.5e}\n Spin: {pv_spin:.5e}\n')
 
 fig.tight_layout()
 plt.show()
@@ -208,22 +208,97 @@
     pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))
 
     # Plot null dist
-    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))
+    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
     axs[k].axvline(r_obs, lw=2, ls='--', color='k')
-    axs[k].set_xlabel('Correlation with {}'.format(fn))
+    axs[k].set_xlabel(f'Correlation with {fn}')
     if k == 0:
         axs[k].set_ylabel('Density')
 
-    print('{}:\n Obs  : {:.5e}\n Moran: {:.5e}\n'.
-          format(fn.capitalize(), pv_obs, pv_rand))
+    print(f'{fn.capitalize()}:\n Obs  : {pv_obs:.5e}\n Moran: {pv_rand:.5e}\n')
 
 fig.tight_layout()
 plt.show()
 
 
 ###############################################################################
-# There are some scenarios where MSR results do not follow a normal
-# distribution. It is relatively simple to check whether this occurs in our
-# data by visualizing the null distributions. Check this interesting paper
-# for more information `(Burt et al., 2020) <https://www.biorxiv.org/content/
-# 10.1101/2020.02.18.955054v1>`_.
+#
+#
+#
+#
+# Variogram Matching
+# ------------------
+#
+# Here, we will repeat the same analysis using the variogram matching approach
+# presented in `(Burt et al., 2020)
+# <https://www.sciencedirect.com/science/article/pii/S1053811920305243>`_,
+# which generates novel brainmaps with similar spatial autocorrelation to the
+# input data.
+#
+#
+# We will need a distance matrix that tells us what the spatial distance
+# between our datapoints is. For this example, we will use geodesic distance.
+
+
+from brainspace.mesh.mesh_elements import get_immediate_distance
+from scipy.sparse.csgraph import dijkstra
+
+# Compute geodesic distance
+gd = get_immediate_distance(surf_lh, mask=mask_tl)
+gd = dijkstra(gd, directed=False)
+
+idx_sorted = np.argsort(gd, axis=1)
+
+
+###############################################################################
+# Now we've got everything we need to generate our surrogate datasets. By
+# default, BrainSpace will use all available data to generate surrogate maps.
+# However, this process is extremely computationally and memory intensive. When
+# using this method with more than a few hundred regions, we recommend
+# subsampling the data. This can be done using SampledSurrogateMaps instead of
+# the SurrogateMaps.
+
+from brainspace.null_models import SampledSurrogateMaps
+
+
+n_surrogate_datasets = 1000
+
+# Note: number samples must be greater than number neighbors
+num_samples = 100
+num_neighbors = 50
+
+ssm = SampledSurrogateMaps(ns=num_samples, knn=num_neighbors, random_state=0)
+ssm.fit(gd, idx_sorted)
+
+t1wt2w_surrogates = ssm.randomize(t1wt2w_tl, n_rep=n_surrogate_datasets)
+curv_surrogates = ssm.randomize(curv_tl, n_rep=n_surrogate_datasets)
+
+
+###############################################################################
+# Similar to the previous case, we can now plot the results:
+
+import matplotlib.pyplot as plt
+from scipy.stats import spearmanr
+fig, axs = plt.subplots(1, 2, figsize=(9, 3.5))
+
+feats = {'t1wt2w': t1wt2w_tl, 'curvature': curv_tl}
+rand = {'t1wt2w': t1wt2w_surrogates, 'curvature': curv_surrogates}
+
+for k, (fn, data) in enumerate(rand.items()):
+    r_obs, pv_obs = spearmanr(feats[fn], embedding_tl, nan_policy='omit')
+
+    # Compute perm pval
+    r_rand = np.asarray([spearmanr(embedding_tl, d)[0] for d in data])
+    pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))
+
+    # Plot null dist
+    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
+    axs[k].axvline(r_obs, lw=2, ls='--', color='k')
+    axs[k].set_xlabel(f'Correlation with {fn}')
+    if k == 0:
+        axs[k].set_ylabel('Density')
+
+    print(f'{fn.capitalize()}:\n Obs      : {pv_obs:.5e}\n '
+          f'Variogram: {pv_rand:.5e}\n')
+
+fig.tight_layout()
+plt.show()
diff --git a/brainspace/gradient/embedding.py b/brainspace/gradient/embedding.py
index 05d4e16..2b34276 100644
--- a/brainspace/gradient/embedding.py
+++ b/brainspace/gradient/embedding.py
@@ -305,7 +305,7 @@ class DiffusionMaps(Embedding):
 
     """
 
-    def __init__(self, n_components=10, alpha=0.5, diffusion_time=1,
+    def __init__(self, n_components=10, alpha=0.5, diffusion_time=0,
                  random_state=None):
         super().__init__(n_components=n_components)
         self.alpha = alpha
diff --git a/brainspace/mesh/mesh_io.py b/brainspace/mesh/mesh_io.py
index b1e3f76..5c0e2d6 100755
--- a/brainspace/mesh/mesh_io.py
+++ b/brainspace/mesh/mesh_io.py
@@ -5,7 +5,9 @@
 # Author: Oualid Benkarim <oualid.benkarim@mcgill.ca>
 # License: BSD 3 clause
 
-
+import tempfile
+import gzip
+import os
 from vtk import (vtkPLYReader, vtkPLYWriter, vtkXMLPolyDataReader,
                  vtkXMLPolyDataWriter, vtkPolyDataReader, vtkPolyDataWriter)
 
@@ -70,7 +72,8 @@ def read_surface(ipth, itype=None, return_data=True, update=True):
     ipth : str
         Input filename.
     itype : {'ply', 'vtp', 'vtk', 'fs', 'asc', 'gii'}, optional
-        Input file type. If None, it is deduced from `ipth`. Default is None.
+        Input file type. If None, it is deduced from `ipth`. Files compressed
+        with gzip (.gz) are also supported. Default is None.
     return_data : bool, optional
         Whether to return data instead of filter. Default is False
     update : bool, optional
@@ -96,6 +99,16 @@ def read_surface(ipth, itype=None, return_data=True, update=True):
     if itype is None:
         itype = ipth.split('.')[-1]
 
+    if itype == 'gz':
+        extension = ipth.split('.')[-2]
+        tmp = tempfile.NamedTemporaryFile(suffix='.' + extension, delete = False)
+        tmp_name = tmp.name
+        tmp.close() # Close and reopen because windows throws permission errors when both reading and writing.
+        _uncompress(ipth, tmp_name)
+        result = read_surface(tmp_name, extension, return_data=return_data, update=update)
+        os.unlink(tmp_name)
+        return result
+
     reader = _select_reader(itype)
     reader.filename = ipth
 
@@ -165,3 +178,27 @@ def convert_surface(ipth, opth, itype=None, otype=None, oformat=None):
     """
     reader = read_surface(ipth, itype=itype, return_data=False, update=False)
     write_surface(reader, opth, oformat=oformat, otype=otype)
+
+def _uncompress(ipth, opth, block_size=65536):
+    """Uncompresses files. Currently only supports gzip.
+
+    Parameters
+    ----------
+    ipth : str
+        Input filename.
+    opth : str
+        Output filename.
+    block_size : int, optional
+        Size of blocks of the input that are read at a time, by default 65536
+    
+    """
+    if ipth.split('.')[-1] == 'gz':
+        with gzip.open(ipth, 'rb') as i_file, open(opth, 'wb') as o_file:
+            while True:
+                block = i_file.read(block_size)
+                if not block:
+                    break
+                else:
+                    o_file.write(block)
+    else:
+        ValueError('Unknown file format.')
\ No newline at end of file
diff --git a/brainspace/null_models/__init__.py b/brainspace/null_models/__init__.py
index e1b666c..406a8a6 100644
--- a/brainspace/null_models/__init__.py
+++ b/brainspace/null_models/__init__.py
@@ -1,10 +1,12 @@
 from .moran import (compute_mem, moran_randomization,
                     MoranRandomization)
 from .spin import spin_permutations, SpinPermutations
-
+from .variogram import SurrogateMaps, SampledSurrogateMaps
 
 __all__ = ['SpinPermutations',
            'MoranRandomization',
            'compute_mem',
            'moran_randomization',
-           'spin_permutations']
+           'spin_permutations',
+           'SurrogateMaps',
+           'SampledSurrogateMaps']
diff --git a/brainspace/null_models/variogram.py b/brainspace/null_models/variogram.py
new file mode 100644
index 0000000..a32fadb
--- /dev/null
+++ b/brainspace/null_models/variogram.py
@@ -0,0 +1,2005 @@
+"""
+Implementation of variogram-matching procedure.
+"""
+
+# Author: Joshua Burt <joshua.burt@yale.edu>
+# License: BSD 3 clause
+
+import numpy as np
+import numpy.lib.format
+from pathlib import Path
+from sklearn.base import BaseEstimator
+from sklearn.linear_model import LinearRegression
+
+
+# ----------------------
+# ------ Checks --------
+# ----------------------
+
+def is_string_like(obj):
+    """ Check whether `obj` behaves like a string. """
+    try:
+        obj + ''
+    except (TypeError, ValueError):
+        return False
+    return True
+
+
+def check_map(x):
+    """
+    Check that brain map is array_like and one dimensional.
+
+    Parameters
+    ----------
+    x : 1D ndarray
+        Brain map
+
+    Returns
+    -------
+    None
+
+    Raises
+    ------
+    TypeError : `x` is not a ndarray object
+    ValueError : `x` is not one-dimensional
+
+    """
+    if not isinstance(x, np.ndarray):
+        e = "Brain map must be array-like\n"
+        e += "got type {}".format(type(x))
+        raise TypeError(e)
+    if x.ndim != 1:
+        e = "Brain map must be one-dimensional\n"
+        e += "got shape {}".format(x.shape)
+        raise ValueError(e)
+
+
+def check_pv(pv):
+    """
+    Check input argument `pv`.
+
+    Parameters
+    ----------
+    pv : int
+        Percentile of the pairwise distance distribution at which to truncate
+        during variogram fitting.
+
+    Returns
+    -------
+    int
+
+    Raises
+    ------
+    ValueError : `pv` lies outside range (0, 100]
+
+    """
+    try:
+        pv = int(pv)
+    except ValueError:
+        raise ValueError("parameter 'pv' must be an integer in (0,100]")
+    if pv <= 0 or pv > 100:
+        raise ValueError("parameter 'pv' must be in (0,100]")
+    return pv
+
+
+def check_deltas(deltas):
+    """
+    Check input argument `deltas`.
+
+    Parameters
+    ----------
+    deltas : 1D ndarray or List[float]
+        Proportions of neighbors to include for smoothing, in (0, 1]
+
+    Returns
+    -------
+    None
+
+    Raises
+    ------
+    TypeError : `deltas` is not a List or ndarray object
+    ValueError : One or more elements of `deltas` lies outside (0,1]
+
+    """
+    if not isinstance(deltas, list) and not isinstance(deltas, np.ndarray):
+        raise TypeError("Parameter `deltas` must be a list or ndarray")
+    for d in deltas:
+        if d <= 0 or d > 1:
+            raise ValueError("Each element of `deltas` must lie in (0,1]")
+
+
+def count_lines(filename):
+    """
+    Count number of lines in a file.
+
+    Parameters
+    ----------
+    filename : filename
+
+    Returns
+    -------
+    int
+        number of lines in file
+
+    """
+    with open(filename, 'rb') as f:
+        lines = 0
+        buf_size = 1024 * 1024
+        read_f = f.raw.read
+        buf = read_f(buf_size)
+        while buf:
+            lines += buf.count(b'\n')
+            buf = read_f(buf_size)
+        return lines
+
+# ----------------------
+# ------ Data I/O ------
+# ----------------------
+
+
+def dataio(x):
+    """
+    Data I/O for core classes.
+
+    To facilitate flexible user inputs, this function loads data from:
+        - txt files
+        - npy files (memory-mapped arrays)
+        - array_like data
+
+    Parameters
+    ----------
+    x : filename or ndarray or np.memmap
+
+    Returns
+    -------
+    ndarray or np.memmap
+
+    Raises
+    ------
+    FileExistsError : file does not exist
+    RuntimeError : file is empty
+    ValueError : file type cannot be determined by file extension
+    TypeError : input is not a filename or array_like object
+
+    """
+    if is_string_like(x):
+        if not Path(x).exists():
+            raise FileExistsError("file does not exist: {}".format(x))
+        elif Path(x).stat().st_size == 0:
+            raise RuntimeError("file is empty: {}".format(x))
+        elif Path(x).suffix == ".npy":  # memmap
+            return np.load(x, mmap_mode='r')
+        elif Path(x).suffix == ".txt":  # text file
+            return np.loadtxt(x).squeeze()
+        else:
+            raise ValueError(
+                "expected npy or txt file, got {}".format(Path(x).suffix))
+    else:
+        if not isinstance(x, np.ndarray):
+            raise TypeError(
+                "expected filename or array_like obj, got {}".format(type(x)))
+        return x
+
+
+def txt2memmap(dist_file, output_dir, maskfile=None, delimiter=' '):
+    """
+    Export distance matrix to memory-mapped array.
+
+    Parameters
+    ----------
+    dist_file : filename
+        Path to `delimiter`-separated distance matrix file
+    output_dir : filename
+        Path to directory in which output files will be written
+    maskfile : filename or ndarray or None, default None
+        Path to a neuroimaging/txt file containing a mask, or a mask
+        represented as a numpy array. Mask scalars are cast to boolean, and
+        all elements not equal to zero will be masked.
+    delimiter : str
+        Delimiting character in `dist_file`
+
+    Returns
+    -------
+    dict
+        Keys are 'D' and 'index'; values are absolute paths to the
+        corresponding binary files on disk.
+
+    Notes
+    -----
+    Each row of the distance matrix is sorted before writing to file. Thus, a
+    second mem-mapped array is necessary, the i-th row of which contains
+    argsort(d[i]).
+    If `maskfile` is not None, a binary mask.txt file will also be written to
+    the output directory.
+
+    Raises
+    ------
+    IOError : `output_dir` doesn't exist
+    ValueError : Mask image and distance matrix have inconsistent sizes
+
+    """
+    op = Path(output_dir)
+    nlines = count_lines(dist_file)
+    if not op.exists():
+        raise IOError("Output directory does not exist: {}".format(output_dir))
+
+    # Load mask if one was provided
+    if maskfile is not None:
+        mask = dataio(maskfile).astype(bool)
+        if mask.size != nlines:
+            e = "Incompatible input sizes\n"
+            e += "{} rows in {}\n".format(nlines, dist_file)
+            e += "{} elements in {}".format(mask.size, maskfile)
+            raise ValueError(e)
+        mask_fileout = str(op.joinpath("mask.txt"))
+        np.savetxt(  # Write to text file
+            fname=mask_fileout, X=mask.astype(int), fmt="%i", delimiter=',')
+        nv = int((~mask).sum())  # number of non-masked elements
+        idx = np.arange(nlines)[~mask]  # indices of non-masked elements
+    else:
+        nv = nlines
+        idx = np.arange(nlines)
+
+    # Build memory-mapped arrays
+    with open(dist_file, 'r') as fp:
+
+        npydfile = str(op.joinpath("distmat.npy"))
+        npyifile = str(op.joinpath("index.npy"))
+        fpd = numpy.lib.format.open_memmap(
+            npydfile, mode='w+', dtype=np.float32, shape=(nv, nv))
+        fpi = numpy.lib.format.open_memmap(
+            npyifile, mode='w+', dtype=np.int32, shape=(nv, nv))
+
+        ifp = 0  # Build memory-mapped arrays one row of distances at a time
+        for il, l in enumerate(fp):  # Loop over lines of file
+            if il not in idx:  # Keep only CIFTI vertices
+                continue
+            else:
+                line = l.rstrip()
+                if line:
+                    data = np.array(line.split(delimiter), dtype=np.float32)
+                    if data.size != nlines:
+                        raise RuntimeError(
+                            "Distance matrix is not square: {}".format(
+                                dist_file))
+                    d = data[idx]
+                    sort_idx = np.argsort(d)
+                    fpd[ifp, :] = d[sort_idx]  # sorted row of distances
+                    fpi[ifp, :] = sort_idx  # sort indexes
+                    ifp += 1
+        del fpd  # Flush memory changes to disk
+        del fpi
+
+    return {'distmat': npydfile, 'index': npyifile}  # Return filenames
+
+# ----------------------
+# - Smoothing kernels --
+# ----------------------
+
+
+def gaussian(d):
+    """
+    Gaussian kernel which truncates at one standard deviation.
+
+    Parameters
+    ----------
+    d : ndarray, shape (N,) or (M,N)
+        one- or two-dimensional array of distances
+
+    Returns
+    -------
+    ndarray, shape (N,) or (M,N)
+        Gaussian kernel weights
+
+    Raises
+    ------
+    TypeError : `d` is not array_like
+
+    """
+    try:  # 2-dim
+        return np.exp(-1.25 * np.square(d / d.max(axis=-1)[:, np.newaxis]))
+    except IndexError:  # 1-dim
+        return np.exp(-1.25 * np.square(d/d.max()))
+    except AttributeError:
+        raise TypeError("expected array_like, got {}".format(type(d)))
+
+
+def exp(d):
+    """
+    Exponentially decaying kernel which truncates at e^{-1}.
+
+    Parameters
+    ----------
+    d : ndarray, shape (N,) or (M,N)
+        one- or two-dimensional array of distances
+
+    Returns
+    -------
+    ndarray, shape (N,) or (M,N)
+        Exponential kernel weights
+
+    Notes
+    -----
+    Characteristic length scale is set to d.max(axis=-1), i.e. the maximum
+    distance within each row.
+
+    Raises
+    ------
+    TypeError : `d` is not array_like
+
+    """
+    try:  # 2-dim
+        return np.exp(-d / d.max(axis=-1)[:, np.newaxis])
+    except IndexError:  # 1-dim
+        return np.exp(-d/d.max())
+    except AttributeError:
+        raise TypeError("expected array_like, got {}".format(type(d)))
+
+
+def invdist(d):
+    """
+    Inverse distance kernel.
+
+    Parameters
+    ----------
+    d : ndarray, shape (N,) or (M,N)
+        One- or two-dimensional array of distances
+
+    Returns
+    -------
+    ndarray, shape (N,) or (M,N)
+        Inverse distance, i.e. d^{-1}
+
+    Raises
+    ------
+    ZeroDivisionError : `d` includes zero value
+    TypeError : `d` is not array_like
+
+    """
+    try:
+        return 1. / d
+    except ZeroDivisionError as e:
+        raise ZeroDivisionError(e)
+    except AttributeError:
+        raise TypeError("expected array_like, got {}".format(type(d)))
+
+
+def uniform(d):
+    """
+    Uniform (i.e., distance independent) kernel.
+
+    Parameters
+    ----------
+    d : ndarray, shape (N,) or (M,N)
+        One- or two-dimensional array of distances
+
+    Returns
+    -------
+    ndarray, shape (N,) or (M,N)
+        Uniform kernel weights
+
+    Notes
+    -----
+    Each element is normalized to 1/N such that columns sum to unity.
+
+    Raises
+    ------
+    TypeError : `d` is not array_like
+
+    """
+    try:  # 2-dim
+        return np.ones(d.shape) / d.shape[-1]
+    except IndexError:  # 1-dim
+        return np.ones(d.size) / d.size
+    except AttributeError:
+        raise TypeError("expected array_like, got {}".format(type(d)))
+
+
+def check_kernel(kernel):
+    """
+    Check that a valid kernel was specified and return callable.
+
+    Parameters
+    ----------
+    kernel : 'exp' or 'gaussian' or 'invdist' or 'uniform'
+        Kernel selection
+
+    Returns
+    -------
+    Callable
+
+    Raises
+    ------
+    NotImplementedError : kernel is not implemented
+
+    """
+    kernels = {'exp': exp,
+               'gaussian': gaussian,
+               'invdist': invdist,
+               'uniform': uniform}
+    if kernel not in kernels.keys():
+        e = "'{}' is not a valid kernel\n".format(kernel)
+        e += "Valid kernels: {}".format(", ".join([k for k in kernels.keys()]))
+        raise NotImplementedError(e)
+    return kernels[kernel]
+
+# ----------------------
+# ---- Core classes ----
+# ----------------------
+
+
+# class Base:
+#     """ Base implementation of map generator.
+#
+#     Parameters
+#     ----------
+#     x : filename or 1D ndarray
+#         Target brain map
+#     D : filename or ndarray, shape (N,N)
+#         Pairwise distance matrix
+#     deltas : 1D ndarray or List[float], default [0.1,0.2,...,0.9]
+#         Proportion of neighbors to include for smoothing, in (0, 1]
+#     kernel : str, default 'exp'
+#         Kernel with which to smooth permuted maps:
+#           'gaussian' : Gaussian function.
+#           'exp' : Exponential decay function.
+#           'invdist' : Inverse distance.
+#           'uniform' : Uniform weights (distance independent).
+#     pv : int, default 25
+#         Percentile of the pairwise distance distribution at which to
+#         truncate during variogram fitting
+#     nh : int, default 25
+#         Number of uniformly spaced distances at which to compute variogram
+#     resample : bool, default False
+#         Resample surrogate maps' values from target brain map
+#     b : float or None, default None
+#         Gaussian kernel bandwidth for variogram smoothing. If None, set to
+#         three times the spacing between variogram x-coordinates.
+#
+#     Notes
+#     -----
+#     Passing resample=True preserves the distribution of values in the target
+#     map, with the possibility of worsening the simulated surrogate maps'
+#     variograms fits.
+#
+#     """
+#
+#     def __init__(self, x, D, deltas=np.linspace(0.1, 0.9, 9),
+#                  kernel='exp', pv=25, nh=25, resample=False, b=None):
+#
+#         self.x = x
+#         self.D = D
+#         n = self._x.size
+#         self.resample = resample
+#         self.nh = nh
+#         self.deltas = deltas
+#         self.pv = pv
+#         self.nmap = n
+#         self.kernel = kernel  # Smoothing kernel selection
+#         self._ikn = np.arange(n)[:, None]
+#         self._triu = np.triu_indices(self._nmap, k=1)  # upper triangular inds
+#         self._u = self._D[self._triu]  # variogram X-coordinate
+#         self._v = self.compute_variogram(self._x)  # variogram Y-coord
+#
+#         # Get indices of pairs with u < pv'th percentile
+#         self._uidx = np.where(self._u < np.percentile(self._u, self._pv))[0]
+#         self._uisort = np.argsort(self._u[self._uidx])
+#
+#         # Find sorted indices of first `kmax` elements of each row of dist. mat.
+#         self._disort = np.argsort(self._D, axis=-1)
+#         self._jkn = dict.fromkeys(deltas)
+#         self._dkn = dict.fromkeys(deltas)
+#         for delta in deltas:
+#             k = int(delta*n)
+#             # find index of k nearest neighbors for each area
+#             self._jkn[delta] = self._disort[:, 1:k+1]  # prevent self-coupling
+#             # find distance to k nearest neighbors for each area
+#             self._dkn[delta] = self._D[(self._ikn, self._jkn[delta])]
+#
+#         # Smoothed variogram and variogram _b
+#         utrunc = self._u[self._uidx]
+#         self._h = np.linspace(utrunc.min(), utrunc.max(), self._nh)
+#         self.b = b
+#         self._smvar = self.smooth_variogram(self._v)
+#
+#         # Linear regression model
+#         self._lm = LinearRegression(fit_intercept=True)
+#
+#     def __call__(self, n=1):
+#         """
+#         Randomly generate new surrogate map(s).
+#
+#         Parameters
+#         ----------
+#         n : int, default 1
+#             Number of surrogate maps to randomly generate
+#
+#         Returns
+#         -------
+#         ndarray, shape (n,N)
+#             Randomly generated map(s) with matched spatial autocorrelation
+#
+#         Notes
+#         -----
+#         Chooses a level of smoothing that produces a smoothed variogram which
+#         best approximates the true smoothed variogram. Selecting resample='True'
+#         preserves the original map's value distribution at the expense of
+#         worsening the surrogate maps' variogram fit.
+#
+#         """
+#         print("Generating {} maps...".format(n))
+#         surrs = np.empty((n, self._nmap))
+#         for i in range(n):  # generate random maps
+#
+#             xperm = self.permute_map()  # Randomly permute values
+#             res = dict.fromkeys(self._deltas)
+#
+#             for delta in self.deltas:  # foreach neighborhood size
+#                 # Smooth the permuted map using delta proportion of
+#                 # neighbors to reintroduce spatial autocorrelation
+#                 sm_xperm = self.smooth_map(x=xperm, delta=delta)
+#
+#                 # Calculate empirical variogram of the smoothed permuted map
+#                 vperm = self.compute_variogram(sm_xperm)
+#
+#                 # Calculate smoothed variogram of the smoothed permuted map
+#                 smvar_perm = self.smooth_variogram(vperm)
+#
+#                 # Fit linear regression btwn smoothed variograms
+#                 res[delta] = self.regress(smvar_perm, self._smvar)
+#
+#             alphas, betas, residuals = np.array(
+#                 [res[d] for d in self._deltas]).T
+#
+#             # Select best-fit model and regression parameters
+#             iopt = np.argmin(residuals)
+#             dopt = self._deltas[iopt]
+#             aopt = alphas[iopt]
+#             bopt = betas[iopt]
+#
+#             # Transform and smooth permuted map using best-fit parameters
+#             sm_xperm_best = self.smooth_map(x=xperm, delta=dopt)
+#             surr = (np.sqrt(np.abs(bopt)) * sm_xperm_best +
+#                     np.sqrt(np.abs(aopt)) * np.random.randn(self._nmap))
+#             surrs[i] = surr
+#
+#         if self._resample:  # resample values from empirical map
+#             sorted_map = np.sort(self._x)
+#             for i, surr in enumerate(surrs):
+#                 ii = np.argsort(surr)
+#                 np.put(surr, ii, sorted_map)
+#
+#         return surrs.squeeze()
+#
+#     def compute_variogram(self, x):
+#         """
+#         Compute variogram values (i.e., one-half squared pairwise differences).
+#
+#         Parameters
+#         ----------
+#         x : 1D ndarray
+#             Brain map scalar array
+#
+#         Returns
+#         -------
+#         v : ndarray, shape (N(N-1)/2,)
+#            Variogram y-coordinates, i.e. 0.5 * (x_i - x_j) ^ 2
+#
+#         """
+#         diff_ij = np.subtract.outer(x, x)
+#         v = 0.5 * np.square(diff_ij)[self._triu]
+#         return v
+#
+#     def permute_map(self):
+#         """
+#         Return randomly permuted brain map.
+#
+#         Returns
+#         -------
+#         1D ndarray
+#             Random permutation of target brain map
+#
+#         """
+#         perm_idx = np.random.permutation(np.arange(self._x.size))
+#         mask_perm = self._x.mask[perm_idx]
+#         x_perm = self._x.data[perm_idx]
+#         return np.ma.masked_array(data=x_perm, mask=mask_perm)
+#
+#     def smooth_map(self, x, delta):
+#         """
+#         Smooth `x` using `delta` proportion of nearest neighbors.
+#
+#         Parameters
+#         ----------
+#         x : 1D ndarray
+#             Brain map scalars
+#         delta : float
+#             Proportion of neighbors to include for smoothing, in (0, 1)
+#
+#         Returns
+#         -------
+#         1D ndarray
+#             Smoothed brain map
+#
+#         """
+#         # Values of k nearest neighbors for each brain area
+#         xkn = x[self._jkn[delta]]
+#         weights = self._kernel(self._dkn[delta])  # Distance-weight kernel
+#         # Kernel-weighted sum
+#         return (weights * xkn).sum(axis=1) / weights.sum(axis=1)
+#
+#     def smooth_variogram(self, v, return_h=False):
+#         """
+#         Smooth a variogram.
+#
+#         Parameters
+#         ----------
+#         v : 1D ndarray
+#             Variogram values, i.e. 0.5 * (x_i - x_j) ^ 2
+#         return_h : bool, default False
+#             Return distances at which the smoothed variogram was computed
+#
+#         Returns
+#         -------
+#         1D ndarray, shape (nh,)
+#             Smoothed variogram values
+#         1D ndarray, shape (nh,)
+#             Distances at which smoothed variogram was computed (returned only if
+#             `return_h` is True)
+#
+#         Raises
+#         ------
+#         ValueError : `v` has unexpected size.
+#
+#         """
+#         u = self._u[self._uidx]
+#         v = v[self._uidx]
+#         if len(u) != len(v):
+#             raise ValueError(
+#                 "argument v: expected size {}, got {}".format(len(u), len(v)))
+#         # Subtract each h from each pairwise distance u
+#         # Each row corresponds to a unique h
+#         du = np.abs(u - self._h[:, None])
+#         w = np.exp(-np.square(2.68 * du / self._b) / 2)
+#         denom = w.sum(axis=1)
+#         wv = w * v[None, :]
+#         num = wv.sum(axis=1)
+#         output = num / denom
+#         if not return_h:
+#             return output
+#         return output, self._h
+#
+#     def regress(self, x, y):
+#         """
+#         Linearly regress `x` onto `y`.
+#
+#         Parameters
+#         ----------
+#         x : 1D ndarray
+#             Independent variable
+#         y : 1D ndarray
+#             Dependent variable
+#
+#         Returns
+#         -------
+#         alpha : float
+#             Intercept term (offset parameter)
+#         beta : float
+#             Regression coefficient (scale parameter)
+#         res : float
+#             Sum of squared residuals
+#
+#         """
+#         self._lm.fit(X=np.expand_dims(x, -1), y=y)
+#         beta = self._lm.coef_
+#         alpha = self._lm.intercept_
+#         y_pred = self._lm.predict(X=np.expand_dims(x, -1))
+#         res = np.sum(np.square(y-y_pred))
+#         return alpha, beta, res
+#
+#     @property
+#     def x(self):
+#         """ 1D ndarray : brain map scalar array """
+#         return self._x
+#
+#     @x.setter
+#     def x(self, x):
+#         x_ = dataio(x)
+#         check_map(x=x_)
+#         brain_map = np.ma.masked_array(data=x_, mask=np.isnan(x_))
+#         self._x = brain_map
+#
+#     @property
+#     def D(self):
+#         """ ndarray, shape (N,N) : Pairwise distance matrix """
+#         return self._D
+#
+#     @D.setter
+#     def D(self, x):
+#         d_ = dataio(x)
+#         if not np.allclose(d_, d_.T):
+#             raise ValueError("Distance matrix must be symmetric")
+#         n = self._x.size
+#         if d_.shape != (n, n):
+#             e = "Distance matrix must have dimensions consistent with brain map"
+#             e += "\nDistance matrix shape: {}".format(d_.shape)
+#             e += "\nBrain map size: {}".format(n)
+#             raise ValueError(e)
+#         self._D = d_
+#
+#     @property
+#     def nmap(self):
+#         """ int : length of brain map """
+#         return self._nmap
+#
+#     @nmap.setter
+#     def nmap(self, x):
+#         self._nmap = int(x)
+#
+#     @property
+#     def pv(self):
+#         """ int : percentile of pairwise distances at which to truncate """
+#         return self._pv
+#
+#     @pv.setter
+#     def pv(self, x):
+#         pv = check_pv(x)
+#         self._pv = pv
+#
+#     @property
+#     def deltas(self):
+#         """ 1D ndarray or List[float] : proportions of nearest neighbors """
+#         return self._deltas
+#
+#     @deltas.setter
+#     def deltas(self, x):
+#         check_deltas(deltas=x)
+#         self._deltas = x
+#
+#     @property
+#     def nh(self):
+#         """ int : number of variogram distance intervals """
+#         return self._nh
+#
+#     @nh.setter
+#     def nh(self, x):
+#         self._nh = x
+#
+#     @property
+#     def h(self):
+#         """ 1D ndarray : distances at which smoothed variogram is computed """
+#         return self._h
+#
+#     @property
+#     def kernel(self):
+#         """ Callable : smoothing kernel function """
+#         return self._kernel
+#
+#     @kernel.setter
+#     def kernel(self, x):
+#         kernel_callable = check_kernel(x)
+#         self._kernel = kernel_callable
+#
+#     @property
+#     def resample(self):
+#         """ bool : whether to resample surrogate maps from target map """
+#         return self._resample
+#
+#     @resample.setter
+#     def resample(self, x):
+#         if not isinstance(x, bool):
+#             e = "parameter `resample`: expected bool, got {}".format(type(x))
+#             raise TypeError(e)
+#         self._resample = x
+#
+#     @property
+#     def b(self):
+#         """ numeric : Gaussian kernel bandwidth """
+#         return self._b
+#
+#     @b.setter
+#     def b(self, x):
+#         if x is not None:
+#             try:
+#                 self._b = float(x)
+#             except (ValueError, TypeError):
+#                 e = "bandwidth b: expected numeric, got {}".format(type(x))
+#                 raise ValueError(e)
+#         else:   # set bandwidth equal to 3x bin spacing
+#             self._b = 3.*np.mean(self._h[1:] - self._h[:-1])
+#
+#
+# class Sampled:
+#     """
+#     Sampling implementation of map generator.
+#
+#     Parameters
+#     ----------
+#     x : 1D ndarray
+#         Target brain map
+#     D : ndarray or memmap, shape (N,N)
+#         Pairwise distance matrix between elements of `x`. Each row of `D` should
+#         be sorted. Indices used to sort each row are passed to the `index`
+#         argument. See :func:`brainsmash.mapgen.memmap.txt2memmap` or the online
+#         documentation for more details (brainsmash.readthedocs.io)
+#     index : filename or ndarray or memmap, shape(N,N)
+#         See above
+#     ns : int, default 500
+#         Take a subsample of `ns` rows from `D` when fitting variograms
+#     deltas : ndarray or List[float], default [0.3, 0.5, 0.7, 0.9]
+#         Proportions of neighbors to include for smoothing, in (0, 1]
+#     kernel : str, default 'exp'
+#         Kernel with which to smooth permuted maps
+#         - 'gaussian' : gaussian function
+#         - 'exp' : exponential decay function
+#         - 'invdist' : inverse distance
+#         - 'uniform' : uniform weights (distance independent)
+#     pv : int, default 70
+#         Percentile of the pairwise distance distribution (in `D`) at
+#         which to truncate during variogram fitting
+#     nh : int, default 25
+#         Number of uniformly spaced distances at which to compute variogram
+#     knn : int, default 1000
+#         Number of nearest regions to keep in the neighborhood of each region
+#     b : float or None, default None
+#         Gaussian kernel bandwidth for variogram smoothing. if None,
+#         three times the distance interval spacing is used.
+#     resample : bool, default False
+#         Resample surrogate map values from the target brain map
+#     verbose : bool, default False
+#         Print surrogate count each time new surrogate map created
+#
+#     Notes
+#     -----
+#     Passing resample=True will preserve the distribution of values in the
+#     target map, at the expense of worsening simulated surrogate maps'
+#     variograms fits. This worsening will increase as the empirical map
+#     more strongly deviates from normality.
+#
+#     Raises
+#     ------
+#     ValueError : `x` and `D` have inconsistent sizes
+#
+#     """
+#
+#     def __init__(self, x, D, index, ns=500, pv=70, nh=25, knn=1000, b=None,
+#                  deltas=np.arange(0.3, 1., 0.2), kernel='exp', resample=False,
+#                  verbose=False):
+#
+#         self._verbose = verbose
+#         self.x = x
+#         n = self._x.size
+#         self.nmap = int(n)
+#         self.knn = knn
+#         self.D = D
+#         self.index = index
+#         self.resample = resample
+#         self.nh = int(nh)
+#         self.deltas = deltas
+#         self.ns = int(ns)
+#         self.b = b
+#         self.pv = pv
+#         self._ikn = np.arange(self._nmap)[:, None]
+#
+#         # Store k nearest neighbors from distance and index matrices
+#         self.kernel = kernel  # Smoothing kernel selection
+#         self._dmax = np.percentile(self._D, self._pv)
+#         self.h = np.linspace(self._D.min(), self._dmax, self._nh)
+#         if not self._b:
+#             self.b = 3 * (self.h[1] - self.h[0])
+#
+#         # Linear regression model
+#         self._lm = LinearRegression(fit_intercept=True)
+#
+#     def __call__(self, n=1):
+#         """
+#         Randomly generate new surrogate map(s).
+#
+#         Parameters
+#         ----------
+#         n : int, default 1
+#             Number of surrogate maps to randomly generate
+#
+#         Returns
+#         -------
+#         ndarray, shape (n,N)
+#             Randomly generated map(s) with matched spatial autocorrelation
+#
+#         Notes
+#         -----
+#         Chooses a level of smoothing that produces a smoothed variogram which
+#         best approximates the true smoothed variogram. Selecting resample='True'
+#         preserves the map value distribution at the expense of worsening the
+#         surrogate maps' variogram fits.
+#
+#         """
+#         if self._verbose:
+#             print("Generating {} maps...".format(n))
+#         surrs = np.empty((n, self._nmap))
+#         for i in range(n):  # generate random maps
+#             if self._verbose:
+#                 print(i+1)
+#
+#             # Randomly permute map
+#             x_perm = self.permute_map()
+#
+#             # Randomly select subset of regions to use for variograms
+#             idx = self.sample()
+#
+#             # Compute empirical variogram
+#             v = self.compute_variogram(self._x, idx)
+#
+#             # Variogram ordinates; use nearest neighbors because local effect
+#             u = self._D[idx, :]
+#             uidx = np.where(u < self._dmax)
+#
+#             # Smooth empirical variogram
+#             smvar, u0 = self.smooth_variogram(u[uidx], v[uidx], return_h=True)
+#
+#             res = dict.fromkeys(self._deltas)
+#
+#             for d in self._deltas:  # foreach neighborhood size
+#
+#                 k = int(d * self._knn)
+#
+#                 # Smooth the permuted map using k nearest neighbors to
+#                 # reintroduce spatial autocorrelation
+#                 sm_xperm = self.smooth_map(x=x_perm, k=k)
+#
+#                 # Calculate variogram values for the smoothed permuted map
+#                 vperm = self.compute_variogram(sm_xperm, idx)
+#
+#                 # Calculate smoothed variogram of the smoothed permuted map
+#                 smvar_perm = self.smooth_variogram(u[uidx], vperm[uidx])
+#
+#                 # Fit linear regression btwn smoothed variograms
+#                 res[d] = self.regress(smvar_perm, smvar)
+#
+#             alphas, betas, residuals = np.array(
+#                 [res[d] for d in self._deltas]).T
+#
+#             # Select best-fit model and regression parameters
+#             iopt = np.argmin(residuals)
+#             dopt = self._deltas[iopt]
+#             self._dopt = dopt
+#             kopt = int(dopt * self._knn)
+#             aopt = alphas[iopt]
+#             bopt = betas[iopt]
+#
+#             # Transform and smooth permuted map using best-fit parameters
+#             sm_xperm_best = self.smooth_map(x=x_perm, k=kopt)
+#             surr = (np.sqrt(np.abs(bopt)) * sm_xperm_best +
+#                     np.sqrt(np.abs(aopt)) * np.random.randn(self._nmap))
+#             surrs[i] = surr
+#
+#         if self._resample:  # resample values from empirical map
+#             sorted_map = np.sort(self._x)
+#             for i, surr in enumerate(surrs):
+#                 ii = np.argsort(surr)
+#                 np.put(surr, ii, sorted_map)
+#
+#         if self._ismasked:
+#             return np.ma.masked_array(
+#                 data=surrs, mask=np.isnan(surrs)).squeeze()
+#         return surrs.squeeze()
+#
+#     def compute_variogram(self, x, idx):
+#         """
+#         Compute variogram of `x` using pairs of regions indexed by `idx`.
+#
+#         Parameters
+#         ----------
+#         x : 1Dndarray
+#             Brain map
+#         idx : ndarray[int], shape (ns,)
+#             Indices of randomly sampled brain regions
+#
+#         Returns
+#         -------
+#         v : ndarray, shape (ns,ns)
+#             Variogram y-coordinates, i.e. 0.5 * (x_i - x_j) ^ 2, for i,j in idx
+#
+#         """
+#         diff_ij = x[idx][:, None] - x[self._index[idx, :]]
+#         return 0.5 * np.square(diff_ij)
+#
+#     def permute_map(self):
+#         """
+#         Return a random permutation of the target brain map.
+#
+#         Returns
+#         -------
+#         1D ndarray
+#             Random permutation of target brain map
+#
+#         """
+#         perm_idx = np.random.permutation(self._nmap)
+#         if self._ismasked:
+#             mask_perm = self._x.mask[perm_idx]
+#             x_perm = self._x.data[perm_idx]
+#             return np.ma.masked_array(data=x_perm, mask=mask_perm)
+#         return self._x[perm_idx]
+#
+#     def smooth_map(self, x, k):
+#         """
+#         Smooth `x` using `k` nearest neighboring regions.
+#
+#         Parameters
+#         ----------
+#         x : 1D ndarray
+#             Brain map
+#         k : float
+#             Number of nearest neighbors to include for smoothing
+#
+#         Returns
+#         -------
+#         x_smooth : 1D ndarray
+#             Smoothed brain map
+#
+#         Notes
+#         -----
+#         Assumes `D` provided at runtime has been sorted.
+#
+#         """
+#         jkn = self._index[:, :k]  # indices of k nearest neighbors
+#         xkn = x[jkn]  # values of k nearest neighbors
+#         dkn = self._D[:, :k]  # distances to k nearest neighbors
+#         weights = self._kernel(dkn)  # distance-weighted kernel
+#         # Kernel-weighted sum
+#         return (weights * xkn).sum(axis=1) / weights.sum(axis=1)
+#
+#     def smooth_variogram(self, u, v, return_h=False):
+#         """
+#         Smooth a variogram.
+#
+#         Parameters
+#         ----------
+#         u : 1D ndarray
+#             Pairwise distances, ie variogram x-coordinates
+#         v : 1D ndarray
+#             Squared differences, ie variogram y-coordinates
+#         return_h : bool, default False
+#             Return distances at which smoothed variogram is computed
+#
+#         Returns
+#         -------
+#         ndarray, shape (nh,)
+#             Smoothed variogram samples
+#         ndarray, shape (nh,)
+#             Distances at which smoothed variogram was computed (returned if
+#             `return_h` is True)
+#
+#         Raises
+#         ------
+#         ValueError : `u` and `v` are not identically sized
+#
+#         """
+#         if len(u) != len(v):
+#             raise ValueError("u and v must have same number of elements")
+#
+#         # Subtract each element of h from each pairwise distance `u`.
+#         # Each row corresponds to a unique h.
+#         du = np.abs(u - self._h[:, None])
+#         w = np.exp(-np.square(2.68 * du / self._b) / 2)
+#         denom = w.sum(axis=1)
+#         wv = w * v[None, :]
+#         num = wv.sum(axis=1)
+#         output = num / denom
+#         if not return_h:
+#             return output
+#         return output, self._h
+#
+#     def regress(self, x, y):
+#         """
+#         Linearly regress `x` onto `y`.
+#
+#         Parameters
+#         ----------
+#         x : 1D ndarray
+#             Independent variable
+#         y : 1D ndarray
+#             Dependent variable
+#
+#         Returns
+#         -------
+#         alpha : float
+#             Intercept term (offset parameter)
+#         beta : float
+#             Regression coefficient (scale parameter)
+#         res : float
+#             Sum of squared residuals
+#
+#         """
+#         self._lm.fit(X=np.expand_dims(x, -1), y=y)
+#         beta = self._lm.coef_.item()
+#         alpha = self._lm.intercept_
+#         ypred = self._lm.predict(np.expand_dims(x, -1))
+#         res = np.sum(np.square(y-ypred))
+#         return alpha, beta, res
+#
+#     def sample(self):
+#         """
+#         Randomly sample (without replacement) brain areas for variogram
+#         computation.
+#
+#         Returns
+#         -------
+#         ndarray, shape (ns,)
+#             Indices of randomly sampled areas
+#
+#         """
+#         return np.random.choice(
+#             a=self._nmap, size=self._ns, replace=False).astype(np.int32)
+#
+#     @property
+#     def x(self):
+#         """1D ndarray : brain map scalars """
+#         if self._ismasked:
+#             return np.ma.copy(self._x)
+#         return np.copy(self._x)
+#
+#     @x.setter
+#     def x(self, x):
+#         self._ismasked = False
+#         x_ = dataio(x)
+#         check_map(x=x_)
+#         mask = np.isnan(x_)
+#         if mask.any():
+#             self._ismasked = True
+#             brain_map = np.ma.masked_array(data=x_, mask=mask)
+#         else:
+#             brain_map = x_
+#         self._x = brain_map
+#
+#     @property
+#     def D(self):
+#         """ndarray or memmap, shape (N,N) : Pairwise distance matrix """
+#         return np.copy(self._D)
+#
+#     @D.setter
+#     def D(self, x):
+#         x_ = dataio(x)
+#         n = self._x.size
+#         if x_.shape[0] != n:
+#             raise ValueError(
+#                 "D size along axis=0 must equal brain map size")
+#         self._D = x_[:, 1:self._knn + 1]  # prevent self-coupling
+#
+#     @property
+#     def index(self):
+#         """ndarray or memmap : indexes used to sort each row of dist. matrix """
+#         return np.copy(self._index)
+#
+#     @index.setter
+#     def index(self, x):
+#         x_ = dataio(x)
+#         n = self._x.size
+#         if x_.shape[0] != n:
+#             raise ValueError(
+#                 "index size along axis=0 must equal brain map size")
+#         self._index = x_[:, 1:self._knn+1].astype(np.int32)
+#
+#     @property
+#     def nmap(self):
+#         """ int : length of brain map """
+#         return self._nmap
+#
+#     @nmap.setter
+#     def nmap(self, x):
+#         self._nmap = int(x)
+#
+#     @property
+#     def pv(self):
+#         """ int : percentile of pairwise distances at which to truncate """
+#         return self._pv
+#
+#     @pv.setter
+#     def pv(self, x):
+#         pv = check_pv(x)
+#         self._pv = pv
+#
+#     @property
+#     def deltas(self):
+#         """ 1D ndarray or List[float] : proportions of nearest neighbors """
+#         return self._deltas
+#
+#     @deltas.setter
+#     def deltas(self, x):
+#         check_deltas(deltas=x)
+#         self._deltas = x
+#
+#     @property
+#     def nh(self):
+#         """ int : number of variogram distance intervals """
+#         return self._nh
+#
+#     @nh.setter
+#     def nh(self, x):
+#         self._nh = x
+#
+#     @property
+#     def kernel(self):
+#         """ Callable : smoothing kernel function
+#
+#         Notes
+#         -----
+#         When setting kernel, use name of kernel as defined in ``config.py``.
+#
+#         """
+#         return self._kernel
+#
+#     @kernel.setter
+#     def kernel(self, x):
+#         kernel_callable = check_kernel(x)
+#         self._kernel = kernel_callable
+#
+#     @property
+#     def resample(self):
+#         """ bool : whether to resample surrogate map values from target maps """
+#         return self._resample
+#
+#     @resample.setter
+#     def resample(self, x):
+#         if not isinstance(x, bool):
+#             raise TypeError("expected bool, got {}".format(type(x)))
+#         self._resample = x
+#
+#     @property
+#     def knn(self):
+#         """ int : number of nearest neighbors included in distance matrix """
+#         return self._knn
+#
+#     @knn.setter
+#     def knn(self, x):
+#         if x > self._nmap:
+#             raise ValueError('knn must be less than len(X)')
+#         self._knn = int(x)
+#
+#     @property
+#     def ns(self):
+#         """ int : number of randomly sampled regions used to construct map """
+#         return self._ns
+#
+#     @ns.setter
+#     def ns(self, x):
+#         self._ns = int(x)
+#
+#     @property
+#     def b(self):
+#         """ numeric : Gaussian kernel bandwidth """
+#         return self._b
+#
+#     @b.setter
+#     def b(self, x):
+#         self._b = x
+#
+#     @property
+#     def h(self):
+#         """ 1D ndarray : distances at which variogram is evaluated """
+#         return self._h
+#
+#     @h.setter
+#     def h(self, x):
+#         self._h = x
+
+
+class SurrogateMaps(BaseEstimator):
+    """ Spatial autocorrelation-preserving surrogate brain maps.
+
+    Parameters
+    ----------
+    deltas : 1D ndarray or List[float], optional
+        Proportion of neighbors to include for smoothing, in (0, 1]
+        Default is [0.1,0.2,...,0.9].
+    kernel : str, optional
+        Kernel with which to smooth permuted maps:
+          'gaussian' : Gaussian function.
+          'exp' : Exponential decay function.
+          'invdist' : Inverse distance.
+          'uniform' : Uniform weights (distance independent).
+        Default is 'exp'.
+    pv : int, optional
+        Percentile of the pairwise distance distribution at which to
+        truncate during variogram fitting. Default is 25.
+    nh : int, optional
+        Number of uniformly spaced distances at which to compute variogram.
+        Default is 25.
+    resample : bool, optional
+        Resample surrogate maps' values from target brain map.
+        Default is False.
+    b : float or None, optional
+        Gaussian kernel bandwidth for variogram smoothing. If None, set to
+        three times the spacing between variogram x-coordinates.
+        Default is None.
+    n_rep : int, optional
+        Number of randomizations (i.e., surrogate maps). Default is 100.
+    random_state : int or None, optional
+        Random state. Default is None.
+
+    See Also
+    --------
+    :class:`.SampledSurrogateMaps`
+    :class:`.SpinPermutations`
+    :class:`.MoranRandomization`
+
+    Notes
+    -----
+    Passing resample=True preserves the distribution of values in the target
+    map, with the possibility of worsening the simulated surrogate maps'
+    variograms fits.
+
+    """
+
+    def __init__(self, deltas=None, kernel='exp', pv=25, nh=25, resample=False,
+                 b=None, n_rep=100, random_state=None):
+
+        self.deltas = np.linspace(0.1, 0.9, 9) if deltas is None else deltas
+        check_deltas(deltas=self.deltas)
+
+        self.resample = resample
+        self.nh = nh
+        self.pv = check_pv(pv)
+        self.kernel = check_kernel(kernel)  # Smoothing kernel selection
+        self.b = b
+        self.n_rep = n_rep
+        self.random_state = random_state
+
+        self._rs = np.random.RandomState(self.random_state)
+
+        # Linear regression model
+        self._lm = LinearRegression(fit_intercept=True)
+
+    def _check_distance(self, dist):
+        d = dataio(dist)
+        if not np.allclose(d, d.T):
+            raise ValueError("Distance matrix must be symmetric")
+        return d
+
+    def _check_map(self, x):
+        x_ = dataio(x)
+        check_map(x=x_)
+        brain_map = np.ma.masked_array(data=x_, mask=np.isnan(x_))
+        if self._dist.shape[0] != brain_map.size:
+            e = "Brain map must have dimensions consistent with distance matrix"
+            e += "\nDistance matrix shape: {}".format(self._dist.shape)
+            e += "\nBrain map size: {}".format(brain_map.size)
+            raise ValueError(e)
+
+        return brain_map
+
+    def fit(self, dist):
+        """ Prepare data for sorrugate map generation..
+
+        Parameters
+        ----------
+        dist : filename or ndarray, shape (N,N)
+            Pairwise (geodesic) distance matrix.
+
+        Returns
+        -------
+        self : object
+            Returns self.
+
+        """
+
+        self._dist = self._check_distance(dist)
+        self.nmap = n = dist.shape[0]
+
+        self._ikn = np.arange(n)[:, None]
+        self._triu = np.triu_indices(n, k=1)  # upper triangular inds
+        self._u = self._dist[self._triu]  # variogram X-coordinate
+        # self._v = self.compute_variogram(self._x)  # variogram Y-coord
+
+        # Get indices of pairs with u < pv'th percentile
+        self._uidx = np.where(self._u < np.percentile(self._u, self.pv))[0]
+        self._uisort = np.argsort(self._u[self._uidx])
+
+        # Find sorted indices of first `kmax` elements of each row of dist.
+        self._disort = np.argsort(self._dist, axis=-1)
+        self._jkn = dict.fromkeys(self.deltas)
+        self._dkn = dict.fromkeys(self.deltas)
+        for delta in self.deltas:
+            k = int(delta*n)
+            # find index of k nearest neighbors for each area
+            self._jkn[delta] = self._disort[:, 1:k+1]  # prevent self-coupling
+            # find distance to k nearest neighbors for each area
+            self._dkn[delta] = self._dist[(self._ikn, self._jkn[delta])]
+
+        # Smoothed variogram and variogram _b
+        utrunc = self._u[self._uidx]
+        self._h = np.linspace(utrunc.min(), utrunc.max(), self.nh)
+        if self.b is None:
+            self.b = 3.*np.mean(self._h[1:] - self._h[:-1])
+
+        return self
+
+    def randomize(self, x, n_rep=None):
+        """ Generate surrogate maps from `x`.
+
+        Parameters
+        ----------
+        x : filename or 1D ndarray
+            Target brain map
+        n_rep : int or None, optional
+            Number of surrogates maps to randomly generate. If None, use
+            the default `n_rep`.
+
+        Returns
+        -------
+        output : ndarray, shape = (n_rep, n_verts)
+            Randomly generated map(s) with matched spatial autocorrelation.
+
+        """
+
+        x = self._check_map(x)
+        nmap = x.size
+
+        n_rep = self.n_rep if n_rep is None else n_rep
+
+        v = self.compute_variogram(x)  # variogram Y-coord
+        smvar = self.smooth_variogram(v)
+
+        surrs = np.empty((n_rep, nmap))
+        for i in range(n_rep):  # generate random maps
+
+            xperm = self.permute_map(x)  # Randomly permute values
+            res = dict.fromkeys(self.deltas)
+
+            for delta in self.deltas:  # foreach neighborhood size
+                # Smooth the permuted map using delta proportion of
+                # neighbors to reintroduce spatial autocorrelation
+                sm_xperm = self.smooth_map(x=xperm, delta=delta)
+
+                # Calculate empirical variogram of the smoothed permuted map
+                vperm = self.compute_variogram(sm_xperm)
+
+                # Calculate smoothed variogram of the smoothed permuted map
+                smvar_perm = self.smooth_variogram(vperm)
+
+                # Fit linear regression btwn smoothed variograms
+                res[delta] = self.regress(smvar_perm, smvar)
+
+            alphas, betas, residuals = np.array(
+                [res[d] for d in self.deltas]).T
+
+            # Select best-fit model and regression parameters
+            iopt = np.argmin(residuals)
+            dopt = self.deltas[iopt]
+            aopt = alphas[iopt]
+            bopt = betas[iopt]
+
+            # Transform and smooth permuted map using best-fit parameters
+            sm_xperm_best = self.smooth_map(x=xperm, delta=dopt)
+            surr = (np.sqrt(np.abs(bopt)) * sm_xperm_best +
+                    np.sqrt(np.abs(aopt)) * self._rs.randn(nmap))
+            surrs[i] = surr
+
+        if self.resample:  # resample values from empirical map
+            sorted_map = np.sort(x)
+            for i, surr in enumerate(surrs):
+                ii = np.argsort(surr)
+                np.put(surr, ii, sorted_map)
+
+        return surrs.squeeze()
+
+    def compute_variogram(self, x):
+        """
+        Compute variogram values (i.e., one-half squared pairwise differences).
+
+        Parameters
+        ----------
+        x : 1D ndarray
+            Brain map scalar array
+
+        Returns
+        -------
+        v : ndarray, shape (N(N-1)/2,)
+           Variogram y-coordinates, i.e. 0.5 * (x_i - x_j) ^ 2
+
+        """
+        diff_ij = np.subtract.outer(x, x)
+        v = 0.5 * np.square(diff_ij)[self._triu]
+        return v
+
+    def permute_map(self, x):
+        """
+        Return randomly permuted brain map.
+
+        Parameters
+        ----------
+        x : 1D masked ndarray
+            Brain map scalars
+
+        Returns
+        -------
+        1D ndarray
+            Random permutation of target brain map
+
+        """
+
+        pidx = self._rs.permutation(x.size)
+        return np.ma.masked_array(data=x.data[pidx], mask=x.mask[pidx])
+
+    def smooth_map(self, x, delta):
+        """
+        Smooth `x` using `delta` proportion of nearest neighbors.
+
+        Parameters
+        ----------
+        x : 1D ndarray
+            Brain map scalars
+        delta : float
+            Proportion of neighbors to include for smoothing, in (0, 1)
+
+        Returns
+        -------
+        1D ndarray
+            Smoothed brain map
+
+        """
+        # Values of k nearest neighbors for each brain area
+        xkn = x[self._jkn[delta]]
+        weights = self.kernel(self._dkn[delta])  # Distance-weight kernel
+        return (weights * xkn).sum(axis=1) / weights.sum(axis=1)
+
+    def smooth_variogram(self, v, return_h=False):
+        """
+        Smooth a variogram.
+
+        Parameters
+        ----------
+        v : 1D ndarray
+            Variogram values, i.e. 0.5 * (x_i - x_j) ^ 2
+        return_h : bool, default False
+            Return distances at which the smoothed variogram was computed
+
+        Returns
+        -------
+        1D ndarray, shape (nh,)
+            Smoothed variogram values
+        1D ndarray, shape (nh,)
+            Distances at which smoothed variogram was computed (returned only
+            if `return_h` is True)
+
+        Raises
+        ------
+        ValueError : `v` has unexpected size.
+
+        """
+        u = self._u[self._uidx]
+        v = v[self._uidx]
+        if len(u) != len(v):
+            raise ValueError(
+                "argument v: expected size {}, got {}".format(len(u), len(v)))
+        # Subtract each h from each pairwise distance u
+        # Each row corresponds to a unique h
+        du = np.abs(u - self._h[:, None])
+        w = np.exp(-np.square(2.68 * du / self.b) / 2)
+        denom = w.sum(axis=1)
+        wv = w * v[None, :]
+        num = wv.sum(axis=1)
+        output = num / denom
+        if not return_h:
+            return output
+        return output, self._h
+
+    def regress(self, x, y):
+        """
+        Linearly regress `x` onto `y`.
+
+        Parameters
+        ----------
+        x : 1D ndarray
+            Independent variable
+        y : 1D ndarray
+            Dependent variable
+
+        Returns
+        -------
+        alpha : float
+            Intercept term (offset parameter)
+        beta : float
+            Regression coefficient (scale parameter)
+        res : float
+            Sum of squared residuals
+
+        """
+        self._lm.fit(X=np.expand_dims(x, -1), y=y)
+        beta = self._lm.coef_.item()
+        alpha = self._lm.intercept_
+        y_pred = self._lm.predict(X=np.expand_dims(x, -1))
+        res = np.sum(np.square(y-y_pred))
+        return alpha, beta, res
+
+    @property
+    def h(self):
+        """ 1D ndarray : distances at which smoothed variogram is computed """
+        return self._h
+
+
+class SampledSurrogateMaps(BaseEstimator):
+    """
+    Spatial autocorrelation-preserving surrogate brain maps wih sampling.
+
+    Parameters
+    ----------
+    ns : int, optional
+        Take a subsample of `ns` rows from `D` when fitting variograms.
+        Default is 500.
+    deltas : 1D ndarray or List[float], optional
+        Proportion of neighbors to include for smoothing, in (0, 1]
+        Default is [0.1,0.2,...,0.9].
+    kernel : str, optional
+        Kernel with which to smooth permuted maps:
+          'gaussian' : Gaussian function.
+          'exp' : Exponential decay function.
+          'invdist' : Inverse distance.
+          'uniform' : Uniform weights (distance independent).
+        Default is 'exp'.
+    pv : int, optional
+        Percentile of the pairwise distance distribution at which to
+        truncate during variogram fitting. Default is 25.
+    nh : int, optional
+        Number of uniformly spaced distances at which to compute variogram.
+        Default is 25.
+    knn : int, optional
+        Number of nearest regions to keep in the neighborhood of each region.
+        Default is 1000.
+    b : float or None, default None
+        Gaussian kernel bandwidth for variogram smoothing. if None,
+        three times the distance interval spacing is used.
+    resample : bool, optional
+        Resample surrogate maps' values from target brain map.
+        Default is False.
+    n_rep : int, optional
+        Number of randomizations (i.e., surrogate maps). Default is 100.
+    random_state : int or None, optional
+        Random state. Default is None.
+    verbose : bool, default False
+        Print surrogate count each time new surrogate map created
+
+    See Also
+    --------
+    :class:`.SurrogateMaps`
+    :class:`.SpinPermutations`
+    :class:`.MoranRandomization`
+
+    Notes
+    -----
+    Passing resample=True will preserve the distribution of values in the
+    target map, at the expense of worsening simulated surrogate maps'
+    variograms fits. This worsening will increase as the empirical map
+    more strongly deviates from normality.
+
+    """
+
+    def __init__(self, ns=500, pv=70, nh=25, knn=1000, b=None, deltas=None,
+                 kernel='exp', resample=False, n_rep=100, random_state=None,
+                 verbose=False):
+
+        self.deltas = np.arange(0.3, 1., 0.2) if deltas is None else deltas
+        check_deltas(deltas=self.deltas)
+
+        self.ns = int(ns)
+        self.resample = resample
+        self.nh = nh
+        self.pv = check_pv(pv)
+        self.kernel = check_kernel(kernel)  # Smoothing kernel selection
+        self.knn = knn
+        self.b = b
+        self.verbose = verbose
+        self.n_rep = n_rep
+        self.random_state = random_state
+
+        self._rs = np.random.RandomState(self.random_state)
+
+        # Linear regression model
+        self._lm = LinearRegression(fit_intercept=True)
+
+    def _check_distance(self, dist):
+        dist = dataio(dist)
+        return dist[:, 1:self.knn + 1]  # prevent self-coupling
+
+    def _check_index(self, index):
+        index = dataio(index)
+        return index[:, 1:self.knn + 1].astype(np.int32)
+
+    def _check_map(self, x):
+        self._ismasked = False
+        x = dataio(x)
+        check_map(x=x)
+        mask = np.isnan(x)
+        if mask.any():
+            self._ismasked = True
+            return np.ma.masked_array(data=x, mask=mask)
+        return x
+
+    def fit(self, dist, index):
+        """ Prepare data for surrogate map generation.
+
+        Parameters
+        ----------
+        dist : ndarray or memmap, shape (N,N)
+            Pairwise distance matrix. Each row of `dist` should be sorted.
+            Indices used to sort each row are passed to through the `index`
+            argument. See :func:`brainspace.variogram.txt2memmap`.
+        index : filename or ndarray or memmap, shape(N,N)
+            Indices used to sort each row of `dist`.
+
+        Returns
+        -------
+        self : object
+            Returns self.
+
+        """
+
+        self._dist = self._check_distance(dist)
+        self._index = self._check_index(index)
+        if self._dist.shape != self._index.shape:
+            raise ValueError('Dimensions of distanace and index matrices are '
+                             'inconsistent: distance shape = {}, index '
+                             'shape = {}.'.format(self._dist.shape,
+                                                  self._index.shape))
+
+        n = self._dist.shape[0]
+        self._ikn = np.arange(n)[:, None]
+        self._dmax = np.percentile(self._dist, self.pv)
+        self._h = np.linspace(self._dist.min(), self._dmax, self.nh)
+        if not self.b:
+            self.b = 3 * (self._h[1] - self._h[0])
+
+        return self
+
+    def randomize(self, x, n_rep=None):
+        """ Generate surrogate maps from `x`.
+
+        Parameters
+        ----------
+        x : filename or 1D ndarray
+            Target brain map
+        n_rep : int or None, optional
+            Number of surrogates maps to randomly generate. If None, use
+            the default `n_rep`.
+
+        Returns
+        -------
+        output : ndarray, shape = (n_rep, n_verts)
+            Randomly generated map(s) with matched spatial autocorrelation.
+
+        """
+
+        x = self._check_map(x)
+        n = x.size
+
+        if n != self._dist.shape[0]:
+            raise ValueError("Size of distance matrix along axis=0 must equal "
+                             "brain map size.")
+
+        n_rep = self.n_rep if n_rep is None else n_rep
+
+        if self.verbose:
+            print("Generating {} maps...".format(n_rep))
+
+        surrs = np.empty((n_rep, n))
+        for i in range(n_rep):  # generate random maps
+            if self.verbose:
+                print(i+1)
+
+            # Randomly permute map
+            x_perm = self.permute_map(x)
+
+            # Randomly select subset of regions to use for variograms
+            idx = self.sample(n)
+
+            # Compute empirical variogram
+            v = self.compute_variogram(x, idx)
+
+            # Variogram ordinates; use nearest neighbors because local effect
+            u = self._dist[idx, :]
+            uidx = np.where(u < self._dmax)
+
+            # Smooth empirical variogram
+            smvar, u0 = self.smooth_variogram(u[uidx], v[uidx], return_h=True)
+
+            res = dict.fromkeys(self.deltas)
+            for d in self.deltas:  # foreach neighborhood size
+
+                k = int(d * self.knn)
+
+                # Smooth the permuted map using k nearest neighbors to
+                # reintroduce spatial autocorrelation
+                sm_xperm = self.smooth_map(x=x_perm, k=k)
+
+                # Calculate variogram values for the smoothed permuted map
+                vperm = self.compute_variogram(sm_xperm, idx)
+
+                # Calculate smoothed variogram of the smoothed permuted map
+                smvar_perm = self.smooth_variogram(u[uidx], vperm[uidx])
+
+                # Fit linear regression btwn smoothed variograms
+                res[d] = self.regress(smvar_perm, smvar)
+
+            alphas, betas, residuals = np.array(
+                [res[d] for d in self.deltas]).T
+
+            # Select best-fit model and regression parameters
+            iopt = np.argmin(residuals)
+            dopt = self.deltas[iopt]
+            kopt = int(dopt * self.knn)
+            aopt = alphas[iopt]
+            bopt = betas[iopt]
+
+            # Transform and smooth permuted map using best-fit parameters
+            sm_xperm_best = self.smooth_map(x=x_perm, k=kopt)
+            surr = (np.sqrt(np.abs(bopt)) * sm_xperm_best +
+                    np.sqrt(np.abs(aopt)) * self._rs.randn(n))
+            surrs[i] = surr
+
+        if self.resample:  # resample values from empirical map
+            sorted_map = np.sort(x)
+            for i, surr in enumerate(surrs):
+                ii = np.argsort(surr)
+                np.put(surr, ii, sorted_map)
+
+        if self._ismasked:
+            return np.ma.masked_array(
+                data=surrs, mask=np.isnan(surrs)).squeeze()
+        return surrs.squeeze()
+
+    def compute_variogram(self, x, idx):
+        """
+        Compute variogram of `x` using pairs of regions indexed by `idx`.
+
+        Parameters
+        ----------
+        x : 1D ndarray
+            Brain map
+        idx : ndarray[int], shape (ns,)
+            Indices of randomly sampled brain regions
+
+        Returns
+        -------
+        v : ndarray, shape (ns,ns)
+            Variogram y-coordinates, i.e. 0.5 * (x_i - x_j) ^ 2,
+            for i,j in idx
+
+        """
+
+        diff_ij = x[idx][:, None] - x[self._index[idx, :]]
+        return 0.5 * np.square(diff_ij)
+
+    def permute_map(self, x):
+        """
+        Return a random permutation of `x`.
+
+        Parameters
+        ----------
+        x : 1D ndarray
+            Brain map
+
+        Returns
+        -------
+        1D ndarray
+            Random permutation of target brain map
+
+        """
+
+        pidx = self._rs.permutation(x.size)
+        if self._ismasked:
+            return np.ma.masked_array(data=x.data[pidx], mask=x.mask[pidx])
+        return x[pidx]
+
+    def smooth_map(self, x, k):
+        """
+        Smooth `x` using `k` nearest neighboring regions.
+
+        Parameters
+        ----------
+        x : 1D ndarray
+            Brain map
+        k : float
+            Number of nearest neighbors to include for smoothing
+
+        Returns
+        -------
+        x_smooth : 1D ndarray
+            Smoothed brain map
+
+        Notes
+        -----
+        Assumes `dist` provided at runtime has been sorted.
+
+        """
+        jkn = self._index[:, :k]  # indices of k nearest neighbors
+        xkn = x[jkn]  # values of k nearest neighbors
+        dkn = self._dist[:, :k]  # distances to k nearest neighbors
+        weights = self.kernel(dkn)  # distance-weighted kernel
+        return (weights * xkn).sum(axis=1) / weights.sum(axis=1)
+
+    def smooth_variogram(self, u, v, return_h=False):
+        """
+        Smooth a variogram.
+
+        Parameters
+        ----------
+        u : 1D ndarray
+            Pairwise distances, ie variogram x-coordinates
+        v : 1D ndarray
+            Squared differences, ie variogram y-coordinates
+        return_h : bool, default False
+            Return distances at which smoothed variogram is computed
+
+        Returns
+        -------
+        ndarray, shape (nh,)
+            Smoothed variogram samples
+        ndarray, shape (nh,)
+            Distances at which smoothed variogram was computed (returned if
+            `return_h` is True)
+
+        Raises
+        ------
+        ValueError : `u` and `v` are not identically sized
+
+        """
+        if len(u) != len(v):
+            raise ValueError("u and v must have same number of elements")
+
+        # Subtract each element of h from each pairwise distance `u`.
+        # Each row corresponds to a unique h.
+        du = np.abs(u - self._h[:, None])
+        w = np.exp(-np.square(2.68 * du / self.b) / 2)
+        denom = w.sum(axis=1)
+        wv = w * v[None, :]
+        num = wv.sum(axis=1)
+        output = num / denom
+        if not return_h:
+            return output
+        return output, self._h
+
+    def regress(self, x, y):
+        """
+        Linearly regress `x` onto `y`.
+
+        Parameters
+        ----------
+        x : 1D ndarray
+            Independent variable
+        y : 1D ndarray
+            Dependent variable
+
+        Returns
+        -------
+        alpha : float
+            Intercept term (offset parameter)
+        beta : float
+            Regression coefficient (scale parameter)
+        res : float
+            Sum of squared residuals
+
+        """
+        self._lm.fit(X=np.expand_dims(x, -1), y=y)
+        beta = self._lm.coef_.item()
+        alpha = self._lm.intercept_
+        ypred = self._lm.predict(np.expand_dims(x, -1))
+        res = np.sum(np.square(y-ypred))
+        return alpha, beta, res
+
+    def sample(self, n):
+        """
+        Randomly sample (without replacement) brain areas for variogram
+        computation.
+
+        Returns
+        -------
+        ndarray, shape (ns,)
+            Indices of randomly sampled areas
+
+        """
+        return self._rs.choice(a=n, size=self.ns,
+                               replace=False).astype(np.int32)
+
+    @property
+    def h(self):
+        """ 1D ndarray : distances at which variogram is evaluated """
+        return self._h
diff --git a/docs/conf.py b/docs/conf.py
index 9acc3f0..ec5eb29 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -109,6 +109,8 @@
     'image_scrapers': ('matplotlib', _get_sg_image_scraper()),
     'within_subsection_order': FileNameSortKey,
     'download_all_examples': False,
+    'remove_config_comments': True,
+    # 'run_stale_examples': True,
 }
 
 
diff --git a/docs/index.rst b/docs/index.rst
index bcdca89..c30de17 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -6,8 +6,11 @@
 Welcome to BrainSpace's documentation!
 ======================================
 
-.. image:: https://img.shields.io/travis/MICA-MNI/BrainSpace/master?style=flat&logo=travis
-   :target: https://travis-ci.org/MICA-MNI/BrainSpace
+.. image:: https://github.com/MICA-LAB/brainspace/workflows/Python%20package/badge.svg
+   :target: https://github.com/MICA-LAB/brainspace/actions
+
+.. image:: https://codecov.io/gh/mica-mni/brainspace/branch/master/graph/badge.svg
+   :target: https://codecov.io/gh/mica-mni/brainspace
 
 .. image:: https://img.shields.io/appveyor/build/OualidBenkarim/brainspace/master?logo=appveyor
    :alt: AppVeyor branch
diff --git a/docs/pages/install.rst b/docs/pages/install.rst
index 2632f33..d045629 100644
--- a/docs/pages/install.rst
+++ b/docs/pages/install.rst
@@ -65,8 +65,11 @@ This toolbox has been tested with MATLAB versions R2018b, although we expect it
 to work with versions R2018a and newer. Operating systems used during testing were OSX Mojave (10.14.6)
 and Linux Xenial Xerus (16.04.6).
 
-To install the MATLAB toolbox simply `download
-<https://github.com/MICA-MNI/BrainSpace/releases>`_ and unzip the GitHub toolbox and run
+The recommended installation approach is to download the code as a toolbox from FileExchange or Github,
+and open this toolbox with MATLAB. 
+
+If you prefer installing BrainSpace by simply adding it to a path, you can also `download
+<https://github.com/MICA-MNI/BrainSpace/releases>`_ and unzip the code from Github and run
 the following in MATLAB: ::
 
     addpath(genpath('/path/to/BrainSpace/matlab/'))
diff --git a/docs/pages/matlab_doc/data_loaders/load_group_fc.rst b/docs/pages/matlab_doc/data_loaders/load_group_fc.rst
index 6c1b3f9..13646a9 100644
--- a/docs/pages/matlab_doc/data_loaders/load_group_fc.rst
+++ b/docs/pages/matlab_doc/data_loaders/load_group_fc.rst
@@ -26,3 +26,8 @@ Usage
 .. _DK: https://surfer.nmr.mgh.harvard.edu/ftp/articles/desikan06-parcellation.pdf
 
 __ DK_
+
+
+.. note ::
+    Data matrices were drived from the discovery (main) and validation (holdout) groups of `(Vos de Wael et
+    al., 2018) <https://www.pnas.org/content/115/40/10154.short>`_. 
\ No newline at end of file
diff --git a/docs/pages/matlab_doc/data_loaders/load_metric.rst b/docs/pages/matlab_doc/data_loaders/load_metric.rst
index 9c3a56e..05d6a5e 100644
--- a/docs/pages/matlab_doc/data_loaders/load_metric.rst
+++ b/docs/pages/matlab_doc/data_loaders/load_metric.rst
@@ -21,3 +21,7 @@ Usage
 - *name*: The type of surface. Either 'thickness' for cortical thickness, 't1wt2w' for t1w/t2w intensity or 'curvature' for curvature. 
 - *metric_lh*: Data on left hemisphere.
 - *metric_rh*: Data on right hemisphere. 
+
+.. note ::
+    Data matrices were drived from the discovery (main) group of `(Vos de Wael et
+    al., 2018) <https://www.pnas.org/content/115/40/10154.short>`_. 
\ No newline at end of file
diff --git a/docs/pages/matlab_doc/examples/tutorial_1_gradient_creation.rst b/docs/pages/matlab_doc/examples/tutorial_1_gradient_creation.rst
index bdb7920..ccb69c2 100644
--- a/docs/pages/matlab_doc/examples/tutorial_1_gradient_creation.rst
+++ b/docs/pages/matlab_doc/examples/tutorial_1_gradient_creation.rst
@@ -10,8 +10,6 @@ data for computational efficiency.
 
 .. code-block:: matlab
 
-    addpath(genpath('/path/to/BrainSpace/matlab'));
-
     % First load mean connectivity matrix and Schaefer parcellation
     conn_matrix = load_group_fc('schaefer',400);
     labeling = load_parcellation('schaefer',400);
diff --git a/docs/pages/matlab_doc/examples/tutorial_2_gradient_customization.rst b/docs/pages/matlab_doc/examples/tutorial_2_gradient_customization.rst
index e7d2e31..41a6f3f 100644
--- a/docs/pages/matlab_doc/examples/tutorial_2_gradient_customization.rst
+++ b/docs/pages/matlab_doc/examples/tutorial_2_gradient_customization.rst
@@ -14,8 +14,6 @@ Customizing Gradient Computation
 As before, we'll start by loading the sample data.
 
 .. code-block:: matlab    
-
-    addpath(genpath('/path/to/BrainSpace/matlab'));
     
     % First load mean connectivity matrix and Schaefer parcellation
     conn_matrix = load_group_fc('schaefer',400);
@@ -113,7 +111,9 @@ angle kernel and laplacian eigenmap approach.
 
 Here, ``Gp`` contains the Procrustes aligned data and ``Gj`` contains the joint
 aligned data. Let's plot them, but in separate figures to keep things organized.
-
+Note that the gradients in the first element of the gradients/aligned property
+correspond to the first data matrix provided in fit() and the gradients in the
+second element correspond to the second data matrix. 
 
 .. code-block:: matlab    
     
@@ -180,8 +180,6 @@ First, let's load the example data of microstructural profile covariance and fun
 
 .. code-block:: matlab
 
-    addpath(genpath('/path/to/BrainSpace/matlab'));
-
     % First load mean connectivity matrix and parcellation
     mpc = load_group_mpc('vosdewael',200);
     fc = load_group_fc('vosdewael',200);
diff --git a/docs/pages/matlab_doc/examples/tutorial_3_null_models.rst b/docs/pages/matlab_doc/examples/tutorial_3_null_models.rst
index 53c93ba..f337c91 100644
--- a/docs/pages/matlab_doc/examples/tutorial_3_null_models.rst
+++ b/docs/pages/matlab_doc/examples/tutorial_3_null_models.rst
@@ -9,7 +9,7 @@ bias the test statistic. In this tutorial we will show two approaches for null
 hypothesis testing: spin permutations and Moran spectral randomization. 
 
 .. note:: 
-    When using either approach to compare gradients to non-gradient markers, we
+    When using these approaches to compare gradients to non-gradient markers, we
     recommend randomizing the non-gradient markers as these randomizations need not
     maintain the statistical independence between gradients.
 
@@ -26,8 +26,6 @@ cortical thickness data, and a template functional gradient.
 
 .. code-block:: matlab
 
-    addpath(genpath('/path/to/BrainSpace/matlab')); 
-
     % load the conte69 hemisphere surfaces and spheres
     [surf_lh, surf_rh] = load_conte69;
     [sphere_lh, sphere_rh] = load_conte69('spheres');
@@ -143,8 +141,6 @@ gradient.
 
 .. code-block:: matlab
 
-    addpath(genpath('/path/to/BrainSpace/matlab')); 
-
     % load the conte69 hemisphere surfaces and spheres
     [surf_lh, surf_rh] = load_conte69();
 
@@ -264,3 +260,80 @@ Indeed, our histograms appear to be normally distributed. This concludes the
 third and last tutorial. You should now be familliar with all the functionality
 of the BrainSpace toolbox. For more details on any specific function, please see
 :ref:`matlab_package`.
+
+
+Variogram Matching
+------------------
+Here, we use the variogram matching presented in `(Burt et al., 2020)
+<https://www.sciencedirect.com/science/article/pii/S1053811920305243>`_, which
+generates novel brainmaps with similar spatial autocorrelation to the input
+data. We will start by loading the conte69 surfaces for left and right
+hemispheres, a left temporal lobe mask, t1w/t2w intensity as well as cortical
+thickness data. 
+
+.. code-block:: matlab
+
+    % load the conte69 hemisphere surfaces and spheres
+    [surf_lh, surf_rh] = load_conte69();
+
+    % Load the data 
+    t1wt2w_lh = load_marker('t1wt2w');
+    curv_lh = load_marker('curvature');
+
+    % Load mask
+    temporal_mask_tmp = load_mask('temporal');
+
+    % There's a one vertex overlap between the HCP midline mask (i.e. nans) and
+    % our temporal mask.
+    temporal_mask_lh = temporal_mask_tmp & ~isnan(t1wt2w_lh);
+
+    % Load the embedding
+    embedding = load_gradient('fc',1);
+    embedding_lh = embedding(1:end/2);
+
+    % Keep only the temporal lobe. 
+    t1wt2w_tl = t1wt2w_lh(temporal_mask_lh);
+    curv_tl = curv_lh(temporal_mask_lh);
+
+Next, we will need a distance matrix that tells us what the spatial distance
+between our datapoints is. For this example, we will use geodesic distance. 
+
+.. code-block :: matlab
+
+    G = surface_to_graph(surf_lh, 'geodesic', ~temporal_mask_lh);
+    geodesic_distance = distances(G);
+
+Now we've got everything we need to generate our surrogate datasets. By default,
+BrainSpace will use all available data to generate surrogate maps. However, this
+process is extremely computationally and memory intensive. When using this method with more
+than a few hundred regions, we recommend subsampling the data. This can be done using
+the 'ns' and 'knn' name-value pairs. 'ns' determines how many data points to sample
+for the generation of the variogram; 'knn' determines how many neighbors to use
+in the smoothing step. 
+
+.. code-block :: matlab
+
+    random_initialization = 0; 
+    n_surrogate_datasets = 10;
+
+    % Note: num_samples must be greater than num_neighbors
+    num_samples = 200;
+    num_neighbors = 100;
+
+    obj_subsample = variogram(geodesic_distance, 'ns', num_samples, ...
+        'knn', num_neighbors, 'random_state', random_initialization);
+    surrogates_subsample = obj_subsample.fit(t1wt2w_lh, n_surrogate_datasets);
+
+.. note:: 
+    The variogram class also supports parallel processing with the 'num_workers'
+    Name-Value pair. This functionality requires the Parallel Computing Toolbox.
+
+The correlation between the real data, surrogate data, and the marker of interest 
+can then be compared as follows. Note that, for this example, we only generated 
+10 surogate datasets. For any practical usage we recommend generating at least 1000.
+
+.. code-block :: matlab
+
+    r_real = corr(t1wt2w_tl, curv_tl);
+    r_surrogate = corr(surrogates_subsample, curv_tl);
+    prctile_rank = mean(r_real > r_surrogate); 
\ No newline at end of file
diff --git a/docs/pages/matlab_doc/main_functionality/gradientmaps.rst b/docs/pages/matlab_doc/main_functionality/gradientmaps.rst
index 3e0910f..7e49f69 100644
--- a/docs/pages/matlab_doc/main_functionality/gradientmaps.rst
+++ b/docs/pages/matlab_doc/main_functionality/gradientmaps.rst
@@ -30,8 +30,7 @@ Properties
 
 The **method** property is a structure array which itself consists of four
 fields. Each of these is set at the initalization of the object (see below) and
-cannot be modifed afterwards. The fields are: "kernel", "approach", "alignment",
-and "n_components". 
+cannot be modifed afterwards. The fields are: "kernel", "approach", "alignment", "n_components", and "verbose". 
 
 The **gradients** property is a cell array containing the (unaligned) gradients
 of each input matrix. Each cell is an n-by-m matrix where n is the number of
@@ -76,9 +75,11 @@ be provided to alter its behavior.
    - Any input accepted by MATLAB's ``rng`` function, or nan to use the current random state. (default: nan)
 'n_components'
    - Any natural number in double format. (default: 10)
+'verbose'
+   - Determines whether non-warning/error messages are diplayed (default: false).
 
 Putting it all together, an example initialization could be: ``gm =
-GradientMaps('kernel','g','approach','pca','alignment','','random_state',10);``
+GradientMaps('kernel','g','approach','pca','alignment','','random_state',10,'verbose',true);``
 
 Public Methods
 ---------------
diff --git a/docs/pages/matlab_doc/main_functionality/moran_randomization.rst b/docs/pages/matlab_doc/main_functionality/moran_randomization.rst
index aa5c4c5..3adce56 100644
--- a/docs/pages/matlab_doc/main_functionality/moran_randomization.rst
+++ b/docs/pages/matlab_doc/main_functionality/moran_randomization.rst
@@ -22,7 +22,7 @@ Usage
 
 - *Y*: An n-by-m data matrix to randomize where n is number of datapoints and m are different modalities. 
 - *MEM*: Moran eigenvectors as returned by :ref:`compute_mem_matlab`.
-- *n_rep*: Number of perutations
+- *n_rep*: Number of pertubations
 - *varargin*: See name-value pairs below. 
 - *Y_rand*: Randomized data. 
 
diff --git a/docs/pages/matlab_doc/main_functionality/variogram.rst b/docs/pages/matlab_doc/main_functionality/variogram.rst
new file mode 100644
index 0000000..fdaf4cb
--- /dev/null
+++ b/docs/pages/matlab_doc/main_functionality/variogram.rst
@@ -0,0 +1,52 @@
+.. _variogram_matlab:
+
+====================
+variogram
+====================
+
+------------------
+Synopsis
+------------------
+
+Generates null model data correcting for spatial autocorrelation. 
+
+------------------
+Usage
+------------------
+
+::
+
+    obj = variogram(D,varargin);
+    surrogates = obj.fit(x,n);
+
+- *D*: A symmetric distance matrix.
+- *n*: Number of surrogate datasets to generate (default: 1000).
+- *varargin*: See name-value pairs below. 
+- *obj*: The variogram object.
+- *surrogates*: Surrogate datasets.
+
+------------------ 
+Description 
+------------------ 
+
+Implementation of the variogram matching procedure as presented by  `(Burt et al., 2020)
+<https://www.sciencedirect.com/science/article/pii/S1053811920305243>`_.
+This class generates surrogate data by permuting the input data, smoothing it with different 
+kernels, and determining which smoothed permuted map best fits the variogram of the 
+empirical data. 
+
+-----------------------
+Name-Value Pairs
+-----------------------
+
+- *deltas*: Proportion of neighbours to include for smoothing (default: 0.1:0.1:0.9).
+- *kernel*: Kernel with which to smooth permuted maps. Valid options are 'gaussian', 'exp' (default), 'invdist', 'uniform', and a function handle.
+- *pv*: Percentile of pairwise distance distribution at which to truncate during variogram fitting (default: 25).
+- *nh*: Number of uniformly spaced distances at which to compute variograms (default: 25).
+- *resample*: Resample surrogate maps' values from target brain map (default: false).
+- *b*:  Gaussian kernel bandwidth for variogram smoothing (default: three times the spacing between variogram x-coordinates).
+- *random_state*: any valid input for the rng() function. Set to nan for no random initialization (default: nan).
+- *ns*: Number of samples to use when subsampling the brainmap. Set to inf to use the entire brainmap (default: inf).
+- *knn*: Number of nearest neighbours to use when smoothing the map. knn must be smaller than ns (default: 1000). 
+- *num_workers*: Number of workers in the parallel pool. Requires the parallel processing toolbox. Set to 0 for no parallelization (default: 0). 
+
diff --git a/docs/pages/matlab_doc/visualization/plot_hemispheres.rst b/docs/pages/matlab_doc/visualization/plot_hemispheres.rst
index cebb6aa..8be4ff2 100644
--- a/docs/pages/matlab_doc/visualization/plot_hemispheres.rst
+++ b/docs/pages/matlab_doc/visualization/plot_hemispheres.rst
@@ -9,7 +9,7 @@ Synopsis
 ------------------
 
 Plots data on the cortical surface (`source code
-<https://github.com/MICA-MNI/BrainSpace/blob/master/matlab/plot_data/plot_hemispheres.m>`_).
+<https://github.com/MICA-MNI/BrainSpace/blob/master/matlab/plot_data/%40plot_hemispheres/plot_hemispheres.m>`_).
 
 
 ------------------
@@ -18,12 +18,12 @@ Usage
 
 ::
 
-   handles = plot_hemispheres(data,surface,varargin);
+   obj = plot_hemispheres(data,surface,varargin);
 
 - *data*: an n-by-m vector of data to plot where n is the number of vertices or parcels, and m the number of markers to plot.
 - *surface*: a surface readable by :ref:`convert_surface_matlab` or a two-element cell array containing left and right hemispheric surfaces in that order. 
 - *varargin*: Name-Value Pairs (see below).
-- *handles*: a structure containing the handles of the graphics objects. 
+- *obj*: an object allowing for further modificatio of the figure (see below). 
 
 ------------------
 Description
@@ -51,5 +51,14 @@ for basic graphic object property modifications.
 
 Name-Value Pairs
 ^^^^^^^^^^^^^^^^^
-- *parcellation*: an k-by-1 vector containing the parcellation scheme, where k is the number of vertices. 
-- *labeltext*: A cell array of m elements containing labels for each column of data. These will be printed next to the hemispheres. 
+- *parcellation*: an k-by-1 vector containing the parcellation scheme, where k is the number of vertices. Default [].
+- *labeltext*: A cell array of m elements containing labels for each column of data. These will be printed next to the hemispheres. Default: [].
+- *views*: A character vector containing the requested view angles. Options are: l(ateral), m(edial), i(nferior), s(uperior), a(nterior), and p(osterior). Default: 'lm'.
+
+Public methods
+^^^^^^^^^^^^^^^
+Public methods can be used with obj.(method) e.g. obj.colorlimits to use the colorlimits method.
+
+- *colorlimits(limits)*: Sets the color limits of each row. Limits must be a 2-element numeric vector or n-by-2 matrix where n is the number of columns in obj.data. The first column of limits denotes the minimum color limit and the second the maximum. When limits is a 2-element vector, then the limits are applied to all axes, with limits(1) as minimum and limits(2) as maximum. 
+- *colorMaps(maps)*: Maps must either be an n-by-3 color map, or a cell array with the same number of elements as columns in obj.data, each containing n-by-3 colormaps.
+- *labels(varargin)*: Modifies the properties of the labeltext. Varargin are valid name-value pairs for MATLAB's text function. e.g. `obj.labels('FontSize',25)`
diff --git a/docs/python_doc/api_doc/brainspace.null_models.rst b/docs/python_doc/api_doc/brainspace.null_models.rst
index 7680817..ce7caff 100644
--- a/docs/python_doc/api_doc/brainspace.null_models.rst
+++ b/docs/python_doc/api_doc/brainspace.null_models.rst
@@ -6,6 +6,7 @@ Null models
 
 - :ref:`pysec-nullmodels-spin`
 - :ref:`pysec-nullmodels-moran`
+- :ref:`pysec-nullmodels-variogram`
 
 
 
@@ -39,3 +40,19 @@ Moran spectral randomization
    MoranRandomization
    compute_mem
    moran_randomization
+
+
+
+.. _pysec-nullmodels-variogram:
+
+Surrogate maps
+-----------------------------
+
+.. currentmodule:: brainspace.null_models.variogram
+
+.. autosummary::
+   :toctree: ../../generated/
+
+
+   SurrogateMaps
+   SampledSurrogateMaps
\ No newline at end of file
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_001.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_001.png
index 120c016..e39d147 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_001.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_001.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_002.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_002.png
index c880f54..1e3d78b 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_002.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial0_002.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_002.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_002.png
index 4fe2eee..a52eeb8 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_002.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_002.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_003.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_003.png
index c750862..a40870c 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_003.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial1_003.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_001.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_001.png
index 91225b8..5686e95 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_001.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_001.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_002.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_002.png
index 28c7b3b..b960d7a 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_002.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_002.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_003.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_003.png
index 4b66f6a..00e38a1 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_003.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_003.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_004.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_004.png
index e997be6..ff1f4f3 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_004.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_004.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_005.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_005.png
index df70981..cdd46a3 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_005.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_005.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_006.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_006.png
index 642a997..f31165d 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_006.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_006.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_007.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_007.png
index 7af6eb9..4d3b098 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_007.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial2_007.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_001.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_001.png
index 19ddcdc..8c79bd2 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_001.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_001.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_002.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_002.png
index 52fd63f..afe96e4 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_002.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_002.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_003.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_003.png
index 10dd0dd..c297c96 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_003.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_003.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_004.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_004.png
index eb14cd8..e6396ae 100644
Binary files a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_004.png and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_004.png differ
diff --git a/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_005.png b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_005.png
new file mode 100644
index 0000000..3d9c419
Binary files /dev/null and b/docs/python_doc/auto_examples/images/sphx_glr_plot_tutorial3_005.png differ
diff --git a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial0_thumb.png b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial0_thumb.png
index 04ec19f..b46de32 100644
Binary files a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial0_thumb.png and b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial0_thumb.png differ
diff --git a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial1_thumb.png b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial1_thumb.png
index d0ffdc7..634030a 100644
Binary files a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial1_thumb.png and b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial1_thumb.png differ
diff --git a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial2_thumb.png b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial2_thumb.png
index a49aa67..35e2f5d 100644
Binary files a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial2_thumb.png and b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial2_thumb.png differ
diff --git a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial3_thumb.png b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial3_thumb.png
index 01638f8..a839196 100644
Binary files a/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial3_thumb.png and b/docs/python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial3_thumb.png differ
diff --git a/docs/python_doc/auto_examples/index.rst b/docs/python_doc/auto_examples/index.rst
index 369f364..e596e35 100644
--- a/docs/python_doc/auto_examples/index.rst
+++ b/docs/python_doc/auto_examples/index.rst
@@ -21,6 +21,7 @@ Tutorials
 .. only:: html
 
  .. figure:: /python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial0_thumb.png
+     :alt: Tutorial 0: Preparing your data for gradient analysis
 
      :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial0.py`
 
@@ -41,6 +42,7 @@ Tutorials
 .. only:: html
 
  .. figure:: /python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial1_thumb.png
+     :alt: Tutorial 1: Building your first gradient
 
      :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial1.py`
 
@@ -61,6 +63,7 @@ Tutorials
 .. only:: html
 
  .. figure:: /python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial2_thumb.png
+     :alt: Tutorial 2: Customizing and aligning gradients
 
      :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial2.py`
 
@@ -81,6 +84,7 @@ Tutorials
 .. only:: html
 
  .. figure:: /python_doc/auto_examples/images/thumb/sphx_glr_plot_tutorial3_thumb.png
+     :alt: Tutorial 3: Null models for gradient significance
 
      :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial3.py`
 
diff --git a/docs/python_doc/auto_examples/plot_tutorial0.ipynb b/docs/python_doc/auto_examples/plot_tutorial0.ipynb
index 815ee4a..ad4db34 100644
--- a/docs/python_doc/auto_examples/plot_tutorial0.ipynb
+++ b/docs/python_doc/auto_examples/plot_tutorial0.ipynb
@@ -15,14 +15,14 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nTutorial 0: Preparing your data for gradient analysis\n=====================================================\nIn this example, we will introduce how to preprocess raw MRI data and how\nto prepare it for subsequent gradient analysis in the next tutorials.\n\nRequirements\n------------\nFor this tutorial, you will need to install the Python package\n`load_confounds <https://github.com/SIMEXP/fmriprep_load_confounds>`_. You can\ndo it using ``pip``::\n\n    pip install load_confounds\n  \n\nPreprocessing\n-------------\nBegin with an MRI dataset that is organized in `BIDS\n<https://bids.neuroimaging.io/>`_ format. We recommend preprocessing your data\nusing `fmriprep <http://fmriprep.readthedocs.io/>`_, as described below, but\nany preprocessing pipeline will work.\n\nFollowing is example code to run `fmriprep <http://fmriprep.readthedocs.io/>`_\nusing docker from the command line::\n\n    docker run -ti --rm \\\n      -v <local_BIDS_data_dir>:/data:ro \\\n      -v <local_output_dir>:/out poldracklab/fmriprep:latest \\\n      --output-spaces fsaverage5 \\\n      --fs-license-file license.txt \\\n      /data /out participant\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>For this tutorial, it is crucial to output the data onto a cortical surface\n    template space.</p></div>\n"
+        "\n# Tutorial 0: Preparing your data for gradient analysis\nIn this example, we will introduce how to preprocess raw MRI data and how\nto prepare it for subsequent gradient analysis in the next tutorials.\n\n## Requirements\nFor this tutorial, you will need to install the Python package\n`load_confounds <https://github.com/SIMEXP/fmriprep_load_confounds>`_. You can\ndo it using ``pip``::\n\n    pip install load_confounds\n  \n\n## Preprocessing\nBegin with an MRI dataset that is organized in `BIDS\n<https://bids.neuroimaging.io/>`_ format. We recommend preprocessing your data\nusing `fmriprep <http://fmriprep.readthedocs.io/>`_, as described below, but\nany preprocessing pipeline will work.\n\nFollowing is example code to run `fmriprep <http://fmriprep.readthedocs.io/>`_\nusing docker from the command line::\n\n    docker run -ti --rm \\\n      -v <local_BIDS_data_dir>:/data:ro \\\n      -v <local_output_dir>:/out poldracklab/fmriprep:latest \\\n      --output-spaces fsaverage5 \\\n      --fs-license-file license.txt \\\n      /data /out participant\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>For this tutorial, it is crucial to output the data onto a cortical surface\n    template space.</p></div>\n"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Import the dataset as timeseries\n++++++++++++++++++++++++++++++++\nThe timeseries should be a numpy array with the dimensions: nodes x timepoints  \n\nFollowing is an example for reading in data::  \n\n   import nibabel as nib\n   import numpy as np\n\n   filename = 'filename.{}.mgz' # where {} will be replaced with 'lh' and 'rh'\n   timeseries = [None] * 2\n   for i, h in enumerate(['lh', 'rh']):\n       timeseries[i] = nib.load(filename.format(h)).get_fdata().squeeze()\n   timeseries = np.vstack(timeseries)\n\n"
+        "### Import the dataset as timeseries\nThe timeseries should be a numpy array with the dimensions: nodes x timepoints  \n\nFollowing is an example for reading in data::  \n\n   import nibabel as nib\n   import numpy as np\n\n   filename = 'filename.{}.mgz' # where {} will be replaced with 'lh' and 'rh'\n   timeseries = [None] * 2\n   for i, h in enumerate(['lh', 'rh']):\n       timeseries[i] = nib.load(filename.format(h)).get_fdata().squeeze()\n   timeseries = np.vstack(timeseries)\n\n"
       ]
     },
     {
@@ -47,7 +47,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Confound regression\n++++++++++++++++++++++++\nTo remove confound regressors from the output of the fmriprep pipeline, first\nextract the confound columns. For example::\n\n   import load_confounds\n   confounds_out = load_confounds(\"path to confound file\",\n                              strategy='minimal',\n                              n_components=0.95,\n                              motion_model='6params')\n\n"
+        "### Confound regression\nTo remove confound regressors from the output of the fmriprep pipeline, first\nextract the confound columns. For example::\n\n   import load_confounds\n   confounds_out = load_confounds(\"path to confound file\",\n                              strategy='minimal',\n                              n_components=0.95,\n                              motion_model='6params')\n\n"
       ]
     },
     {
@@ -108,7 +108,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Calculate functional connectivity matrix\n++++++++++++++++++++++++++++++++++++++++\nThe following example uses\n`nilearn <https://nilearn.github.io/auto_examples/03_connectivity/plot_\nsignal_extraction.html#compute-and-display-a-correlation-matrix/>`_:\n\n"
+        "### Calculate functional connectivity matrix\nThe following example uses\n`nilearn <https://nilearn.github.io/auto_examples/03_connectivity/plot_\nsignal_extraction.html#compute-and-display-a-correlation-matrix/>`_:\n\n"
       ]
     },
     {
@@ -144,7 +144,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Run gradient analysis and visualize\n+++++++++++++++++++++++++++++++++++\n\nRun gradient analysis\n\n"
+        "### Run gradient analysis and visualize\n\nRun gradient analysis\n\n"
       ]
     },
     {
@@ -173,7 +173,7 @@
       },
       "outputs": [],
       "source": [
-        "from brainspace.datasets import load_fsa5\nfrom brainspace.plotting import plot_hemispheres\nfrom brainspace.utils.parcellation import map_to_labels\n\n# Map gradients to original parcels\ngrad = [None] * 2\nfor i, g in enumerate(gm.gradients_.T):\n    grad[i] = map_to_labels(g, labeling, mask=mask, fill=np.nan)\n\n\n# Load fsaverage5 surfaces\nsurf_lh, surf_rh = load_fsa5()\n\n# sphinx_gallery_thumbnail_number = 2\nplot_hemispheres(surf_lh, surf_rh, array_name=grad, size=(1200, 400), cmap='viridis_r',\n                 color_bar=True, label_text=['Grad1', 'Grad2'], zoom=1.5)"
+        "from brainspace.datasets import load_fsa5\nfrom brainspace.plotting import plot_hemispheres\nfrom brainspace.utils.parcellation import map_to_labels\n\n# Map gradients to original parcels\ngrad = [None] * 2\nfor i, g in enumerate(gm.gradients_.T):\n    grad[i] = map_to_labels(g, labeling, mask=mask, fill=np.nan)\n\n\n# Load fsaverage5 surfaces\nsurf_lh, surf_rh = load_fsa5()\n\nplot_hemispheres(surf_lh, surf_rh, array_name=grad, size=(1200, 400), cmap='viridis_r',\n                 color_bar=True, label_text=['Grad1', 'Grad2'], zoom=1.5)"
       ]
     },
     {
@@ -200,7 +200,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.6"
+      "version": "3.7.10"
     }
   },
   "nbformat": 4,
diff --git a/docs/python_doc/auto_examples/plot_tutorial0.rst b/docs/python_doc/auto_examples/plot_tutorial0.rst
index 888701c..acbf985 100644
--- a/docs/python_doc/auto_examples/plot_tutorial0.rst
+++ b/docs/python_doc/auto_examples/plot_tutorial0.rst
@@ -1,7 +1,18 @@
-.. note::
-    :class: sphx-glr-download-link-note
 
-    Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial0.py>` to download the full example code
+.. DO NOT EDIT.
+.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
+.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
+.. "python_doc/auto_examples/plot_tutorial0.py"
+.. LINE NUMBERS ARE GIVEN BELOW.
+
+.. only:: html
+
+    .. note::
+        :class: sphx-glr-download-link-note
+
+        Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial0.py>`
+        to download the full example code
+
 .. rst-class:: sphx-glr-example-title
 
 .. _sphx_glr_python_doc_auto_examples_plot_tutorial0.py:
@@ -42,6 +53,8 @@ using docker from the command line::
     For this tutorial, it is crucial to output the data onto a cortical surface
     template space.
 
+.. GENERATED FROM PYTHON SOURCE LINES 40-54
+
 Import the dataset as timeseries
 ++++++++++++++++++++++++++++++++
 The timeseries should be a numpy array with the dimensions: nodes x timepoints  
@@ -57,8 +70,11 @@ Following is an example for reading in data::
        timeseries[i] = nib.load(filename.format(h)).get_fdata().squeeze()
    timeseries = np.vstack(timeseries)
 
+.. GENERATED FROM PYTHON SOURCE LINES 57-58
+
 As a **working example**, simply fetch timeseries:
 
+.. GENERATED FROM PYTHON SOURCE LINES 58-62
 
 .. code-block:: default
 
@@ -73,6 +89,8 @@ As a **working example**, simply fetch timeseries:
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 63-73
+
 Confound regression
 ++++++++++++++++++++++++
 To remove confound regressors from the output of the fmriprep pipeline, first
@@ -84,8 +102,11 @@ extract the confound columns. For example::
                               n_components=0.95,
                               motion_model='6params')
 
+.. GENERATED FROM PYTHON SOURCE LINES 76-77
+
 As a **working example**, simply read in confounds
 
+.. GENERATED FROM PYTHON SOURCE LINES 77-81
 
 .. code-block:: default
 
@@ -100,8 +121,11 @@ As a **working example**, simply read in confounds
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 82-83
+
 Do the confound regression
 
+.. GENERATED FROM PYTHON SOURCE LINES 83-88
 
 .. code-block:: default
 
@@ -117,8 +141,11 @@ Do the confound regression
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 89-90
+
 And extract the cleaned timeseries onto a set of labels
 
+.. GENERATED FROM PYTHON SOURCE LINES 90-117
 
 .. code-block:: default
 
@@ -153,15 +180,34 @@ And extract the cleaned timeseries onto a set of labels
 
 
 
+.. rst-class:: sphx-glr-script-out
+
+ Out:
+
+ .. code-block:: none
+
+    /home/oualid/Apps/anaconda3/envs/py3/lib/python3.7/site-packages/nilearn/datasets/__init__.py:89: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
+      "Numpy arrays.", FutureWarning)
+
+    Dataset created in /home/oualid/nilearn_data/destrieux_surface
+
+    Downloading data from https://www.nitrc.org/frs/download.php/9343/lh.aparc.a2009s.annot ...
+     ...done. (0 seconds, 0 min)
+    Downloading data from https://www.nitrc.org/frs/download.php/9342/rh.aparc.a2009s.annot ...
+     ...done. (0 seconds, 0 min)
 
 
 
+
+.. GENERATED FROM PYTHON SOURCE LINES 118-123
+
 Calculate functional connectivity matrix
 ++++++++++++++++++++++++++++++++++++++++
 The following example uses
 `nilearn <https://nilearn.github.io/auto_examples/03_connectivity/plot_
 signal_extraction.html#compute-and-display-a-correlation-matrix/>`_:
 
+.. GENERATED FROM PYTHON SOURCE LINES 123-130
 
 .. code-block:: default
 
@@ -179,8 +225,11 @@ signal_extraction.html#compute-and-display-a-correlation-matrix/>`_:
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 131-132
+
 Plot the correlation matrix:
 
+.. GENERATED FROM PYTHON SOURCE LINES 132-148
 
 .. code-block:: default
 
@@ -204,17 +253,21 @@ Plot the correlation matrix:
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial0_001.png
+    :alt: plot tutorial0
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 149-153
+
 Run gradient analysis and visualize
 +++++++++++++++++++++++++++++++++++
 
 Run gradient analysis
 
+.. GENERATED FROM PYTHON SOURCE LINES 153-160
 
 .. code-block:: default
 
@@ -236,13 +289,15 @@ Run gradient analysis
  .. code-block:: none
 
 
-    GradientMaps(alignment=None, approach='dm', kernel=None, n_components=2,
-                 random_state=0)
+    GradientMaps(n_components=2, random_state=0)
+
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 161-162
 
 Visualize results
 
+.. GENERATED FROM PYTHON SOURCE LINES 162-179
 
 .. code-block:: default
 
@@ -259,7 +314,6 @@ Visualize results
     # Load fsaverage5 surfaces
     surf_lh, surf_rh = load_fsa5()
 
-    # sphinx_gallery_thumbnail_number = 2
     plot_hemispheres(surf_lh, surf_rh, array_name=grad, size=(1200, 400), cmap='viridis_r',
                      color_bar=True, label_text=['Grad1', 'Grad2'], zoom=1.5)
 
@@ -268,19 +322,22 @@ Visualize results
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial0_002.png
+    :alt: plot tutorial0
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 181-183
+
 This concludes the setup tutorial. The following tutorials can be run using
 either the output generated here or the example data.
 
 
 .. rst-class:: sphx-glr-timing
 
-   **Total running time of the script:** ( 0 minutes  4.417 seconds)
+   **Total running time of the script:** ( 0 minutes  3.871 seconds)
 
 
 .. _sphx_glr_download_python_doc_auto_examples_plot_tutorial0.py:
@@ -293,13 +350,13 @@ either the output generated here or the example data.
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-python
 
      :download:`Download Python source code: plot_tutorial0.py <plot_tutorial0.py>`
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-jupyter
 
      :download:`Download Jupyter notebook: plot_tutorial0.ipynb <plot_tutorial0.ipynb>`
 
diff --git a/docs/python_doc/auto_examples/plot_tutorial1.ipynb b/docs/python_doc/auto_examples/plot_tutorial1.ipynb
index 42c6b30..65efd47 100644
--- a/docs/python_doc/auto_examples/plot_tutorial1.ipynb
+++ b/docs/python_doc/auto_examples/plot_tutorial1.ipynb
@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nTutorial 1: Building your first gradient\n=================================================\nIn this example, we will derive a gradient and do some basic inspections to\ndetermine which gradients may be of interest and what the multidimensional\norganization of the gradients looks like.\n"
+        "\n# Tutorial 1: Building your first gradient\nIn this example, we will derive a gradient and do some basic inspections to\ndetermine which gradients may be of interest and what the multidimensional\norganization of the gradients looks like.\n"
       ]
     },
     {
@@ -132,7 +132,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.6"
+      "version": "3.7.10"
     }
   },
   "nbformat": 4,
diff --git a/docs/python_doc/auto_examples/plot_tutorial1.rst b/docs/python_doc/auto_examples/plot_tutorial1.rst
index a6b14a5..e65ac72 100644
--- a/docs/python_doc/auto_examples/plot_tutorial1.rst
+++ b/docs/python_doc/auto_examples/plot_tutorial1.rst
@@ -1,7 +1,18 @@
-.. note::
-    :class: sphx-glr-download-link-note
 
-    Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial1.py>` to download the full example code
+.. DO NOT EDIT.
+.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
+.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
+.. "python_doc/auto_examples/plot_tutorial1.py"
+.. LINE NUMBERS ARE GIVEN BELOW.
+
+.. only:: html
+
+    .. note::
+        :class: sphx-glr-download-link-note
+
+        Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial1.py>`
+        to download the full example code
+
 .. rst-class:: sphx-glr-example-title
 
 .. _sphx_glr_python_doc_auto_examples_plot_tutorial1.py:
@@ -13,9 +24,12 @@ In this example, we will derive a gradient and do some basic inspections to
 determine which gradients may be of interest and what the multidimensional
 organization of the gradients looks like.
 
+.. GENERATED FROM PYTHON SOURCE LINES 11-13
+
 We’ll first start by loading some sample data. Note that we’re using
 parcellated data for computational efficiency.
 
+.. GENERATED FROM PYTHON SOURCE LINES 13-25
 
 .. code-block:: default
 
@@ -38,8 +52,11 @@ parcellated data for computational efficiency.
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 26-27
+
 Let’s first look at the parcellation scheme we’re using.
 
+.. GENERATED FROM PYTHON SOURCE LINES 27-34
 
 .. code-block:: default
 
@@ -54,14 +71,18 @@ Let’s first look at the parcellation scheme we’re using.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial1_001.png
+    :alt: plot tutorial1
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 35-36
+
 and let’s construct our gradients.
 
+.. GENERATED FROM PYTHON SOURCE LINES 36-44
 
 .. code-block:: default
 
@@ -84,16 +105,18 @@ and let’s construct our gradients.
  .. code-block:: none
 
 
-    GradientMaps(alignment=None, approach='dm', kernel=None, n_components=10,
-                 random_state=0)
+    GradientMaps(random_state=0)
+
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 45-49
 
 Note that the default parameters are diffusion embedding approach, 10
 components, and no kernel (use raw data). Once you have your gradients, a
 good first step is to simply inspect what they look like. Let’s have a look
 at the first two gradients.
 
+.. GENERATED FROM PYTHON SOURCE LINES 49-65
 
 .. code-block:: default
 
@@ -117,12 +140,15 @@ at the first two gradients.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial1_002.png
+    :alt: plot tutorial1
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 66-72
+
 But which gradients should you keep for your analysis? In some cases you may
 have an a priori interest in some previously defined set of gradients. When
 you do not have a pre-defined set, you can instead look at the lambdas
@@ -130,6 +156,7 @@ you do not have a pre-defined set, you can instead look at the lambdas
 in Laplacian eigenmaps) are more important, so one can choose a cut-off based
 on a scree plot.
 
+.. GENERATED FROM PYTHON SOURCE LINES 72-81
 
 .. code-block:: default
 
@@ -146,12 +173,15 @@ on a scree plot.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial1_003.png
+    :alt: plot tutorial1
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 82-85
+
 This concludes the first tutorial. In the next tutorial we will have a look
 at how to customize the methods of gradient estimation, as well as gradient
 alignments.
@@ -159,7 +189,7 @@ alignments.
 
 .. rst-class:: sphx-glr-timing
 
-   **Total running time of the script:** ( 0 minutes  3.136 seconds)
+   **Total running time of the script:** ( 0 minutes  1.397 seconds)
 
 
 .. _sphx_glr_download_python_doc_auto_examples_plot_tutorial1.py:
@@ -172,13 +202,13 @@ alignments.
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-python
 
      :download:`Download Python source code: plot_tutorial1.py <plot_tutorial1.py>`
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-jupyter
 
      :download:`Download Jupyter notebook: plot_tutorial1.ipynb <plot_tutorial1.ipynb>`
 
diff --git a/docs/python_doc/auto_examples/plot_tutorial2.ipynb b/docs/python_doc/auto_examples/plot_tutorial2.ipynb
index 459c020..26fde5d 100644
--- a/docs/python_doc/auto_examples/plot_tutorial2.ipynb
+++ b/docs/python_doc/auto_examples/plot_tutorial2.ipynb
@@ -15,14 +15,14 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nTutorial 2: Customizing and aligning gradients\n=================================================\nIn this tutorial you\u2019ll learn about the methods available within the\nGradientMaps class. The flexible usage of this class allows for the\ncustomization of gradient computation with different kernels and dimensionality\nreductions, as well as aligning gradients from different datasets. This\ntutorial will only show you how to apply these techniques.\n"
+        "\n# Tutorial 2: Customizing and aligning gradients\nIn this tutorial you\u2019ll learn about the methods available within the\nGradientMaps class. The flexible usage of this class allows for the\ncustomization of gradient computation with different kernels and dimensionality\nreductions, as well as aligning gradients from different datasets. This\ntutorial will only show you how to apply these techniques.\n"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Customizing gradient computation\n+++++++++++++++++++++++++++++++++\nAs before, we\u2019ll start by loading the sample data.\n\n"
+        "## Customizing gradient computation\nAs before, we\u2019ll start by loading the sample data.\n\n"
       ]
     },
     {
@@ -69,14 +69,14 @@
       },
       "outputs": [],
       "source": [
-        "# PCA, Laplacian eigenmaps and diffusion mapping\nembeddings = ['pca', 'le', 'dm']\n\ngradients_embedding = [None] * len(embeddings)\nfor i, emb in enumerate(embeddings):\n    gm = GradientMaps(kernel='normalized_angle', approach=emb, random_state=0)\n    gm.fit(conn_matrix)\n\n    gradients_embedding[i] = map_to_labels(gm.gradients_[:, 0], labeling, mask=mask,\n                                           fill=np.nan)\n\n\n# sphinx_gallery_thumbnail_number = 2\nlabel_text = ['PCA', 'LE', 'DM']\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_embedding, size=(1200, 600),\n                 cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.45)"
+        "# PCA, Laplacian eigenmaps and diffusion mapping\nembeddings = ['pca', 'le', 'dm']\n\ngradients_embedding = [None] * len(embeddings)\nfor i, emb in enumerate(embeddings):\n    gm = GradientMaps(kernel='normalized_angle', approach=emb, random_state=0)\n    gm.fit(conn_matrix)\n\n    gradients_embedding[i] = map_to_labels(gm.gradients_[:, 0], labeling, mask=mask,\n                                           fill=np.nan)\n\n\nlabel_text = ['PCA', 'LE', 'DM']\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_embedding, size=(1200, 600),\n                 cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.45)"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Gradient alignment\n+++++++++++++++++++\n\nA more principled way of increasing comparability across gradients are\nalignment techniques. BrainSpace provides two alignment techniques:\nProcrustes analysis, and joint alignment. For this example we will load\nfunctional connectivity data of a second subject group and align it with the\nfirst group.\n\n"
+        "## Gradient alignment\n\nA more principled way of increasing comparability across gradients are\nalignment techniques. BrainSpace provides two alignment techniques:\nProcrustes analysis, and joint alignment. For this example we will load\nfunctional connectivity data of a second subject group and align it with the\nfirst group.\n\n"
       ]
     },
     {
@@ -152,7 +152,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "The gradients in `galign.aligned_` are now aligned to the reference\ngradients.\n\nGradient fusion\n+++++++++++++++++++\nWe can also fuse data across multiple modalities and build mutli-modal\ngradients. In this case we only look at one set of output gradients,\nrather than one per modality.\n\nFirst, let's load the example data of microstructural profile covariance\n`(Paquola et al., 2019) <https://journals.plos.org/plosbiology/article?\nid=10.1371/journal.pbio.3000284>`_ and functional connectivity.\n\n"
+        "The gradients in `galign.aligned_` are now aligned to the reference\ngradients.\n\n## Gradient fusion\nWe can also fuse data across multiple modalities and build mutli-modal\ngradients. In this case we only look at one set of output gradients,\nrather than one per modality.\n\nFirst, let's load the example data of microstructural profile covariance\n`(Paquola et al., 2019) <https://journals.plos.org/plosbiology/article?\nid=10.1371/journal.pbio.3000284>`_ and functional connectivity.\n\n"
       ]
     },
     {
@@ -226,7 +226,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.6"
+      "version": "3.7.10"
     }
   },
   "nbformat": 4,
diff --git a/docs/python_doc/auto_examples/plot_tutorial2.rst b/docs/python_doc/auto_examples/plot_tutorial2.rst
index 4062630..699c884 100644
--- a/docs/python_doc/auto_examples/plot_tutorial2.rst
+++ b/docs/python_doc/auto_examples/plot_tutorial2.rst
@@ -1,7 +1,18 @@
-.. note::
-    :class: sphx-glr-download-link-note
 
-    Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial2.py>` to download the full example code
+.. DO NOT EDIT.
+.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
+.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
+.. "python_doc/auto_examples/plot_tutorial2.py"
+.. LINE NUMBERS ARE GIVEN BELOW.
+
+.. only:: html
+
+    .. note::
+        :class: sphx-glr-download-link-note
+
+        Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial2.py>`
+        to download the full example code
+
 .. rst-class:: sphx-glr-example-title
 
 .. _sphx_glr_python_doc_auto_examples_plot_tutorial2.py:
@@ -15,10 +26,13 @@ customization of gradient computation with different kernels and dimensionality
 reductions, as well as aligning gradients from different datasets. This
 tutorial will only show you how to apply these techniques.
 
+.. GENERATED FROM PYTHON SOURCE LINES 13-16
+
 Customizing gradient computation
 +++++++++++++++++++++++++++++++++
 As before, we’ll start by loading the sample data.
 
+.. GENERATED FROM PYTHON SOURCE LINES 16-30
 
 .. code-block:: default
 
@@ -43,9 +57,12 @@ As before, we’ll start by loading the sample data.
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 31-33
+
 The GradientMaps object allows for many different kernels and dimensionality
 reduction techniques. Let’s have a look at three different kernels.
 
+.. GENERATED FROM PYTHON SOURCE LINES 33-56
 
 .. code-block:: default
 
@@ -76,18 +93,22 @@ reduction techniques. Let’s have a look at three different kernels.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial2_001.png
+    :alt: plot tutorial2
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 57-62
+
 It seems the gradients provided by these kernels are quite similar although
 their scaling is quite different. Do note that the gradients are in arbitrary
 units, so the smaller/larger axes across kernels do not imply anything.
 Similar to using different kernels, we can also use different dimensionality
 reduction techniques.
 
+.. GENERATED FROM PYTHON SOURCE LINES 62-80
 
 .. code-block:: default
 
@@ -104,7 +125,6 @@ reduction techniques.
                                                fill=np.nan)
 
 
-    # sphinx_gallery_thumbnail_number = 2
     label_text = ['PCA', 'LE', 'DM']
     plot_hemispheres(surf_lh, surf_rh, array_name=gradients_embedding, size=(1200, 600),
                      cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.45)
@@ -114,12 +134,15 @@ reduction techniques.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial2_002.png
+    :alt: plot tutorial2
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 82-90
+
 Gradient alignment
 +++++++++++++++++++
 
@@ -129,6 +152,7 @@ Procrustes analysis, and joint alignment. For this example we will load
 functional connectivity data of a second subject group and align it with the
 first group.
 
+.. GENERATED FROM PYTHON SOURCE LINES 90-99
 
 .. code-block:: default
 
@@ -152,15 +176,17 @@ first group.
  .. code-block:: none
 
 
-    GradientMaps(alignment='joint', approach='dm', kernel='normalized_angle',
-                 n_components=10, random_state=None)
+    GradientMaps(alignment='joint', kernel='normalized_angle')
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 100-103
+
 Here, `gp` contains the Procrustes aligned data and `gj` contains the joint
 aligned data. Let’s plot them, but in separate figures to keep things
 organized.
 
+.. GENERATED FROM PYTHON SOURCE LINES 103-115
 
 .. code-block:: default
 
@@ -180,12 +206,14 @@ organized.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial2_003.png
+    :alt: plot tutorial2
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 116-128
 
 .. code-block:: default
 
@@ -205,12 +233,14 @@ organized.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial2_004.png
+    :alt: plot tutorial2
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 129-141
 
 .. code-block:: default
 
@@ -230,12 +260,15 @@ organized.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial2_005.png
+    :alt: plot tutorial2
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 142-155
+
 Although in this example, we don't see any big differences, if the input data
 was less similar, alignments may also resolve changes in the order of the
 gradients. However, you should always inspect the output of an alignment;
@@ -250,6 +283,7 @@ The first alignment iteration will then be to the reference. For purposes of
 this example, we will use the gradient of the hold-out group as the
 reference.
 
+.. GENERATED FROM PYTHON SOURCE LINES 155-163
 
 .. code-block:: default
 
@@ -272,11 +306,12 @@ reference.
  .. code-block:: none
 
 
-    GradientMaps(alignment='procrustes', approach='le', kernel='normalized_angle',
-                 n_components=10, random_state=None)
+    GradientMaps(alignment='procrustes', approach='le', kernel='normalized_angle')
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 164-176
+
 The gradients in `galign.aligned_` are now aligned to the reference
 gradients.
 
@@ -290,6 +325,7 @@ First, let's load the example data of microstructural profile covariance
 `(Paquola et al., 2019) <https://journals.plos.org/plosbiology/article?
 id=10.1371/journal.pbio.3000284>`_ and functional connectivity.
 
+.. GENERATED FROM PYTHON SOURCE LINES 176-195
 
 .. code-block:: default
 
@@ -316,15 +352,19 @@ id=10.1371/journal.pbio.3000284>`_ and functional connectivity.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial2_006.png
+    :alt: plot tutorial2
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 196-198
+
 In order to fuse the matrices, we simply pass the matrices to the fusion
 command which will rescale and horizontally concatenate the matrices.
 
+.. GENERATED FROM PYTHON SOURCE LINES 198-225
 
 .. code-block:: default
 
@@ -362,9 +402,12 @@ command which will rescale and horizontally concatenate the matrices.
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 226-228
+
 We then use this output in the fit function. This will convert the long
 horizontal array into a square affinity matrix, and then perform embedding.
 
+.. GENERATED FROM PYTHON SOURCE LINES 228-243
 
 .. code-block:: default
 
@@ -387,12 +430,15 @@ horizontal array into a square affinity matrix, and then perform embedding.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial2_007.png
+    :alt: plot tutorial2
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 244-254
+
 .. note::
   The mpc matrix presented here matches the subject cohort of `(Paquola et
   al., 2019) <https://journals.plos.org/plosbiology/article?id=10.1371/
@@ -407,7 +453,7 @@ null hypothesis testing of comparisons between gradients and other markers.
 
 .. rst-class:: sphx-glr-timing
 
-   **Total running time of the script:** ( 0 minutes  8.340 seconds)
+   **Total running time of the script:** ( 0 minutes  6.461 seconds)
 
 
 .. _sphx_glr_download_python_doc_auto_examples_plot_tutorial2.py:
@@ -420,13 +466,13 @@ null hypothesis testing of comparisons between gradients and other markers.
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-python
 
      :download:`Download Python source code: plot_tutorial2.py <plot_tutorial2.py>`
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-jupyter
 
      :download:`Download Jupyter notebook: plot_tutorial2.ipynb <plot_tutorial2.ipynb>`
 
diff --git a/docs/python_doc/auto_examples/plot_tutorial3.ipynb b/docs/python_doc/auto_examples/plot_tutorial3.ipynb
index a3e6544..80c9585 100644
--- a/docs/python_doc/auto_examples/plot_tutorial3.ipynb
+++ b/docs/python_doc/auto_examples/plot_tutorial3.ipynb
@@ -15,14 +15,14 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nTutorial 3: Null models for gradient significance\n==================================================\nIn this tutorial we assess the significance of correlations between the first\ncanonical gradient and data from other modalities (curvature, cortical\nthickness and T1w/T2w image intensity). A normal test of the significance of\nthe correlation cannot be used, because the spatial auto-correlation in MRI\ndata may bias the test statistic. In this tutorial we will show two approaches\nfor null hypothesis testing: spin permutations and Moran spectral\nrandomization.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>When using either approach to compare gradients to non-gradient markers,\n    we recommend randomizing the non-gradient markers as these randomizations\n    need not maintain the statistical independence between gradients.</p></div>\n"
+        "\n# Tutorial 3: Null models for gradient significance\nIn this tutorial we assess the significance of correlations between the first\ncanonical gradient and data from other modalities (curvature, cortical\nthickness and T1w/T2w image intensity). A normal test of the significance of\nthe correlation cannot be used, because the spatial auto-correlation in MRI\ndata may bias the test statistic. In this tutorial we will show three\napproaches for null hypothesis testing: spin permutations, Moran spectral\nrandomization, and autocorrelation-preserving surrogates based on variogram\nmatching.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>When using either approach to compare gradients to non-gradient markers,\n    we recommend randomizing the non-gradient markers as these randomizations\n    need not maintain the statistical independence between gradients.</p></div>\n"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Spin Permutations\n------------------------------\n\nHere, we use the spin permutations approach previously proposed in\n`(Alexander-Bloch et al., 2018)\n<https://www.sciencedirect.com/science/article/pii/S1053811918304968>`_,\nwhich preserves the auto-correlation of the permuted feature(s) by rotating\nthe feature data on the spherical domain.\nWe will start by loading the conte69 surfaces for left and right hemispheres,\ntheir corresponding spheres, midline mask, and t1w/t2w intensity as well as\ncortical thickness data, and a template functional gradient.\n\n"
+        "## Spin Permutations\n\nHere, we use the spin permutations approach previously proposed in\n`(Alexander-Bloch et al., 2018)\n<https://www.sciencedirect.com/science/article/pii/S1053811918304968>`_,\nwhich preserves the auto-correlation of the permuted feature(s) by rotating\nthe feature data on the spherical domain.\nWe will start by loading the conte69 surfaces for left and right hemispheres,\ntheir corresponding spheres, midline mask, and t1w/t2w intensity as well as\ncortical thickness data, and a template functional gradient.\n\n"
       ]
     },
     {
@@ -87,7 +87,7 @@
       },
       "outputs": [],
       "source": [
-        "# sphinx_gallery_thumbnail_number = 2\n# Plot some rotations\nplot_hemispheres(surf_lh, surf_rh, array_name=t1wt2w_rotated[:3], size=(1200, 600),\n                 cmap='viridis', nan_color=(0.5, 0.5, 0.5, 1), color_bar=True,\n                 zoom=1.55, label_text=['Rot0', 'Rot1', 'Rot2'])"
+        "# Plot some rotations\nplot_hemispheres(surf_lh, surf_rh, array_name=t1wt2w_rotated[:3], size=(1200, 600),\n                 cmap='viridis', nan_color=(0.5, 0.5, 0.5, 1), color_bar=True,\n                 zoom=1.55, label_text=['Rot0', 'Rot1', 'Rot2'])"
       ]
     },
     {
@@ -105,14 +105,14 @@
       },
       "outputs": [],
       "source": [
-        "from matplotlib import pyplot as plt\nfrom scipy.stats import spearmanr\n\nfig, axs = plt.subplots(1, 2, figsize=(9, 3.5))\n\nfeats = {'t1wt2w': t1wt2w, 'thickness': thickness}\nrotated = {'t1wt2w': t1wt2w_rotated, 'thickness': thickness_rotated}\n\nr_spin = np.empty(n_rand)\nmask = ~np.isnan(thickness)\nfor k, (fn, feat) in enumerate(feats.items()):\n    r_obs, pv_obs = spearmanr(feat[mask], embedding[mask])\n\n    # Compute perm pval\n    for i, perm in enumerate(rotated[fn]):\n        mask_rot = mask & ~np.isnan(perm)  # Remove midline\n        r_spin[i] = spearmanr(perm[mask_rot], embedding[mask_rot])[0]\n    pv_spin = np.mean(np.abs(r_spin) >= np.abs(r_obs))\n\n    # Plot null dist\n    axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))\n    axs[k].axvline(r_obs, lw=2, ls='--', color='k')\n    axs[k].set_xlabel('Correlation with {}'.format(fn))\n    if k == 0:\n        axs[k].set_ylabel('Density')\n\n    print('{}:\\n Obs : {:.5e}\\n Spin: {:.5e}\\n'.\n          format(fn.capitalize(), pv_obs, pv_spin))\n\nfig.tight_layout()\nplt.show()"
+        "from matplotlib import pyplot as plt\nfrom scipy.stats import spearmanr\n\nfig, axs = plt.subplots(1, 2, figsize=(9, 3.5))\n\nfeats = {'t1wt2w': t1wt2w, 'thickness': thickness}\nrotated = {'t1wt2w': t1wt2w_rotated, 'thickness': thickness_rotated}\n\nr_spin = np.empty(n_rand)\nmask = ~np.isnan(thickness)\nfor k, (fn, feat) in enumerate(feats.items()):\n    r_obs, pv_obs = spearmanr(feat[mask], embedding[mask])\n\n    # Compute perm pval\n    for i, perm in enumerate(rotated[fn]):\n        mask_rot = mask & ~np.isnan(perm)  # Remove midline\n        r_spin[i] = spearmanr(perm[mask_rot], embedding[mask_rot])[0]\n    pv_spin = np.mean(np.abs(r_spin) >= np.abs(r_obs))\n\n    # Plot null dist\n    axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))\n    axs[k].axvline(r_obs, lw=2, ls='--', color='k')\n    axs[k].set_xlabel(f'Correlation with {fn}')\n    if k == 0:\n        axs[k].set_ylabel('Density')\n\n    print(f'{fn.capitalize()}:\\n Obs : {pv_obs:.5e}\\n Spin: {pv_spin:.5e}\\n')\n\nfig.tight_layout()\nplt.show()"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "It is interesting to see that both p-values increase when taking into\nconsideration the auto-correlation present in the surfaces. Also, we can see\nthat the correlation with thickness is no longer statistically significant\nafter spin permutations.\n\n\n\nMoran Spectral Randomization\n------------------------------\n\nMoran Spectral Randomization (MSR) computes Moran's I, a metric for spatial\nauto-correlation and generates normally distributed data with similar\nauto-correlation. MSR relies on a weight matrix denoting the spatial\nproximity of features to one another. Within neuroimaging, one\nstraightforward example of this is inverse geodesic distance i.e. distance\nalong the cortical surface.\n\nIn this example we will show how to use MSR to assess statistical\nsignificance between cortical markers (here curvature and cortical t1wt2w\nintensity) and the first functional connectivity gradient. We will start by\nloading the left temporal lobe mask, t1w/t2w intensity as well as cortical\nthickness data, and a template functional gradient\n\n"
+        "It is interesting to see that both p-values increase when taking into\nconsideration the auto-correlation present in the surfaces. Also, we can see\nthat the correlation with thickness is no longer statistically significant\nafter spin permutations.\n\n\n\n## Moran Spectral Randomization\n\nMoran Spectral Randomization (MSR) computes Moran's I, a metric for spatial\nauto-correlation and generates normally distributed data with similar\nauto-correlation. MSR relies on a weight matrix denoting the spatial\nproximity of features to one another. Within neuroimaging, one\nstraightforward example of this is inverse geodesic distance i.e. distance\nalong the cortical surface.\n\nIn this example we will show how to use MSR to assess statistical\nsignificance between cortical markers (here curvature and cortical t1wt2w\nintensity) and the first functional connectivity gradient. We will start by\nloading the left temporal lobe mask, t1w/t2w intensity as well as cortical\nthickness data, and a template functional gradient\n\n"
       ]
     },
     {
@@ -177,14 +177,61 @@
       },
       "outputs": [],
       "source": [
-        "fig, axs = plt.subplots(1, 2, figsize=(9, 3.5))\n\nfeats = {'t1wt2w': t1wt2w_tl, 'curvature': curv_tl}\nrand = {'t1wt2w': t1wt2w_rand, 'curvature': curv_rand}\n\nfor k, (fn, data) in enumerate(rand.items()):\n    r_obs, pv_obs = spearmanr(feats[fn], embedding_tl, nan_policy='omit')\n\n    # Compute perm pval\n    r_rand = np.asarray([spearmanr(embedding_tl, d)[0] for d in data])\n    pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))\n\n    # Plot null dist\n    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))\n    axs[k].axvline(r_obs, lw=2, ls='--', color='k')\n    axs[k].set_xlabel('Correlation with {}'.format(fn))\n    if k == 0:\n        axs[k].set_ylabel('Density')\n\n    print('{}:\\n Obs  : {:.5e}\\n Moran: {:.5e}\\n'.\n          format(fn.capitalize(), pv_obs, pv_rand))\n\nfig.tight_layout()\nplt.show()"
+        "fig, axs = plt.subplots(1, 2, figsize=(9, 3.5))\n\nfeats = {'t1wt2w': t1wt2w_tl, 'curvature': curv_tl}\nrand = {'t1wt2w': t1wt2w_rand, 'curvature': curv_rand}\n\nfor k, (fn, data) in enumerate(rand.items()):\n    r_obs, pv_obs = spearmanr(feats[fn], embedding_tl, nan_policy='omit')\n\n    # Compute perm pval\n    r_rand = np.asarray([spearmanr(embedding_tl, d)[0] for d in data])\n    pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))\n\n    # Plot null dist\n    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))\n    axs[k].axvline(r_obs, lw=2, ls='--', color='k')\n    axs[k].set_xlabel(f'Correlation with {fn}')\n    if k == 0:\n        axs[k].set_ylabel('Density')\n\n    print(f'{fn.capitalize()}:\\n Obs  : {pv_obs:.5e}\\n Moran: {pv_rand:.5e}\\n')\n\nfig.tight_layout()\nplt.show()"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "There are some scenarios where MSR results do not follow a normal\ndistribution. It is relatively simple to check whether this occurs in our\ndata by visualizing the null distributions. Check this interesting paper\nfor more information `(Burt et al., 2020) <https://www.biorxiv.org/content/\n10.1101/2020.02.18.955054v1>`_.\n\n"
+        "## Variogram Matching\n\nHere, we will repeat the same analysis using the variogram matching approach\npresented in `(Burt et al., 2020)\n<https://www.sciencedirect.com/science/article/pii/S1053811920305243>`_,\nwhich generates novel brainmaps with similar spatial autocorrelation to the\ninput data.\n\n\nWe will need a distance matrix that tells us what the spatial distance\nbetween our datapoints is. For this example, we will use geodesic distance.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "from brainspace.mesh.mesh_elements import get_immediate_distance\nfrom scipy.sparse.csgraph import dijkstra\n\n# Compute geodesic distance\ngd = get_immediate_distance(surf_lh, mask=mask_tl)\ngd = dijkstra(gd, directed=False)\n\nidx_sorted = np.argsort(gd, axis=1)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Now we've got everything we need to generate our surrogate datasets. By\ndefault, BrainSpace will use all available data to generate surrogate maps.\nHowever, this process is extremely computationally and memory intensive. When\nusing this method with more than a few hundred regions, we recommend\nsubsampling the data. This can be done using SampledSurrogateMaps instead of\nthe SurrogateMaps.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "from brainspace.null_models import SampledSurrogateMaps\n\n\nn_surrogate_datasets = 1000\n\n# Note: number samples must be greater than number neighbors\nnum_samples = 100\nnum_neighbors = 50\n\nssm = SampledSurrogateMaps(ns=num_samples, knn=num_neighbors, random_state=0)\nssm.fit(gd, idx_sorted)\n\nt1wt2w_surrogates = ssm.randomize(t1wt2w_tl, n_rep=n_surrogate_datasets)\ncurv_surrogates = ssm.randomize(curv_tl, n_rep=n_surrogate_datasets)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Similar to the previous case, we can now plot the results:\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import matplotlib.pyplot as plt\nfrom scipy.stats import spearmanr\nfig, axs = plt.subplots(1, 2, figsize=(9, 3.5))\n\nfeats = {'t1wt2w': t1wt2w_tl, 'curvature': curv_tl}\nrand = {'t1wt2w': t1wt2w_surrogates, 'curvature': curv_surrogates}\n\nfor k, (fn, data) in enumerate(rand.items()):\n    r_obs, pv_obs = spearmanr(feats[fn], embedding_tl, nan_policy='omit')\n\n    # Compute perm pval\n    r_rand = np.asarray([spearmanr(embedding_tl, d)[0] for d in data])\n    pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))\n\n    # Plot null dist\n    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))\n    axs[k].axvline(r_obs, lw=2, ls='--', color='k')\n    axs[k].set_xlabel(f'Correlation with {fn}')\n    if k == 0:\n        axs[k].set_ylabel('Density')\n\n    print(f'{fn.capitalize()}:\\n Obs      : {pv_obs:.5e}\\n '\n          f'Variogram: {pv_rand:.5e}\\n')\n\nfig.tight_layout()\nplt.show()"
       ]
     }
   ],
@@ -204,7 +251,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.6"
+      "version": "3.7.10"
     }
   },
   "nbformat": 4,
diff --git a/docs/python_doc/auto_examples/plot_tutorial3.py b/docs/python_doc/auto_examples/plot_tutorial3.py
index 04df1f6..f09e8c9 100644
--- a/docs/python_doc/auto_examples/plot_tutorial3.py
+++ b/docs/python_doc/auto_examples/plot_tutorial3.py
@@ -5,9 +5,10 @@
 canonical gradient and data from other modalities (curvature, cortical
 thickness and T1w/T2w image intensity). A normal test of the significance of
 the correlation cannot be used, because the spatial auto-correlation in MRI
-data may bias the test statistic. In this tutorial we will show two approaches
-for null hypothesis testing: spin permutations and Moran spectral
-randomization.
+data may bias the test statistic. In this tutorial we will show three
+approaches for null hypothesis testing: spin permutations, Moran spectral
+randomization, and autocorrelation-preserving surrogates based on variogram
+matching.
 
 .. note::
     When using either approach to compare gradients to non-gradient markers,
@@ -117,14 +118,13 @@
     pv_spin = np.mean(np.abs(r_spin) >= np.abs(r_obs))
 
     # Plot null dist
-    axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))
+    axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
     axs[k].axvline(r_obs, lw=2, ls='--', color='k')
-    axs[k].set_xlabel('Correlation with {}'.format(fn))
+    axs[k].set_xlabel(f'Correlation with {fn}')
     if k == 0:
         axs[k].set_ylabel('Density')
 
-    print('{}:\n Obs : {:.5e}\n Spin: {:.5e}\n'.
-          format(fn.capitalize(), pv_obs, pv_spin))
+    print(f'{fn.capitalize()}:\n Obs : {pv_obs:.5e}\n Spin: {pv_spin:.5e}\n')
 
 fig.tight_layout()
 plt.show()
@@ -208,22 +208,97 @@
     pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))
 
     # Plot null dist
-    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))
+    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
     axs[k].axvline(r_obs, lw=2, ls='--', color='k')
-    axs[k].set_xlabel('Correlation with {}'.format(fn))
+    axs[k].set_xlabel(f'Correlation with {fn}')
     if k == 0:
         axs[k].set_ylabel('Density')
 
-    print('{}:\n Obs  : {:.5e}\n Moran: {:.5e}\n'.
-          format(fn.capitalize(), pv_obs, pv_rand))
+    print(f'{fn.capitalize()}:\n Obs  : {pv_obs:.5e}\n Moran: {pv_rand:.5e}\n')
 
 fig.tight_layout()
 plt.show()
 
 
 ###############################################################################
-# There are some scenarios where MSR results do not follow a normal
-# distribution. It is relatively simple to check whether this occurs in our
-# data by visualizing the null distributions. Check this interesting paper
-# for more information `(Burt et al., 2020) <https://www.biorxiv.org/content/
-# 10.1101/2020.02.18.955054v1>`_.
+#
+#
+#
+#
+# Variogram Matching
+# ------------------
+#
+# Here, we will repeat the same analysis using the variogram matching approach
+# presented in `(Burt et al., 2020)
+# <https://www.sciencedirect.com/science/article/pii/S1053811920305243>`_,
+# which generates novel brainmaps with similar spatial autocorrelation to the
+# input data.
+#
+#
+# We will need a distance matrix that tells us what the spatial distance
+# between our datapoints is. For this example, we will use geodesic distance.
+
+
+from brainspace.mesh.mesh_elements import get_immediate_distance
+from scipy.sparse.csgraph import dijkstra
+
+# Compute geodesic distance
+gd = get_immediate_distance(surf_lh, mask=mask_tl)
+gd = dijkstra(gd, directed=False)
+
+idx_sorted = np.argsort(gd, axis=1)
+
+
+###############################################################################
+# Now we've got everything we need to generate our surrogate datasets. By
+# default, BrainSpace will use all available data to generate surrogate maps.
+# However, this process is extremely computationally and memory intensive. When
+# using this method with more than a few hundred regions, we recommend
+# subsampling the data. This can be done using SampledSurrogateMaps instead of
+# the SurrogateMaps.
+
+from brainspace.null_models import SampledSurrogateMaps
+
+
+n_surrogate_datasets = 1000
+
+# Note: number samples must be greater than number neighbors
+num_samples = 100
+num_neighbors = 50
+
+ssm = SampledSurrogateMaps(ns=num_samples, knn=num_neighbors, random_state=0)
+ssm.fit(gd, idx_sorted)
+
+t1wt2w_surrogates = ssm.randomize(t1wt2w_tl, n_rep=n_surrogate_datasets)
+curv_surrogates = ssm.randomize(curv_tl, n_rep=n_surrogate_datasets)
+
+
+###############################################################################
+# Similar to the previous case, we can now plot the results:
+
+import matplotlib.pyplot as plt
+from scipy.stats import spearmanr
+fig, axs = plt.subplots(1, 2, figsize=(9, 3.5))
+
+feats = {'t1wt2w': t1wt2w_tl, 'curvature': curv_tl}
+rand = {'t1wt2w': t1wt2w_surrogates, 'curvature': curv_surrogates}
+
+for k, (fn, data) in enumerate(rand.items()):
+    r_obs, pv_obs = spearmanr(feats[fn], embedding_tl, nan_policy='omit')
+
+    # Compute perm pval
+    r_rand = np.asarray([spearmanr(embedding_tl, d)[0] for d in data])
+    pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))
+
+    # Plot null dist
+    axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
+    axs[k].axvline(r_obs, lw=2, ls='--', color='k')
+    axs[k].set_xlabel(f'Correlation with {fn}')
+    if k == 0:
+        axs[k].set_ylabel('Density')
+
+    print(f'{fn.capitalize()}:\n Obs      : {pv_obs:.5e}\n '
+          f'Variogram: {pv_rand:.5e}\n')
+
+fig.tight_layout()
+plt.show()
diff --git a/docs/python_doc/auto_examples/plot_tutorial3.py.md5 b/docs/python_doc/auto_examples/plot_tutorial3.py.md5
index 243ef0b..4f24987 100644
--- a/docs/python_doc/auto_examples/plot_tutorial3.py.md5
+++ b/docs/python_doc/auto_examples/plot_tutorial3.py.md5
@@ -1,1 +1,1 @@
-60018ffdc2814a6791cc0edfb57ff6ac
\ No newline at end of file
+cc3447eed3c70fb93955827cb7353b4f
\ No newline at end of file
diff --git a/docs/python_doc/auto_examples/plot_tutorial3.rst b/docs/python_doc/auto_examples/plot_tutorial3.rst
index 0fc90ed..38542c0 100644
--- a/docs/python_doc/auto_examples/plot_tutorial3.rst
+++ b/docs/python_doc/auto_examples/plot_tutorial3.rst
@@ -1,7 +1,18 @@
-.. note::
-    :class: sphx-glr-download-link-note
 
-    Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial3.py>` to download the full example code
+.. DO NOT EDIT.
+.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
+.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
+.. "python_doc/auto_examples/plot_tutorial3.py"
+.. LINE NUMBERS ARE GIVEN BELOW.
+
+.. only:: html
+
+    .. note::
+        :class: sphx-glr-download-link-note
+
+        Click :ref:`here <sphx_glr_download_python_doc_auto_examples_plot_tutorial3.py>`
+        to download the full example code
+
 .. rst-class:: sphx-glr-example-title
 
 .. _sphx_glr_python_doc_auto_examples_plot_tutorial3.py:
@@ -13,15 +24,18 @@ In this tutorial we assess the significance of correlations between the first
 canonical gradient and data from other modalities (curvature, cortical
 thickness and T1w/T2w image intensity). A normal test of the significance of
 the correlation cannot be used, because the spatial auto-correlation in MRI
-data may bias the test statistic. In this tutorial we will show two approaches
-for null hypothesis testing: spin permutations and Moran spectral
-randomization.
+data may bias the test statistic. In this tutorial we will show three
+approaches for null hypothesis testing: spin permutations, Moran spectral
+randomization, and autocorrelation-preserving surrogates based on variogram
+matching.
 
 .. note::
     When using either approach to compare gradients to non-gradient markers,
     we recommend randomizing the non-gradient markers as these randomizations
     need not maintain the statistical independence between gradients.
 
+.. GENERATED FROM PYTHON SOURCE LINES 22-33
+
 Spin Permutations
 ------------------------------
 
@@ -34,6 +48,7 @@ We will start by loading the conte69 surfaces for left and right hemispheres,
 their corresponding spheres, midline mask, and t1w/t2w intensity as well as
 cortical thickness data, and a template functional gradient.
 
+.. GENERATED FROM PYTHON SOURCE LINES 33-53
 
 .. code-block:: default
 
@@ -64,8 +79,11 @@ cortical thickness data, and a template functional gradient.
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 54-55
+
 Let’s first generate some null data using spintest.
 
+.. GENERATED FROM PYTHON SOURCE LINES 55-69
 
 .. code-block:: default
 
@@ -90,8 +108,11 @@ Let’s first generate some null data using spintest.
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 70-71
+
 As an illustration of the rotation, let’s plot the original t1w/t2w data
 
+.. GENERATED FROM PYTHON SOURCE LINES 71-77
 
 .. code-block:: default
 
@@ -105,19 +126,22 @@ As an illustration of the rotation, let’s plot the original t1w/t2w data
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial3_001.png
+    :alt: plot tutorial3
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 78-79
+
 as well as a few rotated versions.
 
+.. GENERATED FROM PYTHON SOURCE LINES 79-86
 
 .. code-block:: default
 
 
-    # sphinx_gallery_thumbnail_number = 2
     # Plot some rotations
     plot_hemispheres(surf_lh, surf_rh, array_name=t1wt2w_rotated[:3], size=(1200, 600),
                      cmap='viridis', nan_color=(0.5, 0.5, 0.5, 1), color_bar=True,
@@ -128,12 +152,15 @@ as well as a few rotated versions.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial3_002.png
+    :alt: plot tutorial3
     :class: sphx-glr-single-img
 
 
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 88-99
+
 .. warning::
 
    With spin permutations, midline vertices (i.e,, NaNs) from both the
@@ -146,6 +173,7 @@ as well as a few rotated versions.
 Now we simply compute the correlations between the first gradient and the
 original data, as well as all rotated data.
 
+.. GENERATED FROM PYTHON SOURCE LINES 100-132
 
 .. code-block:: default
 
@@ -170,14 +198,13 @@ original data, as well as all rotated data.
         pv_spin = np.mean(np.abs(r_spin) >= np.abs(r_obs))
 
         # Plot null dist
-        axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))
+        axs[k].hist(r_spin, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
         axs[k].axvline(r_obs, lw=2, ls='--', color='k')
-        axs[k].set_xlabel('Correlation with {}'.format(fn))
+        axs[k].set_xlabel(f'Correlation with {fn}')
         if k == 0:
             axs[k].set_ylabel('Density')
 
-        print('{}:\n Obs : {:.5e}\n Spin: {:.5e}\n'.
-              format(fn.capitalize(), pv_obs, pv_spin))
+        print(f'{fn.capitalize()}:\n Obs : {pv_obs:.5e}\n Spin: {pv_spin:.5e}\n')
 
     fig.tight_layout()
     plt.show()
@@ -186,6 +213,7 @@ original data, as well as all rotated data.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial3_003.png
+    :alt: plot tutorial3
     :class: sphx-glr-single-img
 
 
@@ -207,6 +235,8 @@ original data, as well as all rotated data.
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 133-155
+
 It is interesting to see that both p-values increase when taking into
 consideration the auto-correlation present in the surfaces. Also, we can see
 that the correlation with thickness is no longer statistically significant
@@ -230,6 +260,7 @@ intensity) and the first functional connectivity gradient. We will start by
 loading the left temporal lobe mask, t1w/t2w intensity as well as cortical
 thickness data, and a template functional gradient
 
+.. GENERATED FROM PYTHON SOURCE LINES 155-168
 
 .. code-block:: default
 
@@ -253,10 +284,13 @@ thickness data, and a template functional gradient
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 169-172
+
 We will now compute the Moran eigenvectors. This can be done either by
 providing a weight matrix of spatial proximity between each vertex, or by
 providing a cortical surface. Here we’ll use a cortical surface.
 
+.. GENERATED FROM PYTHON SOURCE LINES 172-186
 
 .. code-block:: default
 
@@ -285,13 +319,15 @@ providing a cortical surface. Here we’ll use a cortical surface.
  .. code-block:: none
 
 
-    MoranRandomization(joint=False, n_rep=1000, n_ring=1, procedure='singleton',
-                       random_state=0, spectrum='nonzero', tol=1e-06)
+    MoranRandomization(n_rep=1000, random_state=0, tol=1e-06)
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 187-188
+
 Using the Moran eigenvectors we can now compute the randomized data.
 
+.. GENERATED FROM PYTHON SOURCE LINES 188-193
 
 .. code-block:: default
 
@@ -307,10 +343,13 @@ Using the Moran eigenvectors we can now compute the randomized data.
 
 
 
+.. GENERATED FROM PYTHON SOURCE LINES 194-197
+
 Now that we have the randomized data, we can compute correlations between
 the gradient and the real/randomised data and generate the non-parametric
 p-values.
 
+.. GENERATED FROM PYTHON SOURCE LINES 197-223
 
 .. code-block:: default
 
@@ -328,14 +367,13 @@ p-values.
         pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))
 
         # Plot null dist
-        axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(0.8, 0.8, 0.8))
+        axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
         axs[k].axvline(r_obs, lw=2, ls='--', color='k')
-        axs[k].set_xlabel('Correlation with {}'.format(fn))
+        axs[k].set_xlabel(f'Correlation with {fn}')
         if k == 0:
             axs[k].set_ylabel('Density')
 
-        print('{}:\n Obs  : {:.5e}\n Moran: {:.5e}\n'.
-              format(fn.capitalize(), pv_obs, pv_rand))
+        print(f'{fn.capitalize()}:\n Obs  : {pv_obs:.5e}\n Moran: {pv_rand:.5e}\n')
 
     fig.tight_layout()
     plt.show()
@@ -345,6 +383,7 @@ p-values.
 
 
 .. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial3_004.png
+    :alt: plot tutorial3
     :class: sphx-glr-single-img
 
 
@@ -366,16 +405,146 @@ p-values.
 
 
 
-There are some scenarios where MSR results do not follow a normal
-distribution. It is relatively simple to check whether this occurs in our
-data by visualizing the null distributions. Check this interesting paper
-for more information `(Burt et al., 2020) <https://www.biorxiv.org/content/
-10.1101/2020.02.18.955054v1>`_.
+.. GENERATED FROM PYTHON SOURCE LINES 224-236
+
+Variogram Matching
+------------------
+
+Here, we will repeat the same analysis using the variogram matching approach
+presented in `(Burt et al., 2020)
+<https://www.sciencedirect.com/science/article/pii/S1053811920305243>`_,
+which generates novel brainmaps with similar spatial autocorrelation to the
+input data.
+
+
+We will need a distance matrix that tells us what the spatial distance
+between our datapoints is. For this example, we will use geodesic distance.
+
+.. GENERATED FROM PYTHON SOURCE LINES 240-252
+
+.. code-block:: default
+
+
+
+    from brainspace.mesh.mesh_elements import get_immediate_distance
+    from scipy.sparse.csgraph import dijkstra
+
+    # Compute geodesic distance
+    gd = get_immediate_distance(surf_lh, mask=mask_tl)
+    gd = dijkstra(gd, directed=False)
+
+    idx_sorted = np.argsort(gd, axis=1)
+
+
+
+
+
+
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 253-259
+
+Now we've got everything we need to generate our surrogate datasets. By
+default, BrainSpace will use all available data to generate surrogate maps.
+However, this process is extremely computationally and memory intensive. When
+using this method with more than a few hundred regions, we recommend
+subsampling the data. This can be done using SampledSurrogateMaps instead of
+the SurrogateMaps.
+
+.. GENERATED FROM PYTHON SOURCE LINES 259-276
+
+.. code-block:: default
+
+
+    from brainspace.null_models import SampledSurrogateMaps
+
+
+    n_surrogate_datasets = 1000
+
+    # Note: number samples must be greater than number neighbors
+    num_samples = 100
+    num_neighbors = 50
+
+    ssm = SampledSurrogateMaps(ns=num_samples, knn=num_neighbors, random_state=0)
+    ssm.fit(gd, idx_sorted)
+
+    t1wt2w_surrogates = ssm.randomize(t1wt2w_tl, n_rep=n_surrogate_datasets)
+    curv_surrogates = ssm.randomize(curv_tl, n_rep=n_surrogate_datasets)
+
+
+
+
+
+
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 277-278
+
+Similar to the previous case, we can now plot the results:
+
+.. GENERATED FROM PYTHON SOURCE LINES 278-305
+
+.. code-block:: default
+
+
+    import matplotlib.pyplot as plt
+    from scipy.stats import spearmanr
+    fig, axs = plt.subplots(1, 2, figsize=(9, 3.5))
+
+    feats = {'t1wt2w': t1wt2w_tl, 'curvature': curv_tl}
+    rand = {'t1wt2w': t1wt2w_surrogates, 'curvature': curv_surrogates}
+
+    for k, (fn, data) in enumerate(rand.items()):
+        r_obs, pv_obs = spearmanr(feats[fn], embedding_tl, nan_policy='omit')
+
+        # Compute perm pval
+        r_rand = np.asarray([spearmanr(embedding_tl, d)[0] for d in data])
+        pv_rand = np.mean(np.abs(r_rand) >= np.abs(r_obs))
+
+        # Plot null dist
+        axs[k].hist(r_rand, bins=25, density=True, alpha=0.5, color=(.8, .8, .8))
+        axs[k].axvline(r_obs, lw=2, ls='--', color='k')
+        axs[k].set_xlabel(f'Correlation with {fn}')
+        if k == 0:
+            axs[k].set_ylabel('Density')
+
+        print(f'{fn.capitalize()}:\n Obs      : {pv_obs:.5e}\n '
+              f'Variogram: {pv_rand:.5e}\n')
+
+    fig.tight_layout()
+    plt.show()
+
+
+
+.. image:: /python_doc/auto_examples/images/sphx_glr_plot_tutorial3_005.png
+    :alt: plot tutorial3
+    :class: sphx-glr-single-img
+
+
+.. rst-class:: sphx-glr-script-out
+
+ Out:
+
+ .. code-block:: none
+
+    T1wt2w:
+     Obs      : 0.00000e+00
+     Variogram: 0.00000e+00
+
+    Curvature:
+     Obs      : 6.63802e-05
+     Variogram: 3.26000e-01
+
+
+
+
 
 
 .. rst-class:: sphx-glr-timing
 
-   **Total running time of the script:** ( 3 minutes  22.475 seconds)
+   **Total running time of the script:** ( 3 minutes  29.601 seconds)
 
 
 .. _sphx_glr_download_python_doc_auto_examples_plot_tutorial3.py:
@@ -388,13 +557,13 @@ for more information `(Burt et al., 2020) <https://www.biorxiv.org/content/
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-python
 
      :download:`Download Python source code: plot_tutorial3.py <plot_tutorial3.py>`
 
 
 
-  .. container:: sphx-glr-download
+  .. container:: sphx-glr-download sphx-glr-download-jupyter
 
      :download:`Download Jupyter notebook: plot_tutorial3.ipynb <plot_tutorial3.ipynb>`
 
diff --git a/docs/python_doc/auto_examples/sg_execution_times.rst b/docs/python_doc/auto_examples/sg_execution_times.rst
index e796ed6..9cd8e52 100644
--- a/docs/python_doc/auto_examples/sg_execution_times.rst
+++ b/docs/python_doc/auto_examples/sg_execution_times.rst
@@ -5,14 +5,14 @@
 
 Computation times
 =================
-**00:04.417** total execution time for **python_doc_auto_examples** files:
+**03:41.329** total execution time for **python_doc_auto_examples** files:
 
 +------------------------------------------------------------------------------------+-----------+--------+
-| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial0.py` (``plot_tutorial0.py``) | 00:04.417 | 0.0 MB |
+| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial3.py` (``plot_tutorial3.py``) | 03:29.601 | 0.0 MB |
 +------------------------------------------------------------------------------------+-----------+--------+
-| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial1.py` (``plot_tutorial1.py``) | 00:00.000 | 0.0 MB |
+| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial2.py` (``plot_tutorial2.py``) | 00:06.461 | 0.0 MB |
 +------------------------------------------------------------------------------------+-----------+--------+
-| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial2.py` (``plot_tutorial2.py``) | 00:00.000 | 0.0 MB |
+| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial0.py` (``plot_tutorial0.py``) | 00:03.871 | 0.0 MB |
 +------------------------------------------------------------------------------------+-----------+--------+
-| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial3.py` (``plot_tutorial3.py``) | 00:00.000 | 0.0 MB |
+| :ref:`sphx_glr_python_doc_auto_examples_plot_tutorial1.py` (``plot_tutorial1.py``) | 00:01.397 | 0.0 MB |
 +------------------------------------------------------------------------------------+-----------+--------+
diff --git a/matlab/@GradientMaps/GradientMaps.m b/matlab/@GradientMaps/GradientMaps.m
index 94c2c2a..cb90844 100644
--- a/matlab/@GradientMaps/GradientMaps.m
+++ b/matlab/@GradientMaps/GradientMaps.m
@@ -27,6 +27,8 @@
 %   - random_state (default: nan)
 %       - Any valid input for MATLAB's "rng" function or nan for no
 %           initialization.
+%   - verbose (default: false)
+%       - Determines wheter non-warning/error messages will be displayed.
 %
 % For complete documentation, including descriptions of this object's
 % properties and methods please consult our <a
@@ -38,7 +40,7 @@
         method
         gradients
         aligned
-        lambda        
+        lambda
     end
     
     properties (Access = private)
@@ -47,10 +49,7 @@
     
     methods
         %% Constructor
-        function obj = GradientMaps(varargin)
-            disp('Launching BrainSpace, the gradient connectivity toolbox.');
-            disp('');
-            
+        function obj = GradientMaps(varargin)           
             % Parse input
             in_fun = @(x) isa(x,'char') || isa(x,'function_handle');
             p = inputParser;
@@ -59,17 +58,23 @@
             addParameter(p, 'alignment', 'none', in_fun);
             addParameter(p, 'n_components', 10, @isnumeric);
             addParameter(p, 'random_state', nan);
+            addParameter(p, 'verbose', false, @islogical); 
             
             parse(p, varargin{:});
             R = p.Results;
             
             % Set the properties
+            if R.verbose
+                disp('Launching BrainSpace, the gradient connectivity toolbox.');
+                disp('');
+            end
             obj = obj.set( ...
                 'kernel',       R.kernel, ...
                 'approach',     R.approach, ...
                 'alignment',    R.alignment, ...
                 'random_state', R.random_state, ...
-                'n_components', R.n_components);
+                'n_components', R.n_components, ...
+                'verbose',      R.verbose);
         end
     end
     methods(Access = private)
@@ -152,15 +157,21 @@
                     case 'n_components'
                         obj.method.n_components = varargin{ii+1};
                         change_string{end+1} = ['Set the number of requested components to: ' num2str(varargin{ii+1}) '.'];
+                    
+                    case 'verbose'
+                        obj.method.verbose = varargin{ii+1};
+                        change_string{end+1} = ['Verbose display was set to: ' mat2str(obj.method.verbose) '.'];
                         
                     otherwise
                         error('Unknown property. Valid properties are: ''connectivitymatrix'', ''kernel'', ''approach'', and ''nullmodel''.');
                 end
             end
-            for ii = 1:numel(change_string)
-                disp(change_string{ii})
+            if obj.method.verbose
+                for ii = 1:numel(change_string)
+                    disp(change_string{ii})
+                end
+                disp(' ')
             end
-            disp(' ')
         end
         % -------------------------------------
         % -------------------------------------
@@ -182,7 +193,9 @@
             end
             
             % Sparsify input data. 
-            disp(['Running with sparsity parameter: ' num2str(p.Results.sparsity)]);
+            if obj.method.verbose
+                disp(['Running with sparsity parameter: ' num2str(p.Results.sparsity)]);
+            end
             sparse_data = data;
             sparse_data(data < prctile(data,p.Results.sparsity)) = 0; 
             
@@ -201,7 +214,9 @@
                 case {'Pearson','Spearman'}
                     kernel_data = corr(sparse_data,'type',kernel);
                 case 'Gaussian'
-                    disp(['Running with gamma parameter: ' num2str(p.Results.gamma) '.']);
+                    if obj.method.verbose
+                        disp(['Running with gamma parameter: ' num2str(p.Results.gamma) '.']);
+                    end
                     kernel_data = exp(-p.Results.gamma .* squareform(pdist(sparse_data').^2));
                 case {'Cosine Similarity','Normalized Angle'}
                     cosine_similarity = 1-squareform(pdist(sparse_data','cosine'));
@@ -217,7 +232,9 @@
                         
             % Check for negative numbers.
             if any(kernel_data(:) < 0)
-                disp('Found negative numbers in the kernel matrix. These will be set to zero.');
+                if obj.method.verbose
+                    disp('Found negative numbers in the kernel matrix. These will be set to zero.');
+                end
                 kernel_data(kernel_data < 0) = 0; 
             end
             
@@ -288,11 +305,12 @@
                     [~, embedding, ~, ~, lambda] = pca(data);
                     embedding = embedding(:,1:obj.method.n_components);
                 case 'Laplacian Eigenmap'
-                    disp(['Requested ' num2str(obj.method.n_components) ' components.']);
                     [embedding, lambda] = laplacian_eigenmaps(data, obj.method.n_components);
                 case 'Diffusion Embedding'
-                    disp(['Running with alpha parameter: ' num2str(in.alpha)]);
-                    disp(['Running with diffusion time: ' num2str(in.diffusion_time)]);
+                    if obj.method.verbose
+                        disp(['Running with alpha parameter: ' num2str(in.alpha)]);
+                        disp(['Running with diffusion time: ' num2str(in.diffusion_time)]);
+                    end
                     [embedding, lambda] = diffusion_mapping(data, obj.method.n_components, in.alpha, in.diffusion_time);
                 otherwise
                     error('Unknown manifold technique.');
diff --git a/matlab/@GradientMaps/fit.m b/matlab/@GradientMaps/fit.m
index 7b4426f..7371ee7 100644
--- a/matlab/@GradientMaps/fit.m
+++ b/matlab/@GradientMaps/fit.m
@@ -56,22 +56,24 @@
 for ii = 1:numel(connectivity_matrix)
     connectivity_matrix{ii} = connectivity_matrix{ii}';
 end
-    
-disp('Running gradient analysis...');
-if isa(obj.method.kernel,'char')
-    disp(['Kernel: ' obj.method.kernel]);
-else
-    disp(['Kernel: ' func2str(obj.method.kernel)]);
-end
-if isa(obj.method.approach,'char')
-    disp(['Approach: ' obj.method.approach]);
-else
-    disp(['Approach: ' func2str(obj.method.approach)]);
-end
-if isa(obj.method.alignment,'char')
-    disp(['Alignment: ' obj.method.alignment]);
-else
-    disp(['Alignment: ' func2str(obj.method.alignment)]);
+
+if obj.method.verbose
+    disp('Running gradient analysis...');
+    if isa(obj.method.kernel,'char')
+        disp(['Kernel: ' obj.method.kernel]);
+    else
+        disp(['Kernel: ' func2str(obj.method.kernel)]);
+    end
+    if isa(obj.method.approach,'char')
+        disp(['Approach: ' obj.method.approach]);
+    else
+        disp(['Approach: ' func2str(obj.method.approach)]);
+    end
+    if isa(obj.method.alignment,'char')
+        disp(['Alignment: ' obj.method.alignment]);
+    else
+        disp(['Alignment: ' func2str(obj.method.alignment)]);
+    end
 end
 
 % Concatenate matrices for joint alignment. 
@@ -107,13 +109,17 @@
     else
         obj.gradients{ii} = obj.method.approach(kernel_data); 
     end
-    disp('Stored (unaligned) results in the gradients field.');
+    if obj.method.verbose
+        disp('Stored (unaligned) results in the gradients field.');
+    end
 end
 
 %Run the alignment
 if strcmp(obj.method.alignment,'Procrustes Analysis')
     obj.aligned = procrustes_alignment(obj.gradients,alignment_arg{:});
-    disp('Stored aligned results in the aligned field.');
+    if obj.method.verbose
+        disp('Stored aligned results in the aligned field.');
+    end
 elseif strcmp(obj.method.alignment,'Joint Alignment')
     for ii = 1:numel(size_connectivity)
         if ii == 1
diff --git a/matlab/analysis_code/@variogram/variogram.m b/matlab/analysis_code/@variogram/variogram.m
new file mode 100644
index 0000000..d2c135f
--- /dev/null
+++ b/matlab/analysis_code/@variogram/variogram.m
@@ -0,0 +1,339 @@
+classdef variogram
+% VARIOGRAM Generates brain maps with similar spatial autocorrelation.
+%
+%   obj = VARIOGRAM(D,varargin) initializes the variogram object based on
+%   symmetric distance matrix D. Valid name-value pairs are:
+%   deltas : 1-dimensional vector, default [0.1,0.2,...,0.9]
+%       Proportion of neighbors to include for smoothing, in (0, 1]
+%   kernel : character vector, default 'exp'
+%       Kernel with which to smooth permuted maps:
+%           'gaussian' : Gaussian function.
+%           'exp' : Exponential decay function.
+%           'invdist' : Inverse distance.
+%           'uniform' : Uniform weights (distance independent).
+%   pv : scalar (float), default 25
+%       Percentile of the pairwise distance distribution at which to
+%       truncate during variogram fitting.
+%   nh : scalar (integer), default 25
+%       Number of uniformly spaced distances at which to compute
+%       variograms.
+%   resample : logical, default false
+%       Resample surrogate maps' values from target brain map
+%   b : scalar (float) or nan, default nan
+%       Gaussian kernel bandwidth for variogram smoothing. If None, set to
+%       three times the spacing between variogram x-coordinates.
+%   random_state : any valid input for the rng() function. Set to nan for
+%       no random initialization.
+%   ns : scalar (integer) or inf, default inf. 
+%       Number of samples to use when subsampling the brainmap. Set to inf
+%       to use the entire brainmap. 
+%   knn : scalar (integer), default 1000
+%       Number of nearest neighbours to use when smoothing the map. knn
+%       must be smaller than ns. 
+%   num_workers: number of workers in the parallel pool, default 0.
+%       Sets the number of workers in a parallel pool. Only available if
+%       the parallel processing toolbox is installed. Note that if another
+%       pool is running with a different number of workers, then this pool
+%       will be shut down. The parallel pool is not closed at the end of
+%       the script. 
+%
+%   Public Methods: surrogates = obj.fit(x,n) generates n surrogate maps
+%   for brain map x. x must have the same length as D. n is set to 1000 by
+%   default. 
+%
+%   Example usage:
+%   x = rand(100,1);
+%   D = rand(100);
+%   D = D + D'; % make distance matrix symmetric.
+%   obj = variogram(D);
+%   surrogates = obj.fit(x); 
+%
+%   ADD A READTHEDOCS LINK!
+
+    %% Properties 
+    properties(SetAccess = private)
+        D
+        deltas
+        kernel
+        pv
+        nh
+        resample
+        b
+        random_state
+        ns
+        knn
+        num_workers
+    end
+    
+    properties(Hidden, Access = private)
+        disort
+    end
+
+    %% Public Methods
+    methods
+        %%% Constructor %%%
+        function obj = variogram(varargin)
+            % Deal with input.
+            is_square_numeric = @(x)size(x,1)==size(x,2) && isnumeric(x) && numel(size(x))==2;
+            p = inputParser;
+            addRequired( p, 'D'                         , is_square_numeric);
+            addParameter(p, 'deltas'    , 0.1:0.1:0.9   , @(x) all(x>0 & x<=1));
+            addParameter(p, 'kernel'    , 'exp'         , @valid_kernel);
+            addParameter(p, 'pv'        , 25            , @(x) x>0 && x<=100 && isscalar(x));
+            addParameter(p, 'nh'        , 25            , @(x)isscalar(x) && isnumeric(x));
+            addParameter(p, 'resample'  , false         , @islogical);
+            addParameter(p, 'b'         , nan           , @(x)isscalar(x) && isnumeric(x));
+            addParameter(p, 'random_state', nan);
+            addParameter(p, 'ns'        , inf           , @(x)isscalar(x) && isnumeric(x));
+            addParameter(p, 'knn'       , 1000          , @(x)isscalar(x) && isnumeric(x));
+            addParameter(p, 'num_workers', 0             , @(x)isscalar(x) && isnumeric(x));
+            
+            % Assign input to object properties. 
+            parse(p, varargin{:});
+            f = fieldnames(p.Results);
+            for ii = 1:numel(f)
+                obj.(f{ii}) = p.Results.(f{ii});
+            end  
+            
+            if ~isinf(obj.ns) && obj.ns >= size(obj.D,1)
+                error('The number of samples must be smaller than the number of nodes or infinite.');
+            end
+            
+            if obj.num_workers ~= 0 && ~isnan(obj.random_state)
+                warning('Due to the way that random number initialization works on parallel pools, the results of this script will not be identical on separate runs even with the random state defined.');
+            end
+        end
+        
+        %%% Fitting function %%%
+        function surrs = fit(obj,x,n)
+            % Start the parallel pool.
+            if obj.num_workers ~= 0
+                if exist('parpool.m', 'file')  
+                    p = gcp('nocreate');
+                        if isempty(p)
+                            parpool(obj.num_workers);   
+                        elseif p.NumWorkers ~= obj.num_workers
+                            delete(p);
+                            parpool(obj.num_workers);
+                        end
+                else
+                    warning('Could not find the parallel processing toolbox. Continuing without parallel processing.')
+                end
+            end
+            
+            % Set default n to 1000.
+            if ~exist('n','var')
+                n = 1000;
+            end
+            
+            % Set random state.
+            if ~isnan(obj.random_state)
+                rng(obj.random_state)
+            end
+                        
+            % Initialize variables.
+            surrs = zeros(size(x,1),n);
+            
+            % Only compute the sorted distance matrix once because sorting
+            % large matrices takes a while.
+            [~,obj.disort] = sort(obj.D,2);
+            
+            % Compute true variogram (dense). 
+            % Put inside a temporary structure to deal with parfor issues;
+            % parfor doesn't accept variables by the same name being
+            % defined both inside and outside the loop, even if they are
+            % defined in mutually exclusive if statements. 
+            if isinf(obj.ns)
+                v_out = obj.compute_variogram(x);
+                [tmp.utrunc,tmp.uidx,tmp.h] = obj.prepare_smooth_variogram(); 
+                tmp.smvar = obj.smooth_variogram(tmp.utrunc,tmp.uidx,v_out,tmp.h);
+            else
+                % Par-for wants a variable called "tmp" even if its unused.
+                tmp = struct();
+            end
+           
+            % Generate surrograte maps.
+            parfor (ii = 1:n, obj.num_workers)
+            %for ii = 1:n
+                % Initialize variables
+                aopt = nan;
+                bopt = nan;
+                rsquared_best = -inf; 
+                idxopt = nan;
+                
+                % Compute true variogram and permuted map (sampled).
+                if ~isinf(obj.ns) %#ok<PFBNS>
+                    perm = randperm(size(x,1),obj.ns);
+                    v = obj.compute_variogram(x(perm));
+                    [utrunc,uidx,h] = obj.prepare_smooth_variogram(perm);
+                    smvar = obj.smooth_variogram(utrunc,uidx,v,h);                 
+                else
+                    perm = nan; % This line prevents MATLAB from throwing a par-for warning - the line itself doesn't do anything. 
+                    utrunc = tmp.utrunc; %#ok<PFBNS>
+                    uidx = tmp.uidx;
+                    h = tmp.h;
+                    smvar = tmp.smvar;
+                end
+                
+                % Compute permuted map.
+                [x_perm,~] = obj.permute_map(x,~isnan(x));
+                
+                for jj = 1:numel(obj.deltas)
+                    % Smooth permuted map.
+                    sm_x_perm = obj.smooth_map(x_perm,jj);
+                    
+                    % Calculate empirical variogram
+                    if ~isinf(obj.ns)
+                        vperm = obj.compute_variogram(sm_x_perm(perm));
+                    else
+                        vperm = obj.compute_variogram(sm_x_perm);
+                    end
+                    smvar_perm = obj.smooth_variogram(utrunc,uidx,vperm,h);
+                    
+                    % Fit linear regression between smoothed variograms
+                    [alpha,beta,rsquared] = obj.local_regression(smvar_perm,smvar);
+                    
+                    if rsquared > rsquared_best
+                        aopt = alpha; 
+                        bopt = beta;
+                        idxopt = jj;
+                    end      
+                end
+                
+                % Transform and smooth permuted map using best-fit parameters.
+                sm_xperm_best = obj.smooth_map(x_perm,idxopt);
+                surrs(:,ii) = sqrt(abs(bopt)) * sm_xperm_best + sqrt(abs(aopt)) * randn(numel(x),1);
+            end
+            
+            % Resample to the input values. 
+            if obj.resample
+                sorted_map = sort(x);
+                for ii = 1:n
+                    [~,idx_resample] = ismember(surrs(:,ii),sort(surrs(:,ii)));
+                    surrs(:,ii) = sorted_map(idx_resample);
+                end
+            end
+        end
+    end
+    
+    %% Private methods
+    methods(Access = private)
+        function smvar = smooth_variogram(obj,utrunc,uidx,v,h)  
+            % Smooths the variograms. 
+            
+            vtrunc = v(uidx);
+                        
+            % Auto-estimate b if it's nan. 
+            if isnan(obj.b)
+                b_valid = 3 * (h(end)-h(1)) / (numel(h)-1);
+            else
+                b_valid = obj.b;
+            end
+            
+            du = abs(utrunc - h);
+            w = exp(-(2.68 * du / b_valid).^2/2);
+            smvar = sum(w.*vtrunc) ./ sum(w);
+        end
+        
+        function [utrunc,uidx,h] = prepare_smooth_variogram(obj,idx)
+            % Originally part of smooth_variogram. 
+            % Set to a separate function as this only needs to run once
+            % per permutation. 
+            if nargin > 1
+                D_subsample = obj.D(idx,idx);
+            else
+                D_subsample = obj.D;
+            end
+            
+            upper_triangle = triu(ones(size(D_subsample),'logical'),1);
+            u = D_subsample(upper_triangle);
+            uidx = u < prctile(u,obj.pv);
+            utrunc = u(uidx);
+            h = linspace(min(utrunc(:)),max(utrunc(:)),obj.nh);
+        end
+        
+        function sm_map = smooth_map(obj,x,idx_delta)
+            % Smooth x using delta proportion of nearest neighbors.     
+            
+            % Get distances and values of k nearest neighbours
+            if isinf(obj.ns)
+                k = floor(obj.deltas(idx_delta) * size(obj.D,1));
+            else
+                k = floor(obj.deltas(idx_delta) * obj.knn);
+            end
+            jkn = obj.disort(:,2:k+1);
+            jkn_idx = jkn + (0:size(obj.D,1):(size(jkn,1)-1)*size(obj.D,1))'; % Convert row indices to matrix indices. 
+            dkn = obj.D(jkn_idx);            
+            xkn = x(jkn);
+            
+            % Compute kernel weights and smoothed map.
+            weights = obj.smoothing_kernel(obj.kernel, dkn);
+            sm_map = sum(weights .* xkn,2) ./ sum(weights,2);
+        end
+    end
+    
+    methods(Static, Access = private)
+        function [data_perm,mask_perm] = permute_map(data,mask)
+            % Return randomly permuted brain map.
+            perm_idx = randperm(numel(data));
+            data_perm = data(perm_idx);
+            if nargin > 1
+                mask_perm = mask(perm_idx);
+            end
+        end
+        
+        function v = compute_variogram(x)
+            % Computes variograms based on input data x.
+            diff_ij = x-x';
+            upper_triangle = triu(ones(size(diff_ij),'logical'),1);
+            v = 0.5 * diff_ij(upper_triangle).^2;
+        end
+        
+        function smooth_d = smoothing_kernel(type,d)
+            % Applies a smoothing kernel to distance matrix d.
+            if isa(type, 'function_handle')
+                smooth_d = type(d); 
+            else
+                switch lower(type)
+                    case 'exp'
+                        smooth_d = exp(-d ./ max(d,[],2));
+                    case 'gaussian'
+                        smooth_d = exp(-1.25 * (d ./ max(d,[],2)));
+                    case 'invdist'
+                        smooth_d = d.^-1;
+                    case 'uniform'
+                        smooth_d = ones(size(d)) / size(d,1);
+                    otherwise
+                        error('Unknown kernel type')
+                end
+            end
+        end
+        
+        function [intercept,slope,rsquared] = local_regression(x,y)
+            % Computes the intercept, beta and resquared of a regression of
+            % x on y. 
+            
+            % Ascertain column vectors.
+            if ~isvector(x) || ~isvector(y)
+                error('x and y must be vectors');
+            end              
+            x = x(:); y = y(:);
+            
+            % Compute betas
+            x = [ones(size(x,1),1),x];
+            B = (x' * x) ^ -1 * x' * y; 
+            intercept = B(1); 
+            slope = B(2);
+            
+            % Compute R^2
+            rsquared = 1 - sum((y - (intercept + slope * x(:,2))).^2) / sum((y-mean(y)).^2); 
+        end
+    end
+end
+
+function is_valid = valid_kernel(x)
+    is_valid = isa(x, 'function_handle');
+    if ~is_valid
+        is_valid = ismember(lower(x), {'exp','gaussian','uniform','invdist'});
+    end
+end
diff --git a/matlab/analysis_code/procrustes_alignment.m b/matlab/analysis_code/procrustes_alignment.m
index 7b38a2e..922b37d 100644
--- a/matlab/analysis_code/procrustes_alignment.m
+++ b/matlab/analysis_code/procrustes_alignment.m
@@ -1,4 +1,4 @@
-function aligned = procrustes_alignment(gradients,varargin)
+function [aligned, xfms] = procrustes_alignment(gradients,varargin)
 % PROCRUSTES_ALIGNMENT   Performs a Procrustes alignment between gradients.
 %
 %   aligned = PROCRUSTES_ALIGNMENT(gradients,varargin) performs a singular
@@ -28,7 +28,7 @@
     first_target = p.Results.reference;
 end
 
-realigned = procrustes_analysis(embs, p.Results.nIterations, first_target);
+[realigned, xfms] = procrustes_analysis(embs, p.Results.nIterations, first_target);
 
 % Store aligned data in cell format. 
 for ii = 1:size(realigned,3)
diff --git a/matlab/analysis_code/spin_permutations.m b/matlab/analysis_code/spin_permutations.m
index d3f2b3e..6560fc8 100644
--- a/matlab/analysis_code/spin_permutations.m
+++ b/matlab/analysis_code/spin_permutations.m
@@ -8,8 +8,6 @@
 %   the aforementioned. n_rep denotes the number of spins to perform.
 %
 %   Name-Value pairs
-%       'parcellation': an n-by-1 vector containing parcellations on the
-%                       surface.
 %       'surface_algorithm': The program used to construct the surface; 'FreeSurfer'
 %               (default) and 'CIVET' are supported. 
 %
@@ -21,11 +19,9 @@
 
 % Make sure input is in the correct format. 
 p = inputParser;
-addParameter(p, 'parcellation', []);
 addParameter(p, 'surface_algorithm', 'freesurfer', @ischar);
 addParameter(p, 'random_state', nan);
 parse(p, varargin{:});
-parcellation = p.Results.parcellation;
 type = lower(p.Results.surface_algorithm);
 
 if ~isnan(p.Results.random_state)
@@ -44,12 +40,6 @@
     error('Found a different number of data cells than sphere cells.');
 end
 
-if ~isempty(parcellation)
-    if ~iscell(parcellation)
-        parcellation = {parcellation};
-    end
-end
-
 if numel(spheres) > 2
     error('Does not support more than two surfaces.')
 end
@@ -63,19 +53,6 @@
 for ii = 1:numel(spheres)
     S{ii} = convert_surface(spheres{ii});
     vertices{ii} = S{ii}.coord';    
-
-    % If parcellated data on sphere, get centroids of vertices.
-    if ~isempty(parcellation)
-        % Get Euclidean mean of points within each parcel.
-        euclidean_mean = labelmean(vertices{ii}, parcellation{ii},'ignorewarning')';
-        
-        % parcellationmean adds nan-columns if the parcellations are not 1:N. 
-        euclidean_mean(any(isnan(euclidean_mean),2),:) = [];
-        
-        % Find the centroid i.e. closest point to the Euclidean mean. 
-        idx = knnsearch(vertices{ii},euclidean_mean);
-        vertices{ii} = vertices{ii}(idx,:);   
-    end
     
     % Check for correct data dimensions. 
     if size(Y{ii},1) ~= size(vertices{ii},1)
@@ -122,7 +99,6 @@
         end
         rotated_vertices = vertices{ii}*rotation;
         nearest_neighbour = knnsearch(tree{ii}, rotated_vertices); % added 2019-06-18 see home page
-        %Y_rand{ii}= cat(3,Y_rand{ii}, Y{ii}(nearest_neighbour,:));
         Y_rand{ii}(:,:,j) = Y{ii}(nearest_neighbour,:);
     end
 end
diff --git a/matlab/plot_data/@plot_hemispheres/colorlimits.m b/matlab/plot_data/@plot_hemispheres/colorlimits.m
new file mode 100644
index 0000000..62feae8
--- /dev/null
+++ b/matlab/plot_data/@plot_hemispheres/colorlimits.m
@@ -0,0 +1,28 @@
+function colorlimits(obj,limits)
+% Sets color limits for plot_hemispheres. Limits must be a 2-element
+% numeric vector or n-by-2 matrix where n is the number of columns in
+% obj.data. The first column of limits denotes the minimum color limit and
+% the second the maximum. When limits is a 2-element vector, then the
+% limits are applied to all axes, with limits(1) as minimum and limits(2)
+% as maximum. 
+%
+% For more information, please consult our <a
+% href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
+
+% Check for correct input.
+if size(limits,1) ~= size(obj.data,2) && numel(limits) == 1
+    error('The number of color limits must be equal to the number of columns in the data matrix..');
+end
+
+% Set color limits for the axes and colorbar. 
+if numel(limits)==2
+    set(obj.handles.axes,'CLim',limits);
+    set(obj.handles.colorbar,'Limits',limits, 'Ticks',limits);
+else
+    for ii = 1:size(limits,1)
+        set(obj.handles.axes(ii,:),'CLim',limits(ii,:));
+        set(obj.handles.colorbar(ii),'Limits',limits(ii,:), 'Ticks',limits(ii,:));
+    end
+end
+drawnow
+end
\ No newline at end of file
diff --git a/matlab/plot_data/@plot_hemispheres/colormaps.m b/matlab/plot_data/@plot_hemispheres/colormaps.m
new file mode 100644
index 0000000..ac8d438
--- /dev/null
+++ b/matlab/plot_data/@plot_hemispheres/colormaps.m
@@ -0,0 +1,30 @@
+function colormaps(obj,maps)
+% Sets colormaps for plot_hemispheres. Maps must either be a n-by-3 color
+% map, or a cell array with the same number of elements as columns in
+% obj.data, each containing n-by-3 colormaps.
+%
+% For more information, please consult our <a
+% href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
+
+
+% If a single map is provided, set the colormap for the entire figure. 
+if numel(maps) == 1 || isnumeric(maps)
+    if iscell(maps)
+        maps = maps{1};
+    end
+    colormap(obj.handles.figure,maps);
+    return
+end
+
+% If multiple maps are provided, set a map for each axis. 
+if numel(maps) ~= size(obj.handles.axes,1)
+    error('The number of maps must be one or equal to the number of axes.');
+end
+
+for ii = 1:size(obj.handles.axes,1)
+    for jj = 1:size(obj.handles.axes,2)
+        colormap(obj.handles.axes(ii,jj),maps{ii});
+    end
+end
+drawnow
+end
\ No newline at end of file
diff --git a/matlab/plot_data/@plot_hemispheres/labels.m b/matlab/plot_data/@plot_hemispheres/labels.m
new file mode 100644
index 0000000..8122111
--- /dev/null
+++ b/matlab/plot_data/@plot_hemispheres/labels.m
@@ -0,0 +1,51 @@
+function labels(obj,varargin)
+% Sets hemisphere labels for plot_hemispheres. Labeltext is taken from
+% obj.labeltext. Text name-value pairs can be provided. If none are
+% provided, then the defaults are: 
+%   Rotation: 90
+%   Units: Normalized
+%   HorizontalAlignment: 'Center'
+%   FontName: DroidSans
+%   FontSize 18
+%   otherwise default text() parameters. 
+%
+% For more information, please consult our <a
+% href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
+
+% Get default properties
+properties = struct( ...
+    'Rotation',90, ...
+    'Units','Normalized', ...
+    'HorizontalAlignment','center', ...
+    'FontName', 'DroidSans', ...
+    'FontSize',18);
+
+for ii = 1:2:numel(varargin)
+    properties.(varargin{ii}) = varargin{ii+1}; 
+end
+
+% Convert property structure to cell array.
+C=[fieldnames(properties).'; struct2cell(properties).'];
+C=C(:).';
+
+% Set new labels. 
+for ii = 1:numel(obj.labeltext)
+    % If a label exists already modify only the properties.
+    % Rather complex if-statement because there's many exception cases that
+    % can error.
+    if ismember('text',fieldnames(obj.handles))
+        if numel(obj.handles.text) >= ii 
+            if ishandle(obj.handles.text(ii))
+                set(obj.handles.text(ii)        , ...
+                    'String'    , obj.labeltext{ii}, ...
+                    C{:});
+                continue
+            end
+        end
+    end
+    % If no labels exist yet, then simply build some. 
+    obj.handles.text(ii) = text(obj.handles.axes(ii,1), ...
+        -.2, .5, obj.labeltext{ii}, C{:});
+end
+drawnow
+end
diff --git a/matlab/plot_data/@plot_hemispheres/plot_hemispheres.m b/matlab/plot_data/@plot_hemispheres/plot_hemispheres.m
new file mode 100644
index 0000000..46a9136
--- /dev/null
+++ b/matlab/plot_data/@plot_hemispheres/plot_hemispheres.m
@@ -0,0 +1,177 @@
+classdef plot_hemispheres < handle
+    % PLOT_HEMISPHERES   Plots data on the cortical surface.
+    %
+    %   obj = PLOT_HEMISPHERES(data,surface,varargin) plots column vectors of the
+    %   n-by-m data matrix on surfaces provided in one or two element cell
+    %   arary surface. Vertices included in the surfaces must sum to n or the
+    %   length of the parcellation vector (see below). All handles to the
+    %   graphics objects are in the structure h.
+    %
+    %   Valid name-value pairs:
+    %       - 'labeltext': A length m cell array containing text labels to add to
+    %       the surfaces of each column vector.
+    %
+    %       - 'parcellation': A v-by-1 column vector containing a parcellation
+    %       scheme, with identical numbers denoting same parcel.
+    %
+    %       - 'views': A character vector containing the views you wish to
+    %       see. Valid options are  l(ateral), m(edial), i(nferior),
+    %       s(uperior), a(nterior), and p(osterior). When supplying a
+    %       single surface, the lateral/medial designations must be
+    %       inverted for the right hemisphere. Default: 'lm'.
+    %
+    %   For legacy usage: the handles formerly stored in obj are now stored
+    %   in obj.handles.
+    %
+    %   For more information, please consult our <a
+    %   href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
+    
+    properties
+        handles
+        labeltext
+    end
+    
+    properties(SetAccess = private)
+        data
+        surface
+        parcellation
+        views
+    end
+    
+    properties(Dependent)
+        plotted_data
+    end
+    
+    methods
+        function obj = plot_hemispheres(varargin)
+            % Initialize the object.
+            p = inputParser;
+            addRequired(p,'data')
+            addRequired(p,'surface')
+            addParameter(p,'labeltext',[]);
+            addParameter(p,'parcellation',[],@isnumeric);
+            addParameter(p,'views','lm', @ischar);
+            parse(p,varargin{:})
+            
+            % Set properties.
+            obj.handles = struct('figure',gobjects(0),'axes',gobjects(0),'trisurf',gobjects(0), ...
+                'camlight',gobjects(0),'colorbar',gobjects(0),'text',gobjects(0));
+            obj.surface = p.Results.surface;
+            obj.data = p.Results.data;
+            obj.parcellation = p.Results.parcellation;
+            obj.labeltext = p.Results.labeltext;
+            obj.views = p.Results.views;
+            
+            % Check whether input is correct.
+            if ~isempty(obj.parcellation)
+                obj.check_data_parcellation();
+            end
+            obj.check_surface();
+            
+            % Plot figure.
+            obj.plotter()
+        end
+        
+        %% Set methods
+        % Set surface.
+        function set.surface(obj,surface)
+            surface = convert_surface(surface,'Format','SurfStat');
+            if ~iscell(surface)
+                surface = {surface};
+            end
+            if numel(obj.surface) > 3
+                error('More than two surfaces are not accepted.');
+            end
+            obj.surface = surface;
+        end
+        
+        % Set parcellation.
+        function set.parcellation(obj,parcellation)
+            if ~isempty(parcellation)
+                if parcellation ~= round(parcellation)
+                    error('parcellation may only consist of zeros and positve integers.')
+                end
+                if ~isvector(parcellation) 
+                    error('parcellation must be a vector.')
+                end
+                parcellation = double(parcellation(:));
+            end
+            obj.parcellation = parcellation;
+        end
+        
+        % Set labeltext.
+        function set.labeltext(obj,labeltext)
+            if ~isempty(labeltext)
+                if ischar(labeltext)
+                    labeltext = {labeltext};
+                end
+                if ~iscell(labeltext) && ~isstring(labeltext)
+                    error('labeltext must be a char, cell, or string array');
+                end
+                if numel(labeltext) ~= size(obj.data,2)
+                    error('Must have as many labels as columns in the data matrix.')
+                end
+                obj.labeltext = labeltext;
+                if ismember('figure',fieldnames(obj.handles))
+                    if ishandle(obj.handles.figure)
+                        obj.labels();
+                    end
+                end
+            end
+        end
+        
+        function set.views(obj,views)
+            views = lower(views);
+            if numel(unique(views)) ~= numel(views)
+                warning('Detected duplicate characters in views; removing duplicates.');
+                views = unique(views);
+            end
+            for ii = 1:numel(views)
+                if ~ismember(views(ii),{'l','m','s','i','a','p'})
+                    error('Valid characters or views are: l(ateral), m(edial), i(nferior), s(uperior), a(nterior), p(osterior)');
+                end
+            end
+            obj.views = views;
+        end
+        
+        % Get to-plot data
+        function plotted_data = get.plotted_data(obj)
+            if isempty(obj.parcellation)
+                plotted_data = obj.data;
+            else
+                plotted_data = parcel2full(obj.data,obj.parcellation);
+                plotted_data(isnan(plotted_data)) = -Inf;
+            end
+        end
+        
+    end
+    
+    % Private methods for input checks. 
+    methods(Access=private)
+        function check_data_parcellation(obj)
+            % Check for correct size of parcellation.
+            sz = size(obj.data,1);
+            if max(obj.parcellation) ~= sz
+                error('The highest number in the parcellation scheme must be equivalent to length of the first dimension of the parcellated data.');
+            end
+            
+            % Check if all numbers are included.
+            if ~all(ismember(1:sz,obj.parcellation))
+                warning('Some parcel numbers are missing, some data may not be included in the output.')
+            end
+        end
+        
+        function check_surface(obj)
+            N1 = sum(cellfun(@(x)size(x.coord,2),obj.surface));
+            N2 = size(obj.plotted_data,1);
+            
+            if N1 ~= N2
+                error('Number of vertices on the surface and number of data points do not match.');
+            end
+            if size(obj.data,2) > 5
+                warning('Plotting more than five data vectors may result in off-screen axes');
+            end
+        end
+    end
+end
+
diff --git a/matlab/plot_data/@plot_hemispheres/private/make_surface_plot.m b/matlab/plot_data/@plot_hemispheres/private/make_surface_plot.m
new file mode 100644
index 0000000..919c1c0
--- /dev/null
+++ b/matlab/plot_data/@plot_hemispheres/private/make_surface_plot.m
@@ -0,0 +1,27 @@
+function [ax,srf,cam] = make_surface_plot(~,data,surf,pos,view,clim)
+% Constructs a single axis with a hemisphere in it. Called exclusively by
+% plotter. 
+%
+%   For more information, please consult our <a
+%   href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
+
+
+ax = axes('Position',pos);
+srf = trisurf(surf.tri, ...
+            surf.coord(1,:), ...
+            surf.coord(2,:), ...
+            surf.coord(3,:), ...
+            data, ...
+            'EdgeColor', 'None');
+material dull; lighting phong;
+colormap(parula(256));
+
+set(ax                                  , ...
+    'Visible'           , 'off'         , ...
+    'DataAspectRatio'   , [1 1 1]       , ...
+    'PlotBoxAspectRatio', [1 1 1]       , ...
+    'CLim'              , clim          , ...
+    'View'              , view          );
+
+% Add a camlight.
+cam = camlight();
\ No newline at end of file
diff --git a/matlab/plot_data/@plot_hemispheres/private/plotter.m b/matlab/plot_data/@plot_hemispheres/private/plotter.m
new file mode 100644
index 0000000..e8e6e54
--- /dev/null
+++ b/matlab/plot_data/@plot_hemispheres/private/plotter.m
@@ -0,0 +1,68 @@
+function plotter(obj)
+% Workhorse of plot hemispheres. Relies exclusively on properties of the
+% plot_hemispheres class and should never be called directly. 
+%
+% For more information, please consult our <a
+%   href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
+
+% Get data.
+data = obj.plotted_data;
+
+% Determine color limits
+for ii = 1:size(data,2)
+    clim(ii,:) = [min(data(~isinf(data(:,ii)),ii)); ...
+        max(data(~isinf(data(:,ii)),ii))];
+end
+
+% Account for data consisting of only one value
+clim = clim + (clim(:,1) == clim(:,1)) .* [-eps,eps];
+
+% Grab views.
+[views,surf_id] = obj.process_views;
+
+% Plot the surface(s).
+obj.handles.figure = figure('Color','white','Units','normalized','Position',[0 0 .9 .9]);
+for jj = 1:size(obj.data,2)
+    for ii = 1:numel(surf_id)
+        % Get the surface and associated data. 
+        if numel(obj.surface) == 1
+            D = obj.plotted_data(:,jj);
+            surf = obj.surface{1};
+        else
+            if surf_id{ii} == 1
+                D = obj.plotted_data(1:end/2,jj);
+                surf = obj.surface{surf_id{ii}};
+            elseif surf_id{ii} == 2
+                D = obj.plotted_data(end/2+1:end,jj);
+                surf = obj.surface{surf_id{ii}};
+            elseif surf_id{ii} == 3
+                D = obj.plotted_data(:,jj);
+                surf = combine_surfaces(obj.surface{1}, obj.surface{2}, 'surfstat');
+            end
+        end
+              
+        position = [-.07+ii*.11 .99-.18*jj .17 .17];
+        [obj.handles.axes(jj,ii), obj.handles.trisurf(jj,ii), obj.handles.camlight(jj,ii)] = ...
+            obj.make_surface_plot(D, surf, position, views{ii}, clim(jj,:));
+    end
+end
+
+% Add colorbars.
+for ii = 1:size(data,2)
+    obj.handles.colorbar(ii) = colorbar(obj.handles.axes(ii,end));
+    obj.handles.colorbar(ii).Position = [obj.handles.axes(ii,end).Position(1:2) + [.14 .045] ...
+        .007 .08];
+    obj.handles.colorbar(ii).Ticks = clim(ii,:);
+    obj.handles.colorbar(ii).FontName = 'DroidSans';
+    obj.handles.colorbar(ii).FontSize = 14;
+end
+
+% Adding labels before drawing makes the process look a bit weird. Better
+% to do this first.
+drawnow
+
+% Add labels
+if ~isempty(obj.labeltext)
+    obj.labels();
+end
+end
\ No newline at end of file
diff --git a/matlab/plot_data/@plot_hemispheres/private/process_views.m b/matlab/plot_data/@plot_hemispheres/private/process_views.m
new file mode 100644
index 0000000..ab08c66
--- /dev/null
+++ b/matlab/plot_data/@plot_hemispheres/private/process_views.m
@@ -0,0 +1,50 @@
+function [view_angles,surf_id] = process_views(obj)
+% Converts the input views character arrays to actual view angles and
+% defines which surfaces to include. Valid views contains the view angle(s)
+% in the top row, and the surface index (1,2, 3 for both surfaces) in the
+% bottom row.
+%
+% View priority is as follows:
+% Lateral/Medial for surface 1.
+% Medial/Lateral for surface 2.
+% Inferior
+% Superior
+% Anterior
+% Posterior
+%
+%   For more information, please consult our <a
+%   href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
+
+% First take care of the lateral/medial views as these get
+% mixed.
+view_angles = {};
+if contains(obj.views,'l') && contains(obj.views,'m')
+    view_angles(1:4) = {[-90 0],[90 0],[-90 0], [90 0]};
+    surf_id = {1,1,2,2};
+elseif contains(obj.views,'l')
+    view_angles(1:2) = {[-90 0],[90 0]};
+    surf_id = {1,2};
+elseif contains(obj.views,'m')
+    view_angles(1:2) = {[90 0],[-90 0]};
+    surf_id = {1,2};
+end
+
+% If only one surface is provided, remove the calls to surface
+% 2.
+if numel(obj.surface) == 1
+    view_angles(:,cellfun(@(x)x==2,surf_id)) = [];
+    surf_id(cellfun(@(x)(x==2),surf_id)) = []; 
+end
+
+% Now deal with the remaining views
+valid_views = {'i', [180 -90]; ...
+    's', [0 90]; ...
+    'a', [180 0]; ...
+    'p', [0 0]};
+for ii = 1:size(valid_views,1)
+    if contains(obj.views,valid_views{ii,1})
+        view_angles{end+1} = valid_views{ii,2};
+        surf_id{end+1} = 3; 
+    end
+end
+end
\ No newline at end of file
diff --git a/matlab/plot_data/plot_hemispheres.m b/matlab/plot_data/plot_hemispheres.m
deleted file mode 100644
index cb0fd1d..0000000
--- a/matlab/plot_data/plot_hemispheres.m
+++ /dev/null
@@ -1,127 +0,0 @@
-function h = plot_hemispheres(data,surface,varargin)
-% PLOT_HEMISPHERES   Plots data on the cortical surface.
-%
-%   h = PLOT_HEMISPHERES(data,surface,varargin) plots column vectors of the
-%   n-by-m data matrix on surfaces provided in one or two element cell
-%   arary surface. Vertices included in the surfaces must sum to n or the
-%   length of the parcellation vector (see below). All handles to the
-%   graphics objects are in the structure h. 
-%
-%   Valid name-value pairs:
-%       - 'labeltext': A length m cell array containing text labels to add to
-%       the surfaces of each column vector. 
-%       
-%       - 'parcellation': A v-by-1 column vector containing a parcellation
-%       scheme, with identical numbers denoting same parcel.
-%
-%   For more information, please consult our <a
-%   href="https://brainspace.readthedocs.io/en/latest/pages/matlab_doc/visualization/plot_hemispheres.html">ReadTheDocs</a>.
-
-for ii = 1:2:numel(varargin)
-    switch lower(varargin{ii})
-        case 'labeltext'
-            if ischar(varargin{ii+1})
-                label_text = {varargin{ii+1}};
-            else
-                label_text = varargin{ii+1};
-            end
-        case 'parcellation'
-            parcellation = varargin{ii+1}; 
-    otherwise
-        error('Unknown name-value pair.')
-    end
-end
-
-% Check if data is float.
-if ~isfloat(data)
-    data = double(data);
-end
-
-% If parcellated data is provided, bring it to the full mesh.
-if exist('parcellation','var')
-    data = parcel2full(data,parcellation);
-end
-data(isnan(data)) = -Inf;
-
-% Make sure surfaces are in SurfStat format.
-S = convert_surface(surface);
-if ~iscell(S); S = {S}; end
-if numel(S) > 3; error('More than two surfaces are not accepted.'); end
-
-% Check correct dimensions of data and surfaces
-N1 = sum(cellfun(@(x)size(x.coord,2),S));
-N2 = size(data,1);
-if N1 ~= N2; error('Number of vertices on the surface and number of data points do not match.'); end
-if size(data,2) > 4; warning('Plotting more than four data vectors may result in off-screen axes'); end
-
-% Split data over surfaces.
-for ii = 1:size(data,2)
-    D{1,ii} = data(1:size(S{1}.coord,2),ii);
-    if numel(S) == 2
-        D{2,ii} = data(size(S{1}.coord,2)+1:end,ii);
-    end
-end
-
-% Determine color limits
-for ii = 1:size(data,2)
-    clim(:,ii) = [min(data(~isinf(data(:,ii)),ii)); ...
-                  max(data(~isinf(data(:,ii)),ii))];
-end
-
-% Account for data consisting of only one value
-clim = clim + (clim(1,:) == clim(2,:)) .* [-eps;eps];
-
-% Plot the surface(s).
-h.figure = figure('Color','white','Units','normalized','Position',[0 0 .9 .9]);
-colormap(parula(256));
-for jj = 1:size(D,2)
-    for ii = 1:numel(S)*2
-        idx = ceil(ii/2);
-        h.axes(ii,jj) = axes('Position',[-.1+ii*.133 .9-.2*jj .2 .2]);
-        h.trisurf(ii,jj) = trisurf(S{idx}.tri, ...
-            S{idx}.coord(1,:), ...
-            S{idx}.coord(2,:), ...
-            S{idx}.coord(3,:), ...
-            D{idx,jj}, ...
-            'EdgeColor', 'None');
-        material dull; lighting phong;
-    end
-    set(h.axes(:,jj)                    , ...
-    'Visible'           , 'off'         , ...
-    'DataAspectRatio'   , [1 1 1]       , ...
-    'PlotBoxAspectRatio', [1 1 1]       , ...
-    'CLim'              , clim(:,jj)     );
-end
-
-% Add axes embelishments
-set(h.axes(1,:),'View',[-90 0]);
-set(h.axes(2,:),'View',[90 0]);
-if size(h.axes,1) == 4
-    set(h.axes(3,:),'View',[-90 0]);
-    set(h.axes(4,:),'View',[90 0]);
-end
-
-% Add a camlight. 
-for ii = 1:numel(h.axes)
-    axes(h.axes(ii));
-    h.camlight(ii) = camlight();
-end
-
-% Add colorbars. 
-for ii = 1:size(D,2)
-    h.cb(ii) = colorbar(h.axes(end,ii)); 
-    h.cb(ii).Position = [h.axes(end,ii).Position(1:2) + [.175 .065] ...
-                         .007 .08];
-    h.cb(ii).Ticks = clim(:,ii);
-    h.cb(ii).FontName = 'DroidSans';
-    h.cb(ii).FontSize = 14;
-end
-
-% Add labels   
-if exist('label_text','var')
-    for ii = 1:numel(label_text)
-        h.text(ii) = text(h.axes(1,ii),-.2,.5,label_text{ii},'Rotation',90, ...
-            'Units','Normalized','HorizontalAlignment','center', 'FontName', ...
-            'DroidSans','FontSize',18);
-    end
-end
diff --git a/matlab/surface_manipulation/convert_surface.m b/matlab/surface_manipulation/convert_surface.m
index af28c67..05f3274 100644
--- a/matlab/surface_manipulation/convert_surface.m
+++ b/matlab/surface_manipulation/convert_surface.m
@@ -4,8 +4,8 @@
 %
 %   S2 = CONVERT_SURFACE(S) converts surface S to SurfStat format. S can
 %   either be a file (.gii, .mat, .obj, Freesurfer), a loaded variable (in
-%   SurfStat or MATLAB format), or a cell array containing multiple of the
-%   former.
+%   SurfStat or MATLAB format), or a cell/string array containing multiple of 
+%   the former.
 %
 %   S2 = CONVERT_SURFACE(S,'format',F) allows for specifying the output
 %   format F, either 'SurfStat' or 'MATLAB'.
@@ -27,7 +27,7 @@
 path = p.Results.path;
 
 % If multiple surfaces provided in a cell, loop over all.
-if iscell(S)
+if iscell(S) || isstring(S)
     if ~isempty(path) && numel(S) > 1
         error('Multiple inputs are not supported for surface writing');
     end
diff --git a/shared/.DS_Store b/shared/.DS_Store
deleted file mode 100644
index 69609fb..0000000
Binary files a/shared/.DS_Store and /dev/null differ
diff --git a/shared/data/.DS_Store b/shared/data/.DS_Store
deleted file mode 100644
index d83c7cb..0000000
Binary files a/shared/data/.DS_Store and /dev/null differ


Test patch:
diff --git a/brainspace/tests/test_mesh.py b/brainspace/tests/test_mesh.py
index 28b3ac5..53973eb 100644
--- a/brainspace/tests/test_mesh.py
+++ b/brainspace/tests/test_mesh.py
@@ -4,6 +4,8 @@
 
 import os
 import numpy as np
+import shutil
+import gzip
 
 import vtk
 from vtk.util.vtkConstants import VTK_TRIANGLE, VTK_LINE, VTK_VERTEX
@@ -43,10 +45,19 @@ def test_io(ext):
     mio.write_surface(s, io_pth)
     s2 = mio.read_surface(io_pth)
 
+    io_gz_pth = os.path.join(root_pth, 'test_sphere_io.{ext}.gz').format(ext=ext)
+    with open(io_pth, 'rb') as f1:
+        with gzip.open(io_gz_pth, 'wb') as f2:
+            shutil.copyfileobj(f1, f2)
+    s3 = mio.read_surface(io_gz_pth)
+
     assert np.allclose(s.Points, s2.Points)
     assert np.all(s.GetCells2D() == s2.GetCells2D())
-
+    assert np.allclose(s.Points, s3.Points)
+    assert np.all(s.GetCells2D() == s3.GetCells2D())
+    
     os.remove(io_pth)
+    os.remove(io_gz_pth)
 
 
 @pytest.mark.skipif(nb is None, reason="Requires nibabel")
diff --git a/brainspace/tests/test_null_models.py b/brainspace/tests/test_null_models.py
index 9e86a89..0219d3a 100644
--- a/brainspace/tests/test_null_models.py
+++ b/brainspace/tests/test_null_models.py
@@ -3,7 +3,9 @@
 import pytest
 
 import numpy as np
-
+from scipy.spatial.distance import pdist, squareform
+from tempfile import gettempdir
+from os.path import join, exists
 import vtk
 
 from brainspace.vtk_interface import wrap_vtk
@@ -13,6 +15,9 @@
                                           MoranRandomization)
 from brainspace.null_models.spin import (_generate_spins, spin_permutations,
                                          SpinPermutations)
+# from brainspace.null_models.variogram import Base, Sampled,
+from brainspace.null_models.variogram import (txt2memmap, SurrogateMaps,
+                                              SampledSurrogateMaps)
 
 
 def test_moran():
@@ -147,3 +152,84 @@ def test_spin():
     assert np.all(r2[0] == r2bis[0])
     assert np.all(r2[1] == r2bis[1])
 
+
+def test_variogram_base():
+    # Base class
+    # Sphere with points as locations to build distance matrix
+    sphere = wrap_vtk(vtk.vtkSphereSource, radius=20, thetaResolution=20,
+                      phiResolution=10)
+    sphere = to_data(sphere)
+    points = sphere.GetPoints()
+    npoints = sphere.n_points
+    distmat = squareform(pdist(points))  # pairwise distance matrix
+
+    rs = np.random.RandomState(0)
+    brainmap = rs.randn(npoints, 1).flatten()
+    brainmap[0] = np.nan
+
+    # Generate surrogates
+    vsur = SurrogateMaps(random_state=672436)
+    vsur.fit(distmat)
+    surrs = vsur.randomize(brainmap, n_rep=10)
+
+    # np.random.seed(672436)
+    # gen = Base(brainmap, distmat)
+    # surrs2 = gen(n=10)
+    #
+    # assert np.allclose(surrs, surrs2, equal_nan=True)
+    assert surrs.shape == (10, npoints)
+    assert np.allclose(vsur._dist, distmat)
+    # assert np.allclose(gen.x.data[1:], brainmap[1:])
+    # assert np.isnan(gen.x.data[0])
+    assert not np.isnan(surrs).any()
+
+
+def test_variogram_sampled():
+    # Sampled class
+    # Sphere with points as locations to build distance matrix
+    sphere = wrap_vtk(vtk.vtkSphereSource, radius=20, thetaResolution=50,
+                      phiResolution=25)
+    sphere = to_data(sphere)
+    points = sphere.GetPoints()
+    npoints = sphere.n_points
+    distmat = squareform(pdist(points))  # pairwise distance matrix
+
+    rs = np.random.RandomState(0)
+    brainmap = rs.randn(npoints, 1).flatten()
+    brainmap[0] = np.nan
+
+    # Create a mask file
+    mask = np.zeros(npoints, dtype=int)
+    mask[[3, 4, 5]] = 1
+
+    # Save distmat and mask file
+    temp = gettempdir()
+    dist_file = join(temp, 'distmat.txt')
+    mask_file = join(temp, 'mask.txt')
+    np.savetxt(dist_file, distmat)
+    np.savetxt(mask_file, mask)
+    assert exists(dist_file) and exists(mask_file)
+
+    # Convert to memmap
+    files = txt2memmap(dist_file=dist_file, output_dir=temp, maskfile=mask_file)
+    assert exists(files['index'])
+    assert exists(files['distmat'])
+
+    masked_map = brainmap[np.where(mask == 0)]
+
+    # Generate surrogates
+    vsur = SampledSurrogateMaps(knn=100, ns=100, random_state=43)
+    vsur.fit(files['distmat'], files['index'])
+    surrs = vsur.randomize(masked_map, n_rep=5)
+
+    # np.random.seed(43)
+    # gen = Sampled(masked_map, files['distmat'], files['index'],
+    #               knn=100, ns=100)
+    # surrs2 = gen(n=5)
+    #
+    # assert np.allclose(surrs, surrs2, equal_nan=True)
+    assert vsur._dist.shape == (npoints-3, 100)
+    assert surrs.shape == (5, npoints-3)
+    # assert np.allclose(gen.x.data[1:], masked_map[1:])
+    # assert np.isnan(gen.x.data[0])
+    assert not np.isnan(surrs).any()