{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "from functools import partial\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "\n",
    "from alns import ALNS\n",
    "from alns.accept import HillClimbing\n",
    "from alns.select import RouletteWheel\n",
    "from alns.stop import MaxIterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SEED = 5432"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The cutting-stock problem\n",
    "\n",
    "The [cutting-stock problem (CSP)](https://en.wikipedia.org/wiki/Cutting_stock_problem) is well-known problem in operations research. In simple terms, it asks how many pieces of material are needed to cut an ordered amount of differently-sized beams, such that the wastage is minimised. Many different formulations exist, e.g. for the one-, two-, and even three-dimensional case. Here we solve an instance of the one-dimensional problem, obtained from the data [here](http://www.math.tu-dresden.de/~capad/cpd-ti.html#1D). It is known that the optimal solution for this problem requires the use of only 74 beams."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "OPTIMAL_BEAMS = 74"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The first line lists the number of lines for beam orders.\n",
    "# The second line is the length of the available beams. Each\n",
    "# following line is an order of (length, amount) tuples.\n",
    "with open(\"data/640.csp\") as file:\n",
    "    data = file.readlines()\n",
    "\n",
    "NUM_LINES = int(data[0])\n",
    "BEAM_LENGTH = int(data[1])\n",
    "\n",
    "# Beams to be cut from the available beams\n",
    "BEAMS = [\n",
    "    int(length)\n",
    "    for datum in data[-NUM_LINES:]\n",
    "    for length, amount in [datum.strip().split()]\n",
    "    for _ in range(int(amount))\n",
    "]\n",
    "\n",
    "print(\"Each available beam is of length:\", BEAM_LENGTH)\n",
    "print(\"Number of beams to be cut (orders):\", len(BEAMS))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the ALNS meta-heuristic, we need to have destroy and repair operators that work on a proposed solution, and a way to describe such a solution in the first place. \n",
    "Let's start with the solution state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CspState:\n",
    "    \"\"\"\n",
    "    Solution state for the CSP problem. It has two data members, assignments\n",
    "    and unassigned. Assignments is a list of lists, one for each beam in use.\n",
    "    Each entry is another list, containing the ordered beams cut from this\n",
    "    beam. Each such sublist must sum to at most BEAM_LENGTH. Unassigned is a\n",
    "    list of ordered beams that are not currently assigned to one of the\n",
    "    available beams.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, assignments, unassigned=None):\n",
    "        self.assignments = assignments\n",
    "        self.unassigned = []\n",
    "\n",
    "        if unassigned is not None:\n",
    "            self.unassigned = unassigned\n",
    "\n",
    "    def copy(self):\n",
    "        \"\"\"\n",
    "        Helper method to ensure each solution state is immutable.\n",
    "        \"\"\"\n",
    "        return CspState(\n",
    "            copy.deepcopy(self.assignments), self.unassigned.copy()\n",
    "        )\n",
    "\n",
    "    def objective(self):\n",
    "        \"\"\"\n",
    "        Computes the total number of beams in use.\n",
    "        \"\"\"\n",
    "        return len(self.assignments)\n",
    "\n",
    "    def plot(self):\n",
    "        \"\"\"\n",
    "        Helper method to plot a solution.\n",
    "        \"\"\"\n",
    "        _, ax = plt.subplots(figsize=(12, 6))\n",
    "\n",
    "        ax.barh(\n",
    "            np.arange(len(self.assignments)),\n",
    "            [sum(assignment) for assignment in self.assignments],\n",
    "            height=1,\n",
    "        )\n",
    "\n",
    "        ax.set_xlim(right=BEAM_LENGTH)\n",
    "        ax.set_yticks(np.arange(len(self.assignments), step=10))\n",
    "\n",
    "        ax.margins(x=0, y=0)\n",
    "\n",
    "        ax.set_xlabel(\"Usage\")\n",
    "        ax.set_ylabel(\"Beam (#)\")\n",
    "\n",
    "        plt.draw_if_interactive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wastage(assignment):\n",
    "    \"\"\"\n",
    "    Helper method that computes the wastage on a given beam assignment.\n",
    "    \"\"\"\n",
    "    return BEAM_LENGTH - sum(assignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Destroy operators\n",
    "\n",
    "We will consider two simple destroy operators, **random_removal** and **worst_removal**. Random removal randomly removes currently assigned beams, whereas worst removal removes those beams that are currently cut with the most waste. Both remove a fixed percentage of the current solution state, controlled by a degree of destruction parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degree_of_destruction = 0.25\n",
    "\n",
    "\n",
    "def beams_to_remove(num_beams):\n",
    "    return int(num_beams * degree_of_destruction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_removal(state, rng):\n",
    "    \"\"\"\n",
    "    Iteratively removes randomly chosen beam assignments.\n",
    "    \"\"\"\n",
    "    state = state.copy()\n",
    "\n",
    "    for _ in range(beams_to_remove(state.objective())):\n",
    "        idx = rng.integers(state.objective())\n",
    "        state.unassigned.extend(state.assignments.pop(idx))\n",
    "\n",
    "    return state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def worst_removal(state, rng):\n",
    "    \"\"\"\n",
    "    Removes beams in decreasing order of wastage, such that the\n",
    "    poorest assignments are removed first.\n",
    "    \"\"\"\n",
    "    state = state.copy()\n",
    "\n",
    "    # Sort assignments by wastage, worst first\n",
    "    state.assignments.sort(key=wastage, reverse=True)\n",
    "\n",
    "    # Removes the worst assignments\n",
    "    for _ in range(beams_to_remove(state.objective())):\n",
    "        state.unassigned.extend(state.assignments.pop(0))\n",
    "\n",
    "    return state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repair operators\n",
    "\n",
    "We define two equally simple repair operators, **greedy_insert** and **minimal_wastage**. The first considers each currently unassigned ordered beam, and finds the first beam this order may be inserted into. The second does something similar, but finds a beam where its insertion would result in the smallest beam wastage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def greedy_insert(state, rng):\n",
    "    \"\"\"\n",
    "    Inserts the unassigned beams greedily into the first fitting\n",
    "    beam. Shuffles the unassigned ordered beams before inserting.\n",
    "    \"\"\"\n",
    "    rng.shuffle(state.unassigned)\n",
    "\n",
    "    while len(state.unassigned) != 0:\n",
    "        beam = state.unassigned.pop(0)\n",
    "\n",
    "        for assignment in state.assignments:\n",
    "            if beam <= wastage(assignment):\n",
    "                assignment.append(beam)\n",
    "                break\n",
    "        else:\n",
    "            state.assignments.append([beam])\n",
    "\n",
    "    return state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minimal_wastage(state, rng):\n",
    "    \"\"\"\n",
    "    For every unassigned ordered beam, the operator determines\n",
    "    which beam would minimise that beam's waste once the ordered\n",
    "    beam is inserted.\n",
    "    \"\"\"\n",
    "\n",
    "    def insertion_cost(assignment, beam):  # helper method for min\n",
    "        if beam <= wastage(assignment):\n",
    "            return wastage(assignment) - beam\n",
    "\n",
    "        return float(\"inf\")\n",
    "\n",
    "    while len(state.unassigned) != 0:\n",
    "        beam = state.unassigned.pop(0)\n",
    "\n",
    "        assignment = min(\n",
    "            state.assignments, key=partial(insertion_cost, beam=beam)\n",
    "        )\n",
    "\n",
    "        if beam <= wastage(assignment):\n",
    "            assignment.append(beam)\n",
    "        else:\n",
    "            state.assignments.append([beam])\n",
    "\n",
    "    return state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = rnd.default_rng(SEED)\n",
    "\n",
    "state = CspState([], BEAMS.copy())\n",
    "init_sol = greedy_insert(state, rng)\n",
    "\n",
    "print(\"Initial solution has objective value:\", init_sol.objective())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_sol.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heuristic solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alns = ALNS(rng)\n",
    "\n",
    "alns.add_destroy_operator(random_removal)\n",
    "alns.add_destroy_operator(worst_removal)\n",
    "\n",
    "alns.add_repair_operator(greedy_insert)\n",
    "alns.add_repair_operator(minimal_wastage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accept = HillClimbing()\n",
    "select = RouletteWheel([3, 2, 1, 0.5], 0.8, 2, 2)\n",
    "stop = MaxIterations(5_000)\n",
    "\n",
    "result = alns.iterate(init_sol, select, accept, stop)\n",
    "solution = result.best_state\n",
    "objective = solution.objective()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = plt.subplots(figsize=(12, 6))\n",
    "result.plot_objectives(ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having obtained a reasonable solution, we now want to investigate each operator's performance. This may be done via the `plot_operator_counts()` method on the `result` object, like below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure = plt.figure(\"operator_counts\", figsize=(12, 6))\n",
    "figure.subplots_adjust(bottom=0.15, hspace=0.5)\n",
    "result.plot_operator_counts(figure, title=\"Operator diagnostics\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several conclusions follow,\n",
    "\n",
    "* The destroy operators are about equally good, so there is little to be gained changing these around.\n",
    "* The repair operators, however, are not created equal: the `minimal_wastage` operator is significantly better than the simple `greedy_insert`. This makes intuitive sense, as minimising wastage exploits the assignment structure much better than a simple insertion heuristic does. It follows that the `greedy_insert` operator is actually not all that useful, and may be dropped from the ALNS heuristic altogether."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Heuristic solution has objective value:\", solution.objective())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj = solution.objective()\n",
    "print(\n",
    "    \"Number of beams used is {0}, which is {1} more than the optimal value {2}.\".format(\n",
    "        obj, obj - OPTIMAL_BEAMS, OPTIMAL_BEAMS\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "In the code above we implement a simple ALNS-based heuristic for the CSP. Although the operators are simple and few, we do find very good results in limited time: we require just one more beam than optimal solution does, a result that is **1.35%** above optimal.\n",
    "\n",
    "This notebook offers another case where the ALNS library may be put to use to construct powerful, efficient heuristic pipelines from simple, locally greedy operators."
   ]
  }
 ],
 "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": 4
}