Repository navigation
Expand file tree
/
Copy pathutils.py
More file actions
331 lines (256 loc) · 10 KB
/
Copy pathutils.py
File metadata and controls
331 lines (256 loc) · 10 KB
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
# Copyright 2016 Intel Corporation, Princeton University
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import numpy as np
import logging
import re
import warnings
import os.path
"""
Some utility functions that can be used by different algorithms
"""
def from_tri_2_sym(tri, dim):
"""convert a upper triangular matrix in 1D format
to 2D symmetric matrix
Parameters
----------
tri: 1D array
Contains elements of upper triangular matrix
dim : int
The dimension of target matrix.
Returns
-------
symm : 2D array
Symmetric matrix in shape=[dim, dim]
"""
symm = np.zeros((dim, dim))
symm[np.triu_indices(dim)] = tri
return symm
def from_sym_2_tri(symm):
"""convert a 2D symmetric matrix to an upper
triangular matrix in 1D format
Parameters
----------
symm : 2D array
Symmetric matrix
Returns
-------
tri: 1D array
Contains elements of upper triangular matrix
"""
inds = np.triu_indices_from(symm)
tri = symm[inds]
return tri
def fast_inv(a):
"""to invert a 2D matrix
Parameters
----------
a : 2D array
Returns
-------
inva: 2D array
inverse of input matrix a
Raises
-------
LinAlgError
If a is singular or not square
"""
if a.ndim != 2:
raise TypeError("Input matrix should be 2D array")
identity = np.identity(a.shape[1], dtype=a.dtype)
inva = None
try:
inva = np.linalg.solve(a, identity)
return inva
except np.linalg.linalg.LinAlgError:
logging.exception('Error from np.linalg.solve')
raise
def sumexp_stable(data):
"""Compute the sum of exponents for a list of samples
Parameters
----------
data : array, shape=[features, samples]
A data array containing samples.
Returns
-------
result_sum : array, shape=[samples,]
The sum of exponents for each sample divided by the exponent
of the maximum feature value in the sample.
max_value : array, shape=[samples,]
The maximum feature value for each sample.
result_exp : array, shape=[features, samples]
The exponent of each element in each sample divided by the exponent
of the maximum feature value in the sample.
..note::
This function is more stable than computing the sum(exp(v)).
It useful for computing the softmax_i(v)=exp(v_i)/sum(exp(v)) function.
"""
max_value = data.max(axis=0)
result_exp = np.exp(data - max_value)
result_sum = np.sum(result_exp, axis=0)
return result_sum, max_value, result_exp
def concatenate_list(l, axis=0):
"""Construct a numpy array by stacking arrays in a list
Parameters
----------
data : list of arrays, arrays have same shape in all but one dimension or
elements are None
The list of arrays to be concatenated.
axis : int, default = 0
Axis for the concatenation
Returns
-------
data_stacked : array
The resulting concatenated array.
"""
# Get the indexes of the arrays in the list
mask = []
for i in range(len(l)):
if l[i] is not None:
mask.append(i)
# Concatenate them
l_stacked = np.concatenate([l[i] for i in mask], axis=axis)
return l_stacked
def cov2corr(cov):
"""Calculate the correlation matrix based on a
covariance matrix
Parameters
----------
cov: 2D array
Returns
-------
corr: 2D array
correlation converted from the covarince matrix
"""
assert cov.ndim == 2, 'covariance matrix should be 2D array'
inv_sd = 1 / np.sqrt(np.diag(cov))
corr = cov * inv_sd[None, :] * inv_sd[:, None]
return corr
class ReadDesign:
""" A class which has the ability of reading in design matrix in .1D file,
generated by AFNI's 3dDeconvolve.
Parameters
----------
fname: string, the address of the file to read.
include_orth: Boollean, whether to include "orthogonal" regressors in
the design matrix which are usually head motion parameters. All
the columns of design matrix are still going to be read in, but
the attribute cols_used will reflect whether these orthogonal
regressors are to be included for furhter analysis.
include_pols: Boollean, whether to include polynomial regressors in
the design matrix which are used to capture slow drift of signals.
This will be reflected in the indices in the attribute cols_used.
Attributes
----------
design: 2d array. The design matrix read in from the csv file.
n_col: number of total columns in the design matrix.
column_types: 1d array. the types of each column in the design matrix.
0 for orthogonal regressors (usually head motion parameters),
-1 for polynomial basis (capturing slow drift of signals),
values > 0 for stimulus conditions
n_basis: scalar. The number of polynomial bases in the designn matrix.
n_stim: scalar. The number of stimulus conditions.
n_orth: scalar. The number of orthogoanal regressors (usually head
motions)
StimLabels: list. The names of each column in the design matrix.
"""
def __init__(self, fname=None, include_orth=False, include_pols=False):
if fname is None:
# fname is the name of the file to read in the design matrix
self.design = np.zeros([0, 0])
self.n_col = 0
# number of columns (conditions) in the design matrix
self.column_types = np.ones(0)
self.n_basis = 0
self.n_stim = 0
self.n_orth = 0
self.StimLabels = []
else:
# isAFNI = re.match(r'.+[.](1D|1d|txt)$', fname)
filename, ext = os.path.splitext(fname)
# We assume all AFNI 1D files have extension of 1D or 1d or txt
if ext in ['.1D', '.1d', '.txt']:
self.read_afni(fname=fname)
self.include_orth = include_orth
self.include_pols = include_pols
self.cols_used = np.where(self.column_types == 1)[0]
if self.include_orth:
self.cols_used = np.sort(
np.append(self.cols_used, np.where(self.column_types == 0)[0]))
if self.include_pols:
self.cols_used = np.sort(np.append(
self.cols_used, np.where(self.column_types == -1)[0]))
self.design_used = self.design[:, self.cols_used]
if not self.include_pols:
# baseline is not included, then we add a column of all 1's
self.design_used = np.insert(self.design_used, 0, 1, axis=1)
self.n_TR = np.size(self.design_used, axis=0)
def read_afni(self, fname):
# Read design file written by AFNI
self.n_basis = 0
self.n_stim = 0
self.n_orth = 0
self.StimLabels = []
self.design = np.loadtxt(fname, ndmin=2)
with open(fname) as f:
all_text = f.read()
find_n_column = re.compile(
r'^#[ ]+ni_type[ ]+=[ ]+"(?P<n_col>\d+)[*]', re.MULTILINE)
n_col_found = find_n_column.search(all_text)
if n_col_found:
self.n_col = int(n_col_found.group('n_col'))
if self.n_col != np.size(self.design, axis=1):
warnings.warn(
'The number of columns in the design matrix'
+ 'does not match the header information')
self.n_col = np.size(self.design, axis=1)
else:
self.n_col = np.size(self.design, axis=1)
self.column_types = np.ones(self.n_col)
# default that all columns are conditions of interest
find_ColumnGroups = re.compile(
r'^#[ ]+ColumnGroups[ ]+=[ ]+"(?P<CGtext>.+)"', re.MULTILINE)
CG_found = find_ColumnGroups.search(all_text)
if CG_found:
CG_text = re.split(',', CG_found.group('CGtext'))
curr_idx = 0
for CG in CG_text:
split_by_at = re.split('@', CG)
if len(split_by_at) == 2:
# the first tells the number of columns in this condition
# the second tells the condition type
n_this_cond = int(split_by_at[0])
self.column_types[curr_idx:curr_idx + n_this_cond] = \
int(split_by_at[1])
curr_idx += n_this_cond
elif len(split_by_at) == 1 and \
not re.search('..', split_by_at[0]):
# Just a number, and not the type like '1..4'
self.column_types[curr_idx] = int(split_by_at[0])
curr_idx += 1
else: # must be a single stimulus condition
split_by_dots = re.split('\..', CG)
n_this_cond = int(split_by_dots[1])
self.column_types[curr_idx:curr_idx + n_this_cond] = 1
curr_idx += n_this_cond
self.n_basis = np.sum(self.column_types == -1)
self.n_stim = np.sum(self.column_types > 0)
self.n_orth = np.sum(self.column_types == 0)
find_StimLabels = re.compile(
r'^#[ ]+StimLabels[ ]+=[ ]+"(?P<SLtext>.+)"', re.MULTILINE)
StimLabels_found = find_StimLabels.search(all_text)
if StimLabels_found:
self.StimLabels = \
re.split(r'[ ;]+', StimLabels_found.group('SLtext'))
else:
self.StimLabels = []