{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating and fitting multiple voxels with varying initial conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme\n",
    "acq_scheme = wu_minn_hcp_acquisition_scheme()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating N-dimensional datasets with varying parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For real applications, we may need to simulate data with many different initial conditions, or fit a real dataset where we want to give a voxel-dependent initial condition to avoid local minima. The toolbox allows you to do both these things by generating a multi-dimensional parameter vector.\n",
    "\n",
    "As an example, we will generate a simple Ball & Stick model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OrderedDict([('G1Ball_1_lambda_iso', 1),\n",
       "             ('C1Stick_1_mu', 2),\n",
       "             ('C1Stick_1_lambda_par', 1),\n",
       "             ('partial_volume_0', 1),\n",
       "             ('partial_volume_1', 1)])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from dmipy.signal_models import cylinder_models, gaussian_models\n",
    "from dmipy.core.modeling_framework import MultiCompartmentModel\n",
    "stick = cylinder_models.C1Stick()\n",
    "ball = gaussian_models.G1Ball()\n",
    "ball_and_stick = MultiCompartmentModel(models=[ball, stick])\n",
    "ball_and_stick.parameter_cardinality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now say we want to simulate the signal in each voxel to have the same lambdas and partial_volume, but varying 10x10 grid of $\\mu$. This can be done by giving parameters_to_parameter vector just the single value you want for the lambdas and partial_volume, but the mu as a 10x10x2 array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ND parameter vector shape: (10, 10, 6)\n",
      "Parameters of first two voxels:\n",
      "[[  9.67584402e-10   6.70524705e-01   1.04730408e-01   4.83792201e-10\n",
      "    5.00000000e-01   5.00000000e-01]\n",
      " [  9.67584402e-10   2.12527713e-01   5.88270976e-01   4.83792201e-10\n",
      "    5.00000000e-01   5.00000000e-01]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "lambda_par = np.random.rand() * 1e-9\n",
    "lambda_iso = lambda_par * 2.\n",
    "partial_volume = 0.5\n",
    "mu_array = np.random.rand(10, 10, 2)\n",
    "\n",
    "gt_parameter_vector = (\n",
    "    ball_and_stick.parameters_to_parameter_vector(\n",
    "        C1Stick_1_lambda_par=lambda_par,\n",
    "        G1Ball_1_lambda_iso=lambda_iso,\n",
    "        C1Stick_1_mu=mu_array,\n",
    "        partial_volume_0=partial_volume,\n",
    "        partial_volume_1=partial_volume)\n",
    ")\n",
    "\n",
    "print 'ND parameter vector shape:', gt_parameter_vector.shape\n",
    "print 'Parameters of first two voxels:'\n",
    "print gt_parameter_vector[0, :2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10, 10, 288)\n"
     ]
    }
   ],
   "source": [
    "# the data can now be generated exactly as before.\n",
    "E_array = ball_and_stick.simulate_signal(acq_scheme, gt_parameter_vector)\n",
    "print E_array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting N-dimensional data with varying initial condition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To give an initial guess for only some parameters we can use a model's \"set_initial_guess_parameter\" function. It allows to avoid doing any kind of brute force optimization for that parameter, but still do local optimization starting from that value.\n",
    "\n",
    "By giving an initial condition, the optimizer will try to find a global optimum of all other parameters *given* the initial guess parameters, and then polish the result for a local minimum. First, to see how this works we can give an initial guess that our Stick model is aligned with some random orientation, and actually the ground truth orientation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using parallel processing with 8 workers.\n",
      "Cannot estimate signal grid with voxel-dependent x0_vector.\n",
      "Setup brute2fine optimizer in 0.000231027603149 seconds\n",
      "Fitting of 100 voxels complete in 1.22214388847 seconds.\n",
      "Average of 0.0122214388847 seconds per voxel.\n",
      "Using parallel processing with 8 workers.\n",
      "Cannot estimate signal grid with voxel-dependent x0_vector.\n",
      "Setup brute2fine optimizer in 0.000200033187866 seconds\n",
      "Fitting of 100 voxels complete in 1.06655097008 seconds.\n",
      "Average of 0.0106655097008 seconds per voxel.\n"
     ]
    }
   ],
   "source": [
    "random_initial_mu = np.random.rand(10, 10, 2)\n",
    "ball_and_stick.set_initial_guess_parameter('C1Stick_1_mu', random_initial_mu)\n",
    "BAS_fit_random = ball_and_stick.fit(acq_scheme, E_array, Ns=5)\n",
    "\n",
    "true_initial_mu = mu_array\n",
    "ball_and_stick.set_initial_guess_parameter('C1Stick_1_mu', true_initial_mu)\n",
    "BAS_fit_true = ball_and_stick.fit(acq_scheme, E_array, Ns=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how only the mu parameters have a value and the rest is None.  \n",
    "We can then fit the data with our standard brute2fine optimizer and sample the remaining parameters at 5 equal points between their respective parameter_ranges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we then visualize the fitting error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f81efc46290>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4kAAAD5CAYAAACK5HAJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xm8XWV97/HPlwxMMiZVgaBQidU4\nAJqiLa0DqARrAW+lBhywF0sHqAP1FvBexUvFSmuBWnCISqGoBKRQYosiKK1yZQqDQEAkBoQAijFh\nkDnn/O4f6znJys4e1j55zt5Pcr7v12u/cvYafs+zd06+Wc+z11pbEYGZmZmZmZkZwGbD7oCZmZmZ\nmZmVw4NEMzMzMzMzW8ODRDMzMzMzM1vDg0QzMzMzMzNbw4NEMzMzMzMzW8ODRDMzMzMzM1vDg0Qz\nM7NJRNJZkh6SdFumen8vaYmkOyR9VpJy1DUzs+HxINHMzGxyORuYl6OQpN8F9gVeCbwc+G3g9Tlq\nm5nZ8HiQaGZmNolExPeBlfVlkl4k6duSbpD0A0kvaVoO2AKYDmwOTAN+kbXDZmY2cB4kmpmZ2QLg\nryLi1cBHgM812SkirgauBB5Mj8si4o4J66WZmQ3E1GF3wMzMzIZH0nOA3wW+UbuccPO07n8AJ7XZ\n7f6IOEDSHsBLgVlp+eWSXpc+rTQzs42UB4nWlaTdgLuBaRGxeri9MTOzCbAZ8HBE7NW6IiIuAi7q\nsu/bgWsi4tcAkr4FvBbwINEGzscsZvn4dNMNJOkeSc9Imtmy/GZJkQILSbMk/ZukFZIekXSrpPel\ndbulbX/d8nhnhzb/S9JTaZsVki6StNMEv1Qz24QNOsvS3TDH1o/UMu3Xkj46gJdsSUQ8Ctwt6VAA\nVfZsuPu9wOslTZU0jeqmNT7dtFA+ZhkcSWen9+mgluWnp+XvS8+nS/pHScvTe3S3pNNq298j6cmW\n9/qMNu19obb+GUnP1p5/a8JfsG1yPEjM427gsLEnkl4BbNmyzbnAfcALgRnAe1n/4v7tI+I5tcf5\nXdo8JiKeA+wBPAf4zAa+BjOzgWVZRLxsbD3wA1KmpcenWreX5DNfMpF0HnA18FvpwPRI4F3AkZJ+\nBCwBDm5Y7kLgp8CtwI+AH0XENyeg25aPj1kG5yfAEWNPUo4dSvVvZswJwFxgH2Ab4I3ATS11/rDl\nvT6mtaGI+PNapn4KOL+2/YGt2ztTrRcPEvM4lypAxxwB/GvLNr8NnB0Rj0fE6oi4KSI2eGYnIh4G\n/h1Yc5qQpH0kXS3pYUkPSjpD0vTa+pD055LukrRK0plSdSGKpCmSPpNm+5YBf1BvT9LOkhZJWilp\nqaQ/ra37hKRvSPqqpMfSzOOLJZ2g6ju57pP0lg19zWY2YYaWZa0kvV/S91V9795K4P9I+qSks2vb\n7CEpas+3l/QvKfeWSzpJkv+faxERh0XEThExLSJmRcRXIuLuiJgXEXtGxJyIaHcdYrtaIxHxZxHx\n0rTfsRPdf9tgPmZhYMcs3wT2lbRDej4PuAX4eW2b3wYujogHonJPRLT+fWywsbyU9CeS7gW+I+lN\nku5p2W65pDeknzeT9FFJP03v8cLaa7FNnP/zzOMaYFtJL5U0BXgn8NU225wpab6kF+RqWNIM4H8A\nS2uLR4APAzOB3wH2B/6yZde3UQXTnsAfAwek5X+a1u1NNbP1jpb9zgOWAzundZ+StH9t/R9S/Qe0\nA9VM2GVUv2e7UN384Ivje6VmNgBDy7IOfpfq1MXfAE5psP1XgSeBF1Hl1x8AfzJhvTPbOPmYZa2J\nPmZ5ClgEzE/P38v6A/JrgGMl/aWkV4wNgCfQ64CX0DKg7uDYtN3rqG5O9Tjw2YnrmpXEg8R8xmbm\n3gz8GLi/Zf2hVKdUfYzq2o+bJf12yzYr0kza2OOlXdr7rKRHgBVUwfpXYysi4oaIuCbN/t1DFXKt\nX2786Yh4OCLupbp9+dis3h8Dp0fEfRGxEvi7sR0k7Qr8HnBcRDwVETcDXwbeU6v7g4i4LF0w/g2q\ng7tPR8SzwEJgN0nbd3ldZjZcg86ybu6NiM+nT6ue7LahpF2oDi4/HBFPRMTPgdNZe3BmZmv5mKUy\niGOWfwXeK2m79Lr+vWX931FNgr0LWAzcL+mIlm3+veW9/lPG78SUkV0zNfkz4KMRcX9EPAV8Avhj\nn6ExOfh85HzOpbqb2+6sP0tERKwCjgeOV3XB+Geo/tHPqm02s4+7cX0gIr6s6lqC/6Ca4bkXQNKL\ngVOpZtW2ovp7vqFl//qpDk9QXSMA1WzbfbV1P6v9vDOwMiIea1k/t/a8fs3Ck8CKiBipPSe19XCv\nFziZHfDGreNXK0d6bwjccMvTl0XEvAnukk0eg86ybu7rvckaL6T62oZf1CbiNwPuydAPS5xNmwwf\ns1TGdcwi6V2s/ZTxB+2u+RsTEVdJ+g3g/wD/ERFP1j8sTO2dSfXJ7ZbA/wTOknRdrP3O0UMi4opO\nbfSpn1x9AfBNSaO1ZQE8l3X/TorgfMrLg8RMIuJnku4G3goc2WPbFZI+Q3UdwI4b2O6tkj5JFS6v\niogAPk912sRhEfGYpA+x/ikYnTwI7Fp7Xj/N5AFgR0nb1EL3Baw/A2kbaMXKEa69bFbvDYFpO/10\nZu+tzJoZVpZ1aqLl+eNUB5Fjnl/7+T6qg8cdI2IUmxDOpk2Dj1k2TER8DfhaH7t8Ffg41U1putV9\nkuq9+b/AHCbgTsHpPR+zTqaqupnNjNr65cDhEXFt7n5MBOdTXv64OK8jgf0i4vHWFZJOkfRyVbcJ\n3wb4C2BpRPwqQ7vnUM3qjN1meRvgUeDXkl6S2mrqAuADqm5/vQPVTCIAEXEf8EPg7yRtIemVVK+5\nn6C0RoKRGG30MJsAw8qyXm6m+rqFXdMpYK359N/AZyRtm264sIek1w2gX5OIs2kT4mOWwfks1am9\n631/qKQPSXqDpC3T+30E1XvSeofTifBjYBtJB6j6CpsTgWm19V+guo7zBamvz1XLV3qUxfmUkweJ\nGUXETyNicYfVWwEXU52ysIzq1KjWf2gPa93vwWl0l7iIeIYqgD6WFn0EOBx4DPgS0O221K2+RHXh\n9o+AG1n/S5QPA3ajmqG7mOrc9sv7qG8NBDBKNHqY5TasLGvg26ntW4HrqG4IUfduYGvgdmAV1TVG\nz8eycTZtOnzMMjgRsTIivtvyKd6YJ4F/pDp9cwVwNPBHEbGsts03W97rizP1axXV9aHnUH3CupJ1\nTyM9lSp3vyvpMapBd+u1qcVwPuWl9r+vZjZMe+85Pf77W82Obbfb5b4bImJu7y3NzDaMs8nMSuV8\nysvXJJoVyjNdZlYiZ5OZlcr5lI8HiWYFCmDEQWdmhXE2mVmpnE95eZBoVijPhplZiZxNZlYq51M+\nHiSaFSiAEV8vbGaFcTaZWamcT3l5kGhWoCB41rNhZlYYZ5OZlcr5lNeEDBKna/PYgq0norRZRy9+\n5RPZav3klq16b9SHx1i1IiJ+o/EOASPOueymbL11TNs+33e+b/7LZ7LViul54/jpHfJ+w9Hm96/3\nVWrjpunTs9UC4NnVeevtMSVrufjJs9lqaWre35NHV//S2VQAHzfZMGxSx03gfMpsQgaJW7A1r9H+\nE1HarKPLLrs5W60Ddt4rWy2AK+LCn/WzffVdP5bbtO13ZNejP5yt3h5fuDdbrdU75xu8Aiydn/eA\nc/bfdPo6tf5NmTUrWy2A0Z8/lLWeFmyftd7I/r/IVmvKjnl/Ty576PPOpgJswda8ZrM35Suogr8G\nO/cXmZf8WnMbHclaLutx06xXZ6sFcMXI+X1lEzifcvPppmZFEiNo2J0wM2vhbDKzUjmfcvIg0axA\nAYz6lAkzK4yzycxK5XzKy4NEs0J5NszMSuRsMrNSOZ/y8SDRrEABPBuT6DoLM9soOJvMrFTOp7wa\nvZOS5km6U9JSScdPdKfMJrugmg1r8pjMnE1mg+Vsas75ZDZYzqe8en6SKGkKcCbwZmA5cL2kRRFx\n+0R3zmyyCsRIszmcScvZZDZ4zqZmnE9mg+d8yqvJO7kPsDQilkXEM8BC4OCJ7ZaZjYYaPSYxZ5PZ\nEDibGnE+mQ2B8ymfJtck7gLcV3u+HHjNxHTHzGDtKRPWlbPJbMCcTY05n8wGzPmUV5NPEtu92+vd\nYFbSUZIWS1r8LE9veM/MJrFAPBtTGz166XVdjKTNJZ2f1l8rabfauhPS8jslHdCrpqTdU427Us3p\n3dqQNEPSlZJ+LemMWp1tJN1ce6yQdHpa9z5JvwT+AXi7pPev89at//rWZNPI44/3fL/MrLOc2bSJ\n63ns5OMms7x87NT72Km2rn7s1FaTQeJyYNfa81nAA60bRcSCiJgbEXOnsXmDsmbWTY6Lr2vXxRwI\nzAEOkzSnZbMjgVURsQdwGnBK2ncOMB94GTAP+JykKT1qngKcFhGzgVWpdsc2gKeAjwEfqXcoIh6L\niL3GHsDPgItqm5wPHA4sjogvp2U9s2nK1lt3fb/MrLdcN4aQdJakhyTd1mH9/6od0NwmaUTSjmnd\nPZJuTesWZ36JOfQ8dvJxk1l+PnbqfuxUW/9lemgySLwemJ1GudPTC1/UYD8zG6cIMRKbNXr00OS6\nmIOBc9LPFwL7S1JavjAino6Iu4GlqV7bmmmf/VINUs1DurUREY9HxFVUgdeWpNnAc4EftKxyNpkN\nWMZsAjib6iCqQ1vxD7WDnROA/46IlbVN3pjWz92gFzUxnE9mA+Zjp7W6HDs11vNdiojVwDHAZcAd\nwAURsWS8DZpZM6Oo0QOYOXbKUnocVSvT7rqYXVqaWrNN+vf+CDCjy76dls8AHk41Wtvq1EYTh1HN\nftVP1foj4Eaqmfnv4mwyG5g+sqmriPg+sLLnhpXDgPM2pN+D5GMns+HwsdMabY+dJN0i6UJJu3ba\ncUyjiwYi4lLg0oadMrMNVF183fg2ziu6zKQ3uaa40zadlrfrWLftm/ajk/nAe2rPvwmcFxFPS/pz\n4I8jYr+GtcxsA/SZTVlI2orqE8djWrryHUkBfDEiFgy0Uw342MlssHzstI5ex07nUH2C2dGkv7Lc\nrExqerpWL02uKR7bZrmkqcB2VLP73fZtt3wFsL2kqWnGq759pza6krQnMDUibhhbFhG/qm3yJdae\no29mE66vbJrZcr3ggnEO5v4Q+H8tp5ruGxEPSHoucLmkH6dPJs1s0vKxE+Q7dvI3TpoVKIBnY0qj\nRw9NrotZBByRfn4H8L10esIiYH66u9buwGzguk410z5Xphqkmpf0aKOX9U4xk7RT7elBVKdymdkA\n9JlNK8ZuzJIe4/20bz4tORARD6Q/HwIuprrex8wmMR87rZHl2MmfJJoVKFCWU7oiYrWksetipgBn\nRcQSSSdR3Rl0EfAV4FxJS6lmqOanfZdIugC4HVgNHB0RIwDtaqYmjwMWSvokcFOqTac2Uq17gG2B\n6ZIOAd4SEben1X8MvLXlZX1A0kGpTyuB923g22RmDeXKpqYkbQe8Hnh3bdnWwGYR8Vj6+S3ASQPr\nlJkVycdOeY+dPEg0K9RonlMm2l4XExEfr/38FHBoh31PBk5uUjMtX0abGf0ebezWpe+/2WbZCVR3\nOjSzIciVTZLOA95AdVrqcuBEYBpARHwhbfZ24DsRUf+S0+cBF1c3BWQq8PWI+HaWTpnZRs3HTvmO\nnTxINCvQMG4OYWbWS85siojDGmxzNtVXZdSXLQP2zNIJM9tk+NgpLw8SbZNxwM57DbsL2QRiJHrf\nQt76E5uPsvqFHb9WqG/3vPeF2Wo9+ZJ8/QJ46Qn3Zq238//L92XfX9r137PVAjjwt34/a73RNz+U\ntd7qN+bLJl3V9nvnB8bZtJGI0WH3oDNlPogfHclbL6fNel77NtR6B8x6ddZ6w+Z8ysuDRLNCjXo2\nzMwK5Gwys1I5n/LxINGsQBFqcvctM7OBcjaZWamcT3l5kGhWoIBc3/VjZpaNs8nMSuV8ysuDRLNC\n+eJrMyuRs8nMSuV8yseDRLMCBWLUF1+bWWGcTWZWKudTXh4kmhXKs2FmViJnk5mVyvmUjweJZgUK\nfPG1mZXH2WRmpXI+5eVBolmBAhj1xddmVhhnk5mVyvmUlweJZoUawefVm1l5nE1mVirnUz4eJJoV\nKEKeDTOz4jibzKxUzqe8PEg0K5S/68fMSuRsMrNSOZ/y8SDRrEABjPqUCTMrjLPJzErlfMrLg0Sz\nAgXi2VHfocvMyuJsMrNSOZ/y8iDRrFD+rh8zK5GzycxK5XzKx4NEswIFYjR8yoSZlcXZZGalcj7l\n5UGiWaFGPRtmZgVyNplZqZxP+fidNCtQBIyEGj3MzAbF2WRmpcqZT5LmSbpT0lJJx7dZv7mk89P6\nayXtVlt3Qlp+p6QDetWUtHuqcVeqOb1bG5JmSLpS0q8lndHSr/9KbdycHs/t1d9OPEg0K9RoqNGj\nFwedmeWUK5vMzHLLkU+SpgBnAgcCc4DDJM1p2exIYFVE7AGcBpyS9p0DzAdeBswDPidpSo+apwCn\nRcRsYFWq3bEN4CngY8BHOryEd0XEXunxUI9aHfl0U7MCBeLZ2PA7dNVC6c3AcuB6SYsi4vbaZmuC\nQ9J8quB4Z0vQ7QxcIenFaZ9ONceCbqGkL6Tan+/UBmuD7uXp0epdEbG4ZVmnWj1NeXwztrluyyab\nNvLcxY9nqzXl1Duz1QIYGRnNWu/+t22Xrda8h/fJVgtg2d+2+9UZv1f9ft6/i322uzJbratWvihb\nLQB+v7/Nc2UTgKSzgLcBD0XEen+Jkt4AXALcnRZdFBEnpXXzgH8CpgBfjohPZ+nUpkKT6DOAzQq+\nm2XkzWEi8tbbxGTMp32ApRGxDEDSQuBgoH7sdDDwifTzhcAZkpSWL4yIp4G7JS1N9WhXU9IdwH7A\n4Wmbc1Ldz3dqIyIeB66StEcfr6lTrY6/VJMoRcw2HkG22fo1QRcRzwBjQVd3MFUoQRUc+7cGXUTc\nDYwFXduaaZ/9Ug1SzUO6tRERj0fEVVSDxaY69dfMJljGbAI4m2qmvZsf1GbExwaITWb5zWyS6TOf\nZkpaXHscVSu1C3Bf7fnytIx220TEauARYEaXfTstnwE8nGq0ttWpjV7+JZ2B9bHa8VHftfxJolmR\nxGg0nsOZKan+aduCiFiQfm4XSq9p2X+d4JBUD7prWvYdC652NRsHXa2NFT1e279IGgH+DfhkmvEa\nby0z22B9ZVNXEfH9cZ4u3mSW38wmnb7yaUVEzO1YaH2tn7h12qbT8nYd67Z90360eldE3C9pG6pj\np/cA/zqeWv4k0axQo6jRgxR0tceCWpmJCLp+lzftR6t3RcQrqE6I+32qoBtvLTPLpI9s6jZT39Tv\nSPqRpG9Jella1mSW38wmoT7yqZvlwK6157OABzptI2kqsB2wssu+nZavALZPNVrb6tRGRxFxf/rz\nMeDrrD3Vte9aHiSaFSjjHbocdGaWTZ/Z1G0Cq4kbgRdGxJ7APwP/npZ7osjM1pPx2Ol6YHa6Gd90\nqvszLGrZZhFwRPr5HcD30tlOi4D56SZ7uwOzges61Uz7XJlqkGpe0qONtiRNlTQz/TyN6prv28ZT\nC3y6qVmRArF6NMvF12tCCbifKpQOb9lmLDiuphYckhYBX5d0KtWNa8aCTu1qpn3Ggm4h7YNunTY6\ndToN/raPiBW1oLtiPLXMLJ+M2dS7rYhHaz9fKulz6QCoyeSXmU0yufIpXcpyDHAZ1c2xzoqIJZJO\nAhZHxCLgK8C56cY0K6mOhUjbXUB1+vtq4OiIGAFoVzM1eRywUNIngZtSbTq1kWrdA2wLTJd0CPAW\n4GfAZem4aQrVcdOXetXqxINEs0I1OB2iJwedmeWWI5uakPR84BdpAmofqrOffgU8TO/JLzObhHLl\nU0RcClzasuzjtZ+fAg7tsO/JwMlNaqbly1h7tlR9ebc2duvQ9Vd32L5jrU56DhIl7Up1wePzgVGq\nm2L8Uz+NmFl/xu7QlaXWJhp0tWz6T0nOJrMByJlNks4D3kB17eJy4ERgGkBEfIHqTIG/kLQaeBKY\nn84aaDv5laVTmfjYyWzwcuaTNfskcTXw1xFxY7pTzg2SLm/5njUzyyzXHQQ3Yc4msyHIeHfTw3qs\nPwM4o8O6thNVBXE+mQ2Bj53y6TlIjIgHgQfTz4+lL33cBd9q2mziNP+esUnL2WQ2BM6mRpxPZkPg\nfMqqr2sS0/cZ7Q1c22bdUcBRAFuwVYaumU1eweCu+9kUNM2madvsMNB+mW1qnE3965RPPm4yy8v5\nlFfjQaKk51B9KeOH6nccG5Nubb0AYFvt6DsNmm2AAFaP+pSJJvrJpq2et6uzyWwDOJv60y2ffNxk\nlpfzKa9Gg8R0h8F/A74WERdNbJfMDHzxdRPOJrPBczY143wyGzznUz5N7m4qqlvO3xERp058l8ws\n8Hn1vTibzAbP2dSM88ls8JxPeTX5THZf4D3AfpJuTo+3TnC/zCa9UdToMYk5m8yGwNnUiPPJbAic\nT/k0ubvpVeB302ygwqdM9OJsMhsCZ1MjziezIXA+ZdXX3U3NbDD8hbBmViJnk5mVyvmUlweJZgUK\n5Dt0mVlxnE1mVirnU14eJJoVKjwbZmYFcjaZWamcT/l4kGhWKF9YbWYlcjaZWamcT/l4kGhWoPDF\n1xNi6kOP87wzrs5WT1OnZauV/Zu0N8v7+zP66KO9N2po1WGvzlYLYI9PLcla78HXvSRrvcsufSRj\ntZUZa/XP2TSB5NPkhi5G89bL/ndaev8yGul/F+dTXh4kmhXKp0yYWYmcTWZWKudTPh4kmhVJjPji\nazMrjrPJzErlfMrJg0SzAvk2zmZWImeTmZXK+ZSXB4lmJYrq3Hozs6I4m8ysVM6nrPyZrFmhRlGj\nh5nZIDmbzKxUufJJ0jxJd0paKun4Nus3l3R+Wn+tpN1q605Iy++UdECvmpJ2TzXuSjWnd2tD0gxJ\nV0r6taQzanW2kvSfkn4saYmkT9fWvU/SLyXdnB7v7/UeeJBoVqCguvi6yaMXB52Z5ZIzm8zMcsqV\nT5KmAGcCBwJzgMMkzWnZ7EhgVUTsAZwGnJL2nQPMB14GzAM+J2lKj5qnAKdFxGxgVardsQ3gKeBj\nwEfadP8zEfESYG9gX0kH1tadHxF7pceXu74JeJBoVigxGs0eXas46MwsqzzZBCDpLEkPSbqtw/p3\nSbolPX4oac/aunsk3ZomihZnfIFmttHKlk/7AEsjYllEPAMsBA5u2eZg4Jz084XA/pKUli+MiKcj\n4m5gaarXtmbaZ79Ug1TzkG5tRMTjEXEV1THUGhHxRERcmX5+BrgRmNX7fWvPg0SzQo2OqtGjBwed\nmWWVKZsAzqaagOrkbuD1EfFK4G+BBS3r35gmiuaO64WY2SYnUz7tAtxXe748LWu7TUSsBh4BZnTZ\nt9PyGcDDqUZrW53a6EnS9sAfAt+tLf6jNOl2oaRde9XwINGsQBHZTuly0JlZNhmziYj4PrCyy/of\nRsSq9PQaPFFkZl30mU8zJS2uPY6qlWoXYK23xOm0Ta7lTfuxHklTgfOAz0bEsrT4m8BuadLtCtZO\n3Hfku5uaFaqP2zjPbDndakFEjM24T0TQtZtcGnTQnRcRT0v6c6qg269XLTPLY0i3mD8S+FbteQDf\nkRTAF2uZZ2aTWB/5tKLLWQjLgfoE9CzggQ7bLE/HKttRTXp127fd8hXA9pKmpkn0+vad2uhlAXBX\nRJw+tiAiflVb/yXWXvbTkT9JNCtURLMHKehqj/rBUj9BR8Og67R8TdC1aatTG720DbqIeDo9/RLw\n6gZ1zCyTPrKp20x9Y5LeSDVIPK62eN+IeBXVtdFHS3rdhr4uM9v49ZFP3VwPzE4345tOdX+GRS3b\nLAKOSD+/A/heRERaPj/dsG93YDZwXaeaaZ8rUw1SzUt6tNGRpE9SHWN9qGX5TrWnBwF39HgP/Emi\nWaky3R1wTSgB91OF0uEt24yF0NXUQkjSIuDrkk4FdmZt0KldzbTPWNAtpH3QrdNGt47Xgu79Lct3\niogH09NGQWdm+fSRTd1m6huR9Ergy8CB9ZnwiHgg/fmQpIuprpX+/oa0ZWYbvxzHThGxWtIxwGXA\nFOCsiFgi6SRgcUQsAr4CnCtpKdWk9/y07xJJFwC3A6uBoyNiBKBdzdTkccDCdNxzU6pNpzZSrXuA\nbYHpkg4B3gI8Cvxv4MfAjdWtIjgj3eDvA5IOSn1aCbyv1/vgQaJZgYJmdwfsWcdBZ2YZ5cqmJiS9\nALgIeE9E/KS2fGtgs4h4LP38FuCkgXTKzIqVM58i4lLg0pZlH6/9/BRwaId9TwZOblIzLV9GNdHV\nurxbG7t16HrbNyAiTgBO6LBPWx4kmpUosn2S6KAzs3wyZpOk84A3UJ2Wuhw4EZgGEBFfAD5OdYOr\nz6WJotXpk8nnARenZVOBr0fEt7N0ysw2XhnzyTxINCtX73PmzcwGL1M2RcRhPda/n5bTzdPyZcCe\n6+9hZpOej52y8SDRrFCeDTOzEjmbzKxUzqd8PEg0K1SDu2+ZmQ2cs8nMSuV8yseDRLMCBZ4NM7Py\nOJvMrFTOp7w8SDQrUUCMOujMrDDOJjMrlfMpKw8SzUrlUyaK9/AlL8hWa/qCHbPVAtjykuuz1puy\n3bbZaq14Vd5f7pnf2TJrvae33yxrvS0y1orfeUXGasBV54+jE3m7YBNgdCRfLWU+6Fbef182ftos\n399tjBYSDIV0Y1PgQaJZkeRTJsysQM4mMyuV8yknDxLNSuXZMDMrkbPJzErlfMrGg0SzEvkLYc2s\nRM4mMyuV8ykrDxLNSuXZMDMrkbPJzErlfMrGg0SzUnk2zMxK5Gwys1I5n7JpfIspSVMk3STpPyay\nQ2aWRMPHJOdsMhswZ1NjziezAXM+ZdPPJ4kfBO4A8t0H3czaCzwb1pyzyWxQnE39cj6ZDYrzKatG\nnyRKmgX8AfDlie2OmY2JaPaYzJxNZoPnbGrG+WQ2eM6nfJqebno68DfAaKcNJB0labGkxc/ydJbO\nmU1qPmWiCWeT2aA5m5rqmk/OJrMJ4HzKpucgUdLbgIci4oZu20XEgoiYGxFzp7F5tg6aTVYaVaPH\nZOVsMhsOZ1NvTfLJ2WSWn/PVsb47AAAgAElEQVQpnybXJO4LHCTprcAWwLaSvhoR757YrplNYp7p\nasLZZDZozqamnE9mg+Z8yqrnJ4kRcUJEzIqI3YD5wPcccmYTTdXF100ek5SzyWwYnE1NOJ/MhsH5\nlFPjr8AwswHzefVmViJnk5mVKlM+SZon6U5JSyUd32b95pLOT+uvlbRbbd0Jafmdkg7oVVPS7qnG\nXanm9G5tSJoh6UpJv5Z0Rku/Xi3p1rTPZyUpLd9R0uWpjcsl7dDrPehrkBgR/xURb+tnHzMbJwdd\n46BzNpkNkAeJfXE+mQ1QhnySNAU4EzgQmAMcJmlOy2ZHAqsiYg/gNOCUtO8cqrMHXgbMAz6Xvi+1\nW81TgNMiYjawKtXu2AbwFPAx4CNtuv954ChgdnrMS8uPB76b2vhuet6VP0k0K5WDLlvQmVlG+Saw\nzpL0kKTbOqxXmiBaKukWSa+qrTsiTRTdJemIDX5NZrZpyJNP+wBLI2JZRDwDLAQObtnmYOCc9POF\nwP5pMvtgYGFEPB0RdwNLU722NdM++6UapJqHdGsjIh6PiKuojqHWkLQTsG1EXB0RAfxrh1r1Njry\nINGsRJHtDl0OOjPLJ182AZzN2smfdg5k7STRUVQTR0jaETgReA1VHp3Y5NQpM9vE9ZdPM5W+giY9\njqpV2gW4r/Z8eVpGu20iYjXwCDCjy76dls8AHk41Wtvq1EYnu6T92/X7eRHxYKr1IPDcLnWAZnc3\nNbNhaH661kxJi2vPF0TEgvRzu1B6Tcv+64SQpHrQXdOy71jYtKvZOOhqbazo8JoaB52knkFnZhll\nOpU0Ir5fP729jYOBf00TRddI2j5NIL0BuDwiVgJIupxqsHlenp6Z2UareT6tiIi5Hda1m+Vqrdxp\nm07L230w1237pv1o0qdx8SDRbOPnoDOzEnWbwGqi3xl5M7MclgO71p7PAh7osM1ySVOB7YCVPfZt\nt3wFsL2kqWmSvb59pza69XtWh7Z/IWmnNLm+E/BQlzqAB4lmxVKeIZGDrkYSm22e70urf3Xrb2Sr\n9aL/7Pid2+OSe0Q98uivs9Xa48PXZqsFMDJ1WtZ6M364ZdZ69x/d+uH9+D3/6kez1RqvPrKp2wRW\no6baLOs1IWUAyniLf02iK5NKf62Z+xcjI/mKFfLeZTp2uh6YLWl34H6q+zMc3rLNIuAI4GrgHVRf\ncxOSFgFfl3QqsDPV6fLXUeXWejXTPlemGgtTzUu6tdGp0+m46DFJrwWuBd4L/HNLrU+3tNFRGX+j\nZra+PN/1sybo0p1G51MFRd1YcMC6IbQImJ/uTLo7a4Oubc20z1jQQfuga22j/UuvTid9TNJr07WO\n7+1Qq1HQmVlGg/sesk4TVU0mv8xsMsqQT2mi+xjgMuAO4IKIWCLpJEkHpc2+AsyQtBQ4lnQTvYhY\nAlwA3A58Gzg6IkY61Uy1jgOOTbVmpNod2wCQdA9wKvA+SctrNxD8C+DLVPeR+CnwrbT808CbJd0F\nvDk978qfJJqVKNMt5NP1f2OhNAU4ayzogMURsYgqhM5NIbSSatBH2m4s6FaTgg6gXc3U5HHAQkmf\nBG5i3aBbr41U6x5gW2C6pEOAt0TE7VRBdzawJVXI1YPuAklHAvcCh274O2VmjQz26y0WAcdIWkh1\n3fMjaab8MuBTtZvVvAU4YWC9MrMyZcyniLgUuLRl2cdrPz9Fh+OPiDgZOLlJzbR8GdVNuFqXd2tj\ntw7LFwMvb7P8V8D+7fbpxINEs0JpNE8dB52Z5ZQrmySdR3UTmpmSllPdsXQaQER8gSpj3ko1I/4E\n8Cdp3UpJf0t1VgPASWM3sTGzyS1XPpkHiWbl8hU2ZlaifDP1h/VYH8DRHdadBZyVpydmtsnwsVM2\nHiSalcpBZ2YlcjaZWamcT9l4kGhWIEW2O3SZmWXjbDKzUjmf8vIg0axUee4OaGaWl7PJzErlfMrG\ng0SzQvniazMrkbPJzErlfMrHg0SzUvmUCTMrkbPJzErlfMrGg0SzEvm8ejMrkbPJzErlfMrKg0Sz\nUjnozKxEziYzK5XzKRsPEs1K5aAzsxI5m8ysVM6nbDxINCuUT5kwsxI5m8ysVM6nfDxINCuVg87M\nSuRsMrNSOZ+y8SDRrES++NrMSuRsMrNSOZ+y8iDRrFQOOjMrkbPJzErlfMrGg0SzUjnozKxEziYz\nK5XzKRsPEs0KJHzKhJmVx9lkZqVyPuXlQaJZqRx0ZlYiZ5OZlcr5lI0HiWYlCtDosDuxCZJgs82y\nlXvxGfdlq8WuO+erBYw+8POs9eJlL8pX7Ed35qsFxLPP5K336GNZ6+1y4bKs9YbK2bRxUL6cs03Y\npvZ74nzKahP77TDbhETDh5nZIDmbzKxUmfJJ0jxJd0paKun4Nus3l3R+Wn+tpN1q605Iy++UdECv\nmpJ2TzXuSjWnj6cNSb8l6eba41FJH0rrPiHp/tq6t/Z6DzxINCuUotnDzGyQnE1mVqoc+SRpCnAm\ncCAwBzhM0pyWzY4EVkXEHsBpwClp3znAfOBlwDzgc5Km9Kh5CnBaRMwGVqXafbcREXdGxF4RsRfw\nauAJ4OJan08bWx8Rl/Z6Lz1INCuVZ8OyzYaZWUYZP0lskE+n1f6t/0TSw7V1I7V1izb4dZnZxi9P\nPu0DLI2IZRHxDLAQOLhlm4OBc9LPFwL7S1JavjAino6Iu4GlqV7bmmmf/VINUs1DxtlG3f7ATyPi\nZz1fbQceJJqVqGnIeTas0WyYmWWSKZugWT5FxIdrWfDPwEW11U/WcuCgDXxlZraxy5dPuwD1mw4s\nT8vabhMRq4FHgBld9u20fAbwcKrR2la/bdTNB85rWXaMpFsknSVph3YvvM6DRLNCabTZowfPhplZ\nVpmyCZrlU91hrH/QY2a2Rh/5NFPS4trjqHqZNqVbh5adtsm1fDxtVDtVZ3EdBHyjtv7zwIuAvYAH\ngX9sU2MdjQaJkraXdKGkH0u6Q9LvNNnPzMavj/PquwXdpj4b9lVJlzibzAYn4zWJTf79V21KLwR2\nB75XW7xFyrxrJB3Sbr9h8rGT2eD1kU8rImJu7bGgVmY5sGvt+SzggZam1mwjaSqwHbCyy76dlq8A\ntk81Wtvqt40xBwI3RsQvxhZExC8iYiQiRoEvsf6E/HqafpL4T8C3I+IlwJ7AHQ33M7Pxan7KRLeg\n29Rnw14JzHA2mQ1Q82zqNoEFzfJpzHzgwogYqS17QUTMBQ4HTpeU8XtasvCxk9mg5Tnd9HpgdrrP\nwnSq/Gm97nkRcET6+R3A9yIi0vL56V4MuwOzges61Uz7XJlqkGpeMs42xqx31oWknWpP3w7c1utN\n6Pk9iZK2BV4HvA8gnRKS90upzGxdDa/paaCf2bDlfcxUdZ0NS58WtpsN66cN6DAbBmuyaQfSZJez\nyWwA+sumFWkQ10mTfBozHzh6na5EPJD+XCbpv4C9gZ827t0E8rGT2RBkOnaKiNWSjgEuA6YAZ0XE\nEkknAYsjYhHwFeBcSUupjmfmp32XSLoAuB1YDRw9NrnVrmZq8jhgoaRPAjel2oyzja2ANwN/1vKy\n/l7SXukduqfN+vX0HCQCvwn8EvgXSXsCNwAfjIjHG+xrZuMg2k+xj8OamSvgfqqAObxlm7GZqqup\nzVSluwV+XdKpwM6snalSu5ppn7HZsIW0nw1r2saYtrNhEfEgVTaNAltLuglnk9mEy5hN0CyfkPRb\nVBNCV9eW7QA8ERFPS5oJ7Av8fb6ubTAfO5kNWM58SjfFu7Rl2cdrPz8FHNph35OBk5vUTMuX0eb0\nz3G28QTV5Tyty9/Trk43TU43nQq8Cvh8ROwNPA60u031UWOnlDzL0/32w8xaZThlIn2iNzZzdQdw\nwdhsmKSxuwF+BZiRZqqOJf37TjNcYzNV3ybNVHWqmWodBxybas1g3dmwxm3AOrNh9bsZQjUbdivV\nKai7pn0aZdMz8VT3N8zMest0d9OG+QTVZNHCdJrVmJcCiyX9iOpUrU9HxO0b+Mpy6nns5OMmswmQ\nKZ+s2SeJy4HlEXFten4hbQ7E0nVQCwC21Y5++802UMO7A/a0qc6GSXo+cE3tKzB6ZtN2U2Y6m8w2\nUK5sgt75lJ5/os1+PwReka8n2fU8dvJxk1l+OfNpsuv5SWJE/By4L53uAdXt6EuarTPbNHk2rCtn\nk9mQOJt6cj6ZDYnzKZsmnyQC/BXwtXQ3nmXAn0xcl8yMtbdotu6cTWaD5Gzqh/PJbJCcT1k1GiRG\nxM1AtzuUmVluDrqenE1mQ+BsasT5ZDYEzqdsmn6SaGYD5tkwMyuRs8nMSuV8yseDRLNSOejMrETO\nJjMrlfMpGw8SzUoUvkOXmRXI2WRmpXI+ZeVBolmpPBtmZiVyNplZqZxP2XiQaFYg4fPqzaw8ziYz\nK5XzKS8PEs1K5aDL7qldtuQnf/3KbPWmPKlstfjNx/PVAl6/+y+z1rvv9+7KVitG8/5ya2re/8q0\n+eZZ6935D8/NVmuP996Wrda4OZsmRkyS8+Qi8y/QZlPy1ptMcv7OqedXrw+G8ykbDxLNCqXc/5Ga\nmWXgbDKzUjmf8vEg0axEvvjazErkbDKzUjmfsvIg0axUngwzsxI5m8ysVM6nbDxINCuUL742sxI5\nm8ysVM6nfDxINCuVg87MSuRsMrNSOZ+y8SDRrETh2TAzK5CzycxK5XzKyoNEs1I56MysRM4mMyuV\n8ykbDxLNCiRAmb9LzsxsQzmbzKxUzqe8CvnmSzNrpWj2MDMbJGeTmZUqVz5JmifpTklLJR3fZv3m\nks5P66+VtFtt3Qlp+Z2SDuhVU9LuqcZdqeb0DWjjHkm3SrpZ0uLa8h0lXZ7auFzSDr3eAw8SzUoU\nfTx6cNCZWTYZs8nMLKtM+SRpCnAmcCAwBzhM0pyWzY4EVkXEHsBpwClp3znAfOBlwDzgc5Km9Kh5\nCnBaRMwGVqXafbdR69sbI2KviJhbW3Y88N3UxnfT8648SDQrlEabPbrWcNCZWWY5smlNrd6TWO+T\n9Ms0WXSzpPfX1h2RJovuknREvldoZhurTPm0D7A0IpZFxDPAQuDglm0OBs5JP18I7C9JafnCiHg6\nIu4GlqZ6bWumffZLNUg1DxlnG93Ua9Xb6MiDRLNS5Zmtd9CZWV75znJoMokFcH6aLNorIr6c9t0R\nOBF4DVVmnOizCsysj3yaKWlx7XFUrcouwH2158vTMtptExGrgUeAGV327bR8BvBwqtHaVr9tjL0D\n35F0Q8trel5EPJhqPQg8lx584xqzEkW2i6/bBclrOm0TEasl1UPompZ9x0KoXc3GQdewjbGgC+CL\nEbEgLV8n6CT1DDozyyRfNkFtwglA0tgk1u0N9j0AuDwiVqZ9L6c6G+G8XJ0zs41Mf/m0ouUspTq1\nr95om07L230w12378bQBsG9EPJCOjS6X9OOI+H6b7XvyJ4lmherj4utus2ETEXTjCa3xBt2rqD5l\nOFrS69psa2YDlimboNlsPcAfSbpF0oWSdu1zXzObRDLduGY5sGvt+SzggU7bSJoKbAes7LJvp+Ur\ngO1Tjda2+m2DiBj78yHgYtaenfULSTulWjsBD/V6EzxINCtV81MmVkTE3NpjQa2Kg87M8sqTTdBs\nEuubwG4R8UrgCtaeat5kXzObbPKcDn89MDvdjG861b0TFrVsswgYuxb6HcD3IiLS8vnphn27A7OB\n6zrVTPtcmWqQal4ynjYkbS1pGwBJWwNvAW5rU6veRkceJJoVSGSbDXPQmVk2GbMJGkxiRcSvIuLp\n9PRLwKub7mtmk0uufEqXzRwDXAbcAVwQEUsknSTpoLTZV4AZkpYCx5JuohcRS4ALqE6b/zZwdESM\ndKqZah0HHJtqzUi1+24DeB5wlaQfUR2v/WdEfDvV+jTwZkl3AW9Oz7vyNYlmJYqoHhtcJlZLGgul\nKcBZY0EHLI6IRVQhdG4KoZVUgz7SdmMhtJq1IUS7mqnJ44CFkj4J3MS6Qde4DUnPAy6u7m3DVODr\nLUF3gaQjgXuBQzf4jTKzZjJlU7Jmwgm4nyoXDq9vIGmnsWuQgYOoDq6gyp9P1W5W8xbghFwdM7ON\nUMZ8iohLgUtbln289vNTdDj+iIiTgZOb1EzLl9Hmpn39tpHq7Nlh+18B+7db14kHiWaFyvVl1A46\nM8spYzY1mcT6QJq5X001wfS+tO9KSX9LNdAEOGnsJjZmNnnlyifzINGsWE2/Z8zMbJByZlODSawT\n6PAJYUScBZyVrzdmtrHzsVM+HiSalSiAfLeZt2Tz+x5nj2OvzVbv93/0ZLZa59316t4b9WH5G/P+\n/jx2yN7Zam1/1c+y1QJ48hWzstabdvkNWeuN/DrfTTd/9tFeXyPap/+7sL/tnU0TRxlvEzE6kq9W\nbptNyVsvMo4Kcv4dTDY5/x7G3QecTxl5kGhWKuecmZXI2WRmpXI+ZeNBolmhfF69mZXI2WRmpXI+\n5eNBolmp8t1B0MwsH2eTmZXK+ZSNB4lmhfJsmJmVyNlkZqVyPuXjQaJZgRQgX3xtZoVxNplZqZxP\neTW6jZOkD0taIuk2SedJ2mKiO2Y26Y02fExiziazIXA2NeJ8MhsC51M2PQeJknYBPgDMjYiXU33h\n7fyJ7pjZZKeIRo/JytlkNhzOpt6cT2bD4XzKp+npplOBLSU9C2wFPDBxXTIzAt/GuRlnk9kgOZv6\n4XwyGyTnU1Y9P0mMiPuBzwD3Ag8Cj0TEd1q3k3SUpMWSFj/L0/l7ajapRHWHriaPScrZZDYMzqYm\nmuSTs8ksN+dTTk1ON90BOBjYHdgZ2FrSu1u3i4gFETE3IuZOY/P8PTWbZDQajR6TlbPJbDicTb01\nySdnk1l+zqd8mty45k3A3RHxy4h4FrgI+N2J7ZbZJBeg0WaPSczZZDZozqamnE9mg+Z8yqrJNYn3\nAq+VtBXwJLA/sHhCe2VmPh2iN2eT2TA4m5pwPpkNg/MpmybXJF4LXAjcCNya9lkwwf0ys2j4mKSc\nTWZD4mzqyflkNiTOp2wa3d00Ik4ETpzgvphZjW/R3JuzyWzwnE3NOJ/MBs/5lE/Tr8Aws0Fz0JlZ\niZxNZlYq51M2TW5cY2YDpgg00uxhZjYoziYzK1XOfJI0T9KdkpZKOr7N+s0lnZ/WXytpt9q6E9Ly\nOyUd0KumpN1TjbtSzenjaUPSrpKulHSHpCWSPljb/hOS7pd0c3q8tdd74EGiWan8XT9mViJnk5mV\nKkM+SZoCnAkcCMwBDpM0p2WzI4FVEbEHcBpwStp3DjAfeBkwD/icpCk9ap4CnBYRs4FVqXbfbQCr\ngb+OiJcCrwWObun3aRGxV3pc2uut9CDRrFSZDsQ8G2ZmWWUcJDbIp2Ml3S7pFknflfTC2rqRWg4s\nyvgKzWxjlSef9gGWRsSyiHgGWEj1vad1BwPnpJ8vBPaXpLR8YUQ8HRF3A0tTvbY10z77pRqkmoeM\np42IeDAibqzehngMuAPYpfF718LXJFpjr7xRWevd8irPNHcUQIbv8anNXL0ZWA5cL2lRRNxe22zN\nTJWk+VQzVe9smanaGbhC0ovTPp1qjs2GLZT0hVT78+NoY2w27EZJ2wA3SLq81u/TIuIz/b4fIztu\nzaMHvqbf3Tq64Ox882yzTr82Wy0A7bFb1nrPbJPvtcYzz2SrBTDtO5m/WWCzKVnLbbtkWrZaz//n\nvL8nP+l3h0zZBI3z6SZgbkQ8IekvgL8H3pnWPRkRe+XpTQEi45e3KeP/1/LnCeP1ysV5v5DvlrkF\n/12U8HvSXz7NlFT/z2NBRIzdgXgX4L7auuVA68HDmm0iYrWkR4AZafk1LfuODdTa1ZwBPBwRq9ts\nP542AEiT8XsD9f80jpH0Xqqv4/nriFhFFwX8jZpZO4po9OjBs2FmllWmbIIG+RQRV0bEE+npNcCs\nrC/GzDYpfeTTioiYW3vUv6Km3SxLa6h12ibX8vG0Ue0kPQf4N+BDEfFoWvx54EXAXsCDwD+2qbEO\nDxLNStX8lImZkhbXHkfVqrSbDWsdbK0zUwXUZ6ra7dtpeePZsAZtrNFlNuwWSWdJ2gEzG5x8p5s2\nyae6I4Fv1Z5vkTLvGkmHdNrJzCaRPPm0HNi19nwW8ECnbSRNBbYDVnbZt9PyFcD2qUZrW/22gaRp\nVAPEr0XERWvflvhFRIxExCjwJapJuq48SDQrUQSMjjZ7eDas52yYmWXSXzZ1m8CCZvlUbSi9G5gL\n/ENt8QsiYi5wOHC6pBdleIVmtrHqL5+6uR6Yne6zMJ3qspjW654XAUekn98BfC8iIi2fn+7FsDsw\nG7iuU820z5WpBqnmJeNpI52h9RXgjog4td5ZSTvVnr4duK3Xm+BrEs1KledShn5mw5Y3nanqsHzN\nbFj6tLDdbFjjNrrNho39LOlLwH/0fBfMLJ/m2bQiDeI6aZJPSHoT8L+B10fE02PLI+KB9OcySf9F\ndcbBTxv3zsw2PRmOndL1f8cAlwFTgLMiYomkk4DFEbGIajB2rqSlVMcz89O+SyRdANxOdX+FoyNi\nBKBdzdTkccBCSZ+kug77K2l5X21I+j3gPcCtkm5ONT6a7mT695L2opqIuwf4s17vgweJZoVqeE1P\nL2tmroD7qQLm8JZtxmaqrqY2U5XuFvh1SadS3VRmbDZM7WqmfcZmwxbSfjasURu9ZsMi4sH0tNFs\nmJnlkymboEE+Sdob+CIwLyIeqi3fAXgiIp6WNBPYl+qmNmY2ieXKpzSwurRl2cdrPz8FHNph35OB\nk5vUTMuX0eb0z37biIiraH+GBhHxnnbLu/Eg0axUGYLOs2Fmll2+g7Am+fQPwHOAb1RzR9wbEQcB\nLwW+KGmU6tKZT7fcFdXMJiN/R2s2HiSalSiAUc+Gddi+79kwM8skYzZBo3x6U4f9fgi8IltHzGzj\nlzmfJjsPEs2KFE0urDYzGzBnk5mVyvmUkweJZqXyKRNmViJnk5mVyvmUjQeJZiXyKRNmViJnk5mV\nyvmUlQeJZkUKCJ8yYWalcTaZWamcTzl5kGhWKp8yYWYlcjaZWamcT9l4kGhWIp8yYWYlcjaZWamc\nT1l5kGhWKt+hy8xK5Gwys1I5n7LxINGsSOFTJsysQM4mMyuV8yknDxLNShR4NszMyuNsMrNSOZ+y\n8iDRrFSeDTOzEjmbzKxUzqdsPEg0K5WDzsxK5Gwys1I5n7LxINGsSOE7dJlZgZxNZlYq51NOEzJI\nfIxVK66IC3/WY7OZwIqJaD+Tkvs3lL5dsXfjTf3ere+FfW0dECMjE9SVyeuJlctXXPe1j/TKJhjC\n78mS5ps269ud4+/LBtbb+P/95/6nd/r5TbZq1LfbNrgz63E2FeAxVq24YvQbRWZTn0ru38D7tokc\nN8Fw+tdfNoHzKbMJGSRGxG/02kbS4oiYOxHt51By/0ruG5Tdv5L7th6fMpFdk2yCsn9PSu4blN0/\n9y0TZ1N2m0I2Qdn9c9/Gr/T+rcP5lI1PNzUrUYTv0GVm5XE2mVmpnE9ZeZBoVirPhplZiZxNZlYq\n51M2mw2x7QVDbLuJkvtXct+g7P6V3Ld1xOhoo4dNiJJ/T0ruG5TdP/ctA2fTUJX+e1Jy/9y38Su9\nf2vkyidJ8yTdKWmppOPbrN9c0vlp/bWSdqutOyEtv1PSAb1qSto91bgr1Zw+qDa6vgfhEbdZcbbb\nbEa8dvO3Ntr2O0999YaN5loBM9uoOZvMrFS58knSFOAnwJuB5cD1wGERcXttm78EXhkRfy5pPvD2\niHinpDnAecA+wM7AFcCL025ta0q6ALgoIhZK+gLwo4j4/CDa6PYeDfOTRDPrJkabPczMBsnZZGal\nypNP+wBLI2JZRDwDLAQObtnmYOCc9POFwP6SlJYvjIinI+JuYGmq17Zm2me/VINU85ABttHRUAaJ\nvT7CHRZJu0q6UtIdkpZI+uCw+9SOpCmSbpL0H8PuS52k7SVdKOnH6T38nWH3aYykD6e/09sknSdp\ni2H3qZsAYjQaPSyfUrMJNo58KjWbwPmUi7NpeErNJ2fThnE25ZMxn3YB7qs9X56Wtd0mIlYDjwAz\nuuzbafkM4OFUo7WtQbTR0cAHiekj3DOBA4E5wGHpY9MSrAb+OiJeCrwWOLqgvtV9ELhj2J1o45+A\nb0fES4A9KaSPknYBPgDMjYiXA1OA+cPtVQ8Rnq0fsMKzCTaOfCo1m8D5lIezaSgKzydn04ZxNuXS\nXz7NlLS49jiqVkntqrc877RNruWDaqOjYXyS2OQj3KGIiAcj4sb082NU/1B7jrQHSdIs4A+ALw+7\nL3WStgVeB3wFICKeiYiHh9urdUwFtpQ0FdgKeGDI/enJs/UDV2w2Qfn5VGo2gfMpN2fTUBSbT86m\n8XM25ddHPq2IiLm1R/3mPMuBXWvPZ7H+a1+zTXp/tgNWdtm30/IVwPapRmtbg2ijo2F8BUa7j0Jf\nM4R+dKXqDkJ7A9cOtyfrOR34G2CbYXekxW8CvwT+RdKewA3AByPi8eF2CyLifkmfAe4FngS+ExHf\nGXK3unqMVZddMXrBzIabr5jQzkweG0U2QbH5VGo2gfMpG2fT0GwU+eRs6puzKaOM+XQ9MFvS7sD9\nVJ+gHt6yzSLgCOBq4B3A9yIiJC0Cvi7pVKqbyswGrqP6NG+9mmmfK1ONhanmJQNso6NhDBLH9ZHn\nIEl6DvBvwIci4tFh92eMpLcBD0XEDZLeMOz+tJgKvAr4q4i4VtI/AccDHxtut0DSDlQzrrsDDwPf\nkPTuiPjqcHvWWUTMG3YfJqHiswnKzKfCswmcT9k4m4am+HxyNo2LsymjXPkUEaslHQNcRnWa7VkR\nsUTSScDiiFhE9envuZKWUn26Nz/tu0TVnURvpzoV++iIGAFoVzM1eRywUNIngZtSbQbURkfDGCQ2\n+Qh3aCRNowq5r0XERcPuT4t9gYMkvRXYAthW0lcj4t1D7hdUf6/LI2Js9vBCqqArwZuAuyPilwCS\nLgJ+Fyg26Gwois4mKMpYesMAAAFJSURBVDqfSs4mcD7Zxq/ofHI2jZuzqVARcSlwacuyj9d+fgo4\ntMO+JwMnN6mZli+jOqW8dfmEt9HNMK5JXPMRrqovcpxP9XHq0EkS1cj6jog4ddj9aRURJ0TErIjY\njep9+14pQRcRPwfuk/RbadH+VDMcJbgXeK2krdLf8f4UcmG4FaXYbIKy86nkbALnk20Sis0nZ9P4\nOZusZAP/JLHTR7iD7kcH+wLvAW6VdHNa9tE0Krfe/gr4WvoPbBnwJ0PuDwDpFI4LgRupPpa/CVjQ\nfS+bbArPJnA+bSjnk220Cs8nZ9OGcTZZkRRR1CntZmZmZmZmNkTDON3UzMzMzMzMCuVBopmZmZmZ\nma3hQaKZmZmZmZmt4UGimZmZmZmZreFBopmZmZmZma3hQaKZmZmZmZmt4UGimZmZmZmZreFBopmZ\nmZmZma3x/wH2IIuuxm03GAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f81ef8f04d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "mse_random = BAS_fit_random.mean_squared_error(E_array)\n",
    "mse_true = BAS_fit_true.mean_squared_error(E_array)\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=[15, 5])\n",
    "cf = axs[0].imshow(mse_random)\n",
    "axs[0].set_title('MSE Random')\n",
    "fig.colorbar(cf, ax=axs[0], shrink=0.8)\n",
    "cf = axs[1].imshow(mse_true)\n",
    "axs[1].set_title('MSE True')\n",
    "fig.colorbar(cf, ax=axs[1], shrink=0.8)\n",
    "cf = axs[2].imshow(mse_random - mse_true)\n",
    "axs[2].set_title('MSE Random - MSE True')\n",
    "\n",
    "fig.colorbar(cf, ax=axs[2], shrink=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that some voxels find the wrong local minimum with the random orientation guess - resulting in high fitting errors, and the true guess always finds the global minimum (to the solver tolerance)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter Cascading - Preparing a complex model with a simple model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea of parameter cascading is simple: we find some good initial guesses for parameters using some simpler models, and then put these guesses into a complex model to help it better find a global minimum.\n",
    "\n",
    "As an example, we use dipy's DTI module to get the initial guess for orientation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 10, 3)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from dipy.reconst import dti\n",
    "from dmipy.core.acquisition_scheme import gtab_dmipy2dipy\n",
    "\n",
    "# convert dmipy acquisition scheme to dipy gradient table\n",
    "gtab = gtab_dmipy2dipy(acq_scheme)\n",
    "# initialize the dti model\n",
    "tenmod = dti.TensorModel(gtab)\n",
    "# fit the dti model\n",
    "tenfit = tenmod.fit(E_array)\n",
    "# get the estimated eigenvectors of the DTI tensor\n",
    "evecs = tenfit.evecs\n",
    "# the principal eigenvector is our initial guess for the Stick model orientation\n",
    "evecs_principal = evecs[..., 0]\n",
    "evecs_principal.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then estimate the spherical angles from the cartesian ones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from dmipy.utils.utils import cart2mu\n",
    "mu_dti = cart2mu(evecs_principal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and do the same as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using parallel processing with 8 workers.\n",
      "Cannot estimate signal grid with voxel-dependent x0_vector.\n",
      "Setup brute2fine optimizer in 8.32080841064e-05 seconds\n",
      "Fitting of 100 voxels complete in 1.0415418148 seconds.\n",
      "Average of 0.010415418148 seconds per voxel.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f81efb33590>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAScAAAEICAYAAAAdoDKiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGgNJREFUeJzt3X+UHWWd5/H3h85vIIQxOmISDY4s\noqwKRuTHqAjOEVHhrOKKO/7AcSZHVxRH5gdwdnAcZ3bWOS5HPbiyWUQUGVADOqgBdFcZZVeRgOFH\niDgBlASjJEASQoAk3Z/9oyqeS5O+tzpdlVud+3mdUyf33qr7vd/udH/7eZ56qh7ZJiKibfbpdwIR\nEbuS4hQRrZTiFBGtlOIUEa2U4hQRrZTiFBGtlOIUrSPpeElrO56vlHR8H1OKPkhx6kHSLyVtkzR3\n1OsrJFnSwvL5fElXSdogaZOkOySdUe5bWB67ZdT29jE+8wZJT0h6VNJmSbdIOkfS9HL/RR0xtkna\n3vH82o7Pm9Ll63q9pB+Wn7Fe0r9KOqWmb1utbL/Y9g39ziP2rBSnau4D3rHziaR/D8wcdcxlwBrg\necAzgHcDvx11zBzb+3VsX+3ymWfa3h84CDgbOB1YJkm2378zBvBfga92xHxDry9G0mnA14EvA/OB\n3wfOB97c670Re0qKUzWXURSbnd5D8Yvd6RXApbYfs73D9s9sXzvRDy7j3QCcAhwDvHEi8SQJuAD4\nhO2LbW+yPWL7X23/WXnMH0j6vqSHypbg5ZLmdMT4a0kPlK2uuyWdWL4+JOk8SfeU+26RtKDc9xlJ\nazpagq/qiDdT0qWSHpF0F8X3sjPnX0p6Xfn4byV9TdKXy89YKWlRx7FHSvpZue/rkr4q6e8n8j2L\n/khxquYnwGxJh0kaAt4OfGUXx3xO0umSnlt3ArbvB5YDr+p1bA+HAguApV2OEfCPwHOAw8rj/xZA\n0qHAmcArypbd64Fflu/7KEUL82RgNvAnwNZy383Ay4DfA/4Z+LqkGeW+jwF/UG6vpyj+3ZwCXAnM\nAa4BLixzmwZ8A7i0/JwrgP/QI1a0VIpTdTtbT38E/Bx4YNT+twE/Av4GuK8ck3rFqGM2SNrYsR02\nzhx+TfFLNxHPKP9dN9YBtlfb/p7tJ22vp2hpvabcPQxMB14kaartX9q+p9z3p8B/sX23C7fZfqiM\n+RXbD5Wtyv9exji0fN9/BP7B9sO21wCf7fE13Gh7me1hiv+Xl5avHw1MAT5re7vtq4GfVvu2TE6S\nLpH0oKQ7a4r3T2VrdJWkz5Yt7b5IcaruMuA/AWfw9C4dth+xfY7tF1OM4awAvjnqP3eu7Tkd26px\n5jAPeHj30v+dh8p/DxrrAEnPknRl2XXbTNFKnAtF4QI+QtGSerA87jnlWxcA94wR8+zyB36TpI3A\nATtjUrTQ1nQc/qseX8NvOh5vBWaUg//PAR7wU69mX8Pe7VLgpDoCSToWOA54CXA4Rff6NV3f1KAU\np4ps/4piYPxk4Ooex24APkXxyzLRlg4A5djNyylaZxNxN8Uv7Fu7HPOPgIGX2J4NvJOiqweA7X+2\n/YcUg/8GPlnuWkPRNRud+6uAv6ZoIR1oew6wqSPmOorCttPudovXAfNG/UFYMNbBewPbP2TUH6xy\nzPC6cmzvR5JeWDUcMAOYRtGyncrTT+rsMSlO4/M+4ATbj43eIemTkg6XNEXS/sAHgNU7uzW7S9Is\nSa8B/oWii7JsIvHKVsVHgb+R9F5JsyXtI+kPJS0pD9sf2AJslDQP+MuOfA6VdEI5reEJ4HGKrh7A\nxcAnJB2iwkskPaOMtwNYD0yRdD7FmNROXwPOlXSgpPnAh3bzy/txmcuZ5f/DqcBRuxlrMlsCfMj2\ny4G/AP5HlTfZ/jHwA4oivw64fjda97VJcRoH2/fYXj7G7lkUg7EbgXspWhWj5w1t1FPnOX20y8dd\nKOlRir9cnwauAk6yPTKxrwJsL6UY1P8TinGs3wJ/T1EAAT4OHEnRuvkOT20pTgf+G7CBonv1LOC8\nct8FFIXmu8Bm4AsUUy6uB64FfkHRZXuCp3a3Pl6+fl/53st28+vaBryF4o/IRooW37eBJ3cn3mQk\naT/gWIoTDiuA/0nZhZf0Fkl37mK7vtz/AooTIPMphhBOkPTq/nwloNxsLvZmkm4CLrL9xX7n0hQV\nE4G/bftwSbOBu22POabYJc5fAjNsf6J8fj7whO1/qjPfqtJyir2KpNdIenbZrXsPxeDudf3Oa0+x\nvZnibPHboJjXJumlPd620/3Aa8rv3VSKwfB06yJqcihwG0WX9GzgNNtjTpuY7CRdQTHWdqiktZLe\nB/wx8D5JtwErgVMrhltKcbb1Dorv4W22v9VA2pWkWxcRrZSWU0S00phXrU/EtKGZnjnlgNrjjsxo\nJF20Y8InwJ4e84lttccE8Ej9uQI0NhF4SjP/ZzTQ4n/BizbVHhNg9cr9a4/5+MgWto08MaH/tNe/\ndl8/9PBw7wOBW25/8nrbtUz2rKqRn5yZUw7g2Hl/XHvcx//ds2qPCTDt4Sdqj6lV99UeE2Bk69be\nB+2GfaZPbybu3Gf0Pmg3ePv22mN+67pmxs3ffNjxtcf88ZZ/6X1QDw89PMxPr68233XooH+b2/uo\nejX0Zy0i2s7ACM20xOuQ4hQxoIzZ7mrdun5IcYoYYGk5RUTrGDPc4qlEKU4RA2yEFKeIaBkDwy0u\nTpUmYUo6ScW9oldLOqfppCJizxjBlbZ+6NlyKu+Z/TmK29OuBW6WdI3tu5pOLiKaY2B7i8ecqrSc\njqK4adq95f1yrqT6hYQR0VLGDFfc+qHKmNM8nnpjsLXAK0cfJGkxsBhgxlD90/UjomaG4fY2nCq1\nnHZ1/c7TviTbS2wvsr1o2tCsiWcWEY0qZohX2/qhSstpLU+9Sfx8ilu7RsSkJoZ32fZohyrF6Wbg\nEEkHU6zVdjrFEkkRMYkVA+KTuDjZ3iHpTIqb1A8Bl9he2XhmEdGoYp7TJC5OALaXMcEliSKifUYm\nc8spIvZOe0XLKSL2PkYMt/hO3SlOEQMs3bqIaB0jtnmo32mMKcUpYkAVkzAHsVvXwGoe029s6Frj\nBlY00cyZtccEGHrevEbijsyc2kjcHbc0NOtkn/r/4r/5xa+tPSbA9m/Mrj2m31/P158B8YhoHVsM\nu70tp/ZmFhGNG0GVtm4kzZD0U0m3SVop6eO7OGa6pK+W94S7SdLCXrml5RQxoIoB8VpKwJPACba3\nSJoK3CjpWts/6TjmfcAjtl8g6XTgk8DbuwVNyyliQO0cEK+ydY1T2FI+nVpuo+9ccirwpfLxUuBE\n9VhmOsUpYoANW5U2YK6k5R3b4s44koYkrQAeBL5n+6ZRH/W7+8LZ3gFsArouB51uXcSAGucM8Q22\nF40Zyx4GXiZpDvANSYfbvrPjkEr3heuUllPEABvxPpW2qmxvBG4AThq163f3hZM0BTgAeLhbrBSn\niAFVXPi7T6WtG0nPLFtMSJoJvA74+ajDrgHeUz4+Dfi+3X11hXTrIgaUEdvruXzlIOBL5UpN+wBf\ns/1tSX8HLLd9DfAF4DJJqylaTKf3CpriFDGgbGqZhGn7duCIXbx+fsfjJ4C3jSduilPEwOo9wbKf\nUpwiBpSpp+XUlBSniAGWm81FROsY5WZzEdE+xdJQ7S0B7c0sIho2+RfVjIi9kGFcs7/3tBSniAGW\nllNEtI6ttJwion2KAfGsvhIRrdPue4g3Upw8fQpPLux6H6ndMm3r47XHBBhe/1D9Qbdtrz8mwMaN\nzcRVMz+kQ4cd0khcHt5Uf8zh4fpjAlNP29L7oHHSpomvGFQMiGfMKSJaKDPEI6J1MkM8IlprMFf8\njYhWs2H7SIpTRLRM0a1LcYqIFsoM8YhonbZPJejZppO0QNIPJK0q10E/a08kFhFNU+1LQ9WpSstp\nB3C27Vsl7Q/cIul7tu9qOLeIaNikvoe47XXAuvLxo5JWUSwtnOIUMYkVZ+v2kmvrJC2kWAJm9Dro\nlGunLwaYPn1ODalFRJPaPgmzcmdS0n7AVcBHbG8evd/2EtuLbC+aNm3fOnOMiIaMlMtD9dq6qTIu\nLel4SZskrSi383cVq1OllpOkqRSF6XLbV1d5T0S0W41n66qOS//I9puqBu1ZnCSJYinhVbYvGFfK\nEdFqdZyJa2pcukpmxwHvAk7oaJKdPJEPjYj+s8UO71NpA+ZKWt6xLd5VzG7j0sAxkm6TdK2kF/fK\nr8rZuhuhxecbI2K3jaNbt8H2om4H9BiXvhV4nu0tZePmm0DXm32198KaiGjUzjGnKlsvvcalbW+2\nvaV8vAyYKmlut5i5fCVigNUxIF5lXFrSs4Hf2rakoygaRl1vQZviFDGgapzntHNc+g5JK8rXzgOe\nC2D7IuA04AOSdgCPA6fbdregKU4RA6yOy1eqjEvbvhC4cDxxmylOj25l6Ae31h939uz6YzZEQ81c\nFuBm7sGPpjbzozCy+leNxG2EJ75owJ7ikYn/INiwIzebi4g2avPlKylOEQOq7dfWpThFDDCnOEVE\nG03q+zlFxN7JzphTRLSSGM7Zuohoo4w5RUTrtH31lRSniEHlYtyprVKcIgZYztZFROs4A+IR0Vbp\n1kVEK+VsXUS0jp3iFBEtlakEEdFKGXOKiNYxYiRn6yKijVrccEpxihhYGRCPiNZqcdMpxSligA1e\ny0mgKfWH1jMOrD0mwJQD6l/VZccD62qPCfCLz728kbgv/IuVjcTVzBmNxGVO/f9nr/vmbbXHBPjf\nxx9ce0w9MvHVfQyMjLS3OLV3qD4immXAqrZ1IWmBpB9IWiVppaSzdnGMJH1W0mpJt0s6sld66dZF\nDLCa5jntAM62fauk/YFbJH3P9l0dx7wBOKTcXgl8vvx3TGk5RQwyV9y6hbDX2b61fPwosAqYN+qw\nU4Evu/ATYI6kg7rFTcspYmBpPAPicyUt73i+xPaSp0WUFgJHADeN2jUPWNPxfG352piDsylOEYOs\nerdug+1F3Q6QtB9wFfAR25tH7x7vp6c4RQwqg2s6WydpKkVhutz21bs4ZC2woOP5fODX3WJmzCli\noKni1iWCJOALwCrbF4xx2DXAu8uzdkcDm2x3nW9TueUkaQhYDjxg+01V3xcRLVbP2brjgHcBd0ha\nUb52HvBcANsXAcuAk4HVwFbgvb2CjqdbdxbFKHz9s98ioj9qKE62b6RH88q2gQ+OJ26lbp2k+cAb\ngYvHEzwiWqymSZhNqdpy+jTwV8D+Yx0gaTGwGGAGsyaeWUQ0rs03m+vZcpL0JuBB27d0O872EtuL\nbC+aqum1JRgRDRpRta0PqrScjgNOkXQyMAOYLekrtt/ZbGoR0TRN5paT7XNtz7e9EDgd+H4KU8Re\noOqlK30qYJmEGTGw+jfYXcW4ipPtG4AbGskkIva8Fnfr0nKKGGQj/U5gbClOEYNq5zynlkpxihhg\nbT5bl+IUMchaXJxyV4KIaKWGWk6CoYmvDjGat2ytPSbA8Pr1tce85/Ijao8JsPTYCxuJe95/7no7\n5933+BONhB1q4LqL//P6F9YeE2D4ofpX4vHIcC1x0q2LiPYxfbs0pYoUp4hBlpZTRLRRunUR0U4p\nThHRSilOEdE2crp1EdFWOVsXEW2UllNEtFOKU0S0TsvHnHJtXcQgq+k2vZIukfSgpDvH2H+8pE2S\nVpTb+b1ipuUUMcBU383mLgUuBL7c5ZgfjWe18LScImLCbP8QeLjOmClOEYOserdurqTlHdvi3fi0\nYyTdJulaSS/udXC6dRGDanwD4htsL5rAp90KPM/2lnINzG8Ch3R7Q1pOEYNsD61bZ3uz7S3l42XA\nVElzu70nxSlikO2h4iTp2ZJUPj6KovY81O096dZFDChR39k6SVcAx1OMTa0FPgZMBbB9EXAa8AFJ\nO4DHgdPt7rczTXGKGFQ1TsK0/Y4e+y+kmGpQWYpTxCBr8QzxFKeIQTZwxWnWDEZeVv9KFv5/t9Ue\nE2Bo9uzaYx7y/tW1xwQ4880fbiTubP+kkbi4nlVCRhvZ8ljtMX9xwcLaYwK88Kz6c9XmelY3avO1\ndWk5RQyyFKeIaB3Xem1d7VKcIgZZWk4R0UYZc4qIdkpxiojWqenSlKakOEUMKNHubl2lC38lzZG0\nVNLPJa2SdEzTiUVE83auXddr64eqLafPANfZPk3SNGBWgzlFxJ7S4pZTz+IkaTbwauAMANvbgG3N\nphURe0SLi1OVbt3zgfXAFyX9TNLFkvYdfZCkxTtv4bl9R/3T9SOiZhW7dP3q1lUpTlOAI4HP2z4C\neAw4Z/RBtpfYXmR70dQpT6tdEdFGe+hmc7ujSnFaC6y1fVP5fClFsYqISU4j1bZ+6FmcbP8GWCPp\n0PKlE4G7Gs0qIvaINnfrqp6t+xBweXmm7l7gvc2lFBF7xN4wCdP2CmAiy8JERBtN9uIUEXufts8Q\nT3GKGGAaaW91SnGKGFQtH3PKopoRA6yus3WSLpH0oKQ7x9gvSZ+VtFrS7ZJ6TkdKcYoYZPVNwrwU\nOKnL/jcAh5TbYuDzvQI20q3T9h1MXdN1peHd4jkH1B4TgLm/V3tIbdxce0yA2Vfe3Ejcy9b830bi\nXvxIM/N1v3jny2qPeeif3V17TIDhrVtrj+mRela1qXFRzR9KWtjlkFOBL5er/P6kvNPJQbbXjfWG\ntJwiBln1ltPcndfOltvicX7SPGBNx/O15WtjyoB4xKAa3+orG2xPZK6jdp3B2FKcIgbUHp7ntBZY\n0PF8PvDrbm9Ity5ikNnVtom7Bnh3edbuaGBTt/EmSMspYqDV1XKSdAVwPMXY1FrgY8BUANsXAcuA\nk4HVwFYqXJ+b4hQxqGqchGn7HT32G/jgeGKmOEUMsCxHHhGtlOIUEe1j6hrsbkSKU8QAyy1TIqKd\nUpwiom1ys7mIaCc7N5uLiJZqb21KcYoYZOnWRUT7GEi3LiJaqb21KcUpYpClWxcRrZSzdRHRPi1f\nGqqR4nTgIVt5y1XLa4979WtfWntMALbvqD3k8KZHa48JMOX3n9lI3Lee9dFG4u77nRWNxH3+9jvq\nDzpjev0xgaEDZtceU5uHJh4DUK6ti4hWyl0JIqKN0nKKiPYZxDGniJgMcm1dRLRVunUR0TrjW1Rz\nj0txihhkLW45VVpUU9KfS1op6U5JV0ia0XRiEbEHuOLWBz2Lk6R5wIeBRbYPB4aA05tOLCKap5GR\nSlvPONJJku6WtFrSObvYf4ak9ZJWlNuf9opZtVs3BZgpaTswix5rnEfEJGBqmYQpaQj4HPBHwFrg\nZknX2L5r1KFftX1m1bg9W062HwA+BdwPrKNY4/y7u0hwsaTlkpY/+sj2qp8fEX0ijFxt6+EoYLXt\ne21vA64ETp1oflW6dQeWH3Qw8BxgX0nvHH2c7SW2F9letP+BUyeaV0TsCXa1DebubHyU2+KOKPOA\nNR3P15avjfZWSbdLWippQa/UqnTrXgfcZ3s9gKSrgWOBr1R4b0S0WfWzdRtsLxpjn3YVedTzbwFX\n2H5S0vuBLwEndPvAKmfr7geOljRLkoATgVUV3hcRbbZzzKnK1t1aoLMlNJ9R49K2H7L9ZPn0fwEv\n7xW0ypjTTcBS4FbgjvI9S3qmGxGtV9PZupuBQyQdLGkaxdn8a57yOdJBHU9PoUIDp9LZOtsfAz5W\n5diImCxcyyRM2zsknQlcTzHV6BLbKyX9HbDc9jXAhyWdAuwAHgbO6BU3M8QjBpWpbYa47WXAslGv\nnd/x+Fzg3PHETHGKGGS5ti4i2ig3m4uIdkpxiojWsWG4vf26RorTI/82i6tPekX9gWdWuonCuD3y\nyoN6HzROc773WO0xAR4+fmEjcQ/8zujLoOqx7L6bGol70nPHmg84ARUucN2tsFvq/1nwyHBNgdJy\niog2SnGKiNYxkHuIR0T7GDxgY04RMQmYwRsQj4hJImNOEdFKKU4R0T71XPjblBSniEFlGpvbVYcU\np4hBlpZTRLTPAF6+EhGTgMGZ5xQRrZQZ4hHRShlziojWsXO2LiJaKi2niGgf4+Ga7gvVgBSniEGV\nW6ZERGtlKkFEtI0Bp+UUEa3j3GwuIlqqzQPicgOnEiWtB35V4dC5wIbaE2jOZMp3MuUKkyvfNuT6\nPNvPnEgASddRfC1VbLB90kQ+b7waKU6VP1xabruBNX6aMZnynUy5wuTKdzLlOpk1sxBcRMQEpThF\nRCv1uzgt6fPnj9dkyncy5QqTK9/JlOuk1dcxp4iIsfS75RQRsUspThHRSn0rTpJOknS3pNWSzulX\nHr1IWiDpB5JWSVop6ax+51SFpCFJP5P07X7n0o2kOZKWSvp5+T0+pt85dSPpz8ufgzslXSFpRr9z\n2lv1pThJGgI+B7wBeBHwDkkv6kcuFewAzrZ9GHA08MEW59rpLGBVv5Oo4DPAdbZfCLyUFucsaR7w\nYWCR7cOBIeD0/ma19+pXy+koYLXte21vA64ETu1TLl3ZXmf71vLxoxS/PPP6m1V3kuYDbwQu7ncu\n3UiaDbwa+AKA7W22N/Y3q56mADMlTQFmAb/ucz57rX4Vp3nAmo7na2n5LzyApIXAEcBN/c2kp08D\nfwW096rOwvOB9cAXyy7oxZL27XdSY7H9APAp4H5gHbDJ9nf7m9Xeq1/FSbt4rdVzGiTtB1wFfMT2\n5n7nMxZJbwIetH1Lv3OpYApwJPB520cAjwFtHn88kKKFfzDwHGBfSe/sb1Z7r34Vp7XAgo7n82lx\n81jSVIrCdLntq/udTw/HAadI+iVFd/kESV/pb0pjWgustb2zJbqUoli11euA+2yvt70duBo4ts85\n7bX6VZxuBg6RdLCkaRSDitf0KZeuJIliTGSV7Qv6nU8vts+1Pd/2Qorv6/dtt/Kvu+3fAGskHVq+\ndCJwVx9T6uV+4GhJs8qfixNp8QD+ZNeX+znZ3iHpTOB6ijMel9he2Y9cKjgOeBdwh6QV5Wvn2V7W\nx5z2Jh8CLi//SN0LvLfP+YzJ9k2SlgK3UpzF/Rm5lKUxuXwlIlopM8QjopVSnCKilVKcIqKVUpwi\nopVSnCKilVKcIqKVUpwiopX+P7+gTfpeYhkAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f81efc2d650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ball_and_stick.set_initial_guess_parameter('C1Stick_1_mu', mu_dti)\n",
    "\n",
    "BAS_fit_dti = ball_and_stick.fit(acq_scheme, E_array, Ns=5)\n",
    "mse_dti_mu = BAS_fit_dti.mean_squared_error(E_array)\n",
    "\n",
    "plt.imshow(mse_dti_mu)\n",
    "plt.title(\"MSE DTI Cascading\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see we that using DTI cascading we have similar results as giving the ground truth orientations (as we did before)\n",
    "\n",
    "This process is the same for any other paramater, e.g. using CSD to initialize multiple orientations.\n",
    "\n",
    "NOTE: currently dmipy does not accept initial conditions for volume fractions, but everything else is fine."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}