howto optimal inventory-update.ipynb

howto optimal inventory.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, 0, 1, 1, 1, 1, 1, 0, 2])"
      ]
     },
     "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..."
     ]
    },
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-21-d911517b07fd>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mJ_0\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_x\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mJ\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdpsolv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalue_iteration\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mJ_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mJ\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdpsolv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalue_iteration\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mJ\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/deheeger/Github/stodynprog/stodynprog/stodynprog.py\u001b[0m in \u001b[0;36mvalue_iteration\u001b[1;34m(self, J_next, rel_dp, report_time)\u001b[0m\n\u001b[0;32m    510\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    511\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mind_x\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx_k\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mitertools\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mizip\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstate_ind\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstate_grid\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 512\u001b[1;33m             \u001b[0mJ_xk_opt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu_xk_opt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_value_at_state_vect\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx_k\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mJ_next_interp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    513\u001b[0m             \u001b[1;31m# Save the optimal value:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    514\u001b[0m             \u001b[0mJ_k\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mind_x\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mJ_xk_opt\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/deheeger/Github/stodynprog/stodynprog/stodynprog.py\u001b[0m in \u001b[0;36m_value_at_state_vect\u001b[1;34m(self, x_k, J_next_interp, t_k)\u001b[0m\n\u001b[0;32m    675\u001b[0m         \u001b[1;31m# Compute a grid of costs:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    676\u001b[0m         \u001b[0mg_k_grid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msys\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcost\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0msys_params\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# instant cost\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 677\u001b[1;33m         \u001b[0mJ_k_grid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mg_k_grid\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mJ_next_interp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mx_next\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# add the cost-to-go\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    678\u001b[0m         \u001b[1;31m# Expected (weighted mean) cost:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    679\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mnb_perturb\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/deheeger/Github/stodynprog/stodynprog/stodynprog.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, *x_interp)\u001b[0m\n\u001b[0;32m    277\u001b[0m         \u001b[0moutput\u001b[0m \u001b[0mshape\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mshape\u001b[0m \u001b[0mof\u001b[0m \u001b[0mbroadcasted\u001b[0m \u001b[0mcoordinate\u001b[0m \u001b[0minputs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    278\u001b[0m         '''\n\u001b[1;32m--> 279\u001b[1;33m         \u001b[1;32massert\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx_interp\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    280\u001b[0m         \u001b[1;31m# Prepare the interpolated coordinates array\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    281\u001b[0m         \u001b[0mx_mesh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbroadcast_arrays\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mx_interp\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "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": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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
}

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?