import nibabel as ni
import numpy as np
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
import matplotlib.pyplot as plt
%matplotlib inline
We set the data-path to point to the installation-specific location of the data-files
import os
import osmosis as oz
import osmosis.io as oio
oio.data_path = os.path.join(oz.__path__[0], 'data')
We generate strings pointing to all of the data-files for this subject:
subject = 'SUB1'
data_1k_1, data_1k_2 = oio.get_dwi_data(1000, subject)
data_2k_1, data_2k_2 = oio.get_dwi_data(2000, subject)
data_4k_1, data_4k_2 = oio.get_dwi_data(4000, subject)
Wewill only look at white-matter voxels. This loads the white-matter mask:
wm_mask = np.zeros(ni.load(data_1k_1[0]).shape[:3])
wm_nifti = ni.load(oio.data_path + '/%s/%s_wm_mask.nii.gz'%(subject, subject)).get_data()
wm_mask[np.where(wm_nifti==1)] = 1
# Start working on a small chunk:
#wm_mask = np.zeros(wm_mask.shape)
#wm_mask[40:42, 40:42, 40:42] = 1
The following are the parameters defining the performance of the Elastic Net solver , based on the calculations in AppendixA
Initialize the class instances of the SFM objects, using AD/RD derived from the corpus callosum. Parameters are not saved (hence temp
ad_rd = oio.get_ad_rd(subject, 2000)
SD_2k_1_nnls_demean = ssd.SparseDeconvolutionModel(*data_2k_1, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'],
solver='nnls')
SD_2k_2_nnls_demean = ssd.SparseDeconvolutionModel(*data_2k_2, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'],
solver='nnls')
SD_2k_1_nnls = ssd.SparseDeconvolutionModel(*data_2k_1, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'],
solver='nnls', demean=False)
SD_2k_2_nnls = ssd.SparseDeconvolutionModel(*data_2k_2, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'],
solver='nnls', demean=False)
SD_2k_1_EN_demean = ssd.SparseDeconvolutionModel(*data_2k_1, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'])
SD_2k_2_EN_demean = ssd.SparseDeconvolutionModel(*data_2k_2, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'])
SD_2k_1_EN = ssd.SparseDeconvolutionModel(*data_2k_1, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'],
demean=False)
SD_2k_2_EN = ssd.SparseDeconvolutionModel(*data_2k_2, mask=wm_mask, params_file ='temp',
axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'],
demean=False)
Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.nii.gz Loading from file: /home/arokem/usr/local/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs
Calculate rRMSEfor each b value , from a cross-prediction: fit the model to one data-set in each b value and predict the other data set. Calculations average across both directions (1=>2)
rrmse_2k_nnls = oza.cross_predict(SD_2k_1_nnls, SD_2k_2_nnls)
rrmse_2k_nnls_demean = oza.cross_predict(SD_2k_1_nnls_demean, SD_2k_2_nnls_demean)
rrmse_2k_EN = oza.cross_predict(SD_2k_1_EN, SD_2k_2_EN)
rrmse_2k_EN_demean = oza.cross_predict(SD_2k_1_EN, SD_2k_2_EN_demean)
SparseDeconvolutionModel.model_params [******** 22% ] 14910 of 68694 complete
#print(rrmse_2k_nnls[np.where(wm_mask==1)])
#print(rrmse_2k_nnls_demean[np.where(wm_mask==1)])
#print(rrmse_2k_EN_demean[np.where(wm_mask==1)])
#print(rrmse_2k_EN[np.where(wm_mask==1)])
rrmse_2k_nnls_demean[np.isfinite(rrmse_2k_nnls_demean)]
array([ 0.91143276, 0.82941019, 0.89216963, ..., 0.93824783, 0.85796653, 0.96076982])
fig = mpl.probability_hist(rrmse_2k_nnls[np.isfinite(rrmse_2k_nnls)], label='NNLS')
fig = mpl.probability_hist(rrmse_2k_nnls_demean[np.isfinite(rrmse_2k_nnls_demean)], fig=fig, label='NNLS (demean)')
fig.axes[0].plot([1,1], fig.axes[0].get_ylim(), '--k')
fig.axes[0].plot([1/np.sqrt(2),1/np.sqrt(2)], fig.axes[0].get_ylim(), '--k')
fig.axes[0].set_xlim([0.6,1.4])
plt.legend()
fig = mpl.probability_hist(rrmse_2k_EN[np.isfinite(rrmse_2k_EN)], label='Elastic Net')
fig = mpl.probability_hist(rrmse_2k_EN_demean[np.isfinite(rrmse_2k_EN_demean)], fig=fig, label='Elastic Net (demean)')
fig.axes[0].plot([1,1], fig.axes[0].get_ylim(), '--k')
fig.axes[0].plot([1/np.sqrt(2),1/np.sqrt(2)], fig.axes[0].get_ylim(), '--k')
fig.axes[0].set_xlim([0.6,1.4])
plt.legend()
<matplotlib.legend.Legend at 0x108dbd90>
whos
a = 1
fig = mpl.probability_hist(rrmse_2k_nnls[np.isfinite(rrmse_2k_nnls)], label='NNLS')
fig = mpl.probability_hist(rrmse_2k_EN[np.isfinite(rrmse_2k_EN)], fig=fig, label='Elastic Net ')
fig.axes[0].plot([1,1], fig.axes[0].get_ylim(), '--k')
fig.axes[0].plot([1/np.sqrt(2),1/np.sqrt(2)], fig.axes[0].get_ylim(), '--k')
fig.axes[0].set_xlim([0.6,1.4])
plt.legend()
<matplotlib.legend.Legend at 0x10556850>