{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7abfb7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "from dataclasses import dataclass\n",
    "from typing import Optional, List, Tuple\n",
    "\n",
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "from alns import ALNS\n",
    "from alns.accept import SimulatedAnnealing\n",
    "from alns.select import AlphaUCB\n",
    "from alns.stop import MaxIterations, MaxRuntime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e92a956d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b9b14f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "SEED = 2345"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c74062fc",
   "metadata": {},
   "source": [
    "# The permutation flow shop problem\n",
    "This notebook implements ideas of [Stützle and Ruiz (2018)](https://iridia.ulb.ac.be/IridiaTrSeries/link/IridiaTr2018-006.pdf).\n",
    "\n",
    "In the permutation flow shop problem (PFSP), a set of jobs $N=\\{1, \\dots, n\\}$ has to be processed on a set of machines $M = \\{1, \\dots, m\\}$. Every job $j$ requires a processing time $p_{ij}$ on machine $i$. Moreover, all $n$ jobs have to be processed in the same order on every machine, and, the jobs follow the same machine order, starting from machine $1$ and ending at machine $m$. The objective is to find a permutation of the jobs, describing the order in which they are processed, such that the maximum completion time is minimized. This is also known as minimizing the *makespan*.\n",
    "\n",
    "In this notebook, we demonstrate how to use ALNS to solve PFSP, which is known to be NP-hard. We will also cover some advanced features of the `alns` package, including autofitting the acceptance criterion, adding local search to a repair operator, and using the `**kwargs` argument in `ALNS.iterate` to test different destroy rates.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe6fc311",
   "metadata": {},
   "source": [
    "## Data\n",
    "The [Taillard instances](http://mistic.heig-vd.ch/taillard/problemes.dir/ordonnancement.dir/ordonnancement.html) are the most used benchmark instances for the permutation flow shop problem. We will use the `tai_50_20_8` instance, which consists of 50 jobs and 20 machines.\n",
    "\n",
    "We use the [dataclass](https://docs.python.org/3/library/dataclasses.html#dataclasses.dataclass) decorator to simplify our class representation a little."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e8b657b",
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass\n",
    "class Data:\n",
    "    n_jobs: int\n",
    "    n_machines: int\n",
    "    bkv: int  # best known value\n",
    "    processing_times: np.ndarray\n",
    "\n",
    "    @classmethod\n",
    "    def from_file(cls, path):\n",
    "        with open(path, \"r\") as fi:\n",
    "            lines = fi.readlines()\n",
    "\n",
    "            n_jobs, n_machines, _, bkv, _ = [\n",
    "                int(num) for num in lines[1].split()\n",
    "            ]\n",
    "            processing_times = np.genfromtxt(lines[3:], dtype=int)\n",
    "\n",
    "            return cls(n_jobs, n_machines, bkv, processing_times)\n",
    "\n",
    "\n",
    "DATA = Data.from_file(\"data/tai50_20_8.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30a3416f-a8c5-4e5e-b808-94211c6629f8",
   "metadata": {},
   "source": [
    "Let's plot a Gantt chart of a random schedule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79f56d8c-e8c6-4436-aede-1b4aabf14c7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_completion_times(schedule):\n",
    "    \"\"\"\n",
    "    Compute the completion time for each job of the passed-in schedule.\n",
    "    \"\"\"\n",
    "    completion = np.zeros(DATA.processing_times.shape, dtype=int)\n",
    "\n",
    "    for idx, job in enumerate(schedule):\n",
    "        for machine in range(DATA.n_machines):\n",
    "            prev_job = completion[machine, schedule[idx - 1]] if idx > 0 else 0\n",
    "            prev_machine = completion[machine - 1, job] if machine > 0 else 0\n",
    "            processing = DATA.processing_times[machine, job]\n",
    "\n",
    "            completion[machine, job] = max(prev_job, prev_machine) + processing\n",
    "\n",
    "    return completion\n",
    "\n",
    "\n",
    "def compute_makespan(schedule):\n",
    "    \"\"\"\n",
    "    Returns the makespan, i.e., the maximum completion time.\n",
    "    \"\"\"\n",
    "    return compute_completion_times(schedule)[-1, schedule[-1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "835f20ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot(schedule, name):\n",
    "    \"\"\"\n",
    "    Plots a Gantt chart of the schedule for the permutation flow shop problem.\n",
    "    \"\"\"\n",
    "    n_machines, n_jobs = DATA.processing_times.shape\n",
    "\n",
    "    completion = compute_completion_times(schedule)\n",
    "    start = completion - DATA.processing_times\n",
    "\n",
    "    # Plot each job using its start and completion time\n",
    "    cmap = plt.colormaps[\"rainbow\"].resampled(n_jobs)\n",
    "    machines, length, start_job, job_colors = zip(\n",
    "        *[\n",
    "            (i, DATA.processing_times[i, j], start[i, j], cmap(j - 1))\n",
    "            for i in range(n_machines)\n",
    "            for j in range(n_jobs)\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    _, ax = plt.subplots(1, figsize=(12, 6))\n",
    "    ax.barh(machines, length, left=start_job, color=job_colors)\n",
    "\n",
    "    ax.set_title(f\"{name}\\n Makespan: {compute_makespan(schedule)}\")\n",
    "    ax.set_ylabel(f\"Machine\")\n",
    "    ax.set_xlabel(f\"Completion time\")\n",
    "    ax.set_yticks(range(DATA.n_machines))\n",
    "    ax.set_yticklabels(range(1, DATA.n_machines + 1))\n",
    "    ax.invert_yaxis()\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c98304b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(\n",
    "    rnd.choice(range(DATA.n_jobs), size=DATA.n_jobs, replace=False),\n",
    "    \"A random schedule\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90dc2e14-5785-4188-82a8-10649ab192c4",
   "metadata": {},
   "source": [
    "## Solution state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9d39e8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Solution:\n",
    "    def __init__(\n",
    "        self, schedule: List[int], unassigned: Optional[List[int]] = None\n",
    "    ):\n",
    "        self.schedule = schedule\n",
    "        self.unassigned = unassigned if unassigned is not None else []\n",
    "\n",
    "    def objective(self):\n",
    "        return compute_makespan(self.schedule)\n",
    "\n",
    "    def insert(self, job: int, idx: int):\n",
    "        self.schedule.insert(idx, job)\n",
    "\n",
    "    def opt_insert(self, job: int):\n",
    "        \"\"\"\n",
    "        Optimally insert the job in the current schedule.\n",
    "        \"\"\"\n",
    "        idcs_costs = all_insert_cost(self.schedule, job)\n",
    "        idx, _ = min(idcs_costs, key=lambda idx_cost: idx_cost[1])\n",
    "        self.insert(job, idx)\n",
    "\n",
    "    def remove(self, job: int):\n",
    "        self.schedule.remove(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77d6cb77-e0a1-4348-9541-b2409325bec4",
   "metadata": {},
   "source": [
    "A key component in ruin-and-recreate algorithms is a fast way to compute the best insertion place for an unassigned job. In the PFSP, a naive approach is to try all $O(n)$ insertion positions and compute for each the resulting makespan in $O(nm)$ time. This has total time complexity $O(n^2m)$. A more efficient way is to use [Taillard's acceleration](https://scholar.google.com/citations?view_op=view_citation&hl=nl&user=vj-4SjYAAAAJ&citation_for_view=vj-4SjYAAAAJ:IjCSPb-OGe4C), which is a labeling algorithm that reduces this entire procedure to only $O(nm)$ time. We use this to implement `all_insert_cost`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6144cb0e-284b-4abd-981c-722002925d2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def all_insert_cost(schedule: List[int], job: int) -> List[Tuple[int, float]]:\n",
    "    \"\"\"\n",
    "    Computes all partial makespans when inserting a job in the schedule.\n",
    "    O(nm) using Taillard's acceleration. Returns a list of tuples of the\n",
    "    insertion index and the resulting makespan.\n",
    "\n",
    "    [1] Taillard, E. (1990). Some efficient heuristic methods for the\n",
    "    flow shop sequencing problem. European Journal of Operational Research,\n",
    "    47(1), 65-74.\n",
    "    \"\"\"\n",
    "    k = len(schedule) + 1\n",
    "    m = DATA.processing_times.shape[0]\n",
    "    p = DATA.processing_times\n",
    "\n",
    "    # Earliest completion of schedule[j] on machine i before insertion\n",
    "    e = np.zeros((m + 1, k))\n",
    "    for j in range(k - 1):\n",
    "        for i in range(m):\n",
    "            e[i, j] = max(e[i, j - 1], e[i - 1, j]) + p[i, schedule[j]]\n",
    "\n",
    "    # Duration between starting time and final makespan\n",
    "    q = np.zeros((m + 1, k))\n",
    "    for j in range(k - 2, -1, -1):\n",
    "        for i in range(m - 1, -1, -1):\n",
    "            q[i, j] = max(q[i + 1, j], q[i, j + 1]) + p[i, schedule[j]]\n",
    "\n",
    "    # Earliest relative completion time\n",
    "    f = np.zeros((m + 1, k))\n",
    "    for l in range(k):\n",
    "        for i in range(m):\n",
    "            f[i, l] = max(f[i - 1, l], e[i, l - 1]) + p[i, job]\n",
    "\n",
    "    # Partial makespan; drop the last (dummy) row of q\n",
    "    M = np.max(f + q, axis=0)\n",
    "\n",
    "    return [(idx, M[idx]) for idx in np.argsort(M)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f9f2b2",
   "metadata": {},
   "source": [
    "## Destroy operators\n",
    "We implement two destroy operators: a random job removal operator, and an adjacent job removal operator. Both destroy operators rely on the `n_remove` parameter. We set a default value of 2 in case it's not provided. We will show later how we can experiment with different values for `n_remove`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c1ea3da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_removal(state: Solution, rng, n_remove=2) -> Solution:\n",
    "    \"\"\"\n",
    "    Randomly remove a number jobs from the solution.\n",
    "    \"\"\"\n",
    "    destroyed = deepcopy(state)\n",
    "\n",
    "    for job in rng.choice(DATA.n_jobs, n_remove, replace=False):\n",
    "        destroyed.unassigned.append(job)\n",
    "        destroyed.schedule.remove(job)\n",
    "\n",
    "    return destroyed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff11a5d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def adjacent_removal(state: Solution, rng, n_remove=2) -> Solution:\n",
    "    \"\"\"\n",
    "    Randomly remove a number adjacent jobs from the solution.\n",
    "    \"\"\"\n",
    "    destroyed = deepcopy(state)\n",
    "\n",
    "    start = rng.integers(DATA.n_jobs - n_remove)\n",
    "    jobs_to_remove = [state.schedule[start + idx] for idx in range(n_remove)]\n",
    "\n",
    "    for job in jobs_to_remove:\n",
    "        destroyed.unassigned.append(job)\n",
    "        destroyed.schedule.remove(job)\n",
    "\n",
    "    return destroyed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44207f2f",
   "metadata": {},
   "source": [
    "## Repair operator\n",
    "We implement a greedy repair operator and use a specific ordering for the unassigned jobs: the jobs with the highest total processing times are inserted first. This is also known as the *NEH ordering*, after [Nawaz, Enscore and Ham (1983)](https://www.sciencedirect.com/science/article/abs/pii/0305048383900889)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bb06bab",
   "metadata": {},
   "outputs": [],
   "source": [
    "def greedy_repair(state: Solution, rng, **kwargs) -> Solution:\n",
    "    \"\"\"\n",
    "    Greedily insert the unassigned jobs back into the schedule. The jobs are\n",
    "    inserted in non-decreasing order of total processing times.\n",
    "    \"\"\"\n",
    "    state.unassigned.sort(key=lambda j: sum(DATA.processing_times[:, j]))\n",
    "\n",
    "    while len(state.unassigned) != 0:\n",
    "        job = state.unassigned.pop()  # largest total processing time first\n",
    "        state.opt_insert(job)\n",
    "\n",
    "    return state"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "494a071c",
   "metadata": {},
   "source": [
    "## Initial solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33eb4baf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def NEH(processing_times: np.ndarray) -> Solution:\n",
    "    \"\"\"\n",
    "    Schedules jobs in decreasing order of the total processing times.\n",
    "\n",
    "    [1] Nawaz, M., Enscore Jr, E. E., & Ham, I. (1983). A heuristic algorithm\n",
    "    for the m-machine, n-job flow-shop sequencing problem. Omega, 11(1), 91-95.\n",
    "    \"\"\"\n",
    "    largest_first = np.argsort(processing_times.sum(axis=0)).tolist()[::-1]\n",
    "    solution = Solution([largest_first.pop(0)], [])\n",
    "\n",
    "    for job in largest_first:\n",
    "        solution.opt_insert(job)\n",
    "\n",
    "    return solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fc0f4d1-ffa2-4f71-ab09-cc4e6d9cc0a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "init = NEH(DATA.processing_times)\n",
    "objective = init.objective()\n",
    "pct_diff = 100 * (objective - DATA.bkv) / DATA.bkv\n",
    "\n",
    "print(f\"Initial solution objective is {objective}.\")\n",
    "print(f\"This is {pct_diff:.1f}% worse than the best known value {DATA.bkv}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3c771d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(init.schedule, \"NEH\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f2fc12",
   "metadata": {},
   "source": [
    "## Heuristic solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67376595",
   "metadata": {},
   "outputs": [],
   "source": [
    "alns = ALNS(rnd.default_rng(SEED))\n",
    "\n",
    "alns.add_destroy_operator(random_removal)\n",
    "alns.add_destroy_operator(adjacent_removal)\n",
    "\n",
    "alns.add_repair_operator(greedy_repair)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52a42a27-e982-4b22-ba98-c6cff32733e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ITERS = 8000\n",
    "\n",
    "init = NEH(DATA.processing_times)\n",
    "select = AlphaUCB(\n",
    "    scores=[5, 2, 1, 0.5],\n",
    "    alpha=0.05,\n",
    "    num_destroy=len(alns.destroy_operators),\n",
    "    num_repair=len(alns.repair_operators),\n",
    ")\n",
    "stop = MaxIterations(ITERS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7870678-289f-40ad-bea7-ab61a43956de",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Autofitting acceptance criteria\n",
    "To use the simulated annealing criterion, we need to determine three parameters: 1) the start temperature, 2) the final temperature and 3) the updating step. These parameters are often calculated using the following procedure:\n",
    "\n",
    "- Start temperature is set to a specific value, such that a first candidate solution is accepted with 50% probability if it is 5% worse than the initial solution.\n",
    "- Final temperature is set to 1.\n",
    "- The updating step is set to the linear/exponential growth rate that is needed to decrease from the start temperature to the final temperature in the specified number of iterations.\n",
    "\n",
    "Because this is such a common procedure, the `SimulatedAnnealing` class exposes the `autofit` method that can determine these parameters automatically. See the [documentation](https://alns.readthedocs.io/en/latest/accept.html#alns.accept.SimulatedAnnealing.SimulatedAnnealing.autofit) for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50d07c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "accept = SimulatedAnnealing.autofit(init.objective(), 0.05, 0.50, ITERS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ae1785-9a55-4289-848a-b419adac4b57",
   "metadata": {},
   "outputs": [],
   "source": [
    "result = alns.iterate(init, select, accept, stop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98099f7e-582e-4664-b3ed-2b32cc251dfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = plt.subplots(figsize=(12, 6))\n",
    "result.plot_objectives(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35137d49",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = result.best_state\n",
    "objective = solution.objective()\n",
    "pct_diff = 100 * (objective - DATA.bkv) / DATA.bkv\n",
    "\n",
    "print(f\"Best heuristic objective is {objective}.\")\n",
    "print(f\"This is {pct_diff:.1f}% worse than the best known value {DATA.bkv}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab663489",
   "metadata": {},
   "outputs": [],
   "source": [
    "name = \"ALNS with random & adjacent removal + greedy repair\"\n",
    "plot(result.best_state.schedule, name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6f4df42-124d-4c76-9ad9-c7c9ccfcf830",
   "metadata": {},
   "source": [
    "### Adding local search to greedy repair"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baaf095e-e377-43bd-a842-0844dd532339",
   "metadata": {},
   "source": [
    "Some ALNS heuristics include a local search after the repair step. This is also the case for the Iterated Greedy algorithm by [Ruiz and Stützle (2007)](https://www.sciencedirect.com/science/article/abs/pii/S0377221705008507?casa_token=V1QrcrUlY_QAAAAA:Vg3ADRaLM7-AaOu2DAo7P0sLatwnNF4LWEnB6K9OqjtQtYNqEW4YTnOAmIz227byWa2LEE3939s), which is a crossover between large neighborhood search and iterated local search. Iterated Greedy is one of the best performing metaheuristics for the PFSP.\n",
    "\n",
    "There are various ways to implement this. As explained in the [documentation](https://alns.readthedocs.io/en/latest/other%20single-trajectory%20heuristics.html#ils), one could implement a perturbation and local search operator and pass this to ALNS as destroy and repair operator, respectively. We do it differently: we make a new repair operator function that adds a local search step after the repair step is completed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e73142",
   "metadata": {},
   "outputs": [],
   "source": [
    "def greedy_repair_then_local_search(state: Solution, rng, **kwargs):\n",
    "    \"\"\"\n",
    "    Greedily insert the unassigned jobs back into the schedule (using NEH\n",
    "    ordering). Apply local search afterwards.\n",
    "    \"\"\"\n",
    "    state = greedy_repair(state, rng, **kwargs)\n",
    "    local_search(state, **kwargs)\n",
    "    return state\n",
    "\n",
    "\n",
    "def local_search(solution: Solution, **kwargs):\n",
    "    \"\"\"\n",
    "    Improves the current solution in-place using the insertion neighborhood.\n",
    "    A random job is selected and put in the best new position. This continues\n",
    "    until relocating any of the jobs does not lead to an improving move.\n",
    "    \"\"\"\n",
    "    improved = True\n",
    "\n",
    "    while improved:\n",
    "        improved = False\n",
    "        current = solution.objective()\n",
    "\n",
    "        for job in rnd.choice(\n",
    "            solution.schedule, len(solution.schedule), replace=False\n",
    "        ):\n",
    "            solution.remove(job)\n",
    "            solution.opt_insert(job)\n",
    "\n",
    "            if solution.objective() < current:\n",
    "                improved = True\n",
    "                current = solution.objective()\n",
    "                break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4228e899-6b60-4789-8a25-46a85fe1450a",
   "metadata": {},
   "outputs": [],
   "source": [
    "alns = ALNS(rnd.default_rng(SEED))\n",
    "\n",
    "alns.add_destroy_operator(random_removal)\n",
    "alns.add_destroy_operator(adjacent_removal)\n",
    "\n",
    "alns.add_repair_operator(greedy_repair_then_local_search)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "389f79d0-94ec-405e-b696-8a74d62d1069",
   "metadata": {},
   "outputs": [],
   "source": [
    "ITERS = 600  # fewer iterations because local search is expensive\n",
    "\n",
    "init = NEH(DATA.processing_times)\n",
    "select = AlphaUCB(\n",
    "    scores=[5, 2, 1, 0.5],\n",
    "    alpha=0.05,\n",
    "    num_destroy=len(alns.destroy_operators),\n",
    "    num_repair=len(alns.repair_operators),\n",
    ")\n",
    "accept = SimulatedAnnealing.autofit(init.objective(), 0.05, 0.50, ITERS)\n",
    "stop = MaxIterations(ITERS)\n",
    "\n",
    "result = alns.iterate(init, select, accept, stop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97de15f6-3560-4c76-a5ad-c253b79d27a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = plt.subplots(figsize=(12, 6))\n",
    "result.plot_objectives(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3374c6f-3782-4f85-8b45-9bbfdf91ea88",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = result.best_state\n",
    "objective = solution.objective()\n",
    "pct_diff = 100 * (objective - DATA.bkv) / DATA.bkv\n",
    "\n",
    "print(f\"Best heuristic objective is {objective}.\")\n",
    "print(f\"This is {pct_diff:.1f}% worse than the best known value {DATA.bkv}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95d33e4b-5ec4-4f2a-bb3c-a01be4eb11e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "name = \"ALNS with random & adjacent removal + greedy repair with local search\"\n",
    "plot(result.best_state.schedule, name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de4737b9-a3d7-4e41-87c6-2ecf14b4c87b",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Tuning the destroy rate\n",
    "Besides the selection scheme and criteria, `ALNS.iterate` can also take keyword-arguments. This can be helpful, for instance, when we want to tune parameters of the destroy/repair operators. In our case, we want to test for which values of `n_destroy` we obtain the best results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77907e0b-a63e-464d-b311-dc4c111bb400",
   "metadata": {},
   "outputs": [],
   "source": [
    "alns = ALNS(rnd.default_rng(SEED))\n",
    "\n",
    "alns.add_destroy_operator(random_removal)\n",
    "alns.add_destroy_operator(adjacent_removal)\n",
    "\n",
    "alns.add_repair_operator(greedy_repair)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f4509d1-4742-4215-836a-8fcb9a46dda5",
   "metadata": {},
   "outputs": [],
   "source": [
    "objectives = {}\n",
    "\n",
    "for n_remove in range(2, 10):\n",
    "    select = AlphaUCB(\n",
    "        scores=[5, 2, 1, 0.5],\n",
    "        alpha=0.05,\n",
    "        num_destroy=len(alns.destroy_operators),\n",
    "        num_repair=len(alns.repair_operators),\n",
    "    )\n",
    "\n",
    "    iters = 2000 / (n_remove * 2)  # higher destroy rates use less iterations\n",
    "    accept = SimulatedAnnealing.autofit(init.objective(), 0.05, 0.50, iters)\n",
    "    stop = MaxRuntime(30)\n",
    "    result = alns.iterate(init, select, accept, stop, n_remove=n_remove)\n",
    "    objectives[n_remove] = result.best_state.objective()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e57a6cc4-5e2c-40ac-8c54-a366cafd5a4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=[8, 6])\n",
    "\n",
    "ax.plot(*zip(*sorted(objectives.items())))\n",
    "ax.set_title(\"Final objective value for various n_destroy values\")\n",
    "ax.set_ylabel(\"Objective value\")\n",
    "ax.set_xlabel(\"n_destroy\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d93bf45e-90fe-41ad-b78e-7aa7fdbe4bf6",
   "metadata": {},
   "source": [
    "From this simple experiment, it looks like removing 3 works best. Our experiment is clearly too simple to draw serious conclusions from, but this example can be easily extended to include more instances and also more random seeds."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b3ced50",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "In this notebook, we implemented several variants of adaptive large neighborhood search heuristic to solve the permutation flow shop problem. We obtained a solution that is only 1-2% worse than the best known solution. Furthermore, we showed several advanced features in ALNS, including:\n",
    "\n",
    "* Autofitting acceptance criteria\n",
    "* Adding local search to repair operators\n",
    "* Using the `**kwargs` argument in `ALNS.iterate` to tune the `n_remove` parameter in destroy operators\n",
    "\n"
   ]
  }
 ],
 "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.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}