"array_correlation",

"center_mass_exp",

"concatenate_not_none",

"cov2corr",

"from_tri_2_sym",

"from_sym_2_tri",

"gen_design",

"phase_randomize",

"p_from_null",

"ReadDesign",

"sumexp_stable",

"usable_cpu_count",

"""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()

"""Randomize phase of time series across subjects

# Check if input is 2-dimensional

data_ndim = data.ndim

# Get basic shape of data

data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data)

# Random seed to be deterministically re-randomized at each iteration

if isinstance(random_state, np.random.RandomState):

prng = random_state

else:

prng = np.random.RandomState(random_state)

# Get randomized phase shifts

if n_TRs % 2 == 0:

# Why are we indexing from 1 not zero here? n_TRs / -1 long?

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

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

else:

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

neg_freq = np.arange(data.shape[0] - 1,

(data.shape[0] - 1) // 2, -1)

if not voxelwise:

phase_shifts = (prng.rand(len(pos_freq), 1, n_subjects)

* 2 * np.math.pi)

else:

phase_shifts = (prng.rand(len(pos_freq), n_voxels, n_subjects)

* 2 * np.math.pi)

# Fast Fourier transform along time dimension of data

fft_data = fft(data, axis=0)

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

fft_data[pos_freq, :, :] *= np.exp(1j * phase_shifts)

fft_data[neg_freq, :, :] *= np.exp(-1j * phase_shifts)

# Inverse FFT to put data back in time domain

shifted_data = np.real(ifft(fft_data, axis=0))

# Go back to 2-dimensions if input was 2-dimensional

if data_ndim == 2:

shifted_data = shifted_data[:, 0, :]

side='two-sided', exact=False,

axis=None):

"""Compute p-value from null distribution

if side not in ('two-sided', 'left', 'right'):

raise ValueError("The value for 'side' must be either "

"'two-sided', 'left', or 'right', got {0}".

format(side))

n_samples = len(distribution)

logger.info("Assuming {0} resampling iterations".format(n_samples))

if side == 'two-sided':

# Numerator for two-sided test

numerator = np.sum(np.abs(distribution) >= np.abs(observed), axis=axis)

elif side == 'left':

# Numerator for one-sided test in left tail

numerator = np.sum(distribution <= observed, axis=axis)

elif side == 'right':

# Numerator for one-sided test in right tail

numerator = np.sum(distribution >= observed, axis=axis)

# If exact test all possible permutations and do not adjust

if exact:

p = numerator / n_samples

# If not exact test, adjust number of samples to account for

# observed statistic; prevents p-value from being zero

else:

p = (numerator + 1) / (n_samples + 1)

"""Checks response time series input data (e.g., for ISC analysis)

# Convert list input to 3d and check shapes

if type(data) == list:

data_shape = data[0].shape

for i, d in enumerate(data):

if d.shape != data_shape:

raise ValueError("All ndarrays in input list "

"must be the same shape!")

if d.ndim == 1:

data[i] = d[:, np.newaxis]

data = np.dstack(data)

# Convert input ndarray to 3d and check shape

elif isinstance(data, np.ndarray):

if data.ndim == 2:

data = data[:, np.newaxis, :]

elif data.ndim == 3:

pass

else:

raise ValueError("Input ndarray should have 2 "

"or 3 dimensions (got {0})!".format(data.ndim))

# Infer subjects, TRs, voxels and log for user to check

n_TRs, n_voxels, n_subjects = data.shape

logger.info("Assuming {0} subjects with {1} time points "

"and {2} voxel(s) or ROI(s) for ISC analysis.".format(

n_subjects, n_TRs, n_voxels))

"""Column- or row-wise Pearson correlation between two arrays

# Accommodate array-like inputs

if not isinstance(x, np.ndarray):

x = np.asarray(x)

if not isinstance(y, np.ndarray):

y = np.asarray(y)

# Check that inputs are same shape

if x.shape != y.shape:

raise ValueError("Input arrays must be the same shape")

# Transpose if axis=1 requested (to avoid broadcasting

# issues introduced by switching axis in mean and sum)

if axis == 1:

x, y = x.T, y.T

# Center (de-mean) input variables

x_demean = x - np.mean(x, axis=0)

y_demean = y - np.mean(y, axis=0)

# Compute summed product of centered variables

numerator = np.sum(x_demean * y_demean, axis=0)

# Compute sum squared error

denominator = np.sqrt(np.sum(x_demean ** 2, axis=0) *

np.sum(y_demean ** 2, axis=0))