Skip to content

Commit

Permalink
Merge pull request #16 from tsalo/selcomps
Browse files Browse the repository at this point in the history 
[FIX] Fix errors in tedana workflow and selcomps
  • Loading branch information
@emdupre
emdupre authored Aug 24, 2018
2 parents 7d4ac90 + 7db535f commit ed0f9a5
Show file tree
Hide file tree
Showing 2 changed files with 90 additions and 73 deletions.
110 changes: 56 additions & 54 deletions tedana/selection/select_comps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,20 +76,20 @@ def selcomps(seldict, mmix, manacc, n_echos):
varex = seldict['varex']
Z_maps = seldict['Z_maps']
Z_clmaps = seldict['Z_clmaps']
F_S0_maps = seldict['F_S0_maps']
F_R2_maps = seldict['F_R2_maps']
F_S0_clmaps = seldict['F_S0_clmaps']
F_R2_clmaps = seldict['F_R2_clmaps']
Br_clmaps_S0 = seldict['Br_clmaps_S0']
Br_clmaps_R2 = seldict['Br_clmaps_R2']
Br_S0_clmaps = seldict['Br_clmaps_S0']
Br_R2_clmaps = seldict['Br_clmaps_R2']

n_vols = mmix.shape[1]
n_comps, n_vols = mmix.shape
if len(kappas) != n_comps:
raise ValueError('Number of components does not match between seldict '
'({0}) and mmix ({1}).'.format(len(kappas), n_comps))

# List of components
all_comps = np.arange(len(kappas))
acc = np.arange(len(kappas))
midk = []
ign = []
all_comps = np.arange(n_comps)
acc = np.arange(n_comps)

# If user has specified
if manacc:
Expand All @@ -114,26 +114,26 @@ def selcomps(seldict, mmix, manacc, n_echos):
"""
countsigFS0 = F_S0_clmaps.sum(axis=0)
countsigFR2 = F_R2_clmaps.sum(axis=0)
countnoise = np.zeros(len(all_comps))
countnoise = np.zeros(n_comps)

"""
Make table of dice values
"""
dice_table = np.zeros([all_comps.shape[0], 2])
for i_comp in acc:
dice_FR2 = utils.dice(Br_clmaps_R2[:, i_comp], F_R2_clmaps[:, i_comp])
dice_FS0 = utils.dice(Br_clmaps_S0[:, i_comp], F_S0_clmaps[:, i_comp])
dice_FR2 = utils.dice(Br_R2_clmaps[:, i_comp], F_R2_clmaps[:, i_comp])
dice_FS0 = utils.dice(Br_S0_clmaps[:, i_comp], F_S0_clmaps[:, i_comp])
dice_table[i_comp, :] = [dice_FR2, dice_FS0] # step 3a here and above
dice_table[np.isnan(dice_table)] = 0

"""
Make table of noise gain
"""
tt_table = np.zeros([len(all_comps), 4])
counts_FR2_Z = np.zeros([len(all_comps), 2])
tt_table = np.zeros([n_comps, 4])
counts_FR2_Z = np.zeros([n_comps, 2])
for i_comp in all_comps:
comp_noise_sel = utils.andb([np.abs(Z_maps[:, i_comp]) > 1.95,
Z_clmaps[:, i_comp] == 0]) == 2
comp_noise_sel = ((np.abs(Z_maps[:, i_comp]) > 1.95) &
(Z_clmaps[:, i_comp] == 0))
countnoise[i_comp] = np.array(comp_noise_sel, dtype=np.int).sum()
noise_FR2_Z = np.log10(np.unique(F_R2_maps[comp_noise_sel, i_comp]))
signal_FR2_Z = np.log10(np.unique(F_R2_maps[Z_clmaps[:, i_comp] == 1,
Expand All @@ -147,28 +147,28 @@ def selcomps(seldict, mmix, manacc, n_echos):
Assemble decision table
"""
d_table_rank = np.vstack([
len(all_comps)-stats.rankdata(kappas, method='ordinal'),
len(all_comps)-stats.rankdata(dice_table[:, 0], method='ordinal'),
len(all_comps)-stats.rankdata(tt_table[:, 0], method='ordinal'),
n_comps-stats.rankdata(kappas, method='ordinal'),
n_comps-stats.rankdata(dice_table[:, 0], method='ordinal'),
n_comps-stats.rankdata(tt_table[:, 0], method='ordinal'),
stats.rankdata(countnoise, method='ordinal'),
len(all_comps)-stats.rankdata(countsigFR2, method='ordinal')]).T
n_comps-stats.rankdata(countsigFR2, method='ordinal')]).T
d_table_score = d_table_rank.sum(1)

"""
Step 1: Reject anything that's obviously an artifact
a. Estimate a null variance
"""
rej = acc[utils.andb([rhos > kappas, countsigFS0 > countsigFR2]) > 0]
rej = np.union1d(rej, acc[utils.andb([dice_table[:, 1] > dice_table[:, 0],
varex > np.median(varex)]) == 2])
rej = np.union1d(rej,
acc[utils.andb([tt_table[acc, 0] < 0,
varex[acc] > np.median(varex)]) == 2])
rej = acc[(rhos > kappas) | (countsigFS0 > countsigFR2)]
rej = np.union1d(
acc[(dice_table[:, 1] > dice_table[:, 0]) & (varex > np.median(varex))],
rej)
rej = np.union1d(
acc[(tt_table[acc, 0] < 0) & (varex[acc] > np.median(varex))],
rej)
acc = np.setdiff1d(acc, rej)
varex_ub_p = np.median(varex[kappas > kappas[getelbow_mod(kappas)]])

"""
Step 2: Make a guess for what the good components are, in order to
Step 2: Make a guess for what the good components are, in order to
estimate good component properties
a. Not outlier variance
b. Kappa>kappa_elbow
Expand All @@ -178,19 +178,23 @@ def selcomps(seldict, mmix, manacc, n_echos):
f. Estimate a low and high variance
"""
# Step 2a
varex_ub_p = np.median(varex[kappas > kappas[getelbow_mod(kappas)]])
ncls = acc.copy()
# NOTE: We're not sure why this is done, nor why it's specifically done
# three times. Need to look into this deeper, esp. to make sure the 3
# isn't a hard-coded reference to the number of echoes.
for nn in range(3):
ncls = ncls[1:][(varex[ncls][1:] - varex[ncls][:-1]) < varex_ub_p]

# Compute elbows
kappas_lim = kappas[kappas < utils.getfbounds(n_echos)[-1]]
rhos_lim = np.array(sorted(rhos[ncls])[::-1])
rhos_sorted = np.array(sorted(rhos)[::-1])
kappa_elbow = min(kappas_lim[getelbow_mod(kappas_lim)],
kappas[getelbow_mod(kappas)])
rho_elbow = np.mean([rhos_lim[getelbow_mod(rhos_lim)],
rhos_sorted[getelbow_mod(rhos_sorted)],
kappa_elbow = min(getelbow_mod(kappas_lim, return_val=True),
getelbow_mod(kappas, return_val=True))
rho_elbow = np.mean([getelbow_mod(rhos[ncls], return_val=True),
getelbow_mod(rhos, return_val=True),
utils.getfbounds(n_echos)[0]])
good_guess = ncls[utils.andb([kappas[ncls] >= kappa_elbow,
rhos[ncls] < rho_elbow]) == 2]

good_guess = ncls[(kappas[ncls] >= kappa_elbow) & (rhos[ncls] < rho_elbow)]

if len(good_guess) == 0:
return [], sorted(rej), [], sorted(np.setdiff1d(all_comps, rej))
Expand All @@ -205,23 +209,22 @@ def selcomps(seldict, mmix, manacc, n_echos):
Step 3: Get rid of midk components; i.e., those with higher than
max decision score and high variance
"""
max_good_d_score = EXTEND_FACTOR*len(good_guess) * d_table_rank.shape[1]
midkadd = acc[utils.andb([
d_table_score[acc] > max_good_d_score,
varex[acc] > EXTEND_FACTOR*varex_ub]) == 2]
midk = np.union1d(midkadd, midk)
max_good_d_score = EXTEND_FACTOR * len(good_guess) * d_table_rank.shape[1]
midk_add = acc[(d_table_score[acc] > max_good_d_score) &
(varex[acc] > EXTEND_FACTOR * varex_ub)]
midk = midk_add.astype(int)
acc = np.setdiff1d(acc, midk)

"""
Step 4: Find components to ignore
"""
good_guess = np.setdiff1d(good_guess, midk)
loaded = np.union1d(good_guess, acc[varex[acc] > varex_lb])
igncand = np.setdiff1d(acc, loaded)
igncand = np.setdiff1d(igncand,
igncand[d_table_score[igncand] < max_good_d_score])
igncand = np.setdiff1d(igncand, igncand[kappas[igncand] > kappa_elbow])
ign = np.array(np.union1d(ign, igncand), dtype=np.int)
ign_cand = np.setdiff1d(acc, loaded)
ign_cand = np.setdiff1d(
ign_cand, ign_cand[d_table_score[ign_cand] < max_good_d_score])
ign_cand = np.setdiff1d(ign_cand, ign_cand[kappas[ign_cand] > kappa_elbow])
ign = ign_cand.astype(int)
acc = np.setdiff1d(acc, ign)

"""
Expand All @@ -237,29 +240,28 @@ def selcomps(seldict, mmix, manacc, n_echos):
len(acc) - stats.rankdata(countsigFR2[acc], method='ordinal')]).T
d_table_score = d_table_rank.sum(1)
num_acc_guess = np.mean([
np.sum(utils.andb([kappas[acc] > kappa_elbow,
rhos[acc] < rho_elbow]) == 2),
np.sum((kappas[acc] > kappa_elbow) & (rhos[acc] < rho_elbow)),
np.sum(kappas[acc] > kappa_elbow)])
conservative_guess = num_acc_guess * d_table_rank.shape[1] / RESTRICT_FACTOR
candartA = np.intersect1d(acc[d_table_score > conservative_guess],
acc[kappa_ratios[acc] > EXTEND_FACTOR * 2])
midkadd = np.union1d(midkadd, np.intersect1d(
midk_add = np.union1d(midk_add, np.intersect1d(
candartA, candartA[varex[candartA] > varex_ub * EXTEND_FACTOR]))
candartB = acc[d_table_score > num_acc_guess * d_table_rank.shape[1] * HIGH_PERC / 100.]
midkadd = np.union1d(midkadd, np.intersect1d(
midk_add = np.union1d(midk_add, np.intersect1d(
candartB, candartB[varex[candartB] > varex_lb * EXTEND_FACTOR]))
midk = np.union1d(midk, midkadd)
midk = np.union1d(midk, midk_add)

# Find comps to ignore
new_varex_lb = stats.scoreatpercentile(varex[acc[:num_acc_guess]],
LOW_PERC)
candart = np.setdiff1d(
acc[d_table_score > num_acc_guess * d_table_rank.shape[1]], midk)
ignadd = np.intersect1d(candart[varex[candart] > new_varex_lb],
candart)
ignadd = np.union1d(ignadd, np.intersect1d(
ign_add = np.intersect1d(candart[varex[candart] > new_varex_lb],
candart)
ign_add = np.union1d(ign_add, np.intersect1d(
acc[kappas[acc] <= kappa_elbow], acc[varex[acc] > new_varex_lb]))
ign = np.setdiff1d(np.union1d(ign, ignadd), midk)
ign = np.setdiff1d(np.union1d(ign, ign_add), midk)
acc = np.setdiff1d(acc, np.union1d(midk, ign))

acc = list(sorted(acc))
Expand Down
53 changes: 34 additions & 19 deletions tedana/workflows/tedana.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -312,6 +312,10 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None,
dn_ts_e[echo].nii Denoised timeseries for echo number ``echo``
====================== =================================================
"""
METHOD = 'kundu_v2_5'
if METHOD not in ['kundu_v2_5', 'kundu_v3_2']:
raise ValueError('Component selection must be either "kundu_v2_5" or '
'"kundu_v3_2"')

# ensure tes are in appropriate format
tes = [float(te) for te in tes]
Expand Down Expand Up @@ -395,34 +399,45 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None,
t2s, t2sG, stabilize, ref_img,
tes=tes, kdaw=kdaw, rdaw=rdaw,
ste=ste, wvpca=wvpca)
mmix_orig, fixed_seed = decomposition.tedica(n_components, dd, conv, fixed_seed,
cost=initcost, final_cost=finalcost,
mmix_orig, fixed_seed = decomposition.tedica(n_components, dd, conv,
fixed_seed, cost=initcost,
final_cost=finalcost,
verbose=debug)
np.savetxt(op.join(out_dir, '__meica_mix.1D'), mmix_orig)
LGR.info('Making second component selection guess from ICA results')
seldict, comptable, betas, mmix = model.fitmodels_direct(catd, mmix_orig,
mask, t2s, t2sG,
tes, combmode,
ref_img,
reindex=True)
(seldict, comptable,
betas, mmix) = model.fitmodels_direct(catd, mmix_orig, mask, t2s,
t2sG, tes, combmode, ref_img,
reindex=True)
np.savetxt(op.join(out_dir, 'meica_mix.1D'), mmix)

acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask, ref_img, manacc,
n_echos, t2s, s0, strict_mode=strict,
filecsdata=filecsdata)
if METHOD == 'kundu_v2_5':
acc, rej, midk, empty = selection.selcomps(seldict, mmix, manacc,
n_echos)
elif METHOD == 'kundu_v3_2':
acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask,
ref_img, manacc,
n_echos, t2s, s0,
strict_mode=strict,
filecsdata=filecsdata)
else:
LGR.info('Using supplied mixing matrix from ICA')
mmix_orig = np.loadtxt(op.join(out_dir, 'meica_mix.1D'))
seldict, comptable, betas, mmix = model.fitmodels_direct(catd, mmix_orig,
mask, t2s, t2sG,
tes, combmode,
ref_img)
(seldict, comptable,
betas, mmix) = model.fitmodels_direct(catd, mmix_orig, mask, t2s,
t2sG, tes, combmode, ref_img)
if ctab is None:
acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask,
ref_img, manacc,
n_echos, t2s, s0,
filecsdata=filecsdata,
strict_mode=strict)
if METHOD == 'kundu_v2_5':
(acc, rej,
midk, empty) = selection.selcomps(seldict, mmix, manacc,
n_echos)
elif METHOD == 'kundu_v3_2':
(acc, rej,
midk, empty) = selection.selcomps(seldict, mmix, mask,
ref_img, manacc,
n_echos, t2s, s0,
filecsdata=filecsdata,
strict_mode=strict)
else:
acc, rej, midk, empty = utils.ctabsel(ctab)

Expand Down

0 comments on commit ed0f9a5

Please sign in to comment.