File: fiac_util.py

package info (click to toggle)
nipy 0.6.1-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 7,352 kB
  • sloc: python: 39,115; ansic: 30,931; makefile: 210; sh: 93
file content (414 lines) | stat: -rw-r--r-- 13,390 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""Support utilities for FIAC example, mostly path management.

The purpose of separating these is to keep the main example code as readable as
possible and focused on the experimental modeling and analysis, rather than on
local file management issues.

Requires matplotlib
"""

#-----------------------------------------------------------------------------
# Imports
#-----------------------------------------------------------------------------

# Stdlib
import csv
import os
from io import StringIO  # Python 3
from os import listdir, makedirs
from os.path import abspath, exists, isdir, splitext
from os.path import join as pjoin

# Third party
import numpy as np
import pandas as pd

# From NIPY
from nipy.io.api import load_image


def csv2rec(fname):
    return pd.read_csv(fname).to_records()


def rec2csv(recarr, fname):
    pd.DataFrame.from_records(recarr).to_csv(fname, index=None)


#-----------------------------------------------------------------------------
# Globals
#-----------------------------------------------------------------------------

# We assume that there is a directory holding the data and it's local to this
# code.  Users can either keep a copy here or a symlink to the real location on
# disk of the data.
DATADIR = 'fiac_data'

# Sanity check
if not os.path.isdir(DATADIR):
    e=f"The data directory {DATADIR} must exist and contain the FIAC data."
    raise OSError(e)

#-----------------------------------------------------------------------------
# Classes and functions
#-----------------------------------------------------------------------------

# Path management utilities
def load_image_fiac(*path):
    """Return a NIPY image from a set of path components.
    """
    return load_image(pjoin(DATADIR, *path))


def subj_des_con_dirs(design, contrast, nsub=16):
    """Return a list of subject directories with this `design` and `contrast`

    Parameters
    ----------
    design : {'event', 'block'}
    contrast : str
    nsub : int, optional
        total number of subjects

    Returns
    -------
    con_dirs : list
        list of directories matching `design` and `contrast`
    """
    rootdir = DATADIR
    con_dirs = []
    for s in range(nsub):
        f = pjoin(rootdir, "fiac_%02d" % s, design, "fixed", contrast)
        if isdir(f):
            con_dirs.append(f)
    return con_dirs


def path_info_run(subj, run):
    """Construct path information dict for current subject/run.

    Parameters
    ----------
    subj : int
        subject number (0..15 inclusive)
    run : int
        run number (1..4 inclusive).

    Returns
    -------
    path_dict : dict
        a dict with all the necessary path-related keys, including 'rootdir',
        and 'design', where 'design' can have values 'event' or 'block'
        depending on which type of run this was for subject no `subj` and run no
        `run`
    """
    path_dict = {'subj': subj, 'run': run}
    if exists(pjoin(DATADIR, "fiac_%(subj)02d",
                    "block", "initial_%(run)02d.csv") % path_dict):
        path_dict['design'] = 'block'
    else:
        path_dict['design'] = 'event'
    rootdir = pjoin(DATADIR, "fiac_%(subj)02d", "%(design)s") % path_dict
    path_dict['rootdir'] = rootdir
    return path_dict


def path_info_design(subj, design):
    """Construct path information dict for subject and design.

    Parameters
    ----------
    subj : int
        subject number (0..15 inclusive)
    design : {'event', 'block'}
        type of design

    Returns
    -------
    path_dict : dict
        having keys 'rootdir', 'subj', 'design'
    """
    path_dict = {'subj': subj, 'design': design}
    rootdir = pjoin(DATADIR, "fiac_%(subj)02d", "%(design)s") % path_dict
    path_dict['rootdir'] = rootdir
    return path_dict


def results_table(path_dict):
    """ Return precalculated results images for subject info in `path_dict`

    Parameters
    ----------
    path_dict : dict
        containing key 'rootdir'

    Returns
    -------
    rtab : dict
        dict with keys given by run directories for this subject, values being a
        list with filenames of effect and sd images.
    """
    # Which runs correspond to this design type?
    rootdir = path_dict['rootdir']
    runs = filter(lambda f: isdir(pjoin(rootdir, f)),
                  ['results_%02d' % i for i in range(1,5)] )

    # Find out which contrasts have t-statistics,
    # storing the filenames for reading below

    results = {}

    for rundir in runs:
        rundir = pjoin(rootdir, rundir)
        for condir in listdir(rundir):
            for stat in ['sd', 'effect']:
                fname_effect = abspath(pjoin(rundir, condir, 'effect.nii'))
                fname_sd = abspath(pjoin(rundir, condir, 'sd.nii'))
            if exists(fname_effect) and exists(fname_sd):
                results.setdefault(condir, []).append([fname_effect,
                                                       fname_sd])
    return results


def get_experiment_initial(path_dict):
    """Get the record arrays for the experimental/initial designs.

    Parameters
    ----------
    path_dict : dict
        containing key 'rootdir', 'run', 'subj'

    Returns
    -------
    experiment, initial : Two record arrays.

    """
    # The following two lines read in the .csv files
    # and return recarrays, with fields
    # experiment: ['time', 'sentence', 'speaker']
    # initial: ['time', 'initial']

    rootdir = path_dict['rootdir']
    if not exists(pjoin(rootdir, "experiment_%(run)02d.csv") % path_dict):
        e = "can't find design for subject=%(subj)d,run=%(subj)d" % path_dict
        raise OSError(e)

    experiment = csv2rec(pjoin(rootdir, "experiment_%(run)02d.csv") % path_dict)
    initial = csv2rec(pjoin(rootdir, "initial_%(run)02d.csv") % path_dict)

    return experiment, initial


def get_fmri(path_dict):
    """Get the images for a given subject/run.

    Parameters
    ----------
    path_dict : dict
        containing key 'rootdir', 'run'

    Returns
    -------
    fmri : ndarray
    anat : NIPY image
    """
    fmri_im = load_image(
        pjoin("%(rootdir)s/swafunctional_%(run)02d.nii") % path_dict)
    return fmri_im


def ensure_dir(*path):
    """Ensure a directory exists, making it if necessary.

    Returns the full path."""
    dirpath = pjoin(*path)
    if not isdir(dirpath):
        makedirs(dirpath)
    return dirpath


def output_dir(path_dict, tcons, fcons):
    """Get (and make if necessary) directory to write output into.

    Parameters
    ----------
    path_dict : dict
        containing key 'rootdir', 'run'
    tcons : sequence of str
        t contrasts
    fcons : sequence of str
        F contrasts
    """
    rootdir = path_dict['rootdir']
    odir = pjoin(rootdir, "results_%(run)02d" % path_dict)
    ensure_dir(odir)
    for n in tcons:
        ensure_dir(odir,n)
    for n in fcons:
        ensure_dir(odir,n)
    return odir


def test_sanity():
    import nipy.modalities.fmri.fmristat.hrf as fshrf
    from nipy.algorithms.statistics import formula
    from nipy.modalities.fmri import design, hrf
    from nipy.modalities.fmri.fmristat.tests import FIACdesigns
    from nipy.modalities.fmri.fmristat.tests.test_FIAC import matchcol

    """
    Single subject fitting of FIAC model
    """

    # Based on file
    # subj3_evt_fonc1.txt
    # subj3_bloc_fonc3.txt

    for subj, run, design_type in [(3, 1, 'event'), (3, 3, 'block')]:
        nvol = 191
        TR = 2.5
        Tstart = 1.25

        volume_times = np.arange(nvol)*TR + Tstart
        volume_times_rec = formula.make_recarray(volume_times, 't')

        path_dict = {'subj':subj, 'run':run}
        if exists(pjoin(DATADIR, "fiac_%(subj)02d",
                        "block", "initial_%(run)02d.csv") % path_dict):
            path_dict['design'] = 'block'
        else:
            path_dict['design'] = 'event'

        experiment = csv2rec(pjoin(DATADIR, "fiac_%(subj)02d", "%(design)s", "experiment_%(run)02d.csv")
                             % path_dict)
        initial = csv2rec(pjoin(DATADIR, "fiac_%(subj)02d", "%(design)s", "initial_%(run)02d.csv")
                                % path_dict)

        X_exper, cons_exper = design.event_design(experiment,
                                                  volume_times_rec,
                                                  hrfs=fshrf.spectral)
        X_initial, _ = design.event_design(initial,
                                           volume_times_rec,
                                           hrfs=[hrf.glover])
        X, cons = design.stack_designs((X_exper, cons_exper), (X_initial, {}))

        # Get original fmristat design
        Xf = FIACdesigns.fmristat[design_type]
        # Check our new design can be closely matched to the original
        for i in range(X.shape[1]):
            # Columns can be very well correlated negatively or positively
            assert abs(matchcol(X[:,i], Xf)[1]) > 0.999


def rewrite_spec(subj, run, root = "/home/jtaylo/FIAC-HBM2009"):
    """
    Take a FIAC specification file and get two specifications
    (experiment, begin).

    This creates two new .csv files, one for the experimental
    conditions, the other for the "initial" confounding trials that
    are to be modelled out.

    For the block design, the "initial" trials are the first
    trials of each block. For the event designs, the
    "initial" trials are made up of just the first trial.

    """

    if exists(pjoin("%(root)s", "fiac%(subj)d", "subj%(subj)d_evt_fonc%(run)d.txt") % {'root':root, 'subj':subj, 'run':run}):
        designtype = 'evt'
    else:
        designtype = 'bloc'

    # Fix the format of the specification so it is
    # more in the form of a 2-way ANOVA

    eventdict = {1:'SSt_SSp', 2:'SSt_DSp', 3:'DSt_SSp', 4:'DSt_DSp'}
    s = StringIO()
    w = csv.writer(s)
    w.writerow(['time', 'sentence', 'speaker'])

    specfile = pjoin("%(root)s", "fiac%(subj)d", "subj%(subj)d_%(design)s_fonc%(run)d.txt") % {'root':root, 'subj':subj, 'run':run, 'design':designtype}
    d = np.loadtxt(specfile)
    for row in d:
        w.writerow([row[0]] + eventdict[row[1]].split('_'))
    s.seek(0)
    d = csv2rec(s)

    # Now, take care of the 'begin' event
    # This is due to the FIAC design

    if designtype == 'evt':
        b = np.array([(d[0]['time'], 1)], np.dtype([('time', np.float64),
                                                    ('initial', np.int_)]))
        d = d[1:]
    else:
        k = np.equal(np.arange(d.shape[0]) % 6, 0)
        b = np.array([(tt, 1) for tt in d[k]['time']], np.dtype([('time', np.float64),
                                                                 ('initial', np.int_)]))
        d = d[~k]

    designtype = {'bloc':'block', 'evt':'event'}[designtype]

    fname = pjoin(DATADIR, "fiac_%(subj)02d", "%(design)s", "experiment_%(run)02d.csv") % {'root':root, 'subj':subj, 'run':run, 'design':designtype}
    rec2csv(d, fname)
    experiment = csv2rec(fname)

    fname = pjoin(DATADIR, "fiac_%(subj)02d", "%(design)s", "initial_%(run)02d.csv") % {'root':root, 'subj':subj, 'run':run, 'design':designtype}
    rec2csv(b, fname)
    initial = csv2rec(fname)

    return d, b


def compare_results(subj, run, other_root, mask_fname):
    """ Find and compare calculated results images from a previous run

    This script checks that another directory containing results of this same
    analysis are similar in the sense of numpy ``allclose`` within a brain mask.

    Parameters
    ----------
    subj : int
        subject number (0..4, 6..15)
    run : int
        run number (1..4)
    other_root : str
        path to previous run estimation
    mask_fname:
        path to a mask image defining area in which to compare differences
    """
    # Get information for this subject and run
    path_dict = path_info_run(subj, run)
    # Get mask
    msk = load_image(mask_fname).get_fdata().copy().astype(bool)
    # Get results directories for this run
    rootdir = path_dict['rootdir']
    res_dir = pjoin(rootdir, 'results_%02d' % run)
    if not isdir(res_dir):
        return
    for dirpath, dirnames, filenames in os.walk(res_dir):
        for fname in filenames:
            froot, ext = splitext(fname)
            if froot in ('effect', 'sd', 'F', 't'):
                this_fname = pjoin(dirpath, fname)
                other_fname = this_fname.replace(DATADIR, other_root)
                if not exists(other_fname):
                    print(this_fname, 'present but ', other_fname, 'missing')
                    continue
                this_arr = load_image(this_fname).get_fdata()
                other_arr = load_image(other_fname).get_fdata()
                ok = np.allclose(this_arr[msk], other_arr[msk])
                if not ok and froot in ('effect', 'sd', 't'): # Maybe a sign flip
                    ok = np.allclose(this_arr[msk], -other_arr[msk])
                if not ok:
                    print('Difference between', this_fname, other_fname)


def compare_all(other_root, mask_fname):
    """ Run results comparison for all subjects and runs """
    for subj in range(5) + range(6, 16):
        for run in range(1, 5):
            compare_results(subj, run, other_root, mask_fname)