https://github.com/cvnlab/GLMsingle
Tip revision: 3a1a580eae6bc7e221cbe01101be3c669687dc6a authored by Kendrick Kay on 01 November 2022, 22:14:32 UTC
Merge branch 'main' of https://github.com/kendrickkay/GLMsingle into main
Merge branch 'main' of https://github.com/kendrickkay/GLMsingle into main
Tip revision: 3a1a580
example_singletrial.py
import os
import glob
import numpy as np
import nibabel as nib
import pandas as pd
from glmsingle.design.make_design_matrix import make_design
from glmsingle.glmsingle import GLM_single
import time
sub = 2
ses = 1
stimdur = 0.5
tr = 2
proj_path = os.path.join(
'/home',
'adf',
'charesti',
'data',
'arsa-fmri',
'BIDS')
data_path = os.path.join(
proj_path,
'derivatives',
'fmriprep',
'sub-{}',
'ses-{}',
'func')
design_path = os.path.join(
proj_path,
'sub-{}',
'ses-{}',
'func')
runs = glob.glob(
os.path.join(data_path.format(sub, ses), '*preproc*nii.gz'))
runs.sort()
runs = runs[:-1]
eventfs = glob.glob(
os.path.join(design_path.format(sub, ses), '*events.tsv'))
eventfs.sort()
runs = runs[:3]
eventfs = eventfs[:3]
data = []
design = []
for i, (run, eventf) in enumerate(zip(runs, eventfs)):
print(f'run {i}')
y = nib.load(run).get_fdata().astype(np.float32)
dims = y.shape
# y = np.moveaxis(y, -1, 0)
# y = y.reshape([y.shape[0], -1])
n_volumes = y.shape[-1]
# Load onsets and item presented
onsets = pd.read_csv(eventf, sep='\t')["onset"].values
items = pd.read_csv(eventf, sep='\t')["stimnumber"].values
n_events = len(onsets)
# Create design matrix
events = pd.DataFrame()
events["duration"] = [stimdur] * n_events
events["onset"] = np.round(onsets)
events["trial_type"] = items
# pass in the events data frame. the convolving of the HRF now
# happens internally
design.append(
make_design(events, tr, n_volumes)
)
data.append(y)
opt = {'wantlss': 0}
outputdir = 'GLMestimatesingletrialoutputs'
start_time = time.time()
gst = GLM_single(opt)
results = gst.fit(
design,
data,
stimdur,
tr,
outputdir=outputdir)
elapsed_time = time.time() - start_time
print(
'elapsedtime: ',
f'{time.strftime("%H:%M:%S", time.gmtime(elapsed_time))}'
)