example_qiba_v6.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:a2dabd85c3f66109a8b9bc0216b9a789fc76be60331b689ffb24236773581e86"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 0. Parameters and Initialization ##"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.io import loadmat, savemat\n",
      "from util import mosaic, qimshow, status_check\n",
      "import matplotlib.pyplot as plt\n",
      "from numpy.ma import masked_array\n",
      "from skimage.filter import threshold_otsu\n",
      "import time\n",
      "import dcemri\n",
      "import dicom\n",
      "\n",
      "\n",
      "# SCRIPT FLAGS\n",
      "plotting = True\n",
      "Rel = 4.5   # Relaxivity of Gd-DTPA at 3 T [s^-1 [mmol Gd-DTPA]^{-1}]   \n",
      "flip_angle = 30 * pi / 180.0 # rad\n",
      "TR = 5e-3   # sec\n",
      "nx = 80\n",
      "ny = 50\n",
      "nt = 1321\n",
      "\n",
      "data_dir = 'qiba/v6'\n",
      "file_ext = 'QIBA_v06_Tofts_beta1'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# add cubehelix3 color map\n",
      "mpl.cm.register_cmap(name='cubehelix3', data=mpl._cm.cubehelix(gamma=1.0, s=0.4, r=-0.5, h=1.5))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 1. Load data ##"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data_dicom = zeros((nt, nx, ny))\n",
      "T1map = ones((nx, ny)) # s\n",
      "\n",
      "t = 0.5*arange(nt)  # ms\n",
      "\n",
      "for k in range(nt):\n",
      "    file_name = '%s/%s_%04d.dcm' % (data_dir, file_ext, k+1)\n",
      "    dcm = dicom.read_file(file_name)\n",
      "    data_dicom[k,:,:] = dcm.pixel_array.astype('float')\n",
      "    #print k, dcm.pixel_array.shape, file_name\n",
      "    \n",
      "data_dce = data_dicom[:,10:70,:]\n",
      "nt, nx, ny = data_dce.shape\n",
      "S0map = ones((nx, ny)) * 50000.0 #\n",
      "data_aif = mean(mean(data_dicom[:,70:,:], axis=2), axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# subsample data to speed up the run\n",
      "deltat = 2\n",
      "data_dce = data_dce[::deltat,:,:]\n",
      "data_aif = data_aif[::deltat]   # TODO: do this better\n",
      "t = t[::deltat]\n",
      "nt = len(t)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(data_dce[:,:,:].max(axis=0), interpolation='nearest')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<matplotlib.image.AxesImage at 0x10f67b550>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAANcAAAD+CAYAAACp3pBLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACo5JREFUeJzt3W1sXmUdx/FvN0YAESfBdMu2rBhhEITwoKisuEKYguL0\nFWACmYT4wmAgEhHGC2ExuqkvwBfGGAUlQBASAtlAZDzdmEFEiJuMQXmIdA6EYQJMjDGiu31xndve\nbU97Ttf+e53W7yc5OQ+97t7/dPv1us7pnfxBkiRJkiRJkiQ12NnAIPAScFXmWqQ5Yz7wMtAHLAC2\nA8fmLEhqknlTeO2ppHANAe8BvwK+OA01SXPCAVN47RJgd9f5q8AnRg5Z3oZdU3gLaTY4DXiiZ/TV\nqYSrXT1kF7CKNLn1dW05LSn2m4A1OQsZZSlwK3Bh7kK6LAF+AnwtdyGj3EreW/ytwONd5z8sHTWV\ncL0GLOs6X0aavUYZAFrFXpoL+outozxcU7nneho4ijQVHQicT5oOJDG1mevfwNeBB0hPDm8Eni8f\n2jeFt4myIncBJU7IXUCJj+UuoMTK3AXUMuYmbJq14drgt5isJdVDsliau4ASTf1ZNa2uI6AkS1NZ\nFkqagOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIY\nLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoLUCddNwB5gR9e1w4EHgReB\nLcDC6S9Nmt3qhOsXwNmjrl1NCtfRwMPFuaQudVsI9QGbgeOL80FSP9Y9wCJS68hjSl7Xhh9NrcJp\n18QeWNC8tjiQehs2UH/1kBm1tQemsYVQLylYFPve/fw+0pw1HQ802tRqPi79f9nftq2d5eAbwGLg\nzfGH3t91/BEau9SQ6nqnBXtblcP2N1ybgLXA94v9PeMPPWc/30JqqIUDaevYvb50WJ1l4e3AE6QO\n3buBi4GNwGrSo/gzi3NJXerMXF8e5/pZ01mINNf4CQ0piOGSghguKYjhkoIYLimI4ZKCGC4piOGS\nghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKC\nGC4piOGSghguKYjhkoIYLilInXAtAx4FdgLPApcV1w8nNR1/EdgCLIwoUJqt6oTrPeAbwHHAJ4FL\ngWOBq0nhOhp4uDiXVKjTn+uNYgP4O/A8qfX8GmBVcf1moEVpwJrWpnVl7gLK9S3IXcFYA7kLGMc3\ncxcwykfLL0/2nqsPOAl4Eugl9Uam2PfuX2XS3DSZnsiHAncBlwPvjvpau9hK3Np1fEKxSbPY71vw\nVKtyWE/Nb7cAuBe4H7ihuDZIWji8ASwmPfQ4ZtTr2vDrmm8xUxrabdZlYX2NWxb2QEmW6iwLe4Ab\ngecYDhbAJmBtcbwWuGdqFUpzS51l4UrgQuAZYFtxbR2wEbgTuAQYAs4LqE+ateqEayvjz3ANXWNJ\n+fkJDSmI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGS\nghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSglSF6yBSF8ntpBZCG4rr\nNhuXKlSF65/AGcCJpJaQZwD92GxcqlSnhdA/iv2BwHzgbWo3GwdYPpX6pl9/Azs4QvqV1TRfyV1A\nuXUrvp27hBE2jHO9zj3XPNKycA+pNetObDYuVaozc+0jLQs/ADxAWhp2m6DZOMCPu44/Dpw6qQKl\nptnVGuLPrV2V4+qEq2MvcB9wCmm2WsRws/E3x3/ZpZN4C6n5lg/0sXyg73/nW9f/tnRc1bLwCIaf\nBB4MrCb1RbbZuFShauZaTHpgMa/YbiE9HdyGzcalCVWFawdwcsn1t7DZuDQhP6EhBTFcUhDDJQUx\nXFIQwyUFMVxSEMMlBTFcUhDDJQUxXFIQwyUFMVxSEMMlBTFcUhDDJQUxXFIQwyUFMVxSEMMlBTFc\nUhDDJQUxXFIQwyUFMVxSEMMlBTFcUhDDJQUxXFKQuuGaT+pssrk4t+G4VKFuuC4HnmO4g6QNx6UK\nPTXGLAV+CXwXuAL4AjBIajje6TDZAo4peW07NaRskHsPy11BqQ9/fmfuEsY4nztyl1Dqez//Tu4S\nRuj5atqNvl5n5roeuJLUG7nDhuNSharmd+eS+h1vAwbGGVPRcHxD13E/cHrt4qQmar2QtipV4ToN\nWAN8DjgIOIzUunUSDcfX1ShXmj0GVqStY/3m8nFVy8JrgGXAkcAFwCPARdhwXKo02b9zdZZ/G4HV\npEfxZxbnkrpULQu7PVZsYMNxqZKf0JCCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjh\nkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGS\nghguKUjdRgxDwN+A/wDvAaeSmo7fASwvvn4e8M60VyjNUnVnrjaps+RJpGCBTcelCU1mWTi6ofIa\n4Obi+GbgS9NSkTRHjOlAPo4/AXtJy8KfAj8D3gY+2PV93uo672iP7FOe38ntx3OXUGo91+YuYYxz\ntzySu4RS1302dwUjrU+7MVmqe8+1Engd+BBpKTg46usTNB2/rut4gPH7lkuzw1CxVakbrteL/V+B\nu0n3XTWbjl9X8y2k2aGv2DoeKx9W657rEOD9xfH7gM8AO7DpuDShOjNXL2m26oy/DdgCPA3cCVzC\n8KN4SYU64XoFOLHkuk3HpQn4CQ0piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI\n4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKYjhkoIYLimI4ZKCGC4piOGSghguKcgM\nhas1M28zCe+2tuUuYYwdrbdzlzBG64+5KxhrKHcBNRmuBmlkuJ7JXcFYQ7kLqMlloRTEcElB6jYc\n318tYFXwe0i5PYbNviVJkiRpTjkbGAReAq7KVMNNwB5gR9e1w4EHgReBLcDCGa5pGfAosBN4Fris\nIXUdBDwJbAeeAzY0pC6A+cA2YHODaspmPvAy0AcsIP2DHZuhjtOBkxgZrh8A3yqOrwI2znBNi4AT\ni+NDgRdIP5vcdQEcUuwPAH4H9DekriuA24BNxXkTasrmU8Bvus6vLrYc+hgZrkGgtzheVJzndA9w\nFs2q6xDgKeA48te1FHgIOIPhmSt3TZUi/4i8BNjddf5qca0JeklLRYp97wRjo/WRZtYnaUZd80ir\njD0ML11z13U9cCWwr+ta7poqRYarHfi9p1ObfLUeCtwFXA68O+prueraR1qyLgU+TZotus10XecC\nb5Lut8b70EPOf8NxRYbrNdKNe8cy0uzVBHtISwmAxaR/vJm2gBSsW0jLwqbU1bEXuA84hbx1nQas\nAV4BbgfOJP3MmvSzKhUZrqeBo0jLngOB8xm+Gc1tE7C2OF7L8H/umdID3Eh6IndDg+o6guGnbgcD\nq0kzRs66riH9Yj4SuAB4BLgoc02NcA7pSdjLwLpMNdwO/AX4F+ke8GLSY9yHyPcYt5+0/NpO+s+7\njfRni9x1HQ/8oajrGdJ9Dg2oq2MVw7+gm1KTJEmSJEmSJEmSNKf8FxxBtEfMBrJEAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1037cf250>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(data_dce[:,-1,-1],'k.-')\n",
      "title('example tissue signal curve')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "<matplotlib.text.Text at 0x10f693f10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEKCAYAAAD0Luk/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8lPWB7/FPkgmCEHJZiCBiQlWiWEAkK65iGespxZwe\nRS2VSot4KfVlV311z7GVnvaY1nNsbbcXu2d1tbZyOUAP6i5HC2GFrklvr0pludVUSNwkSoQkSBLC\nxRBIzh+/32SeTGYyk2QmM8883/frNa955jfPzPwml+c7v9szICIiIiIiIiIiIiIiIiIiIiIiIiIi\nIhJXlcB9cXquZcC/xum5EqkDKB6B16kHbhqB1xEXy0x2BcTzeuxlsIqBbvr+Da8HPh2HOiVaDuYA\nnWhD/dmKhygExO0ykl0Bj8tKdgVkeBQC6e1C4BWgGfgP4CFbXgC8D3zG3h4H1AJfsLf/M7AbaAfe\nAx53PGcx5hP4Cnvfh8ADwF8D+4BW4B8c+68Afm/L2oC/AJ8coM73AtXAMWAbcHGE/X5jr9uA48C1\n9rV+a8szgB8DTfZ97AOutPeVAW/bxx0C/qujroHHB3QDH7Pb5wF/DzQAR4BngdER6ncpUGXr1wL8\nMsJz/hXwmq3jTuB/htShG/gycBDzs/3fjvsuAf4NOGpf4/8AuRHqE2oM8ENMi6TNvuZowI/523Cq\nJ/g7KwdeBtbZOn8DOAXkO/afY+sTCIhYf6ciEkeZwC7gm4APmAa8Cyy0938KOAxMBH4GbHI8dgHB\nA+ZMzAHvVnu7GHNgegYYZZ+nE/gXYAImeJqAT9j9VwBdwCOYg8LnMAedPHv/G5iDBPY1aoASW///\njgmQcIro3x20guAB9NPAW8B4e7sEmGS3DwPX2+1czEEr9PEBzgP2j4HNtu7jgFeBJyPUbyOwym6P\nAq6L8Jy/BDZgDsBXYIL1NyH7vmrfx1RMoAe6vC7B9PlnY372VbaOAXVEDtx/xATIZMzP8FpbTz/9\nQ8D5POXAGeAWe3s08Gvgfsf+P8D8fcDgfqciEkfzMJ9YnVYBv3Dc/imwH/NPn09kPwF+ZLeLMQem\nyY77jwJLHLdfxhz0wRxYG0Oe702CrQ5nCFQ4tsEcNE5iDn6hAvWIFAKfBA5gfg6hLd4GYCXBgAj3\n+IDAATsDOEHw4A3wN5gWVjhrgOeAKWHuCzxnFuaAepnjvifo3xJwBsj/Bb4e4TUXA//uuB0pBDIx\nn95nhrnPT/QQqAy5/z5MEID5Ob0HzLe3B/M7lSRQd1D6KsJ8Km91XFYBhY59fob5xL/a3h8wD3Nw\nbsZ8av8yptvCqcmxfTrM7bGO26Eh0EDfEHHW+WlHfT+05eEOpNH8G6br5B9t3Z7DDMgC3IHpEqrH\nHNCujeH5JgLnY1pXgfpVYD6Bh/M1zAFxJ/Bn4J4Iz+mj70H3UJj9jji2T2FaIQAXYFoShzBdM+vo\n/3sKZwLmE/y7MewbTmgd/xkTiJMwLcBu4Hf2vnj+TiUBFALp6z3MJ7h8x2U8wXGALOB5YC3wFUzX\nQsAGTLfHRZiuj39ieH8rof/wRcAHEeq8MqTOY4E/htk3llkv/wCUAjOA6cCjtvwtzKfmiZj3GegK\nO4k50AdMcmwfxYTbDEfd8ujfmghosu9lCiZEn6FvKwJMv/lZ+n4qjuUTcuC9PwmcAz6O6db6IrH9\nno4CH2HGLUKF/gyyMD+ncK8f0Aq8DtwJ3IXpCgsYzO9UkkAhkL52Yuajfw0zCJiFOViU2vu/gTmA\n3IPpw11L8O9hHOYf+wxwDeYfe7BTDZ2zdgqBhzF910uAy4GtYR7zT7ZeM+ztXPp2Mzm1YD5xXhLh\n/lJMiyYb8+n5I8z7zcasJ8i1tzvsNcBeTMtoNuaTcrnj+boxLaefEDwoTiE4xhJqCSZEwbSmeuxz\nOJ3DfIoux/yOLsccyAf6WTt/ruMwB+3jti6Phn1Ef92YbsEfYVpkWZhP8qMwA9CjMS2lbMyY0nkx\nPOcG4G5MK2uDo3wwv1NJAoVA+urGfOq/CtNv3YL55D8emAt8FViOOeA8Za8Dfc0PAt/BHFy+hemH\ndoolEJz7vInp927B9HnfQd/up4DNti6/xHRv7CfyvP9TwP/CDDIewxzwnfPix2Pe7zFMt89RTNiB\nGY+os6+xEhMKYA6A3wF2YMYTfhvyPr6OmUX1R/vY7ZgWRjildr8O4P9hQrDe3ud8zr/FHBiPYMYR\nNmLClzD7ht7+NnC1rctrmJlgsYb1f8P8fP+E6aL5LiZg2jG//xcw3T4n6NtdFWntwauYlsVh+7wB\ng/mdSgoajfkH3oOZ4vVdW16O+QPZbS83Ox6zCjMb4B36fkqai/kDqMH0EYo3rKD/YKtE9hTwYrIr\nIeIU6B/0YT7ZzMfMG/+7MPvOwARGNmb2Ri3B5utOTNcCmK6ARYmprqSYFSgEBlICzML8n1yDaS3d\nMuAjROIolu6gU/Z6FKbvMNCMD7dS81ZMc7YL0/StxTTTJ2NmZuy0+63FDMxJ+tOpCwaWg+nGOYHp\nMvl7TNeKyIiIJQQyMZ/umzDTBt+25Q9hBtJ+TnDhz4X0nT52CDNgFVreiKaIecUaggvHpL+3MOMl\nYzGzh55KbnXEa2IJgW7M4OJFmH9mP2a5/DRbfhiz/FxERFzGN4h924EtmFkPlY7yFzAzE8B8wnfO\nc74I0wJoJDhdLlAeuoCISy65pOfdd4e6fkVExLPeJfy6j6iitQQmEOzqGYM5T8xu+i6iuY3glLBX\ngaWY8YNpmGbuTsz0t+OY8YEMzFzozaEv9u6779LT0+Pay+OPPz7g/XffffdQfke9fD4fe/fuTVr9\nU/2i+qvuXq0/kdfLRD+uRLl/MqZPN9Ne1mHOEbIW0xXUg5lv/WW7fzVm9WU1ZiXkgwQHBR/EnJ5g\nDGZ20LahVtqNli9fzrp16/qVZ2RkkJ+fT2ZmMI/b29vp6urqt+/Zs2eZPXs2+fn57N69m6KiooTW\nWUTSX7QQ2I9ZjBJq+QCPeZLwZ1bcRfgTVqW9u+++O2wAfOpTn2LTpk3k5eX1KW9ra2PZsmW89dZb\ntLS0BJK+V2trK8XFxSxYsIDNmzf3e7yISKy0YjiO/H5/v7J77rmHtWvX9inz+XzU19fz+uuvhz2A\n5+XlsWXLFpqamqirq6OwsLDfPgBVVVXk5+fj9/tpa2tLSP3dRPVPHjfXHdxf/+FItW9l6gn91Otm\n9957Ly++2Hfxp8/no7a2dtBdOQ0NDVx77bU0NTX1axkAZGdnU1NToy4iEQ/KyMiAIR7P1RJIkHgG\nAEBRURGHDx+mrq6O887rfz6vrq4uiouL49YqEBFvUAgkQLwDwKmoqIgjR45QVlbGqFGj+t1fVVVF\nYWEhDQ2h3ycjItKfuoPirKSkhIMHD/Ypi1cAhAoMIFdUVPTrIsrIyKCurk7dQyIeMJzuIIVAHI1k\nADg1NDRQUlJCZ2dnn3IFgYg3aEwgRdTU1PS5PX78+IQHAPTtInLq6elh2rRp6hoSkYgUAnGyaNGi\nPl0yEydOpKGhYcQ+hQemlS5YsKBPeU9PjwaMRSQidQfFQVdXF6NHj6a723x7YE5ODu+9915SFnEN\nNE5QWFjIgQMHtLhMJM2oOyjJFi5c2BsAADfccEPSDrSBFkFdXV3gD6NXc3MzJSUlahGISC+1BIap\np6eH7Oxszp0z31U+fvx4GhoaUuLTdqQBY7UIRNKLWgJJdMcdd/QGAMD8+fNT5uAaGDAOPe2EWgQi\nEqCWwDCNHTuWU6fMN3CmUivAqa2tjZKSEpqbm/uUq0Ugkh7UEkiS7u5uTp8+3Xs7lVoBTnl5eRw4\ncEAtAhHpRyEwDLfffnvv9syZM1m/fn0SazMwBYGIhKMQGIYdO3b0TsOcOnVqSrYCnAYKgmXLliWp\nViKSTAqBYXB2BWVnZyexJrGLFAQVFRVaWSziQQqBIVq6dGnv2oCZM2eyevXq5FZoEAJB4AwunWJC\nxJsUAkO0a9eu3u2PfexjKd8VFCovL4+ampo+C8oUBCLeoymiQzRmzBg++ugjcnNz2bt3r2vP1NnQ\n0MC0adP6nGJCZx8VcRdNER1h586d48yZMwC0t7fz6KOPJrlGQ1dUVNTvFBNqEYh4h0JgCD772c/2\nbs+ZM4fnn38+ibUZPgWBiHcpBIZgz549vYPCF198sevGA8JREIh4k8YEhmD06NF0dna6fjwgHI0R\niLiPxgRGUDqNB4QTqUWgVcUi6UkhMEh33nln73Y6jAeEEy4IOjs7FQQiaUghMEj79u3r7SpJl/GA\ncMIFgc4zJJJ+FAKD1NjYCEBubi5PP/10kmuTWEVFRdx00019ynSeIZH0Ei0ERgNvAnuAauC7trwA\n2A4cBF4HnB+HVwE1wDvAQkf5XGC/vc+1R8/AJ+N0HA8I56WXXtJ5hkTSWLQQ+Ai4EbgKmGW35wOP\nYUJgOvBrextgBnCnvV4EPENwxPpZ4D7gMntZFK83MVK+9KUvcfLkSSB9xwNC6TxDIuktlu6gU/Z6\nFJAFtAK3AGts+Rpgsd2+FdgIdAH1QC0wD5gM5AA77X5rHY9xjf379/dup/N4QKhI5xnS+ICI+8US\nApmY7qAm4A3gbeACext7fYHdvhA45HjsIWBKmPJGW+4qge8SLi0tddVZQ+NBM4ZE0lMsIdCN6Q66\nCPgEpkvIqcde0t7ixYu59NJL2b59u2daAU6aMSSSfnyD2Lcd2IIZ4G0CJgFHMF09gW8wbwSmOh5z\nEaYF0Gi3neWN4V6kvLy8d9vv9+P3+wdRxcRat24dAHfddRcbNmzwbBDcdNNN7Nixo7csEAT60nqR\nkVFZWUllZWVcnivaMuMJwFmgDRgD/CvwbeDTwIfAU5hB4Tx7PQPYAFyD6e7ZAVyKaSm8CTyMGRfY\nAvwU2Bbyeil92oj8/PzeT7xLlixh06ZNSa5RcrS1tVFSUkJzc3Of8sLCQgWBSBIM57QR0R40EzPw\nm2kv64AfYKaIbgIuxgwAfw4TFADfAO7FhMcjmOAA04JYjQmTrZhACJXSITBu3DhOnjxJaWmpZ7uE\nAhQEIqkjkSEw0lI6BHJzc7nxxhtZvXq1DnIoCERShU4gNwJWrFhBR0cHnZ2dya5Kyoj0pfXNzc2s\nWLEiOZUSkUFRCMSourqanp4etm3bxsqVK5NdnZQRKQi2bt2qxWQiLqAQGKTS0lJPrBQejHBB0NXV\npVXFIi6gEIjRihUruPjiiz0/IBxJIAhGjRrVW6ZVxSKpTyEQo9bWVu666y4FwADy8vI4ePBgnzKt\nKhZJbZodFKMZM2Zw5swZpk+f7tmFYrHy+/1UVVX1KdOMIZHE0RTRETBx4kSOHj0KeHuhWCw0dVRk\nZGmK6AgIhJMGhqMbaOqouoZEUotaAjGaPXs2hYWFvPTSS/okGyO1CERGhloCI+D48eM899xzOnAN\ngloEIqlPIRCjY8eOUVBQkOxquI6CQCS1KQRicP/993P8+HE+//nP66A1BAoCkdSlEIhBdXU1gE4Z\nMQwKApHUpBCIgc9nvntHM4OGR0EgknoUAjH45je/SUFBgU4ZEQcDBcGyZcuSVCsR71IIxKCrq4t5\n8+YpAOIkUhBUVFTohHMiI0whEAPNDIq/QBBkZ2f3lvX09OjMoyIjTCEQA4VAYuTl5VFTUxNY6AIo\nCERGmlYMx+Dqq6/m2LFjzJgxQyePS4CGhgamTZuG83efkZFBXV0dRUVFSayZiDvoBHIJNmXKFD74\n4ANAJ49LFAWByNDptBEJFuiu0BTRxCkqKqKurq5f11BxcTF+v1/TR0USRC2BGCxevJimpiYqKirU\nFZRg4VoEoJPOiQxELYEE6+zs5Fvf+pYOQCMgXIsAtKBMJFEUAjE4ceIE48aNS3Y1PCMQBOedd16f\ncgWBSPwpBGLQ0dFBTk5OsqvhKUVFRRw5ciTsyuLCwkJNIRWJE4VADNQSSI5IK4u7urrUIhCJE4VA\nDNQSSJ5IQdDZ2akgEIkDhUAUK1eupKWlheXLl+uAkySBIJg4cWKfcnUNiQyfQiCKAwcO0NPTw/bt\n2/VdAkmUl5fHwYMHw3YN6TQTIkMXLQSmAm8AbwN/Bh625eXAIWC3vdzseMwqoAZ4B1joKJ8L7Lf3\nPT3Meo+YUaNGAVoolgoidQ3pfEMiQxdtccEke9kDjAN2AYuBzwEdwI9C9p8BbAD+GpgC7AAuA3qA\nncDf2uutwE+BbSGPT7nFYtXV1cydO5fDhw9rnUCKaGtrY9myZVRUVPRbVLZgwQI2b96s35V4SiIX\nix3BBADACeAvmIN7pBe8FdgIdAH1QC0wD5gM5GACAGAtJkxSXlZWFlOnTtVBJYXk5eWxZcuWsIvK\nqqqqNGAsMgiDGRMoBuYAf7S3HwL2Aj8HAkfICzHdRAGHMKERWt5IMExSmmYGpa6BVhdrwFgkNr4Y\n9xsHvAw8gmkRPAt8x973BPBD4L54VKi8vLx32+/34/f74/G0Q9bR0aE1AiksEAQlJSV0dnb2lnd1\ndVFcXKzuIUlLlZWVVFZWxuW5YulDygZ+BVQAPwlzfzHwGjATeMyWfc9ebwMeBxowA8xX2PLPAwuA\nB0KeK+XGBF577TWee+45fvWrXyW7KjKAtrY2SkpKaG5u7ndfWVkZW7ZsSUKtREZGIscEMjDdPdX0\nDYDJju3bMLN+AF4FlgKjgGmYQeGdmLGF45jxgQzgi8DmoVR4pKkl4A6BmUNlZWX9uoe2bt2q01GL\nRBAtBK4HvgDcSN/poE8B+zBjAguAr9r9q4FN9roCeBAzMwi7/QJmimgt/WcGpaQTJ05oTMAlog0Y\na5xApL9oYwK/I3xQVAzwmCftJdQuTJeRq2hg2H0C4wSh30sQWFimbysTCdKK4Sg2bNjAK6+8QllZ\nmboTXCTS6aj1bWUifSkEomhubua9996joqJCp41wmcDpqMONE6h7SMRQCESRlZUF6LQRbjXQOEFg\nGqlaBeJl+o7hKJYvX051dTU7duzQXHOXa2ho6LeeICA7O5uamhqNFYgr6TuGE6i7u5tHHnlEAZAG\nBuoeCrQK9u3bl6TaiSSHQiCK06dPM2bMmGRXQ+LE2T0UOmgMMHv2bHUPiacoBKJQCKQnZ6sgVFVV\nFfn5+QoD8QSFQBQKgfQVaBXs3bs37P2aQSReoBCIQiGQ/mbNmkV9fX3Y7iHNIJJ0pxCIQiHgDc7u\nocC3yTmpi0jSlUIgCoWAdwS6h5qamsLOIAJ1EUn6UQhEoRDwnmgziAJdRPPnz1erQFxPIRBFS0sL\nS5cu1bmDPChaF9Hvf/97dRGJ62nFcBSZmZm9Z6JcsmQJmzZtSnKNJBkCX26/devWiPvoW8wkWYaz\nYlghMICenh4yM01jqbS0lO3bt+sf3OP27dvHVVddxUB/pwoDGWk6bUSCdHZ24vP5WLJkiQJAADOd\n9NixYxG7iEAzicRd1BIYQGtrK8XFxbS3tye7KpKCAl1EO3bs4MyZMxH3U8tAEk0tgQTRzCAZSOiU\nUrUMxI0UAgM4ffo0559/frKrISlOYSBuphAYgFoCMhixLDaDYBgUFBRo0ZkknUJgAAoBGQrnYrNJ\nkyaRnZ0ddr/AmJPP51PrQJJGITAAhYAMR1FREYcPH6a5uXnAbqJz586pdSBJoxAYgEJA4iG0m6iw\nsDBiV5FaBzLSFAIDUAhIPDnDIFpXkVoHMlIUAgNQCEiihHYVRQoDCLYOMjIyFAgSd1osNoAFCxZw\n4MABrr76ajZs2KDFPpIwgYVnb731Fq2trXR1dUV9jM/nY9euXcyaNWsEaiipTOcOSpDLLruM2tpa\nQCePk5HjDISWlpYBz1MUkJ+fz+7duykqKhqBGkqqUQgkyOWXX86BAwd08jhJmoaGBq699lrOnj3L\nhx9+GFMggELBaxJ52oipwBvA28CfgYdteQGwHTgIvA44j46rgBrgHWCho3wusN/e9/RQKjvSbr/9\ndq644goFgCRNYOygpaWldzC5oKAg6uOc4wiaaSQDiRYCXcBXgSuBa4GvAFcAj2FCYDrwa3sbYAZw\np71eBDxDMJ2eBe4DLrOXRfF6E4m0bNkyBYCkhEAgfPjhh9TX1zNp0iQmTJgQcbppgHOmkUJBQkUL\ngSPAHrt9AvgLMAW4BVhjy9cAi+32rcBGTHjUA7XAPGAykAPstPutdTwmZZ06dUqzgyQlhWshTJgw\nQaEggzaYKaLFwBzgTeACoMmWN9nbABcChxyPOYQJjdDyRlue0jRFVNzAGQihoRAtEKB/KGRkZJCf\nn6+pqB7hi3G/ccArwCNAR8h9PfYSF+Xl5b3bfr8fv98fr6ceNIWAuFEgFKDvwDJAe3t7TNNP29ra\nKC4u7lOmwebUUVlZSWVlZVyeK5YQyMYEwDpgsy1rAiZhuosmA822vBEzmBxwEaYF0Gi3neWN4V7M\nGQLJphAQt3MGAvSdftrd3R1zKEBwsNlJaxWSI/QD8re//e0hP1e07qAM4OdANfATR/mrwN12+26C\n4fAqsBQYBUzDDADvxITFccz4QAbwRcdjUpZCQNKN89QVLS0tvSuWCwsLKSgoIDs7O6YupICzZ88y\ne/bs3m6kwEXjDO4RrSVwPfAFYB+w25atAr4HbMLM9qkHPmfvq7bl1cBZ4EGCXUUPAquBMcBWYFsc\n6p9Q+lIZSXeBUHAK7UICBrVGAfqOM4RSt1Jq0WKxAVx33XX84Ac/4Prrr092VUSSKjQYuru7OXbs\nWEJeKysri/nz5+t7mQdhOIvFYh0Y9iR1B4kYoWMLEL7FALEPPkcyUCsiQEERPwqBASgERCILFwzQ\nf/C5u7ubjo4OcnJyaG1tHVS3UiSxBIWTuqAiU3fQAMaNG8fHP/5xCgoKdBZRkTiI1HoIGOzYQyK5\naeaTTiCXINnZ2b1/rDqLqEjiRQsJGH5300gaqRaIQiBBfD4f586d01lERVJIoLtp586dtLe3k5OT\nQ2Zm5NnuqdS6cMrMzOSGG26Iy7iGQiBBsrKyuO2223jhhRcUACIu1dDQwDXXXENra2vUwIDEznwK\nJx69DAqBBOjq6mLMmDEDNktFJD3F0i0F9Bn0Dg2XWFogs2bNoqqqSi0Bh5QJgePHjzNlyhQ6OkJP\nlSQiEl24tRWBwPD5fJSWlrJ+/fq49DJonUACaHqoiAxHpCm0qWYwp5L2FIWAiHiBQiACfaGMiHiB\nQiACtQRExAsUAhEoBETECxQCESgERMQLFAIRKARExAsUAhHoC2VExAsUAhE8++yzVFVVUVZWpq/I\nE5G0pRCIoLGxkebmZioqKli5cmWyqyMikhAKgQh8PrOYurS0lOeffz7JtRERSQyFQAR33HEHl19+\nuU4hLSJpTSEwgLvuuksBICJpTSEQgaaIiogXKAQiUAiIiBcoBCJQCIiIFygEIlAIiIgXKAQiUAiI\niBcoBCJQCIiIFygEItCXyoiIF8QSAr8AmoD9jrJy4BCw215udty3CqgB3gEWOsrn2ueoAZ4eco1H\niFoCIuIFsYTAi8CikLIe4EfAHHupsOUzgDvt9SLgGSDD3vcscB9wmb2EPmdKUQiIiBfEEgK/BVrD\nlGeEKbsV2Ah0AfVALTAPmAzkADvtfmuBxYOs64h6//33ue+++3QWURFJa8MZE3gI2Av8HAicW+FC\nTDdRwCFgSpjyRluesjo7O9m5c6fOIioiac03xMc9C3zHbj8B/BDT1TNs5eXlvdt+vx+/3x+Ppx20\n7u5uQGcRFZHUU1lZSWVlZVyeK1yXTjjFwGvAzCj3PWbLvmevtwGPAw3AG8AVtvzzwALggZDn6unp\n6YmxSonT09NDVlYWt99+Oy+88IJOIiciKS0jIwNiP573MdTuoMmO7dsIzhx6FVgKjAKmYQaAdwJH\ngOOY8YEM4IvA5iG+dsKdOXMGn8/Hyy+/rAAQkbQWS3fQRsyn9gnA+5hP9n7gKswsoTrgy3bfamCT\nvT4LPGj3wW6vBsYAWzGthJR08uRJxo4dm+xqiIgk3JCaDwmUEt1Bhw4dYt68eTQ2Nia7KiIiUSWj\nOyitqSUgIl6hEAjj1KlTnH/++cmuhohIwikEwlBLQES8QiEQhloCIuIVCoEw1BIQEa9QCIShloCI\neIWmiIZxww03UFtby5w5c9iwYYMWjIlIStMU0Tg7fPgwR44c0cnjRCTtKQTCyMrKAnTyOBFJfwqB\nMD7zmc9w5ZVXsn37dnUFiUhaUwiE0d3dzb333qsAEJG0pxAIQ1NERcQrFAJhaIqoiHiFQiAMtQRE\nxCsUAmGoJSAiXqEQCEMtARHxCoVAGGoJiIhX6LQRYeTm5jJ9+nQmTpyo00aISMobzmkjFAJhjBo1\niq6uLgCWLFnCpk2bklwjEZHIFAJxlpWVRXd3N6WlpVo1LCIpTyeQiyNnC0ABICLpTi2BEMeOHeOS\nSy6htbU1qfUQEYmVWgJx1N7eTm5ubrKrISIyIhQCIY4fP8748eOTXQ0RkRGhEAihloCIeIlCIMST\nTz5JdXU1ZWVltLW1Jbs6IiIJpRAIUV9fz7Fjx/TVkiLiCQqBEHaUXV8tKSKeEEsI/AJoAvY7ygqA\n7cBB4HXAOZl+FVADvAMsdJTPtc9RAzw99ConVllZmb5aUkQ8I5YQeBFYFFL2GCYEpgO/trcBZgB3\n2utFwDME564+C9wHXGYvoc+ZEk6cOMFDDz2kABART4glBH4LhK6cugVYY7fXAIvt9q3ARqALqAdq\ngXnAZCAH2Gn3W+t4TEo5evQoEyZMSHY1RERGxFDHBC7AdBFhry+w2xcChxz7HQKmhClvtOUp5w9/\n+ANPPPGEZgeJiCfEY2C4x17SQltbG3v37tXsIBHxBN8QH9cETAKOYLp6mm15IzDVsd9FmBZAo912\nljeGe+Ly8vLebb/fj9/vH2IVB6+7u5vOzk5As4NEJHVVVlZSWVkZl+eK9YRDxcBrwEx7+/vAh8BT\nmEHhPHs9A9gAXIPp7tkBXIppKbwJPIwZF9gC/BTYFvI6ST2B3AcffMBVV12F3+/n+eef1+CwiLjC\ncE4gF0tLYCOwAJgAvA/8D+B7wCbMbJ964HN232pbXg2cBR4k2FX0ILAaGANspX8AJN0DDzzA6dOn\nOXHiRLIKVqGRAAAFPklEQVSrIiIyInQqaYfp06dTU1MD6BvFRMQ99M1icTJ27FhOnTpFbm4ue/fu\npaioKGl1ERGJlb5PIA7uv/9+Tp06BZgziT766KNJrpGISOKlXEsg2RUAyMvLo66uTgPDIuIKagnE\n2XXXXacAEBFPUEsgxJVXXsnvfvc7hYCIuEaip4iOqMLCQrq7uwGzeKujo4OcnBwyMwdutAx338zM\nTEpLS1m/fr0CQEQ8I+VaAsmcHSQi4kYaExARkSFRCIiIeJhCQETEwxQCIiIephAQEfEwhYCIiIcp\nBEREPEwhICLiYQoBEREPUwiIiHiYQkBExMMUAiIiHqYQEBHxMIWAiIiHKQRERDxMISAi4mEKARER\nD1MIiIh4mEJARMTDFAIiIh6mEBAR8bDhhkA9sA/YDey0ZQXAduAg8DqQ59h/FVADvAMsHOZri4jI\nMA03BHoAPzAHuMaWPYYJgenAr+1tgBnAnfZ6EfBMHF4/pVRWVia7CsOi+ieXm+vv5rqD++s/HPE4\nCGeE3L4FWGO31wCL7fatwEagC9OCqCUYHGnB7X9Iqn9yubn+bq47uL/+wxGPlsAO4C3gS7bsAqDJ\nbjfZ2wAXAoccjz0ETBnm64uIyDD4hvn464HDwERMF9A7Iff32EskA90nIiIJFtqVMxyPAycwLQI/\ncASYDLwBXE5wbOB79nqbfcybjueoBS6JY51ERLzgXeDSkX7R84Ecuz0W+D1mxs/3ga/b8scIHvRn\nAHuAUcA0TKXjGUIiIjKCpmEO6nuAP2Omf4KZIrqD8FNEv4H5tP8O8OkRq6mIiIiIiKSuRZgWQg3B\n7qRU8gvMbKf9jjI3LYybihmfeRvTcnvYlrvlPYzGjB/tAaqB79pyt9Q/IAuzuPI1e9tN9a/HvYtD\n84CXgb9g/n7m4Z66l2B+5oFLO+b/1y31j0kWppuoGMjG/KNfkcwKhXEDZlGcMwS+D3zNbn+d/uMf\n2Zj3VEvyF8ZNAq6y2+OAA5ifsZvew/n22gf8EZiPu+oP8HfAeuBVe9tN9a/DHHic3FL/NcC9dtsH\n5OKeujtlYmZkTsWd9Y/obzCzhQIeIzibKJUU0zcE3iG4DmISwSmyq+jbmtkGXJvoyg3SZuA/4c73\ncD7wJ+BK3FX/izDjZTcSbAm4qf51wF+FlLmh/rnAf4Qpd0PdQy0Efmu341L/VEmHKcD7jttuWUjm\n1oVxxZhWzZu46z1kYj7hNBHs2nJT/X8MPAp0O8rcVH+3Lg6dBrQALwL/DvwMM6PRDXUPtRRz5gWI\nU/1TJQTSYdGYWxbGjQNeAR4BOkLuS/X30I3p0roI+ATmE7VTKtf/M0Azpk830tToVK4/mMWhc4Cb\nga9gukidUrX+PuBqzPnKrgZO0r+nIVXr7jQK+C/AS2HuG3L9UyUEGjF9XAFT6ZtkqaoJ0wwDszCu\n2W6Hvp+LbFmyZWMCYB2mOwjc9x7ADIxtAebinvpfhzmvVh3mk9wnMb8Ht9QfTF80mE/V/4I595cb\n6n/IXv5kb7+MCYMjpH7dnW4GdmF+/uCOn33MfJjFY8WYtEvFgWHoPybgpoVxGcBaTJeEk1vewwSC\nsx/GAL8BbsI99XdaQHBMwC31d/vi0N9gzmwMUI6pt1vqHvBL4G7HbbfVP6qbMTNWagkuPEslG4EP\ngDOY8Yt7cNfCuPmY7pQ9BKeaLcI972Empj93D2aa4qO23C31d1pAcHaQW+rv9sWhszEtgb3AP2MG\ni91SdzDBe5RgEIO76i8iIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIon0/wEOxUdzQSKB8wAAAABJ\nRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10f681650>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 2. Derive the AIF ##"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# turn Sb into Cp\n",
      "T1p = 1.440\n",
      "R1p = 1 / T1p\n",
      "Hct = 0.45\n",
      "S0 = data_aif[:4].mean()\n",
      "R1_eff_aif = dcemri.dce_to_r1eff(data_aif, S0, R1p, TR, flip_angle)\n",
      "Cb = dcemri.r1eff_to_conc(R1_eff_aif.flatten(), R1p, Rel)  \n",
      "Cp = Cb.flatten() / (1.0 - Hct)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "converting DCE signal to effective R1\n",
        "converting effective R1 to tracer tissue concentration\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "if plotting:\n",
      "    figure(1)\n",
      "    clf()\n",
      "    plot(t, Cp, 'k.-')\n",
      "    xlabel('time (s)')\n",
      "    ylabel('$C_p$ (mM)')\n",
      "    title('derived AIF')\n",
      "    savefig('qiba/qiba_v6_aif.pdf')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEZCAYAAACJjGL9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH59JREFUeJzt3XuUlPWd5/F30Vf6AtUMNxGhiRozccUILRhDsDCa0dZ1\nNjg9S4LrddNnZ7xkMqsT455dcc7R3GZnxnV3zDBjEseISWuiIUc72kwsYgAVowhEQEQkXBvoi1yb\nrr7sH7+nup+urqqueqqqn0t/XufUoeqpp6u+hVKf/l0fEBERERERERERERERERERERERERERCYwi\ntwsQKaAfAfOBVx3+/FbrtidfBVmimH977+T5dUUKYpzbBYgUUL91c+o/AL/JUy12mdR1K9AH/HnC\n8Qiw1/Y4CpwGjttuC/NQo8gABYUEXcjBzxTnvYrs3QJsAW4e4bx+4E6g2nZ7o7ClyVijoJAguQR4\nGzgG/AQoT3j+emAT0AGsAy6yPfcR8DfAZsxv5UXWsSuBGcApoCbhvY4w2H17O/Ae0A78CphlO/dq\nYDvQCTyGCa90ATYb+Bxwm/Wz09KcK1JwCgoJilLgBeBJzBf6s8CNDHbxXAI8AXwVmAT8M7AaKLG9\nxjLgWiAM9Np+9gCwwXq9uK9Y79EL/CnwTeBLwGTgNeAZ67zJwM+AB4A/AnZhQiBd19PNwFpM6L0F\nLB/hsztpNYmIjDmLgf0Jx9YBf2vdf9x2P2478Hnr/m7MuIDdbkyLAuAO4N+t+yHgD8Ai63EzpkUR\nNw44iWlV3AysT3jdvQnnJ9oJNFr3/wrTCoqLMHyM4iSmldSBCRaRvFKLQoJiBsODwj5baTbw3xn8\nQu0AZlo/F2f/Ak70c+CzwHRMKPUBv7W99qO2122zjp8NnAXsS3itdO/zOaDWej+A5zBdZBenOL8f\nuBvTiqoB6tK8togjCgoJioOYL2a72bb7fwAeZvALtQaoAn5qOyddd1AH8ArwnzHdTs/YnvsDpgVg\nf+1KTHfVQeAc27mhhMeJbrHO2WL97EbbcRERyUEJpgVxj3V/KdDNYHfTfMwX+gLMF3ElcB0mLGBo\nNxMpjv05ZtzgKEMHwv8T5ov909bjiUCDdX8yZnD9S5jZVF8DYiTveirHDHjfBky13f4SOIQZOI8w\ntEXyKqZbTEREMjCfobOenmHouMSfAG9iWgcHMK2JSuu5TIKi3HrtLUne+ybMjKmPMYH0rwnvu4PB\nWU+vkjwolmG6zxIXwo7HhFM9cIX1+nGpXkvE934AtDL0H9wkoAV4H9PED7tQl4iIJHBrjOKHwDUJ\nx+7HBMUnMbNL7h/tokRExFtqGdqi2M7gwqLp1mMREXGZl2Y9TcN0R2H9qdWoIiIe4KWgsMt1MzcR\nEckTL2x+FteK6XI6hFmkdDjZSeeee27/rl27RrMuEZEg2AWc5+QHvdSiWM3goqJbMPv2DLNr1y76\n+/t9e3vwwQddr2Gs1u/n2lW/+ze/1w+c6/TL2a2geAaz/80FmMVDtwHfxuyU+T5m7vq3XapNRERs\n3Op6+nKK41eNahUiIjIiL3U9jQmRSMTtEnLi5/r9XDuofrf5vf5c+HEf+36rv811d9xxBxs3bmTm\nzJmsWrWKcFiLyUXEm0KhEDj8zldQ5GDu3Lls2WLWDDY0NNDU1ORyRSIiyeUSFOp6ykFpaSkAFRUV\ndHR00NnZ6XJFIiL5p6DIwW233QbAqVOnWLNmDY2NjSP8hIiI/ygocnDy5MmB+3V1daxcudLFakRE\nCkNBkYPDh83i8c985jO0tLRoMFtEAklBkYPVq1dTUlIypGUhIhI0CooctLe3E4vF2Llzp8YnRCSw\nND02B1OmTOHo0aNMmTKF999/X11PIuJZuUyP9dLusb5zwQUXUFVVxfz58xUSIhJY6nrKQW9vL8uW\nLaOrq8vtUkRECkZBkYPTp08TjUZ57bXXqK+v14I7EQkkBUUOurq66Ozs5NixYzQ3N2tAW0QCSUGR\ng9OnT1NZWQlowZ2IBJeCIgddXV08/PDDVFZWasGdiASWgiIHp0+fZtasWUycOFEhISKBpaDIwenT\npwmHw5w+fdrtUkRECkZB4VBPTw99fX1MmDBB02NFJNAUFA51dXUxfvx4ysvL6erqwiurxUVE8k1B\n4VBXVxfl5eUUFRVRXFxMLBZzuyQRkYJQUDh0zz33cPz4cerr6ykrK9M4hYgEloLCoV27dtHd3U1z\nczPd3d0apxCRwFJQOBS/XnZdXR1Tp05VUIhIYCkoHHrkkUeYOHEiLS0tVFRUKChEJLAUFA6Vl5dz\n/vnnEw6HOXr0KMuXL9fGgCISSAoKh7q7uwe6n7q7u/nd736njQFFJJAUFA7Zg6KkpATQxoAiEkxe\nvMLdN4GbgD5gC3AbcMbVipKIxWIDAXHJJZcQi8VYvXq19nwSkcDxWouiFvgqMA+4CCgClrlZUCr2\nFkV1dTVf//rXFRIiEkhea1EcA2JABdBr/bnf1YpSsAdFfBsPEZEg8lqLoh3438AfgANAJ7DG1YpS\nsHc9jR8/XkEhIoHltaA4F/grTBfUDKAKWO5mQamoRSEiY4XXup7qgPVAm/X458DlwNP2k1asWDFw\nPxKJEIlERqc6m8Sg0F5PIuIl0WiUaDSal9fyWlBsB/4nMB7oAq4C3kw8yR4UbrF3PalFISJek/hL\n9EMPPeT4tbzW9fQu8G/AW8Bm65gnFybYWxQaoxCRIPNaiwLgu9bN07q7u4e0KI4cOeJyRSIiheG1\nFoVvxGIxDWaLyJigoHBIs55EZKxQUDhk73rSGIWIBJmCwqHEridNjxWRoFJQOKSuJxEZKxQUDtm7\nnr7//e/z29/+VhcuEpFACrldgAP9/f39btfApz71Kfr6+jjvvPPYu3cvW7duBaChoYGmpiaXqxMR\nGSoUCoHD73wFhUPTp0+ntbUVgEmTJtHe3k5dXR0tLS3ablxEPCeXoFDXk0Pjxpm/urq6On784x8z\nYcIEhYSIBJIXV2b7wsKFCzl06BDNzc0cPnyYadOmKSREJJDUonBo3Lhx3HvvvYTDYcrKyjTrSUQC\nS0HhUCwWo7jYNMjKy8s5c8Zzl/UWEckLBYVDPT09A9Nj1aIQkSBTUDiU2KJQUIhIUCkoHOrp6RkI\nirKyMmKxGF6Ytisikm8KCofsXU+hUIiSkhKNU4hIICkoHLJ3PYEGtEUkuBQUDtlbFKABbREJLgWF\nQ2pRiMhYoaBwyD6YDWpRiEhwKSgcSux6UotCRIJKQeFQsq4ntShEJIgUFA4ltij27t1LY2OjLl4k\nIoGjoHAosUVx5swZNm3aRHNzM42NjS5WJiKSXwoKhxJbFPH7dXV1rFy50q2yRETyTtejcChx1tPC\nhQtpa2vjpZde0nUpRCRQFBQOJXY9VVdXs2zZMoWEiASOup4c0vRYERkrFBQOJbYotOBORILKi0ER\nBp4DtgHvAZe5W85wfX199PX1UVRUNHBMLQoRCSovjlE8CrwE/Bmmvkp3yxmut7eX4uJiQqHQwDEt\nuBORoPJaUEwEPg/cYj3uAT52r5zkErudwHQ9qUUhIkHkta6nOcAR4IfA28C/ABWuVpRE4kA2qEUh\nIsHltRZFMTAPuAvYCPwjcD/wv+wnrVixYuB+JBIhEomMWoEAd911F6dOnaK+vp5Vq1YRDofVohAR\nT4lGo0Sj0by8VmjkU0bVdGADpmUBsAgTFNfbzul3+9rUl19+ORs2bACgoaGBpqYmHn/8cTZv3szj\njz/uam0iIslYY6qOvvO91vV0CNgLfNJ6fBXwe/fKSa6srAwYul2HWhQiElReCwqAu4GngXeBucAj\n7pYz3KOPPsr48eNpaWkZWImtMQoRCSqvjVGACYhL3S4inYqKCmbMmDFkuw4FhYgElRdbFJ6XbHrs\nypUriUajuh6FiASOgsKBZNNjDx48SEdHh65HISKBo6BwIFmLYvz48YCuRyEiwaOgcCDxWhQA3/rW\nt6ipqRkywC0iEgQKCgeSdT1NnjyZs88+WyEhIoGjoHAgWdeTZj2JSFApKBxI1qLQgjsRCSoFhQNq\nUYjIWKKgcCDV7rFqUYhIECkoHEg260mXQhWRoFJQOJCs66m0tJRYLEZfX59LVYmIFIaCwoFkXU+h\nUIjS0lK6u7tdqkpEpDAUFA4ka1GABrRFJJiyDYpyoKwQhfhJshYFaIqsiATTSEExDlgKPAvsB3YD\ne6z7zwFfwntXySs4tShEZCwZKSiiwHzg74BPAGdhLlf6CevYpcDaAtbnSclmPQF0dnbS0NCgrcZF\nJFBGunDR1UCyvpQzwOvWbcx1RaXqeurp6WHjxo0ANDY20tTUNNqliYjk3Ugtikw63Mdcp3yqrqd4\neGircREJkpFaFL8E+kk+DtEP3JD3inwgVYvi4osvBuAXv/iFdpEVkcAYKSguA/YBzwBvWMfiodFf\nqKK8LlWLIhwOc/vttyskRCRQRgqKszDjFF+2bi9iQuP3Ba7L01INZldWVnLy5EkXKhIRKZyRxih6\ngGbgZkzr4gPMLKe7ClyXp6XqeqqoqODUqVMuVCQiUjgjtSjALLK7DlgG1AKPAs8XsCbPS9X1pBaF\niATRSEHxFHAh8BLwt8CWglfkA83NzZw6dYq1a9eyatWqgTEJBYWIBNFIXU/LgfOBrwHrgeO227HC\nluZd7e3tfPjhhzQ3N9PY2DhwvKKiQkEhIoEzUotCmwYmUVRUBAxfL1FZWUlbW5tbZYmIFEQmYxRg\ntup4ADNGEf+ZfmBuAWryvMWLF7N9+3ZaWlqGTIVV15OIBFGmQfE0cC+wFRjzV+YZN24cd99997D1\nEgoKEQmiTLuWjgCrgQ+Bj2y3QikC3sGsDPecWCyWdHrsj370I15++WVtCigigZJpUDwEPIFZdHej\ndVtaqKIwg+fv4dHV36mCorW1lba2tmGD3CIifpZp19MtwAXW+faup5/nvSKYCdQDDwN/XYDXz1mq\noKisrAS0KaCIBEumQVEHfIrR+Q3/H4D7gAmj8F6OpAqK73znO9xwww3DBrlFRPws06BYD3yawu/x\ndD1wGDM+EUl10ooVKwbuRyIRIpGUpxZEqqCYPn06U6dOVUiIiOui0SjRaDQvr5XpZUy3A+diLoUa\nv/5EIabHPgL8F8weU+WYVsXPMHtNxfX397s7dBGJRHjwwQdZsmTJkOP79+9nwYIF7N+/36XKRESS\nC4VC4PDS1Zm2KK5x8uIOPGDdAK7ATMm9OfXp7kg3RqHpsSISNCMFRQjTcvgog3MKwVeznrSFh4gE\n0UjTY6OYgeVPJnnuAuAbmG3HC2EtHr2CXqqgKC0tHXheRCQoRgqKLwJtwP8DDgLvAzut+/8XaAWu\nKmSBXpQqKEDdTyISPCN1PZ0BfmDdioDJ1vGjQG8B6/K0dEHR3d3NtddeS01NzZAtyEVE/Cqb3WF7\nMS2IVsZwSED6oOjv7+f111/X6mwRCQxtI+5AuqAoKysDtDpbRIIj0+mxYpMuKObNm0dPTw+rV69W\nt5OIBIKCwoF0QTFlyhSWLl2qkBCRwMil6+ka4FagND+l+Ee6oAiHw9piXEQCJZeg+Bh4HbOz7Jii\noBCRsSSXrqdrgE7gGDAJaM9LRT6goBCRsWSkFsUtwBeAiUme+yXmEqmngP+T57o8q7+/n97eXoqL\nk2fsiy++yJNPPqmr3IlIYIwUFMcwV7JrAKqsY1djWiJvYbYEfwb4r4Uq0GtisRjFxcXxnRiHaW9v\n5+DBg1pHISKBMVLX0wTgzoRjazGXRH0Js70HQFee6/KsdN1OANXV1YDWUYhIcIzUokjW5dQNPIW5\nXOmYM1JQPPbYY1RVVekqdyISGCMFxRTMQHUyZXmuxRfuvvtuTp06lXIM4vzzz6eoqEghISKBMVJQ\n/BPwU8yAtl0I+OOCVORxH3zwAT09PSnHIMLhMMeOHWPx4sUa0BaRQMjksnifAH4MVGOuT3EauAz4\ne+CFglWWmquXQo1EIqxdu5a6urqU3UtlZWV0d3cD0NDQQFNT02iXKSIyRC6XQs3mhy4HPou5nvWL\nwAdO3jAPXA2Kt99+m8WLF7Nv376U3UuTJk2io6MjbZiIiIym0bhmNsB66zamjR8/npkzZ6b98r/u\nuuvYtm2bQkJEAkHbjGdppFlPANu3b6e1tZWvfOUrGqMQEd9TUGQpk6Do6upi3759WnQnIoGgoMhS\nJkExbdo0AKqqqujo6FCrQkR8TUGRpUyC4rnnniMUCnHixAnWrFnDrbfeOjrFiYgUgIIiS5kERTgc\npqioaODx+vXr1aoQEd9SUGQpk6AA0+0Ud+TIEbUqRMS3FBRZyjQo6urqhjzesGGDWhUi4ksKiixl\nGhTPPvssU6dOHXh8+PBhtSpExJcUFFnKNCjC4TA7duwYcq5aFSLiRwqKLGUaFGDC4oorrhh4fPjw\nYa2rEBHf8VpQnAO8Cvwe2Arc4245w2UTFGC6oEpLSwGYMGEC3/ve9wpVmohIQXgtKGLA14ELMTvU\n3onHtjPPNijC4TDz588H4NixY9x3332FKk1EpCC8FhSHgE3W/RPANmCGe+UMl21QAAMbA5aVlWml\ntoj4jteCwq4WuAR4w+U6hnASFKtWraK4uJgzZ85opbaI+E4224yPpirgOeBrmJbFECtWrBi4H4lE\niEQio1WX4xZFVVXVQEvC2hdeRKRgotEo0Wg0L6/lxaAoAX6Guape0ivo2YNitDkJCjAL8NasWUNZ\nWRknTpygs7NT16oQkYJJ/CX6oYcecvxaXut6CgFPAO8B/+hyLUk5DYpnn31W3U8i4kteC4rPATcB\nS4B3rNs1rlaU4Pnnn+epp56ivr4+q0HpcDhMZWXlwOP44rs77riD0tJSQqHQwG3RokUa8BYRz/Bj\nZ7mr18yeNWsWe/fuBaChoYGmpqaMf/bqq69mzZo1A49DoRCpPsvUqVPZsWOHuqdEJC9yuWa2giJL\ntbW17Nmzh7q6uqyvid3Z2cmUKVPo6enJ6PySkhJ27tzJ7NmznZYrIgLkFhRe63ryvCuvvJJ58+Zl\nHRJgup9SzdCKdzvZxWIx5syZw549e5yWKyKSMwVFlvr7+7nzzjsddwkl7ioLsGTJEtrb22lvbx/2\nXH9/v8JCRFyloMhSV1cX5eXljn8+vqtsfX09M2bM4KOPPuLXv/414XB4yHP21kV/fz+1tbVEIhEN\ncovIqNMYRZaWLl3KTTfdxNKlSwv6Pnv27GHOnDnDBrs1biEiTmiMYhR1dXVRVlZW8PeZPXs2u3fv\nTjpuodaFiIwmBUWWcu16ykY8LJIF09q1a6mpqdGaCxEpOAVFls6cOTNqQQEmLA4dOjRs3CJu3bp1\n1NTUUFxcrFaGiBSEgiJLo9X1ZBcOh3nxxRdTti4Aent7B1oZoVBIwSEieaOgyNJodj0lircurr76\n6hF3oE0MjpKSEjZv3jxKlYpIkCgosjTaXU+JwuEwr7zyCu3t7dTX1zN16tSMNins6enh4osvJhQK\nMWnSJK3LEJGMKSiy5EbXUzLx7qjW1lYOHz48EBqTJ08eMTg6Ojqora1VF5WIZETrKLI0efJktm3b\nxpQpU1yrYSSdnZ0sX76ct956i56eHtrb2zP+2ZqaGt555x2t0xAJGG0KOIqqq6s5cOAA1dXVrtWQ\nrT179nDZZZfR09NDW1tbyh1rExUVFbFo0SJeeOEF7WIr4nNacDdKGhsbOXHiBA0NDb7qqpk9ezYH\nDx7kyJEj7N69m+nTp2fURZU4IK7xDZGxSUGRhR07dgDw8ssv09jY6HI1zthDI3FsIxP28Q1daElk\nbFBQZKG0tBQw179euXKly9Xkzj4gfuTIEd59911KS0uZNGlSVpd7jS/6s4eHWh4iwaExiixs3bqV\nSy+9lIMHDwa+z94+IN7X1weQ1fhGOhowFxl9GsweJe+99x433ngj27Ztc+X93WYfFO/r66OjoyMv\nwZGMwkQkvzSYPUpOnDjhq9lO+WYf32hraxuy6G/y5MkDt5FWjWci2VhIqptWnYsUVrHbBfjJ8ePH\nqaqqcrsMz4iPcSSytzwA+vr6OH78ONXV1QVphcRXnedKrRiR5BQUWRjrLYpMxVseySSGSKKPP/6Y\nWCxWyPJSirdivCAcDrNu3TrOO+88iouLGTdOjX9xj4IiC2pR5C5diEDyQfRU+vr6slp17iednZ1c\neOGFbpchLvPKolcFRRbUoii8VN1ZqYzUQomzd38l++08XzO6RPIpvui1sbGRpqYm1+rQrKcsLFy4\nkAMHDnDRRRexatWqwE+RHUsSZ3SlCxW7Qpyr0BK7uXPnsnbt2py/bzQ9dpTMmjWLvXv3AtDQ0OBq\nwktw7dmzhwULFtDR0eFaWOlc988dN24cdXV1PP3003n5pVRBMUrOPvtsDhw4QF1dHS0tLWpRiIhv\nKChGyZIlS+jt7WX16tUKCRHxlaAtuLsG2A7sBL7hci1DHD16lMcee0whISJjitdaFEXADuAqYD+w\nEfgyYN8zw7UWRXl5OfPmzSMcDmswW0R8JUgtigXAB8BHQAz4CfCnbhYUt2/fPnp7e9mwYQPNzc2+\n3WZcRCRbXguKs4G9tsf7rGOuu/TSSwfm6s+dOzcQ24yLiGTCawvuMupTysemc7k4evSoup1EZMzw\nWlDsB86xPT4H06rwlPXr17tdgohIWtFolGg0mpfX8tpgdjFmMPsLwAHgTZIMZrtQ14DXXnuNRYsW\nuVmCiEjWchnM9lqLoge4C3gZMwPqCYaGBMCI13cuxErK4uJiXn/9dW1BLSJjjtdaFJlwbXqsiIhf\nBWl6rIiIeIyCQkRE0lJQiIhIWgoKERFJS0EhIiJpKShERCQtBYWIiKSloBARkbQUFCIikpaCQkRE\n0lJQiIhIWgoKERFJS0EhIiJpKShERCQtBYWIiKSloBARkbQUFCIikpaCQkRE0lJQiIhIWgoKERFJ\nS0EhIiJpKShERCQtBYWIiKSloBARkbQUFCIikpaCQkRE0lJQiIhIWgoKERFJy0tB8T1gG/Au8HNg\norvliIgIeCsoXgEuBC4G3ge+6W45hRGNRt0uISd+rt/PtYPqd5vf68+Fl4KiBeiz7r8BzHSxloLx\n+/9sfq7fz7WD6neb3+vPhZeCwu524CW3ixARESge5fdrAaYnOf4A8Evr/v8AuoFVo1WUiIikFnK7\ngAS3Al8FvgB0pTjnA+Dc0SpIRCQgdgHnuV1Erq4Bfg9MdrsQEREZ5KUWxU6gFGi3Hm8A/tK9ckRE\nREREJFCuAbZjWh7fcLmWVH4AtAJbbMcmYQbx38esFQnbnvsm5vNsB744SjWmcw7wKqYLcCtwj3Xc\nL5+hHDO1ehPwHvAt67hf6gcoAt5hcHKHn2r/CNiMqf9N65if6g8Dz2EW/r4HLMQ/9V+A+XuP3z7G\n/Pv1S/15UYQZxK4FSjBfBH/sZkEpfB64hKFB8V3gb6z73wC+bd3/NOZzlGA+1we4P115OvAZ634V\nsAPz9+ynz1Bh/VkMvA4swl/1/zXwNLDaeuyn2ndjvpjs/FT/k5ip+WD+/5mIv+qPGwccxPzi58f6\nHfss8Cvb4/utmxfVMjQotgPTrPvTrcdg0tzeMvoVcFmhi8vSC8BV+PMzVAAbMav9/VL/TGANsITB\nFoVfagcTFH+UcMwv9U8EPkxy3C/1230ReM26n5f6/ZIgZwN7bY/3Wcf8YBqmOwrrz/h/tBmYzxHn\ntc9Ui2kdvYG/PsM4zG9KrQx2o/ml/n8A7mNwhwLwT+0A/ZigewszzR38U/8c4AjwQ+Bt4F+ASvxT\nv90y4Bnrfl7q90tQ9LtdQJ70k/6zeOVzVgE/A74GHE94zuufoQ/TfTYTWIz57dzOq/VfDxzG9C+n\nmo3o1drjPof55eJa4E5MV6ydl+svBuYB/2T9eZLhvRZerj+uFPiPwLNJnnNcv1+CYj+mvy3uHIam\noZe1Mrga/SzMlwEM/0wzrWNuK8GExFOYrifw32cAM5j3IjAff9R/OXADpvvmGeBKzH8DP9Qed9D6\n8wjwPLAA/9S/z7pttB4/hwmMQ/ij/rhrgd9h/huAf/7+86IYs6qwFpOYXh3MhuFjFN9lsC/wfoYP\nJpVimr27cH9dSwj4N0wXiJ1fPsNkBmd1jAd+g1nl75f6465gcIzCL7VXANXW/UpgHaav3C/1g/n/\n5ZPW/RWY2v1UP8BPgFtsj/1Wf86uxczC+QDvbkH+DHAAs1fVXuA2zCyQNSSfnvYA5vNsB/5kVCtN\nbhGm62YTg9PsrsE/n+EiTP/yJsw0zfus436pP+4KBmc9+aX2OZi/902YqdXxf6N+qR/MJQ42MvSa\nOH6qvxI4ymBgg7/qFxEREREREREREREREREREREREREREZHCmwj8he3xDJJvf5AP12MWcqUyF3ii\nQO8tIiIO1TJ0JX0hvcrgpmypRIGphS9FJHN+2etJpFC+DZyLWYX+HWA2g8FxK2a/q1cwezDdBdyL\nWf29AaixzjsXaMbsmvobzEVkEp2D2S4hvpNng/U+m4C1tvOaredERMQj7MEAQ1sYt2KuAFaJ2Ufq\nY6DReu7vMbvrAvw7cJ51f6H1ONEy4DHb482YTdoAJtiOLwF+mt1HECmsYrcLEHHZSBuhvYrZcvok\n0MngZn1bMGMKlZidX+3jGqVJXmcWg7urgtk070mgCbOvUNxBTFiJeIaCQiS9M7b7fbbHfZh/P+OA\nDsx1GEZiD6W/wGzDfR1mW+j5QLt1jheuayAyQGMUMtYdZ+hum5mKf+kfx4xf/Jnt+Nwk5+9h8LoA\nYMY13gQexFw7YKZ1/CzrXBHPUFDIWNeG6QbaghnMtl8FLPGKYIn344+XA3cwuMX2DUneZx3mQjhx\n38WMU2yxnttsHV+AGRAXEZEx6NcMDmCnEkXTY8Vj1KIQGT1/B/y3NM/PxVxI5nCac0RERERERERE\nRERERERERERERERERETEuf8PJdrWoGSAFtkAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10f681ad0>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 3. Reduce the problem size averaging 10x10 ROIs to single pixels. ##"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data_dce = transpose(data_dce, (1,2,0)).copy('C')\n",
      "nxx = 6\n",
      "nyy = 5\n",
      "data_dce_red = zeros((nxx, nyy, nt))\n",
      "t1_map_red = zeros((nxx, nyy))\n",
      "s0_map_red = zeros((nxx, nyy))\n",
      "for x in range(nxx):\n",
      "    for y in range(nyy):\n",
      "        tmp = data_dce[10*x:10*(x+1), 10*y:10*(y+1), :]\n",
      "        data_dce_red[x,y,:] = mean(mean(tmp, axis=0), axis=0)\n",
      "        tmp = T1map[10*x:10*(x+1), 10*y:10*(y+1)]\n",
      "        t1_map_red[x,y] = mean(tmp)\n",
      "data_dce = data_dce_red\n",
      "T1map = t1_map_red\n",
      "nx = nxx\n",
      "ny = nyy"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 4. Convert signal to tissue concentration $C_t$ ##"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "R1map = 1.0 / T1map"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data_dce = reshape(data_dce, (nx, ny, nt))\n",
      "S0 = data_dce[:,:,:4].mean(axis=2)\n",
      "R1_eff = dcemri.dce_to_r1eff(data_dce, S0, R1map.T, TR, flip_angle)\n",
      "data_dce = reshape(data_dce, (-1, nt))\n",
      "R1_eff_old = dcemri.dce_to_r1eff_old(data_dce, S0map, range(nx*ny), TR, flip_angle)\n",
      "data_dce = reshape(data_dce, (nx, ny, nt))\n",
      "R1_eff_old = reshape(R1_eff_old, (nx, ny, nt))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "converting DCE signal to effective R1\n"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "R1_eff_old.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 34,
       "text": [
        "(6, 5, 661)"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# show map of R1_eff \n",
      "R1_eff = reshape(R1_eff, (nx, ny, nt))\n",
      "if plotting:\n",
      "    figure(7)\n",
      "    clf()\n",
      "    imshow(mosaic(R1_eff.max(axis=2)), interpolation='nearest', cmap='cubehelix3')\n",
      "    title('max $R_1(t) (s^{-1})$')\n",
      "    xlabel('ro')\n",
      "    ylabel('pe')\n",
      "    colorbar()\n",
      "    figure(8)\n",
      "    clf()\n",
      "    imshow(mosaic(R1_eff_old.max(axis=2)), interpolation='nearest', cmap='cubehelix3')\n",
      "    title('max $R_1(t) (s^{-1})$ old method')\n",
      "    xlabel('ro')\n",
      "    ylabel('pe')\n",
      "    colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEdCAYAAAASMeKuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGDtJREFUeJzt3XuUHGWZx/HvzCQzyRCSSQiXkAsxBgRUIBEIoEAALzEL\n6LIuCior6q67R5SzCiiePbpZRY+7CqwLiwhydQX3gAocEGGBhFtIDAkJEaKEXEzIhUBmCJAQMt29\nfzzVMz3VXV3V1VVdVd2/T06d6Z6peuvtzPTT76XqfUBEREREREREREREREREREREREREIvcOj+9P\nALprPMbPGOAs4NIaziVSVXvSFWhR04DjPH62DbikxmP8vAY8DXQGPFct1gGnefzsJuC7dZZfj3V4\n161WN5Hsa0kdBY9kfAm4zeNn/cC9wHk1HBOW17lqUXC2Wn8WtXXAqTGev5GvJROGJV2BjJoJ/CvW\nHbgF6AKOBH4JLPA59khgo88+fwC+4pQd9Jh9KX/zvAI8VOO5otYWU7luhQacq1GvJROy3vJYB1wE\nrABeB34O7A/8DmuqPwj0OPt+E1gN7AD+CHy8pJx3Aq8CM5znB2JN+pM8zrvUOd/lzjn/2znnTwLU\n+XTgYdf3HqI8kG8Dplc5xm0b8CvXVho4qv3hl56rksOA+UAvsBI4w2O/Gdj/zQ7gdmCET53XEfz3\nB/Z7uRN4GViDBT2AW4EpwD1OORe56rQc6HPq1BXwNdX6WiRj1gJPYp+6BwJbsV/4kdgfyUPAt519\nPwEc4Dw+G3ij5DnAF7GgMhL4PfDvPudew+BgYydwB/B3Aer8W4a+kSdSOTCcB3zS45hajQK+DtwN\nvMfnXG7DsaD7TSzAnYK9oQ52fr4Wa/F0AuuBC4EO4G+At4F/q1KvWn5/7di4zb849XgH8CLwYVc9\nSq0DnsJ+z2OB57Dun9drOsQ5LsxrkYxZC5xT8vwO4OqS5xcAv/E4dhlwput7dwHPAs9gf2BeDnOO\n/wjwj9gn5pSAdX6g5PGHgP/Fujufde13BvaGdx9TNBp7fXOBrwU8t5fSc7mdCGx2fe+XwHecx8U3\n7UnAS679nsA/eAT9/c3C3tClLgVucNXDXf65Jc9/CFyD/2sK81paTta7LWCfVkW7XM/fwj51wT5d\nl2HN1F7sE3gfV1nXA+8G/gvYU+Wcp2KB5vfAT7EpzwkB69tR8vhBbNDyx1jTu9QuBmdHOij3cWAy\n9sl9eMBzeyk9l9uBwAbX99Y733fv537Drce/xRT093eQc47eku1SYD+f8reUPN7plDeByq9povM4\n7GtpKc0QPNwq/YKnANcBXwbGYU3Yla59RwFXYgFknrOPl9nAwpLn4xh6DYbXdRVgwaK0rjOw5rjb\nGGB7hWOKfgeMx1pKlY6vRem53DZhQar0/+ogyt9cmxh885XuV+sMhdcbdAPWkhhbso3GxoOo8Txe\nr6k4KL2ZaF5LU2vG4FHJXkAem31oB86nvO//n8Bi4B+w6cufepTVhjVrF5V8773Ym6/Y+vC6rgLs\nk7D4aXo48Lzz+FOu/SZg/XL3MWBN+MuALwDvA072qGtQpedyewr7xL4E68rNxt6wt7v2W4gFua86\n+50FHFNnvUotxgZDL8HGpTqw3+HRzs+3YgPf1RSDxSKqv6a4X0tTaMbgUXA9LmBv0B9jfxRbsD+6\nx0v2+xg28PZPzvOvYdOxpf1xgCOA72N/vGeVfP8G7AKuD+NvAXCs8/hVLNCcQ/kU71FYP9t9DNhs\nw9PYmM25eI9XBFV6Lrc92JjIR7FZmauw8Zk/V9jvLOBz2Os6G5sZqVWl3x9ADnuDH4UNVm8Dfoa1\nPgB+gA2m9uI9BlQsz+81vR3RaxEJ5SAGB+BK9QDf8zl2BDYNXMsxYbnPJRJIM7Y80sKr796HdZ/G\nVzn2U8C1NR4TlvtcIpIgv+sq2oC/9zh2MtaNquWYsLzOJSIiIiIiIpJ1SV8xp4tupOl1tXWzu7Cz\nrvfa2LFjC729vUF378UuXIyVgodLd/so/50Ceju/m872Lv8dA+hqj+amzl25NxnZsVckZY3s2DuS\ncgB29Pcyeli1i3qD64ro9fW+vZWxnfvXXc7anSuh/vdaIVfIBdqxo60jivP50noeIhmRL+STrsIQ\nCh4iGZFHwaNlOM3HVBnWVm2lgeRE1S2L0oiIuj9R6S9Uuj8yOQoeMepoS99/7/B2rzvvk9XVPjLp\nKpQZ2RHd+FcU8ikbIkzfX7eIVKRui4iEogFTEQlF3RYRCUXdFhEJRd0WEQmln2BXmDaKgodIRuQL\n6RrziHslsTnAKuAF4Bsxn0ukqeUD/muUOFseHdjCsh/Elun/A7ay1vPVDhKRylqp5XEstpz/Omy1\n6tvRkncioaWt5RFn8JjI0KxcGylPpCMiAeUoBNqq6MCyJt7j8fOfYEMMyxlM+u4pzm5LutpYIhkX\nQbflQizZd6WFWOYC07EE5rOwnL7HVSsszuDxErY6d9FkBtP5iTStXbk3eCv3ZuTl1jlVOwkLEJdR\nOSnWmcDNzuNFWK6g/RmaO3iIOIPHEiyKTcVyg36S8gxsIk1nZMeoIXfk9vVvi6TcOlseVwAXM5hh\nz63SMMMkqgSPOMc8+oELsEzyzwG/QjMtIqHlKQTaKjgdS1G6jOrLE7p/VjVaxX2R2O+cTUTq5HV5\n+pJHF7Hk0cXVDj0B65bMxdKLjgZuAc4r2cc9zDDJ+Z4nLYDsEuUCyFFK40pbUS6AHKWoFkCOSlQL\nIC/ZuSrQjkd3H1rtfCcDF2GJvkvNxXoKc7GB0itJcMBURCLkMw1bi2JBX3K+XgvchwWO1cCbwPl+\nhSh4iGRERFeYLnA2KE9wfkEtBSl4iGREv27JF5EwtJKYiISSthvjFDxEMkItDxEJRcFDRELJpyt2\nKHiIZMUezbaISBjqtohIKJptEZFQ0tVpUfAQyQy1PFyGt3UmXYUhulN6p2j3sDFJV6HMXp3jk65C\nRaP2fkfSVRjCuau2bgoeIhKKui0iEkq/Wh4iEoa6LSISirotIhKKWh4iEorubRGRUNJ2eXqceVtE\nJEL5QrCtghFYFrhnsBxKP/A4xWwst8tKYL5ffdTyEMmIOqZq3wJOAXZi7/nHgQ84X4t6gKuBj2DZ\n4nyvAFTwEMmIOsc8djpfO4EOYLvr5+cCdzKYT/oVvwLVbRHJiELAzUM71m3ZCjyCdV9KHQyMc362\nBPisX33U8hDJiDqnavPAUcAYLH/0bIaOawwHZgKnAd3AQuAp4AWvAhU8RDLCq9uy+smlvLhwWdBi\nXgPuBY5maPDYgHVVdjnbo8CRJBg8bgD+CsvQ/d6YzyXS1LyCx7TjZzLt+JkDzx+4/Eb3LuOBfqAP\nGAl8CJjn2ucu4CpsPKQLmAVcXq0+cY953AjMifkcIi0hH3CrYALwMDbmsQi4B3gIy1VbzFe7Crgf\nWOHscx3l4yJDxN3yeAyYGvM5RFpCLvx0y7PYeIabO1ftj5wtEI15iGSELk8XkVB0V61LrtA/8LiN\ndtrbdOmJZNuOt7bw+u6tkZerlodLR1viVRCJ1OgRBzB6xAEDzze9viKSclN2R37ssy23AU8Ch2Dz\nyOfHfD6RplXHjXGxiPtj/5yYyxdpGeq2iEgoOQUPEQkjbWMeCh4iGaFui4iEouAhIqGo2yIioajl\nISKh5FJ2fbqCh0hGqNsiIqEoeIhIKBrzEJFQ1PIQkVDS1vLQ4hkiGVEoBNsqCJJu8tPAcmwN0yeA\nI/zqo5aHSEbk8m1hDw2SbnINcBKWmmEO8DPguGqFKniIZESdYx5+6SYXljxeBEzyKzDx4DG6oyfp\nKgwxrsv3/ywRe3VPTLoKZUaMnpJ0FSrqnH5U0lUYauOtkRRTZ/BoB5YC7wSuoXpahS8A9/kVmHjw\nEJFg6gwefukmi04BPg+8369ABQ+RjPAKHlueXsqWpXWnmwQbJL0OG/Po9StIwUMkI/Ie97bsN2Mm\n+80YzOm04vpQ6SanAL8GPgOsDlIfBQ+RjKij2zIBuBkb92gHbmUw3SRY5rhvA2Ox8RCAPcCx1QpV\n8BDJiHz4qdog6Sa/6GyBKXiIZIQuTxeRUBQ8RCSUghYDEpEw1PIQkVAKhdADprFQ8BDJCLU8RCSU\nfC7pGgyl4CGSEWlrecS9GNBk4BHgj8BK4Ksxn0+keRUCbg0Sd8tjD/DP2ApGo4CngQeB52M+r0jT\nabWp2i3OBvAGFjQORMFDpGZp67Y0csxjKjADW6VIRGrUqsFjFHAHcCHWAhGRGhXC3xgXi0YEj+HA\nncAvgN+6f7gzNxhLhrd1Mry9swFVEolP36ur6Nu+KvJyW23Mow34ObZe4pWVdujuGBVzFUQaq2ef\nQ+nZ59CB5+tfvDuSctMWPOKeqn0/tjLRKcAyZ5sT8zlFmlOLTdU+jhJLiUSiVQdMRaROaeu2KHiI\nZEXKWh7qUohkRCEXbKsgyG0i44H7savBVwKf86uPgodIRtSR6Lp4m8i7sfyzXwYOc+1zATahcRSW\nEOrH+PRMFDxEsiIfcCu3BWtRwNDbREptBkY7j0cDr2K5XjxpzEMkIyKabZlK5dtErgMeBjYBewNn\n+xWklodIVoRveRRVu03kW1jr5ECs63I1FkQ8qeUhkhFeLY/tq5bSu2qp3+FVbxMBTgAucx6/CKwF\n3gUs8SpQwUMkKzxaFeMOmcm4QwYTwq296+fuXXxvEwFWAR8EngD2xwLHmmrVUfAQyQiPadggireJ\nrMBmVMC6KVOcx9cC3wduBJZjwxmXANurFRokeByANWcmYvelHA4cj0UyEWmU8COmQW4TeQU4o5ZC\ngwyY3gQ8wODUzgvYnLGINFAhH2xrlCDBYzzwK6DYaNqDz/yviMQgg3fVvoEFkKLjgNfiqY6IeMrg\njXFfB+4CpgFPAvsCn4iqAt3DeqIqKhJjeg713ykBXWMmJV2FMh1T3pV0FSo78Z1J12Co+6MpJou3\n5D+HzQvvAnY4j/8UZ6VEpIJ8uqJHkOBxCxY0LsPmi88FbgX+NsZ6iYhbBtNNvhubni16GGuNiEgD\npa3bEmS2ZSl2XUfRcVjmNxFppHwh2NYgQVoeR2OXrG7AJoKmYGMezzrPj4itdiIyKGUtjyDBQ6ud\ni6RBBqdq18VdCREJIGWDHroxTiQjtHq6iISTU8tDRMJQy0NEQtGYh4iEopaHiITSYi2PEcACoAvo\nxO7OvTTmc4o0pxZrebwFnALsdM71OPAB56uI1CJlLY9G5G3Z6XztBDrwWVRVRDzkCsG2ckFy1RYd\ng60UeJZfdRoRPNqxZDJbsRegO3JFwgh/Y1yQXLVgH+4/xJYvavOrTiOCRx7LQDUJOAlLoisiNSrk\nC4G2CoLkqgX4CpZRbluQ+jRytuU14F7sLt35A9/c88rADl3t3Yzo6G5glUSi17dmGX1rlvnvWKt4\nc9VOBD4GnIp1XXzPFnfwGI/1n/qAkcCHgHmlO4wZPr7CYSLZ1TNtBj3TZgw8X//wjdEUXP9aHdVy\n1V4JfBMLGm0E6LbEHTwmADdj3aN2bPnCh2I+p0hz8pht6du4nL6Ny/2O9stV+z7gdufxeOCj2FjJ\n3V4Fxh08ngVm+u4lIv48bozrmXAEPRMG1+T6y+JfuHcJkqt2WsnjG4F7qBI4QFeYimRGW/huS5Bc\ntTVT8BDJinhz1ZY6P8hOCh4iWZHBvC0ikgYKHiISSsrubVHwEMkKtTxEJJRcuu7JV/AQyQp1W0Qk\nFHVbRCQUtTxEJBS1PEQkFLU8RCQUtTxEJJRcLukaDJF48Ojs2CvpKgwxctqspKtQ2dTJSdeg3OHj\nkq5BRcf8dbo+oRd8K6KC1G0RkVDUbRGRUAq6wlREwlDLQ0RC0ZiHiISSsm5LI5I+iUgUcvlgW7mg\n6SZ/ArwALMdyu1SllodIVuRDtzyK6SafwXK3PA08iGWOK5oLTAcOBmYB12CpKT2p5SGSFYVCsK1c\nkHSTZ2I5lsCyyfUA+1erjloeIlkRzZjHVLzTTW4oeb4Ryy+91asgBQ+RrIg33SSUp5isekIFD5Gs\n8Gh59PX+mb6+P/sd7Zdu8iVsYLVokvM9TwoeIlnhMWDaM2Y6PWOmDzxfv/4+9y5B0k3eDVyA5as9\nDktO79llAQUPkezIh76rNki6yfuwGZfVwJsEyBrXiODRASzBBmDOaMD5RJpT+KnaoOkmL6il0EYE\njwux5tLeDTiXSNMqpOzy9Liv85iENYWup3wkV0RqUcgH2xok7pbHFcDFwOiYzyPS/MJ3W2IRZ8vj\ndOBlbIBGrQ6ReoW/wjQWcbY8TsAueZ0LjMBaH7cA55XutP3tzQOPR3aMYmSHhkYk2/6yaBkbFi/z\n37FW+f7oy6xDnMHjW84GcDJwEa7AATCuc0KMVRBpvCmzZjBl1uBNqQuvuimScgsp67Y08jqPdA0V\ni2RNytbzaFTwWOBsIhJWyqZqdYWpSFa0aMtDROpUUPAQkVBaeMBUROpQaKGpWhGJkrotIhKKgoeI\nhKEBUxEJRwOmIhKGWh4iEo6Ch4iEkbapWmWME8mIQiEfaKvgBmwl9GerFD8bW3tnJTA/SH0UPESy\nIvwyhDcCc6qU3ANcjS1Q/h7gE0Gqo26LSEbUMWD6GJZm0su5WEKojc7zV4IUqpaHSFYUcsG22h0M\njAMewdKkfDbIQYm3PIa1dyVdhaEOnZZ0DSrqOrEn6SqUmX5k6CREsfqPQ05MugpD/CiicmKcqh0O\nzAROA7qBhcBTwAvVDko8eIhIMF7BY8furby+++V6it6AdVV2OdujwJEoeIg0B6+p2r2H78Pew/cZ\neL75jZW1Fn0XcBWW3bELmAVc7neQgodIRtTRbbkNW4R8PNbK+A7WVQHLU7sKuB/LZZsHrsOyPFal\n4CGSEQVCjzGdE2CfH1Hj8IyCh0hW6PJ0EQlDN8aJSCgKHiISioKHiISSz+9JugpDKHiIZEQBtTxE\nJAR1W0QklEK4m95i04jgsQ7YAeSAPcCxDTinSNNpxZZHAVulaHsDziXStFp1zKOtQecRaVqt2vL4\nP6zbci12042I1ChfaL2p2vcDm4F9gQexO/gea8B5RZpKK7Y8NjtftwG/wQZMB4LHK7s3DOzY3TGa\n7mFjGlAlkfjMnz+f+fPnR15uvsWCRze2wMjrwF7Ah4F5pTuM75occxVEGmv27NnMnj174Pm8efO8\nd65Bqw2Y7o+1Norn+h/ggZjPKdKUWu06j7XAUTGfQ6QltFq3RUQiki+kK92kgodIRuRTNuahpE8i\nGVEo5AJtFfjlqv00sBxbAPkJ4Igg9VHwEMmIfMB/Ffjlql0DnIQFje8CPwtSH3VbRDIixly1C0se\nLwImBSlUwUMkIxo05vEF4L4gOyp4iGSEV8tjT/5t9kRz38spwOexW0p8KXiIZES/R4Boa2ujs61z\n4PlbuZ1hij8Cu2l1DtAb5AAFD5GMiPHy9CnAr4HPAKuDHqTgIZIRdVxh6per9tvAWOAa53uBVvxT\n8BDJiDoGTP1y1X7R2WrSNNd57Ox/LekqlOnbtCLpKpTZ/vzSpKtQ0Ut/WJZ0FcrEcVt9PfKFfKCt\nUZoneOR2JF2FMmkMHr0pDR6bFDx8FcgH2hpF3RaRjEjbvS0KHiIZ0Z+yu2qTXtV8PjYKLNLMFmDp\nR+pRqGHfXmBcnecTERERERGRdJiD5YN5AfhGwnUB/8VXkjAZeAT4I7AS+Gqy1QFgBHYL+DPAc8AP\nkq3OEB3AMuCepCsi8enArsefil1y+wxwWJIVAk4EZpCu4HEAg4tRjwL+RPL/T2DpOcBm/p4CPpBg\nXUp9DVvt/+6kK5JWzXCR2LFY8FiHXZN/O/CxJCuELb4S6M7EBtqCBVaAN4DngQOTq86A4i2gndgH\nQRoSok8C5gLXk/yMZGo1Q/CYiN3sU7TR+Z54m4q1jBYlXA+wv8FnsG7eI1j3JWlXABdDyq7KSplm\nCB61zH+LdVnuAC7EWiBJy2PdqUnYOpqzE60NnA68jI13qNVRRTMEj5ewwcCiyVjrQ8oNB+4EfgH8\nNuG6uL0G3AscnXA9TgDOxBKW3QacCtySaI0kNsOAF7GmeCfpGDAFq0+aBkzbsDfBFUlXpMR4oMd5\nPBJ4FDgtueqUORnNtjS9j2KzB6uBSxOuC9gn1iZgNzYec36y1QFsFiOPBddlzlZtOf5GeC+wFKvT\nCmycIU1ORrMtIiIiIiIiIiIiIiIiIiIiIiKSIW3ofg0RCWgqdvXtzdgiQDdgl82vAM5OrloiknZT\ngRy27slZwANY62M/YD22SJBIzZrhrlrxtx5YjK1w9ktsGYOXsZQAxyRYL8kwBY/W8KbztUD5mIfW\nQ5FQFDxay2PAJ7Hf+77Y4juLE62RZJbSTbaGYuviN8DxwHLnexdj3RcRERERERERERERERERERER\nERERkUH/D7N/M9N33amDAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x11360d350>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEdCAYAAAASMeKuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFk9JREFUeJzt3XmUXFWdwPFvdyfdSROSJoQJgSQTIyCiAonsigTcMCPg\ncBwVVGZQZzuinFFxm3Oc4cygxxkFxoHDIAICrnNwAQ6IMEDYDWJCIEKUXcISliRsgZB01fxxq9PV\n1bW8elWv3ntd30/OO11V/d59N5D+9V3evT+QJEmSJEmSJEmSJEmSJEkT3utqfD4HGGzymkZmAMcC\nX2niXsqY3rQroExYCBxU43vPAF9s8ppGngd+B/RHvJekjPpmg+/vD5zQ5DWN/DnwLxHvpQyalHYF\n1BaLgX8ldAcuBgaAfYAfATc2uHYfYG2Dc34LfKZUdtRrdgKOqPjsWeC6Ju8ldYVHgC8AdwMvAucD\ns4FfEZrq1wJDpXO/DDwAvAD8HvhAWTmvB54DFpXe70Jo0r+jzr1/CBxT9v4DwKoIdf5n4E0Vn13H\n+F8s3wZ2q3NNsxZQveVReS+pKzwM3Eb4rbsLsA5YQfhNPUD4ofxa6dwPAjuXXn8IeKnsPcCnCEFl\nKvBr4D8a3PshRgcb+4FLgb+OUOdfAj1l73cFrq9y3gnAh2tc06xpwOeBy4E3N7iX1BUeBo4re38p\ncHbZ+5OAX9S4diVwdMVnlwH3AHcBk+vc942l698L/AOhxTM/Yp2vKXv9buB/Cd2dj1ecdxThB77y\nmhHTCX+/pcDnIt67lvJ7KaOcbWm/dWWvX6l4/yrhty6E364rgQ2l483AjhVlfY/QPfhvYEudex5B\nCDS/Bv6HMOU5J2J9+8peXwtsJXQbLqk47xVGZ0f6GO8DwDxCy2uviPeupfxeyiiDR/KqNe/nA+cB\nnwZmAjsAqyvOnQacSQggp5bOqWUJcHvZ+5mMfQaj1nMVEIJFeV0XEaZRK80A1le5ZsSvgFmEllK1\n65tRfi9llMEjHdsBBcLsQy9wIuP7/v8F3AH8HXAloUVRTQ9hIHV52WdvIfzwjbQ+aj1XAfAUo62h\nvYD7Sq8/UnHeHMIAb+U1AAcCpwGfBN4KHFajrlGV30sZZfBIXrHidZHwA/ptQmvhKULguKXsvGOA\n9wD/WHr/OcJ0bPl4CsDewNcJg6rHln1+AeEBrvdEqN+NwAGl188RAs1xjJ/i3Re4tco1AE8TgtPR\nwPG0Pl5Rfi9JKav1UNYQ8O8Nrp0CnN7kNXFV3ksZZcuje9SaWt1I6D7NqnPtR4Bzm7wmrsp7SUpR\no+cqeoC/rXHtPMY+fBblmrhq3UuSJEmSpGa1sj6hHYqNT5HybaBnkM3FTS39rO2www7FDRs2RD19\nA+FBwUQZPCoM9k5rfFJErxU209870JayBnqntKWcV4ZfZmrfdm0pa2rf9m0pB+CFrRuYPqneQ7TR\nDbTp77fhtXXs0D+75XIe3rQaWv9ZKw4XhyOd2NfT1477NeR+HlJOFIqFtKswhsFDyokCBo+uUWo+\nZsqknnor+9PTrm5ZO01pU/enXbYWq61HTI/BI0F9Pdn7zzu5N5sr3Qd6p6ZdhXGm9rVv/KsdChkb\nIszev25JVdltkRSLA6aSYrHbIikWuy2SYrHbIimWrUR7wrRTDB5SThSK2RrzSHonsSOBNcD9wJcS\nvpc0oRUi/umUJFsefcBZwLuAxwk5SC9ndHduSU3oppbHAYTt8x8hJCz6CW4xJ8WWtZZHksFjV+Cx\nsvdrS59JimGYYqSjjj5ClsIranz/O4QhhlWMJlmvKcluS7baWFLOtaHbcjJwL1BtI5alwG7A7oQk\nXucQcv/UlGTweJywG/aIeYTWhzShvTL8Eq8Ov9z2clucqp1LCBCnUT0R+dHARaXXywm5eWYzNtfy\nGEkGjzsJUWwB8ATwYcZnPJMmnKl908asyN249Zm2lNtiy+MM4BRgeo3vVxtmmEud4JHkmMdW4CRC\n5vZ7gZ/iTIsUW4FipKOK9xNSgq6k/vaEld+rG62SfkjsV6VDUotqPZ5+503LufOmO+pdegihW7KU\nkM5zOnAxcELZOZXDDHNLn9XkBsgV2rkBcjtlcaetdm6A3E7t2gC5Xdq1AfKdm9ZEOnG/wT3r3e8w\n4AvAURWfLyX0FJYSBkrPJMUBU0lt1GAathkjBf196eu5wFWEwPEA8DJwYqNCDB5STrTpCdMbSweM\nTyh+UjMFGTyknNjqknxJcbiTmKRYsrYwzuAh5YQtD0mxGDwkxVLIVuwweEh5scXZFklx2G2RFIuz\nLZJiyVanxeAh5YYtjwqTe/rTrsIYgxldKTo4aUbaVRhnu/5ZaVehqmnbvy7tKoxRWlXbMoOHpFjs\ntkiKZastD0lx2G2RFIvdFkmx2PKQFItrWyTFkrXH05PM2yKpjQrFaEcVUwhZ4O4i5FD6Ro1bLCHk\ndlkNLGtUH1seUk60MFX7KnA4sInwM38L8PbS1xFDwNnAewnZ4ho+AWjwkHKixTGPTaWv/UAfsL7i\n+8cDP2M0n/SzjQq02yLlRDHiUUMvoduyDriB0H0ptzsws/S9O4GPN6qPLQ8pJ1qcqi0A+wIzCPmj\nlzB2XGMysBh4JzAI3A78Bri/VoEGDyknanVbHrhtBQ/evjJqMc8DVwL7MTZ4PEboqrxSOm4C9iHF\n4HEB8BeEDN1vSfhe0oRWK3gsPHgxCw9evO39NadfWHnKLGArsBGYCrwbOLXinMuAswjjIQPAgcDp\n9eqT9JjHhcCRCd9D6gqFiEcVc4DrCWMey4ErgOsIuWpH8tWuAa4G7i6dcx7jx0XGSLrlcTOwIOF7\nSF1hOP50yz2E8YxKlblqv1U6InHMQ8oJH0+XFIuraisMF7due91DL709PnqifHvh1ad4cfO6tpdr\ny6NCX0/qVZDaavqUnZk+Zedt75948e62lJuxFfmJz7b8GLgN2IMwj3xiwveTJqwWFsYlIulf+8cl\nXL7UNey2SIpl2OAhKY6sjXkYPKScsNsiKRaDh6RY7LZIisWWh6RYhjP2fLrBQ8oJuy2SYjF4SIrF\nMQ9JsdjykBSLLQ9JsWSt5eHOO1JODBd6Ih1VRMlV+1FgFWED5FuBvRvVx5aHlBMttDyi5Kp9CHgH\nIa/LkcB3gYPqFZp68JjeN5R2FcaYOTA37SpUtd3grmlXYZwp0+enXYWq+nfbN+0qjLX2krYU02K3\npVGu2tvLXi8HGv4g2G2RcqJYjHbU0ChXbblPAlc1qk/qLQ9J0bTY8miUq3bE4cAngLc1KtDgIeVE\nocbalqdWrGDdipZz1UIYJD2PMOaxoVFBBg8pJ2q1PGYvWszsRaMJ4e45P1au2vnAz4GPAQ9EqY/B\nQ8qJQvVp2CjmABcRxj16gUsYzVULIe3k14AdgHNKn20BDqhXqMFDyokWxjyi5Kr9VOmIzOAh5UTW\nnjA1eEg5UXQzIElx2PKQFEuxGHvANBEGDyknbHlIiqUwnHYNxjJ4SDmRtZZH0gvj5hEW4fweWA18\nNuH7SRNXMeLRIUm3PLYA/0RYzTcN+B1wLXBfwveVJpxum6p9qnQAvEQIGrtg8JCalrVuSyfHPBYA\niwgbjUhqUrcGj2nApcDJhBaIpCYV4y+MS0Qngsdk4GfAD4BfVn5z0/BoLJnc08/k3v4OVElKzsbn\n1rBx/Zq2l9ttYx49wPmELc/OrHbCYN+0hKsgddbQjnsytOOe294/+uDlbSk3a8Ej6anatxE2Fzkc\nWFk6jkz4ntLE1GVTtbfgJstSW3TrgKmkFmWt22LwkPIiYy0PuxRSThSHox1VRFkmMgu4mvA0+Grg\nbxrVx+Ah5UQLSZ9Glom8iZBC8tPAGyvOOYkwobEvIafLt2nQMzF4SHlRiHiM9xShRQFjl4mUexKY\nXno9HXiOkK6hJsc8pJxo02zLAqovEzkPuB54Atge+FCjgmx5SHkRv+Uxot4yka8SWie7ELouZxOC\nSE22PKScqNXyWL9mBRvWrGh0ed1lIsAhwGml1w8CDwNvAO6sVaDBQ8qLGq2KmXssZuYeozmdHr7s\n/MpTGi4TAdYA7wJuBWYTAsdD9apj8JByosY0bBQjy0TuJsyoQOimzC+9Phf4OnAhsIownPFFYH29\nQqMEj50JzZldCetS9gIOJkQySZ0Sf8Q0yjKRZ4Gjmik0yoDp94FrGJ3auZ8wZyypg4qFaEenRAke\ns4CfAiONpi00mP+VlIAcrqp9iRBARhwEPJ9MdSTVlMOFcZ8HLgMWArcBOwEfbFcFBicNtauotpgx\ntGfjk1IwMGNu2lUYp2/+G9KuQnWHvj7tGox1dXuKyeOS/HsJ88KvAC+UXv8hyUpJqqKQregRJXhc\nTAgapxHmi48HLgH+KsF6SaqUw3STbyJMz464ntAakdRBWeu2RJltWUF4rmPEQYTMb5I6qVCMdnRI\nlJbHfoRHVh8jTATNJ4x53FN6v3ditZM0KmMtjyjBw93OpSzI4VTtI0lXQlIEGRv0cGGclBPuni4p\nnmFbHpLisOUhKRbHPCTFYstDUixd1vKYAtwIDAD9hNW5X0n4ntLE1GUtj1eBw4FNpXvdAry99FVS\nMzLW8uhE3pZNpa/9QB8NNlWVVMNwMdoxXpRctSP2J+wUeGyj6nQiePQSksmsI/wFXJErxRF/YVyU\nXLUQfrl/k7B9UU+j6nQieBQIGajmAu8gJNGV1KRioRjpqCJKrlqAzxAyyj0TpT6dnG15HriSsEp3\n2bYPtzy77YSB3kGm9A12sEpS+218aCUbH1rZ+MRmJZurdlfgGOAIQtel4d2SDh6zCP2njcBU4N3A\nqeUnzJg8q8plUn4NLVzE0MJF294/ev2F7Sm49b066uWqPRP4MiFo9BCh25J08JgDXEToHvUSti+8\nLuF7ShNTjdmWjWtXsXHtqkZXN8pV+1bgJ6XXs4D3EcZKLq9VYNLB4x5gccOzJDVWY2Hc0Jy9GZoz\nuifXn+74QeUpUXLVLix7fSFwBXUCB/iEqZQbPfG7LVFy1TbN4CHlRbK5asudGOUkg4eUFznM2yIp\nCwwekmLJ2NoWg4eUF7Y8JMUynK01+QYPKS/stkiKxW6LpFhseUiKxZaHpFhseUiKxZaHpFiGh9Ou\nwRipB4/+vu3SrsIYUxcemHYVqlswL+0ajLfXzLRrUNX+f5mt39A3frVNBdltkRSL3RZJsRR9wlRS\nHLY8JMXimIekWDLWbelE0idJ7TBciHaMFzXd5HeA+4FVhNwuddnykPKiELvlMZJu8i5C7pbfAdcS\nMseNWArsBuwOHAicQ0hNWZMtDykvisVox3hR0k0eTcixBCGb3BAwu151bHlIedGeMY8F1E43+VjZ\n+7WE/NLrahVk8JDyItl0kzA+xWTdGxo8pLyo0fLYuOGPbNz4x0ZXN0o3+ThhYHXE3NJnNRk8pLyo\nMWA6NGM3hmbstu39o49eVXlKlHSTlwMnEfLVHkRITl+zywIGDyk/CrFX1UZJN3kVYcblAeBlImSN\n60Tw6APuJAzAHNWB+0kTU/yp2qjpJk9qptBOBI+TCc2l7TtwL2nCKmbs8fSkn/OYS2gKfY/xI7mS\nmlEsRDs6JOmWxxnAKcD0hO8jTXzxuy2JSLLl8X7gacIAja0OqVXxnzBNRJItj0MIj7wuBaYQWh8X\nAyeUn7T+tSe3vZ7aN42pfQ6NKN/+tHwlj92xsvGJzSpsbX+ZLUgyeHy1dAAcBnyBisABMLN/ToJV\nkDpv/oGLmH/g6KLU28/6flvKLWas29LJ5zyyNVQs5U3G9vPoVPC4sXRIiitjU7U+YSrlRZe2PCS1\nqGjwkBRLFw+YSmpBsYumaiW1k90WSbEYPCTF4YCppHgcMJUUhy0PSfEYPCTFkbWpWjPGSTlRLBYi\nHVVcQNgJ/Z46xS8h7L2zGlgWpT4GDykv4m9DeCFwZJ2Sh4CzCRuUvxn4YJTq2G2RcqKFAdObCWkm\nazmekBBqben9s1EKteUh5UVxONrRvN2BmcANhDQpH49yUeotj0m9A2lXYaw9F6Zdg6oGDh1Kuwrj\n7LZP7CREifrPPQ5NuwpjfKtN5SQ4VTsZWAy8ExgEbgd+A9xf76LUg4ekaGoFjxc2r+PFzU+3UvRj\nhK7KK6XjJmAfDB7SxFBrqnb7yTuy/eQdt71/8qXVzRZ9GXAWIbvjAHAgcHqjiwweUk600G35MWET\n8lmEVsa/ELoqEPLUrgGuJuSyLQDnEbI81mXwkHKiSOwxpuMinPMtmhyeMXhIeeHj6ZLicGGcpFgM\nHpJiMXhIiqVQ2JJ2FcYweEg5UcSWh6QY7LZIiqUYb9FbYjoRPB4BXgCGgS3AAR24pzThdGPLo0jY\npWh9B+4lTVjdOubR06H7SBNWt7Y8/o/QbTmXsOhGUpMKxe6bqn0b8CSwE3AtYQXfzR24rzShdGPL\n48nS12eAXxAGTLcFj2c3P7btxMG+6QxOmtGBKknJWbZsGcuWLWt7uYUuCx6DhA1GXgS2A94DnFp+\nwqyBeQlXQeqsJUuWsGTJkm3vTz311NonN6HbBkxnE1obI/f6IXBNwveUJqRue87jYWDfhO8hdYVu\n67ZIapNCMVvpJg0eUk4UMjbmYdInKSeKxeFIRxWNctV+FFhF2AD5VmDvKPUxeEg5UYj4p4pGuWof\nAt5BCBr/Bnw3Sn3stkg5kWCu2tvLXi8H5kYp1OAh5USHxjw+CVwV5USDh5QTHXg8/XDgE4QlJQ0Z\nPKSc2FpjYdzW4ha2tj6Nuzdh0eqRwIYoFxg8pJyo9Xh6X08ffT19296/Nry52aLnAz8HPgY8EPUi\ng4eUEy08YdooV+3XgB2Ac0qfRdrxz+Ah5UQLA6aNctV+qnQ0ZcI857Fp6/NpV2GcjU/cnXYVxll/\n34q0q1DV479dmXYVxkliWX0rCsVCpKNTJk7wGH4h7SqMk8XgsSGjweMJg0dDRQqRjk6x2yLlRNbW\nthg8pJxow3RsW6W9q/kywiiwNJHdSEg/0opiE+duAGa2eD9JkiRJUjYcScgHcz/wpZTrAo03X0nD\nPOAG4PfAauCz6VYHgCmEJeB3AfcC30i3OmP0ASuBK9KuiJLTR3gefwHhkdu7gDemWSHgUGAR2Qoe\nOzO6GfU04A+k/98JQnoOCDN/vwHenmJdyn2OsNv/5WlXJKsmwkNiBxCCxyOEZ/J/AhyTZoUIm69E\nWpnYQU8RAivAS8B9wC7pVWebTaWv/YRfBFlIiD4XWAp8j/RnJDNrIgSPXQmLfUasLX2m2hYQWkbL\nU64HhH+DdxG6eTcQui9pOwM4BTL2VFbGTITg0cz8t0KX5VLgZEILJG0FQndqLmEfzSWp1gbeDzxN\nGO+w1VHHRAgejxMGA0fMI7Q+NN5k4GfAD4BfplyXSs8DVwL7pVyPQ4CjCQnLfgwcAVycao2UmEnA\ng4SmeD/ZGDCFUJ8sDZj2EH4Izki7ImVmAUOl11OBm4B3pledcQ7D2ZYJ732E2YMHgK+kXBcIv7Ge\nADYTxmNOTLc6QJjFKBCC68rSUW87/k54C7CCUKe7CeMMWXIYzrZIkiRJkiRJkiRJkiRJypEeXK8h\nKaIFhKdvLyJsAnQB4bH5u4EPpVctSVm3ABgm7HtyLHANofXxZ8CjhE2CpKZNhFW1auxR4A7CDmc/\nImxj8DQhJcD+KdZLOWbw6A4vl74WGT/m4X4oisXg0V1uBj5M+P++E2HznTtSrZFyy3ST3WGkdfEL\n4GBgVemzUwjdF0mSJEmSJEmSJEmSJEmSpFH/D+1Lnq7gAKYKAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x113627990>"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# convert effecitve R1 to tissue concentration Ct\n",
      "Ct = dcemri.r1eff_to_conc(transpose(R1_eff, (2,0,1)), R1map, Rel)\n",
      "Ct = transpose(Ct, (1,2,0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "converting effective R1 to tracer tissue concentration\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# show map of Ct\n",
      "Ct = reshape(Ct, (nx, ny, nt))\n",
      "if plotting:\n",
      "    figure(9)\n",
      "    clf()\n",
      "    imshow(Ct.max(axis=2), interpolation='nearest', cmap='cubehelix3')\n",
      "    title('max $C_t$ [mM Gd-DTPA]')\n",
      "    colorbar()\n",
      "    savefig('qiba/qiba_v6_Ct.pdf')\n",
      "Ct = reshape(Ct, (-1, nt))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAEMCAYAAAAmrhhaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGQ5JREFUeJzt3XuUHGWZx/HvzCRzI2SGJEqMIYSAMbjLJVkQDhAJ6Coq\nYnblLAFclPXK2ZWIAiuucg67e9YVlwOuF1YJRLmIZwWMoKCCEggCIZhJQgjRBcIlhCSQZHIPZKZr\n/3iqZ2qquqrfnupUdXd+H06d9KUu7wxTT73vW/W+D4iIiIiIiIiIiIiIiIiIiIjI/qQA7AD+Le+C\nOPoRtV3W3wO7gUV5F6RRNeddgP3I0cDXU2z/AvAGMDb0eQ8WeCYlbDsHWIwFpw3A48BFCet7/pJU\nll3ANmAL8Afgc0CT//0OYLu/FPx1i+/PDe1jO7AemA8cEDrOQmAz0Br6/HTg8wnlk5QUGNx9BLuK\nXgR8EjgHuA5oz+j4HvA8gycWwFFAB8kn8Zexcn4TONhfPg+cTPSEC2pK+M4DzgRGYwHpP4F/Bm70\nvx8FHOgvL/rrFt/fHtrHgcAM4Djga4FjTAbeDWwEzqqwfJJSvQWGF4BLgRXYleZG7A/9PmArcD/Q\n7a/7FeBZ7Kr2NDA7sJ/DgU3AdP/9BOA14D0ljtkMzAPeiV3xr8eq2m9gJ+aeDH6WoluBCwLvPwHc\nTPxJ0gVchQWzu4Cd/ufLgI8Db/rvpwNLsd/VT6ks2G0H7sEC5SeAv6hg26J1wK+Bvwx8dgHwAHCL\nv1+RWGuAR4G3YCfzBuwP+higDfgdcKW/7tnAeP/132HV2/GBfX0aCxgdwG+Aq2OOeRUWGMIOBi53\nLHcBmJLiZymu/15gNTANaAFexq7YcU2JM4C9JF8AWrGr+lx/nx/DAsa/JmyzBqvOh71ItIoft27x\n5wE4BFiJ/a6LngXOB97hl+etoe0/ifoY9pkReRdgGL6DXd3B/jA2AMv99z9n8I/tjsA2/wtcgVVN\n7/Y/m4c1D54A+oF/KXGsscBlWG0hbCPw3WH9BINcf5agW7Cr6cPAKuCVhP2PA17HAkfRo8CRWPD5\nAFbbGAF82//+TmBJhT9H0TrgIMd1m4AFQB9WQ/ol8B/+d6cAb8f+X23Hfs7zsCaRZKDemhJgJ0/R\n7tD7PVj7Fuzk6cE6x7Zg1dRwx908rOr7HezKGjYTeAm7Mod5WOdZGq4/S/CYt2BX0nLNCLDm0jiG\n/n8+CTt5N/nbTiAaXF4MvD6fwY7DXyUcC+xk3lxmnSIP+KhflsnAP2HNM7Cf7bf+MQF+hpoTmarH\nwBBW6sSYBNwA/CMwBvvjWxladxR2BZqHVWFLXen6if9DPz/w+hKGNlOGy6VD7SWsE/KDWL9Bksew\nk212wjrrsBM66NDA69sY7Dj8cMJ+jvf380iZMpXTgTX9Tgde9ZcvY02so1PuWxw1QmAo5QCs+vw6\n9jNeyNCOLbCq8xPAZ7Er4f+U2M/vsCtu8KRvBj6DdZYVTcVuuWXlU9iJs7vMer1Y0Ps+1ndwIFb+\nYxm8NfgYVp2/GBgJ/C12kpdTDGKjsbsLt2O1maddf4gYs/3yHIkFg2P814sY2vGadCdGUqrHPoYw\nL/TaA54BrsH+6AtYlTt4Jfso8H7srgLAl7Ce+nMZvJ0G1lQ4Czu5/oRVvwvYlXqTv85MrCp8IvZ8\nQByX2kCpn6WU5xO2C/sW1lS4HPs97PS3vxz7/ezFgsENwL8D92L9DOXcg53ABSwYXEPp4FqpC4Cb\ngLWhz7+LBfPL/WPqdqXUtEOxe/hJdjN49Zb07sdurd6fd0FE4pyH9aInPXkoUlcatY8hS7uAt6Hf\npYiIiIiIiEhanhYtjb60NXV6pHTQQQdVcsxSD9adgY2V+T9K3wm7FHvatwd4Crud3I2NRXkQu628\nEntmJVE17gWn/oVVW2dz+Eni4Xuz8AatzW1V2Vdbc3VGaO/u30lHS3jqguHpaDmwKvsB2Na3hdEj\nXIdKJGur0s+35c0NHNR6cOr9rNm1EtKfL16/1++0YktTS/h4LdizNO/DnktZgj1380zMLs4Evuiv\nP95flmFP/P4Re5Asblv1pItkqeAVnJYS3o2NOH0Beyjtp9iDenHOY/BhvfVYUAAbZfwMNkYmlgKD\nSIYKjv+V8HaGDuZbS3SMS1EnNnK21BOsk7H5NxYnlbMRHonep/wqXU0Z0TQy7yKUVK2mUjW1V6lJ\nUi19Xt9wN62kyf4RbAhAb+jzUdh0BHOxmkMsBYYyWppq71c0sjlpRrb8tDV35F2EiI6W6vU3VUMh\n5vxetPBhFj2UOO/MK1gnYtEhRMeTFM1h6JgfsAFyd2KzgC0oV051PmaoFq+o1ex8rKZqdT5WS7U6\nH7f2b3NasatldPh4I7DOx/diQ+WfoHTnYxc2SG4ig6Nvm4AfYwP/LnE5fu1dDkUaWEzHoos+bDKb\n32B3KG7EgsLn/O9/4P87218nOCT/ZGyOzxXYrUywGc2CUwcMoRpDhlRjcNeoNYZNfVucVhxrt31z\nG1quGoNIhmLuONQcBQaRDKVoSmRKgUEkQ324PfmYNwUGkQwVvJrrkitJgUEkQ/XSx+DySHS5EV0i\n4qjgeU5L3srVGFqw2XmDI7ruJmFUlojEa5QaQ6UjukQkQT+e05K3cjWGUiO6Tth3xRFpbLXQTHBR\nLjDUx08hUmW7+3ewp39n1ffbKLcrKxnRJdIwOlpGDRmZ2dv3WsLa7uqlxlCuj+FJ4B3Y5A6twDkM\nppEXkQoV8JyWvJWrMcSN6BKRYWikR6Lv8xcRSakWagMu9OSjSIZq4VakCwUGkQw1SuejiFRRn1dw\nWmK4DE+Yhc3StBJYGPquxf/unnLlVI1BJEMp+hhchid0A9/Dpo5fC4wL7WMusAooO22XagwiGUox\niMpleMJ52EzQxWeNXg98NxH4EDAPhynjFBhEMpTiOQaXhDPvAMZgeSqfBP4+8N21wGXgNopLTQmR\nDKVoSrhsOBKYgU0x3wk8BjwOvBPYiPUvzHI5mAKDSIYKMaf38kVLWLHoyaRNXYYnvIw1H3b7y8PA\nMViwOAtrSrQDo4GbgQviDqbp4zOk6ePdNer08b/cutRpxTO7ZoSP55JwZhrWQfkBoA3LT3kO1uFY\ndCpwKZbGLpZqDCIZStGUcEk4sxpLIrMC60u4gaFBoahsIVRjyJBqDO4atcawoDexuTBgdvdx1Tje\nsKnGIJKh+hhCpcAgkql6eSS6KoFhZFNtpWXvrNHqceeIrryLEHFAa/jhuNow6sDD8i7CEH5TIrX9\nKjCIiBs1JUQkok81BhEJU1NCRCLUlBCRCNUYRCQibqxErVFgEMmQJoMVkQjVGEQkQrcrRSRCNQYR\niaiTuKDAIJIl3a4UkYh6aUpolmiRDBU8tyVGuYQzs4Ct2KSvPcDXAt91A3dgsz6tAk5MKqdLjeEm\n4MPYLLNHOawvIjFSPBLtknAG4CFs4tewbwP3Amdj533iFFkuNYb5WKQSkZT6C57TUoJLwhkoPR1c\nFzATu8iDzR+5NamcLoFhEbDFYT0RKSNFU8Il4YwHnAQsx2oH7/I/Pwx4DbvIL8Umie1MKqc6H0Uy\nFNeUWPNYD2se70na1KXbcimWb2IX8EFgATAVO89nYLNMLwGuA74CXBm3o6oEhn6vb+B1E800N6lP\nU+rbtj3r2f7GhqrvN65j8dATp3PoidMH3i+8bn54FZeEM9sDr+8Dvo+lrFvrL0v87+7AAkOsqgSG\nliZVPKSxjG4fz+j28QPv121fUZX9pniM4UksN+VkLOHMOVjCmaCDsZsEHtYn0QRs9r97Gas9/Bnr\nwHw66WA6o0UylOI5BpeEM2cDF/nr7gLmBLb/AnAb0Ao8B1yYdDCXwHA7ltZqLBZ1rsQ6MUSkQikf\ncLrPX4J+EHj9PX8pZTlwvOuBXAJDuLoiIsPUXydPPqopIZKhOhkqocAgkqV6GSuhwCCSIQUGEYlQ\nU0JEIlRjEJGI/jrJOKPAIJIhNSVEJEKBQUQi1McgIhGqMYhIhGoMIhKhGoOIRPQXSk3JWHsUGEQy\ntF/VGEa3dFdjN1Uzpm1i3kUo6YDO8Nyd+WsfPSnvIpTUesSxeRdhqLW3VGU39RIYNDmjSIY8z22J\nUS7hTNHx2CxOHwt8dgU2ndtTwE+AtqRyKjCIZChFYCgmnDkDmxb+XODImPW+Cfw68Nlk4DPYTNFH\n+evMiWwZoMAgkqFCwW0pwTXhzBewWaBfC3y2zd+mE+s+6MRmnY6lwCCSoRQ1BpeEM2/HgsX1xcP5\n/24GrgFewmaY7gUeSCqn7kqIZKgQc7tyw9KlbOxJnXCmmEjGw6aOLx7scOCLWJNiK/Az4Hxs1uiS\nFBhEMhTXsfjW6TN46/QZA+9Xzr8pvIpLwpm/wpoYAOOwbFR9WEfjo8Am/7u7sFR2sYFBTQmRDKVo\nSgQTzrRiCWfuDq0zBctTeRjWz3AR8AvgT1ja+w6sFvE+YFVSOVVjEMmQN/yJWlwSzsRZDtyMBZcC\nluPyh0kHU2AQyVDKB5zKJZwJCmeautpfnCgwiGTI8zRWQkRC6uWRaAUGkQwV+vMugRsFBpEM1UuN\nweV25SHAg9gAjJXAxfu0RCKNzHNccuZSY9gLXAIsA0YBfwTux26ViEgFUtyuzJRLYFjvLwA7sIAw\nAQUGkYrVS1Oi0j6GycB0YHH1iyLS+BoxMIzCHrOci9UcRKRCXoPN+TgSuBO4FVgQ/nJX/2CcGNnU\nysjm1qoUTiQvvZtW07t5ddX320h9DE3Yc9mrsGGdEZ0to6pZJpHcdY+dRvfYaQPvX3wuPF5peOol\nMLjcrjwZ+DhwGtDjL2fsy0KJNKwGul35CBqeLVIVjdj5KCIp1UtTQoFBJEuqMYhImFcng6jUdyCS\noRwTzrhuCygwiGSr4LhEpUk447rtAAUGkQylqDGkSTjjuu0ABQaRLA2/xpAm4YzLtkOo81EkQ3H9\nB5tXL2XL6qWJmzrsPi7hTMX3QhQYRLIU8xzDmKkzGDN1MOHMml/cGF5luAln9jpuO4QCg0iGUtyu\nDCacWYclnDk3tM6UwOv5wD1YUpoRDtsOocAgkqXhPxOdJuFM3LaxFBhEMpTykeg0CWdKbRtLgUEk\nS3okWkQi9qdBVJ0juquxm6rp6p5WfqUctHVNzLsIES2T3pl3EUqbeXjeJRjq1+VXcaFh1yISVaiP\nyKDAIJKlOhldqcAgkiE1JUQkSk0JEYmoj7igwCCSqf3pdqWIOKqTTgYFBpEMaZZoEYnqV41BRMJU\nYxCRCPUxiEiEagwiElEnNQaXWaLbgcXAMmAV8I19WiKRRjb8WaKhfNKYjwLLsYz0fwRO9z8/BHgQ\neBpYCVxcrpguNYY9wGnALn/9R4BT/H9FpBLDrzEUk8a8D5vcdQk2n2NwirYHgF/4r48Cfg4cgU0I\newl2cR+FBY37SZjezTWvxC7/31a/gJsdtxORoH7PbYlySRqzM/B6FPC6/3o9FhQAdmABYUJSMV0D\nQ7O/4w1YlWSV43YiElTw3JYo16Qxs7ET/z5KNxkmA9Ox7oFYrp2PBeBYoAubaXYWsNBxWxHxeTGj\nK3tfXkbvy8sTN3U8xAJ/mQncAgSn6BqFpa+bi9UcYlV6V2Ir8CvgOAKBYeve1wdWaGvupL2ls8Ld\nitSW3ud76H2+p/o7jjm9uyceS/fEYwfev/T4LeFVKk0aswg7v8cCm4CRwJ3ArVjgSOQSGMZh89L3\nAh3AXwNXBVfoGjnOYTci9aN7ynS6p0wfeP/i7+dXZ8fDn4/BJeHM4cDzWPgpprXahKWquxHrArjO\n5WAugeFtwI+xfoZmrHryO5edi0jIvk048zHgAqxzcgcwx//uZODjwArsVibAFSRMcesSGJ5iMPqI\nSBrpBlGVSzhztb+EPUKFme315KNIhpo0tZuIRNTJI9EKDCJZUo1BRCIUGEQkQk0JEYlQjUFEIvrr\nY6YWBQaRLKkpISIRakqISIRqDCISoRqDiESoxiAiEaoxiEhEf3/eJXBSlcDQ2nJANXZTNR1TTsi7\nCKVNPqT8Oll715i8S1DS8X9TW1fWh75apR3VSVOiojHaIpLS8CeDhfJ5Jc7H8kqsAP4AHB36vgWb\nqOWecsVUU0IkS96wn3x0ySvxPPAebG7WM4AfAicGvp+LTe92YLmDqcYgkqXh1xhc8ko8hgUFsOnh\nJwa+mwh8CJiHzQGZSIFBJEue57ZEueaVKPoUcG/g/bXAZTim1VVTQiRLw29KVNJreRrwD9gksABn\nAhux/oVZLjtQYBDJUszoyt7Nq+ndvDppS9e8EkcDN2B9DFv8z04CzsKaEu3AaOBmbEbpkhQYRLJU\nKB0Yurun0t09deD9i8/dHV7FJa/EJOAubKr4ZwOff9VfAE4FLiUhKIACg0i29m1eiSuBg4Dr/c/2\nYp2WkVKUO5gCg0iWht/HAOXzSnzaX5I85C+JFBhEsqSxEiISka7GkBkFBpEsxXQ+1hoFBpEsFfaj\n0ZUi4qhOagyuj0Q7j8oSkXie5zkteXOtMTiPyhKRBHXS+ehSY6hoVJaIJCgU3JacudQYiqOyRu/j\nsog0vhpoJrgoFxicRmVtfvPVgdcdLaPoaFGLQ+rbS4t7ePmJnurvuNBX/X3uA+UCg9OorDGtb9sn\nhRPJy6QTpjPphOkD7x/77o+qsl+vBpoJLsr1MXwVG955GDAH+D1lRmWJSAKv4LbkrNLnGOqjgSRS\nqxqkjyHIaVSWiCSogdqACz35KJIhT4FBRCIapPNRRKrIK/Q5LTHKJZyZhk0hvwf4cui7buAObNan\nVQzNNxGhGoNIlvZtwplNwBeA2SW2/zY2nfzZ2HmfmFdSNQaRLA3/dqVLwpnXsElj94Y+7wJmAjf5\n7/sYTExTkgKDSIY8r+C0lFBpwpmgw7CgMR9Yik0v35m0gQKDSJaGP4gqzQMQI4AZwPf9f3cCXym3\ngYhkJO525dbda9m2+5WkTV0TzpSy1l+W+O/vQIFBpIbEBIau9gl0tU8YeL92y5LwKi4JZ4rC0yOs\nx5ohU4E/Yx2YTycVU4FBJEMJtyLLcUk4Mx6rFYzGktfOBd4F7MDuVtwGtALPARcmHUyBQSRDKZ98\nLJdwZj1DmxtBy4HjXQ+kwCCSJT0SLSJhGishIlHefpRXYkRzWzV2Uz3TpuRdgpLaZnbnXYSII46p\nzT/Ub02dmXcRhvivKu1HNQYRiVBgEJGIFLcrM6XAIJIh1RhEJMKjNvt0whQYRLKkGoOIhKkpISIR\nCgwiEqHAICIRhUJ41rXapMAgkiEP1RhEJERNCRGJ8OpkEJUmgxXJUIpZoqF8whmA//a/Xw5MD3x+\nBTad21PAT4DEkY+ugeEFYAXQAzzhuI2IhHgUnJYSiglnzsCmazsXODK0zoeAI7C5IT8LXO9/Phn4\nDDZD9FH+vuYkldO1KeEBs4DNjuuLSAkp+hiCCWdgMOFMMBPVWcCP/deLsbR0BwPbsCQ0nUC//2/i\nlNSVNCXCM8+KSIUK3l6npQSXhDNx62wGrgFewmaY7gUeSCqna2Dw/B09iVVJRGQYUvQxuCacKXUB\nPxz4ItakmACMAs5P2olrU+Jk4FXgLcD9WAfIouKXr78xGKQ6W0bTOaLLcbcitWnhwoUsXLiw6vst\nxDQldvfvYE//zqRNXRLOhNeZ6H82C3gUS3oLcBdwEjadfEmugeFV/9/XgJ9j7Z2BwDCuLW7GapH6\nNGvWLGbNmjXw/qqrrqrKfuMecGpv6aS9ZTCdZG/fxvAqLgln7sZyT/wUS3PfC2wA/gR8HegA9mAJ\nZxJvIrgEhk6sF3M7ljr7/UB1fksi+5kUzzG4JJy5F7sz8SyWn7KYVGYZcDMWXApYYtsfJh3MJTAc\njNUSiuvfBvzW6UcRkSHimhKOyiWcAQsepVztL05cAsMa4FjXHYpIvIKnOR9FJKSgQVQiElYvYyUU\nGEQypBqDiERo2LWIRKjGICIRqjGISERf6QFSNUeBQSRDmvNRRCJSPvmYGQUGkQzVS+djTc35uKtv\na95FiOhdtyLvIkRsfmZp3kUo6ZUlPXkXIWJfDJ1Oo+AVnJa81VZg6N+WdxEiajEwbKnRwLBOgaGs\nFHM+ZkpNCZEM1UtTQoFBJEN9dTK6shoTvC4ETq3CfkRq2UPYFGlpuM7bCLAFGJPyeCIiIiIiIvsr\nl7x8WboJm2H3qbwLEnAI8CCWg3AlcHG+xQGgHct6tAxYBXwj3+IM0YKlVbwn74LI8LRgM9tOBkZi\nf2ThvHxZm4klBa2lwDCewfk3R2HTguf9ewKbSRzsLtfjwCk5liXoS9jkxXfnXZB6UysPOAXz8u1l\nMC9fnhZhPcO1ZD0WNAF2YNOHT8ivOAN2+f+2YkG+FnKcTsSmUp+H0itWrFYCg0tePhlqMlajWZxz\nOcD+jpZhTa8HsSZF3q4FLoM6eaKoxtRKYKjk/q5YM+IOYC5Wc8hbAWviTATeQ/r7/WmdCWzE+hdU\nWxiGWgkMLnn5xIwE7gRuBRbkXJawrcCvgONyLsdJWEr4NcDtwOlYJiapMyOA57DqcSu10fkIVp5a\n6nxswv7Ar827IAHjgG7/dQfwMPDe/IoTcSq6K1HXPoj1sj8LXJFzWcCuNOuAN7D+jwuTV8/EKVi1\nfRlWTe7BbvPm6SgsF+IyYAXWrq8lp6K7EiIiIiIiIiIiIiIiIiIiIiIiIo3j/wGmD6CqLieFpgAA\nAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10fa04c10>"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## 5. Fit voxel-wise Ct to the model to get $K^{trans}$ and $v_e$. ##"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "params, covs =  dcemri.fit_tofts_model(Ct, Cp, t, extended=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "fitting perfusion parameters\n",
        "using Standard Tofts-Kety\n",
        "fitting 30 voxels\n",
        "10% complete, 68 of 76 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "20% complete, 66 of 83 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "30% complete, 57 of 81 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "40% complete, 47 of 78 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "50% complete, 40 of 80 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "60% complete, 30 of 76 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "70% complete, 22 of 75 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "80% complete, 14 of 74 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "90% complete, 7 of 72 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100% complete, 0 of 72 s remain"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "72 s elapsed\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Kt = zeros(nx*ny)\n",
      "idxs_dce = range(nx*ny)\n",
      "Kt[idxs_dce] = params[0]\n",
      "ve = zeros(nx*ny)\n",
      "ve[idxs_dce] = params[1]\n",
      "Kt = reshape(Kt, (nx, ny))\n",
      "ve = reshape(ve, (nx, ny))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# convert Kt to min$^{-1}$\n",
      "Kt = Kt * 60.0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# PLOT RESULTS\n",
      "if plotting:\n",
      "    figure(13)\n",
      "    clf()\n",
      "    imshow(Kt, interpolation='nearest', cmap='cubehelix3', vmin=0, vmax=0.35)\n",
      "    #plot(Kt.T,'o-')\n",
      "    colorbar()\n",
      "    title('$K^{trans}$ (min$^{-1}$)')\n",
      "    savefig('qiba/qiba_v6_Kt.pdf')\n",
      "\n",
      "    figure(14)\n",
      "    clf()\n",
      "    imshow(ve, interpolation='nearest', cmap='cubehelix3', vmin=0, vmax=0.5)\n",
      "    #plot(ve,'o-')\n",
      "    title('$v_e$')\n",
      "    colorbar()\n",
      "    savefig('qiba/qiba_v6_ve.pdf')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAEOCAYAAABrZrlRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFZ5JREFUeJzt3X20HVV5x/HvvTcJSUjg8lYgELwUQRIKTSKaGAggIAIi\nVbFVC6XS1lpXkUiVpelatsuii9qlC98q8iJUhGIrKIQSUFAC4SUhkISQBAJRAgokIOQSkhCSe2f6\nxzP33nPOnJd9zszZM3Py+7Bm5dx75mUbc57z7L1n9gMiIiIiIiIiIiIiIiIiIiIiIiLS0S4GDsi6\nEZ7sCXwEmJd1Q6R9urNuQIc4Atjg6VqXAXNbPHYVcELC678OPAaMqfLeEmBqwvNLDozKugEdYA7Q\nB7wH2A24AHgwev1d4EJgHXA2cCPwYeA24FTgWuA54BzgTWAz8ECda+0H/BVwWItt/ZMWj3P1DeDf\ngI+2+TrSZgoMjc3E/rFPBP4T+3AD/BwYBH4J3As8HO3zCeCqaJ+/AF7CPuxzgPVYKr4IOImRb92P\nAf8HLGjQlk8CdwBvJfkf5GA/4OSK3/0B+FWD424HfgDsD2xsQ7tEcuU27MM75DQsYAD8JXA8cAiw\nN/C1kv2+D+wT7f+RaJ9/j967BjgI+AyWAXwNOKZBO34VXW/IeuALwErgDeCH2IfyTizlvxvordj/\n5JLXnwceB/qBn2BZjos+4F9rvPdL4HzH80hOaYyhsR7s234hMA44D/swLYne3wYciP1dzsCyhyG/\nAN6HpfB/HP059P4r2AfsWeBw4GlgTYO2HA2sLfk5xALOKcA7gLOwoPAl4I+iNl1UsX+pPwfeDxyK\nBaVPNrg+wASs6/NOqndNngT+1OE8IoX2LuwD+3Zs7GBOhm3ZgQ10DnkW67oMuRnr7gy5EOvylO5/\ncsnr0uzj68AVKbTxq1jmIgWmjKGxk7FU+wBgPvDZDNuyCRvHKFXal3+z4uft2Dd8LaUzKW822NfV\nHlg7pcAUGBo7CfgWNoB4NTZeMDmjtqzEugz1dLV47spuRqumYF0tKTAFhvpGA8cBv45+fg24CfjH\njNqzADixTeduNaCUGouNs9ydwrkkQwoMtU3H+t0h8IHodxOB3YF/AP46gzZdD5yJfQBrCSteu2YC\nzexbywexwVVfN3uJSORrtH7nY7stRnc+ioiIiIiIiIiIiEgtieeuu+kOA4I02iKSW910ExAk+rzs\ntdde4aZNzjeFbsIeystEGje1hBO690jhNPBWsJ3duutN0funNrnLY7vSatOWYDMk/7yEg+Gg0449\nXT1pXK9lWo9BxKMgLEZ2rcAg4lFRut25Cgw9XblqDqA2NSOP7cpbmwbCgayb4CRXf2ujcvZ/IqhN\nzchju/LWpiC1h1jbK19/ayIdTl0JEYnR4KOIxKgrISIx6kqISIy6EiISM4DbnY9ZU2AQ8SgIizHG\n4LLm4+nAU8AzwBfb2xyRzhY4/pe1RoGhB/geFhymYsVNprS7USKdKghDp62GRl/Sf4Yt3b8cq0g+\nVFxoMrZI72qs4vlFVY4t06gr8W6sUvP66OefRBd/stGJRSQuQTYw9CV9KvACsBQrgFT6WbwHq7MK\nVs7w51gFtZ3AxcAKrKjQY9gS/zU/x40yhoOA35X8/PvodyLSgkFCp62K0i/pnYx8SZfaWvJ6Alah\nHGw5/xXR6y1YQJhUr52NMoZijJSIFESCwcdqX9Izq+z3IeAyrNDyaVXe78Nqpiyp8t6wRoHhBcrL\nsU2OGlTmrWD78OuerlG5e3BFpFkD4QCDbXgSMsF0pWtEuTXa5gA/pryk4QSs8PFcLHOoqdEn+FGs\nRHsf8CLwMcqrKwPkbtUekaRGVXzB7Qx2pHLeWhnD0vsXs/T+ul/iTl/SJRZhn+99gFexcou3ADdg\ngaMul6WjzsCKuvZg5c0vq3g/taXdRPIqraXdVrz5jNOO08YdXnm9UcBa4BTsS/oR7Eu6dADxMOC3\nWHYxA/hp9Lsu4EdYgLjY5fouOf+d0SYiCSW4JXoAuBD4BSNf0k8Cn47evxI4BzgfG5zcAnw8eu84\n4DysWvry6HfzgLtqXSxXi8GK5FVaGcOj255y2vHY8Uemcb2WaZRQxKMaU5G5o8Ag4lFRnpVQYBDx\naECPXYtIJa3gJCIx6kqISIwyBhGJUWAQkZigGHFBgUHEp52alRCRSupKiEiMZiVEJKYYHQkFBhGv\ndqmMYfdRe6ZxGpHc2rJjcyrn2aUCg4i4UVdCRGIGCpIxuFSiEpGUtLngzLlYwZmVwIPAMSXvzcMK\nzjwB/DewW712KjCIeBQ4blW4VIX7LXACFhAuBa6Kft8HfApbB/Lo6Fwfpw51JUQ8SjD46FIV7uGS\n10uAg6PXm7F1IMcDg9GfL9S7mDIGEY+C0G2rotmqcH8LLIhevwZ8E3geW2G6HytnV5MCg4hHAaHT\nVkUzqcZ7gb9hZBziMOBzWJdiElZ45tx6J1BXQsSjWk9Xrn1wGU8/tKzeoa4FZ44BrsbGIjZFvzsW\neAirKwHwM2A2cGOtiykwiHhUa7rysNnTOWz29OGf7/jGtZW7uFSFOwT70J+HjUcMeQr4MjAO2I5V\nzH6kXjsVGEQ8SrAeg0vBmX8B9gKuiH63Exu0fBy4HgsuAbCMkRmLqlIpOLP/mMmN9xIpsI07fgcp\nFJz59ksPOu0498Dj0rhey5QxiHikZyVEJKYoS7u5TFdeC2zEbqUUkQQS3MfglUtguA6b+hCRhBLc\nEu2VS1diETZFIiIJDeYhHXCgMQYRjwoSFxQYRHzKQzfBRSqBYcvA68Ovx3TvxpjusWmcViQzO4Lt\n7AjeSv28u1TGMEFrPkqHGdM9tuwLbmtKaz4W5DYGp1mJm7AHMI7AHvu8oK0tEulgRZmudMkYKh/U\nEJEW5eFD70KDjyIeDSowiEiloowxKDCIeKSuhIjEKDCISIy6EiISU5SMQatEi3g0GLhtNSSpRAW2\nJNxy4PZG7VTGIOJRgq7EUCWqU7EVo5cC8ykvODNUiep1LIhcBcwqeX8usAaY2OhiyhhEPApDt62K\n0kpUOxmpRFXqYSwoQHklKqLXZwLX4LCWpAKDiEcZVaICuBy4BMcHPNWVEPEoQVeilUpUx0U/nwW8\njI0vnORyAgUGEY9qzUq8uHQ5Lz66vN6hSSpRzQbOxroSY4E9sDoT59e6mOpKiDhIq67Ep5Yvctrx\n6ulzKq83ClgLnIJVonoEe8CxdPDxEODXWCWqxTVOfSLwBeCD9a6vjEHEo8Gg5diSpBJVpYbdEmUM\nIg7SyhguePQBpx2vO/b4NK7XslQyhgP3nN54J5EC2/jK7xrv5EC3RItIjAKDiMQoMIhITFCQ9eMV\nGEQ8UsYgIjFB69OVXikwiHikjEFEYhQYRCQm1OCjiFRSxiAiMWGowUcRqaCMQURigsGsW+BGgUHE\no6JkDC5rPk4G7gVWA6uAi9raIpFOFjpuGXPJGHYCFwMrgAnAY8DdlK8cIyIOijJd6ZIxbMCCAsAW\nLCBMaluLRDpYguXjIVnBmUbHlml2+fg+YDq2Zr2INClBYBgqOHM6MBVb73FKxT5DBWeOAS7FCs64\nHlummcAwAbgZq2azpYnjRCQSBl1OWxVJCs64HFvGdVZiNHALcANwa+WbL219avj1hNH7MnHMvo6n\nFcmnN3b8gS07/5D6eROMMVQrODOzzv6lBWeaPdYpMHRhK9KuAb5VbYcDdz/S4TQixTFxTPkX3MZt\nT6dy3gSBIUnBmabnOVwCw3HYOvUrsUo2APOAu5q9mMgur8ZH9NVVy3h19bJ6RyYpOON67DCXwPAA\nqnEpkopaMw57HzWDvY+aMfzzuv/9YeUujwKHYxMALwIfwwYRSx0C/Az7Il/X5LFldOejiEcJuhJJ\nCs7UOrYmBQYRn5Ld1XhntJW6suT130Wb67E1KTCIeBTqISoRqVSUh6gUGER8KsizEgoMIh4pYxCR\nOGUMIlJJGYOIxCljEJFKmq4UkbiC9CUUGEQ8KsrSbgoMIj4VI2FQYBDxalfKGPaYcloap0lPVzHK\ngEmB3Dc/ldMUZIhBGYOIV0ExIoMCg4hPmq4UkUrqSohIXEG6ElrLUcSnZLUrG1WTOhKrLbEd+HzF\ne71YXZgnsRXfZ9VrpjIGEZ9an64cqiZ1Krbq81JgPuVrN74KfBb4UJXjv43Vmfgo9rnfvd7FlDGI\n+NR6jTqXalKvYCtC76z4/Z7AHODa6OcBRipWVaXAIOJRGLhtVVSrJnWQ42UPxYLGdcAyrO7E+HoH\nKDCI+DQYum1xSUYtRwEzgO9Hf24FvtToABHxpcYYQ/+zy+lfv7z6m6bpalIlfh9tS6Ofb0aBQSRH\natzI0Ns3jd6+acM/P3fff1Xu0kw1qcpnAjZg3ZAjgKexAczV9ZqpwCDiU3srUR2AZQV7RFeaC0wF\ntmCzFTcCY4DfABfUu5gCg4hPyW59bFSJagPl3Y1SjwPvcr2QS2AYC9wH7IZFm9uwatci0qwOeux6\nO/BeYFu0/wPA8dGfItKMgjws4dqV2Bb9OQbr37zWnuaIdLjqU5G543ofQzewAtgI3Ivday0izQpC\nty1jroEhAKYBBwMnACe1q0EinSwMQqcta83OSrwO3AEcCywc+uX69XcM79Dbezi9vUek0TaRzPT3\nP01//zPpnzj7z7wTl8CwLzaH2g+MA94HfKV0h76+D6TfMpEM9fYeUfYF9/xzC9I5cQ6yARcugeFA\n4EdYt6Mb+DHwq3Y2SqRjddCsxBPYgxciklRBZiV056OIR10d1JUQkbR0UFdCRNKijEFEYhQYRCRG\nXQkRiVHGICIxg8V47lqLwYr41Pry8dB6wZnJ2MOPq4FVwEWNmqmMQcSn1rsSSQrO7AQuxp6QngA8\nBtxdcWwZZQwiPmVTcGYDFhTA1n98EphUr5kKDCI+tb4eQ5KCM6X6gOnAkno7qSsh4lPr05VpTGdM\nwGpKzMUyh5oUGER8qjHG0P/yavpfqVvqIUnBGYDRwC3ADcCtjXZWYBDxaXCw6q979zmS3n2OHP75\nuTU/rdwlScGZLqwOxRrgWy7NTCUwjPvMO9M4TWq6RxfjJhIpkPtSOk/rXYkkBWemAecBK4GhOnjz\ngLtqXUwZg4hPye58bLXgzAM0OdGgwCDiU40a93mjwCDik56VEJEYPV0pIjHqSohITEGerlRgEPEp\nUGAQkUoaYxCRGI0xiEiMpitFJEYZg4jEaPBRRGKC6k9X5o1rYOjBHvv8PfDB9jVHpMN1WMYwF3uW\ne2Ib2yLS8cKCTFe6PIp5MHAmcA3xBSBEpBlh4LZlzCVjuBy4BFv8QUSSKEhXolHGcBbwMrbqi7IF\nkaSSFZzxplHGMBs4G+tKjMWyhuuB80t3WnfzNcOv9546g72nzki3lSKevbpqGa+tXpb+iYOBJEef\njq3Z2IN17b9eZZ/vAGcA24BPUr6U23nYkm9PABcAb9W6UDNZwInAF4jPSoTvv+nhJk7TflrzUdJ2\n50dnQ/KsOZwz/VKnHRct/3Ll9XqAtZRXovoE5dWkzsTWhTwTmAl8G5iFLSD7a2AKFgz+B1gA/KjW\n9ZstOKNPnEgSrQ8+ulSiOpuRD/sSoBfYH9gcHTMe6yWMx4JLTc0EhvuiC4tIq1ofY3CpRFVrn9eA\nbwLPY0vP9wP31GumStSJ+NR6xuCarVfr7hwGfA7rUkzCKlKdW+8kuiVaxKOwxj0K/VvW8/rW5+od\n6lKJqnKfg6PfnQQ8hFXDBvgZNrFwY62LKWMQ8SkIqm694w/hbfvNGd6qKK1ENQarRDW/Yp/5jMwY\nzsK6DBuxQctZwDgsozgVu5O5JmUMIh6FrU9XulSiWoDNSKwDtmJTkgArsNsMHsWmK5cBV9W7mAKD\niE/JbnduVIkKLHhU8x/R5kSBQcSnHDwH4UKBQcSjWoOPeaPAIOJTQR6iUmAQ8UgZg4jEKTCISKUE\n05VeKTCIeKSuhIjEKTCISCVlDCISF3ZWXYm6vvvhA9M4TWoO3+1tWTdBOkxaC54qYxCRGAUGEYnR\ndKWIxChjEJGYkF1o8FFEHBUkY9DSbiIehWHgtNVwOvAU8AzwxRr7fCd6/3FgesV7PVgBmtsbtVOB\nQcSjBIGhB/geFhymYsVmplTscybwdmxtyL8Hrqh4f6hqfcMVpxUYRDxKEBiSFJyBJqvWa4xBxKMg\n2NnqodWKycx02OcgbKXopqrWK2MQ8SgkcNqqHuqmMhvoooWq9coYRDyqNbC4dWc/Wwf66x2apODM\nOThUrS+ljEHEozAcrLqNHzWR/cZOHt6qaLXgzAbgn7GAcSjwcazydc2gAO4Zw3qsYu4gNvDxbsfj\nRKREgjsfkxSciTWj0cVcHxp7FngnVjU3dpGnt693PI0ferpS0tbV1QXJH7IM37HHLKcd125enMb1\nWtbMGENmjRTpFJ32rEQI3IN1Ja4Erm5bi0Q6WBC2PF3plWtgOA54CdgPuBu7LXNRuxol0qk6LWN4\nKfrzFeDn2ODjcGD4zqWXD+8484RZzDzxPWm1TyQTCxcuZOHChamfNyhIYHAZNxiPjYK+AewO/BL4\nSvQnaPBRdgFpDT72jT/Kacf121ancb2WuWQM+2NZwtD+NzISFESkCWEHLQb7LDCt3Q0R2RUUpSuh\nW6JFPApCrfkoIhWC6g9I5Y4Cg4hHnTTGICIpUcYgIjGddoOTiKRAGYOIxChjEJGYgQ57iEpEUlBj\nPcfcUWAQ8agodz5qzUcRjwLH/2pIUonK5dhhuQoMS+57OOsmxLTj0duk8tgmyGe78tamIAyctiqS\nVKJyObZMvgLD/YuzbkJM3v5hQT7bBPlsV97alKCuRKuVqA5wPLZMrgKDSKdL0JWoVWXKZZ9JDseW\n0eCjiEcDrT9d2WolqswsxBqtTVsnbwtJrpnrba44dhZwV8nP84gPIv4AKygz5ClsoSWXY0WkgEYB\nv2GkEtUKqg8+LohezwIWN3GsiBTUGcBabCBxXvS7TzNSjQps9mEdNl05o8GxIiLF1NQNGB5cC2wE\nnsi6ISUmA/cCq4FVwEXZNgew6slLsPR0DXBZts0p04OVfr8964ZIa3qwFKcPGE0++kBzsDvH8hQY\nDmBkYd4JWGqY9d8TWIkBsL7sYuD4DNtS6p+wVc0rq0JLA3m5j6HpGzA8WARsyrgNlTZgQRNgC1bt\neFJ2zRm2LfpzDBbkqxU/9u1gbDDuGnI0hVcUeQkMLjdvSLk+LKNZknE7wP4drcC6XvdiXYqsXQ5c\nAgV5nDFn8hIYwqwbUDATgJuBuVjmkLUA6+IcDJwAnJRpa+As4GVsfEHZQgvyEhhewAbWhkzGsgaJ\nGw3cAtwA3JpxWyq9DtwBHJtxO2Zjzw08C9wEnAxcn2mLpCV5vQGjj3wNPnZh/8Avb7SjR/tiD+sA\njAPuB07JrjkxJ6JZiULL2w0YNwEvAm9h4x8XZNscwEb7AyxwLo+20zNtERwNLMPatBLr1+fJiWhW\nQkREREREREREREREREREREREpHP8P4qJvKtvuLIyAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c652d10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAEICAYAAAC9P1pMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE1JJREFUeJzt3XusHFd9wPHv9doGO4/aJCnGcaq0AUqCksYhNWkMyUJp\nGhAFJCo1qC0VtCmqWohoi1qqqsrtH61aVQK1RWkKDeKRgipIQhAJTdJmwbzywHYcCKSEYASOgTxs\nJ8GAr+9O/zjr692Znd0zO3tnZi/fDxrdfcyZOULZn3/nzJnfgCRJkiRJkiRJkiRJkiRJJWwBXgf8\nV+/924E/r687WulW1d0BRXk+sAvY1Hv/UaBbX3ckNcXfAm/uvX4hcF6NfdEKZ8YwO14E3NV7vRXY\nU2NftMK16u6AorWAFwOnAfcAB4HfAzYALwC+WV/XtNLM1d0BTeyPgcPAN4AHgUfr7Y6kJvgX4KTe\n61Pr7IgqcznwdcI/Bn8x5Ps2cIgwUb0L+OsCbbVCbAfeAFwGnFFzX7T8WsBDwJnAGmA3cHZqnzZw\n84RtB6wu01PV6vO9TT8dthF+3Ht77z8KvBb4Wmq/YdMDsW2XeFVCmg2nA9/pe//d3mf9EuBi4D7g\nFuCcAm0HmDFIsyGJ2GcnYVh5GHglcBNhcVxhpTOGOeYSQqfd3FbsdumllyaUtHHjxiLnfDLVfB+D\nc0lnEP7l7/cUISgA3EqYT3hWb79xbQeUzhgSEtbNrS97GAAWkiOsmVtb+jhnnPDCKfQmeOwn+zj1\nGSOzrminbPqVqRznO49/iTNOuWgqx1pz2aumchyAvXd/iDO3/e5UjvWZ9/z6VI5z9dVXc/XVV5c+\nztxc+Sv7Bw4cYDFZjNq3Ndc6KfXRvcDzCBOIjwC/RZh87vds4AeEwLKNMN/wRGTbAQ4lpAp1k4lv\ncTkK/Anw34SrDP9BmDx8S+/7a4HfBP6ot+9h4IoxbXMZGKQKdcvd+3Zrb+t3bd/r9/S22La5GhUY\nVjVwhfb6TEZXv5PXbam7C0NtOL1593W12+26uzDgaHK07i5EaVRgaM01MDCsPrnuLmT8zPqmBoZf\nqrsLGU0LDF1Kz2FWolGBQVrpSg4lKmNgkCpUYvKxUgYGqUIOJSRlOJSQlOFQQlLGUeJWPtbNwCBV\nqJs4xyApZVbmGGLurrQklDQl3SSJ2uo2LmNoAf8KvIJw2+c9hNJRI2/AkDTcSskY+ktCLXC8JJSk\nCSySRG11GxcYCpeEkpSv5FAidlj/y4RbrV/f99lewkOKdgF3j+vnuKFEVOhaSI4svV5Fq5E3Q0lF\ndDodOp3O1I9b4nJl7LC+BfwD8OnU5wmhivQTMScbFxhiyklNpeqS1CTtdnvgzsz5+fmpHLfExGJs\npee3Ah8jZA1p0WWoxg0l+ktCrSWUhBpWt15ShC5J1DZEzLD+dEKwuKb3vv9ACXAH4Td95bh+jssY\nCpeEkpSvxJLomFTj3cBf9vadYzBD2A7sJzz79HbCXMWOvAPFLHAqVBJKUr68uyvv/exdfPmzI+cE\nY4b1LyIMMSA8tvCVhKuJNxOCAoRnnN5IGJqUCgySpiTvUuTWS7ax9ZJtS+/f+3eZ0o0xlZ5/oe/1\n+4FPEoLCekLG/xRwAuGxhiMnTQwMUoVKTD7GVInOswm4ofd6NXA9cNuokxkYpAodLXfb9bgq0f3e\n1Pf6YeD8IicyMEgVsoKTpIwm3CAVw8AgVciMQVKGgUFSRnc24oKBQarSgsVgJaU5lJCU4VUJSRmz\nMZAwMEiVMmOQlGFgkJThUEJSxlEzBklpszKUiHkSlaQp6UZuOcqUjy/0RDkDg1ShEs+VOFY+/nLg\nHEL1prNz9kuXj49tu8TAIFWom8RtQ8Q+Fe5Y+fhHJ2i7xMAgVaim8vGFnyjn5KNUoby7Kx/8/E7+\n7ws7RzUtUz6+8IyngUGqUN7lyrMu3spZF29dev+pf7ouvUuZ8vFRT5TrZ2CQKlSiHkOZ8vGrI9oO\nMDBIFSqxiqFM+fjCT5QzMEgVKrnAadLy8XltcxkYpApZ2k1SxqwEhph1DNcB3wfuX+a+SCteySXR\nlYkJDO8nLKWUVNJiN4na6hYzlNhBuMwhqaQG/OajOMcgVagJw4QYUwkMC8mRpderaNGaa03jsFJt\nOp0OnU5n6sf9qcoY1sytncZhpMZot9u02+2l9/Pz81M57ozUaXEoIVVpVjKGmKsSHwG+ADyfcOtm\nekWVpEgl6jFUKiZjGHmzhaR4iw340cdwKCFVyDkGSRlNGCbEsLSbVKGScwzjKj2/FrgP2AV8GXh5\n33d7gT297+4e108zBqlCJYYSxyo9v4JQkekeQhGW/roKdwCf6L0+F7gReO6xUwNt4ImYk5kxSBVa\n5irRP+x7fSLwWOr7OSIZGKQKLXbjtiFiKz2/jpBF3Aq8re/zhJBR3AtcOa6fDiWkCpUYSsS2vKm3\nvRT4EPCLvc+3A/uB04DbCXMVO/IOYmCQKpQXGPbds4t99+wa1bRopecdhN/3KcDjhKAA4UE0NxKG\nJgYGqQnyrjg858KtPOfC4+Xj77nm/eldYqpEnwU8TMguLuh99jiwnjB5+RRwAnAZMPLmDwODVKES\nQ4mYKtGvB95ImJx8Grii990m4Ibe69XA9cBto05mYJAqVHKB07gq0f/Y29IeBs4vciIDg1Qhl0RL\nyljsRi8lqJWBQaqQGYOkDAODpAwDg6QMA4OkjO6M1I83MEgVMmOQlNH1cqWkNDMGSRkGBkkZiZOP\nktJmJWOwtJtUoSSZi9pylKkSPa7tADMGqUI1VYmOaTvAjEGqUHcxbhuiTJXomLYDDAxShZIkbhui\nTJXo2LZLYgLDGcCdwFeBrzBYklpSEUnkNrxljJuAs4HfIFSJnmhFVcwcwwLwdmA3IT35MqH8dO74\nRNJweZcrH9uzk8f27BzVdNIq0c/q7VekbVRg+F5vg1Bg8mvAZgwMUmF5k4+nnHsBp5x7wdL7B//z\nuvQuZapEH4poO6DoVYkzga3AXQXbSaK2KtF5bXMVCQwnAh8DruqdVFJBSbmbqCatEp3XNldsYFgD\nfBz4MGFyY8BCcmTp9SpatOZaseeXGqnT6dDpdKZ+3JW0JHqOkHo8ALx72A5r5tZOs09S7drtNu12\ne+n9/PzIBzdFm5XAEHO5cjvwO8DLCEstdxGWV0oqavLLlZWKyRg+hwuhpKmYlZuovFdCqtCsDCUM\nDFKVzBgkpSXDb5BqHAODVCHnGCRlOccgKc2MQVKWGYOkNDMGSVlmDJLSZuVypUudpSqVKPrI+BLw\nv00oH78H+DxwXt93e3uf7wLuHtdNMwapQiWWRMeUgH8YuIRQsely4N+Bi46dGmgDT8SczIxBqtLk\nd1fGlID/IiEoQKiytiX1fXSVGAODVKVu5JZVtAT87wO39L1PCA+kuRe4clw3HUpIFcqbPjjwjZ0c\n/MbIKtFFLnS+DHgzoZbKMduB/cBphCrvXydUkh7KwCBVqTv8973xrK1sPGvr0vu9t2aqRMeWjz8P\neC9hjuFA3+f7e38fJTy6bhsjAoNDCalKi5FbVn/5+LWEEvA3p/b5OeAGQsW1h/o+Xw+c1Ht9AnAZ\ncP+obpoxSBVa5vLxfwNsBK7pfbZAyAw2EQIGhN/89cBto05mYJCqlDOUiDSufPwf9La0h4Hzi5zI\nwCBVyXslJGV4r4SkjBm5vdLAIFXIKtGSshbNGCSlmTFIynCOQVKGGYOkjBWUMTwT+AzwDMIa7U8A\n71zOTkkr1grKGH5MuI3zcG//zwEv6f2VVMQKyhggBAUIGUOLyPJQklJm5HJl7G3Xq4DdwPeBO4EH\nlq1H0krWTeK2msUGhi7h7qwthGKT7eXqkLSSJd0kastRpkr0uLYDil6VOAR8CrgQ6Bz7cCE5srTD\nKlq05loFDys1S6fTodPpTP/AkycDZapEx7QdEBMYTiUUiTgIrAN+DZjv32HN3NqIw0izo91u0263\nl97Pz8/n71zE5MOE/irRcLxKdP+P+4t9r/urRMe0HRATGJ4DfIAw7FgFfAj4n4h2ktImvyoxrEr0\ni0fs318lumjbqMBwP3BBxH6Sxpn8qkSZKtGFT+rKR6lCczlDiYP77uPgI3tGNS1TJTq27RIDg1Sl\nnKHEhs3nsWHz8YsI3773+vQu/VWiHyFUiX5Dap+8KtExbQcYGKQqTT75WKZKdF7bXAYGqUr1VInO\na5vLwCBVaYXdKyFpGhqw3DmGgUGq0uJs3HdtYJCq5FBCUoZDCUkZZgySMswYJGWYMUjKMGOQlLG4\nWHcPohgYpCo5lJCU4VBCUkYyGysfY6tES5qGcuXjx1V6fgGh7uOPgT9LfbeXUD16F3D3uG6aMUhV\nmnyOIabS8+PAW4HXDTsz4bEPUQ+LMmOQqpR047as/krPCxyv9NzvUUK1poWcs8/FdtPAIFVpsRu3\nZQ2r9Hx6gTMnwB2EwHHluJ0dSkhV6k48+Vj2csZ2YD9wGnA7Ya5iR97OBgapSjlzDAefeJCDBx4c\n1bJwpeeU/b2/jwI3EoYmBgapEXIuV27Y+Dw2bHze0vtvf+uT6V2KVHpOzyWsJ0xePgWcAFxG6mly\naQYGqUrLWyV6E+FqxcmEB1FfBZwD/CyhrDyE3/z1wG2jTmZgkKpUboHTuCrR32NwuHHM04Sn1Ucz\nMEhVmnzysVIGBqlKXe+ulJQ2IxlD7AKnFmGNdWaqVFK8JEmitrrFZgxXAQ8AJy1jX6SVbwXdXbkF\neBXwPgqstZY0RLcbt9UsJmN4F/AOwrVRSWU0YJgQY1xgeDXwA8L8Qjtvp4XkyNLrVbRozbWm0Tep\nNp1Oh06nM/0Dd49O/5jLYFxguBh4DWEo8UxC1vBB4I39O62ZW7ssnZPq0m63abfbS+/n50euII6W\nNGCYEGPcHMNfEVZS/TxwBfC/pIKCpAImr8dQqaLrGGZjgCQ11QqZY+j3md4maVINyAZiuPJRqlBi\nYJCUsUImHyVNUdI9GrXlKFM+flzbAWYMUpUmH0qUKR8f03aAGYNUpXrKx8e0HWBgkCqUJN2obYgy\n5eMLt3UoIVUpZ/Lx0OHvcuhHI4s+l1kAUbitgUGqUN7lypPXbebkdZuX3n/3QObxkmXKxxdu61BC\nqtLkcwz95ePXEsrH35xzlnR5hCJtATMGqVIjLkWOU6Z8/NM5bXMZGKQKlVz5OGn5+Ly2uQwMUpVc\nEi0pzXslJGUlPldCUooZg6QMA4OkjBKXKytlYJAqZMYgKSPByUdJaWYMktIcSkjKMDBIyjAwSMro\ndtNV15rJegxShRK6UVuOmErP/9z7/j5ga9/ne4E9hAdUZ6rApJkxSBUqMZSIqfT8KuC5hKIsLwau\nAS46dmrCE+ufiDmZGYNUoSRZjNqGiKn0/BrgA73XdwEbgGf3fZ+u7JTLwCBVaJmrRI/aJwHuIJR5\nu3JcP2OHEnuBJ4FFQrTaFtlOUp+8+YPDRw9x+OiTo5vGycsKXgI8ApwG3E6Yq9iRd5DYwFBofCJp\nuLw5hnWtk1jXOmnp/eNH9qV3ian0nN5nS+8zCEEBwkNpbiT8454bGIoMJaLHJ5KG6yYLUdsQMZWe\nbwbe2Ht9EXAQ+D6wHjgWdU4ALgPuH9XPIhnDHYShxLXAeyPbSepT4qpETJXoWwhXJh4Cfgi8qffd\nJuCG3uvVwPXAbaNOFhsYtgP7yRmfLCRHlnZcRYvWXCvysFIzdTodOp3O1I/bXd4q0RCCR9rDwPlF\nThQbGPb3/g4dn6yZW1vknFLjtdtt2u320vv5+fmpHHfE4qVGiZljKDw+kTRciXUMlYrJGJ5NyBKO\n7T92fCJpuJJDicrEBIZvUXB8Imm4bmLNR0kp3RmZYzAwSBVqwvxBDAODVCEzBkkZVnCSlGHGICnD\njEFSxtHhN0g1joFBqtCsLIk2MEgVmpWVj5Z2kyrUjfxfjjJVomPaLmlUYFhs4OKPMeW2anHocLpw\nTzMc3Hdf3V3IWI5bp8voJt2obYhjVaIvB84B3gCcndqnv0r0HxKqRMe2HdCowNBt4JOADy8+VXcX\nMp78UVMDw566u5DRtMBQ4rkSk1aJ3hTZdkCjAoO00pUYSpSpEr05ou0AJx+lCh2d/O7KslWiK9ch\ndNrNbSVvHcorcr705NZFwKf73r+T7CTivwFX9L3/OqGeSkxbSTNoNfBNjleJ3s3wycdbeq8vAr5U\noK2kGfVK4EHCROI7e5+9heOVoiFcfXiIcLnygjFtJWk2FVqAUYHrCA/raFLh2zOAO4GvAl8B3lZv\ndwB4JuHS2G7gAeDv6+3OgBbhse+frLsjmkyLkOKcCayhGWOglxJWjjUpMGzieP3NEwmpYd3/P0Go\nJA5hLPslwnMSm+BPCcWL009s0hhNWcdQeAFGBXYAB2ruQ9r3CEET4GnCk4g219edJYd7f9cSgnwT\nnnG6hTAZ9z4adAlvVjQlMMQs3tCgMwkZzV019wPCf0e7CUOvOwlDirq9C3gHzMjtjA3TlMCQ1N2B\nGXMi8DHgKkLmULcuYYizBbiE8GT0Or0a+AFhfsFsYQJNCQwxj/hWsAb4OPBh4Kaa+5J2CPgUcGHN\n/biYcN/At4CPAC8HPlhrjzSRpi7AOJNmTT7OEf4Df1fdHelzKuFmHYB1wGeBX62vOxmX4lWJmda0\nBRgfAR4BfkKY/3jT6N0r8RJC2r6bkCbvIlzmrdO5wE5Cn/YQxvVNcilelZAkSZIkSZIkSZIkSZKk\nleP/AfCs0gnR8d0jAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10fe95410>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Kt_truth = array([0.01, 0.02, 0.05, 0.1, 0.2, 0.35])\n",
      "ve_truth = array([0.01, 0.05, 0.1, 0.2, 0.5])\n",
      "Kt_error = ((Kt.T - Kt_truth) / Kt.T).T\n",
      "ve_error = (ve - ve_truth) / ve"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "if plotting:\n",
      "    figure(10)\n",
      "    imshow(Kt_error, interpolation='nearest', cmap='bwr', vmin=-0.02, vmax=0.02)\n",
      "    title('relative error in $K^{trans}$')\n",
      "    colorbar()\n",
      "    savefig('qiba/qiba_v6_error_Kt.pdf')\n",
      "\n",
      "    \n",
      "    figure(11)\n",
      "    imshow(ve_error, interpolation='nearest', cmap='bwr', vmin=-0.01, vmax=0.01)\n",
      "    title('relative error in $v_e$')\n",
      "    colorbar()\n",
      "    savefig('qiba/qiba_v6_error_ve.pdf')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARQAAAENCAYAAADOqZoqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGVJJREFUeJzt3X+8VHWdx/HX5QL+ui1m9hAVEStcMTErAzTxDohFVGpr\nWRb5ox5b265pP1YDU0EfPTazH5bu1rqZ/TK0VsswM2U3hwslpAgqKpituLpgqegCgQt07/7xOcOc\nO/ecme+c+c45Zw7v5+NxHnfOzPc78z1w53O/5/s95/MFERERERERERERERERERERERGRzvFpYHTW\njZD2GpZ1A1K2DjgpYd3VwIn+mpIraRzb4cCzbf4MydjwrBuQsoFgc7EO+Ajw62D/qHY0KCfafWxT\ngXHAccAewLnAb4LH1wLnAU8ApwA/At4D/ByYAdwAPAWcDmwDNgFL29xeSahIPRTfwXEA6PL8nq2I\nOr5mj7ldf0AmA3cBvwU+FHr+Z8AtwATgHuBeYAXwMvBvWDA5A9iABYkXsUA+ClgC7ARGBu/1fmB/\nFExyrdMDyjrgIuAhYDNwCHAr8Cfgv4BP1qk7B/uruAl4BDgt9NoPgbHA7cH7Xhh81vTg9c8B/17z\nft8INoCDmmhHvbK1x9cd8dww7Atbxr6Qq4F313mPqP/zdQw+FVwHfBZ4EHgJuBnrTcRZjgWJb2A9\nDIC3AVcC7wW2YMFmLDACeD5UtxS0/Xgs2HQDLwSvHQpsBT4BzAX+Gji6TjtEWrIOeAA4GNgb+4W8\nBPtLfBjwB+wXu+JJqkHhvVQHCc/AfulHx5St3R8L/BnoCfa7gfXAJOwL26gdFY3Kho9vj5jnRmCB\ncU7wHtOwIDm+znvUijrWZdi/xyuBR4GPx9StHP9G4ABgL2B28LjiNOB92GnPjGCrOBX4APAZ4B+B\nWcDbg9e+CLwVmBlsZ7P7naZLip4EzgkeT8bOtcPmYufg4fLTibYSO4ePK1u7vwT4cPD4ZOxL7doO\nHMuGj4+Y56ZipwxhC4B5dd6jVtSxfjC0/yXgW3XqvwULOq/DxkamNvg8KagiRPung5+HYqcPL4Ze\n6wb6YuqdhU1ljgv2e4BXNfG5C4AzsdOjD1Lt6jfTDpeyTzNU+LmDIso8FTxf7z0aCc/IbKt5v1rT\nsVOj0cBC7LRtSYLPlA5XhIBSmbV5GvvLerhDnUOxQcHp2EDhANZDCQ/CNpoNugX4KnYqcRowJXj+\nv5toh0vZqHaEn1uPjR11hZ4/FFjT4D2a0ah+Cfg6NmD6KNbLOoRkgUw6WKcPyoYtxwYdL8LO47ux\n6dBjI8rug31Jnsf+Dc5l6NTpH4HX1vm857DBxO9hg6lrg+d/10Q7mikbZxk2cHkRNp5SAt6FDaT6\nUm+2awQ2zlGZXt8I3AT8g8fPlw5RpIDSj32RjsG+4M9hvZC/iij7KNa7uBfr2h/F0OnIL2KDpS9i\nMx5Rf6UXYLMjCxK2o5mycXZgszrvCOr/Mza283gT79FI3PU7b8TGVwaAdwbPvQIL2H+HDaKKiIiI\nSJHNxMbGfo9dCxXlmuD1B7EeJNh41j3Y9VargfND5fcDFmE92ruBfb23WkRypxu7NGEcNm61Crug\nMWwW8Mvg8WRsfA1s9u2Y4HEPNt53RLB/FTb+BhakrvTcbhHJoeOAX4X25wRb2L9ityhUrGHwBYYV\nt1G9MjpcZjSDZwebVqRBWZEiO5jB0/DPBM81KjOmpsw47FRoebB/ADajSfAzKgA5U0AR6Qyu1xLV\nTvGH6/Vg109dgN1qEvUZLV2z1HJAmTq1t9IIbdoKu514Yu8ALXplc5+5qab6/2CDqxWHYD2QemXG\nBM+BjbvcCtyInfJU/JHqPWwHYjepZmpg+/YBL9sll8zz8j4D/f3etnmXXebtvXw167LL5nl7r4Gd\nO71t8y691Nt75e3fCvuSt/xdGRg2zGmL+Lzh2I2j47CUDo0GZadQHZTtAn4AXB3RpquozhjNocVB\n2SJcei/SOYY5nhT099c+sxNLRHUXNuPzHeAxqneBX4cFk1nYbNCfsSvAwa5kno2lsFgZPDcXG+S9\nEvgJ8FHszvQzmjmcWgooImlyDSjR7gy2sOtq9s+LqLeU+OGNjQxOJ9GSXAWU3t5S1k0YolQqZd2E\nIfLYJoBSb2/WTRgid/9Wwx2/ctu3t7cdbeIjxeHA9u0+Ti/9GTE8X+2pGMhVRknT1f+XrJsQaWBY\nd9ZNGGTYsC5o/fsyMNDT07gU0LVli4/PS12ueigihdfaKU/uKaCIpEkBRUS8UUAREW8UUETEGwUU\nEfHGddq4QxX76ETyRj0UEfGm4AHF5ehc0s6JiIthw9y2DtWoh9KNZVGfgd0GfR+2kNNjbW6XSDF1\ncLBw0ejoJmF3Lq7Dlmu4GVuLVkSS2M17KFEp5Sa3rzkiBdfBwcJFo4DidJfdFVfM3/W4t7eUy7uG\nRZpRLpcpl8v+37i1aeOZ2JKv3cD12CJrta7BFn3bCpxDNf/JDdhibH8CJtbU+STw98BfgDtoYay0\n0d2MU4D52IGAJWXpZ/CB6G5jR7rb2F1h7zY+qnbF22hdq1fXfl43tvxFeDzzTAaPZ87C8qHMws4k\nvkF1ze2pWB7ZHzA4oEwDLg7q7ABeja1AmUij/tf9wHiqaefejw3KikgSycdQXMYzTwG+Hzxeji3a\nVckXuwRbVrfWJ7Bld3cE+4mDCTQOKOG0c48CP0YzPCLJJQ8oSZfRqC1TazxwIpZ/tgwc28TRDOFy\nQheVdk5EkogZlC1v3kx58+Z6NV3P4+stoxFlOPBK7NToLVh+2dc4flbkm4lIWmICSmnUKEqjRu3a\nv/zZZ2uLtLqMRpxngJ8Gj+/DxkhfBbzQoF6kYs9hieRN8lMel/HMhcBZweMpwEtUVwWMcxswPXh8\nePDeiYIJqIcikq7k08atLKMBcBPQi/U+ngYuA76LTSffADwMbKcakBJRkuoUadrYXWGnjU84walg\n19KlPj4vdeqhiKRpN79SVkR8UkAREW8UUETEGwUUEfFGOWVFxBv1UETEGwUUEfFGAUVEvFFAaSxv\nV6bu2JnPCwzz+LvU35+vK1J36c+6AW2Sx18Cj9RDEUmTAoqIeKNpYxHxpuA9lGIfnUjetLYuj8sq\nntcErz8IvDH0/A1YbpSHa8p/GUuD8CCWaGkULVBAEUlT8oBSWcVzJnAklvF+Qk2ZWcDrsERMHwO+\nFXrtu1RXrwi7G3g98AbgcWxli8QUUETSlL+s94uozqktx9JGJqaAIpKm/GW9D/sIlvUtMQ3KiqQp\nLuv9+vWUN2yoV7NdWe8rPo+lgFzgWD6SAopImmKmjUtjx1IaO3bX/uUrV9YWaVfWe7AlS2cBJzmU\nrUunPCJpyl/W+5nAhdh4zMuJjilEPRSRNCW/DqVdWe+vxQLUoqDcvdjC6Yl4yXo/0J+vGy90L4+7\nnP3X5dbIkZ6y3p9/vlPBrmuu8fF5qVMPRSRNefyr4pECikiaCh5QXI4u7pJdEWlWa5fe555Ly+Mu\n2RWRZg0f7rZ1KJeWL8GmqkSkVR3c+3DRuaFQpBMpoDQ2f/78XY9LpRKlUsnH24pkZvHiMosXl/2/\nccEDius89zjgdmBixGu6DsVRHn+XcvZfl1verkMJ/fGtp8vK5fMXuQ6d8oikKY9/VTxyObqbgN8C\nh2OX7J5bv7iIxCr4tLFLD+XMtrdCZHfRwVPCLop9dCJ508G9DxcKKCJpKnhAKfbRieRNdlnv4+pO\nAn4HrATuA96S9NBAAUUkXdlkva9X9yrgUiz4XBbsJz+8ViqLSJOyyXpfr+4Gqmvx7ItbyshYGkMR\nSVPyWZ6ojPaTHcocDBxUp+4cYCnwFayDcVzSBoJ6KCLpSt5DSZr1vpHvAOcDY4FPY+lKElMPRSRN\ncctorF1Lee3aejWTZr1/BhhRp+4kYEbw+Bbg+nqNaEQBRSRNMQGlNGECpQnVMdbLf/GL2iLhrPfr\nsaz3tRedLsQSWd/M4Kz3L9Sp+wSWvHoxMB1bjjQxBRSRNGWT9T6uLths0L8AewDbgv3ElPU+RXm8\npiln/3W55e1u4+99z6lg1znn+Pi81KmHIpKmPP5V8UgBRSRNujlQRLxRD6Xz7NyZdQs6R8F/v/On\n4P/ghQwoIrmlgCIi3iigiIg3Cigi4o0Cioh4o2ljEfFGPRQR8UYBRUS8UUAREW8UUETEm4IHlGIf\nnUje5G8ZjYrPAv3Afs0fVJV6KCJpSj5tXFkKYwaW6vE+LEPbY6Ey4WU0JmPLaExxqHsIcDLwVNLG\nVaiHIpKm/C2jAfA14CIvh+dQ5hDgHuARYDWWIVtEkkgeUOKWyHApE7WMRqXuqcH+Q8kPqsql/7UD\nS6+/CugBVgCLGNzVEhEXcVnv77+f8ooV9Wq2YxmNvYCLsdOdJPWHcAkozwYbwBYskByEAopI8+Ky\n3k+aRGnSpF37l3/727VF2rGMxmuxTPgPhsqvwE6R/tTgSCI1O4YyDhs5Xp7kw0R2e8lPecLLaIzE\nlsJYWFNmIXBW8Di8jEZc3dXAAcBhwfYM8CYSBhNobpanB1sI6AKspyIizcrfMhphrqdVsVwDygjg\nVuBG4LbaF+fPn7/rcalUolQqtdoukUwtXlymr6/s/41bu9v4zmALu65m/7wm6tZ6TZJGhbkMwHRh\nU1EvYIOztXK3Ls+2lztuOZPMFPzCTW/23NPTujxr1jgV7DriCB+flzqXX6e3ArOBacDKYJvZzkaJ\nFFZrV8rmnkv/aym6AE7Ejw4OFi506b1ImhRQRMQbBRQR8UY5ZUXEG/VQRMQbBRQR8UYBRUS8UUAR\nEW8UUETEGwUUEfGm4NPGxQ6XInmTv6z3+2EZGB8H7sby0CamgCKSpuQBpZK5fiZwJHAmMKGmTDjr\n/cewrPeN6s7BAsrhwH8G+8kPr5XKItKk/GW9D9f5PnBaS4fXSmURaVL+st4fgKWJJPh5QKLjCngZ\nIbr9F/nKAzNjRtYtiJbHAf6c5cbaJY//Vl7EZb0vlymXy/Vq+sx63xXzfgNNfE6kYg85i+TMQMz3\nvbc0jd7StF37l19xRW0Rn1nvxwRlwXolo7GVLQ6khQTVoFMekVTt3Om2RWhH1vtKnbODx2cTkTO6\nGeqhiKSohVPMdmW9vxL4CfBRbND2jMQtxE8S3IGFC1vOvu+VxlDcaQzFja8k1Vu3un1X9t7by+el\nTj0UkRTlNYD7ooAikiIFFBHxRgFFRLxRQBERb2KmhAtDAUUkReqhiIg3Cigi4o0Cioh4o4AiIt4o\noMCewGJgD+zGop8Dc9vZKJGiUkCBl4FpwNag/FLghOCniDRB08Zma/BzJHa34sb2NEek2IreQ3G9\np3MYsArLrXAP8GjbWiRSYP39bluTXDPXJ818PxbYAny2UUNceyj9wDHAKCynQgkoV15csGD+roIT\nJ5aYOLHk+LYi+bR4cZm+vrL3921TD6WSuf4qLFDMYWj2+krm+xlYtrb7sORKjznU/xpwh0tDkuRb\nuBTYBnwl2Fc+FEd5y/EB+e2C5+3fylc+lIcecvuuHH10U5+3Builms6xDBxRU+Y4YB7WS4FqwLiy\nQf3TgOOxhE1bgK/Wa4jLf9v+VLtAewEnAysd6olIjTad8rhkrq+XNT+ufg9wETDftSEupzwHYut1\nDAu2H2ILAolIk1qY5VmE9R5qfb5mPy5zfe1zLpnv5wNXY5MyTr0ll4DyMPAmlzcTkfrieh8rVpRZ\nsaJcr+rJdV5zyVwflRG/Ueb7ScDp2NjKvthY6jbgm3ENUU7ZFOVtXAA0huLK1xjKsmVu35UpU5r6\nvKuAF4AvYWMj+zJ0UHY4sBY4CVgP/A5bkvQxx/rzgM3YAG2snP23iRRbm8ZQrsR6MI8D04N9sBUD\nK7Mz4cz3jwI/ZnDm+6j6TdO9PCIpalOPcCM2HVxrPfDO0P6dweZaP+xyl4YooIikKK+nmL4ooIik\nSAFFRLzRzYEi4o16KCLijQKKiHijgCIi3iigiIg3Cigi4o0Cioh4o2ljB6ec4uNd/HnppaxbEC1v\nN7wBdPm4PVScqYciIt4ooIiINwooIuJN0QNKDs/qRYqrw5bR2BO4CXgIy6FSm3RpCAUUkRS1KaBU\nlsE4HMv3HPXFryyjMRM4EsvWNqFB/Q8EP48G3gx8HFujJ5YCikiKdu5025p0CpZInuDnaRFlJgFP\nAOuAHcDNwKkN6m8A9sGC0T7AdmBTvYYooIikqMOW0bgLCyAbsED0ZaDuRRkalBVJUVywWLu2zOOP\nl+tVzWIZjdnYWlwHYuMsS7BToifjGqmAIpKiuIAyfnyJ8eNLu/bvuGNICtcsltE4HvgZ8BfgOeA3\nwLHUCSg65RFJUZtOeRYCZwePzwZuiyhzPzAeGAeMBN4f1KtXfw2WBR9sDGUK1Uz5kRRQRFLUYcto\nXIcFn4exdXxuAFbXa4iXhb6iT8Wyo3t53OX1Xp7u7qxbMNjee/tZ6OvrX3f7rnzqU14+L3UaQxFJ\nke42FhFvin7pvWtA6cYGdZ4B3t2+5ogUmwKKuQAbyHlFG9siUnhFDyguw4RjgFnA9XTgIJFInrRp\nlic3XHooVwMXAn/V5raIFF4nBwsXjQLKu7Cr5lYCpfhi80OPS/WLinSAvr4yfX1l7+9b9FmeRqcw\n/wR8GLsoZk+sl3IrcFaojK5DcaTrUNwV9TqUSy5x+6584QudeR1Ko1/xi7Hr/w/DciP8msHBRESa\noDGUwfLVFRHpMJ0cLFw0E1AWB5uIJKSAIiLeKKCIiDdFDyg5nHcQKa425ZRtNev9+4BHsERKbw49\nfzJ2y81Dwc9pjRqigCKSopxmvX8YeA/Qx+CJl+ewa9GOxhIv/bBRQxRQRFLUpoDSatb7NVjvptYq\nLC0k2L18ewEj6jVEYygiKWrTGErSrPeTm/iM04EVWDCKpYAikqK4gLJ+fZn168v1qvrOet+M11NN\nE1mXAopIiuICyujRJUaPLu3af+CBtme9PwTrpTQyBvgpdgtObLb7Co2hiKQop1nvw8L3D+2LJbn+\nHHCvS0MUUERS1KZp41az3r8HG1+ZEpS/M3j+POC1wDws48BKYP96DVHW+xTpbmN3Rb3bePZst+/K\njTd25t3GGkMRSVHRr5RVQBFJkQKKiHijgNKBnn8+6xZE27/ucFY28jZWUXQKKCLijQKKiHhT9CTV\nCigiKVIPRUS8UUAREW8UUETEGwUUEfFGAUVEvFFAERFvij5tnMP7X0WKq035UHxmvX9TTZ2jsVwo\nq7Hs93vUa4gCikiKOiDrfdhwLNP9x4CjgF4a5JRVQBFJUYdlvX8b1it5ONh/EajbOgUUkRS1KaAk\nzXp/cIP3HY9lT/sVlvH+wkYNcR2UXQdsws6xdmDRTkSaFBcsNm0qs3lzuV7VLLLejwBOAI4FtmGn\nUyuAX8dVcA0oA0AJ2JigUSISiAsoPT0lenpKu/Y3bMhF1vunsXGVyvf+l9igbWxAaeaUp+PyW4rk\nTZuSVLcr6/1dwERsxcDh2KDsI/Ua4hpQBoD/CBr1t451RKRGm8ZQ2pX1/iXga8B9WMb7FaHXIrme\n8rwV2AC8GjuXWwMsqb48P1S0FGwinauvr0xfX9n7+7bpStmNwIyI59cD7wzt30l0QPhZsEX5UbA5\nSXIaMw/YAnw12M/dMhpPPJF1C6IpBaS7vLXL1zIaRxzh9l1Zs6Yzl9FwOeXZG3hF8HgfbG764fji\nIhKnTac8ueFyynMA1e7QcKz7c3fbWiRSYJ0cLFy4BJQngWPa3RCR3UHRbw7U3cYiKVIPRUS8UUAR\nEW8UUETEGwUUEfFGAUVEvFFAERFvNG0sIt6ohyIi3iigiIg3RQ8oOcspW866AUMsW1bOuglDLFlS\nzroJkfLYrnakIGhFTpfR+DKWG+VB4KfAqNBrc4Pya7Abg+tSQGlg+fJy1k0YYunSctZNiKSA0lhO\nl9G4G3g98AYsKM0Nnj8Sy+x2ZFDvmzSIGTkLKCLFltNlNBZRXR5jOTAmeHwqcFNQfl1Qv26Ceo2h\niKSoTdPGSZfRmBxR7iNYEAFLIbmspk6jpTdaVqaaul+btqJuZVrXzOdtqqm7CEtsVrudgi3AFRa1\nOsXpwLdD+7OBa2vKfB64NbR/LfCh0P71wN/EHRz46aGUPLyHyO6glZSO7V5G4xxgFnBSnTpjgudE\npMCuojprM4dq1vuw4cAfqC6jsYrqoOxMbHmM2qzHRwblRgKHBfU7Ls+tiDRnP2yZm9pp4/AyGgDv\nANZig6tzQ8//HngKWypjJTabU3FxUH4N8PY2tF1EJP/iLrrJyg3YuWmeMvwfAtyDdU9XA+dn2xwA\n9sSmGldhC0h9MdvmDNKN/cW9PeuGSLq6sW7VOGyB5vD5XVamAm8kXwFlNNWE4T1Y9zXrfyewpVbA\nztOXYQts58FnsFUaopbclDbIy4Vt9S66ycoShk7HZe1ZLNiCLbb2GHaenLWtwc+R2B+HqGnLtI3B\nZi2uRwOJqclLQIm66KbtF9B0uHFYD2p5xu0A+z1ahZ0i3oOd+mTtauBCqleASgryElAGsm5Ah+kB\nbgEuwHoqWevHTsXGACeS/bVJ78KuxViJeiepyktAaXTRjVSNwK5mvBG4LeO21PpfbJry2IzbcTx2\nBemT2GXk04EfZNoiSVW9i26yNI58Dcp2YV+Mq7NuSMj+VK972AvoY/DVllnrRbM8u6W4i26ychOw\nHvg/bHzn3GybA9jsST8WcCsXIc3MtEUwEXgAa9ND2LhFnvSiWR4RERERERERERERERERERERkY7y\n/8RRR2fXjRxVAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x110267190>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAARQAAAEMCAYAAAAF9UmPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGHxJREFUeJzt3X+0VWWdx/H3hSuigpBpoiHROLCUjDFigAK9B80ixtJm\nKqvlj5qaWtNYrWoysFZezZVoU2bNVM4QU1OpOVoOZVlWHrmMv/khCF6IggYDf6SYGBYiZ/747sPd\n99yz93nOOc/+cfb9vNY6655z7vPs893o+d5nP3vv7wMiIiIiIiIiIiIiIiIiIiIiIpKUbcDpLfZ9\nCDjVXyi5UuR9E0nMVuA0x7bbmmgrUmgjsg4gA92et1cBujxvsx319q/Zffb9byRSKNuAi4B1wHPA\nccDNwOPAb4AP1bQPj1AWAVuAZ4ANwNmhdt8GXgD2ALuBTzB4xPJJ4L9rtn1N8AA4tkEcYXFta/dv\nZJ33RgAnAmVgF3b48qaYbdT7Y7ONwYeC24CPAw8CTwM3AAfH7INIIWwDVgMvBQ4FVgGfxv4Svxz4\nNfD6UPtwQnkrMCF4/nbg2dDr2ra1rycBfwTGBK9HAjuAWdgXtlEcVY3ahvfv4Ij3DsIS46JgG/Ox\nJDklZhu16u3rPdi/x4uAjcAHIvqKFMZW4N3B89nAb2t+vxhYVtM+al5kDfDmmLa1r/uA84LnZ2Bf\natc4cGwb3j8i3jsF2FnT5jrgkpht1Kq3r+8Kvb4S+FqdfhOxkd2NweuPAv/c4LOkAw2nY+Xtwc+X\nYYcPu0K/GwmsiOh3PvYFmBy8HgO8uInPvQ54J3Z49C7guy3E4dJ2O0OF3zu2TpvfBu/HbaORR0PP\nn6vZXtVULBF/LHh9A/ZvIgUznBJKJfi5HfvLOtWhz8uAf8f+Kt8dbGMNgydhK3X6hd0EfAE7lDgb\nmBO8/39NxOHStl4c4fd2YHNHXaH3Xwb0N9hGM6L6/xK4DPhm8PoI4Odtfpbk0HA8y3MvNoF6EXAI\n9pf+JGBmnbaHYV+S32P/Vu8J2oY9Bhwf83lPYBOh38QmUzcF79/XRBzNtI1yDzZ5fBE2n1ICzsRG\nC77Ene16NfZvD/AqbPJXCmY4JpT92BfpZOwL/gQ2Cjm8TtuN2OjibmxofxKwsqbNFdhk6S7sjEe9\nv9LXYWdHrmsxjmbaRnkeO6vzxqD/v2JzO5ub2EYjFaJHKTdgI7S/Be4K3hsBXAD0AG/wGIdkJE/X\nT8jw80/YqOlX2MjtiWzDEZFO9hVgbPD8yCwD6RALsDmvX2HXONXz5eD3D2KHllXLsMPz9TXtjwBu\nx0aqPwPGe4xXJFVzsbM9r8cmjCXaSOySg8nYHNha7ELFsIXAj4Pns7F5s6pTsARTm1CuwubVwJLU\nEm8Ri0huvQa4LfR6UfAI+zpwTuh1P4MvwpzM0ITSDxwdPJ/A4LN+TRuOk7IineilDL5O6JHgvWbb\n1DoaOxQi+Hl0TNuGlFBEOoPrNUK1J1qaubYo7iydk7YTSs/JJ1eD0EOPwj7mzu2p0KYXNfeZz9R0\n/x2D55mOw0YgcW0mBu/FeYyBw6JjsJtPW9b2lbJ3rl1L5a67Gjd00Lt0Kb3ve1/7Gzqp9tqz1vV+\n7nP0Xnyxn42NHu1lM72XXUbvZz7jZVtPPnOQl+0AXHllL5/8ZK+XbXV7uob7iit6Wby4t+3tjB/f\n/hUWu4DKCLe/4V3794+teesB7EbOydhVz+cw9PaF5cCF2DU/c7A7wB8j3nLsWqArg5+3OAUYQYc8\nImkaMcLtMdQ+LFn8FLvg8nvAw9jd3dU7vH+MXfi4BbgW+GCo//XYBYVTsXmW9wTvL8FuWt2M3WLS\n1lme4XQvj0j2HEcoEX4SPMKurXl9YUTfqJsxnwJe105QYblKKKUZM7IOYYjSKadkHcIQpZ6erEOo\na+7cUtYhDDFvXinrEAZzPZbbuzfZOBLi49L7iq85FG88zqF45WkOxSefcyg++ZpD8SWYQ2n3+1Kp\njBnTuBXQ9eyzPj4vdTn7zyZScO0d8uSeEopImpRQRMQbJRQR8UYJRUS8UUIREW/ydvrKs2LvnUje\naIQiIt4UPKG47J1L2TkRcdH6vTwdodEIZSRWHf112G3Q92N3Jz6ccFwixdTBycJFo72bhd25uA1b\nhuEG4KyEYxIprmE+QqlXUm52cuGIFFwHJwsXjRKKU5Wq3qVLDzwvzZiRy7uGRZrR11dm5cqy/w23\nd9p4AfAlbCpiKVYUqdaXscXc9gDvxpbOjes7C5vWOAirufJBbGqjJY3uZpwD9AbBACzGVrEL74ju\nNnalu42d5e1yDW93Gzv+v9n10EO1nzcSWwwtPJ/5TgbPZy7E6qEsxI4krsG+w3F9y9jqlz/FEtFF\nwPzmd800Gn+Fy86NwsrOLW/1w0SGvdbnUFzmM98MfCt4fi+2aNeEBn13AuOC5+NpXIM2VqO/A+Gy\ncyOBb6AzPCKta30OxWU+M2oZjWNj+i7C1uv+F2yA8ZpWAwS3C9vqlZ0TkVZEJJTy7t2Ud++O6+la\ndb/Zw7JvAB8GfgC8DVuy9Iwmt3FAzo5URQouIqGUxo2jNG7cgdeXPvpobZNWl9F4BJtwjeo7i4Ga\nsjdhE7YtK/Y5LJG8aX0OxWU+czlwfvA8vIxGXN8tQLVI8WlY9fuWaYQikqbWT19FzWdWl9C4FltG\nYyGWJP7IwFIZcXOh7wf+DTgYeC543TIVqU6TThs7K+xp43nznBp2rVzp4/NSl7P/bCIFN8yvlBUR\nn5RQRMQbJRQR8UYJRUS8ydtss2fF3juRvNEIRUS8UUIREW+UUETEGyUUB3m7AnTUqKwjqKvSnb+r\nUg8/POsI6ivs3KUSioh4o4QiIt4Uduhlir13InlT8BFKsfdOJG/aW5fHZRXPLwe/fxB4lWPfD2Hl\nDB6ifiV9ZxqhiKSp9RGKyyqeC4G/xIopzQa+xkDV+6i+87Hi1tOxAtZHtRogaIQikq78Vb3/R2wZ\njeeD10+0tXvtdBaRJrWeUKIq2ru0qVf1vtp3CnAqcA+2Rs/MFvcM0CGPSLqiqt7v2EF55864nklV\nve8GXoQdGv01cCPwF01uY9DGRCQtEaeNS5MmUZo06cDrS9esqW2SVNX7R4DvB8/vx1YGfTHwZOx+\nRNAhj0ia8lf1/has2j3A1OD3LSUT0AhFJF2tn+VJqur9suCxHtjLQEJqiRKKSJrau7Ct3iqe19a8\nvrCJvmBnd85rJ6gwJRSRNBX8SlklFJE0FTyhuOzdMmxiZ33CsYgUX3uX3ueeS+T/id0HICLt6u52\ne3Qol8j7sNNNItKuDh59uOjcVCjSiZRQGuv9+tcPPC/NnElpZlu3A4hkrlwuUy6X/W+44AnF9br/\nycAPgVfW+V2lsnq1t4C8mDYt6wjqqow6OOsQhti3L+sI6svbNMKIEV3Q/H0ytSqV3l6nhl3Wrt3P\nS13O/rOJFFzBRygue3c9cBd2nf92Bi7nFZFmFfy0scsI5Z2JRyEyXOTtWM6zYu+dSN508OjDhRKK\nSJoKnlCKvXcieZPPqvcAH8eKKx3R/E4NUEIRSVPrCaVauX4BMA2b2zyxpk246v37sar3Ln2PA84A\nftv27rW7ARFpQv6q3gN8EbjIx+5pDkUkTa2f5alX0X62Q5uoqvfVvmcFr9e1GliYEopImlqflE2i\n6v0hwMXY4U4r/YdQQhFJU9QyGps2Ud60Ka5nElXvj8duq3kw1H4Vdoj0eFwwUZRQRNIUkVBKJ55I\n6cSBedJLf/Sj2ibhyvU7sMr1tRedLsdqyt7A4Kr3T0b0fRg4OtR/K/Bq4KnmdmqAEopImvJX9T7M\n9bAqkhKKSJryV/U+rOUVA6uUUETSVPArZZVQRNKkmwNFxBuNUBwcfriXzXiT078CeayOtn9/1hHU\nt3dv1hEkRAlFRLxRQhERb5RQRMQbJRQR8UYJRUS8yekJA1+KvXcieaMRioh4o4QiIt4ooYiINwVP\nKMXeO5G8yV/V+89jpQweBL4PjGt110AJRSRd+at6/zPgFcBfAZuBxW3tXjudRaRJ3d1uj6GSqnp/\nO7YeT7XPxHZ2TwlFJE2tj1CiKtq7tKlX9b62L8DfY1XfWuaSUI4D7gA2AA8BH27nA0WGtdYTShJV\n78M+BewFrmuxP+B2lud54KPAWmAMVhX7durXpBSROFFV7x94gPKqVXE9k6h6X/VubP7l9LgAXLSS\nzW4BvgL8InhdqWzZ0m4cfk2enHUEdT2/f2TWIQyR13ooeTN6dBe0uWYNUKmsXu3UsGvGjNrP6wY2\nYV/6HcB9DFSur1qI1ZRdiFW9/1LwM67vAuALQA/w+5b2qibIZkzGTkXd2+4HiwxL+at6/xVgFHbU\nAXA38MFWg2wm444BysDl2CilSiMURxqhdC5vI5R1bit+dk2f7uPzUuc6QjkIuBn4DoOTCQC911xz\n4Hlp9mxKc+Z4CU4kK3feWWbFirL/DRf8bmOXDNiFndt+EpucraURiiONUDqXtxFKf79Tw64TTvDx\nealzOaCbC5wLzAfWBI8FSQYlUljtXXqfey7jr5XoAjgRPzo4Wbgo9gGdSN4ooYiIN0ooIuJNwc/y\nFHvvRPJGIxQR8UYJRUS8UUIREW+UUETEGyUUEfFGCUVEvCn4aeNip0uRvMnfMhpHYLVQNmMV8Me3\numughCKSrvwto7EISyhTsSqMi9ravXY6i0iTWk8oSS2jEe7zLeDstnavnc4i0qT8LaNxNPBY8Pyx\n4HXLvMwQ/Xni8T42482onKbJ/fuyjmAoFVhKWVTV+3KZcrkc19PnMhpdEdurNPE5dRV7ylkkZyoR\n3/ee0nx6SvMPvL70sstqm/hcRmNi0BZsVDIBeBQ4Bni88V5Ey+nfcpFi2rfP7VHHA9hk62SsSv05\nwPKaNsuB84Pnc4CnsYQR13c5cEHw/ALq1IxuhkYoIilq4xAzqWU0lgA3Au/FJm3f3nKE+CmCW/nT\nn9o67PJu1KisI6hv796sIxhKcyhuDj3UT5HqPXvcviuePi91GqGIpKjoCVwJRSRFSigi4o0Sioh4\no4QiIt5EnBIuDCUUkRRphCIi3iihiIg3Sigi4o0SCowG7gQOxu4D+B9gcZJBiRSVEgr8CZgP7Ana\nrwTmBT9FpAlKKGZP8HMUdnPRU8mEI1JsRT9t7Fq+YASwFrsV+g5gY2IRiRTY/v1uj07lmlD2Aydj\nhVlOBUpJBSRSZAklFNfK9c1Wvj8Dq6WyLvg5nwaaPcvzB+BWYCZQrr752c/2Hmhw6qklenpKTW5W\nJF9WrCizYkXZ+3YTGn1UK9dfhSWKRQytXl+tfP86rFrb/VhxpYdj+j8BnIlVc3sFVk9lYlwgLvUW\njsQKtDwNHBJs9FKs5D6oHooz1UPpXL7qoaxb5/ZdmT69qc/rB3oYKOdYBk6oafMa4BJslAIDCWeJ\nY/8u4PfB75+PCsRlhHIMVl5/RPD4NgPJRESakFACd6lcX68i/uwm+v8dsIqYZAJuCWU9MMOhnYg0\n0MZZntux0UGtT9W8jqpcX/teM5XvX4GNZM5oFKSulBVJUdQIZdWqMqtWleO6xn2ZXSrX16uI71L5\nfiLwfeA8YGtcgKCasqnSHErn8jWHcs89bt+VOXOa+ryrgCeBK7G5kfEMnZTtBjYBpwM7gPuwJUkf\njuk/HrtK/hIcq+FrGQ2RFCV02rh6OLIZOC14DbZi4K3B83Dl+43A9xhc+b5e/wuB47GEsiZ4HBkX\niEYoKdIIpXP5GqGsXOn2XZk3T1XvRaSBoidwJRSRFCmhiIg3Rb85UAlFJEUaoYiIN0ooIuKNEoqI\neKOEIiLeKKGIiDdKKCLijU4bOzj3XB9b8WfZsqwjqK9b6dvZCy9kHUEyNEIREW+UUETEm6InFJUv\nEElRh1W9r5oEPAt8vFEgSigiKUoooVSr1k/F6j3XFleCgar3C4BpWHGlEx37f5GBuiqxlFBEUpRQ\nQnkzVkie4OfZddrMArYA27BC0zcAZzn0Pxv4DY6L+ymhiKRo3z63R5NarXr/0gb9xwAXAb2ugWhS\nViRFbUzKZlH1vhe4Glvb3Kl6nBKKSIqiEsqmTWU2by7Hdc2i6v0sbD2eq7CJ2v3Ac8BXowJRQhFJ\nUVRCmTKlxJQppQOvb7310mY2uxy4AKtafwH1K9Q/AEwBJmNV78/BJmbj+p8a6n8JsJuYZAKaQxFJ\nVYdVvW+aRigiKUrowransEXQa+0A/ib0+ifBw7V/mNOQSQlFJEVFv1JWCUUkRbrbWES80QjFjMRm\niR8B3pRcOCLFpoRiPoLNDI9NMBaRwit6QnE5bTwRWAgspQPXWhXJk4ROG+eGywjlauATwOEJxyJS\neJ2cLFw0SihnYpfhrgFKUY02bOg98Pyoo0q85CWRTUU6Ql9fmb6+svftFv0sT6NDmM8B52FX2Y3G\nRik3A+eH2lTe+tZ69xhlRzVlO1/easqOHdsF7R/yVz79abfvyuWXe/m81DWaQ7kYu6Ho5cA7gF8y\nOJmISBM0hzJYvoYiIh2mk5OFi2YSyp3BQ0RapIQiIt4ooYiIN0VPKKqHIpKihGrKJrmMxnTgbuAh\nYB1wcFwgSigiKeqwZTS6gW8D7wdOAnqwivmRlFBEUtRhy2i8HhuVrA9e78LqykZSQhFJUUIJJall\nNKZil4rcBqzCbsGJpUlZkRRFJYsdO8rs2FGO65rFMhrdwDxgJlbt/hdYYvllVJBKKCIpikooEyaU\nmDChdOD16tVDSrhmsYzGdmAFVnMW4MfADGISig55RFKU0CFPdRkMcFtGYxS2jMbyBv1/BrwSOAQb\nfPQAG+IC0QhFJEUJ3W28BLgReC826fr24P1jgf/AKt+Hl9EYCXyDwcto1Ou/C1so/X7sMOhW6lfN\nP0AJRSRFHbiMxneDhxMlFJEUFf1KWSUUkRQpoYiIN0ooDm66ycdW/Ln66qwjqG/cuKwjGKqr42qC\ndTYlFBHxRglFRLwpepFqJRSRFGmEIiLeKKGIiDdKKCLijRKKiHijhCIi3iihiIg3RT9trHooIilK\nqB5KUlXvRwPXY3VlN1K/+PUgSigiKeqwqvfvCH5OB14NfACYFBeIEopIijqs6v1O4DAsGR0G7AWe\niQtECUUkRR1W9f6nWALZiSWizwNPxwXiOim7LdjwC1h2m+XYT0RCopLFM8+U2b27HNc1i6r352L1\nZI/B5ln6sEOirVFBuiaUClBioPq1iLQgKqGMGVNizJjSgdc7d+ai6v1rgR9gA4kngP/FltSITCjN\nHPKocoZImxJa2zipqvf9wGnB88OAOQwUtq7LNaFUgJ8HQf2DYx8RqZHQHMoSbASzGUsAS4L3j8Uq\n1cPgqvcbge8xuOp9vf7XYslnPXAfsAxbND2S6yHPXGxi5ijsWK4fO54K9IaaloKHSOfq6yvT11f2\nvt0Oq3r/Z2wexVkrhzGXAM8CXwheV+rP7WRn+/bGbbKgEpCda+zYLmj/sL9ywglu35X+fi+flzqX\nQ55DgbHB88OwFdnXRzcXkSgJHfLkhsshz9HYTG+1/Xexy3NFpEmdnCxcuCSUrcDJSQciMhwU/eZA\n3W0skiKNUETEGyUUEfFGCUVEvFFCERFvlFBExBslFBHxRqeNRcQbjVBExJuiJxSVgBRJUU6r3r8N\n2IAVUppR02dx0L4fu48vVs4SSjnrAIa4++5y1iEMkcRt9T7kMa68xZTTqvfrgbcAK2r6TMMKMU0L\n+n2VBjlDCaUBJRR3eYwrbzHltOp9Pza6qXUWti7P80G/LTSoJ52zhCJSbDmteh/l2KCdcx9Nyoqk\nqI3Txr6r3rcq8WpqZQZ2Qg89ivoo075mPi92Qa0a/Qwkm2OC17XmALeFXi9m8MQswB0MnpRdxOD5\nmNuA2U3EJSId6CoGksMiBopMh3UDv2ag6v1aBiZlq+7Alhytmha0GwW8POjfcWUpRaQ5R2CrUtSe\nNg5XvQd4I7AJm1xdHHr/Ldj8ynPY2jzhQtYXB+37gTckELuISP5FXXSTlWXYjHmeCnIfhw1LN2Dr\no3w423AAGA3ciw2NNwJXZBvOICOBNcAPsw5E0jUSG1ZNBg6i/vFd2k4BXkW+EsoEBur7jsGGr1n/\nO4GtjAB2nH4PMC/DWMI+hhVVX96oofiRl+tQ4i66yUofsCvjGGo9iiVbsLWRHsaOk7O2J/g5Cvvj\nkIc1sCcCC4GlaCIxNXlJKK1cdDPcTcZGUPdmHAfY/0drsUPEO7BDn6xdDXwCKPjtePmSl4RSyTqA\nDjMGuAn4CDZSydp+7FBsInAq2a9FeybwODZ/otFJivKSUH6HTThWHcfgS35lwEHAzcB3gFsyjqXW\nH7DTlDMzjuO12P0tW7F7UU4D/ivTiCRVLhfdZGEy+ZqU7cK+GFdnHUjIkQxc93AIdsfq6dmFM0QP\nOsszLEVddJOV67HV6/+Mze+8J9twADt7sh9LuGuCx4JMI4JXAquxmNZh8xZ50oPO8oiIiIiIiIiI\niIiIiIiIiIiIiIgMT/8P+HFnfN5/OHMAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10fe93710>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# mean relative error\n",
      "print 'rms relative error in Kt:', sqrt((Kt_error**2).mean())\n",
      "print 'rms relative error in ve:', sqrt((ve_error**2).mean())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "rms relative error in Kt: 0.0149482687766\n",
        "rms relative error in ve: 0.0043530556982\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}