work in progress

example.py

@@ -85,4 +85,4 @@ def format_time(time):
 fit_dur = format_time(time.time()-start)
 print('Fit took {}!'.format(fit_dur)
 
-gd.plot_figures()
+gd.plot_figures(spatialdims=dims[:3])

example_singletrial.py

@@ -0,0 +1,97 @@
+import os
+import glob
+import numpy as np
+import nibabel as nib
+import pandas as pd
+from glmdenoise.design.make_design_matrix import make_design
+from glmdenoise.single_trial import GLM_singletrial
+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_singletrial(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))}'
+)

glmdenoise/check_inputs.py

@@ -0,0 +1,90 @@
+import numpy as np
+
+
+def check_inputs(data, design):
+    """
+    check that the data and design meet the required
+    specifications for glm single trial estimates.
+
+    Arguments
+    _________
+
+        data (list): could be x,y,z,t or XYZ,t
+        design (list of runs or single run): design matrix
+
+    Returns:
+    ________
+
+        data (list): flattened XYZ data format
+        design (list): design matrix with a list entry per run
+    """
+    # massage <design> and sanity-check it
+    if type(design) is not list:
+        design = [design]
+
+    numcond = design[0].shape[1]
+    for p in range(len(design)):
+        np.testing.assert_array_equal(
+            np.unique(design),
+            [0, 1],
+            err_msg='<design> must consist of 0s and 1s')
+        condmsg = \
+            'all runs in <design> should have equal number of conditions'
+        np.testing.assert_equal(
+            design[p].shape[1],
+            numcond,
+            err_msg=condmsg)
+        # if the design happened to be a sparse?
+        # design[p] = np.full(design[p])
+
+    # massage <data> and sanity-check it
+    if type(data) is not list:
+        data = [data]
+
+    # make sure it is single
+    for p in range(len(data)):
+        data[p] = data[p].astype(np.float32, copy=False)
+
+    np.testing.assert_equal(
+        np.all(np.isfinite(data[0].flatten())),
+        True,
+        err_msg='We checked the first run and '
+        'found some non-finite values (e.g. NaN, Inf).'
+        'Please fix and re-run.')
+    np.testing.assert_equal(
+        len(design),
+        len(data),
+        err_msg='<design> and <data> should have '
+        'the same number of runs')
+
+    # reshape data in 2D mode.
+    is3d = data[0].ndim > 2  # is this the X Y Z T case?
+    if is3d:
+        xyz = data[0].shape[:3]
+        n_times = data[0].shape[3]
+        for p in range(len(data)):
+            data[p] = np.reshape(
+                data[p],
+                [np.prod(xyz), n_times])
+            # force to XYZ x T for convenience
+    else:
+        xyz = False
+
+    # check number of time points and truncate if necessary
+    for run_i in np.arange(len(data)):
+        if data[run_i].shape[1] > design[run_i].shape[0]:
+            print(
+                f'WARNING: run {run_i} has more time points'
+                'in <data> than <design>. We are truncating'
+                '<data>.\n')
+            data[run_i] = data[run_i][:, np.arange(
+                design.shape[0])]
+
+        if data[run_i].shape[1] < design[run_i].shape[0]:
+            print(
+                f'WARNING: run {run_i} has more time points in <design>'
+                'than <data>. We are truncating <design>.\n')
+            design[run_i] = design[run_i][np.arange(
+                data[run_i].shape[0]), :]
+
+    return data, design, xyz

glmdenoise/defaults.py

@@ -1,11 +1,33 @@
+import time
+import numpy as np
+
 default_params = {
     'numforhrf': 50,
     'hrfthresh': 0.5,
     'hrffitmask': 1,
     'R2thresh': 0,
     'hrfmodel': 'optimise',
     'n_jobs': 1,
-    'n_pcs': 20,
+    'n_pcs': 10,
     'n_boots': 100,
-    'extra_regressors': False
-}
+    'extra_regressors': False,
+    'wantlibrary': 1,
+    'wantglmdenoise': 1,
+    'wantfracridge': 1,
+    'chunklen': 45000,
+    'wantfileoutputs': [1, 1, 1, 1],
+    'wantmemoryoutputs': [1, 1, 1, 1],
+    'wantpercentbold': 1,
+    'wantlss': 0,
+    'brainthresh': [99.0, 0.1],
+    'brainR2': [],
+    'brainexclude': False,
+    'pcR2cutoff': [],
+    'pcR2cutoffmask': 1,
+    'pcstop': 1.05,
+    'fracs': np.linspace(1, 0.05, 20),
+    'wantautoscale': 1,
+    'seed': time.time(),
+    'suppressoutput': 0,
+    'lambda': 0
+}