openfmri.py

import os
import glob
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import nibabel
from pypreprocess.external.nistats.design_matrix import (
    make_design_matrix, plot_design_matrix)

# meta data
tr = 2.
drift_model = 'Cosine'
hrf_model = 'Canonical With Derivative'
hfcut = 128.

# load stuff
design_matrices = []
for run_path in sorted(glob.glob(
        "/home/elvis/nilearn_data/ds001/sub001/model/model001/onsets/task*")):
    run_id = os.path.basename(run_path)
    run_func = os.path.join(
        "/home/elvis/nilearn_data/ds001/sub001/BOLD/%s/bold.nii.gz" % run_id)
    run_onset_paths = sorted(glob.glob(
        ("/home/elvis/nilearn_data/ds001/sub001/model/model001/onsets/"
         "%s/*" % run_id)))
    onsets = map(np.loadtxt, run_onset_paths)
    for x, p in zip(onsets[:2], run_onset_paths[:2]):
        print x[:5], p
    # break
    conditions = np.hstack([["cond%03i" % (c + 1)] * len(onsets[c])
                            for c in range(len(run_onset_paths))])
    onsets = np.vstack((onsets))
    onsets *= tr
    run_func = nibabel.load(run_func)
    n_scans = run_func.shape[-1]
    onset, duration, amplitude = onsets.T

    frametimes = np.linspace(0, (n_scans - 1) * tr, n_scans // 2)
    paradigm = pd.DataFrame(dict(name=conditions, onset=onset,
                                 duration=duration, amplitude=amplitude))
    design_matrix = make_design_matrix(frametimes, paradigm,
                                       hrf_model=hrf_model,
                                       drift_model=drift_model,
                                       period_cut=hfcut)
    design_matrices.append(design_matrix)
    plot_design_matrix(design_matrix)
    plt.title(run_id)
plt.show()