This takes a (very long) while to run unless the tools in osmosis.parallel are used to parallelize the computation of all the model params and the rRMSE values.
For completeness, the notebook AppendixA will do that for you, but god only knows how long that might take. This here, on the other hand, works FAST!
We import the needed modules:
import nibabel as ni
import os
import osmosis.model.analysis as oza
import osmosis.model.sparse_deconvolution as ssd
import osmosis.model.dti as dti
import osmosis.viz.mpl as mpl
The following points the data-path variable to where the data is stored, relative to the install location
import os
import osmosis as oz
import osmosis.io as oio
oio.data_path = os.path.join(oz.__path__[0], 'data')
Set ht subject to be analyzed:
subject = 'SUB1'
We sweep through particular regularization parameter combinations in each b value:
alphas = [0.0001,0.0005, 0.001, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05]
l1_ratios = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
In the following loop, we collect the median value of the rRMSE from pre-computed rRMSE files for each combination of regularization parameters:
rmse_matrix = {}
rel_matrix = {}
for b in [1000, 2000, 4000]: #
ad_rd = oio.get_ad_rd(subject, b)
rmse_matrix[b] = np.zeros((len(l1_ratios), len(alphas)))
rel_matrix[b] = np.zeros((len(l1_ratios), len(alphas)))
for l1_idx, l1_ratio in enumerate(l1_ratios):
for alpha_idx, this_alpha in enumerate(alphas):
rrmse_fname = "%s/%s/RMSE_b%s_l1ratio%s_alpha%s.nii.gz"%(oio.data_path,
subject,
b,
l1_ratio,
this_alpha)
rrmse_fname = "%s/%s/PDD_angle_b%s_l1ratio%s_alpha%s.nii.gz"%(oio.data_path,
subject,
b,
l1_ratio,
this_alpha)
rrmse = ni.load(oio.data_path + '/%s/'%subject + 'RMSE_b%s_l1ratio%s_alpha%s.nii.gz'%(b, l1_ratio, this_alpha)).get_data()
median_rmse = np.median(rrmse[np.isfinite(rrmse)])
rmse_matrix[b][l1_idx, alpha_idx] = median_rmse
rel = ni.load(oio.data_path + '/%s/'%subject + 'PDD_angle_b%s_l1ratio%s_alpha%s.nii.gz'%(b, l1_ratio, this_alpha)).get_data()
median_rel = np.median(rel[np.isfinite(rel)])
rel_matrix[b][l1_idx, alpha_idx] = median_rel
The results are visualized:
vmax = 1.5
vmin = 0.5
for b in [1000, 2000, 4000]:
fig, ax = plt.subplots(1)
cax = ax.imshow(rmse_matrix[b], interpolation='nearest', cmap=cm.RdYlGn_r, vmax=vmax, vmin=vmin)
cbar = fig.colorbar(cax, ticks=[0.5, 1., 1.5])
cax.axes.set_xticks([0,1,2,3,4,5,6,7,8])
cax.axes.set_xticklabels([str('%1.4f'%this) for this in alphas])
cax.axes.set_xlabel(r'$\lambda$')
cax.axes.set_yticks([0,1,2,3,4,5])
cax.axes.set_yticklabels([str('%1.1f'%(1-this)) for this in l1_ratios])
cax.axes.set_ylabel(r'$\alpha$')
im = ax.matshow(rmse_matrix[b], cmap=matplotlib.cm.RdYlGn_r, vmax=vmax, vmin=vmin)
fig.set_size_inches([8,6])
fig.savefig('figures/AppendixA_acc_b%s.svg'%b)
We also visualize the reliability of the PDD predicted from the SFM model parameters in each combination of regularization parameters:
for b in [1000, 2000, 4000]:
fig, ax = plt.subplots(1)
cax = ax.imshow(rel_matrix[b], interpolation='nearest', cmap=cm.OrRd, vmax=90, vmin=0)
cbar = fig.colorbar(cax)
cax.axes.set_xticks([0,1,2,3,4,5,6,7,8])
cax.axes.set_xticklabels([str('%1.4f'%this) for this in alphas])
cax.axes.set_xlabel(r'$\lambda$')
cax.axes.set_yticks([0,1,2,3,4,5])
cax.axes.set_yticklabels([str('%1.1f'%(1-this)) for this in l1_ratios])
cax.axes.set_ylabel(r'$\alpha$')
#ax.matshow(rmse_matrix[b], cmap=matplotlib.cm.hot)
fig.set_size_inches([8,6])
fig.savefig('figures/AppendixA_rel_b%s.svg'%b)