{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "234feb6e-4c52-443c-a3d6-6da8e903dd5c",
   "metadata": {},
   "source": [
    "Travelling salesman problem\n",
    "\n",
    "TSP es *casi* (multi-depot) VRPTW pero:\n",
    "- Solo 1 camion\n",
    "- Todos los ciudades (aka. almacenes pequeños) tienen al menos 1 pedido\n",
    "- Capacidad infinita\n",
    "- Sin ventanas de tiempo (aka. plazos de entrega)\n",
    "- Solo 1 deposito\n",
    "\n",
    "Cambios identificados, necesarios para adaptar el TSP a nuestro caso:\n",
    "- Tramos no conectados -> distancia grande entre ellos\n",
    "- Distancias no son euclidianas, usar \"geodistance\"\n",
    "- [...]\n",
    "\n",
    "Refs:\n",
    "\n",
    "- [scikit-opt](https://github.com/guofei9987/scikit-opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "934ae28c-ccc1-4e63-832c-4f53ceb13f50",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0fb6e4d9-d0cd-438c-b0a1-21a85a23de89",
   "metadata": {},
   "source": [
    "Differential Evolution scikit-opt example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eda106bc-2411-4ace-acf0-d2b66dc2d127",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "min f(x1, x2, x3) = x1^2 + x2^2 + x3^2\n",
    "s.t.\n",
    "    x1*x2 >= 1\n",
    "    x1*x2 <= 5\n",
    "    x2 + x3 = 1\n",
    "    0 <= x1, x2, x3 <= 5\n",
    "'''\n",
    "\n",
    "\n",
    "def obj_func(p):\n",
    "    x1, x2, x3 = p\n",
    "    return x1 ** 2 + x2 ** 2 + x3 ** 2\n",
    "\n",
    "\n",
    "constraint_eq = [\n",
    "    lambda x: 1 - x[1] - x[2]\n",
    "]\n",
    "\n",
    "# r(x1, x2, x3) >= 0\n",
    "constraint_ueq = [\n",
    "    lambda x: 1 - x[0] * x[1],\n",
    "    lambda x: x[0] * x[1] - 5\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b3dcfaf-ab6b-42b7-9218-e58ea6828605",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sko.DE import DE\n",
    "\n",
    "de = DE(func=obj_func, n_dim=3, size_pop=50, max_iter=800, lb=[0, 0, 0], ub=[5, 5, 5],\n",
    "        constraint_eq=constraint_eq, constraint_ueq=constraint_ueq)\n",
    "\n",
    "best_x, best_y = de.run()\n",
    "print('best_x:', best_x, '\\n', 'best_y:', best_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b49859f-c3a3-49a4-a605-5ec880556dc5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7563715b-c74b-4010-98c5-28c1e5b1410d",
   "metadata": {},
   "source": [
    "Genetic Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03283268-a9af-4e26-a673-4c8225bf5fcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def schaffer(p):\n",
    "    '''\n",
    "    This function has plenty of local minimum, with strong shocks\n",
    "    global minimum at (0,0) with value 0\n",
    "    https://en.wikipedia.org/wiki/Test_functions_for_optimization\n",
    "    '''\n",
    "    x1, x2 = p\n",
    "    part1 = np.square(x1) - np.square(x2)\n",
    "    part2 = np.square(x1) + np.square(x2)\n",
    "    return 0.5 + (np.square(np.sin(part1)) - 0.5) / np.square(1 + 0.001 * part2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24da5ddc-04fe-4536-b5b9-48dca9a6118c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sko.GA import GA\n",
    "\n",
    "ga = GA(func=schaffer, n_dim=2, size_pop=50, max_iter=800, prob_mut=0.001, lb=[-1, -1], ub=[1, 1], precision=1e-7)\n",
    "best_x, best_y = ga.run()\n",
    "print('best_x:', best_x, '\\n', 'best_y:', best_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e4903dc-6eff-4087-a986-1f6cf3470645",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# col: individuos,  row: iterations\n",
    "Y_history = pd.DataFrame(ga.all_history_Y)\n",
    "fig, ax = plt.subplots(2, 1)\n",
    "ax[0].plot(Y_history.index, Y_history.values, '.', color='red')\n",
    "Y_history.min(axis=1).cummin().plot(kind='line')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdfb384b-b2b5-4755-b869-3a1da68cedd9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "730c1714-7f9d-427a-9ecd-eecd416c293d",
   "metadata": {},
   "source": [
    "TSP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fb2ef00-3c16-4986-afa1-428fc43f5d36",
   "metadata": {},
   "source": [
    "\"geodistance\" (using longitude, latitude): https://stackoverflow.com/questions/31632190/measuring-geographic-distance-with-scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "740a307c-7f4d-4978-a561-2424309cc310",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import spatial\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "num_points = 5\n",
    "\n",
    "points_coordinate = np.random.rand(num_points, 2)  # generate coordinate of points\n",
    "distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')\n",
    "\n",
    "\n",
    "def cal_total_distance(routine):\n",
    "    '''The objective function. input routine, return total distance.\n",
    "    cal_total_distance(np.arange(num_points))\n",
    "    '''\n",
    "    num_points, = routine.shape\n",
    "    return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b570dd4-09fa-41d8-8d80-ba0c6ae05fa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sko.GA import GA_TSP\n",
    "\n",
    "ga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)\n",
    "best_points, best_distance = ga_tsp.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4af7656-ace9-4594-ac32-3ccd027c7811",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2)\n",
    "best_points_ = np.concatenate([best_points])\n",
    "# \"path\"\n",
    "best_points_coordinate = points_coordinate[best_points_, :]\n",
    "ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')\n",
    "ax[1].plot(ga_tsp.generation_best_Y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f361519e-2397-4b70-be35-45dd14cc15af",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40f6c413-89f8-4934-a165-1183fc5458eb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "089bbdae-063a-4a5c-8f7a-fb8bb26dc83c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7266973-639c-42f9-86b6-3369ba04d03c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5127a9e8-17a8-4a07-b825-2aaf11a6270e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sko.ACA import ACA_TSP\n",
    "\n",
    "aca = ACA_TSP(func=cal_total_distance, n_dim=num_points,\n",
    "              size_pop=50, max_iter=200,\n",
    "              distance_matrix=distance_matrix)\n",
    "\n",
    "best_x, best_y = aca.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9b939a8-3a5f-41d6-8263-25dacebc4602",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2)\n",
    "best_points_ = np.concatenate([best_x, [best_x[0]]])\n",
    "best_points_coordinate = points_coordinate[best_points_, :]\n",
    "ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')\n",
    "pd.DataFrame(aca.y_best_history).cummin().plot(ax=ax[1])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}