{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#import os\n",
      "#os.environ['DISPLAY'] = 'localhost:99'\n",
      "\n",
      "#from IPython.display import Image, display"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import tempfile\n",
      "from IPython.display import Image, display\n",
      "\n",
      "import nibabel as ni\n",
      "\n",
      "import osmosis.model.analysis as ozm\n",
      "import osmosis.model.dti as dti\n",
      "import osmosis.model.sparse_deconvolution as ssd\n",
      "import osmosis.viz.maya as maya\n",
      "import osmosis.viz.mpl as mpl\n",
      "\n",
      "import osmosis.utils as ozu\n",
      "import osmosis.volume as ozv\n",
      "import osmosis.tensor as ozt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import osmosis as oz\n",
      "import osmosis.io as oio\n",
      "oio.data_path = os.path.join(oz.__path__[0], 'data')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subject = 'SUB1'\n",
      "data_1k_1, data_1k_2 = oio.get_dwi_data(1000, subject)\n",
      "data_2k_1, data_2k_2 = oio.get_dwi_data(2000, subject)\n",
      "data_4k_1, data_4k_2 = oio.get_dwi_data(4000, subject)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vox_idx = (53, 43, 47)\n",
      "mask = np.zeros(ni.load(data_1k_1[0]).shape[:3])\n",
      "mask[vox_idx[0]:vox_idx[0]+2, vox_idx[1], vox_idx[2]]=1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "TM_1k_1 = dti.TensorModel(*data_1k_1, mask=mask, params_file='temp')\n",
      "TM_1k_2 = dti.TensorModel(*data_1k_2, mask=mask, params_file='temp')\n",
      "TM_2k_1 = dti.TensorModel(*data_2k_1, mask=mask, params_file='temp')\n",
      "TM_2k_2 = dti.TensorModel(*data_2k_2, mask=mask, params_file='temp')\n",
      "TM_4k_1 = dti.TensorModel(*data_4k_1, mask=mask, params_file='temp')\n",
      "TM_4k_2 = dti.TensorModel(*data_4k_2, mask=mask, params_file='temp')\n",
      "\n",
      "alpha = 0.0005\n",
      "l1_ratio = 0.6\n",
      "\n",
      "solver_params = dict(alpha=alpha,\n",
      "                     l1_ratio=l1_ratio,\n",
      "                     fit_intercept=False,\n",
      "                     positive=True)\n",
      "\n",
      "\n",
      "ad_rd = oio.get_ad_rd(subject, 1000)\n",
      "SD_1k_1 = ssd.SparseDeconvolutionModel(*data_1k_1, mask=mask, solver_params=solver_params, params_file='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'])\n",
      "SD_1k_2 = ssd.SparseDeconvolutionModel(*data_1k_2, mask=mask, solver_params=solver_params, params_file='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'])\n",
      "ad_rd = oio.get_ad_rd(subject, 2000)\n",
      "SD_2k_1 = ssd.SparseDeconvolutionModel(*data_2k_1, mask=mask, solver_params=solver_params, params_file='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'])\n",
      "SD_2k_2 = ssd.SparseDeconvolutionModel(*data_2k_2, mask=mask, solver_params=solver_params, params_file='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'])\n",
      "ad_rd = oio.get_ad_rd(subject, 4000)\n",
      "SD_4k_1 = ssd.SparseDeconvolutionModel(*data_4k_1, mask=mask, solver_params=solver_params, params_file='temp', axial_diffusivity=ad_rd[0]['AD'], radial_diffusivity=ad_rd[0]['RD'])\n",
      "SD_4k_2 = ssd.SparseDeconvolutionModel(*data_4k_2, mask=mask, solver_params=solver_params, params_file='temp', axial_diffusivity=ad_rd[1]['AD'], radial_diffusivity=ad_rd[1]['RD'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0009_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0011_01_DWI_2mm150dir_2x_b1000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DTI_2mm_150dir_2x_b2000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0005_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvals\n",
        "Loading from file: /home/arokem/usr/lib/python2.7/site-packages/osmosis/data/FP/0007_01_DWI_2mm150dir_2x_b4000_aligned_trilin.bvecs\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ten_rrmse_1k = ozm.cross_predict(TM_1k_1, TM_1k_2)\n",
      "sfm_rrmse_1k = ozm.cross_predict(SD_1k_1, SD_1k_2)\n",
      "\n",
      "ten_rrmse_2k = ozm.cross_predict(TM_2k_1, TM_2k_2)\n",
      "sfm_rrmse_2k = ozm.cross_predict(SD_2k_1, SD_2k_2)\n",
      "\n",
      "ten_rrmse_4k = ozm.cross_predict(TM_4k_1, TM_4k_2)\n",
      "sfm_rrmse_4k = ozm.cross_predict(SD_4k_1, SD_4k_2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r",
        "SparseDeconvolutionModel.model_params [*****************50%                  ]  1 of 2 complete "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%gui wx\n",
      "\n",
      "fig = maya.plot_signal_interp(TM_4k_1.bvecs[:,TM_4k_1.b_idx], TM_4k_1.relative_signal[vox_idx], cmap='hot', colorbar=True, vmin=0, vmax=1)#, figure=fig, offset=-4, roll=90, file_name=fn)\n",
      "fig = maya.plot_signal_interp(TM_4k_1.bvecs[:,TM_4k_2.b_idx], TM_4k_2.relative_signal[vox_idx], cmap='hot', colorbar=True, figure=fig, vmin=0, vmax=1, origin = [0,0,-1])#, figure=fig, offset=-4, roll=90, file_name=fn)\n",
      "fig = maya.plot_signal_interp(TM_4k_1.bvecs[:,TM_4k_1.b_idx], TM_4k_1.fit[vox_idx]/TM_4k_1.S0[vox_idx], cmap='hot', colorbar=True, figure=fig, origin = [0,0,1], roll=90, vmin=0, vmax=1)#, file_name=fn)\n",
      "fig = maya.plot_signal_interp(SD_4k_1.bvecs[:,SD_4k_1.b_idx], SD_4k_1.fit[vox_idx]/SD_4k_1.S0[vox_idx], cmap='hot', colorbar=True, figure=fig, origin = [0,0,-2], roll=90, vmin=0, vmax=1)#, file_name=fn)\n",
      "\n",
      "print(TM_4k_1.fractional_anisotropy[vox_idx])\n",
      "print(ten_rrmse_4k[vox_idx])\n",
      "print(sfm_rrmse_4k[vox_idx])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Predicting signal from TensorModel\n",
        "Predicting signal from SparseDeconvolutionModel with ElasticNet(alpha=0.0005, copy_X=True, fit_intercept=False, l1_ratio=0.6,\n",
        "      max_iter=1000, normalize=False, positive=True, precompute=auto,\n",
        "      rho=None, tol=0.0001, warm_start=False)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "0.247992798372"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1.21245936767\n",
        "0.811665956894\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/python2.7/dist-packages/scipy/interpolate/rbf.py:122: RuntimeWarning: divide by zero encountered in log\n",
        "  result = r**2 * log(r)\n",
        "/usr/lib/python2.7/dist-packages/scipy/interpolate/rbf.py:122: RuntimeWarning: invalid value encountered in multiply\n",
        "  result = r**2 * log(r)\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%gui wx\n",
      "fig = maya.plot_signal_interp(TM_2k_1.bvecs[:,TM_2k_1.b_idx], TM_2k_1.relative_signal[vox_idx], cmap='hot', colorbar=True, vmin=0, vmax=1)#, figure=fig, offset=-4, roll=90, file_name=fn)\n",
      "fig = maya.plot_signal_interp(TM_2k_1.bvecs[:,TM_2k_2.b_idx], TM_2k_2.relative_signal[vox_idx], cmap='hot', colorbar=True, figure=fig, vmin=0, vmax=1, origin = [0,0,-1.1])#, figure=fig, offset=-4, roll=90, file_name=fn)\n",
      "fig = maya.plot_signal_interp(TM_2k_1.bvecs[:,TM_2k_1.b_idx], TM_2k_1.fit[vox_idx]/TM_2k_1.S0[vox_idx], cmap='hot', colorbar=True, figure=fig, origin = [0,0,1.1], roll=90, vmin=0, vmax=1)#, file_name=fn)\n",
      "fig = maya.plot_signal_interp(SD_2k_1.bvecs[:,SD_2k_1.b_idx], SD_2k_1.fit[vox_idx]/SD_2k_1.S0[vox_idx], cmap='hot', colorbar=True, figure=fig, origin = [0,0,-2.1], roll=90, vmin=0, vmax=1)#, file_name=fn)\n",
      "\n",
      "#print(TM_2k_1.fractional_anisotropy[vox_idx])\n",
      "#print(ten_rrmse_2k[vox_idx])\n",
      "#print(sfm_rrmse_2k[vox_idx])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Predicting signal from TensorModel\n",
        "Predicting signal from SparseDeconvolutionModel with ElasticNet(alpha=0.0005, copy_X=True, fit_intercept=False, l1_ratio=0.6,\n",
        "      max_iter=1000, normalize=False, positive=True, precompute=auto,\n",
        "      rho=None, tol=0.0001, warm_start=False)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%gui wx\n",
      "fig = maya.plot_signal_interp(TM_1k_1.bvecs[:,TM_1k_1.b_idx], TM_1k_1.relative_signal[vox_idx], cmap='hot', colorbar=True, vmin=0, vmax=1)#, figure=fig, offset=-4, roll=90, file_name=fn)\n",
      "fig = maya.plot_signal_interp(TM_1k_1.bvecs[:,TM_1k_1.b_idx], TM_1k_2.relative_signal[vox_idx], cmap='hot', colorbar=True, figure=fig, vmin=0, vmax=1, origin = [0,0,-1.1])#, figure=fig, offset=-4, roll=90, file_name=fn)\n",
      "fig = maya.plot_signal_interp(TM_1k_1.bvecs[:,TM_1k_1.b_idx], TM_1k_1.fit[vox_idx]/TM_1k_1.S0[vox_idx], cmap='hot', colorbar=True, figure=fig, origin = [0,0,1.1], roll=90, vmin=0, vmax=1)#, file_name=fn)\n",
      "fig = maya.plot_signal_interp(SD_1k_1.bvecs[:,SD_1k_1.b_idx], SD_1k_1.fit[vox_idx]/SD_1k_1.S0[vox_idx], cmap='hot', colorbar=True, figure=fig, origin = [0,0,-2.1], roll=90, vmin=0, vmax=1)#, file_name=fn)\n",
      "\n",
      "#print(TM_1k_1.fractional_anisotropy[vox_idx])\n",
      "#print(ten_rrmse_1k[vox_idx])\n",
      "#print(sfm_rrmse_1k[vox_idx])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Predicting signal from TensorModel\n",
        "Predicting signal from SparseDeconvolutionModel with ElasticNet(alpha=0.0005, copy_X=True, fit_intercept=False, l1_ratio=0.6,\n",
        "      max_iter=1000, normalize=False, positive=True, precompute=auto,\n",
        "      rho=None, tol=0.0001, warm_start=False)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    }
   ],
   "metadata": {}
  }
 ]
}