howto optimal inventory-update.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A StoDynProg usage example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(\"/home/deheeger/Github/stodynprog/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from stodynprog import SysDescription\n",
    "invsys = SysDescription((1,1,1), name='Shop Inventory')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dyn_inv(x, u, w):\n",
    "    'dynamical equation of the inventory stock `x`. Returns x(k+1).'\n",
    "    return x + u - w\n",
    "# Attach the dynamical equation to the system description:\n",
    "invsys.dyn = dyn_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.stats as stats\n",
    "demand_values = [0,   1,   2,   3]\n",
    "demand_proba  = [0.2, 0.4, 0.3, 0.1]\n",
    "demand_law = stats.rv_discrete(values=(demand_values, demand_proba))\n",
    "demand_law = demand_law.freeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.2,  0.1])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demand_law.pmf([0, 3]) # Probality Mass Function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 2, 0, 0, 1, 1, 2, 1, 1])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demand_law.rvs(10) # Random Variables generation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "invsys.perturb_laws = (demand_law, ) # a list, to support several perturbations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def admissible_orders(x):\n",
    "    'interval of allowed orders U(x_k)'\n",
    "    U1 = (0, 10)\n",
    "    return [U1, ] # tuple, to support several controls\n",
    "# Attach it to the system description.\n",
    "invsys.control_box = admissible_orders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(h,p,c) = 0.5, 3, 1\n",
    "def op_cost(x,u,w):\n",
    "    'operational cost of the shop'\n",
    "    holding = x*h\n",
    "    shortage = -x*p\n",
    "    order = u*c\n",
    "    return np.where(x>0, holding, shortage) + order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "op_cost(1,1,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7. ,  4. ,  1. ,  1.5,  2. ])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Vectorized cost computation capability (required):\n",
    "import numpy as np\n",
    "op_cost(np.array([-2,-1,0,1,2]),1,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "invsys.cost = op_cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dynamical system \"Shop Inventory\" description\n",
      "* behavioral properties: stationnary, stochastic\n",
      "* functions:\n",
      "  - dynamics    : __main__.dyn_inv\n",
      "  - cost        : __main__.op_cost\n",
      "  - control box : __main__.admissible_orders\n",
      "* variables\n",
      "  - state        : x (dim 1)\n",
      "  - control      : u (dim 1)\n",
      "  - perturbation : w (dim 1)\n"
     ]
    }
   ],
   "source": [
    "invsys.print_summary()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from stodynprog import DPSolver\n",
    "dpsolv = DPSolver(invsys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([-3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.,  6.])]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xmin, xmax = (-3,6)\n",
    "N_x = xmax-xmin+1 # number of states\n",
    "dpsolv.discretize_state(xmin, xmax, N_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n",
      "0\n",
      "3\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([array([ 0.,  1.,  2.,  3.])], [array([ 0.2,  0.4,  0.3,  0.1])])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Nw = len(demand_values)\n",
    "print Nw\n",
    "print demand_values[0]\n",
    "print demand_values[-1]\n",
    "dpsolv.discretize_perturb( demand_values[0], demand_values[-1], Nw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dpsolv.control_steps=[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SDP solver for system \"Shop Inventory\"\n",
      "* state space discretized on a 10 points grid\n",
      "  - Δx = 1\n",
      "* perturbation discretized on a 4 points grid\n",
      "  - Δw = 1\n",
      "* control discretization steps:\n",
      "  - Δu = 1\n",
      "    yields 11 possible values\n"
     ]
    }
   ],
   "source": [
    "dpsolv.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "J_0 = np.zeros(N_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "value iteration...\r",
      "value iteration run in 0.00 s\n",
      "[ 9.   6.   3.   0.   0.5  1.   1.5  2.   2.5  3. ]\n",
      "[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]\n",
      "value iteration...\r",
      "value iteration run in 0.00 s\n",
      "[ 4.  3.  2.  1.  0.  0.  0.  0.  0.  0.]\n",
      "value iteration...\r",
      "value iteration run in 0.00 s\n",
      "value iteration...\r",
      "value iteration run in 0.00 s\n",
      "value iteration...\r",
      "value iteration run in 0.00 s\n",
      "[ 6.  5.  4.  3.  2.  1.  0.  0.  0.  0.]\n"
     ]
    }
   ],
   "source": [
    "J_0 = np.zeros(N_x)\n",
    "J,u = dpsolv.value_iteration(J_0)\n",
    "print J\n",
    "print u[...,0]\n",
    "J,u = dpsolv.value_iteration(J)\n",
    "print u[...,0]\n",
    "J,u = dpsolv.value_iteration(J)\n",
    "J,u = dpsolv.value_iteration(J)\n",
    "J,u = dpsolv.value_iteration(J)\n",
    "\n",
    "print u[...,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 6.5)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEhCAYAAACXwKDgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYFOWZ/vHvzRkUGZCggsAAUYgiB5WQGFEiSTTGQ2Lc\nmChGTDYmccVVY1bz0xViskp2s8ZV18TEAxE1yYpnRaNRENQoig4IKEZgRBkUFRhRDjLM8/ujqrFp\n5lA90z31zvTzua6+pqu6uvqeGnin+6l3npKZ4ZxzrrS0SzuAc865lueDv3POlSAf/J1zrgT54O+c\ncyXIB3/nnCtBPvg751wJ8sHfFZykOZJ+30KvVSvp1JZ4rfj15kj6Q4H21aLZm0vS+Dhz33i5PF4+\nLO1sLn8++JcYSXtJulbSSklbJa2VNFPSyCbs61JJK+t46OvABc1PGySLb4WwN3BXgfaVhlVE38P8\ntIO4/PngX0Ik9QdeAD4H/AgYAnwN+Bh4VtLRhXgdM9tgZh8WYl9pkNSxyPvvBGBma81sazFfq5jM\nrDb+HmrSzuLy54N/aflfoD3wRTP7q5m9ZWbPm9mpwBPAdEldACRNlfQPSadKWiFps6RHJQ2MH58E\nXA4MjD/610q6LH5sp9JIvHyjpF/GnzTWS7pckZ9Lejte/8vssPFrPydpg6R3JT0oab98v2lJx0pa\nIGmLpHck/a+kblmPT5f0mKTJkiqBLZI6Sxoo6RFJmyStkjS5jn13jI9V5hgtlnRWzja18b7vkLQB\n+GPW+tNytvuxpBmSPpD0pqSLc/a1p6Q7JX0oaY2kyzL5G/j+M+WZ0yQ9Hn8/yyWdkrPdUEkPSdoY\n3+6XNCTBfg/LWtdH0i3xz3SzpFclnRk/tkLSz3L2sVv8vZ6Wu39XXD74lwhJPYFjgevqeVd+JbAX\n8KWsdfsQfUI4GRgH7AHcHT/2Z+BXwFtEH/33Bn4dP1ZXaeRkol88hxGVhC4FHgY6A4cDFwL/T9Ix\nWc/pRPQLZnScazvwUD7vzCWNAO4H5gAjgDOA44Df5Wz6WWA8cHy8XQ1wD9ATODJefzxwcM7z/kBU\n5joLGBbn/ZWk7+VsNwV4Kv5eLs1an3ucpsRZRxL9TK6QdFTW47cABxF9YpsAlAMn1rGfuvwncGO8\n7zuA2yWNApDUFXiU6JgfEX/PuwOPJD3e8T6ejPOdSnQ8zgYy/95+D3w/52nfJvrkeWeS13AFZGZ+\nK4Eb0eBWC5xYz+O94sd/Ei9PjZcHZ22zX7zui/HypcDKOvY1G/h91vIc4MWcbRYDC3PWVQD/1cD3\nkMn4+ax1tcCpDTxnBvBszroTiH6R9I+XpwPrgG5Z23wp3vens9b1BjZlvjdgULyf/XP2fxnwUk7G\nP9SRbafs8fLVOdssBa6o6/jH6zoQ1d4fbeAYlMfP+3nO+qeBW+P73wc+AnplPd4n/n5Pj5fHx/vp\nm7Pfw7L2sTnzeB05+gBbgQlZ6/4O/Cbt/x+lePN3/q4h75rZisyCmf0DeA84MM/9GLAwZ93bwKI6\n1n0qsyBplKR74nLBB8Ab8UMD83jtA4C5OevmAoofy3jFzDblPO89M3t9xzdh9h6wLGubQ+P9LMgq\nlWwEfgZ8Ouc1k54UrchZriIaNDOZAJ7NylRDdB4nib/nLD/DJz/LA4ElZrYua99rib7fA0jmkHgf\nVXU9GO/vPuAHAJKGA2OJPj25FtYh7QCuxbxONAgfRPQfMFdmEFhWx2OFsC1n2epYB3EpMq7JP0o0\nUE8C3iEaaJcQlSbyoQTbbGp8k132lXnz9Pk6np9bhvko4f4/rmNd7pu03H0n+f7qk72vuvaT774b\n2/53wCxJewL/DDxjZkvzfA1XAP7Ov0TE7+hmAedI6l7HJj8jeuedfeLwU5IGZxYk7U9U+sj8Z/2Y\nqI5fsJhZ9z8Tv9YlZjbXzJYRlX3yHYyWENWwsx0Zv9aSBp63FOgtacc7eEm9gaFZ2yyIvw40sxU5\nt7qmwDZF9jHJHPfsE6wdiN5xJ/H5nOXDsva5BDggHpQz+94L2J+oRJfEC/E++jWwzWyiMtWPgIn4\nu/7U+OBfWv6F6ETmE5KOltRf0hhJdxDVcyfZzlMPNwG3SDpE0qFEs1ReMrMn4sdXAHtL+pyk3vEJ\nP4gG6OxBOnc5ybo3iOrD50oaImkC8D/kP8f+v4CDJV0laVh8Qvla4DYze6u+J5nZ34hKVbfFx2gU\ncDvRLzzF27wO3Az8QdJESZ+WNFLS9yT9W54567PjmMRltweA/5V0hKQDgBuITsQnOS7fk/QdSftL\nupxoyu9V8WO3A+8Cf5E0WtIhRCf13wL+kjDrn4h+bvdLmiBpUPz1W5kNLCr0/57ovEi7PPbtCswH\n/xJiZquI3iU+RzRovE70aaAj0UnUR3OesibebiYwj2jWxklZj99LNEvjIWAt8NPMS7HzYFTX7J8G\n18X19YnAl4neef4n8BOiE4yJmdnLRCd4jyCqp99KNID+qJEsEM3iqSYqPd0PPAi8mLPtWcBvgEuI\n3j3/DTgdWJ5Pzoa+hZzXO5PoeDxMND33LaLy2JYE+7o4zrsQOA04zcwqAMxsC/AVol+4c4lO0m8E\njrGd5/HX9TMj3sdmok9Vi4l+cSwl+kXbJec5t8Rfb49f16VA0S9i53YmaSrR4JD3vHrXciS1B14F\n7jWzn9azTTnRp7TDzeyZlktXN0kHAi8DI+Nfzi4FfsLXuVZE0jiiv8d4CegOnA8MIJquGjRFf9n8\nKaK/X3jCB/50ednH1aeQPWxc4bQnKjFVEJV9yonm/Td08hrC+FmeSnSydyDw45SzlDwv+zjnXAny\nd/7OOVeCgh7858yZkyk9BHPzTJ6pFHJ5plafqVGhD/5pR9iFZ0rGMyUXYi7PlExrzhT04F9ZWZl2\nhF14pmQ8U3Ih5vJMyYSYKamgB3/nnHPF0X7q1KlpZ6hXWVnZ1PLy8rRj7KSsrAzP1DjPlFyIuTxT\nMiFmAigvL/95Y9uEPtUz6HDOOReoRhsgBl32ac0nU1qSZ0omxEwQZi7PlEyImZIKevB3zjlXHF72\ncc65tqd1l32cc84VR9CDf4j1NM+UjGdKLsRcnimZEDMllergL6lM0kxJr0haKulzaeZxzrlSkWrN\nX9IfgSfN7Ob4WqS7mVl11iYGsHHLNhaurubwIb1Tyemcc61MuDV/ST2AcWZ2M4CZ1eQM/EA08F8/\nbwUj+/Vo8YzOOddWpVn2GQS8K+kWSS9K+oOkbtkbzJkzh+vnreDscYPp3qVjSjF3FmKNzzMlE2Im\nCDOXZ0omxExJpXkZxw7AwcA5Zva8pKuJLjB9WWaDmTNn8nrVe2x+ehgQ/Sn1qFGjGD9+PPDJgW/J\n5YqKilRfv67ljFDyhLpcUVERVB7/+eW3HOLPL4TxIHM/02Ru0qRJO7ZpSGo1f0l7A383s0Hx8uHA\nxWZ2XNZmNvbXT3DZMZ/ha8P3SSWnc861QuHW/M3sbeBNSfvHq74E7HId0lqDqQ+/wo3PrCDwP0hz\nzrlWI+15/pOB2yUtBEYAV2Q/OGfOHCYfMQSAG56u5IpHl1GzvbblU+ZkCo1nSibETBBmLs+UTIiZ\nkkqz5o+ZLQTGNLTNd8cOpF9ZV6bMWsq9i6pYU72ZaScexO6dU43unHOtWqvp7bO4qpqf3LOIdZu2\nMbj3bvzmpBH07dE1zWzOOReqRmv+rWbwB6iq3sx5dy1i5fsf0atbR646aSQH7rNHWtmccy5U4Z7w\nTSK3nta3R1duOvVgxgzoybpN2/jhn1/kiWVrU80UAs+UTIiZIMxcnimZEDMlFfTgX5fuXTpyzckj\nOXHEPmytqeWi+xczY/4bPhPIOefy0KrKPjs9YMaM+au4du5yAL4+oi8XfWl/OrRvdb/PnHOu0NpW\nzb8ujy9by5RZS9laU8vYgT19JpBzzrW1mn9dJgztw+9OGU2vbh157o31fP+OBVRVb041U0vzTMmE\nmAnCzOWZkgkxU1JBD/5JDe/bg1smHsqgPXdjxXsfceZtL7BkzQdpx3LOuWC1+rJPto1btnHRfYt5\nftV6Ondox+XHHsBRQ/sUK5tzzoWq+DV/Sb3MbJ2k9ma2vVk721Xe4Wq21zLtsWXc9/IaAM49cggT\nxwxAavRYOOdcW9EiNf8p8df+ko4twP52aEo9rUP7dlxy9DDOiXsCXfPk8oL2BAqxxueZkgkxE4SZ\nyzMlE2KmpJo8+Es6QlJH4A5JRwG/AD5TsGTNIIkzxg5k2gnD6dyhHfcuquK8uxby4daatKM551wQ\nmlz2kXQdsBVoD3wROB941sw2FS5e/mWfXN4TyDlXglpmnr+kPsA44EAzu7zZO/xEQc5Gr96wmfPv\nXsjK9zd5TyDnXCkobs1f0tGS9jaztWZ2V4EH/oLV0/qVdeWmUw8pSE+gEGt8nimZEDNBmLk8UzIh\nZkoq78Ff0t8k3STp20AF8NXCxyq8HT2BDvKeQM45l3fZJz7JOxaYQHTpxRVmdkYRskGByj477dCM\nW+ev4jrvCeSca7taZJ7/P5nZnc3aSf2K9rbcewI559qwwtf8Jf1J0lxJZ0k6ABjSpGgJFLOe1tSe\nQCHW+DxTMiFmgjBzeaZkQsyUVFNqHXcCpwJ9gP8E3itoohY0vG8Pbj7tUAbt2c17AjnnSkpTav79\ngWFm9lhxIu2kRc7Gek8g51wb0/b7+ReK9wRyzrUhbb+ff6Ek7QkUYo3PMyUTYiYIM5dnSibETEml\nOvhLqpS0SNJLkuanmSXO4z2BnHMlIXHZR1I7MytMa8xP9rkSOMTM1tWzSWo1qeyeQHvv0Zn//sYI\n9u/TfcfjG7dsY+Hqag4f0jutiM45V5/ClH0kdQA+ktS52ZHq2H0R9tls2TOB3v5gK9+7fQHz34h+\nR23cso3r561gZL8eKad0zrmmSTT4m1kN8A+g0G9zDfibpBck/SD3wbTradk9gbbW1DL5zgp+/oc7\nuX7eCs4eN5juXTqmmi8j7eNUF8+UXIi5PFMyIWZKKp8/ab0NeEDSNcCbZJVkzOyJJr7+F8xsjaRP\nAY9JetXM5mUenDlzJtOnT6e8vByAsrIyRo0axfjx44FPDnyxl685+QimPbaMGfc+zI1Vr3PxhaPZ\nvXOHFnv9xpYzQskT6nJFRUVQefznl99yiD+/ioqK1PNk7ldWVgIwadKkHds0JJ+af2V8d5cnmNmg\nRDtpeP9TgA/N7L+zd93c/RbKB5s/5ry7FvFy/Edgxx64N/9+9DDvCeScC1G48/wldQPam9lGSbsB\njwI/N7NHszYLYvDP1PjPHjeY+W+s57KHlvDxduOQ/j349TdGek8g51xoCjvPX9JXJN0s6cF4+dD4\nEo5NsRcwT1IF8BzwYM7Av8vH4rQsXF29o8bffs1Sbvj2wZR17cCCN6sT9wQqplCOUzbPlFyIuTxT\nMiFmSirx4C9pMvBbohO/R8SrtwC/bMoLm9lKMxsV34ab2ZVN2U9LOHxI751O7g7v24PpE8d4TyDn\nXKuVT81/BTDBzFZKWm9mPSW1B941s15FyhdE2ac+3hPIOReogpZ9diea5ZOtE9FF3EuSXx3MOdda\n5TP4zwMuzlk3GZhduDg7C7GelpspaU+glswUAs+UXIi5PFMyIWZKKp/BfzLwDUlvALtLeg04BfhJ\nUZK1It4TyDnX2uQ11VNSO2AMMJCoBDTfzLYXKRsEXvOvS3ZPoMG9d+M3J42gb4+uacdyzpWWcOf5\nJxR0uPqs3rCZ8+9eyMr3N9GrW0euOmkkB+6zR9qxnHOlo3knfCX9QtLl8dfM/V1uhcu7sxDraUky\nZfcEWrdpGz/884s8sWxtqplammdKLsRcnimZEDMl1VjNv3982xf4NNEJ3wnx/Qnx8n7FDNha+Uwg\n51zI8pnn/2fgTjO7K2vdScC3zOzbRcrX6kdKM+PW+au4bu5yAL4+oi8XfWl/7wnknCumwtX8JX0A\n9Mw+wSupI/C+mRWroN3qB/+Mx5etZcqspWytqWXswJ5MO/Eg7wnknCuWgv6R1+vAOTnrfhyvL4oQ\n62lNzTRhaB9+d8poenXryHNvrC9oT6C2dJyKKcRMEGYuz5RMiJmSymfw/z5wgaTVkuZLWk00x3+X\ni7C4umVfHcx7Ajnn0pTXNXyJLv7yOaAvsAZ4xsy2FS9e2yn7ZPOeQM65IitMzT++hu9GoMzMWrKX\nT5sc/AFqttcy7bFl3PfyGgDOPXIIE8cMQAryksbOudalMDX/Il7Dt0Eh1tMKlamQPYHa8nEqpBAz\nQZi5PFMyIWZKKu1r+Ja0TE+gfcu6MmXWUu5dVMWa6s0+E8g5V3TBXMO3Hm227JMruyfQoD134+pv\nek8g51yTeW+f1sR7AjnnCqTg1/DdT9IUSTdIukzS/k3P1rgQ62nFzNTUnkCldpyaKsRMEGYuz5RM\niJmSyucavscDC4ChwDpgGPCCpBOLlK0k7egJNMJ7Ajnniiefmv9iYLKZzc5aNx64zsyGFydeaZV9\nspkZM+av4lrvCeScy19Be/usBz4VT/vMrOtIdAH3siZHbFjJDv4Z3hPIOdcEBa35LwQu3LHn6K+R\nLgAq8s+VTIj1tJbOlNsT6Hu379oTyI9TMiFmgjBzeaZkQsyUVD6D/4+Bf5a0RtJ8oAo4Czi7qS8u\nqb2klyQ90NR9lILhfXtwy8RDGbTnbqx833sCOeeaL99r+Hbkk94+VcCzzentI+kC4BCgu5mdUMcm\nJV/2yZbdE6hDO3Hp0cP42vB9dnp84epqDh/Son+I7ZwLT2GneprZNjObZ2Z/ib82Z+DfFzgWuJEE\nQd3OM4Fqao2pD7/Cjc+swMzYuGUb189bwch+PdKO6ZxrBdKcOvIb4KdAvc1sQqynpZ2pQ/t2XPKV\nYUyOewLd8HQlE6+cznVzl3P2uMF079Ix1XwZaR+nuoSYCcLM5ZmSCTFTUqlMG5F0HLDWzF6Kp4vW\naebMmUyfPp3y8nIAysrKGDVqFOPHR0/JHPiWXK6oqEj19TPL3x07kHeXvcgtz1byQq3RrfxDnpo7\nl66d2geRL8TlioqKoPLkDhyh5Al1OcSfXwjjQeZ+ZWUlAJMmTdqxTUNSae8g6QrgdKAG6ALsAdxl\nZt/N2dRr/g3YuGUbv3zkVRa8uZ7qLTUM6NmNa/9ppPcEcs4VtrdP3M7hO0A/4C3gz2b2WpPjRfs8\nErjQzI6v42Ef/OuRqfGfPW4wG7fWMPnOhaxav4myrh24+pujvCeQc6WtcCd84/YOLxC1d3ifwrZ3\nqHOQz/1YHIJQMi1cXb2jxv/aS88xfeIhjN63Bxs21yTuCVRMoRynbCFmgjBzeaZkQsyUVD4nfK8E\nTjSzU83sZ2Z2KnAC8B/NCWBmT9YzzdM14PAhvXc6udu9S0eu/9Zo7wnknEvE2zu0Md4TyDlHa2/v\n4PInie+OHci0E4bTuUM77l1UxXl3LeTDrTWNP9k5VzJSbe/QmBDraa0lU5KeQC2dKW0hZoIwc3mm\nZELMlFTiwd/MXgE+A3wL+O/46zAzW1qkbK6ZvCeQc64++dT8LzSzX9ex/gIzu6rgySJe8y+A7J5A\nnTu04/JjD+CooX3SjuWcK56C9vPfaGbd61i/3sx6NiFcEj74F0jN9lqm/W0Z9y1aA8C5Rw5h4pgB\nRKdunHNtTPNP+Eo6StIEoH18P/v2A6BodYQQ62mtNVNuT6BrnlzOFY8uo2Z7va2Vip6ppYWYCcLM\n5ZmSCTFTUkl6+9xM9A68M3BT1noD3gEmFyGXK4LMTKB+ZV2ZMmsp9y6qYk31Zr86mHMlKJ+yzwwz\nO73IeXJ52adIFldV85N7FrFu0zYG7bkbV39zhPcEcq7tKGxvnxQEHa61W71hM+ffvZCV72+iV7eO\nXHXSSO8J5FzbUNiLubS0EOtpbSlTv7Ku3HTqIYwZ0JN1m7YVtCdQWzpOxRZiLs+UTIiZkgp68HfF\nt+PqYAd5TyDnSomXfRwQ9QS6df4qrvOeQM61BQWd538UUGlmKyTtA/wK2A78zMzeblbM+vng38Ie\nX7aWKbOWsrWmlrEDe/pMIOdap4LW/K8nuvIWwFVE00QN+H3+uZIJsZ7W1jMVqidQWz9OhRRiLs+U\nTIiZkspn8O9rZqviNs5HAz8EfgR8oSjJXGqG9+3BzacdyqA9u3lPIOfaqHzKPm8BhwIHAlPNbJyk\nzkT9/Is1P9DLPinynkDOtVoFLftcC8wH7iAqAUH0rv+V/HO51sBnAjnXduXT0vlXwJeBL5jZn+LV\nbwH/XIxgEGY9rdQydWjfjkuOHsY5efYEKrXj1Bwh5vJMyYSYKam85vGZ2TIzez1r+TUze7nwsVxI\nJHGGXx3MuTYln5p/GXAuMBrYPeshM7OvFCEbeM0/ONk9gQb33o3fnOQ9gZwLUEHn+T9G9EnhHmBL\n1kNmZjfV/axm88E/QN4TyLngFfSE72eBY83sOjO7MetWrIE/yHqaZ0rWE8iPU3Ih5vJMyYSYKal8\nBv9ngGGFemFJXSQ9J6lC0lJJVxZq3674fCaQc61bPmWfvYCHgb8TXcQl87HCzOzyJr241M3MNknq\nADwFXGhmT2Vt4iNJ4HJ7Ah174N78+9HDdvQE2rhlGwtXV3P4kN5pxnSu1BS07HMF0A/YC9gP+HR8\n269J0QAz2xTf7QS0B9Y1dV8uHdkzgTq1F7OWvM05d77Eh1tr2LhlG9fPW8HIfj3Sjumcy5HP4P8t\nYLSZnWxmp2ffmvriktpJqiD6JDHbzJZmPx5iPc0z1W3C0D7c8O2DKevagQVvVvPVS2/gvx5/jbPH\nDaZ7l45pxwPCOE51CTGXZ0omxExJ5dOucSWwrZAvbma1wChJPYC/ShpvZnMyj8+cOZPp06dTXl4O\nQFlZGaNGjWL8+PHAJwe+JZcrKipSff26ljPSzvPeay9xVvlW7lizB0uXb+Whv/6NQVsqOfObXwsi\nX0VFRaqvH/rPL/TlEH9+IYwHmfuVlZUATJo0acc2Dcmn5n8hcBJRm4d3sh8zsycS7aTh/f87sNnM\nfp296+bu17WsjVu28T9zXmfl+x+xqOoDOrUXv/jagd4TyLmWVdB5/pXx3V2eYGaD8ooV7a83UGNm\nGyR1Bf4K/NzMHs/edb77denJ1PjPHjeYrh3b84tHXmHW0uh9wrlHDmHimAFIjf6bdM41X+FO+JpZ\neXwblHtrYrh9gCfimv9zwAM5A/8uH4tD4Jnqt3B19Y4a/1Pz5jL12AP4wWHlQPKeQMUUynHKFWIu\nz5RMiJmSyusSTZK+Anwb6GNmx0k6FNijKWWfuCfQwfk+z4UrdzqnJM76wmCG9N6dKbOWcu+iKtZU\nb/argzkXgHzKPpOB84AbiS7duIek4cDvzeywIuXzsk8b8XJVNRd6TyDnWkpBa/4rgAlmtlLSejPr\nKak90cVcejUzaH188G9DvCeQcy2moH/ktTvwZs66TsDWfBLlI8R6mmdKpq5MSXoCtXSmEISYyzMl\nE2KmpPIZ/OcBF+esmwzMLlwc19Z5TyDnwpBP2acv8ADQG+hL9EdfG4HjzGxNkfL5iNBG5fYE+vqI\nvlz0pf139ARyzjVL4Wr+ELVjAMYAA4FVwPz4r3SLxQf/Nu7xZWuZMmspW2tqGTuwp88Ecq4wClfz\nl3ShmdWa2XNm9n9m9qyZ1Uq6oHkZ6xdiPc0zJZM004ShffjtKaPp1a0jz72xnu/fsYCq6s2pZmpp\nIebyTMmEmCmpfD5jT6ln/b8XIogrXQf17cHNpx3KoD27seK9jzjztgUsWfNB2rGca9MaLftIOoro\nI8QDwHE5Dw8BLjWzgcWJ52WfUrJxyzYuum8xz69aT+cO7bj82AO8J5BzTdP8mn/c08eAAUR1/gwj\navB2pZnd3/SMDfLBv8TUbK9l2mPLuO/laA6B9wRyrkmaX/PP9PMB7sjp6TPYzD5fxIE/yHqaZ0qm\nqZk6tG/HJUcP45wjhgBRT6ArC9QTKMTjBGHm8kzJhJgpqQanVUg6wszmxovT4xLQLgrR0tm5jMzV\nwfYt68qUWUu5Z1EVVd4TyLmCarDsI2mxmQ2P71dSTxmmGZ09G+NlnxLnPYGca5LCzvNPQdDhXMvY\nuSdQJ646aYT3BHKuYQXt7dPiQqyneaZkCplp555AHze5J1CIxwnCzOWZkgkxU1JBD/7OZXhPIOcK\ny8s+rlXJ7Qn0jRF9+TfvCeRcruaVfSSdk3X/04VI5FxzZGYCTTthOJ07tOOeRVWcd9dCPtxak3Y0\n51qVxt4uXZF1/8ViBqlLiPU0z5RMsTM1pSdQiMcJwszlmZIJMVNSjU2aXiHpv4GlQEdJ3yP6OJEp\nxwgwM7u5iBmdq1OmJ9D5dy/c0RPIZwI5l0xj8/yHAv9G1MJ5PNEFXXZhZl8sRji85u8S8J5Azu2i\noNfwfcLM6vwL3yLywd8l4j2BnNtJ4eb5m9lRkvaTNEXSDZIuk7R/8/I1LMR6mmdKpqUzJekJFOJx\ngjBzeaZkQsyUVD4XczkeWAAMBdYBw4AXJJ3YlBeW1F/SbElLJC2WdG5T9uNcRmYm0JU+E8i5RuVT\n9lkMTDaz2VnrxgPXZfr/5PXC0t7A3mZWIWl3ol8sXzezV7I287KPa5LsnkADe3XjmpNH7ugJtHHL\nNhaurubwIb1TTulc0RS0vUM/dj3h+zSwbz6JMszsbTOriO9/CLxCdGF455otMxNoQM+uvLFuE2fM\neIElaz5g45ZtXD9vBSP79Ug7onOpymfwXwhcmFlQdCbtAqCiuSEklQOjgeey14dYT/NMyYSQqV9Z\nV6ZPPJTR+/Zgw+ZtnPKLm7jkgSWcPW4w3bt0TDveDiEcq1yeKZkQMyWVT3P0HwMPSPpX4E2gP7AJ\nOL45AeKSz0zgX+NPADvMnDmT6dOnU15eDkBZWRmjRo1i/PjxwCcHviWXKyoqUn39upYzQskT2vL1\n3zqCyx5awox5/+CvjxtjBvZk4pgBPPnkk0HkywjleIW6XFFREVSeUMaDzP3KykoAJk2atGObhuTV\n20dSR+BzROWZKuBZM9uWeAd17+9B4GEzu7qOTbzm75otU+rp1qk9t86PrkTqPYFcGxduP/+4bPRH\n4H0zO7+hwdVAAAAScklEQVSezXzwd82SGfgzpZ4HFq/hl4+8Qq3B2IE9/epgrq0Kup//F4CJwBcl\nvRTfjsneIPdjcQg8UzKhZFq4unrHwD9nzhyOH74P15w8kt07t9/RE2hNIz2Bii2UY5XNMyUTYqak\nUhv8zewpM2tnZqPMbHR8eyStPK5tOnxI711O7o4t35PbvvtZBu3ZjRXvfcSk2xawZM0HKSV0Lh35\nzPNvZ2a1jW9ZUF72cUXjPYFcG1aYmr+kDsBGoMzMthYgWFI++Luiyu4JJGCy9wRybUNhav5mVgP8\nA2jRP4kMsZ7mmZJpLZmyewIZdfcESiNX2jxTMiFmSiqfaQ63Ec3zv4Zonv+Od+Vm9kShgznXUjI9\ngfqVdWXqrKXcs6iKqurNPhPItWn51Pwr47u7PMHMBhUw0067LtJ+natTdk+gwb134+qTRrBP3BPI\nuVYk3Hn+CQUdzrVNqzds5vy7F7Ly/U306tbJrw7mWqOg5/k3KsR6mmdKpjVn6lfWlZtOPYQxA3qy\nbtPH/PDPL/LEsrWp52pJnimZEDMlldfgL+krkm6W9GC8fKiklr66l3NF171LR645eSQnHrQPW2tq\nufj+xcyY/waBf1J2LrF8av6TgfOAG4GfmdkekoYDvzezw4qUz/+nuVSZGbfOX8V1c5cD3hPItRoF\nvYbvCmCCma2UtN7MekpqD7xrZr2aGbQ+Pvi7IDy+bC1TZi1la02t9wRyrUFBa/67E03xzNYJKNof\nfYVYT/NMybS1TBOG9uG3p4ymV7eOBe8J1NaOVbF4psLKZ/CfB1ycs24yMLuObZ1rczJXB/OeQK4t\nyKfs0xd4gOivfPsCK4laPhxnZmuKlM/LPi443hPItQKFnecf9+AfAwwkKgHNL3KzNx/8XZC8J5AL\nXOFq/pI6A5cDdxBdhOU24HJJXZocrxEh1tM8UzJtPVMhewK19WNVKJ6psPKp+f8W+CJRnX9M/HV8\nvN65kpPpCTTthOF07tCOexZVcd5dC/lwa03a0ZxrVD41/3XAEDNbn7WuF7DczHoWKZ+XfVyr4D2B\nXGAKOtVzDdAtZ11Xogu5O1fSfCaQa20aHPwlTZB0VNzCYQbwsKSzJH1V0g+Bh4FbixUuxHqaZ0qm\nFDPV2RPotcZ7ApXisWoKz1RYjf2J4k3sXHoR8LOc5R8BvypwLudapUxPoMxMoIvvW8zkIz/NxDH9\nfSaQC4q3dHauCLwnkEuZ9/N3Lk3eE8ilpKDz/MskXSbpHkmPZd0ebV7G+oVYT/NMyXimyIShffhd\nIz2B/Fgl45kKK5/PoHcCRwKPA3/JueUtvi7AO5JebsrznWsthvtMIBegfOb5VwN9zKwgXTwljQM+\nBG41s4Pq2czLPq7NyO4J1KGduPSYYXztwH12enzh6moOH9I7xZSujSjoPP9ngGFNz7IzM5sHrG90\nQ+faiOyrg9XUGlNnvcKNz6zEzNi4ZRvXz1vByH490o7pSkQ+g/8k4BZJ/xvX/qfEt8uKlC3Ieppn\nSsYz1S27JxDADU+v5PQrp3Pd3OWcPW4w3bt0TDlhJIRjlcszFVY+0w6uAPoBewF7FCfOzmbOnMn0\n6dMpLy8HoKysjFGjRjF+/HjgkwPfkssVFRWpvn5dyxmh5Al1uaKiIpg8Z4wdyHvLXuSWZyt5vtbo\nWv4hT82dS9dO7YPIF+JySD+/zHII40HmfmVlJQCTJk3asU1D8qn5bwSGmlnB2jlIKgce8Jq/K0Ub\nt2zjl4+8yoI311O9pYaBvbpx7ckjvSeQK4SC1vxXAtuansU5l5Gp8V96zDD+ePoYBvTsyhvrNvHd\nGS/4TCDXIvIZ/G8F7pP0nUy/n6y+P3mT9Ceik8j7S3pT0pm52+SWNULgmZLxTA1buLp6R43/HxXP\nMX3ioYzetwcbNm9L3BOomEI6VhmeqbDyqfmfQ1SGuaKOxwbl+8Jm9p18n+NcW5E7nbN7l45c/63R\n3hPItRhv7+BcQLwnkCsQ7+3jXGvkPYFcMxW0t88vJF0ef83cv1zS5c3LWL8Q62meKRnPlFxduZL0\nBGrpTGnzTIWVz2fJ/vFt3/j2WeBCYEgRcjlX8rwnkCumZpV9JB0DnGpm3y1cpJ142ceVvOyeQJ07\ntOPyYw/gqKF90o7lwlbcmr+k9sB6MyvWX/z64O8cULO9dsdMIAGTjxzCxDEDfCaQq09Ba/6Dc27D\ngV8Cq5qTsCEh1tM8UzKeKbkkubJ7AhlwzZPLufLRZdRsr00tU0vzTIWVz/SB13OWNwEVwBmFi+Oc\nq48kzhg7kH3LujJl1lLuWVRFVfVmnwnkmsSnejrXCr1cVc2F9yxi3aZtDO69G1efNMJ7Arlsza/5\nS5od381smL1TAzCzJrV4SMAHf+fqsXrDZs6/eyEr399Er26duOqkERy4T4s03HXhK0jN//b4dkd8\nux24DZgNjAQ+34yADQqxnuaZkvFMyTU1V7+yrtx06iGMGdCTdZs+jnoCLStMT6AQj5VnKqxGB38z\nuzH7BtwLfAb4CXAXsF+RMzrn6pF9dbCtNbVcfP9iZsx/g8DLuS4A+fTz70H0R12TgQeBKWa2vIjZ\nwMs+ziXiPYFcjoLU/LsB/0o08M8BLjOzJYVIl4AP/s7lwXsCuVhBav4rgQuA/wSuB/YqRD//JEKs\np3mmZDxTcoXMNWFoH36b0xOoqgk9gUI8Vp6psJK8Jcj8y/lRA9vk3c/fOVccB8U9gc6/eyEr3vuI\nM29b4DOB3C58nr9zbZT3BCpp3s/fuVKW3RMI4FzvCVQqCnoB9xYXYj3NMyXjmZIrZq7snkCQvCdQ\niMfKMxVW0IO/c675Mj2Bpp0wnM4d2nHPoirOu2shH26tSTuaS5GXfZwrIbk9gX5z0gj6ek+gtshr\n/s65nXlPoJIQds1f0jGSXpX0D0kX5T4eYj3NMyXjmZJr6VxJegKFeKw8U2GlNvjHVwG7DjgGOAD4\njqTPpJXHuVKS2xPoIu8JVHJSK/tI+jxRf6Bj4uWLAcxsWtZm/i/RuSLK7Qn0tQP35tKjh+3oCbRx\nyzYWrq7m8CG9U8n31PL3GNmvB927dNyxzjMlyhRuzV/SycDRZvaDeHkiMNbMJmdt5oO/cy3g8WVr\nueyhJXy83Tikfw8mH7kfm7bVcNdLq/nm6H7s1imd/kAffbxzhtxlz7RrpjEDekHgg/83gWMaGvzn\nzJlj48ePTyVffebMmYNnapxnSi6UXC9XVXPB3QvZsLmGD5ZXsMeQUWlH2olnSub5nx4FCQb/NNv9\nrQb6Zy33B97K3mDmzJlMnz6d8vJyAMrKyhg1atSO/yiZky0tuVxRUZHq69e1nBFKnlCXKyoqgsoT\n2s/v/dde4qzyrTy5qQ+PLYfu771Kh3bt2PszhwDw9isLAFJZ/nh7LU/Pe51+Pbqw74FjUs8D8NaS\n53m76nVGj/0Cndq3Sy0PwDuvLuDDd9dQa8acMe0SvZlI851/B2AZMAGoAuYD3zGzV7I287KPcy1o\n45ZtXD9vBad/dgAz5q/i7HGDd6pte6bwM1305aEQ8lRPM6sBzgH+CiwF/pIz8DvnWlBm8Dh73GD6\n9ujK2eMGc/28FWzcss0ztaJMiZlZsLfZs2dbaDxTMp4puVByzXv9Xftg88dm9kmmDzZ/bPNef9cz\ntaJMsUbHV7/Ej3MOoM5pit27dExt+iJ4pqSa8tre3sE559qecGv+zjnn0hP04B9i3wzPlIxnSi7E\nXJ4pmRAzJRX04O+cc644vObvnHNtj9f8nXPO7SrowT/EeppnSsYzJRdiLs+UTIiZkgp68M/0YgmJ\nZ0rGMyUXYi7PlEyImZL+Qgp68N+wYUPaEXbhmZLxTMmFmMszJRNipjYx+DvnnCuOoAf/ysrKtCPs\nwjMl45mSCzGXZ0omxExJtZ86dWraGepVVlY2NdPLPxRlZWV4psZ5puRCzOWZkgkxE0B5efnPG9sm\n9Hn+zjnniiDoso9zzrni8MHfOedKUNCDv6RfSFooqULS45L6N/6somf6L0mvxLnultQj7UwAkv5J\n0hJJ2yUdnHKWYyS9Kukfki5KM0uc52ZJ70h6Oe0sGZL6S5od/8wWSzo3gExdJD0X/39bKunKtDNl\nSGov6SVJD6SdJUNSpaRFca75aecBkFQmaWY8Ri2V9Ll6tw255i+pu5ltjO9PBkaa2T+nnOnLwONm\nVitpGoCZXZxmJgBJw4Ba4AbgJ2b2Yko52hNdm/lLwGrgeXa9NnNLZxoHfAjcamYHpZUjm6S9gb3N\nrELS7sAC4OtpHqc4Vzcz2xRfY/sp4EIzeyrNTHGuC4BDgO5mdkLaeQAkrQQOMbN1aWfJkPRH4Ekz\nuzn+Ge5mZtV1bRv0O//MwB/bHXgvrSwZZvaYmdXGi88B+6aZJ8PMXjWz19LOAXwWeN3MKs1sG/Bn\n4MQ0A5nZPGB9mhlymdnbZlYR3/8QeAXom24qMLNN8d1OQHsg9YFN0r7AscCNJGhY1sKCyRNXIcaZ\n2c0QXSe9voEfAh/8AST9h6RVwBnAtLTz5PgeMCvtEIHpB7yZtfxWvM7VQ1I5MJrozUSqJLWTVAG8\nA8w2s6VpZwJ+A/yU6JNtSAz4m6QXJP0g7TDAIOBdSbdIelHSHyR1q2/j1Ad/SY9JermO2/EAZnaJ\nmQ0AphP9I0g9U7zNJcDHZnZHS2RKmisA4dYRAxSXfGYC/xp/AkiVmdWa2SiiT7RHSBqfZh5JxwFr\nzewlAnqXHfuCmY0Gvgr8S1xeTFMH4GDgejM7GPgIqLcknfoF3M3sywk3vYMWepfdWCZJk4g+hk5o\niTwZeRyrNK0Gsk/M9yd69+9ySOoI3AXcZmb3pp0nm5lVS3oIOBSYk2KUw4ATJB0LdAH2kHSrmX03\nxUwAmNma+Ou7ku4hKnnOSzHSW8BbZvZ8vDyTBgb/1N/5N0TSflmLJwIvpZUlQ9IxRB9BTzSzLWnn\nqUea75BeAPaTVC6pE3AKcH+KeYIkScBNwFIzuzrtPACSeksqi+93Bb5Myv/nzOz/mVl/MxsEfBt4\nIoSBX1I3Sd3j+7sBXwFSnU1mZm8Db0raP171JWBJfdun/s6/EVdKGgpsB5YDP045D8C1RCfDHov+\n//J3Mzs73Ugg6RvANUBv4CFJL5nZV1s6h5nVSDoH+CvRCcObApjB8ifgSGBPSW8Cl5nZLWlmAr4A\nTAQWScoMsD8zs0dSzLQP8EdJ7YjeGM4ws8dTzFOXUMqKewH3xGNAB+B2M3s03UgATAZuj994LQfO\nrG/DoKd6OuecK46gyz7OOeeKwwd/55wrQT74O+dcCfLB3znnSpAP/s45V4J88HfOuRLkg79zzpUg\nH/yda6a4r3uLtvpwrrl88HdtlqTDJT0jaYOk9yU9JenQ+LFKSUcV6KWMcP7y1LlEQm/v4FyTSNoD\neBD4IfB/QGdgHJDpx2SE1yXSuRbj7/xdW7U/YGb2F4tsiS/Es1jSDGAA8ICkjZIuBJD0GUlzJK2P\nL62Y3cK7v6LLdq6V9J6ka+t60XgfKySd0iLfpXNN5IO/a6uWAdslTVd0TeGemQfM7HRgFXCcmXU3\ns1/H7ZUfAB4BPsUnDbL2iy9N+SCwEhhIdHGaP+e+oKJrJz8CnGNmfyny9+dcs/jg79qk+BKghxOV\nd/4ArJV0n6Q+9Tzlc0TXO50WX/5uNtGAfypRn/Z9gJ+a2WYz22pmT+c8/0jgPuB0M6v3uhOSTpD0\nNUnTJJ0maUZ8/WXnWpQP/q7Niq9rfKaZ9QeGE10jt77e+X3Z+fKTAG8QvcvfF3gj69rNuUR0buFp\nM5tbXx5JA4j69z9E1Cv/IeAvRJ9CnGtRPvi7kmBmy4A/AgdmVuVsUgX0jy+ykjGQ6OpIbwID4vJP\nnbsnGvwHSrqqgQyrzOx1SXsBG81sg5k9aGabJF0n6YgmfGvONYkP/q5NkjRU0gWS+sXL/YHvAM/G\nm7wDDMl6yrPAJuDfJHWMr117HFFtfz6wBpgWX8Gpi6TDcl5yI3AM0XVvr6wn0zBJI4kuATo3Xnds\n/PBW4O/N+Z6dy4cP/q6t2giMBZ6T9CHRwLoI+En8+JXApfHMngvMbBtwPNHFuN8FriOq378Wl3uO\nBz5NVKJ5E/hW7guaWTVROeerkn5eR6avEP1CEdBF0teJfglBdNWzgyX9j6Rezf/2nWuYX8nLuZTF\nJ6EfA843syfSzuNKg7/zdy5944DzgTPjC9/Xd27BuYLxwd+59B1IdM6hAjjAzLannMeVAC/7OOdc\nCfJ3/s45V4J88HfOuRLkg79zzpUgH/ydc64E+eDvnHMlyAd/55wrQT74O+dcCfLB3znnSpAP/s45\nV4L+Pw+K5jYInNwqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fd35a071390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    ">>> from pylab import *\n",
    ">>> xr = range(xmin, xmax+1)\n",
    ">>> # or equivalent:\n",
    ">>> xr = dpsolv.state_grid[0]\n",
    ">>> plot(xr, u, '-x')\n",
    ">>> # Annotations:\n",
    ">>> title('Optimal ordering policy')\n",
    ">>> xlabel('Stock $x_k$')\n",
    ">>> ylabel('Number of items to order $u_k$')\n",
    ">>> ylim(-0.5, u.max()+0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "IPy2 - statmath",
   "language": "",
   "name": "statmath"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

howto optimal inventory.ipynb

stodynprog-diff.py

################################################################################
# Interpolation class
# TODO : use a nicer n-dim method (like multilinear interpolation)
from scipy.interpolate import RectBivariateSpline, UnivariateSpline
from dolointerpolation.multilinear_cython import multilinear_interpolation

class UnivariateSpline(UnivariateSpline):
    '''extended UnivariateSpline class,
    where spline evaluation works uses input broadcast
    and returns an output with a coherent shape.
    '''
    #@profile
    def __call__(self, *x):
        # flatten the inputs after saving their shape:
        shape = np.array(x).shape
        # Evaluate the spline and reconstruct the dimension:
        z = super(UnivariateSpline, self).__call__(np.ravel(x))
        return z.reshape(shape)


#-----
#-----

    def interp_on_state(self, A):
        '''returns an interpolating function of matrix A, assuming that A
        is expressed on the state grid `self.state_grid`
        
        the shape of A should be (len(g) for g in self.state_grid)
        '''
        # Check the dimension of A:
        expect_shape = self._state_grid_shape
        if A.shape != expect_shape:
            raise ValueError('array `A` should be of shape {:s}, not {:s}'.format(
                             str(expect_shape), str(A.shape)) )
        
        if len(expect_shape) == 1:
            A_interp = UnivariateSpline(self.state_grid[0], A, ext=3)
            return A_interp
        elif len(expect_shape) <= 5:
            A_interp = MlinInterpolator(*self.state_grid)
            A_interp.set_values(A)
            return A_interp
            
#        if len(expect_shape) == 2:
#            x1_grid = self.state_grid[0]
#            x2_grid = self.state_grid[1]
#            A_interp = RectBivariateSplineBc(x1_grid, x2_grid, A, kx=1, ky=1)
#            return A_interp
        else:
            raise NotImplementedError('interpolation for state dimension >5'
                                      ' is not implemented.')
    # end interp_on_state()
    
    

As to adding the UnivariateSpline interpolator, it is a very useful addition, since it removes the need to compile the multilinear_cython module at least for 1D case. However, it would be better to have a way to choose which interpolator to use. What do you think?