arokem/t_mrdTensor.ipynb

## t_mrdTensor.ipynb
{
 "metadata": {
  "name": "t_mrdTensor"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "# Calculating the tensor from diffusion weighted images\n\nThis tutorial script demonstrates the calculation of the diffusion tensor (Basser et al. 2004) from diffusion-weighted data.\n\nWe start by setting up the matlab session. We import the pymatbridge, fire up a matlab session and connect to it:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import os\nimport pymatbridge as pymat\nreload(pymat)\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 1,
       "text": "<module 'pymatbridge' from '/usr/local/lib/python2.7/dist-packages/pymatbridge/__init__.pyc'>"
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "ip = get_ipython()\npymat.load_ipython_extension(ip,matlab='/opt/local/matlab/r2010b/bin/matlab')",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Starting MATLAB on http://localhost:54222\n"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "."
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "."
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "."
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "MATLAB started and connected!\n"
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\naddpath(genpath('/home/arokem/source/vistasoft/trunk'))\naddpath(genpath('/home/arokem/source/vistadata'))\n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "We calculate the tensor (Q) from diffusion weighted images.  We call the tensor Q because the tensor is a positive-definite Quadratic form. These are the key actors and relations.\n\nThe 3x3 positive-definite quadratic form, Q, predicts the apparent diffusion data using the formula: \n\n   $ADC = u^t Q u $                                                               \n                                                                              \nwhere $u$ is a unit vector in 3-space that defines the direction and $u^t$ is the transpose of this vector. \n                                                                              \nAccording to the Stejskal-Tanner equation, the diffusion signal is related to the ADC as:\n                                                                              \n   $dSig = S_0 exp^(-b ADC)$\n\nWhere $S_0$ is the signal in the voxel in non diffusion-weighted images and $b$ is a parameter of the measurement.  \n\nThe diffusion ellipsoid (mean diffusion distance) is related to Q by finding the vectors that solve:\n\n$1 = v^t Q^{-1} v$\n\n \nIn the follwing calculations, we will implement these equations step by step and apply them to diffusion-weighted MRI data. \n\nOK - let's start by loading some data from file. The vistadata diffusion sample data are 40-directions.  The directory\ncontains the dwi data as well as the bvals and bvecs."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\ndataDir = fullfile(mrvDataRootPath,'diffusion','sampleData');\ndwi = dwiLoad(fullfile(dataDir,'raw','dwi.nii.gz'));\ndwi",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "text": "\nFiles loaded: \n  dwi   = /home/arokem/source/vistadata/diffusion/sampleData/raw/dwi.nii.gz \n  bvecs = /home/arokem/source/vistadata/diffusion/sampleData/raw/dwi.bvecs \n  bvals = /home/arokem/source/vistadata/diffusion/sampleData/raw/dwi.bvals\n\n\ndwi = \n\n     name: 'dwi'\n     type: 'dwi'\n    nifti: [1x1 struct]\n    bvecs: [87x3 double]\n    bvals: [87x1 double]\n    files: [1x1 struct]\n\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## Calculate the ADC from the diffusion-weighted data\n\nExtract the non-diffusion bvecs and bvals (b ~= 0):"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nbvecs = dwiGet(dwi,'diffusion bvecs');\nbvals = dwiGet(dwi,'diffusion bvals');",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Pick a coordinate in the volume:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\n% coords = [47 54 43];  % isotropic\ncoords = [44 54 43];  % anisotropic",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "We note that it is in principle possible to pick multiple voxels and extract the data from them:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\n% coords = [47 54 43; 44 54 43]; % Both",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "We extract the non diffusion-weighted image, which is used for normalization:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nS0 = dwiGet(dwi,'b0 image',coords);",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Diffusion data are nCoords x nImages, excludes the b=0 measures"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\ndSig = dwiGet(dwi,'diffusion data image', coords);",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "$dSig = S_0  exp^{-b * ADC}$ \n\n$ADC = \\vec{b}Q\\vec{b}^t$"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nADC = - diag( (bvals).^-1 )*log(dSig(:)/S0);  % um2/ms                          \nmrvNewGraphWin;                                                                  \nplot(ADC)                                                                       \nxlabel('Direction list')                                                        \nylabel('ADC (diffusion weighted)')                                              ",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAARvklEQVR4nO3d25ab\nyhmFUcjw+78yuSAmGNFqHTj8q2rOiwxv79gbIVQfVSB6nKZpAIA0/7l7AwDgEwIGQCQBAyCSgAEQ\nScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAED\nIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQ6c/dG3CkcRznX0zTdO+WAHC21mZg\n0zRN07SUDIBWNRUwEy+AfjS1hDisVhF3/xGAF9WfErQWsHmPr7tV/z0YhmEcR9t5INt5oIiNHGzn\n0SLO/ttZQozY3QAcpZ0Z2PrejYgTHAC+0U7ABt0C6EnGauzHUpabAUqJGDzbuQYGQFcEDIBIAgZA\nJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGDAP8bx7i2A\n1wgYAJEEDIBIAgZAJAEDIJKAAf/nDg6CCBgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkY\nAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAz4yjj6Oc7cQ8AAiCRgVOfs/jLjOEzT\n3RsBLxMw4B/T5KSBDAIGQCQBAyCSgAEQScAAiCRglOZuAuAnf+7egLeNf4e06eGG33E12j3+WwBa\nkhew4W+cxnF8rJRuAXQiL2DPEzVPwtb/n/HfRSiFA3g0Bq7X5wVstjv9GvYmZ4oF8KvNUBnRs7yA\nPc6xFloF8KjVh4RF3oW4G6qI8wVyOb6gmrAZ2FypzY2I84LhNE1PblAEoDFhAXu+cqhbTZqfLeu9\nBTYilxABQMDgdy6AQUECBnzLjxDjFgIGQCQBoy73btC8cTR5/ZyAARBJwIBhMN8lkIABEEnAAIgk\nYAAta/gmEQEDIJKAARBJwACIJGAdaXgpHOiQgAEQScAAiCRgwP2sb5+q1R8XIGAARBIwivJoPuA5\nAetFkwsIQM8EDIBIAgZstXrNn8YIGPzC1bhXaN7H7LqPCRi8xCgD1QgYcDNnBnxGwACIJGCA63xE\nEjAAIgkYEMB1Mh4JGAHcAViWtUduJGDA/Zyj8AEBoyLn9cCvBAyASAIGQCQB64jLDEBLBAzgQ1+e\nEbrW+yUBA2hW240UMAAiCRhQnWu37BKwA/h0AVxPwABu5g7hzwgYcAyjMBcTMAAiCRgAkQSMckp9\nc6XUxlCQVdMbCRgAkf7cvQH/N/57JjOFnPc6/wK4RYmAzenaFGv3N4HDWSYlVImA7VZKugB4okTA\nxr1lOAGDHrw4/5u/ZFZqVJi3x0WEG5UI2NyqcRyXaO0mDQAW7kI8gLMwgOsJGACRSiwhzqZpWlYO\nXQCjQ9Wu8Vys4FUuiisUsEG3qM0IC6XUCpgZGMAZmjz9KnQNbLkLcb2WyCHaO3A5m1uTqK9QwGDQ\nWuBltQJm4gVsdHJOY8r7gUIBW3+L2TUwNr75eI+joQEaVChgc7dmpmIn6fYsr9sXfjH7mSuVuAtx\nyZVucZ4m78KCnpUI2OOzED/jLnye0zBoSYmALdYzsM8idFQLacamWBpGP5o/1GsF7MvqiBav0DBo\nQ62AHWIz/dpcV2smcobgC1y5k13/5V6JtyCUCNjuTRwflGb+45s/2EyxeqDKcJfNUBnRsxIBO7Ax\nF+fqrgE34dDibbesbTppuIs9/70SATvEfL7gRkTWjBHQsEIBe5yxvhUhxQLoSqEncQzDsDyJY0h7\nJv2VDyAwqwAYqgUMIILzyAoEDIBIha6BrdcMl1XEW7eIbM6RuYDvxd+oxAxscwPhEPIVhFv4qNAS\nxzPfKDEDM99iMJYFOvwtM5vhLSUCtvAtLgBeVGIJcbY8wzDrBvorOTmFhvlxoO8qFDCADSdtPFEr\nYCZeHOX5wNfPqW4nL5M+FQrYct3Lj6OEXP2cHHC7Wjdx6NYT1lIA1grNwMZxtIT4Oue5QOcKzcDm\n6Zc76QF4RaEZ2OAmjo4du0B63nLrSRNf68PwgUIzsLleJl67DHAAG4UCZgmRUho4afh1stjAa6Rn\nlhC74+4PoA2FZmBxX//anL16DunZ7OEP2Gk0rNAMLKteANyrUMAAOEoPM28BA/a5XEpxta6Brf/R\niiKw5noeG4UCNojWOep/5lO+xVzzv8v1vNdFWEIEIFKtGdh6FdFsDOjNqcuk7a3BFgqYYgHwOkuI\nQCFufeR1hWZggwchEsuYe4bG1rvWGn5pVyoUsPWjpOIeKwW8yyDOlywh/mIcnVwDVCRgUMiBV4Bu\nnN+4jsU1Ci0hTtPkGhgQp73b01OUmIHN3Vp/CazIDwZzUF7ATgY+U2IGNs+3zLq4Up2z5iKbAXFK\nzMA2xnEsMgMrwgDHGRxXpCsUsGUhschUzMebNrilglYVClgDjBTN+OD0xbsPFxOw3xmYAAqqFbB5\n/bDCKqL1Q8CZa3G1Ajbs3VJfk8Idwm4EPlbiNvrZetZVP2AV1LkRvHPeBbhFuRlYJwQaivMhra/Q\nDKzOrOvsE+rbX6ipG5Tl4/m6QgG7/cYNLuZTSnG7LXHc1lFiCTHlxg0YfK0CyqgyA3ts2F0TMqdX\nNMBhTA9KBGxuVYWvfwE0oJMzmEJLiHTlsw/Yr8t3nXxue3DvW+lAilBiBjZUWkLc5b4guIyrjLyo\nRMAsIQLwrhIBW9JVdgYGbFiT4HYlAuYnMvfG2EcKlw8qK3ETBxHG0ZWJ4xkcC/KmpCgxA9u9C9GE\nrJT5Iy1goR6nEbeP0bdvAA0oEbDHmzjcWF+KseYJO2eXlTcuYAnxYC9OU1JmM/OyYW/DUIcv+Scp\nByp9EjB+NI/jhw/lDefBcA9XKrGEOJumaVk5bPgCWMrwnbKdQLdKBGzu1jRNjz+U+eKSPRm171rT\n7+c/ChVccPD7fB2lRMAev8U8nJau3g6dD15vwYqfpLGXs2jyRd2o1eOkASUCNgtaNnQ0f8yuo0ki\nd4tCAXvL7oMTPYnqez6EQIq8gD3/iphufUO9oALzuRfl3Ua/uddjYxxHX4KmHw72Y8lGliozsPUD\n6b+ZRb3wUA+H5z4f3SKcfT9h53zpyQ5MPPWvErDF/G2wzxq2+6c2vxn4HsE+B/P1Gs7nw1AZcHiV\nWELcFGv9jea3/pJDNwoCNPzsj1Y7wYHKzcA+MPevkwd5UErD5+NQX4kZ2Ac2M7blF89v8fj6P9rO\n2e7mtVwzEBvueVG146TUZ7/OltyuRMA2a4Zf3scB8IGgEywNm5UI2PC3YbNT61XnTOqQLanzcgAu\nViVgbsGgiJTT8GU7G/sRdLzCuzkrcRPH5sHztzyH/kBPv2mRMTiGitu9d21w3I6CXVVmYLs3ZXAN\nwxnEeTIJ6+cTXSVgpLB2MdTeCf0MXiBgAEQqcQ1scBMHrbjlYX2eEHiBgtPu3fe9qyOhRMBc9LpR\n3OFusB4C3zWu0duBYQmRi/T20cJdBmfY7Lr1Tu5wr5aYga0Vf55hwWUEgA7rNdQJ2Poa2F3p6vMI\n4GMOGCroeVG9RMCe/vzJivo8Vqis51GMbt/3EtfAlgch3r0hhzlqpfHsIcmQl8i7doZqe9XVileU\nCNjw9yehzA1rqWQAnKRKwGbLD/TSMK7hPJfXlZqiMVQL2KLmLYiHM3pClj5GphhFA1bZx0dwtUX2\nwaexLc6H6I2A9Uu9HhU8ybidLlKWgHEFYTiKPQmLHgPmjPJLdmBlT/LmjTvVi7vXKciBegwYVHDG\nQGZkpCsCdhZnu+exb7/R6gzAIdEhAQPI0OrJx8cEDIBIAgYfuv5cuI0T8PNehVXE3ggY3KCNFMG9\nBKxlXY2SXb3Y2zW2txt7Of0QME7X5Ojg61ZwOwGDciQQXiFgN/t+qDLYAX0SsBNJC7uaXFOF6wnY\n/5w9phizAI4lYMAvPl5LsAjxDXvvVwIG/2OWnM6I/4qW9pKA8YmWPgM0wMlHnwSMSLkFNdTCUQSs\nLiMda20cD228ilvYdY8EjHP51JXSwHuRO/nmcJ0GzGeAbjXQMJh1GrDLKGV93iMIJWBwHQuqBV35\npjgAjiVgQDvMp7siYABEEjAAIgnY/Sx6dML1D6IEHKwCBkAkAbtCt6fe3b5w4AICBsdrb1m4vVd0\nPedzhxMwUmUNqQavHmQdk0+M4zAMAa9EwE7X6rDVzGd1VjwwxTfvRW28ilvYdbsEDGhKY6dWPCFg\ntMD5KU0S4+cEbBgMf5zMAXYqu7dbAnaFPj9dhhXgVAIG9Ouy0yznc2cQMDiXkQtOImAluFRLqxzb\nnKffgPlcAUTrN2AtEWNYu/ET4cN4JQGDE7kABucRsCo2J27pA1/69gP1CRjBLNcE6fmcpufXfqoG\nAzYa0qjhm2FLm3N57y7TVMDGcVQv2lPk/N24fJeLD4Aix9srmgrYNE1Tyo7fc+wAccFRaEQDbvTn\n7g043WZOFl04TuKggNVQOaUsZbUfMMXiFo67ns2LE1nHwDJUjuMwTRkNa2oJkSLiPrq055X1bQdq\nOgGrJe6qUtwGw8Vk8jwNBsyaIXA753YXaDBgAIOEdEDA+JZhghs1v0DX9qv7koAB72m+GaQQsHJM\naOAo936ayn6Wy27YuwTM6eTB7E8eOSq+Ydf9RMA4wO0/P9AnnLIcnOfpOmDNzKOBn9y+ish5ug5Y\nWQ56FuaXJ7FjGyBgRfloUVDoosXjZqtXGwSMY4QObST68mBTryeydo6AtSnrKKRtDkVOImAcSTip\naZm0OURbImANuusjalygOPVqTPs/0LI3PqKwy+eiPWZgHMnMj13eI84gYE0x/QJe0cZtwwLWDvUC\nuiJgAERyE0cjzL2A3piBARBJwKAu1zWfs386J2AARBIwACIJGPAGS3bUIWAADEPgNcXeAzZNYW8Y\nwCEaeBhH7wGDspxawXMCBkAkAQMgkoAB77G2SRECBkAkAQMgkoBBXRbr4AkBAyCSgAF0Kv27zAIG\nQCQBAyCSgAEQScAAiCRgQKS4n/1RXOL+FDAAIgkYAJEEDIBIAgZAJAED6Ff0wzgEDIBIAgZAJAED\nIJKAARBJwACIJGAARBIwACIJGJAn7rGznEHAAIgkYABdy30Yh4ABELkqK2BApMQBl2MJGACRBAyA\nSAIGQCQBA+hd6AVFASthDLmJ1XYeK2I7IzZysJ1d+nP3Brxtefunh3OG9ZHx+G8BaElYwMZxXMq0\n/vVCtwA60doS4jiOZugAPdiZxFT2fAa2/M76F9dvJEAD6tchbAnxOSuKAP1oZwnRZAugK2FLiMPe\nXYiPC4ZxLwqAd+UFDACGlpYQAehKUzdxbESsKG7uq5x/UWqDd9dsN79TwWaram7kovi69+MzAYpv\nZ+WDc3N5fpqmghu5mLet8pu+1mzAfv3K8+02h3XlDV5/OSFiO4e9S6R1LG99/Z05q7mdm6G27MG5\n2ZM1N3K2Oa8qu50LS4i3maap5jGxEbGRQ852lh0LNlKeCbBJQmUp2xmk2RkYx4r47EUMuCk2zwSo\nqf4UIc4y94r4NAkYv9gsJlRW//Er87at/7emiLc7S/3E7l6PL84SIr8r/sEbcj5v019D4b2asjOh\n+knBN+rfQjOUvwvx8QaqoeR2Du5CPFTE3X1D1HZubuWYf1FqI4e0D9HQdsAAaJglRAAiCRgAkQQM\ngEgCBkAkAYM3jP9afvPdv+SbP7j7p9z7TocEDN6z/i7Xl9/y/uwPunMYZp7EAR9aHrezPEl2+OEL\nNOt/3PyR3WeTPz4Mfu2np5vXf9YDHEvA4DCP07Ldp3r/9Dzy9a9/fSzWulWbvxM6IWBwit32HNWY\n4j9TCq4hYPChJ5Oenx6KetQ8qfLPlILLCBi859f7/R6nR7tXuZZH+n4wl9r9aRdKRm8c8QBEchs9\nAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQS\nMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIj0X9Eyy8LQJ/pIAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "In practice, you don't need to compute as above, but use the following call:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nADC = dwiGet(dwi,'adc data image',coords);",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab \ndwiPlot(dwi,'adc',ADC);",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO2d23bl\nqA5FnTP6/38558FdbpevGJBYEnM+1Egle9u6gIQA45/f398FAAAgGv8bLQAAAEANJDAAAAgJCQwA\nAEJCAgMAgJCQwAAAICQkMAAACAkJDAAAQkICAwCAkJDAAAAgJCQwAAAICQkMAABCQgIDAICQkMAA\nACAkJDAAAAgJCQwAAEJCAgMAgJCQwCA5Pz8/h//uOf++8Jrn7/YSuOTur78p/OLD7z01AqiDBAaZ\n+fn5+f39PcTi3z8sf8L0+rGVxIH70+vXc5sCcvDPaAEAhrHG6DV77X+5/nD4fQlbxN9f5O7W+1vs\nf1j+Th6vlzon6fNXLu9yKfn+8+tlvxoBwA0SGKSlMfg+f3cL7tu/W+7Z3/qQS0rk2bLR/lJ7XZ6T\n4t1XDrc+X5xcBeEggUFmtlj/KTq3hPLzktv282VdVX4pa76KBzAcEhhk5rlqOc/gLW3Z63C7w2U/\n3eW13upOLyMAuEECg5w8hOBzSrhbcyqcRXy+yOVdHoS5vFTJloqHr7xmo8vdLgDiMM4CiI3dChZ1\nGIjDNnqA2FA5wbQwwgIAgJBQgQEAQEhIYAAAEBISGAAAhIQEBgAAISGBAQBASEhgAKlwfrcLwEBI\nYADhuTyKHiA9PAcGkBlO04DE0LgBACAkTCECxKB9cYvlMUgGFRiALnUTgEwbwiRQgQEIcSiS7PIQ\n1RgkgJEawGDWREJPBPgKCQxgAMzyAbRDLwIAgJCwBgbggeeaE4tbMAlUYABWRF/cii4/pIcEBtCT\ndXGLJS4AB+hmAE3MU6aQlUENWiTAZ8STFpkGJoFNHABFHJ4vnjxDsE8EFCCBAfRnbHx3SK6T528Q\ngQQGcI3bqU7RoRqDUTBXDvAf4otbALCHCgxmJ9/i1sCSiGoMPKECg+mgzALIARUYTIFzmUUhsmAE\nsIcKDNLC41A64AuwgAoMcjLz8F9NdzV5IA0kMMgDG981OczZks+gFyQwiM3d4tbY7JX+QWYABUhg\nEA/KrAMP+VKw3DlUY4ISQhRYWYUAsPH9QCyDsIMDjKACA1EaN77nG9fvi5Vng9zpLl7uiIsHgjAy\nAiEYqh/AIAAPUIEB9Ke6knDbSCle7oiLByIwvoORxFrLAQApqMDAG59TnQKN37tXG1+vZl3uNF6c\nagzuoAIDD/zXcsRXj0zFE9cdoBc0dDBh8rnBcwqJYhBxOcXFA2dIYNCNQ9SeuQ5YdSfaApgyb4iB\nLsycpS7BIGfsbIK1J4dNHJCTUcv+CvHUSHfBnRScETw5JDD4RpQT32UFc0BNdzV5IA0kMHinbuN7\n4hHx68bulLqLlzvi4oEFJDC4JkSZ5cmnLJ7eYukVhBCQwOBfus8NRn8jV5TJ0jPO9cfr7fzrIaqx\nSWAPz9Qk3sRVp1qOje+J3Qqwh4YOs2MR7qdNIeKKi4sHX2EKcTrcZlSiTN0Q0Toibkxx8eArjEfy\nk2NaDMKxlTuCdQ+dIgdUYDlZ++fWSyfsqIf6L0o52AURZcVbnbh4UILcyAiquTxAdmb/Tq6+FOK+\noCALinSrgmfEe92QV6gsGgYRj9cAOWAKMRg+b4MMhKZBBoqh//pHffHEJYQNEhhYYfrCxte7TBWD\nHgyiYIeDDAqDjAez6AyD4BUSWACqz4NQCF4WEF9inRIiKF75eZ5ZO1EOmKlXRGctB+owejh6aWsV\ntCtIBglMhXzL/o0ahTZIaOEBosAU4jBizQJVUH0UYfXXdeglvN38lefkWMWNrMVjYjAHjBNdcZ7D\n0a8D7Ayir/slXQyipruaPIukSFABXjRnYFfR7KWaUg1k+NLUcAGeQTy4g1BizrTxeq/4bEZ41Xc2\ngwBYwBpYf/bT98SpldmMcNY3ypKn+MZxcfEWVtd8Ibz2gWkENRSGDqNahYLue9TkWSRFggqowOrR\nPMRoj9qr5WdDs1WU06vc8TmT5RM+fqFHWMMw5AOUWWcYyR5YDYJZABwIXIH9/MHhLuvPFQPqfEOw\n8rWcfLpfEmVxqx1xhyLehAQeJ26D3O6j3Y4XzDESz6FFR6jFPyHYfgRFggoCV2Ah2t9YIVsGfdGr\niu4DXv0lzw2Hwb54PYF4kxB+GNJlJJV1QM0ws4WsraI74s1MXDxoIXAFtq49nN9+VDi6aVzcCsHX\nd6/0HRhGH2a2tIrouoMbNJUWAiew5SpAr0HnLhb7zwKJt85A02KxmMqSe2UFG3wI8VbBBMUTJ2oC\n2/x9mav2aWzsWo5gIHMzSJT1v3xRY6BGl05fe6GUmfM5fU6Szw5PPv29PZO0jE4n/jy4Pr1BRJr9\nbp/wv78REOo/RKwELSR34bRtdFrF73A2CPY/sOYwTAJ9iTqFGAXPmQq1je/DZ2nUDDI3P8oecDgS\noQVx8QaSfJw4fCA8XICpSD83GBTBXqAmkpo8UaACs8W0Ub6Oy2YYtd1tpJxB90v0R+v64vlLWB4o\n9P3rSfK0n29c80mjfOqfmUHHVzACzAkVWACq13KGT586XFkzcDuPkeMeqSxeTyCeOMkHbsNHpnUC\n5FjL6Wj8HAYBuGN4pAoKFZgtn6b79t9K0JobVWg0SL6Raflw++5jowbs5WK7iVdxI/FyR1w8I5Kn\n/cnHNeHU71u0xdL9kuoKPoHuAK8kb+j0ZHGyzg1WN7ysBvmKuB3ExVumCX3JlRzuxYMA+u2+I3fG\nH+4UB75uFl2maRUAHWENzJx8i1uFzKPpmVfdp2oVgmszaiKZyqOmbEdIYP3Zr6aOjU06Ddd/47uO\n7pdYP+He8l1l04mLt9Qa3+29EOLW+8o/owVIwn4WKPdoupDDtBgvsvGcJ2y5i4WEHa/Z61JqzQPq\noAKrR38WyFkk/eeL3fjZva9u/Y1mC3kma7ljini5Iy7eV0hg9ZTEowRN5BmS1oFDLR7aJvrypxEv\nfaAwIkwCOzv4Z8cQkUoQ72AVHAz+oODw/Z9uN1LL4ka6V1/Wxxef7uIfNATDVIJqLMCG5v0kzOH3\nJXu9upwHoW+lS3ptWI+48d1UZvFWEdFfABUEqMAe5jGMyi/9xa1C4kreTnfdP7WKseNZB7+LD9gR\nrwVx8fYESGAPrHGk3dx2s0CBmsKBdoPE1R1eER8bIV4L4uLtCZzA7k55KPx69B1i3emexaObtMUg\n0XX/xGunGzuUuby7zuhq27A6WpAjgiKdiZrA7oy7FmR3U4v+c4PiGxnSTJaeqet+OYY1zqFH3Eoh\nxBMUUlCkM2EWe7d16f0P658eVsge/joDl4v5rPAfmMcgPz+LtaLixkS8ZCS3Fw1iGZfIZY3vYBBB\n3beqzE0uQSNAMqJOIUZBYR457lRYR5wnS2UNrirXABT65sbdc65DhNkEGHj3QkhgtowKZAqNb/j6\nn9rzxW5cev/31zt7bTYfHovP6FeHjDtLIIHl4Ryv1aIGzIlgIFYTSU2eRVKkMySwwLwWGSGaYEd0\n9hBqDh1+fpZRcgU6tUiwXoQ7SGC2dO8JgTa++5+AJ24QN8SNoC+euISwQQILQNC1HOujCK3v0kLz\nCZzfSqWHscLfthqzieNZvCHlTuFN3cRrf2yxOyHKUPWVzEb0l2pNSaN+xcb30LqvoeNB/GeDqOmu\nJs8iKRJUkNyLkzfT0OqHFt6C9sfXxB/tR7wWxMUzInmMCBQEA4laSIVGaYzwqkjhoRhpDAJgAWtg\ntpTP/hOnMvGavbZ/T39yahXie+1mE890t9eoKziQfHynNoCds8x/wM4gaq7fY32qk5ruavIseiKp\nybNIinSGCsycgfu8NcdQbHxfdwNueg+vNqwFaN6T2V+8jq2ui3jWJ3NWfCtExwyQY1sYMoigzDoT\nYjTXkYI1sLkMAmABFVh/9lWFZg3kACt8B6IYZHg5+Iy4eIt8lxcX7yvJh4HDx7nDBfBkKmVLoBZf\nsWsY1VdWa6tq8iySIp2hAqunZCwzPH163kKquY8aaSqs8CUbZcel3BGC7USqO98RPoF59tUos0Cw\njPOOQquwK3fqvmhnk+ozgn3cJH6Esbh4JQROYD6z4fu7hNsyZ7Q9vTCLB+0SXZhZdwA3Aicwu3TS\ncRYoRyBTmBaT4nXwlNJE4gP20OINF/gs3nCRSgicwDpiNzcYN5C1GyTZ+l+gLO4Qer4erGwnySWv\n4o2NzsqNJxbZElj1I3spm9Qna7DCdwCDFPLazASHMjreFNyqugojJdIdATZKPnPe67lvECF2gg4E\n+xzoFU3KDWt9rBRAYrJVYMufcupylcJ/3kB8HlnwABt/As0N6iPudMRLRvIB+PDyXLDEERQpB69v\noXz44iiH+NxasMmpiaQmzyIp0pmEFdiBsT5QmP2fcy0n0GD293f5+bl+u0odhbo/vNXFU4xRSIl3\n7pXDD80KESgC5NgWQgwiLFgVn1b9M+KmqK7e2u8rbBWAF/JXYPNwLrMU6j8RAZxN8VX3/dtVPHG7\n6fDG8MzwcucZcfEGIj0sbWf4uNtUgOErfM/4G1/cIAB3DI9UZwRFOkMFZkv3FhBoy5yPbHcGcTlm\nzPoONQiO1g/yKIj3IIOgAfeIi+cMCSwzWRu6wp4U6x0Qn3gwiEIbOIikMOp6kGHIuLD8jm7iKbjp\nFRJYAKrjtVoT7LXXrkQva92fL2+dNp53lqr5HYygGiOB2VLXvHJsfL/TPU2Pc955sW8V5cPwwyGt\nk8e7ZIgvIjhAAlMh0OJWIWcV1ohfqFnoUGvhvvZWYdSuXgtrQVeqiaQmzyIp0hkSmC0P8cKhzArR\nBPd0NEg43S+pruBLPqNsotnEM93tVUeIMXSAjZItqO0Edd7nrab+mbEb30c9PvwATwIAlEMFZs7A\nucHhz8C9/n7gZKnp+P7T+Lducasjs5U7fREXb8kyG3HJP6MFSE7ipvOJfS2Yqby4exnKq45SBhku\nwDMziGc6WVJ3Zf35m4UKzILhA2oR9ooPqSqeP/BpR0k7njtLW4ZNFvVExwv2Es9uZCk+ZhUX7ysB\ncmwLboMIzaWLEGMoIxxdf5sFR7WKBH7nPZ9QAhVYPSJrOQ+In2BrSoXudc9ZX95nS11DWkWvmw5f\n3XnVQ6rJnUE8a0hgH6iYBRrYRPq+Yur+LtKHGFXQKHVQrS8ZOCwrnN0VHDXuQTxrAiSwnz/c/d40\nZERf0Opum/IsHs5WdatilwYJp3sJ1R3NJ6l/3flpJ0nFHYePey4D7BBJPqE+V76fzT/M7JdM9Nct\nBiRYQljp9brCiAYxlVlzyXNDzV9q8iySIkEFASqwB4bP0cvTtBtt+zliV//9/f06ifr8ef0lzw0H\n2T71O8k9qIqPSIogLt6e2AlsjSPt5vbc4uxJiyLtRgjUDS4pbxVr5turG133Vw5zIQMluQTxWggU\nAAMnsEsrl9dkPotbgq3zjpTl7CevHhbAWlpFoBDwQLrmoItg1xMU6UzUBHZn3K0mu/yA/yxQjzrG\nMI6YGmT4LI3b/Z2fiX6l06O+//37zGs9MXov/oV4OtFZvBoTJ8BK5r5OWnarr4ffv34xIo2Pc16u\nVE+yfF2uZj6DdNFI8KTjdsR9LS6eIMntlaBBtO8kFN8y54+DQXI0vKUtgYkbQV+8hW77RtQpxCi0\nzwlU1167K4zZMic1H6JgkIi0+PB128tYxNsArbQEElge1DZSDpdhoEGG697OqkGv5wj3vhieupbT\n6Gr4DqbXx5zVnrwWQbqIbkd8lsCa2dTfz7rMpvseNd3P8vR6xL4aNRNBHVRggRk+bFTgbm5w2vAk\n2CTOh2QOd87j9KZ0txIXz5nkw5B846xAGvmIqmmQCqlaihJNIwBYQwVmS6cncrQWt4ajb5A6qaob\ny3M9UXlRF8TrCTfx6u5iKp6yXzZIYLbUBrIMR1t1FLvCICG6nxFqjefBF+JbZMW3AoqL5wAJTAWL\nfd6hg/hsG99fz/JoH26PKnfKnyh3E6+iOZmK1+Xo0b7ihehxyafOWRuIRVB/nZ/57aVIUIMA+EAF\nZkvoGqiREt1zTJZ2BIOssDbWiLh4vUg+vhs+gD0IwPEwi4BT7Kg7fil9qxD0uJpIpvLUXVzNRJdQ\ngdny+/frypzXciYZhV0irrtpqxDXHc4IpgpBkc6QwPrDLNAZf4MMMXv5PU3Fa7m4xeSY3X7UakzL\nHaMrV3OeBErDP6MFSMJ+FkgnYw2U5DAtpmMTa1ZFz9MvUeYJxSUUFw+coQKrp2QWKNl455nodefh\niPSWE9O38xg3m0R8EoCtChXUlTtDnh94vqmgbc+QwOopiUfhYtZXHpJWj0db6lNIdfc7fO/rZYY/\nvtbxTSX6zyNWi+cTncW7v7h4JQROYD9/GC2ILhbGSbzCd3iBSOH7RDQNYtEtqpuTj00+iTdiYPEk\n3vBmcxZvuEglBNgoecl+jeFhu2f7TtAoSxemdNxQu3WTkuvJvtX+tVWE2IIMEJ3AFZgdHWeBZi4Q\n23V/PV3Jk+Fzg1KIN2zEa0FcvD0ksGVRnQUaS7tBLr+15qT0Bk7fhD4pWBgQe63ereVxF5GMOFtP\nKmdsW5BGC/JOtgT2aduPww6x4eeAfPpwpize3v3iGsQ59EhFurMsUuKdkV2kEBTpTNSZ+oc1sH1a\nOv8pqL51vOo7m0FekY0mUI7s0il0J3D82ieqhw/EVdAIt6QVJTtayNlyzZZXMw9H3OmIl4zk9hre\nIIYLsIxL5Aq6X6I8svm0S1Mc2QYgiKCtBEU6k20NTA2FFjBqy9zwocPdf9PvIVRb9VGTZxlxyMi6\nRaX8njzkWkKAHNtCiEGEBavi06q/aJdZ1gj6XVAkZzLV1jpQgeXhvGVu5pAhVf85M7PfZZnkARJn\nSGC2mAYy8X3e/kHc2iBfZ4FGITj1FOiNHszdBYIEFoxAazn+J+BZG0TkdV/idz9zyAcK4j2kqEFH\nMMulTEGRzpDAbOl+mLfFCQjhKDGIke6as0APO1bE24B4ueMmXl2gMBVPYZzxSvKV1cmXjtOoP/OO\njEswCMBCBZab0BsZAk2WnrEYF7cffka50wLiCZJkhH7H8BKkXIDhonZn06j8aJ/QRrDYJ31nkNCn\ndcAZwZYvKNIZKjBbnl8ZpbyHsJ3hQ4eBd6+mpFWsf3/Qr0T35wF7382WFb6wrifapwfEG5i4eL0I\nkGNbUBtEOC9dqKl/xs4g+rq38FyBtevOebgQgsydfNGIYgoySDGDQT7pOHxHxnABnkG8RhL3uLSK\nrQzx3L5BJ246r+x1n80Or/rOZhAAC1gD60+4LXO92B9UIfjsqidnfaMseVqs7nS8YC/x7JaI6q7s\ntmRVfqMQq2jJh4GTj3Od1Zc6rlTB9Z6TS/tVKwXd99jJU31lNRNBHVRg9YQYoVjzd1UhdFDFqPA0\n9vG1P2lMwwd/UJNnGXHOmSDi4pUQIIH9/OHu957Vd4hZoA02svuzrX263/e/f40QdKj4GcGIZ80/\nowV44XkjgEOA3s8C1Z19oJ/nPnGYFkumXS+GmAVXxEUwUAiKdCZABfaAUfnVcRZIvwWUUGeQxPWf\n+HOsRrJVD9h9bPVJvBElsnS5Iy7eHbET2BpM280dbm6wkBbLpDRIC2MXt9T4pL6/rV7vONZ957sL\nNidBkc5oJbCfv3n+8M0pOzUhO2s8+qRUsizeLn+LQcaOYd13jkgP2BHvK4H6vtYa2KfDCy4/vBVk\nIj4Qn0fei9ddTnHdL2lc8ry64LJkX50S9zLiJUarAjuzJqSV/f6u8+/3X9n+eria4GBnLIWdZ/+E\ncveLD8dibjCK7t257GI6/W6TREeklVUeKamkhLkj3hj5E1LVmD+XNVBdYRSukohY/+1pMbiz7q+i\nCvpCUCSoQL0Ca4dm2mVxS+cJ5UJm9ruI7lurE5FnT5fNX6aI73QVIfkwZNpxloLiCjJstAhTcUSW\nmO7WDzgvi1l1Hq70B0/yV2Bj8RxDZdpD2IWBBtGx/+vbL9spqc4bO4J1Nxpe7jzffYh4Ieo/Epgt\npoFMfOO7vzx7g3SsgaTOeCzh71YxUJB/afGF0RFZh+gs/iCNuHgDEZrosEBqJqcLnzTKp/6ZOx1n\n0P0S/Y1LB9eoTRKKG1BcPGeSd/J8UaxCI6m3nHRB363u+wBvbxfBVv/+oC0mKMIUoi1d5pG7zBP6\nRwe7OXTxiOzAp9nj4Qs8z4SbodVB3LMOqI/OGhk+/Kx96IpZAjiybxV1VcuodjW8G55RE0lNnkVS\npDNUYLZ8Wq/af6vTOUaBR2fN+9YC635Je6vYn2UDUgimCkGRzgTIsS0oDyIchsPK6l8STmBr7Awi\nXuUjXiOTdKXkSqp5Ub/dd6TE+FkNUt3wnA2i1kEAPsEUoi3nxXbixQwG+aTUvpE4G2SbVFSeV5xB\nPFMF6y6ubPMNEpgtKaNzIX8/61O/kbLiIPxFvvtJZfHtBQ5G12+8soV9Oio73H2viIvXAgmsPzoH\nZOgE8cRdqBCdVnGHXTXWUdle4tnZv048z0njwk8KNtEzyWfA3ab4s67lQB1rwzu0CrUjJx5QWxsL\nZDrwhAqsnpJZIJ0ayIGDslPpfknHySVnW2qujZ1lkRLvDOJZozXO6k7fgWRFmaU2ku1OegW/shmk\nY1H+8NiyQ2nS7mLBRqImkpo8i6RIZ8JXYNaDiMYdYsPPAbG4pvhazorb6PLOIB3Lr7E23k2BahVk\nB5RlW97EG96PzuINF6mEADn2jn1eefhM3UlOcc2yp5ciaQzSi6/Vld0SqfPiEC0BpAhcgenvXh1O\ni33ayyzxEfFXWja+27XVy2Nw6x48KLsd1VgHEK8XgRNYR+ymxQI1hQMMDqJMlt5h/Bbm64f8RAgn\nnpSQ+0Vccf4ZLUARB1M+zxk+f+Dyk3axKVDU6z7NNXz9r1GAQL478Pu7/PyYzyseonD3mXxTBEWC\nCsJ78bIhbrH48NfZWu2rvokNklg1ZcTNjnjJyDmFuD8aJ/QsUCOvj6YlNsjXHRbDEV9VKmTtdwkU\ngRAkT/jDRzTDBVjGnRLioHvdHjyOTfFBofE/oC/eQit9I2cFBi1b5nrhk72KPzzeIJ9IUMRs1Zig\nLmrZ6/IxrOGD74F3L0TLi91Ra6amnJVNr/7DVgUGsH15OA3k4fPrhwV9cV4dX8Qk3CMu3kBi7EKE\nEiZs3weNDztLQwwhS+g4EPmah2rv8rMs/92g+9la7RzEEJHqDnHxBpJ+hJ5ZQamIcMbH+Lld3B27\nBFbuCHGX6XerRVg8Z1gDs6V7ERBuLUeETysxypVbY4taT+4Yfrii5sLYypBuVW4NN/FkHbSHBGZL\nr6MI+17QBztRKwxS3u3Xa8t23j/TcRISVrfMzR2jMlnhTZUT7SIvngPStXw74pMVd/SaJQiq/pm+\n0yavV9N/feKfqOXtX4v5K+bEoJokAe6OQBE8kKg+YBAR3ByxFROafhdPtOLiGcEUoi2fprZNJfGn\nYnLDf7LUfxKm1+zf4SJ9tXBzxN4a4mu6iCcICQysKOlOh/zh3wODdvtztjo/1fTxgoMdcbi19RPQ\njVfuLpvpbq9RV3CABGbL+cliz2ah2QRlN1LKLomvZUp30RQccXdPqYZxprt43ZVtHNAsQeaEeJDZ\nnP0SQog2MS0+3ulyk+eLHIKXeKt7kM5OEXGbQCFUYP1Rm40Zxd4IIgYJSq+Ht3RaZgVjt92XIyhe\nezUmCwmsD/t+tZ9eSNZcyvn5+102ZethEs82XRLaj6GFP7B2rkwarQhqJCjSGRJYPSVLCLEGuY20\nrKlE6CxhuHOE+MJSOe0DRLvoLF7uiIv3FfX58Ub0FwBMcVC//BZpfCGoSMUzQGkeGxJ0B7gRxvfn\nZloym0/j7k6awBcdHHFAvLMjngUBhL5bRymxeLtXGsNE0GZxoLsWaYKvs3+NbvfsDv2DtTZydDco\nJ8Aa2MPEvdGWJIXnY3KTxqpptEijyPrDa1gYu/xzvrvUctRhB5YyARLYA722JNltL9bZyN7yXVMt\n9DdGF2I0lnLe+H52x/DXr1TwK//GltffDGQVRkqkO7Qq7lPPKXqMcf+nw8eev7UEcVI1VwuHt8Fo\nBoOEAEf0RXxeUVw8ccLY7iEznf+0/rCOwmKdSmDK+W28agZJE7u/GlbNESuZ3LGkUAQOKHabS7Ye\nvv9h/ZNyjaUWmNbJbSmR5kShcU6IWn/cWAWTEk9KmDsCiNhCCB9YEzRWBhX7zGHIFVejNL1J0BGC\nthUU6UwAEVsI4QM7Jlcf7AjdtATLnTOCWVaQ2LsQ4Y6t9ctuxPqE8o6yZw5iB9UiAefNnOK5Ic2z\nDaZIj0HaGT7I8hRguLIPnPePtF1NWdPpBs6a7vjqiAkdlwAqsMB8fUho4PA/d1jYO+J14Gz0dkpY\n2o4g2D4/pNwXLM0FRTqjOHTqiObYsIVeGqWxzChFqu/btxhVw98dFnekGosCb2QOwL6L0qlGcQhq\n1Y7I7UCfs0KWHo544HAklWynExfPgSTD8DvS1BlQQnd3OwSI89koaaJSR3eM7ciEEVlYA7MlxDxy\nGgyijO1OsLV1HNpImu1nObRY/ihiujbWfuXu4oWIXSQwW75ughpO3A3rB6q12L54mV0seGgjM7vj\n0wYlB9aBhaw70ox7PpG8NFau/dPMFEXnzhG93oOVe9dGRwL1CH1RlUNfR5IrqeZF/Xa/J5a0D5yb\ngadqvRJYJncsu3m5uBqFFj4HyR0wvIWlCToJGN4YYCWZI/T7eDKD72ENrD9p1i32RFTq8kHvcFrc\nEUiRhyfuI7arA13WxkyNUPvAYgC/pM3MK4mHHnAJHhdhWkeIF2TJ/EIF5kqIQc0rCUbNK3stQp/w\nlMYdCRQRTw/i4n0lVTY+YzrcSDaWiUu/47WWhb2CDdAjDlCNWRO4Avv5g/Mdt/82+j7BYHNlyMmn\n3R8S+vn5WZaf4N15WXzdYeGITBxW+x4+OSQaPIsXIkDFPgvRYSfufgzV9y5pervXZnTbAyFTusOi\na1g7IhnnJIHROhK4ArNOWttd3M4nTYDdYTae3T7HYkxHKANn8cMAAAmtSURBVLPa+f3DctVNhls1\naKINPwf6PMYsHIEGclhucIQdn6oxWUfICvYVZUUCrY3FEPQwYNmPYs4nLFhPoVij3LI/UT56SKCs\nPq92xhFDwOwtxFgDu3Pw+ffbE4Vx20RcyQ9EH0mcCd208rljSaGIbMgKYdsAIl5yWZNdfiB0v9Vs\n2RUkqIxTgiMU2Lo57vhKZntd5jAYwj4T5+ilQccWZ7GDKnImVru6M3sad/gQyeUlXG6Xok0MIVZA\nSQyOECHHPhopJmrZFGQOfI2VaWKrlCLVsS9N0BRxR6M9x7pDxIbPBBCxL6SxvqQJedHBESJYOGKI\nc0lg0qTp8P6KGLXsEB2mhByKZOogCbRYEinSkagW6VJIUY3VMbYjceTuBhFNBE9HELX2xD5K6rft\nPXLb4S6Zjg4yUsT5MKE0Hul+ttaQU51wx/k6o47X8jncLgrhR3Adxz5p5ky6gDVEwBFGVGw4WsQc\nYVqNhajvA4j4gNFh2+sPoS2zkmNPoGDgqKPcvJqOWJnKHcqO2AghpBFRNXfoRWk66jMzt34d4ja2\nZO0nriNWkrnjlajauvkpeoPe2CwWWqPQwu9JE2jSKJKMSfwSUsnXgxCN7hjRViuTtOagpPFODkWi\nd/aNRneE8GaM0+gP+Js1Qe2yJ0TTfCWiO3JYPgEPjkjjoP1ZemmUOpBWMTuixM3ErTYQX1tLlNb1\nilrzqzMs7hAnp1YOCG5WrO6iOirkIE3Uiw6OOJOsv6dSZgiZDtxM07j9FbF7oiOBR5wPqjC6F+4Q\nJI8mA3Guxky7aI724KBImnAWHRxRzXM3CRENAogYCKO+RBdtIUQ/fCVNG8AdUkRXJENjUiN6mwCA\n2Qg6sAh8mK8s+wOCqw8PlTo7deazXAUVn9Mdyiory1ZI0DPNQ2bdQHxaHgs6CEoGBbQIOGIUayAK\nEY4CiJiDy94YsYtGlPmSff+Mq1RcyQ8ECprPRFfB/5yjFmLbOhxpwk0CogeaNCR2RBTV9kkr1iA7\nhn2jc5hIFG8T5cRS5EHaKIHmlRCKvDabWO3qAXF3lC9wyHpE2r4JeA6ad38CSIZ4KHdAxAIPxVZE\nQh7mG4iHJrKvxp4/GYI0+Vgk0LSTRhFoJNay1ido4iqkSQDDIXCfGWITHHGHj2WSFVuX0MI60/4O\nniV+a3PWwu52Odyx+EbMBOaKi1vSEpk6IoF1o6NHCQSvMLqvJscB0Gn6SLsBh1Ram9hjeyJRoDMd\n3UkXNbrO5HSJmDhiOCo1EAksE93dmSaNrXyyTzLdBSl0h7gjxMUr59UdUstaChmUBNYZo/GIQlvp\nQkkXja5jDnCECFJJ68zYdsI2+hgkOPRoZS95jhAZ3SMbOdyxZFFEOW/pWFhFjjQ4uDZfNabTHyYH\nR4zlnLRkh0ciUYj2GhjZxv3KWfIcoTOcR+4EzuGOJYIiJZWWSLYQRN278EqUoKkfSiYBRwxHeXow\nFjTlJGiO0epiZZoIq6NIiyRRRkivDHeHZicNjUoHgy4M7yFpgl10cIQIFFumkMBykiZ+DR8198Jn\nd4/PBqIEHjG1FUnLjSTRAS7xCTdpckx0cEQdHU+K2X72PKt3ZqfzHFhm7Hbi7vt8gkN7nel73pi/\nI7Zb+99UjbHFFk+hzKv5bDSGG6KVCDjClJJkoDZDSAKDWfgU/pQ7Rpo4XhgxE2gaggdTa07ZTd42\nplZ+Wh664uT9YSx74+OI4ewrrUUsby2JxnAt0EmmJkcfQAvoyCFvLapOYYizkMCmJUovnZY04SmK\nIs/TEnd/GoV4dehGjLYFvXjtpdHbQywtnldcomgRl6/bMXCKGvgD/iJWAohIRRBMEzcVFFHbQwgt\njG9PIAhprCMYczgkrayQwOAWwan/FjyH/3b3SpMOeXMetEMCS0uv3ksUKERhfgz8iy38PhBMnxOL\nJ4rSjP2XLEEnjUca3TFqhpCx3XAydGM4Y/RILD0WdBBpjTkGQ0HhMF/4wD4pHn4Tl08BSDZaTVKN\nsR0D9pDAoIbtGOwEcfNZ+IOCsprKCvaVh+eIL/8KM0MCg3oOBVmC4LIO//fqRFQqhzsyVflghOh8\nCLTj3/9Dx03ZucGpuCy2xDMZLWcgmB46EyWNvcqZIzDpu6N8hlBfF3AmQxcFQYg18ADLWtAFEhgY\nkiaN5ajGltGKkLegL0m6JSgzdg0jTe7pjo9lSFpgB30bnHBLY9ZlX5qMaH1g40oOW4EmSboiBMIi\nwaSZqxxFl2Q2auOr5x1BChIYjKE95aSphEIzsNgyOi8NAsGDzDCG7SyP/X+fiXIoRnRekwEzhCAC\nCQxGsj+S6jUUEit9uLMz83WgBgkMxlNRjYmTZk2OYguUIYGBCtshhAmif2jhD0lr/18AKVj5hMHc\njfETpLGVKIq8VsCCimSq2qECEhhIIxg0M1ExQ8h+P9CBtggBSJPGFBRhWQvSQAKDMChE/7gw2wb5\nIIGBN41zUJkCsfV0HMUW5IYEBn70zT0UZJcoJC3WycAHttGDH4fnvTpeLUG4zFGYsucePCGBQWzS\npLEK4RWKrQN9xygAz5DAwIRDCLMOr2nS2MpDNSaYtABGQQIDE4bE1jRHUp0f6L7801hkBYN5IIFB\nNvYHBC+RA6t4hlCTByaEBAbe+AS+/Zui3G7azsMMYSxFABxgtytMgXL0/7SspawIgDMkMJgIqeif\nYK0OYCwkMJiOgWmMPYQAHSGBwaS4FUAkLQAjSGAwO0YF2agZQmYmYR7YhQiz0/EhaJFia9OIHAa5\nIYEBLEtDGhNJWlIyAPjAGA3gSEkaE5+po/yCGaCVA1xzTmNqxdblUR1SjwoAmEICA3hCLWm9Qu0F\n80BbBzhymbRCVDbixycC9IUEBvAXzxVMiDQGMAkkMIDPME0HoAD9EAAAQsJzYDAv4lvhAeCZ/40W\nAGAkv7+/v7+/h70PABACKjCYF6PCi8IOwAcqMJgdix0ZFHYADlCBwRR4HlpB4QXgAwkMpuAuqVjP\nIgKAHWyjh0npcmjF80V4XAzAFCowmJQuqeV8EZIWgBt0NoDOsAsRwAcSGAAAhIRt9AAAEBISGAAA\nhIQEBgAAISGBAQBASEhgAAAQEhIYAACEhAQGAAAhIYEBAEBISGAAABASEhgAAISEBAYAACEhgQEA\nQEhIYAAAEBISGAAAhIQEBgAAISGBAQBASEhgAAAQEhIYAACEhAQGAAAhIYEBAEBISGAAABASEhgA\nAISEBAYAACEhgQEAQEhIYAAAEBISGAAAhIQEBgAAISGBAQBASEhgAAAQEhIYAACEhAQGAAAhIYEB\nAEBISGAAABCS/wNUyNQl+TBNqwAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Solve for the quadratic, Q, that predicts the ADC values.\n\nThe ADC values are measured.  Each measurement has a direction and amplitude, $m = b \\vec{g}$, where $\\vec{g}$ is the direction of the applied diffusion gradients.  We want to find Q such that:\n\n$\\vec{b}t Q \\vec{b} = ADC$                                               \n\nWe express this as a quadratic equation.  Suppose the entries of Q are (sorry about this) $Q_{ij}$. Suppose that for each direction, $\\vec{g}$, we have a particular level of b.  Then we have for each direction:                      \n                                                                               \n$ADC(:) = Q_{11} \\vec{b(:,1)}^2 + ... 2 Q_{ij} \\vec{b(:,i)} \\vec{b(:,j)} + ... Q_{33} \\vec{b(:,3)}^2$\n\nThe coefficients Q_{ij} are the same in every equation.  So, we can pull them out into a column vector.  We have a matrix, V, whose entries are these b values. \n\n$ADC = V Q$                                                                   \n                                                                                                                                                              \nWe compute using only the diffusion-weighted data. Because the $\\vec{b}$ enter the equation on both sides of Q, we take the square root of the bvals and multiply it times the bvecs. \n\nThis can be computed scaled: "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\n% b = bvecs .* repmat(bvals.^0.5,1,3); ",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Or not-scaled. In this case we will scale when we compute the diffusion:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nb = bvecs;                 ",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Here is the big matrix                                                        "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nV = [b(:,1).^2, b(:,2).^2, b(:,3).^2, 2* b(:,1).*b(:,2), 2* b(:,1).*b(:,3), 2*b(:,2).*b(:,3)];",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Now, we divide the matrix V by the measured ADC values to obtain the $Q{ij}$ values in the parameter, tensor                                               "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\ntensor = mldivide(V,ADC);                                                                         ",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "We convert the format from a vector to a 3x3 Quadratic "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nQ = dt6VECtoMAT(tensor);                                                      ",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "To compare the observed and predicted, do this                                "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nADCest = zeros(size(ADC));                                                      \nfor ii=1:size(bvecs,1)                                                          \n    u = bvecs(ii,:);                                                            \n    ADCest(ii) = u(:)'*Q*u(:);\nend",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "We would like to have a measure of how well the model does compared to a repeated measure.  Or, we would like a bootstrap of the estimate.  That is done in dtiRawFitTensor"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\nmrvNewGraphWin;\nplot(ADCest,ADC,'.');\nxlabel('ADC (estimated)');\nylabel('ADC (observed)');",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAOm0lEQVR4nO3d4XKb\nShaFUZjK+79yzw8mGiLJCMkgzu5eq6amcnNzbSJjfTpNC8+ttQkA0vzn6gMAgE8IGACRBAyASAIG\nQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEE\nDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQKQ/Vx/AkeZ5Xn7RWrv2SAA4W28T\nWGuttXYrGQC96ipgBi+AcXS1hDitVhGf/iMAO9UfCXoL2PKIr7tV8Gswz7Oj2qPgIU0lj6rgIU2O\nareChzSFvPrvZwkx4uEG4Cj9TGDrvRsFX84AcKx+AjbpFsBIKq69Hqjm4jJAcRFPnv1cAwNgKAIG\nQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEE\nDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAC/M8zTP\nVx8EDwQMgEh/rj4AgOpau/oIeMYEBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAED\nIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAOea52merz6IHgkYUJFn\nfF4SMKCcpV7dNKy1qbWrD6JHAgYU5UmfbX+uPgCAe9LFHiYwACIJGACRBAyASAIGQCQBAyBS3i7E\n+e97Q9rDRqV59baRx38LQE/yAjb9jdM8z4+V0i2AQeQFbDtRyxC2/jPzv+/mVzhYm2dvumKaHp4q\nI+QFbPF0/JqeDWeKBT+53bHJdwl3T5URPcsL2OOMdaNVjEBvYBG5C/FpqCJeL8AvHXuX2+Ums3JI\nqLAJbKnU3UbEZcGwtbaxQRH60JoJDP4nLGDbK4e6xQic5rCIXEIEAAEDIJKAARBJwIDvmefDtlCC\ngAFfIl0cS8CAL1n2T9pFyVEEDPge9eJAAgZAJAEDIJKAARBJwACIJGAA7/FutiIEDIBIAgZAJAED\neI93sxUhYADvOfbnYvMxAQN4jwmsiLCfyAxQgYZVYAIDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAHH8FOy+DIBAyCSeyECx3B7QL7MBAZAJAEDIJKAARBJwACIJGBALbbjs5OAARDJ\nNnqgFtvx2ckEBkAkExjwidtlKgMTVzGBAZ/QLS5nAgM+pGFcywQGQCQBAyCSgAEQScAAiCRgAEQS\nMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAgL38rGRKETAAIrkbPbCX289TigkMgEgCBkAkS4gAX3Lb\nAmMx9hAmMOBLbGLkWCYw4Kvmedz5Y9i/+EkEDPgST98cyxIiAJG6msDmv+vrzSs9gN51FbDpb7rm\nedYwDmTzGBTU1RKiaMHCfj9GUGgCm//9hvu4Rnfj11EflpE5a+jeHPiSp0TAlgfuLi1Pf/ODj6NY\nDMhZz7uePgMXVyJgTxvzWXjkCmAQJQJ2iOX1go2IAIMoEbCns+q7BVIsgKGUCNjj3veI5VcALtTV\nNnoAxiFgAEQqsYS4aK3ZglHQ8jXxBQGqMYEBEKlQwG6bONajGJdrzfgFVFQoYACwX62AGbwA2KlQ\nwNZvArOJA8pabnXv1SaXK7QLUbcA2K9QwGyjhwi+OymiUMAmq4gA7FYrYCYwAHYqFDBTFwD7VdyF\nCNATmzZPUihg3gQGwH6FlhCnfxtmICvOTX5hp9ZMYKcoFDDFSjTPGgav+TY5Q6GATXYhRlleVPpC\nAVepdQ3M3eizqBdwoUIBAw5k5xvdqxUwgxccwncSIyh0Dey2cugdzfBLO7+BbCUlWq0JrLW2pMso\nBl/ju41QhQJm/IIva+1//4NEhQIGAPsJGACRagVsWT+0igjAS7V2Id79Aq7iJiNQX6GATW4lRQ3L\naahhUFyhgK1XDq0icjknIBRXKGBQhHRBhFqbOABgpxIT2O3SlxtwALBTiYCtL33d/Q5czg0DoaYS\nAVvYxAHAfoUCBjV5KQU12cQB73GhFoooNIHdfh7Y5BoYJUkXlFIoYJNuAbBbiSXEeZ4fN9A//U24\n0PL6yqssKKLEBPb0pzCbxijIWTkg76Moq0TAFooFwH6FAgZQkJfWZZW4BgYA7xIw+jHPdrrDQKoE\nzP18O/blrjiDYBDlroEtb2e2oYMPOGtgKCUmsLtirW/JQQdakxbgeCUCBgDvEjAAIpUI2N2aoWtg\nA7JmDLyryiYOt6Jnnl0qA95QYgKb7J4H4E0lJrClXrfB6+4fGYGvNvCuKhPY3Tb6C48EgAhVAsYX\nuNMS0BMBAyBSiWtgk00cX2Fplt/wcx2ppkTAXPQC4F0lArbm3WAA7FElYOslROmipsHX0Ib9i1NW\niYCt7x2182LY09tNqSBcaPDA830lAna7j9Se6mwXTrf4vZ+eiJ1cG2zD4vuqbKNvrd0y9jJRG5Wa\n59mGRvi+5ZtS4/mmEhPYzXoh8bNZavmvNtYkjWi85Bz5jMctWuJL/1oBu/lNvQ75UABDuXuqjOhZ\nlSXE34t4uAE4StEJ7C3LgqGfKNa9l5vc7IKDoaQG7Ond63WLM+gi1JQaMAb0MiEaA0Pp5xrYmFz4\nA4YlYMGWemnY2Voz20FFAgZAJNfAghkLgJGZwCLNs5VDYHQCBkAkS4iRLB4CmMAAiCRgbHGxDShL\nwPiRdAGVuQbGj1xpq+b2ksKXBiYTGAChTGAQw+AFayawbDZZ/J7HEEIJGEOTLr7MC6YDCVi2768p\ndfbttzyAluYgkWtg8QZ58j3vxyIP8gBShPPtQALGe3z7AUVYQgQgkgmM9+xZyjvj/bYmP+COCYxC\nOtshApzKBMZ79kxCpiXgCwSMQpQP2M8SIgCRBIzTubIFnEHA+BINA47lGlikxJ8LFXSoQAQB43TS\nBZxBwCJJAoBrYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARDJ+8D4UOLdQICemMCCuUkuMDITGB8y\neAHXErBgEgKMzBIiAJEEDIBIAgZAJAHjYPZGAt9hE0c2b8YChmUCO91o40hragp8gwnsXEu95vms\n53SpAIZlAjuXwACcxAR2Og0DOIMJDIBIAgZAJAEDIJKAscW7koGyBIzX+m6YSEMoAWNL91sopQty\nCRhbun9+XwrdfaehSwLGlhGe3/v+20HHBIwXPL8DNQkYAJEEDIBIAgZAJAEDIFKHAZu73/oNQGc/\nTkW6lgfAvkFgBF1NYK211t2T9/4bHQ2fb2AsXU1gT92NZd8s3O0zf+dztjbNs/EL+ETiClb/AUuf\nyd46/PC/K3CZu6fKiJ71H7ALfTknLoABQ+nqGhgA4+hwAktfM/zYqH9vYFAmMAAiCRgAkQQMgEgC\nVtr+dzEDjEbAAIgkYKfbP0I9zlutvd5baEQDxiRg51rqcl5jzv74AGUJ2LmW+WljilqmrttNNN59\nL9fLjw/Qqw7fyFzNnrp8ViCDFzAyAbuY4QngMwIWTPyAkbkGBkAkATvFZ29Adk0LYD8Bq+W2LV7M\nALYJGACRbOI4xWfbK27/ld0ZAC+ZwIqyhAiwTcAqcoMogJcErC4LiQAbXAOrSLoAXjKBARBJwACI\nJGAARBIwACIJWEVuJQXwkoABEMk2+opsowd4yQQGQCQBu54rXgAfELCLSRfAZwTsYsvlrruLXmYy\ngJcE7Hq2bAB8QMAAiGQbfUVmMoCXTGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQS\nMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAjaKeb76CAAOJWBDWOql\nYUBPBGwIrV19BABH+3P1AfAlGgZ0xgQGQCQBO9c8u/IEcAoBAyCSgJ2rtScXn+aSQ1nBoyp4SFPJ\noyp4SJOj2q3gIaXI28Rx+2K3hzKsz4PHfwtAT8ICNs/zrUzrX9/oFsAgeltCnOfZPA4wgidDTGXb\nE9jtd9a/+P5BAnSgfh3ClhC3WVEEGEc/S4iGLYChhC0hTs92IT4uGMb9pQB4V17AAGDqaQkRgKF0\ntYlj2vE25yIT59M3sRU5jMtXYn96cC580IIeqCLn+eUPzvZhXHt425/9kvN8zw0iLj+pHnUVsI1N\n9o877K9SZLPJ9mFc9VhtHNVVj9tPn3fdiToPVJ3zfLruwdl5GNce3k+f/cLnh6eHdO15/pIlxG9r\nrVU4CTYO48LD++moLvzO2f56LW+c//6xFTmLNhQ5vIIn+cZnv/Y83/i3V53nL3U1gW0rsrSSoub5\nWsr2Xc2uUuo8L/LIFFyUvvyzP7V9f76CBzxEwF7eQZG1Us+A09/juf1/nQOrptR5XuQs+ukwrj28\np5/92vO8yNfrXUMEjHeVOo/rPCnzliJfrJRVxMvP8yJfr7d0FbDW2t2GmeVUePx9Hi2P1fpl4HT1\nw1WzWAVPqmqHVOQsenoYt3+86vAej+ry8/ynB6rUSfVUxScIAHjJLkQAIgkYAJEEDIBIAgZAJAGD\nJ+7u6DP/6/H393yQl5/r3dsIvfzzd7vLoDNdbaOHQ9zeUbDeo7t+m87jHzhkJ/TdPu+jdjBX2KgN\nZxAweM/Sg5/ydrO+r+76z6znofUb79aZubt36uMg9fiOvenfxD49Zg2jMwIG/zj2if7pzbzvyrfz\nM9716ae37d/NhUf9RaAgAYN764tS+2O28YefDk/TmwuDGzUSKsYkYHBve4JZr+ztuQZ299E+u3K2\nMVeZtxiWgMH/7ZyiFntmqcc/s/6dPZ934wDufn/j40OXXNeFU5TaNFHqYOAo3gcGpzAJwdm8LgMg\nkgkMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIw\nACIJGACRBAyASAIGQCQBAyCSgAEQ6b+uwV0wES4HdQAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Here is the ellipsoid associated with the tensor.\n\nWe have a prediction ADC = u'Qu.  The diffusion ellipsoid associated with Q are the scaled unit vectors such that v'inv(Q)v = 1.  The vectors v = sqrt(ADC)*u will be close to the ellipsoid defined by Q."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\n% ellipsoidFromTensor(Q);\n% title('Diffusion distance ellipsoid')",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Here are the measured ADC values shown as vectors. These should be that peanut like shape"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\n%dwiPlot(dwi,'adc',ADC,Q)\n\n% pts = diag(ADC(~b0).^0.5)*bvecs(~b0,:);\n% plot3(pts(:,1),pts(:,2),pts(:,3),'b.');  % pts or pts/b?",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "A comparison of the estimated and observed ADC values"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab \nmrvNewGraphWin;\nsubplot(1,2,1)\nerr = 100*(ADCest - ADC) ./ ADC;\nhist(err)\nxlabel('Percent Error')\nylabel('N voxels')\n\nmean(err)\nstd(err)\n\nerr3 = diag(err)*bvecs;\nsubplot(1,2,2)\nplot3(err3(:,1),err3(:,2),err3(:,3),'.')\ngrid on\naxis equal",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "text": "\nans =\n\n    4.8759\n\n\nans =\n\n   23.9028\n\n"
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAWd0lEQVR4nO3d4ZKb\nOrYGUDN1HjiPkjdmfjhNCGAMQhLaYq26dSfpY4NM2vqQ2IhhHMcXAETzv7sbAAApBBgAIQkwAEIS\nYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQ\nBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACE9N/dDUg3DMP7D+M4\nzv86/QSAjgUOsNcsut5/kFsAzxF4CnEdV8MwzMdhAHQs9gjsNRt+vVYDste/84qQxsge2hQ4wN7h\ntEivtdZ6n3m+tkB79jkHgmYFnkJ8zcJJLwPwNFFHYO/EmhciLooSAehbW9M12bU2H0U4foWgWbGn\nEAF4LAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACE\nJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACOm/uxvQrWEYDr5yHMeiLQHokgAr\n6veB1/wq3gqAHplCBCAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAE\nGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAk\nAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIT0390N\nSDcMw/sP4ziu/wpA3wIH2Ovf6JpyaxgGGQbQvcABdjClppHZ8bfA/NcGaFPgAHt7j7d2uhuhRYL5\ngP7elgCfBA6wxcwhAI8SuwpRegE8VtQR2Hv4Na88VIUI8ChRA2ydUnIL4FFiTyEC8FgCDICQBBgA\nIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEG\nQEgCDICQBBgAIQkwAEISYACEJMAACEmAARDSf3c3gNcwDEdeNo5j6ZYABCLAWvD7wGt+FW8FQCim\nEAEISYABgR2cgadLAgwIzLXhJxNgAIQkwAAISYABge1cAxuGwRWyvgkwoCtTaI3j6ApZ3wQY0IqE\nAdM7ouZv3AwtQ7EuCTCgf4ZiXRJgQCtOxczXURfdE2BASOvJw4NMJ3ZDgAExbFYVJoy9DNe6IcCA\nVqzzaR5aJaoKjcZCsxo90Ip1PpUeLRmNhWYEBpSSNr459S5DqCcTYEApx8c3SgpJIMCAViTfyFyN\nAV9TBBhwj2ZHXTsp1VQ7EWBAKfvjlfolGwc10gy+EmBAJSWWhy89p2fOsGUCDDjqbG8+jmPpG7lc\nA3syAQYclbbsRVotYiNMJ7ZMgAGZJedQ42kxDEPjLXwaAQZkNvXyFUZUNZ/ILL1aI8CAo/bTYv3D\nW3p8T2R+DgEGJMp+I1f7T2Ru8Crdk1nM9xy/vgT1/s2d9+0JV3QWr29nfPO1Jbma2s5H5iXAkvw+\n9rJfZVsB1b1P4E5VFZ7q8c8+kTktTvLWYpw9JmQkwOARxvG1mD6I/ijIafLw3lY1dUyexjUw6NCH\nkooim91RoXOvHMP7R8AlhsoEGLD0qSMuHUj7T2Sus8crjMYqE2DwUA2uub65vG/RxlwsY+FeAgxi\nGH7s/OSURV9cYqDz5CcyN9uwnggwCOM9HHn3jO/ihflPEpy9kavcNbBmnw2WrI9P0TgBBjGU6BDb\n6WTX9yOfemM1CS0sfRnvyQQYRHKxanzemd6eFs1263mvDlrRqhwBBjG8s+diV1i5M/VE5rOaDfU2\nhQ+w+enkxWva0Li8l6nqd9OeyKyDyivwShzr34NGTt8gu6lw4/3XRe3G/Df/Sk3HkWacXRpqvtJS\nxMt461rNK3ustmbjQwQOsPVVX4uS0avN3+pPv+rlMqz0she3rwtV0/rDrhdcZl/gAFvbXButv/Jc\nKnjCPE+5tEjecuPf0LxHbDO9OCX8NbDJztnopHKTiOsJvzPlSgqnLd/7RObsyv0+DMPQ9e9aKZ0E\n2BPOl+F2DS7mu1YnNUvsZRzNH57TSYC9J/2z1BkDB2Wfn0++Ne1rS3LlzXShPW8/o9dKEz7Apn/4\nJ8z5QEbXb2Ru5+tWrbqvwkc2n3Rc+AADqjk7dVZ0yjHLCpDX5Z1ONId0SldViMBxBzvK0jdyJdus\nOr6rGdzCCAzYM5+Zb7CIo/StaQv7R+D6UEwcniLAgKVPHXHp7nW939JrL+XduPipTIDBQy367hZu\n+d9c3rdoY85uPO0JZ9ep7NgkwCCqi6tXN1tQV7ROpHF3jX2DEmAQ0uCJzG0rt5hWZ5l9hQCD52on\nKtb3I596YzUJLfQIlXIEGDzIvDO9PS2a7dZ3GpZW9GjZjkIEGHRuHlqVV6vZj6jNko2SzTnqxmac\nDfVmTwLqEGDQlfWE1afQaiQtLor1RObs+vhHTNZ0gJk7hk8+LWBddIyVUMRR+ltcugdfPbWrbI+U\n/HGe2Vu2GGDDz9PTH35yAfvqL2Bd+grQo7rg7PdQfz3O/R3eFgMMqKBcd3bLrWkVhiB5z6rrn6D3\nNyQQYPBQfTyReb7TytOJGV0pDd3cWn+DrU2NBtj7TMcsIjTlbH35Xd/fCkOxQrvIcsSe82TEtgJs\nfZPKQ84jIKLsX89cGyw9VAp0a9f8kPbXnbb1PLCHnDVACxJmOJp9InM1dRaQDH2lraa2RmBAT25/\nftgwvGKNOopeN+kvzBoNsIvLbANfHX8i8y3fRF//7Ee+v0PaYoBN5yBXltkGsrjricyZyhle02ay\n5MH+Fq5vv+g96UZgQP+6fCpV6Dr7LOYlcn2MDRoNsD4OLrRs8S2L8kSuCpe1dnYR+onMU+O7qbNv\nMcDm8xV9HGW4bn1h+OKl4golhSUX+yi04Zt1OfYtp60y+okMg7XpqY+L2/yTvyZn33j29Vm+vFOX\nPo5/b8N6/7xQ37De7LsNp3aX9tmPvEuvOGlxBDadVJpIhEmJPqv9fnDRB0xzX++fZ+khTvUzCZ1S\n9gtOV1asz9iMFrQYYO/f0fdZRvtfMKjp4tn3erGbU+5Yf3bv51mac+xDDa/X38Hf2e1H6cfCJVyL\nU4jTLIGRMo+16EqmW0oSvhHzN0b8Qn1u8vB6bf+3aXbx7DTjpz6n2nErOqn79ZXhfj1aDLB5qYwM\n45mSu9H1SfSnd9270u71vX/awnQApmnGtF2F6Hzab2FRLQbYa/YlfPg/D7wtrgpP0+zTX6dXToUe\nhZpRLniu2GzYOJ6+SLY+kpVdvL71tA6zxQDLUl4FPbl3XqvcHnN9wWc9xuLn6e0p3fnUX7S3v+60\nxSIOoIIro7T2n3tyxDAsSwQPjsCy3EzdwfnH7QQYPFR/3dnrz5BreE8e7sfMfDr21C5KzM5eKQ19\nshanED9N7gONqHYd6+QWXq9V29aTirM655S9JFxaO7DNcfEHjmgxwPqbqIWOZb+OldSGf/68s6V3\neH192b5A/VPfJQUtTiG+R2DdrJcMbWr5OlZa2w42ZPGyUxe0KnRKjSzbEUKLAfb6qRIePQ8MQrkl\nFNf7nI+ufh4JNl2VmP/wUp19Cf0NkopqNMCmEZh/Syikjy/XfKnf12x0tZghnJZz2qyzT77TOeVt\nZ7bvDH6fa2BAE9exNn29P3dRUvGzVP/8J39f+e+W//zwa5M/LUVfennl7Nvvr19tMcD6O8rQuPrB\nc9DBlXY3q/gWo5fFElMPNJ2mdLNsR6NTiEBoJZZg/zSldmRHi5ecmjk8NkorPp14fSPzmO8gvV4C\nDLrUXwHUT/n7nwHE/qfb7Jw3M6y1UkxPZD5FgEGfpttR7m5ITj8jp2VvPpXCrz/uvKxj/vqf/1pp\nmHhQ2gMHHquta2DJ8wPA2nTBI+1LVOJKyaeCiLy7mNfQb+7u/CK/x+cbi/dXV1as76w7bSvAVmvA\n9Ha4ob4x9WGYt3/75j3AagJw/Pev2wOs9U8+faadQ/SpGOTuw/N65e4kw3W5jU4hugkMMrq+LMAt\n01YHe4Cf7Hm/5e//vWZTi1+rNs7WNRx5bcJ0YtHj/PUDhuty2xqBvTqq74TWJA/F0t7yYTt//7x5\nkppw5round9cb/fz2GujYfPNblZ/nG1YIQ/vKtsKsM1nCjz8Xwgyauo2oM1FZo+36sjltLSPWHmt\n/fk/x5XrW1feHlRbAfa0ow+3uF7ckfermrC1TzNt4/h32nD9+iPl9as3Fl+jpNoaKP1dl2krwIBq\nbizuWOw3oRmbk4TTf5rt6PtP6rsrSDpLr5cAgye7q85+VUOYsvevI6f3Xc9HXry7l5w1fnk3iACD\np2uqzj6tJes15jeX9D2yKNR+26585ITrfOxrtIweqKnQ4/cSNphxmb6pdH6+jO+Vj5h0rW5vfxWW\nSplvv78lPIzAgD+u1NlvvrHaUOPTfqarbPPpxNKOF1Wqs79IgAF/JV8V26yJv2h/U8dXxJi1rUb5\nRt+Z0RRTiMDS2YWAFzdutrbE+2yDpabR+pudC8EIDNhwaihWYeZw3ZKaNynvFN8Hqi0M0chTjMCA\njz4NxY6Pz5odjeUqoMhYdTIpPUzs5jk7RmDAns2h2L3n8lkusyUs3ft1sFVh2Y4sm+1mKGYEBhyS\ndtreTV/5OjDYKndjXDmhR2MCDPjiPbBo55ksNacT43buB5WYAq1GgAFfLCYPk8/Zo1zHyvI0jAbT\nur8wDn8NbJp0nv5t4p5NQAjtPJOl8nNPTilU2VH0gIdbrj5wgC2WSClxHyXwycWFgLO7sCTxn/+9\n+ESu3V3UrjpJ20Ij/5THBQ6wxcALOOvTDcjHO7Irq081ZqP9Tx6NhRA4wA7ycGcSPOfEaHPq4lTH\nfdczWTabkbDfiL1CyEYX0H+A+Zcmwbwrv7clpWXMj6aeybK52XkLE1YYydKwEps9vvfO+kNViPBo\n14vjM26qdOXelZLxuMt2dFzg1skILG36Hh6uxJeltati3V/H6q8047jwI7D54iih78iDysrNjiYP\nxXJ9f0NM/BYdjYU4Atd1MgIDzlrPW+Sdxrg4FGt25FRClg97/EGa3RBg8FyLbi57r3elQDGhqnC/\n+25qYpMswk8hAo27svrU6/B6UUcuBdUvoDgoeZC6c5vQE2YRjcCA4koMxS7OTya/d+7eZTv2X/aE\nsaYAAyq5XqB49kaunZYkv/eTOhftEg5gx3OnAgyo5+KyHc90sTqj40PtGhhQ26mrYovLPJ1dx6q8\n5SMfNtAjLgUYcIPpXrHNvvLIEqYtB8/1ttWPkOFHoBtqTSECt6lZZ3/EvTef7UwVlqs6Cb2GkQAL\n4/gZWcRfRJ5sXdxx5Sawiy25vpG5xssdC225GgEWyO9jL/tVthWQ1bwb7a+441RJ4cHPfvEozYdc\n05xh8tbuJcCAO817z+Q6+2ZHTvsqJ8f6Q8WNrjcBBjTk4lCs/Sm769tPmF+9vtM2CTCgOe0MxV45\nQrFcad9O0necWxMBBlRyqktt56pYodq/clUnHVRnHCTAgBrmXfap5f6a6o6Pt3z/lUWnKHOVeOTa\nZjkCDGhaO0OxU2q2dsr4LAWK8+hKOOeoyUocQADXn8mSqxnZ93LxQTPzC2xj0oOwNzfbYFytGYEB\nMTS4bMdi42l7yVtVePGTzodcV7ZThwADIsnyTJagy3Y0dTmwBQIMqGE+u5VllHBXb15uv08uiE8j\nwIBK8va/yTOKCc24+ESuNPVza32S0Xh2CjAgsIszim1WK2SPjeljft3y4ucNHpw5AQbEVqi4I23U\ndTFNF7vLMl05bTbtVryWCTCgB9eHYov31ikp/LT3K1tbbDlLbX2bBBjQietDsbzjkq9ba/wKU/sE\nGNCVU0Oxiw/STFMtt947mv//zggwoDfHh1NFF8PdXMM3y5bPNqDXiUQBBvTp01Dsa4oUukP5+pY/\nVXkc3HLGW/EaIcCAbm3e8nyq7867PO71a2zzweXxqsL6N7HVIcCAnm0WkR+XqyA+eVPXG9MxAQZ0\naxFaRW9+ml1gKnuLVR+3cGXhcSpAt9brShR6Jsv8P5XMyDAPOqnDCAx4kLzLdtyyXGGdHYUgwIDH\nybJsR/Lb07azuJ3r/RDLzqoKzxJgQEiLWbuzXXnaM1nme7l4W9XZlQlLPL4yOgEGRLUZAKemBw/O\nKH6qKnRF6l4C7I/+7lGH7mVcMuPTpvZ3Ib3uJcDmfh94za/irQCOWTzjKsumQjzIkTcBBrRuEVHv\nUCkRLReLO6hMgAGt25zZK5QxeR+nstimGsK8BBgQz7r3bzMM9pcubKqpEQkwIKT1Kht3teSTBpvU\nGUtJARSk1L4cIzCAqzbLTNSDlCbAAK5yo9gtBBhAfpYurECAAeRn6cIKFHEAEJIAAyAkAQZASAIM\ngJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAITU1fPA5k/19twd\ngL51FWAvuQXwGL1NIQ7DMB+HAdCrPkdgwzBMQzHziiRwGgTt6yrANvNJaJFg8wQIaEo/U4g6GoBH\n6WcENo7jlGFGXQDd6yfAXnIL4En6mUIE4FEEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIM\ngJC6WomDt4PLQlq4BAhNgHXp94HX/CreCoCSTCECEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQk\nwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAh\nCTAAQhJgAIQkwAAISYABEJIAAyCk/+5uAJ0YhuH4i8dxLNcS4CEEGBn9PvayX2VbATyDKUQAQhJg\nAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJLFfLnBwaXrLVoP7BBg\n3OLIuvUWrQf2mEIEICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADIKSuAmz4cXdDInG4\ngKD6WUpqGIZp6bz5nwHoUlcjMACeo5+RyuYIzPwY13XzHYHO9DOFuEnXA9ArU4gAhNTPFOJrNmHY\n04cCYFNXAQbAc5hCBCCkDos45pWHi1rEe4ebjTTjtTpE9zZsUT46b8ktDdtsz6uBAwUsdBhgr3/7\nl0ZucG6kGZPNxlRu2OImh3VyVG7Y5k0XLRwoYFOfU4gWlPqqhUM0jmNTMbDZnhYOFLCp5xGY0+Qd\nbvQ+yIGCZoUPsEXP0tpJfZscooMcKGhZ+ABbdDFGXV85RAc5UNC4Dr+i61KxRorHGmnGuiWqEFtu\nD/BJhwEGwBP0WYUIQPcEGAAhCTAAQhJgAIQUvoz+UdbLPObd+PqehPlf1fsATRFgwVRejk9oAc0S\nYLF9uultvXT6zu1x8619TazF9jfXaF8sxQtQggALZp4W85yYZ8/7P+3cjTt/5fSTzfTanLTc2dH6\nNQCFCLBg9jPm01ver0lYkVYOAc0SYLHtzAdO1s/WqtAwgNIEWGA7zwg+cg1sbb8QcfONnlMM3MVa\niACE5EZmAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMg\nJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEj/B3ZPCFB+GtrhAAAAAElFTkSu\nQmCC\n"
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "##  A comparison of the estimated and observed signal values\n\n$d = S_0 exp^{(-b (u^tQu))}$\n        "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%matlab\ndSigEst = S0*exp(-bvals.* diag(bvecs*Q*bvecs'));\n\nmrvNewGraphWin;\nsubplot(1,2,1)\nplot(dSigEst,dSig,'.')\nxlabel('Estimated diffusion sig')\nylabel('Measured diffusion sig')\n\nsubplot(1,2,2),\np = 100*(dSigEst(:) - dSig(:)) ./ dSig(:);\nhist(p,20);\nxlabel('% Error')\n\n%% The error in 3-space\np3 = diag(p)*bvecs;\nplot3(p3(:,1),p3(:,2),p3(:,3),'.')\ngrid on\naxis equal",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAej0lEQVR4nO3d7ZKr\nqroGUD217/+WPT/c7Xb6FVQQXh2jVq3qzuwoSZRHEEg/DEMHANH8X+0CAMAVAgyAkAQYACEJMABC\nEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyA\nkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABC+k/tAqTq+376eRiG9SPT\nr+O/AvBuYQKs20qm6ZG+7zd/BuCtInUh9n0/b3VtPgLAR8Rrgc0bWOu+xAXxxn0a9NCmMAF20H94\n9ol1tdbDqTzHXANBs2J0Ia4rEdUKwMfFaIGtBxmuHzEKEeBT2uquya61/ijCcQhBs2J0IQLAggAD\nICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTA\nAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmC8St93fV+7EMAjBBgAIf2ndgEgp2GoXQLgKVpg\nAIQkwAAIKUwXYj+7NT8Mw/yRzV8BeLcwAdb9m0x930+/jtE1/1WGAbxepC7Evu97Q6QB6LouYgvs\nbANr3fcIP7lUgvaFCbDL2SO0uGDRQQ00KEYXokrk3SyfAVwQowU2DMNikOH6EaMQAT4lRoB1W7G0\neERuxeWjAy6I0YUIAAsCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkw\nAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgC\nDICQBBgAIQkwAEL6T+0CnNP3/TAM08/T48MwTL9OfwDAi0UKsHlijeZhtvkzAG8VJsDGWFpk2Pjr\ncVwtGmqFisfLrK+WgNaECbBNYyAd1zVCiws2e6qBpsQYxDFWIvP/d5IJ4NtiBNjwp0trdQHwelG7\nENfDDrXMAD4lWIDNw2kRVHIL4FNidCECwIIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYAB\nEJIAAyAkAQZASAKMYPq+81UEQCfAAAgq2Gr04FsHgJEWGAAhaYER1XQnTJsMvkkLDICQtMCISsOL\nEvq+9/XuUWiBARCSAAM+qt+aUaj5FYgAAz5KVkUnwAAISYAB/M9mvyJtEmAEYy1EcnEPLDoBBkBI\nAgz4KI2t6AQYwQyDKcz8I+9dK/fAAgkWYPNjq/+z+StAUWqb6iItJbVIr6n5Pz4+/1XPAHzH5fN9\ns65I35p6prowLTCxBOSlSokuUgvsmnm7zfFKIr1D0L4YATbWJtP/T+WQ0OKCRQc1LcvYPbO4H/HY\nfrkmRhfi8KcTSMC/7twDW/xwZ2s8L0YLbG0YhumYGw+4xa8AZ+ndiSVYgM2PmMXR42CC0B7ukZvv\nbnO/egjbF6MLEXi9+z2BWXanLycQAUYwFvPlmuwLHZzamtFAJQgwIlEJsHbQVJqH1jQQbP6vp7Z2\n+S8pJNg9MD5OjfFiJe45HW/wyQSSdiVogQGvUrSzzoKrTRFgQBPuDOLIOJHrOJ/m/ZCSrDoBBmRT\npU7PO3yxUF+ftCtBgBFMyihEIxW/oPryGWY9VyfAiGSqskRUm+5U02fbKPcjQahEJ8CIRIXzYj/j\nZG8ARQvfyKyHsArD6AlGhn3KfHn41hpMPxej2vtjctEC44WGQc7Vkb0hsp59vPk31zZeaCLz/TdB\ney6RAAOyyXIP7LHqu1CTaHOzml8lCDBgqW4LQF3vHUgkwID67twiqjWIQ0dfdQIMWHomS3JN5DKR\n+bOaGIW4/mi1oJkOCsfCW0U/zcuV35DFRK20wIY/48+uVhg5i2M5rnaPV8J9/qxfl7ZcGSzbUUIr\nAQYLKUPhrcex52ZFXKgePx4TX/Ee2HwLD+To5V2MB7xjfiLAgKXPtgBOvXDfyFxdEwE29hmO2pxy\nz/OOrzTHf1q30lyf1jVW0w9X1teqi8WXNa+39vNVXK6mUpbt2HmiSfr/qBxg68PddQo/jceII+Wn\nWneVLidK7uIc+bnGR97LaMt2lFB5FKL2FnsODoph2E0vh9LcnbtKZ5/b/sC5+bKK7bi/bEdrr+hJ\nTXQhTnxdN8fmY+s/fNq2Iu83cpUYxLHoJ0zZhYnMgTQRYFNH4pcvJfhJdfGYxDPxyWHop4ox/6eD\nfzWRObomAgwoIXulGaKPpGgJ94Z+/GR8YwkCjDB0G1ZXdODDnYr4WmfmhdIm9kPe9+VYStfEUlLd\nX//hQS/i+qpnccheuyyCF3tyEMd97Szmm1iScu+SQRyJmgiwxG81Hf9181tQ5w+6lwZfcOoLkau7\nXC+p0A6E6ULc/AhD9MhDOI00ho7dT69C38h8X97hnS/WRAss0d6i9cfnjEOBC9q8MFr0k//sNj91\n8Z7lSv/mPbCfawHf3MVCxk39LHwLS9e/rzEXKcAST9rNZ8Ep6esJPWx+FvzsNk9ZsmizHdNOTTd/\npY0UadNjy3as/6nlt6W0JroQ+39t/sHPR+D1sldVTdV9exl8eeXAa+5PZG7kO1myP71BrbTAfra+\nFw2v9SNGIb7GtFAvm6Yhu9ee2yWcI3WHL57dSN7SFupHzduobXNZrOe1EmA//TzOfJa83tlqa11p\nhjhNQhSyKzyROf2OWjv9vc9rogux+7cXsXZZoF0P96eddacx9HzJm633LWeV6OXR/eVrE7Jo5xBa\n1D4Hk/fv9y8986pz7SV771wjn3iiL3cnttKF6A4Wubz1Flr63ZrLd8h+bjmv6Cf787MUNkV/G+9o\nogtx/GhH72vkkotvW27fz0mZJU7wjBOZT8XMMxOZn3liUE0EGMxtBlX6iWnN30lrd5WOx8RXLO0D\nu14Mu8i45fS1Gt4XbwKMthx/1bJkOiXEilBNKdGoutBzkGV84xc+xCYCbOw5HH25P5fur/20eRQ4\nNKK4VnU+HLd7DZfp8Z+bvZwQtWq599WulQOs//su5sUjQEU3s+ThivLa7vbuuG+urXXHeEEWbtmO\nECoH2LSsxlzdIsFrPHk52P6l58EQkrGXr0pH396yHelbTvnj9j+da5roQpyYyAyBrHvhbk4+y1Cm\n1TbnU3SirEtSaNHL99WuTQTY1JHY7PEEEZW+q5T3hC1R2sROnfmd13C1fAvj/mtpIsCAphwMbGu/\nfi9awvk7kL2jj7MEGJDqZ4OmSjW9GI5xanWMza39+zf/3BubvwMnGzTnFvMVeClaCbDpGyLe18iF\nWlqbyHyswdIelyi9wBde2c0YTvnLF2RkEwE2XdRIL8jo3ROZr7W6TjmePl+uvjq1bMflYrygvm1i\nMd/1Mtu1SgLUcvbEny8EVW41+oMNL/5yqsbWT7lTuowv7Xgdr4gVb0MtMJPA4JS8X6GXpUHz2ACK\nubylvba1cst2PCNo3dtEgPEFbZ+/8aR8h8OpWrXNKux4Ildd/zbUhs1V0KYxIKdKfuoDbTwai2oi\nwPqZ2mWhiPGD9fG2ZhiGvF9Sk30QR9HQKrTl+/XYqWU70iUWLFA93MQ9sKauqj6r6PdADkPX91bj\nfdr6zEpcq/DUHZGDez+nXL6nVfQe2AWJt9DumH+O2ScyB6qQm1vMlxfzOVe3tZzS5qCDmlVYUxXo\nQQv1VOP1eBnGs1rrTa2lfhfilGF6Eesq+j2QPtXH1JhNlWEdpjvFvj+IY3Gn7XJJzhdg73EnTJIm\nhk6WG8EZdGzoKwXtQmz5ENqrc6f+pZYL320MQ3+0tIm9qeWLsXFenOzCTX0hjR8PFzTRhcgXvOvE\nacLe5JOrax39zzMnZvWOyuHBNdr39rL5Hpx6Z77cnagLEXjI8Qn+QNzO+gznY09+Lyi1uYdTldUz\nGXNcpPflXOVRiNMV0PveWYjrTpYcPLfYgk95huGN19IpeTZtL71jfP7OlKvxvlaRttKFqAUG2Tmb\nut3pZRszjrudANgZqPm/J6aMJMwykfnUH3/h02+iBZZ4+3H+lPUjVQYRAXN3+lRytUsWgxpKzoO+\nOLz2Z5EMaktU+fWkr0U2vfX96otXFsfrM011PiLoIXSz2LFeda3SzmdwX14HoG7h43zI2yp3IU7j\nZ/bGU83/8sFyEVXehZE+61r7qURJ0vf+WAFWX3F55Vmzp9cf3BFXE0tJJbr2GcyfJQVJ9NYTPq95\n66HigKzSo+H3Jlpdam8lPevO27jZtHpr1RcpwK4dpm/95F6vbhfH+lZrRKXj5GAK2gU3Szv7yMZf\nL29pd8t/v3brHSXursSnkf6+/XvBkb8kz6s/D6z7WzJgtHcD7PlSEVHRBbHCyT6zqvo44Vp7n4+X\n/vvh3BZOHZlnJzKfK8qLtNICO/4Mhtk3Hk33zIxCbNM7bg6z6eH23NUCjFVBzqLO9jtMR/jYHzj/\ntct05N9Z5urgSe+rG1sJsJ9+dla877P5OJ9nRYUmMj+zzVwF+Hc882IX//v/5h9sbm3M1JTSPTay\nI3q12UQXIm+iE68dVdaGzyLX2KuzfZ4p+x0f7vvu1LYvT2T+uY/LH1b09OoEGLB2uWq7GTZZtjN3\nfqHbYR0H60Hz0/2K1e5SF+b4UYhZmTPGzEHUVb9kuaZyF2L6RGbg7H3fs1f9FVcAeKAC3XxRP+82\nzWcr70lckqPQO5plKcigtW4TS0mtl9UANs2j637erCdyrX/OrqmJSgnL+45/tnjW9l2x9V/ONpXh\nmuAFd63y0oUIYWSvvEoMuMi7wYzSRjlur9s7dQyuX1+hqV07xbNsxz8EGO2Kcx496tTE1fUjRaun\nnwXr+z5xmE/eoi42lbjlg6IuHp+ybVHsVa9vnqN6vpf0dyn5VYdp5DUxjH49qQv+qoOPjmlc1DWb\nPYc/XV4po1xXVZYbNvf3e3bj89U3Fsfk/Nf153N1LtePZz3W39s4LTBKuXmx+eGzsuv+hs+tF7lO\naeIULtpp94t0fx3UEsb22dacsNNto2e8L+qauCU4/4aUvJd+7nlWlDJ8q33tHELrNtlev8WddRxu\nmt6uWidylc+r//sq5yrL0FT8uKtroguRV/rkCVVQen/glG0V07doB+Bj+92zNZRjmP/rzVLUWnwk\nnFa6EFtrawOnzBfjLjpUZD0cY3/p4dNbvrF2Sdd1i5thy42nby37KswXyhBCc90jecvTTv8PQQU9\nhG4WO+KrXo+tGGV8HdN6hpNFYk3rdER786JqpQW2vlkN1PLM8kXZN75e8+ny6/hZ7PmaUnuTmi9v\nPIv3tbfWmgiwL7zR8A61ztYLy87dS+Fhmte13t3mlg92txiRG3EoaZtaGcRRYh1P+LgSPYHnZ1Bl\nW5b31x9k2c+P/f47x6vsrhPtfcqbCz++rHZtogW2N+UFuCPXWIBcDYKz2zm1xkTeZTumCV7j16bc\n3ODeih7Hy3Ykb/z3H7+1SddEC2xz0QGgor11fu9IrGov7PdUCdeztRZTqWarJP/3gfSN31Gu6hu3\n/L6qtYkA+3cKxTuvFCCQEpOR5w7m3j62ZPt8rOD+jLrlAI2z78nxYMhF1fe+gCmtiS5EoIRTl4Pr\nXriiDYJZKyfzpLHku2VJ+130/mXJ1wsDUo43eK9EgTXRAvvyB8CBdyxGFcU71uw4NrVyxv89sMe/\nG2n/m0B2th81o/c18poIsJe9p9CIEMuKFl3KL2ixT7l8p/AFedbWC8h+TLzgE4qiyjKmDwh6CB0U\nu9Yrsl+ya+IeWPWFR+E7stzFKbTfEncT5jfbsm982vJ647lqM8t2HGgiwHiBvZkuVDTVSg+PklgX\nIPuWE/c7XwLjgTEpRQN4z6nZcpf30iYBBu9U6Ku5UvY7/fzYfhdtoIr7fazBVGgB9FiaGMTR/Z1j\nehHhpsXKDkVPqM2b1hlbUWc2NUb1shtgvYG8N2tvjH2/VYxr+018VqB6OExBrwn0ScTy1iEba9EP\noQfKX2s83mK/KcdklokZWV7vhTOonXGP7Wji5EyZPrluL++13E1uf8DB6bf5tUxxP4QXHEJTxRd0\nwHreja+Oz6T35FoB7hd7bwvCbNRKF+LUf/jzNuPmUmnzB19Q47Rv7w2OOZTp5Z65E/bMxu9X3Ks+\nxpDvybUtv69uDDOIY+8yJOjoT7hmfsz3f1KeOF4d3jlfsox2u7nf+RpU2S2KXWVUyPF+F98rtv77\nYy9Lr66dFliixRVEytDVWkfhK/3sD1z/U9C3vM0Lo/nAwm6nB+LAzabYF9aaCrHf6N3yGTUUYMf9\nh+uug8TDS2jl0mSVXsoDs19rOXUz7FraZZkr/fyZW3EgdPIlyMZp+OUqrokuxHnnxsGHsRid8UTJ\nmBnf/rMny2anB5eNZ8rNOmvsiNub5py306LxZTv2JpCVXrbj8sS1ccWA/ZvQqd1R79BEC2w+CWzz\n5Oz/1poaf12cfouOxGavR17Q8A9d+OgWI5Xub3CvKZb3DGptIEPiRhovdpWNN6iJAPtp81MpetZ9\nUKF8vbzBF+R9+9Z31J7sQ8s4FvxUsRuZuNaVfLcXn2lXePxLLa0E2PvatmuvO3h4VIleh/kUsedv\n/zw8tG89/ebF+53v7n25NWkiwKYz833TFAJp7Y1vrTwtyN7rsKjjHriOLNT6CVpvBC12O5oYxNH9\nNW99nPcZNMFlx4M7LltP5Co6g3Oz/Lk6KlOKfW2/99+Qvf2+eL5s5RbP5tuasUgfbNK5dZTXBw+h\nzkpFCZxoLWiuBeacuanZ90/T8GGLsdrrodsH1k2xs5fw6X+ft3Ew39qTLY+MOzpV7HKfSwitXF0W\nGgT/zcvnNgW9Yg16CG0OQpv+qf0Be+nqDSnsuhrHc/ufyJOaGMTRzc60oPUFNGg+yPDaFp4fW39W\nrYLVej+a/SCqaCXAXFa8ns/2YRlnPc+nPB+07e5LTMoL+81S1L2NlJ7ONcq+i4Nit3zJMlc5wNpf\nPoO6+tV37FLFZlOs1mj4C/u1bEfixmO1JVoZxNHP1C4LrRiPBUdEUxYDQ4rtZTnqJ+8oiXJ/XGLu\nwYWSXNjddAcnSnp17QziKCRKQ5g91Yd+xD2EFt0b93s7+n8XYyz9tlT86Gt96M/vd35URDzU45X4\nlIgfCU1xCO2J1dfEwubHF+5ob6ULERbMG2vcA8t2TI+U7Kss1R2aWOxr+71c2uOuwljp1WmB0azp\nDK37ATqEftIUa99bh8u1MoyeKqrfYTrQZqlebHGLa/whpb5bDKkvmmelB6yHK/bPjb/78kKAAcsV\nmC4M1jg15flyrVpoItf9LafsutAuDkbD7/1rusuXNc8QYByZN9Fabq5xR//3Zej3N5XYAsu1NvzZ\n7dRbtuOfteGLFiNjk+v+ZU1pBnF82jAIJDJbDO6Y1JriWWt26cM7nc+jrXV98DwtMI7MD+DmD2au\n6P++S7bLWmetm2JPLtsxfyFP1sIp+82+Xvl6syGyJwstMPi04U+Xe1XD7t91hIsOhV9sPPsLSXR2\nv3fek/lo+MV+77/8xWXNza2VowUG/GOeOperwnlDpPRAuFPjR65t/MCd/V544jMjKRq83bVJgMVQ\ndACF0Rl0/1aIeeusdcCUqBYfqGc3i/1Y/V53QHyWy5rsBBjwhAeaYq9UNzbKXdZkIcBiGIb86ypN\nDa/2Dkvata5P02vY+319tSLw+f1mnMi12E6DDanLBFgY8Q82XmKeQxc6Boedr8dM33WivIMqy218\nveXNPV528/NqmQALJuP9quCHLnVkqfKiLNtxsIXsVf+iYZQrXaJH1DEBBpyWpXo9aIFln8iVHpbP\nTyDbfBOyTxd7ZZIJsGDedxAaAxlL3g6uvdZGC2s+lfbMvah3D5wJE2A/b0W+6c4ktKzECPhn5opt\nerh18nxN9eIqMUyAdYe3Irt/g+3FH9j7+KwCWSzNMK4BcbNGnp+/88Ed5c7ixcbzLttxnJR5d5dS\nni7359WUMAF2c9ztzY3wNS0vn1PR5hl087RaZ0npa9CH23lZAuPaRkp8Xk0JE2CjQF+g8Hrvvne1\n7qnmMYsexSx5ttn6KdokytjQednY94zCBNi7b0UCc/PoKjeavHR0ldsFozAB1jkUGlNicRDe7Wyj\n5ObgjhZGZ2gwFRUjwFJuRb7pzmQU3mnSXesHW/Qi/nzi5YlchRa4Uh0VFSPAUo4MBwq8VXpT7HI9\ncOGJLpqrixFgtO/dYzqobt0UWz+YS8po+O7B3HrZ2PeMBBi7pltcThkasdkUKzGxevPx7EPJ0peK\nl1ubBBi7Tg3TcH7xgKl+f3JwRNHWz8H6DELrJwHGEWcQuWTpB1sM8Js22/5crk0i6iYBVo2bRnxN\n3vq63LIdmz17jy1tRToBVocZVJBFxmU7DtpbWQJmsarLzVludAKslvH2Ut3jVhOQplye2Xl/2Y5n\ngsREsewEWDWOW1hYD5HvkhtV6wZNyqzn+X4f9vql4h8gwD5KHyatuV9lb+bfphb67jTI7hNgH+U0\noU3TgPLLWzi4GVa0fbP+ksJCO2IiwD7NbTCqWA9n6HK3iuY9ig/EyfwVmc71GAEGPG2vTi86zr7c\nV7vdbzVyjQD7NJeGNGI9oiFXs2m+OL320MsIMKC+QlOvsm9zby7X9H8B+SQBhjthkOogaKXX8wTY\n1+m3h7xM53qMAPu6FtYEgRbc/JL3ot/wwiYBhvSC5TdkdvoGI/i/2gUAaEXf9+IqEC2wwAy+gIzW\n62jQOAEW1eIUcx8L0m0uBUI4Aiyq+RnnehFOkVjvIMAANsa+GwrfPgH2Bs2eXzo2CWQRVHKrfUYh\nUsp4/ap7EyhEgD2n7//730eM16+uYoFCBNhzPliVf/AlA48Jdg9s8ztPuziLjzVcNIBgwgTY5tTC\nzbVeTKQH+IIwXYjDMKxjaVz3pUp5qvjULTSAY2FaYJtSvil80dNYvEy8wqcujCCowAF24TsOonvR\nS2mdZfGgfWG6EBdUK+H4xIC8orbArPsSyzSp2YcD5BIswA6+81Rutc9HBGQULMAISnQB2UW9BwbA\nxwkwAEISYG9ggjPwQe6BfdeUeW5QAREJsDeQQMAHCbDWlZs7JfaA0NwDe86FO1W+1BhgjwALQFMJ\nYE0XYtNEF8AeAfacLGk0dicKNgBdiJG4GQYwEWCRjA0vzS+AToCFI70ARgIMgJAEGAAhCTAAQhJg\nAIRkHlgYFo8HmNMCAyAkLbAWbTa2NLwA5rTAAAhJC6xFGlsAP2mBhWEhRIA5ARaD9AJYEGAAhOQe\nWAzuigEsBAuwvu+Hv7q8/+tWGx9Z/ArAu4UJsP7fu0DrJJv/KsMAXi/MPbBhGF4TS31vUAbAXWFa\nYJfNm25NRWDfu7PVrt4lBjTv/QHWVGgRxfpWK9CaMF2IADAXtQU2DEPcUYjNFxAggGABNg+nRVC1\nn1sXjKH8xlcGcJcuRABCCtYC+xptL4A9WmCBbY6PM8kM+AgBFtWYUrIK+CwBFtVe7+Iw6HgEPsE9\nsMAEFfBlWmAAhCTAHpV4y8pADICfBNhzDLsAyEiAPSf9lpWBGAA/GcTxKLEEkIsWGAAhCTAAQhJg\nAIQkwAAISYDVYaYXwE0CrCYZBnCZAKtjHE9vVD3AZQKsGukFcIcAAyAkAQZASAIMgJAEGAAhCTAA\nQhJgAIQkwAAISYABEFLgAOtn5r/WLtcPrZVQeYCgYn8j8/C3mkXf95s/A/BWgVtg3V+rq3YpAKjg\nDS2w4wxrMOFaK5LyABEFDrCUfkJ9iQBvFbUL0UU6wMcFHu8wZdiiIzHuKwIgXeAAA+DLonYhAvBx\ngQdx/NRIp+I0L21dnodL+LPTtW55Ukr4TKmaKg+w57UB1sLU5vlIk3V5qpRwXhE3VZ51AdYlfKA8\nBx9ZlfIAB3QhFjQMQ1PVXFOF6dorj1iCWATY5zRVTVtLBbhMgH3ImBbtpFfXUiN1WhK6+7cjEWiW\nAPuWRtKiay8khj9dS+8ScKCt6/G8Ghkz1sgoxEVgjIMmKpZnc3ctjPozChGieHOAAfBiuhABCEmA\nARCSAAMgJAEGQEgC7FH9v/b+ZvFD+sav/cFij2cLcHNAfGvj6YEoXrsWYrPSh32uh26XHsZ9bbM3\nC2McLHCNAKts3v6YImqxuO18xdtua7WI+RPnD653sbfrxSPzH+YlWU8mWxesW0XvXmH2nt5JNSCB\nAHvaupqeV9aLRdkPrBdKX+fH5hLv0897oTJPlM2SbD64uZj9QewdP/3wpQN0nQB73qL6vtwxeFDL\nFw2Amz2ZD3SEAh8hwCq79i1cB+2V0u2Ym18b1sL3tAHvIMCetuivm7dI5n/zs2b/eWdr6sTb/Mu9\nx3/uLmWDB4XPVR4Al8A0R8sMSKEFRivcGwNOcakLQEhW4gAgJAEGQEgCDICQBBgAIQkwAEISYACE\nJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgA\nIf0/9i12P5KulHIAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    }
   ],
   "metadata": {}
  }
 ]
}