Taxi company car allocation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taxi company car allocation\n",
    "\n",
    "\n",
    "### Optimization Methods, MA06002, 2016\n",
    "\n",
    "* Artem Korenev\n",
    "* Andrei Kvasov\n",
    "* Anton Marin\n",
    "* Oleg Sudakov\n",
    "\n",
    "It is common knowledge, that the financial efficiency, or, equivalently, profit of a taxi company depends on how many customers its drivers will be able to serve, and what resources should be spent to achieve such result. The company should also consider its reputation as a brand, which depends on, among other factors, on how fast its cars are able to pick up clients and to drive them to their destination. \n",
    "\n",
    "Both of these metrics greatly depend on how many cars a company should employ, and where they should be allocated geographically during the day.\n",
    "\n",
    "In this problem, our team will examine historical data of Uber pickups in New York City, collected from April to September 2014.\n",
    "\n",
    "# Problem formulation\n",
    "\n",
    "The dataset contains the following information about Uber pickups in New York City: Date and time of the Uber pickup, the latitude of the Uber pickup, and the longitude of it. The pickup locations will be represented as points of a graph, where the weight of a given edge between vertices will be the distance between two points. We will be solving vehicle allocation problem, assuming that we know the orders' locations and times in advance.\n",
    "\n",
    "Let's examine elements of the problem individually:\n",
    "\n",
    "## Taxi drivers\n",
    "\n",
    "The total number of taxi drivers will be denoted as $N$. Each driver will start his working day at position $X_0^j$ at time $T_0^j$, where $j$ denotes the index of a taxi car.\n",
    "\n",
    "We will assume, that the driver's shift has a following schedule:\n",
    "\n",
    "1. The car $j$ waits for orders at a specific point $X_w^j$, starting from the time $T_w^j$;\n",
    "2. The car drives to pick up the client $k$ at the same point at time $T_p^k$;\n",
    "3. After picking up the client, car immediately drives to its destination $X_d$. It arrives at time $T_a^k$;\n",
    "4. The car drops the client at the destination point and afterwards heads to the next pickup point $X_{w+1}^j$ to wait for the order. It starts waiting for the next client at time $T_{w+1}^j$ (if the car waits for the next client at the drop-off point of the previous order, $T_{w+1}^j = T_a^k$.\n",
    "\n",
    "Obviously, $T_a^k = T_p^k + time(X_p, X_d)$, and $T_{w+1}^j = T_a^k + time(X_d, X_{w+1}^j)$, where the $time(A, B)$ denotes the time needed to get from point $A$ to point $B$. $T_p^k$ is given by the dataset.\n",
    "\n",
    "We will introduce two penalties for each car:\n",
    "\n",
    "* $c_{d}$ - downtime cost per unit of time, induced on each taxi car when it is stationary.\n",
    "* $c_{f}$ - fuel cost per unit of time, induced on each taxi car when it is en route.\n",
    "\n",
    "Therefore, if the car $j$ is assigned to order $k$, the penalty for this individual order equals to:\n",
    "\n",
    "$$P_{j, k} = time(X_w^j, X_p)c_{f} + (T_p^k - T_w^j - time(X_w^j, X_p))c_{d} + time(X_p, X_d)c_{f} + time(X_d, X_{w+1}^j)c_{f}$$\n",
    "\n",
    "This formula does not take into account the downtime at the end of the day. Given end of workshift time $T_{end}$, car index $j$ and the number of orders that this car completed during the day as $n_j$ the penalty for the downtime at the end of the day for this car is $$E_j = (T_{end} - T^j_{n_j})c_{d}$$\n",
    "\n",
    "The penalty for maintaining individual car will be denoted as $c_{a}$ ($a$ stands for auto). Then, the additional penalty for maintaining $N$ cars on a given day equals to $Nc_a$\n",
    "\n",
    "## Clients\n",
    "\n",
    "The total number of clients, or orders will be denoted as $C$. Individual order is represented by a tuple $(t^k, X_p, X_d)$, where $t^k$ denotes the minimum possible time for a pickup, $X_p$ the pickup point, and $X_d$ the destination point. \n",
    "\n",
    "A cost for a late pickup of the client will be denoted as $c_{t}$. Therefore, for individual order $k$ the penalty for the pickup delay equals to $D_k = (T_p^k - t^k)c_{t}$.\n",
    "\n",
    "## Total penalty\n",
    "\n",
    "To sum up, we will aim to minimize taxi car downtime, pickup delays and fuel costs. Thus, for a given formulation of a problem our goal is:\n",
    "\n",
    "$$\\Sigma_{j, k = 1}^{N, C}P_{j, k} + \\Sigma_{j=1}^N E_j + Nc_a + \\Sigma_{k=1}^C D_k \\rightarrow min$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chosen Methods\n",
    "\n",
    "## 1. Greedy Algorithms\n",
    "\n",
    "One of the most straightforward approaches to the problem are greedy algorithms. We can perform the greedy steps in order to make decisions upon pick ups of clients. However, it can be clearly seen almost instantly that applying such algorithms will not produce an optimal solution. However, applying such technique we can reach a sufficient approximation.\n",
    "\n",
    "### 1.1 Straightforward greedy algorithm\n",
    "\n",
    "For a greedy algorithm we choose the following strategy:\n",
    "\n",
    "1. Sort the clients by their pickup time.\n",
    "2. Iterate through a clients. For a given client $i$ assign the driver to $k$ to this client with the minimal possible value of a function: $time(X_k, X_p^i) c^f + (T_p^i - T_k - time(X_k, X_p^i)) c^d + (-T_p^i + T_k + time(X_k, X_p^i)) c_t$. Where $T_k$ is the time for $k$th driver after delivering its previous client and $X_k$ is the position where the driver $k$ delivered his or her last client. After choosing the best possible driver we assign the client to him or her and update the values $T_k$ and $X_k$ according to the data from the client's order.\n",
    "3. We repeat the 2nd step until all clients are assigned to all drivers.\n",
    "\n",
    "### 1.2 Simulated annealing algorithm\n",
    "\n",
    "Since our previous algorithm has an obvious issue of falling into local optima, one possible improvement can be to perform simmulated annealing technique.\n",
    "\n",
    "Consider that our initial solution is a random assignment of all clients to different clients and also in random order. To cover the full space of the solutions we can perform two different optimization steps:\n",
    "1. Take a client $i$ assigned and assign it to a random driver.\n",
    "2. Take a driver $k$, find two clients that are assigned to this driver. Swap the order of the delivering these clients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Greedy\n",
    "\n",
    "depot_position_x = 40.5\n",
    "depot_position_y = -74.0\n",
    "\n",
    "N_drivers = 100\n",
    "N_clients = int(times.shape[0] / 2)\n",
    "\n",
    "C_assignment = 40000\n",
    "speed_limit = 0.001\n",
    "\n",
    "\n",
    "# current drivers state\n",
    "class Solution():\n",
    "    def __init__(self):\n",
    "        self.cur_driver_pos_x = []\n",
    "        self.cur_driver_pos_y = []\n",
    "        self.cur_driver_time_finish = []\n",
    "        self.cur_driver_path = []\n",
    "        self.cur_driver_cost = []\n",
    "        self.cur_driver_assigned = []\n",
    "\n",
    "        # initializing current drivers states\n",
    "        for i in range(N_drivers):\n",
    "            self.cur_driver_pos_x.append(depot_position_x)\n",
    "            self.cur_driver_pos_y.append(depot_position_y)\n",
    "            self.cur_driver_time_finish.append(0)\n",
    "            self.cur_driver_path.append([])\n",
    "            self.cur_driver_cost.append(0)\n",
    "            self.cur_driver_assigned.append(0)\n",
    "\n",
    "\n",
    "greedy_solution = Solution()\n",
    "\n",
    "\n",
    "def assignment_cost(driver_ind):\n",
    "    return C_assignment\n",
    "\n",
    "\n",
    "def get_positions(client_ind, driver_ind):\n",
    "    client_pos = np.zeros(2)\n",
    "    driver_pos = np.zeros(2)\n",
    "    driver_pos[0] = greedy_solution.cur_driver_pos_x[driver_ind]\n",
    "    driver_pos[1] = greedy_solution.cur_driver_pos_y[driver_ind]\n",
    "    client_pos[0] = xp[client_ind]\n",
    "    client_pos[1] = yp[client_ind]\n",
    "    return client_pos, driver_pos\n",
    "\n",
    "\n",
    "def fuel_costs(client_ind, driver_ind, depot=False):\n",
    "    # warning! if changed, please change at return_all_to_depot\n",
    "    client_pos, driver_pos = get_positions(client_ind, driver_ind)\n",
    "    return np.sqrt(np.dot(driver_pos, client_pos))\n",
    "\n",
    "\n",
    "def waiting_costs_and_time_arrival(client_ind, driver_ind):\n",
    "    client_pos, driver_pos = get_positions(client_ind, driver_ind)\n",
    "    distance = np.sqrt(np.dot(driver_pos - client_pos, driver_pos - client_pos))\n",
    "    time_to_travel = 1.0 * distance / speed_limit\n",
    "    time_arrival = greedy_solution.cur_driver_time_finish[driver_ind] + time_to_travel\n",
    "    waiting_costs = max(0, np.sign(time_arrival - td[client_ind])*(time_arrival - td[client_ind]) ** 2)\n",
    "    return waiting_costs, time_arrival\n",
    "\n",
    "\n",
    "def one_client_route_time(client_ind):\n",
    "    client_start_pos = np.zeros(2)\n",
    "    client_finish_pos = np.zeros(2)\n",
    "\n",
    "    client_start_pos[0] = xp[client_ind]\n",
    "    client_start_pos[1] = yp[client_ind]\n",
    "    client_finish_pos[0] = xd[client_ind]\n",
    "    client_finish_pos[1] = yd[client_ind]\n",
    "    distance_to_travel = np.sqrt(np.dot(client_start_pos - client_finish_pos, client_start_pos - client_finish_pos))\n",
    "    time_to_travel = 1.0 * distance_to_travel / speed_limit\n",
    "    return time_to_travel\n",
    "\n",
    "\n",
    "def calculate_costs_driver(client_ind, driver_ind):\n",
    "    costs_driver = 0\n",
    "\n",
    "    # calculating assignment costs\n",
    "    assignment = 0\n",
    "    if not greedy_solution.cur_driver_assigned[driver_ind]:\n",
    "        assignment = assignment_cost(driver_ind)\n",
    "\n",
    "    # calculating fuel costs\n",
    "    fuel = fuel_costs(client_ind, driver_ind)\n",
    "\n",
    "    # calculating waiting costs\n",
    "    waiting = waiting_costs_and_time_arrival(client_ind, driver_ind)[0]\n",
    "    costs_driver = assignment + fuel + waiting\n",
    "    return costs_driver\n",
    "\n",
    "\n",
    "def calculate_costs_total_plus_driver(client_ind, driver_ind):\n",
    "    # calculating previous total costs\n",
    "    total_costs = 0\n",
    "    for i in range(N_drivers):\n",
    "        total_costs += greedy_solution.cur_driver_cost[i]\n",
    "    total_costs += calculate_costs_driver(client_ind, driver_ind)\n",
    "    return total_costs\n",
    "\n",
    "\n",
    "def assign_client(client_ind, driver_ind):\n",
    "    costs = calculate_costs_driver(client_ind, driver_ind)\n",
    "    time_arrival = waiting_costs_and_time_arrival(client_ind, driver_ind)[1]\n",
    "    time_to_travel = one_client_route_time(client_ind)\n",
    "\n",
    "    # assignment\n",
    "    greedy_solution.cur_driver_pos_x[driver_ind] = xd[client_ind]\n",
    "    greedy_solution.cur_driver_pos_y[driver_ind] = yd[client_ind]\n",
    "    greedy_solution.cur_driver_time_finish[driver_ind] = max(time_arrival, td[client_ind]) + time_to_travel\n",
    "    greedy_solution.cur_driver_path[driver_ind].append(client_ind)\n",
    "    greedy_solution.cur_driver_cost[driver_ind] += costs\n",
    "    greedy_solution.cur_driver_assigned[driver_ind] = 1\n",
    "    pass\n",
    "\n",
    "\n",
    "def total_costs():\n",
    "    costs = 0\n",
    "    for i in range(N_drivers):\n",
    "        costs += greedy_solution.cur_driver_cost[i]\n",
    "    return costs\n",
    "\n",
    "\n",
    "def return_all_to_depot():\n",
    "    for i in range(N_drivers):\n",
    "        if greedy_solution.cur_driver_assigned[i]:\n",
    "            cur_pos = np.zeros(2)\n",
    "            cur_pos[0] = greedy_solution.cur_driver_pos_x[i]\n",
    "            cur_pos[1] = greedy_solution.cur_driver_pos_y[i]\n",
    "\n",
    "            depot_pos = np.asarray([depot_position_x, depot_position_y])\n",
    "            distance_to_depot = np.sqrt(np.dot(cur_pos - depot_pos, cur_pos - depot_pos))\n",
    "            time_to_depot = 1.0 * distance_to_depot / speed_limit\n",
    "\n",
    "            greedy_solution.cur_driver_pos_x[i] = depot_position_x\n",
    "            greedy_solution.cur_driver_pos_y[i] = depot_position_y\n",
    "            greedy_solution.cur_driver_time_finish[i] += time_to_depot\n",
    "            greedy_solution.cur_driver_path[i].append(-1)\n",
    "            greedy_solution.cur_driver_cost[i] += distance_to_depot\n",
    "    pass\n",
    "\n",
    "\n",
    "for i in range(N_clients):\n",
    "    min_costs = np.inf\n",
    "    min_costs_driver = -1\n",
    "    for j in range(N_drivers):\n",
    "\n",
    "        cur_costs = calculate_costs_total_plus_driver(i, j)\n",
    "        if cur_costs < min_costs:\n",
    "            min_costs_driver = j\n",
    "            min_costs = cur_costs\n",
    "\n",
    "    # assigning client to driver\n",
    "    assert min_costs_driver != -1\n",
    "    assign_client(i, min_costs_driver)\n",
    "\n",
    "return_all_to_depot()\n",
    "\n",
    "costs = total_costs()\n",
    "print 'Costs = ', costs\n",
    "\n",
    "assigned_drivers = 0\n",
    "for i in range(N_drivers):\n",
    "    if greedy_solution.cur_driver_assigned[i]:\n",
    "        assigned_drivers += 1\n",
    "\n",
    "print (assigned_drivers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Dynamic Programming\n",
    "\n",
    "The second approach to the car allocation problem solution would be to use dynamic programming. Then, it could be solved using the same method, as speed scheduling problem, which was considered in the course. Given some optimal car assignment to orders up to the order $j$ (if the orders are sorted in time), we find the optimal value for each driver $i$ if we decide to assign him or her to order $j+1$. \n",
    "\n",
    "Given the information about the previous orders, the penalty for current assignment can easily be calculated, as we have information about geographical location of previous drop-off point, drop-off time, time and distance between previous drop-off point and current pickup point.\n",
    "\n",
    "To solve such problem, a matrix of size $N \\times C$ will be used, where $C$ denotes the number of clients, or equivalently, the number of orders, and $N$ denotes the number of drivers. Then, the final assignment will be recovered from the minimal value of column $C-1$ (numeration starts from 0) by backtracking."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Linear Programming\n",
    "\n",
    "The third approach to the problem would be to use LP-relaxation. Given the variable matrix $A$ of size $N \\times C$, where $C$ denotes the number of clients, and $N$ denotes the number of drivers, $a_{i,j} = 1$, if driver $i$ was assigned to order $j$, and $0$ otherwise.\n",
    "\n",
    "Formulating the car allocation problem will require implementation of sophisticated restrictions and weights function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Column Generation\n",
    "\n",
    "The fourth approach is to use column generation technique to get routes for individual drivers. Given their start point, and assuming, as before, that the driver's workshift starts and ends at the same point, for any subset of clients we can calculate the minimum cost $c$ for driving them to their destinations by simply rearranging the order, in which we pick them up (and subsequently deliver them to their destinations).\n",
    "\n",
    "The column generation formulation is, therefore, as follows:\n",
    "\n",
    "$$c^Tx \\rightarrow min$$\n",
    "$$ s.t.: Ax \\ge 1 \\text{ (very client has to be delivered \"at least once\")}$$\n",
    "$$\\Sigma_i x_i \\le N \\text{ (the number of cars does not exceed N)}$$\n",
    "$$ x \\in \\{0,1\\} \\text{ (we either take some route, or not)}$$\n",
    "\n",
    "This leads to the dual problem:\n",
    "\n",
    "$$ \\Sigma_j y_j + Uz \\rightarrow max $$\n",
    "$$ s.t.:A^Ty + z \\le c $$\n",
    "$$ y \\ge 0 $$\n",
    "$$ z \\le 0 $$\n",
    "\n",
    "The algorithm is as follows:\n",
    "\n",
    "0. Find initial routes \n",
    "1. Find optimal solution using found routes\n",
    "2. Find the dual variables at the optimal vertex\n",
    "3. Find a new route and add it to the used routes using values of dual variables. Go back to step 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cvxpy as cvx\n",
    "\n",
    "def compute_cost(route):\n",
    "    pass\n",
    "\n",
    "def get_new_route(y):\n",
    "    pass\n",
    "\n",
    "def find_initial_routes(C, N):\n",
    "    '''Creates route matrix A, covering all clients for a given number of taxis'''\n",
    "    if N >= C:\n",
    "        return np.identity(C)\n",
    "    else:\n",
    "        A = np.zeros((C, N))\n",
    "        mult = int(np.ceil(1.0*C/N))\n",
    "        for j in range(N):\n",
    "            for i in range(mult):\n",
    "                row = mult*j + i\n",
    "                if row > C-1:\n",
    "                    return A\n",
    "                else:\n",
    "                    A[row, j] = 1\n",
    "    return A\n",
    "\n",
    "def solve_cg(C, N, maxiters = 100000, verbose = False):\n",
    "    '''Finds a solution of the given problem using column generation\n",
    "       Initial A - identity matrix of size CxC\n",
    "       Uses a limit of iterations to get a solution in reasonable time\n",
    "       C - number of clients, N - number of taxi cars'''\n",
    "    num_routes = C\n",
    "    A = find_initial_routes(C, N)\n",
    "    costs = np.zeros(num_routes)\n",
    "    for i in range(num_routes):\n",
    "        costs[i] = compute_cost(A[:,i])\n",
    "    \n",
    "    #Forming constraints via dual problem:\n",
    "    for i in range(maxiters):\n",
    "        if verbose:\n",
    "            print \"Iteration: {}\".format(i)\n",
    "        y = cvx.Variable(C)\n",
    "        z = cvx.Variable(1)\n",
    "\n",
    "        obj = cvx.Maximize(cvx.sum_entries(y) + N*z)\n",
    "        \n",
    "        constraints = [\n",
    "            (A.T)*y + z <= costs,\n",
    "            y >= 0,\n",
    "            z <= 0\n",
    "        ]\n",
    "\n",
    "        prob = cvx.Problem(obj, constraints)\n",
    "        prob.solve()\n",
    "        \n",
    "        new_route = get_new_route(y.value)\n",
    "        #TODO: Fix route dimensions\n",
    "        if new_route in (A.T).tolist():\n",
    "            break\n",
    "        else:\n",
    "            num_routes += 1\n",
    "            costs = np.append(costs, compute_cost(new_route))\n",
    "            A = np.concatenate((A, new_route), axis = 1)\n",
    "\n",
    "    #Solving initial problem:\n",
    "    x = cvx.Variable(num_routes)\n",
    "    \n",
    "    obj = cvx.Minimize(costs.T*x)\n",
    "    \n",
    "    constraints = [\n",
    "        A*x >= 1,\n",
    "        cvx.sum_entries(x) <= N,\n",
    "        x >= 0\n",
    "    ]\n",
    "    \n",
    "    prob = cvx.Problem(obj, constraints)\n",
    "    prob.solve()\n",
    "    \n",
    "    return x.value\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [ 1.,  0.],\n",
       "       [ 1.,  0.],\n",
       "       [ 1.,  0.],\n",
       "       [ 0.,  1.],\n",
       "       [ 0.,  1.],\n",
       "       [ 0.,  1.]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = find_initial_routes(7,2)\n",
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Flow Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            Date  Hour  Minute      Lat      Lon    Base\n",
      "0     04.12.2014     0       0  40.7480 -73.9870  B02512\n",
      "1087  04.12.2014     0       0  40.6646 -73.7253  B02598\n",
      "1088  04.12.2014     0       0  40.7744 -73.9632  B02598\n",
      "1089  04.12.2014     0       0  40.7280 -73.9930  B02598\n",
      "1090  04.12.2014     0       0  40.7157 -74.0049  B02598\n",
      "\n",
      "Total clients: 9085\n",
      "Desired Times: [   0    0    0 ..., 1048 1048 1048], shape: (9085L,)\n",
      "Pickup X-coordinate: [ 40.748   40.6646  40.7744 ...,  40.7188  40.6449  40.7253], shape: (9085L,)\n",
      "Pickup Y-coordinate: [-73.987  -73.7253 -73.9632 ..., -74.0067 -73.7823 -74.0092], shape: (9085L,)\n",
      "Destination X-coordinate: [ 40.7166  40.7283  40.7606 ...,  40.7363  40.7394  40.7249], shape: (9085L,)\n",
      "Distanation Y-coordinate: [-74.0042 -73.9888 -73.9698 ..., -74.0083 -74.0087 -73.9514], shape: (9085L,)\n"
     ]
    }
   ],
   "source": [
    "# Data preparations\n",
    "\n",
    "import pandas as pd\n",
    "from dateutil.parser import parse\n",
    "\n",
    "df = pd.read_csv('data/04.12.2014.csv')\n",
    "\n",
    "# sorting by time\n",
    "df = df.sort_values(['Hour', 'Minute'])\n",
    "print df.head()\n",
    "\n",
    "# Computing times with ms\n",
    "times = df.as_matrix(columns=['Hour', 'Minute'])\n",
    "# times = 60 * 1000 * (60 * times[:, 0] + times[:, 1])\n",
    "times = (60 * times[:, 0] + times[:, 1])  # TODO: Change to MS (the line above)\n",
    "\n",
    "# Vectors of X and Y positions\n",
    "x_pos = df.as_matrix(columns=['Lat']).reshape(-1)\n",
    "y_pos = df.as_matrix(columns=['Lon']).reshape(-1)\n",
    "\n",
    "# number of clients\n",
    "n = times.shape[0] / 2\n",
    "\n",
    "# Desired times\n",
    "td = times[:n]\n",
    "\n",
    "# Pickup XY coordinates\n",
    "xp = x_pos[:n]\n",
    "yp = y_pos[:n]\n",
    "\n",
    "# Destination XY coordinates\n",
    "xd = x_pos[n:]\n",
    "yd = y_pos[n:]\n",
    "\n",
    "print\n",
    "print 'Total clients: {}'.format(n)\n",
    "print 'Desired Times: {}, shape: {}'.format(td, td.shape)\n",
    "print 'Pickup X-coordinate: {}, shape: {}'.format(xp, xp.shape)\n",
    "print 'Pickup Y-coordinate: {}, shape: {}'.format(yp, yp.shape)\n",
    "print 'Destination X-coordinate: {}, shape: {}'.format(xd, xd.shape)\n",
    "print 'Distanation Y-coordinate: {}, shape: {}'.format(yd, yd.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[        inf,  0.2837062 ,  0.07086494,  0.01598124,  0.00114018,\n",
       "         0.17659179,  0.02203361,  0.0256002 ,  0.01385208,  0.04258568],\n",
       "       [ 0.01978206,         inf,  0.05273111,  0.0042107 ,  0.02044431,\n",
       "         0.15903261,  0.00560803,  0.02067124,  0.01062121,  0.03068827],\n",
       "       [ 0.02132135,  0.26267137,         inf,  0.0400125 ,  0.0569914 ,\n",
       "         0.13711316,  0.03795629,  0.03892364,  0.04749947,  0.04739726],\n",
       "       [ 0.02191347,  0.25941983,  0.01712308,         inf,  0.05637961,\n",
       "         0.13470301,  0.03640494,  0.03987794,  0.04614109,  0.0442685 ],\n",
       "       [ 0.01524106,  0.2897863 ,  0.0454657 ,  0.02441905,         inf,\n",
       "         0.16948891,  0.0310498 ,  0.00875557,  0.0348851 ,  0.05457408],\n",
       "       [ 0.0357817 ,  0.27332984,  0.01145469,  0.05656792,  0.07262816,\n",
       "                inf,  0.05654069,  0.04975862,  0.06543921,  0.0677983 ],\n",
       "       [ 0.07626402,  0.28874946,  0.11011449,  0.05544908,  0.04010012,\n",
       "         0.19884069,         inf,  0.06631327,  0.04738143,  0.06407574],\n",
       "       [ 0.10115063,  0.24997872,  0.12732042,  0.08245635,  0.07449597,\n",
       "         0.18029756,  0.07835388,         inf,  0.07193421,  0.06756249],\n",
       "       [ 0.26470121,  0.53860984,  0.2843657 ,  0.26366481,  0.25609414,\n",
       "         0.41577784,  0.27303284,  0.24928265,         inf,  0.29768851],\n",
       "       [ 0.03950494,  0.27872705,  0.01742297,  0.0603841 ,  0.07593273,\n",
       "         0.14607375,  0.06109992,  0.0519312 ,  0.06966262,         inf]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_adj_graph(xp, yp, xd, yd):\n",
    "    A = np.fromfunction(lambda i,j: np.sqrt((xd[i]-xp[j])**2+(yd[i]-yp[j])**2), (xp.size, xp.size), dtype=int)\n",
    "    for i in range(A.shape[0]):\n",
    "        A[i,i] = np.inf\n",
    "    return A\n",
    "\n",
    "get_adj_graph(xp[:10], yp[:10], xd[:10], yd[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Miscellaneous here\n",
    "\n",
    "# TODO: function to draw solutions to debug more easily\n",
    "\n",
    "import time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Greedy algorithm soluions here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processed: 9085 clients with 400 drivers for 27.3320000172 seconds\n",
      "Assignments of the drivers for clients: [   0.    1.    2. ...,   50.  226.  144.]\n",
      "Times of taxis: [[  0.00000000e+00   3.58022346e-02   3.58022346e-02 ...,   9.21009070e+02\n",
      "    9.21009070e+02   9.21009070e+02]\n",
      " [  0.00000000e+00   0.00000000e+00   5.45758171e-01 ...,   1.01700130e+03\n",
      "    1.01700130e+03   1.01700130e+03]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00 ...,   1.04806089e+03\n",
      "    1.04806089e+03   1.04806089e+03]\n",
      " ..., \n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00 ...,   1.01001970e+03\n",
      "    1.01001970e+03   1.01001970e+03]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00 ...,   9.82007398e+02\n",
      "    9.82007398e+02   9.82007398e+02]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00 ...,   1.04713168e+03\n",
      "    1.04713168e+03   1.04713168e+03]]\n"
     ]
    }
   ],
   "source": [
    "# Dynamic programming solutions here\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "## SIMPLE IMPLEMENTATION OF THE DYNAMIC PROGRAMMING:\n",
    "## Assumptions:\n",
    "## 1. Clients are sorted by the pickup desired time\n",
    "## 2. Clients are picked up in the order of their pickup desired time\n",
    "\n",
    "\n",
    "# Data Mock so far, not real data\n",
    "\n",
    "# clients amount\n",
    "# n = 3\n",
    "\n",
    "# # time desired\n",
    "# td = np.asarray([ 10, 10, 10])\n",
    "\n",
    "# # pickup point\n",
    "# xp = np.asarray([ 10.0, 20.0, 30.0])\n",
    "# yp = np.asarray([ 10.0, 20.0, 30.0])\n",
    "\n",
    "# # destination point\n",
    "# xd = np.asarray([ 10.0, 20.0, 30.0])\n",
    "# yd = np.asarray([ 20.0, 30.0, 40.0])\n",
    "\n",
    "# drivers amount\n",
    "m = 400\n",
    "\n",
    "# driver's initial positions\n",
    "dr_x = np.ones(m) * 40.7480\n",
    "dr_y = np.ones(m) * -73.9870\n",
    "\n",
    "table = np.zeros((m, n + 1))\n",
    "\n",
    "# positions for a driver at each point of the solution\n",
    "d_x_positions = np.zeros((m, n + 1))\n",
    "d_y_positions = np.zeros((m, n + 1))\n",
    "d_x_positions[:, 0] = dr_x\n",
    "d_y_positions[:, 0] = dr_y\n",
    "\n",
    "d_times = np.zeros((m, n + 1))\n",
    "assignments = np.zeros(n)\n",
    "\n",
    "def distance(x1, y1, x2, y2):\n",
    "    return np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "for i_c in range(1, n + 1):  # iterating by clients\n",
    "    d_x_positions[:, i_c] = d_x_positions[:, i_c - 1]\n",
    "    d_y_positions[:, i_c] = d_y_positions[:, i_c - 1]\n",
    "    d_times[:, i_c] = d_times[:, i_c - 1]\n",
    "\n",
    "    # desired time of current client\n",
    "    desired_time = td[i_c - 1]\n",
    "    \n",
    "    # point of pickup of current client\n",
    "    pickup_x = xp[i_c - 1]\n",
    "    pickup_y = yp[i_c - 1]\n",
    "    \n",
    "    # min driver\n",
    "    i_min = 0\n",
    "    min_penalty = np.inf\n",
    "    for i_d in range(m):  # iterating by drivers\n",
    "        distance_to_pickup = distance(d_x_positions[i_d, i_c - 1], d_y_positions[i_d, i_c - 1], pickup_x, pickup_y)\n",
    "        lateness = np.maximum(0.0, d_times[i_d, i_c] + distance_to_pickup - desired_time) # travelling speed 1.0\n",
    "                \n",
    "        total_penalty = distance_to_pickup + lateness ** 2\n",
    "                \n",
    "        if min_penalty > total_penalty:\n",
    "            min_penalty = total_penalty\n",
    "            i_min = i_d\n",
    "    \n",
    "    # assignment index of current client\n",
    "    assignments[i_c - 1] = i_min\n",
    "    \n",
    "    # updating the position of the assigned driver to the destination point of the client\n",
    "    d_x_positions[i_min, i_c] = xd[i_c - 1]\n",
    "    d_y_positions[i_min, i_c] = yd[i_c - 1]\n",
    "    \n",
    "    # updating time of the driver after delivering the client\n",
    "    # it is max from desired time or time of pick up + time that is spent on delivering\n",
    "    d_times[i_min, i_c] = \\\n",
    "                np.maximum(td[i_c - 1], d_times[i_min, i_c - 1] + distance(d_x_positions[i_min, i_c - 1], d_y_positions[i_min, i_c - 1], pickup_x, pickup_y)) + \\\n",
    "                                      distance(pickup_x, pickup_y, xd[i_c - 1], yd[i_c - 1])\n",
    "         \n",
    "end = time.time()\n",
    "\n",
    "print 'Processed: {} clients with {} drivers for {} seconds'.format(n, m, (end - start))\n",
    "print \"Assignments of the drivers for clients: {}\".format(assignments)\n",
    "print 'Times of taxis: {}'.format(d_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [py27]",
   "language": "python",
   "name": "Python [py27]"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}