diff --git a/.circleci/config.yml b/.circleci/config.yml
new file mode 100644
index 0000000..c423da2
--- /dev/null
+++ b/.circleci/config.yml
@@ -0,0 +1,29 @@
+version: 2
+jobs:
+  pytest:
+    docker:
+      - image: circleci/python:3.7.3
+    steps:
+      - checkout
+      - run:
+          name: pytest
+          command: |
+            export PATH=/home/circleci/.local/bin:$PATH
+            mkdir test-results
+            pip install --user -r requirements.txt
+            set +e
+            coverage run -m pytest --junitxml=test-results/junit.xml
+            coverage report --include="./*" --omit="/home/circleci/.local/*"
+            coverage html --include="./*" --omit="/home/circleci/.local/*"
+      - store_test_results:
+          path: test-results
+      - store_artifacts:
+          path: test-results
+      - store_artifacts:
+          path: htmlcov
+
+workflows:
+  version: 2
+  unit_tests:
+    jobs:
+      - pytest
diff --git a/README.md b/README.md
index 4d6a3f2..99a19c6 100644
--- a/README.md
+++ b/README.md
@@ -1 +1,90 @@
 # CPAC_regtest_pack
+
+```bash
+$ python cpac_correlations_wf.py --help
+usage: cpac_correlations_wf.py [-h] [--old_outputs_path OLD_OUTPUTS_PATH]
+                               [--new_outputs_path NEW_OUTPUTS_PATH]
+                               [--s3_creds S3_CREDS]
+                               [--replacements REPLACEMENTS]
+                               [--corr_map CORR_MAP]
+                               [--working_dir WORKING_DIR]
+                               num_cores run_name
+
+positional arguments:
+  num_cores             number of cores to use - will calculate correlations
+                        in parallel if greater than 1
+  run_name              name for the correlations run
+
+optional arguments:
+  -h, --help            show this help message and exit
+  --old_outputs_path OLD_OUTPUTS_PATH
+                        path to a CPAC outputs directory - the folder
+                        containing the participant-ID labeled directories
+  --new_outputs_path NEW_OUTPUTS_PATH
+                        path to a CPAC outputs directory - the folder
+                        containing the participant-ID labeled directories
+  --s3_creds S3_CREDS   path to your AWS S3 credentials file
+  --replacements REPLACEMENTS
+                        text file containing strings you wish to have removed
+                        from the filepaths if they occur - place one on each
+                        line
+  --corr_map CORR_MAP   YAML file with already-calculated correlations, which
+                        can be provided if you only want to generate the box
+                        plots again
+  --working_dir WORKING_DIR
+                        if you are correlating two working directories of a
+                        single participant to check intermediates
+```
+
+```bash
+$ python correlation_matrix.py --help
+usage: correlation_matrix.py [-h] [--old_outputs_path OLD_OUTPUTS_PATH]
+                             [--old_outputs_software {C-PAC,fmriprep}]
+                             [--new_outputs_path NEW_OUTPUTS_PATH]
+                             [--new_outputs_software {C-PAC,fmriprep}]
+                             [--save] [--no-save]
+                             [--subject_list SUBJECT_LIST] [--session SESSION]
+                             [--feature_list FEATURE_LIST]
+                             num_cores run_name
+
+Create a correlation matrix between two C-PAC output directories.
+
+positional arguments:
+  num_cores             number of cores to use - will calculate correlations
+                        in parallel if greater than 1
+  run_name              name for the correlations run
+
+optional arguments:
+  -h, --help            show this help message and exit
+  --old_outputs_path OLD_OUTPUTS_PATH
+                        path to an outputs directory - the folder containing
+                        the participant-ID labeled directories
+  --old_outputs_software {C-PAC,fmriprep}
+                        (default: fmriprep)
+  --new_outputs_path NEW_OUTPUTS_PATH
+                        path to an outputs directory - the folder containing
+                        the participant-ID labeled directories
+  --new_outputs_software {C-PAC,fmriprep}
+                        (default: C-PAC)
+  --save                save matrices & heatmap (default)
+  --no-save             do not save matrices & heatmap
+  --subject_list SUBJECT_LIST
+                        (default: subjects in OLD_OUTPUTS_PATH sorted by
+                        session, subject ID). TODO: handle path to file
+  --session SESSION     limit to a single given session (integer)
+  --feature_list FEATURE_LIST
+                        TODO: handle path to file (default: ['GS', 'CSF',
+                        'WM', 'tCompCor0', 'aCompCor0', 'aCompCor1',
+                        'aCompCor2', 'aCompCor3', 'aCompCor4', 'FD'])
+
+The following features currently have available definitions to calculate Pearson's r between C-PAC and fmriprep:
+
+key       feature name              documentation link
+--------  ------------------------  ----------------------------------------------------------------------------------
+aCompCor  aCompCor                  https://fcp-indi.github.io/docs/user/nuisance.html#acompcor
+CSF       mean cerebrospinal fluid  https://fcp-indi.github.io/docs/user/nuisance.html#mean-white-matter-csf
+FD        framewise displacement    https://fcp-indi.github.io/docs/user/nuisance.html#regression-of-motion-parameters
+GS        global signal regression  https://fcp-indi.github.io/docs/user/nuisance.html#global-signal-regression
+tCompCor  tCompCor                  https://fcp-indi.github.io/docs/user/nuisance.html#tcompcor
+WM        mean white matter         https://fcp-indi.github.io/docs/user/nuisance.html#mean-white-matter-csf
+```
\ No newline at end of file
diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/configs/defaults.py b/configs/defaults.py
new file mode 100644
index 0000000..d342736
--- /dev/null
+++ b/configs/defaults.py
@@ -0,0 +1,53 @@
+feature_headers = {
+    'GS': {
+        'name': 'global signal regression',
+        'link': 'https://fcp-indi.github.io/docs/user/nuisance.html#'
+                'global-signal-regression',
+        'C-PAC': ['GlobalSignalMean0', 'GlobalSignal_mean'],
+        'fmriprep': 'global_signal'
+    },
+    'CSF': {
+        'name': 'mean cerebrospinal fluid',
+        'link': 'https://fcp-indi.github.io/docs/user/nuisance.html#'
+                'mean-white-matter-csf',
+        'C-PAC': ['CerebrospinalFluidMean0', 'CerebrospinalFluid_mean'],
+        'fmriprep': 'csf'
+    },
+    'WM': {
+        'name': 'mean white matter',
+        'link': 'https://fcp-indi.github.io/docs/user/nuisance.html#'
+                'mean-white-matter-csf',
+        'C-PAC': ['WhiteMatterMean0', 'WhiteMatter_mean'],
+        'fmriprep': 'white_matter'
+    },
+    'aCompCor': {
+        'name': 'aCompCor',
+        'link': 'https://fcp-indi.github.io/docs/user/nuisance.html#acompcor',
+        'C-PAC': ['aCompCorPC', 'aCompCor'],
+        'fmriprep': 'aCompCor_comp_cor_0'
+    },
+    'tCompCor': {
+        'name': 'tCompCor',
+        'link': 'https://fcp-indi.github.io/docs/user/nuisance.html#tcompcor',
+        'C-PAC': ['tCompCorPC', 'tCompCor'],
+        'fmriprep': 'tCompCor_comp_cor_0'
+    },
+    'FD': {
+        'name': 'framewise displacement',
+        'link': 'https://fcp-indi.github.io/docs/user/nuisance.html#'
+                'regression-of-motion-parameters'
+    }
+}
+motion_list = ['FD']
+regressor_list = [
+    'GS',
+    'CSF',
+    'WM',
+    'tCompCor0',
+    'aCompCor0',
+    'aCompCor1',
+    'aCompCor2',
+    'aCompCor3',
+    'aCompCor4'
+]
+software = ["C-PAC", "fmriprep"]
diff --git a/configs/subjects.py b/configs/subjects.py
new file mode 100644
index 0000000..6a6093f
--- /dev/null
+++ b/configs/subjects.py
@@ -0,0 +1,162 @@
+import os
+
+from itertools import chain
+from string import ascii_lowercase
+
+
+def cpac_sub(sub):
+    """
+    Function to convert a string from f"sub-{sub_number}{ses_letter}" to
+    f"sub-{sub_number}_ses-{ses_number}"
+
+    Parameter
+    ---------
+    fmriprep_sub: str
+
+    Returns
+    -------
+    sub: str
+
+    Example
+    -------
+    >>> print(cpac_sub("sub-0025427a"))
+    sub-0025427_ses-1
+    """
+    return(f"{sub[:-1]}_ses-{str(ascii_lowercase.find(sub[-1])+1)}")
+
+
+def fmriprep_sub(sub):
+    """
+    Function to convert a string from f"sub-{sub_number}_ses-{ses_number}" to
+    f"sub-{sub_number}{ses_letter}"
+
+    Parameter
+    ---------
+    sub: str
+
+    Returns
+    -------
+    fmriprep_sub: str
+
+    Example
+    -------
+    >>> print(fmriprep_sub("sub-0025427_ses-1"))
+    sub-0025427a
+    """
+    return(f"{sub.split('_')[0]}{ascii_lowercase[int(sub[-1])-1]}")
+
+
+def generate_subject_list_for_directory(path, old_outputs_software="C-PAC"):
+    """
+    Function to take a path and return a subject list.
+
+    Parameter
+    ---------
+    path: str
+
+    old_outputs_software: str, optional, default="C-PAC"
+
+    Returns
+    -------
+    sub_list: list
+    """
+    output = os.path.join(path, "output")
+    sub_ses_list = list(chain.from_iterable([[
+        d for d in os.listdir(
+            os.path.join(output, o)
+        ) if all([
+            os.path.isdir(os.path.join(output, o, d)),
+            d not in ["log", "logs"]
+        ])
+    ] for o in os.listdir(output)]))
+    return(sessions_together([
+        cpac_sub(s) if s[
+            -1
+        ] in ascii_lowercase else s for s in sub_ses_list
+    ]))
+
+
+def generate_subject_list_for_range(
+    subject_start_stop,
+    session_start_stop=None
+):
+    """
+    Function to create a subject list for a given range. All values are
+    inclusive.
+
+    Parameters
+    ----------
+    subject_start_stop: 2-tuple of integers (start, stop) or list of specific
+    values
+
+    session_start_stop: 2-tuple of integers (start, stop) or list of specific
+    values or None
+
+    Returns
+    -------
+    List of strings
+
+    Example
+    -------
+    >>> generate_subject_list_for_range((25427,25428), (1,2))
+    ['sub-0025427_ses-1', 'sub-0025428_ses-1', 'sub-0025427_ses-2', 'sub-0025428_ses-2']
+    """
+    return([
+        f'sub-00{sub}{ses_string}' for ses_string in ([
+            f'_ses-{ses}' for ses in _expand_range(
+                session_start_stop
+            )
+        ] if session_start_stop else [
+            ''
+        ]) for sub in _expand_range(subject_start_stop)
+    ])
+
+
+def sessions_together(sub_list):
+    """
+    Function to sort by session then by subject
+
+    Parameter
+    ---------
+    sub_list: list of str
+
+    Returns
+    -------
+    sub_list: list of str
+
+    Example
+    -------
+    >>> sub_list = [
+    ...    'sub-0025427_ses-1', 'sub-0025427_ses-2', 'sub-0025428_ses-1'
+    ... ]
+    >>> print(sessions_together(sub_list))
+    ['sub-0025427_ses-1', 'sub-0025428_ses-1', 'sub-0025427_ses-2']
+    """
+    sub_list.sort()
+    sub_list.sort(key=lambda x: x.split("ses-")[-1])
+    return(sub_list)
+
+
+def _expand_range(tuple_or_list):
+    """
+    Function to expand an inclusive tuple to a range or return a literal list
+
+    Parameter
+    ---------
+    tuple_or_list: 2-tuple of integers or list
+
+    Returns
+    -------
+    list
+    """
+    return(
+        list(
+            range(
+                tuple_or_list[0],
+                tuple_or_list[1] + 1) if all([
+                isinstance(tuple_or_list, tuple),
+                len(tuple_or_list)==2,
+                *[isinstance(v, int) for v in tuple_or_list]
+            ]) else tuple_or_list
+        )
+    )
diff --git a/correlation_matrix.py b/correlation_matrix.py
new file mode 100644
index 0000000..c0ea67e
--- /dev/null
+++ b/correlation_matrix.py
@@ -0,0 +1,529 @@
+# coding=utf-8
+import sys
+
+if (sys.version_info < (3, 6)):
+    raise EnvironmentError("This module requires Python 3.6 or newer.")
+
+import argparse
+import glob
+import numpy as np
+import os
+import pandas as pd
+import scipy.io as sio
+
+from afnipy.lib_afni1D import Afni1D
+from itertools import chain
+from scipy.stats import pearsonr
+from tabulate import tabulate
+
+try:
+    from configs.defaults import feature_headers, motion_list, regressor_list, \
+                                 software
+    from configs.subjects import fmriprep_sub, \
+                                 generate_subject_list_for_directory
+    from heatmaps import generate_heatmap, reshape_corrs
+except ModuleNotFoundError:
+    from .configs.defaults import feature_headers, motion_list, regressor_list,\
+                                  software
+    from .configs.subjects import fmriprep_sub, \
+                                  generate_subject_list_for_directory
+    from .heatmaps import generate_heatmap, reshape_corrs
+
+sorted_keys = list(feature_headers.keys())
+sorted_keys.sort(key=str.lower)
+feat_def_table = tabulate(
+    [
+        [
+            key,
+            feature_headers[key].get("name"),
+            feature_headers[key].get("link")
+        ] for key in sorted_keys
+    ],
+    headers=["key", "feature name", "documentation link"]
+)
+del(sorted_keys)
+
+def calc_corr(data1, data2):
+    """
+    Function to calculate Pearson's r between two np.ndarrays or lists
+
+    Parameters
+    ----------
+    data1: np.ndarray or list
+
+    data2: np.ndarray or list
+    """
+    if not any([
+        data1 is None,
+        data2 is None
+    ]):
+        if isinstance(data1, np.ndarray) and data1.shape == data2.shape:
+            return(pearsonr(data1.flatten(), data2.flatten())[0])
+        lens = (len(data1), len(data2))
+        if lens[0]==lens[1]:
+            return(pearsonr(data1, data2)[0])
+        if lens[0]>lens[1]:
+            return(pearsonr(data1[lens[0]-lens[1]:], data2)[0])
+        if lens[0]<lens[1]:
+            return(pearsonr(data1, data2[lens[1]-lens[0]:])[0])
+    return(float(np.nan))
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Create a correlation matrix between two C-PAC output "
+                    "directories.",
+        epilog="The following features currently have available definitions to "
+               "calculate Pearson's \x1B[3mr\x1B[23m between C-PAC and "
+               f"fmriprep:\n\n{feat_def_table}",
+        formatter_class=argparse.RawDescriptionHelpFormatter
+    )
+
+    path_help = ("path to an outputs directory - the "
+                 "folder containing the participant-ID "
+                 "labeled directories")
+
+    parser.add_argument("--old_outputs_path", type=str,
+                        help=path_help, default="fmriprep")
+
+    parser.add_argument("--old_outputs_software", type=str,
+                        choices=software, default="fmriprep",
+                        help="(default: %(default)s)")
+
+    parser.add_argument("--new_outputs_path", type=str,
+                        help=path_help)
+
+    parser.add_argument("--new_outputs_software", type=str,
+                        choices=software, default="C-PAC",
+                        help="(default: %(default)s)")
+
+    parser.add_argument("--save", dest="save", action='store_true',
+                        help="save matrices & heatmap (default)")
+
+    parser.add_argument("--no-save", dest="save", action='store_false',
+                        help="do not save matrices & heatmap")
+
+    parser.set_defaults(save=True)
+
+    parser.add_argument("--subject_list", type=str,
+                        help="(default: subjects in OLD_OUTPUTS_PATH sorted by "
+                             "session, subject ID). TODO: handle path to file")
+
+    parser.add_argument("--session", type=int,
+                        help="limit to a single given session (integer)")
+
+    parser.add_argument("--feature_list", type=str,
+                        default=regressor_list + motion_list,
+                        help="TODO: handle path to file (default: %(default)s)")
+
+    parser.add_argument("num_cores", type=int, \
+                            help="number of cores to use - will calculate " \
+                                 "correlations in parallel if greater than 1")
+
+    parser.add_argument("run_name", type=str, \
+                            help="name for the correlations run")
+
+    args = parser.parse_args()
+
+    subject_list = args.subject_list if (
+        "subject_list" in args and args.subject_list is not None
+    ) else generate_subject_list_for_directory(args.old_outputs_path)
+
+    if "session" in args and args.session is not None:
+        subject_list = [
+            sub for sub in subject_list if sub.endswith(str(args.session))
+        ]
+
+    corrs = Correlation_Matrix(
+        subject_list,
+        args.feature_list,
+        [{
+            "software": args.new_outputs_software,
+            "run_path": args.new_outputs_path if args.new_outputs_path.endswith(
+                "/"
+            ) else f"{args.new_outputs_path}/"
+        }, {
+            "software": args.old_outputs_software,
+            "run_path": args.old_outputs_path if args.old_outputs_path.endswith(
+                "/"
+            ) else f"{args.old_outputs_path}/"
+        }]
+    )
+
+    path_table = corrs.print_filepaths(plaintext=True)
+
+    if args.save:
+        output_dir = os.path.join(
+            os.getcwd(), "correlations_{0}".format(args.run_name)
+        )
+
+        if not os.path.exists(output_dir):
+            try:
+                os.makedirs(output_dir)
+            except:
+                err = ("\n\n[!] Could not create the output directory for the "
+                       "correlations. Do you have write permissions?\n "
+                       f"Attempted output directory: {output_dir}\n\n")
+                raise Exception(err)
+
+        path_table.to_csv(os.path.join(output_dir, "filepaths.csv"))
+        sio.savemat(
+            os.path.join(output_dir, "corrs.mat"), {'corrs':corrs.corrs}
+        )
+
+    generate_heatmap(
+        reshape_corrs(corrs.corrs),
+        args.feature_list,
+        subject_list,
+        save_path=os.path.join(
+            output_dir, "heatmap.png"
+        ) if args.save else args.save,
+        title=f"{args.new_outputs_software} "
+        f"{args.new_outputs_path.split('/')[-1]} vs "
+        f"{args.old_outputs_software} {args.old_outputs_path.split('/')[-1]}"
+    )
+
+
+class Subject_Session_Feature:
+    """
+    A class for (subject × session) × feature data
+    """
+    def __init__(self, subject, feature, runs):
+        """
+        Parameters
+        ----------
+        subject: str
+            (subject × session)
+
+        feature: str
+
+        runs: list of dicts
+            [{"software": str, "run_path": str}]
+        """
+        if "_" in subject:
+            self.subject, self.session = subject.split("_", 1)
+        else:
+            self.subject = subject
+            self.session = None
+        self.feature = feature
+        self.paths = (
+            self.get_paths(
+                self.subject,
+                self.feature,
+                runs[0]["run_path"],
+                runs[0]["software"],
+                self.session
+            ),
+            self.get_paths(
+                self.subject,
+                self.feature,
+                runs[1]["run_path"],
+                runs[1]["software"],
+                self.session
+            )
+        )
+        self.data = (
+            self.read_feature(
+                self.paths[0],
+                self.feature,
+                runs[0]["software"]
+            ),
+            self.read_feature(
+                self.paths[1],
+                self.feature,
+                runs[1]["software"]
+            )
+        )
+        if self.data[0] is not None:
+            print(f"{runs[0]['software']} {self.feature}: {len(self.data[0])}")
+        if self.data[1] is not None:
+            print(f"{runs[1]['software']} {self.feature}: {len(self.data[1])}")
+
+    def get_paths(self, subject, feature, run_path, software="C-PAC",
+        session=None):
+        """
+        Method to find a path to specific outputs
+
+        Parameters
+        ----------
+        subject: str or int
+
+        feature: str
+
+        run_path: str
+
+        software: str
+
+        session: str, int or None
+
+        Returns
+        -------
+        paths: list of str
+        """
+        paths = []
+        if software.lower() in ["cpac", "c-pac"]:
+            subject = str(subject)
+            session = f"*{str(session)}*" if session else ""
+            if feature in regressor_list:
+                paths = glob.glob(
+                    f"{run_path}working/"
+                    f"resting_preproc_*{subject}{session}/"
+                    "nuisance_*_0/_*/_*/"
+                    f"{get_feature_label(feature, 'C-PAC')[1][:-1]}*/"
+                    "*1D"
+                ) if "compcor" in feature.lower(
+                ) else list(chain.from_iterable([
+                    glob.glob(
+                        f"{run_path}working/"
+                        f"resting_preproc_*{subject}{session}/"
+                        "nuisance_*_0/_*/*/build*/*1D"
+                    ),
+                    glob.glob(
+                        f"{run_path}working/"
+                        f"resting_preproc_*{subject}{session}/"
+                        "nuisance_*_0/_*/_*/"
+                        f"{get_feature_label(feature, 'C-PAC')[1]}/"
+                        "roi_stats.csv"
+                    )
+                ]))
+            elif feature in motion_list:
+                # frame wise displacement power
+                paths = glob.glob(
+                    f"{run_path}output/*/*{subject}{session}"
+                    "/frame_wise_displacement_power/*/*"
+                )
+        elif software.lower()=="fmriprep":
+            fmriprep_subject = fmriprep_sub("_".join([subject, session]))
+            if feature in regressor_list:
+                paths = [
+                    f"{run_path}output/fmriprep/{fmriprep_subject}/func/"
+                    f"{fmriprep_subject}_task-rest_run-1"
+                    "_desc-confounds_regressors.tsv"
+                ]
+            elif feature in motion_list:
+                paths = [
+                    f"{run_path}working/fmriprep_wf/"
+                    f"single_subject_{fmriprep_subject[4:]}_wf/"
+                    "func_preproc_task_rest_run_1_wf/"
+                    "bold_confounds_wf/fdisp/fd_power_2012.txt"
+                ]
+        return(paths if len(paths) else [])
+
+    def read_feature(self, files, feature, software="C-PAC"):
+        """
+        Method to read a feature from a given file
+
+        Parameters
+        ----------
+        files: list of str
+            paths to files
+
+        feature: str
+
+        software: str
+
+        Returns
+        -------
+        feature: np.ndarray or list or None
+        """
+        if not files:
+            return(None)
+        software = "C-PAC" if software.lower() in [
+            "c-pac",
+            "cpac"
+        ] else software.lower()
+
+        feature_label = get_feature_label(feature, software)
+
+        if software=="C-PAC":
+            for file in files:
+                if file.endswith(".1D"):
+                    data = Afni1D(file)
+                    if "compcor" in file.lower():
+                        return(data.mat[int(feature_label[1][-1])][1:])
+                    header = data.header[-1] if len(data.header) else ""
+                    header_list = header.split('\t')
+                    if isinstance(feature_label, list):
+                        for fl in feature_label:
+                            if(fl in header_list):
+                                return(data.mat[header_list.index(fl)])
+                    else:
+                        return(
+                            data.mat[header_list.index(feature_label)] if (
+                                feature_label in header_list
+                            ) else data.mat[0][1:] if (
+                                len(data.mat[:])==1
+                            ) else ([None] * len(data.mat[0][1:]))
+                        )
+                elif file.endswith('.csv'):
+                    return(list(pd.read_csv(
+                        file,
+                        sep="\t"
+                    )["Sub-brick"][1:].dropna().astype(float).values))
+
+        elif software=="fmriprep":
+            for file in files:
+                if file.endswith(".tsv"):
+                    data = pd.read_csv(file, sep='\t')
+                    if feature_label in data.columns:
+                        return(data[feature_label])
+                elif file.endswith(".txt"):
+                    with open(file) as f:
+                        return([
+                            float(x) for x in [
+                                x.strip() for x in f.readlines()
+                            ][1:]
+                        ])
+
+        return(None)
+
+class Correlation_Matrix:
+    """
+    A class for (subject × session) × feature correlation matrices
+    """
+    def __init__(self, subject_sessions, features, runs):
+        """
+        Parameters
+        ----------
+        subject_sessions: list of strings
+            ["subject_session", ...]
+
+        features: list of strings
+            ["feature", ...]
+
+        runs: list of dicts
+            [{"software": str, "run_path": str}]
+        """
+        self.subjects = subject_sessions
+        self.features = features
+        self.runs = runs
+        self.data = {
+            subject: {
+                feature: Subject_Session_Feature(
+                    subject, feature, runs
+                ) for feature in features
+            } for subject in subject_sessions
+        }
+        self.corrs = np.zeros((len(subject_sessions), len(features)))
+        self.run_pearsonsr()
+
+    def print_filepaths(self, plaintext=False):
+        """
+        Function to print a table
+        """
+        columns = ["\n".join([
+            self.runs[i]["software"], self.runs[i]["run_path"]
+        ]) for i in range(2)]
+        plaintext_columns = ["\n".join([
+            self.runs[i]["software"], wrap(self.runs[i]["run_path"])
+        ]) for i in range(2)]
+        path_table = pd.DataFrame(
+            [[
+                "Not found" if not
+                self.data[sub][feat].paths[i] else (
+                    self._join_paths(self.data[sub][feat].paths, i)
+                ) for i in range(2)
+            ] for sub in self.data for feat in self.data[sub]],
+            columns=columns,
+            index=[
+                f"{sub} {feat}" for sub in self.subjects for
+                feat in self.features
+            ]
+        )
+        if plaintext:
+            plaintext_path_table = pd.DataFrame(
+                [[
+                    f"\u001b[3m\u001b[31mNot found\u001b[0m{' '*13}" if not
+                    self.data[sub][feat].paths[i] else wrap(
+                        self._join_paths(self.data[sub][feat].paths, i)
+                    ) for i in range(2)
+                ] for sub in self.data for feat in self.data[sub]],
+                columns=plaintext_columns,
+                index=[
+                    f"{sub} {feat}" for sub in self.subjects for
+                    feat in self.features
+                ]
+            )
+            print(tabulate(
+                plaintext_path_table,
+                headers=plaintext_columns
+            ))
+        else:
+            stored_options = (
+                pd.options.display.max_rows,
+                pd.options.display.max_colwidth
+            )
+            pd.options.display.max_rows = 999
+            pd.options.display.max_colwidth = 1000
+            try:
+                from IPython.display import display
+                display(path_table)
+            except:
+                print(path_table)
+            (
+                pd.options.display.max_rows,
+                pd.options.display.max_colwidth
+            ) = stored_options
+            del stored_options
+        return(path_table)
+
+    def run_correlation(self, subject, feature, data1, data2):
+        """
+        A method to fill a cell in a correlation matrix with Pearson's r
+
+        Parameters
+        ----------
+        subject: int
+            subject index
+
+        feature: int
+            feature index
+
+        data1: np.ndarray or list
+
+        data2: np.ndarray or list
+        """
+        corr = calc_corr(data1, data2)
+        print(
+            f"Running subject: {subject} {feature} "
+            f"correlation score: {str(corr)}"
+        )
+        self.corrs[subject][feature] = round(corr, 3)
+
+    def run_pearsonsr(self):
+        for i, subject in enumerate(self.data):
+            for j, feature in enumerate(self.data[subject]):
+                self.run_correlation(i, j, *self.data[subject][feature].data)
+
+    def _join_paths(self, data_paths, index):
+        return(
+            "\n".join([
+                data_path.replace(
+                    self.runs[index]["run_path"], "", 1
+                ) if data_path.startswith(
+                    self.runs[index]["run_path"]
+                ) else data_path for data_path in data_paths[index]
+            ])
+        )
+
+
+def get_feature_label(feature, software):
+    return(feature_headers.get(feature, {}).get(software, "") if (
+        "CompCor" not in feature
+    ) else [
+        f"{feature[:-1]}{feature[-1]}",
+        f"{feature[:-1]}_{feature[-1]}"
+    ] if (
+        software=="C-PAC"
+    ) else f"{feature[0]}_comp_cor_0{feature[-1]}" if (
+        software=="fmriprep"
+    ) else "")
+
+
+def wrap(string, at=25):
+    return('\n'.join([
+        string[i:i+at] for i in range(0, len(string), at)
+    ]))
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/heatmaps.py b/heatmaps.py
new file mode 100644
index 0000000..2b2f0d2
--- /dev/null
+++ b/heatmaps.py
@@ -0,0 +1,306 @@
+import argparse
+import matplotlib as mpl
+import numpy as np
+import os
+import pandas as pd
+import sys
+import yaml
+
+from matplotlib import gridspec as GS
+from matplotlib import pyplot as plt
+from scipy import io as sio
+from warnings import filterwarnings
+
+try:
+    from configs import defaults
+    from configs.subjects import generate_subject_list_for_range
+except ModuleNotFoundError:
+    from .configs import defaults
+    from .configs.subjects import generate_subject_list_for_range
+
+filterwarnings(
+    "ignore",
+    "Warning: converting a masked element to nan"
+)
+
+def annotate_heatmap(im, data=None, valfmt="{x:.2f}",
+                     textcolors=["black", "white"],
+                     threshold=None, **textkw):
+    """
+    A function to annotate a heatmap.
+
+    Parameters
+    ----------
+    im
+        The AxesImage to be labeled.
+    data
+        Data used to annotate.  If None, the image's data is used.  Optional.
+    valfmt
+        The format of the annotations inside the heatmap.  This should either
+        use the string format method, e.g. "$ {x:.2f}", or be a
+        `matplotlib.ticker.Formatter`.  Optional.
+    textcolors
+        A list or array of two color specifications.  The first is used for
+        values below a threshold, the second for those above.  Optional.
+    threshold
+        Value in data units according to which the colors from textcolors are
+        applied.  If None (the default) uses the middle of the colormap as
+        separation.  Optional.
+    **kwargs
+        All other arguments are forwarded to each call to `text` used to create
+        the text labels.
+    """
+    import matplotlib
+
+    if not isinstance(data, (list, np.ndarray)):
+        data = im.get_array()
+
+    # Normalize the threshold to the images color range.
+    if threshold is not None:
+        threshold = im.norm(threshold)
+    else:
+        threshold = im.norm(data.max())/2.
+
+    # Set default alignment to center, but allow it to be
+    # overwritten by textkw.
+    kw = dict(horizontalalignment="center",
+              verticalalignment="center")
+    kw.update(textkw)
+
+    # Get the formatter in case a string is supplied
+    if isinstance(valfmt, str):
+        valfmt = matplotlib.ticker.StrMethodFormatter(valfmt)
+
+    # Loop over the data and create a `Text` for each "pixel".
+    # Change the text's color depending on the data.
+    texts = []
+    for i in range(data.shape[0]):
+        for j in range(data.shape[1]):
+            kw.update(color=textcolors[int(im.norm(data[i, j]) < threshold)])
+            text = im.axes.text(j, i, valfmt(data[i, j], None), fontsize=15, **kw)
+            texts.append(text)
+
+    return texts
+
+
+def generate_heatmap(corrs, var_list, sub_list, save_path=None, title=None):
+    """
+    Function to generate a heatmap.
+
+    Parameters
+    ----------
+    corrs: numpy ndarray with shape (number of features, number of subject_sessions)
+        This matrix contains the values to plot
+    var_list: list of strings
+        The labels, in order, of the features (rows)
+    sub_list: list of strings
+        The labels, in order, of the subject_sessions (columns)
+    save_path: string or falsy
+        The path to save the file to, or a falsy value to display in IPython
+    title: str
+        String to use as plot title. Optional.
+
+    Returns
+    -------
+    None
+    """
+    fig, ax = plt.subplots(figsize = (50, 15))
+    im, cbar = heatmap(
+        corrs, var_list, sub_list, ax=ax, vmin=0, vmax=1,
+        cbarlabel="correlation score"
+    )
+    texts = annotate_heatmap(im)
+    if title:
+        plt.title(
+            label=title,
+            fontdict={
+                'fontsize': max(24, len(sub_list)*0.75),
+                'fontweight' : 'bold',
+            }
+        )
+    fig.tight_layout()
+
+    if save_path:
+        plt.savefig(save_path, bbox_inches="tight")
+    else:
+        try:
+            from IPython.display import display
+            plt.show()
+        except:
+            print("No save path or display configured")
+
+
+def heatmap(data, row_labels, col_labels, ax=None,
+            cbar_kw={}, cbarlabel="", **kwargs):
+
+    """
+    Create a heatmap from a numpy array and two lists of labels.
+
+    Parameters
+    ----------
+    data
+        A 2D numpy array of shape (N, M).
+    row_labels
+        A list or array of length N with the labels for the rows.
+    col_labels
+        A list or array of length M with the labels for the columns.
+    ax
+        A `matplotlib.axes.Axes` instance to which the heatmap is plotted.  If
+        not provided, use current axes or create a new one.  Optional.
+    cbar_kw
+        A dictionary with arguments to `matplotlib.Figure.colorbar`.  Optional.
+    cbarlabel
+        The label for the colorbar.  Optional.
+    **kwargs
+        All other arguments are forwarded to `imshow`.
+    """
+
+    if not ax:
+        ax = plt.gca()
+
+    # Plot the heatmap
+    im = ax.imshow(data, **kwargs)
+
+    # Create colorbar
+    cbar = ax.figure.colorbar(
+        im, ax=ax, values=None, boundaries=None, fraction=0.03, pad=0.03
+    )
+    cbar.ax.set_ylabel(cbarlabel, fontsize=20, rotation=-90, va="bottom")
+
+    # We want to show all ticks...
+    ax.set_xticks(np.arange(data.shape[1]))
+    ax.set_yticks(np.arange(data.shape[0]))
+
+    # ... and label them with the respective list entries.
+    ax.set_xticklabels(col_labels)
+    ax.set_yticklabels(row_labels)
+
+    # Let the horizontal axes labeling appear on top.
+    ax.tick_params(top=True, bottom=False,
+                   labeltop=True, labelbottom=False)
+
+    # Rotate the tick labels and set their alignment.
+    plt.setp(ax.get_xticklabels(), fontsize=20, rotation=-30, ha="right",
+             rotation_mode="anchor")
+    plt.setp(ax.get_yticklabels(), fontsize=20, ha="right",
+             rotation_mode="anchor")
+
+    # Turn spines off and create white grid.
+    for edge, spine in list(ax.spines.items()):
+        spine.set_visible(False)
+
+    ax.set_xticks(np.arange(data.shape[1]+1)-.5, minor=True)
+    ax.set_yticks(np.arange(data.shape[0]+1)-.5, minor=True)
+    ax.grid(which="minor", color="w", linestyle='-', linewidth=3)
+    ax.tick_params(which="minor", bottom=False, left=False)
+
+    return im, cbar
+
+
+def reshape_corrs(correlation_matrix):
+    """
+    Function to reshape a given correlation matrix file to the shape expected by matplotlib.
+
+    Parameter
+    ---------
+    correlation_matrix: str or np.ndarray
+        path to matrix file or matrix
+
+    Returns
+    -------
+    matrix: np.ndarray
+        numpy n-dimensional array in the shape of the heatmap
+        [features, subject_sessions]
+    """
+    return(
+        abs(np.transpose(
+            sio.loadmat(
+                correlation_matrix_path
+            )['corrs'] if isinstance(
+                correlation_matrix, str
+            ) else correlation_matrix
+        ))
+    )
+
+
+def parse_args(args):
+    parser = argparse.ArgumentParser(
+        description="generate heatmaps"
+    )
+
+    parser.add_argument(
+        'config',
+        help='path to a YAML configuration file specifying the data, '
+             'features and participants to plot'
+    )
+
+    parser.add_argument(
+        '-o', '--output',
+        dest='save_path',
+        help='path to save heatmap to',
+        required=False
+    )
+
+    parsed = vars(parser.parse_args(args[1:] if len(args)>1 else args))
+    return(parsed.pop('config'), parsed)
+
+
+def main(config_path, save_path=None):
+    with open(config_path, 'r') as config_file:
+        config_settings = yaml.safe_load(config_file)
+    generate_heatmap(
+        reshape_corrs(
+            config_settings['correlation_matrix']
+        ) if 'correlation_matrix' in config_settings else
+        defaults.correlation_matrix,
+        var_list=config_settings[
+            'var_list'
+        ] if 'var_list' in config_settings else (
+            config_settings.get(
+                'regressor_list', []
+            ) + config_settings.get(
+                'motion_list',
+                []
+            )
+        ) if any([
+            l in config_settings for l in [
+                'regressor_list',
+                'motion_list'
+            ]
+        ]) else (
+            defaults.regressor_list + defaults.motion_list
+        ),
+        sub_list=generate_subject_list_for_range(
+            (
+                config_settings['subjects']['start'],
+                config_settings['subjects']['stop']
+            ) if all([
+                'subjects' in config_settings,
+                'start' in config_settings['subjects'],
+                'stop' in config_settings['subjects']
+            ]) else config_settings[
+                'subjects'
+            ] if 'subjects' in config_settings else (
+                defaults.subjects['start'],
+                defaults.subjects['stop']
+            ), (
+                config_settings['sessions']['start'],
+                config_settings['sessions']['stop']
+            ) if all([
+                'sessions' in config_settings,
+                'start' in config_settings['sessions'],
+                'stop' in config_settings['sessions']
+            ]) else config_settings[
+                'sessions'
+            ] if 'sessions' in config_settings else (
+                defaults.sessions['start'],
+                defaults.sessions['stop']
+            )
+        ),
+        save_path=save_path
+    )
+
+
+if __name__ == "__main__":
+    parsed = parse_args(sys.argv)
+    main(parsed[0], **parsed[1])
diff --git a/pytest.ini b/pytest.ini
new file mode 100644
index 0000000..c1fab1c
--- /dev/null
+++ b/pytest.ini
@@ -0,0 +1,2 @@
+[pytest]
+addopts = --continue-on-collection-errors --doctest-ignore-import-errors --doctest-modules --ignore-glob="corr_two_*.py"
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..97440f2
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,11 @@
+git+https://github.com/afni/afni.git@master#egg=afnipy&subdirectory=src/python_scripts
+coverage
+git_python
+matplotlib
+nibabel
+numpy
+pandas
+pytest
+pyyaml
+scipy
+tabulate
\ No newline at end of file
diff --git a/test_cpac_correlations_wf.py b/test_cpac_correlations_wf.py
index 6c1d0e7..19dfd81 100755
--- a/test_cpac_correlations_wf.py
+++ b/test_cpac_correlations_wf.py
@@ -2,7 +2,10 @@
 
 def test_create_unique_file_dict():
 
-    from cpac_correlations_wf import create_unique_file_dict
+    try:
+        from cpac_correlations_wf import create_unique_file_dict
+    except ModuleNotFoundError:
+        from .cpac_correlations_wf import create_unique_file_dict
 
     filepaths = [
       "/path/sub001/centrality_outputs/_scan_rest_1/degree_centrality_weighted.nii.gz",
@@ -39,7 +42,10 @@ def test_create_unique_file_dict():
 
 def test_create_unique_file_dict_with_replacements():
 
-    from cpac_correlations_wf import create_unique_file_dict
+    try:
+        from cpac_correlations_wf import create_unique_file_dict
+    except ModuleNotFoundError:
+        from .cpac_correlations_wf import create_unique_file_dict
 
     filepaths = [
       "/path/sub001_site1/centrality_outputs/_scan_rest_1/degree_centrality_weighted.nii.gz",
@@ -79,7 +85,10 @@ def test_create_unique_file_dict_with_replacements():
 
 def test_match_filepaths():
 
-    from cpac_correlations_wf import match_filepaths
+    try:
+        from cpac_correlations_wf import match_filepaths
+    except ModuleNotFoundError:
+        from .cpac_correlations_wf import match_filepaths
 
     old_files_dict = {
       'alff_to_standard_smooth':