Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
io.filewrite(masksum, op.join(out_dir, 'adaptive_mask.nii'), ref_img)
os.chdir(out_dir)
LGR.info('Computing T2* map')
t2s_limited, s0_limited, t2s_full, s0_full = decay.fit_decay(
catd, tes, mask, masksum, fittype)
# set a hard cap for the T2* map
# anything that is 10x higher than the 99.5 %ile will be reset to 99.5 %ile
cap_t2s = stats.scoreatpercentile(t2s_limited.flatten(), 99.5,
interpolation_method='lower')
LGR.debug('Setting cap on T2* map at {:.5f}'.format(cap_t2s * 10))
t2s_limited[t2s_limited > cap_t2s * 10] = cap_t2s
io.filewrite(t2s_limited, op.join(out_dir, 't2sv.nii'), ref_img)
io.filewrite(s0_limited, op.join(out_dir, 's0v.nii'), ref_img)
if verbose:
io.filewrite(t2s_full, op.join(out_dir, 't2svG.nii'), ref_img)
io.filewrite(s0_full, op.join(out_dir, 's0vG.nii'), ref_img)
# optimally combine data
data_oc = combine.make_optcom(catd, tes, mask, t2s=t2s_full, combmode=combmode)
# regress out global signal unless explicitly not desired
if 'gsr' in gscontrol:
catd, data_oc = gsc.gscontrol_raw(catd, data_oc, n_echos, ref_img)
if mixm is None:
# Identify and remove thermal noise from data
dd, n_components = decomposition.tedpca(catd, data_oc, combmode, mask,
t2s_limited, t2s_full, ref_img,
LGR.warning('Could not save PCA solution')
else: # if loading existing state
voxel_comp_weights = pcastate['voxel_comp_weights']
varex = pcastate['varex']
comp_ts = pcastate['comp_ts']
ct_df = pcastate['comptable']
np.savetxt('mepca_mix.1D', comp_ts.T)
# write component maps to 4D image
comp_maps = np.zeros((OCcatd.shape[0], comp_ts.shape[0]))
for i_comp in range(comp_ts.shape[0]):
temp_comp_ts = comp_ts[i_comp, :][:, None]
comp_map = utils.unmask(model.computefeats2(OCcatd, temp_comp_ts, mask), mask)
comp_maps[:, i_comp] = np.squeeze(comp_map)
io.filewrite(comp_maps, 'mepca_OC_components.nii', ref_img)
# Add new columns to comptable for classification
ct_df['classification'] = 'accepted'
ct_df['rationale'] = ''
# Select components using decision tree
if method == 'kundu':
ct_df = kundu_tedpca(ct_df, n_echos, kdaw, rdaw, stabilize=False)
elif method == 'kundu-stabilize':
ct_df = kundu_tedpca(ct_df, n_echos, kdaw, rdaw, stabilize=True)
elif method == 'mle':
LGR.info('Selected {0} components with MLE dimensionality '
'detection'.format(ct_df.shape[0]))
ct_df.to_csv('comp_table_pca.txt', sep='\t', index=True,
index_label='component', float_format='%.6f')
# Legendre polynomial basis for denoising
bounds = np.linspace(-1, 1, optcom.shape[-1])
Lmix = np.column_stack([lpmv(0, vv, bounds) for vv in range(dtrank)])
# compute mean, std, mask local to this function
# inefficient, but makes this function a bit more modular
Gmu = optcom.mean(axis=-1) # temporal mean
Gmask = Gmu != 0
# find spatial global signal
dat = optcom[Gmask] - Gmu[Gmask][:, np.newaxis]
sol = np.linalg.lstsq(Lmix, dat.T, rcond=None)[0] # Legendre basis for detrending
detr = dat - np.dot(sol.T, Lmix.T)[0]
sphis = (detr).min(axis=1)
sphis -= sphis.mean()
io.filewrite(utils.unmask(sphis, Gmask), 'T1gs', ref_img)
# find time course ofc the spatial global signal
# make basis with the Legendre basis
glsig = np.linalg.lstsq(np.atleast_2d(sphis).T, dat, rcond=None)[0]
glsig = stats.zscore(glsig, axis=None)
np.savetxt('glsig.1D', glsig)
glbase = np.hstack([Lmix, glsig.T])
# Project global signal out of optimally combined data
sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T, rcond=None)[0]
tsoc_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T,
np.atleast_2d(glbase.T[dtrank])) + Gmu[Gmask][:, np.newaxis]
io.filewrite(optcom, 'tsoc_orig', ref_img)
dm_optcom = utils.unmask(tsoc_nogs, Gmask)
io.filewrite(dm_optcom, 'tsoc_nogs', ref_img)
comptable = metrics.kundu_metrics(comptable, metric_maps)
comptable = selection.kundu_selection_v2(comptable, n_echos, n_vols)
else:
comptable = pd.read_csv(ctab, sep='\t', index_col='component')
comptable = selection.manual_selection(comptable, acc=manacc)
# Save decomposition
data_type = 'optimally combined data' if source_tes == -1 else 'z-concatenated data'
comptable['Description'] = 'ICA fit to dimensionally reduced {0}.'.format(data_type)
mmix_dict = {}
mmix_dict['Method'] = ('Independent components analysis with FastICA '
'algorithm implemented by sklearn. Components '
'are sorted by Kappa in descending order. '
'Component signs are flipped to best match the '
'data.')
io.save_comptable(comptable, op.join(out_dir, 'ica_decomposition.json'),
label='ica', metadata=mmix_dict)
if comptable[comptable.classification == 'accepted'].shape[0] == 0:
LGR.warning('No BOLD components detected! Please check data and '
'results!')
mmix_orig = mmix.copy()
if tedort:
acc_idx = comptable.loc[
~comptable.classification.str.contains('rejected')].index.values
rej_idx = comptable.loc[
comptable.classification.str.contains('rejected')].index.values
acc_ts = mmix[:, acc_idx]
rej_ts = mmix[:, rej_idx]
betas = np.linalg.lstsq(acc_ts, rej_ts, rcond=None)[0]
pred_rej_ts = np.dot(acc_ts, betas)
comptable.classification.str.contains('rejected')].index.values
acc_ts = mmix[:, acc_idx]
rej_ts = mmix[:, rej_idx]
betas = np.linalg.lstsq(acc_ts, rej_ts, rcond=None)[0]
pred_rej_ts = np.dot(acc_ts, betas)
resid = rej_ts - pred_rej_ts
mmix[:, rej_idx] = resid
comp_names = [io.add_decomp_prefix(comp, prefix='ica', max_value=comptable.index.max())
for comp in comptable.index.values]
mixing_df = pd.DataFrame(data=mmix, columns=comp_names)
mixing_df.to_csv('ica_orth_mixing.tsv', sep='\t', index=False)
RepLGR.info("Rejected components' time series were then "
"orthogonalized with respect to accepted components' time "
"series.")
io.writeresults(data_oc, mask=mask, comptable=comptable, mmix=mmix,
n_vols=n_vols, ref_img=ref_img)
if 't1c' in gscontrol:
gsc.gscontrol_mmix(data_oc, mmix, mask, comptable, ref_img)
if verbose:
io.writeresults_echoes(catd, mmix, mask, comptable, ref_img)
if not no_png:
LGR.info('Making figures folder with static component maps and '
'timecourse plots.')
# make figure folder first
if not op.isdir(op.join(out_dir, 'figures')):
os.mkdir(op.join(out_dir, 'figures'))
viz.write_comp_figs(data_oc, mask=mask, comptable=comptable,
comp_names = [io.add_decomp_prefix(comp, prefix='ica', max_value=comptable.index.max())
for comp in comptable.index.values]
mixing_df = pd.DataFrame(data=mmix, columns=comp_names)
mixing_df.to_csv('ica_orth_mixing.tsv', sep='\t', index=False)
RepLGR.info("Rejected components' time series were then "
"orthogonalized with respect to accepted components' time "
"series.")
io.writeresults(data_oc, mask=mask, comptable=comptable, mmix=mmix,
n_vols=n_vols, ref_img=ref_img)
if 't1c' in gscontrol:
gsc.gscontrol_mmix(data_oc, mmix, mask, comptable, ref_img)
if verbose:
io.writeresults_echoes(catd, mmix, mask, comptable, ref_img)
if not no_png:
LGR.info('Making figures folder with static component maps and '
'timecourse plots.')
# make figure folder first
if not op.isdir(op.join(out_dir, 'figures')):
os.mkdir(op.join(out_dir, 'figures'))
viz.write_comp_figs(data_oc, mask=mask, comptable=comptable,
mmix=mmix_orig, ref_img=ref_img,
out_dir=op.join(out_dir, 'figures'),
png_cmap=png_cmap)
LGR.info('Making Kappa vs Rho scatter plot')
viz.write_kappa_scatter(comptable=comptable,
out_dir=op.join(out_dir, 'figures'))