"from_tri_2_sym",

"from_sym_2_tri",

"sumexp_stable",

"concatenate_not_none",

"cov2corr",

"ReadDesign",

"gen_design",

"usable_cpu_count",

"center_mass_exp",

"phase_randomize",

"ecdf",

"p_from_null"

"""convert a upper triangular matrix in 1D format

symm = np.zeros((dim, dim))

symm[np.triu_indices(dim)] = tri

"""convert a 2D symmetric matrix to an upper

inds = np.triu_indices_from(symm)

tri = symm[inds]

"""Compute the sum of exponents for a list of samples

max_value = data.max(axis=0)

result_exp = np.exp(data - max_value)

result_sum = np.sum(result_exp, axis=0)

"""Construct a numpy array by stacking not-None arrays in a list

# 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)

"""Calculate the correlation matrix based on a

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]

"""A class which has the ability of reading in design matrix in .1D file,

def __init__(self, fname=None, include_orth=True, include_pols=True):

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

# The two flags above dictates whether columns corresponding to

# baseline drift modeled by polynomial functions of time and

# columns corresponding to other orthogonal signals (usually motion)

# are included in nuisance regressors.

self.cols_task = np.where(self.column_types == 1)[0]

self.design_task = self.design[:, self.cols_task]

if np.ndim(self.design_task) == 1:

self.design_task = self.design_task[:, None]

# part of the design matrix related to task conditions.

self.n_TR = np.size(self.design_task, axis=0)

self.cols_nuisance = np.array([])

if self.include_orth:

self.cols_nuisance = np.int0(

np.sort(np.append(self.cols_nuisance,

np.where(self.column_types == 0)[0])))

if self.include_pols:

self.cols_nuisance = np.int0(

np.sort(np.append(self.cols_nuisance,

np.where(self.column_types == -1)[0])))

if np.size(self.cols_nuisance) > 0:

self.reg_nuisance = self.design[:, self.cols_nuisance]

if np.ndim(self.reg_nuisance) == 1:

self.reg_nuisance = self.reg_nuisance[:, None]

else:

self.reg_nuisance = None

# Nuisance regressors for motion, baseline, etc.

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(r'\..', 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(r'\..', 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:

temp_res=0.01,

hrf_para={'response_delay': 6, 'undershoot_delay': 12,

'response_dispersion': 0.9,

'undershoot_dispersion': 0.9,

'undershoot_scale': 0.035}):

""" Generate design matrix based on a list of names of stimulus

if np.ndim(scan_duration) == 0:

scan_duration = [scan_duration]

scan_duration = np.array(scan_duration)

assert np.all(scan_duration > TR), \

'scan duration should be longer than a TR'

if type(stimtime_files) is str:

stimtime_files = [stimtime_files]

assert TR > 0, 'TR should be positive'

assert style == 'FSL' or style == 'AFNI', 'style can only be FSL or AFNI'

n_C = len(stimtime_files) # number of conditions

n_S = np.size(scan_duration) # number of scans

if n_S > 1:

design = [np.empty([int(np.round(duration / TR)), n_C])

for duration in scan_duration]

else:

design = [np.empty([int(np.round(scan_duration / TR)), n_C])]

scan_onoff = np.insert(np.cumsum(scan_duration), 0, 0)

if style == 'FSL':

design_info = _read_stimtime_FSL(stimtime_files, n_C, n_S, scan_onoff)

elif style == 'AFNI':

design_info = _read_stimtime_AFNI(stimtime_files, n_C, n_S, scan_onoff)

response_delay = hrf_para['response_delay']

undershoot_delay = hrf_para['undershoot_delay']

response_disp = hrf_para['response_dispersion']

undershoot_disp = hrf_para['undershoot_dispersion']

undershoot_scale = hrf_para['undershoot_scale']

# generate design matrix

for i_s in range(n_S):

for i_c in range(n_C):

if len(design_info[i_s][i_c]['onset']) > 0:

stimfunction = generate_stimfunction(

onsets=design_info[i_s][i_c]['onset'],

event_durations=design_info[i_s][i_c]['duration'],

total_time=scan_duration[i_s],

weights=design_info[i_s][i_c]['weight'],

temporal_resolution=1.0/temp_res)

hrf = _double_gamma_hrf(response_delay=response_delay,

undershoot_delay=undershoot_delay,

response_dispersion=response_disp,

undershoot_dispersion=undershoot_disp,

undershoot_scale=undershoot_scale,

temporal_resolution=1.0/temp_res)

design[i_s][:, i_c] = convolve_hrf(

stimfunction, TR, hrf_type=hrf, scale_function=False,

temporal_resolution=1.0 / temp_res).transpose() * temp_res

else:

design[i_s][:, i_c] = 0.0

# We multiply the resulting design matrix with

# the temporal resolution to normalize it.

# We do not use the internal normalization

# in double_gamma_hrf because it does not guarantee

# normalizing with the same constant.

""" Utility called by gen_design. It reads in one or more

design_info = [[{'onset': [], 'duration': [], 'weight': []}

for i_c in range(n_C)] for i_s in range(n_S)]

# Read stimulus timing files

for i_c in range(n_C):

with open(stimtime_files[i_c]) as f:

for line in f.readlines():

tmp = line.strip().split()

i_s = np.where(

np.logical_and(scan_onoff[:-1] <= float(tmp[0]),

scan_onoff[1:] > float(tmp[0])))[0]

if len(i_s) == 1:

i_s = i_s[0]

design_info[i_s][i_c]['onset'].append(float(tmp[0])

- scan_onoff[i_s])

if len(tmp) >= 2:

design_info[i_s][i_c]['duration'].append(float(tmp[1]))

else:

design_info[i_s][i_c]['duration'].append(1.0)

if len(tmp) >= 3:

design_info[i_s][i_c]['weight'].append(float(tmp[2]))

else:

design_info[i_s][i_c]['weight'].append(1.0)

""" Utility called by gen_design. It reads in one or more stimulus timing

design_info = [[{'onset': [], 'duration': [], 'weight': []}

for i_c in range(n_C)] for i_s in range(n_S)]

# Read stimulus timing files

for i_c in range(n_C):

with open(stimtime_files[i_c]) as f:

text = f.readlines()

assert len(text) == n_S, \

'Number of lines does not match number of runs!'

for i_s, line in enumerate(text):

events = line.strip().split()

if events[0] == '*':

continue

for event in events:

assert event != '*'

tmp = str.split(event, ':')

if len(tmp) == 2:

duration = float(tmp[1])

else:

duration = 1.0

tmp = str.split(tmp[0], '*')

if len(tmp) == 2:

weight = float(tmp[1])

else:

weight = 1.0

if (float(tmp[0]) >= 0

and float(tmp[0])

< scan_onoff[i_s + 1] - scan_onoff[i_s]):

design_info[i_s][i_c]['onset'].append(float(tmp[0]))

design_info[i_s][i_c]['duration'].append(duration)

design_info[i_s][i_c]['weight'].append(weight)

""" Calculate the center of mass of negative exponential distribution

assert isinstance(interval, tuple), 'interval must be a tuple'

assert len(interval) == 2, 'interval must be length two'

(interval_left, interval_right) = interval

assert interval_left >= 0, 'interval_left must be non-negative'

assert interval_right > interval_left, \

'interval_right must be bigger than interval_left'

assert scale > 0, 'scale must be positive'

if interval_right < np.inf:

return ((interval_left + scale) * np.exp(-interval_left / scale) - (

scale + interval_right) * np.exp(-interval_right / scale)) / (

np.exp(-interval_left / scale) - np.exp(-interval_right / scale))

else:

"""Get number of CPUs usable by the current process.

try:

result = len(os.sched_getaffinity(0))

except AttributeError:

try:

result = len(psutil.Process().cpu_affinity())

except AttributeError:

result = os.cpu_count()

"""Randomly shift signal phases

random_state = check_random_state(random_state)

F = fft(D, axis=1)

if D.shape[1] % 2 == 0:

pos_freq = np.arange(1, D.shape[1] // 2)

neg_freq = np.arange(D.shape[1] - 1, D.shape[1] // 2, -1)

else:

pos_freq = np.arange(1, (D.shape[1] - 1) // 2 + 1)

neg_freq = np.arange(D.shape[1] - 1, (D.shape[1] - 1) // 2, -1)

shift = random_state.rand(D.shape[0], len(pos_freq),

D.shape[2]) * 2 * math.pi

# Shift pos and neg frequencies symmetrically, to keep signal real

F[:, pos_freq, :] *= np.exp(1j * shift)

F[:, neg_freq, :] *= np.exp(-1j * shift)

"""Empirical cumulative distribution function

xp = np.sort(x)

yp = np.arange(len(xp) + 1) / len(xp)

def ecdf_fun(q):

return yp[np.searchsorted(xp, q, side="right")]

max_null_input=None, min_null_input=None):

"""Compute p value of true result from null distribution

if (min_null_input is None) or (max_null_input is None):

real_data = X[..., 0]

leading_dims = tuple([int(d) for d in np.arange(X.ndim - 1)])

# Compute maximum/minimum in each null dataset

max_null = np.max(X[..., 1:], axis=leading_dims)

min_null = np.min(X[..., 1:], axis=leading_dims)

else:

real_data = X

# maximum & minimum in each null dataset should be provided as input

max_null = max_null_input

min_null = min_null_input

# Compute where the true values fall on the null distribution

max_null_ecdf = ecdf(max_null)

if two_sided:

min_null_ecdf = ecdf(min_null)

p = 2 * np.minimum(1 - max_null_ecdf(real_data),

min_null_ecdf(real_data))

p = np.minimum(p, 1)

else:

p = 1 - max_null_ecdf(real_data)