examples/ga_k_queens.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genetic Algorithm for combinatorial optimisation\n",
    "\n",
    "This notebook showcases how to use some of the building blocks provided by EvoLP to solve a combinatorial problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Statistics\n",
    "using EvoLP\n",
    "using OrderedCollections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example,  we will solve the 8-queen problem.\n",
    "This is a constraint satisfaction problem in which the goal is to place 8 queens in a chess board such that neither of them _check_ each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![One queen on the board](img/8-queens01.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure above, we have placed a queen represented by a blue dot. All conflicting cells have been highlighted. The problem becomes harder when we add more queens to the board:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![2 queens on the board](img/8-queens02.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve this problem using a Genetic Algorithm (GA) that deals with the constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing the solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dealing with constraints\n",
    "\n",
    "To implement the solution, we need to deal with the constraints in some way.\n",
    "For this problem, we can divide the constraints in _vertical_, _horizontal_ and _diagonal_ clashes.\n",
    "\n",
    "Interestingly enough, the vertical and horizontal clashes can be dealt with by using a vector of **permutations** where the genes (index of the vector) represent a column and the alleles (values the index can take) represent the row used in that column:\n",
    "\n",
    "![A possible representation of the 8-queen problem](img/8-queens03.png)\n",
    "\n",
    "For the phenotype above, the genotype representation would look like this:\n",
    "\n",
    "```julia\n",
    "x = [1, 3, 5, 2, 6, 4, 7, 8]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EvoLP contains a **generator for integer permutations**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{verbatim}\n",
       "permutation_vector_pop(n, d, pool; replacement=false, rng=Random.GLOBAL_RNG)\n",
       "\\end{verbatim}\n",
       "Generate a population of \\texttt{n} permutation vector individuals, of size \\texttt{d} and with values sampled from \\texttt{pool}. Usually \\texttt{d} would be equal to \\texttt{length(pool)}.\n",
       "\n",
       "Sampling is \\textbf{without replacement} by default (generating permutations if \\texttt{pool} is a set). When \\texttt{replacement=true} then it generates combinations of (possibly) repeated values.\n",
       "\n",
       "\\section{Examples}\n",
       "\\begin{verbatim}\n",
       "julia> permutation_vector_pop(1, 8, 1:8)\n",
       "1-element Vector{Vector{Int64}}:\n",
       " [7, 3, 8, 1, 5, 6, 4, 2]\n",
       "\n",
       "julia> permutation_vector_pop(2, 5, [\"a\", \"b\", \"c\", \"d\", \"e\"]; replacement=false)\n",
       "2-element Vector{Vector{String}}:\n",
       " [\"e\", \"b\", \"c\", \"d\", \"a\"]\n",
       " [\"b\", \"d\", \"a\", \"e\", \"c\"]\n",
       "\\end{verbatim}\n"
      ],
      "text/markdown": [
       "```\n",
       "permutation_vector_pop(n, d, pool; replacement=false, rng=Random.GLOBAL_RNG)\n",
       "```\n",
       "\n",
       "Generate a population of `n` permutation vector individuals, of size `d` and with values sampled from `pool`. Usually `d` would be equal to `length(pool)`.\n",
       "\n",
       "Sampling is **without replacement** by default (generating permutations if `pool` is a set). When `replacement=true` then it generates combinations of (possibly) repeated values.\n",
       "\n",
       "# Examples\n",
       "\n",
       "```julia\n",
       "julia> permutation_vector_pop(1, 8, 1:8)\n",
       "1-element Vector{Vector{Int64}}:\n",
       " [7, 3, 8, 1, 5, 6, 4, 2]\n",
       "\n",
       "julia> permutation_vector_pop(2, 5, [\"a\", \"b\", \"c\", \"d\", \"e\"]; replacement=false)\n",
       "2-element Vector{Vector{String}}:\n",
       " [\"e\", \"b\", \"c\", \"d\", \"a\"]\n",
       " [\"b\", \"d\", \"a\", \"e\", \"c\"]\n",
       "```\n"
      ],
      "text/plain": [
       "\u001b[36m  permutation_vector_pop(n, d, pool; replacement=false, rng=Random.GLOBAL_RNG)\u001b[39m\n",
       "\n",
       "  Generate a population of \u001b[36mn\u001b[39m permutation vector individuals, of size \u001b[36md\u001b[39m and\n",
       "  with values sampled from \u001b[36mpool\u001b[39m. Usually \u001b[36md\u001b[39m would be equal to \u001b[36mlength(pool)\u001b[39m.\n",
       "\n",
       "  Sampling is \u001b[1mwithout replacement\u001b[22m by default (generating permutations if \u001b[36mpool\u001b[39m\n",
       "  is a set). When \u001b[36mreplacement=true\u001b[39m then it generates combinations of\n",
       "  (possibly) repeated values.\n",
       "\n",
       "\u001b[1m  Examples\u001b[22m\n",
       "\u001b[1m  ≡≡≡≡≡≡≡≡≡≡\u001b[22m\n",
       "\n",
       "\u001b[36m  julia> permutation_vector_pop(1, 8, 1:8)\u001b[39m\n",
       "\u001b[36m  1-element Vector{Vector{Int64}}:\u001b[39m\n",
       "\u001b[36m   [7, 3, 8, 1, 5, 6, 4, 2]\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> permutation_vector_pop(2, 5, [\"a\", \"b\", \"c\", \"d\", \"e\"]; replacement=false)\u001b[39m\n",
       "\u001b[36m  2-element Vector{Vector{String}}:\u001b[39m\n",
       "\u001b[36m   [\"e\", \"b\", \"c\", \"d\", \"a\"]\u001b[39m\n",
       "\u001b[36m   [\"b\", \"d\", \"a\", \"e\", \"c\"]\u001b[39m"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc permutation_vector_pop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{Vector{Int64}}:\n",
       " [2, 7, 3, 5, 1, 6, 4, 8]\n",
       " [5, 4, 1, 6, 8, 7, 2, 3]\n",
       " [6, 5, 2, 1, 4, 3, 7, 8]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pop_size = 100\n",
    "population = permutation_vector_pop(pop_size, 8, 1:8)\n",
    "first(population, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To deal with the _diagonal_ constraints, we can use the **fitness function**.\n",
    "\n",
    "The penalty of a queen is the number of queens she can check.\n",
    "The penalty of a board configuration would then be the sum of all penalties of all queens.\n",
    "\n",
    "This is what we want to **minimise**."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's build our fitness function step by step.\n",
    "\n",
    "Assume a queen $q$ is in a position $(i, j)$. Then, we can define the diagonal neighbourhood as the following:\n",
    "\n",
    "- Top-left: $(i-1, j-1)$\n",
    "- Top-right: $(i-1, j+1)$\n",
    "- Bottom-left: $(i+1, j-1)$\n",
    "- Bottom-right: $(i+1, j+1)$\n",
    "\n",
    "We can then use this information to iterate in all directions and check how many queens are there in the diagonals.\n",
    "\n",
    "If we do this **for every queen** $q$, then we will count some of the clashes twice. It is a good idea to create a set of these clashes so that we can sum them afterwards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "diag_constraints (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function diag_constraints(x)\n",
    "    # rows are values in x\n",
    "    # columns are indices from 1:8\n",
    "    fitness = []\n",
    "    for q in 1:8\n",
    "        tl = collect(zip(x[q]:-1:1, q:-1:1))\n",
    "        tr = collect(zip(x[q]:-1:1, q:1:8))\n",
    "        bl = collect(zip(x[q]:1:8, q:-1:1))\n",
    "        br = collect(zip(x[q]:1:8, q:1:8))\n",
    "\n",
    "        constraints = Set(vcat(tl, tr, bl, br))\n",
    "        delete!(constraints, (x[q], q))\n",
    "        q_fit = sum([(i, j) in constraints ? 1 : 0 for (i, j) in zip(x, 1:8)])\n",
    "        push!(fitness, q_fit)\n",
    "    end\n",
    "\n",
    "    return sum(fitness)\n",
    "end"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To handle the corners and not specify \"emtpy\" diagonals, we consider the position $(i,j)$ of a queen to count as a \"clash\" itself.\n",
    "This means that a queen in $(1,1)$ will consider $(1,1)$ as the top-left diagonal, and $(i,i), i \\in [1,8]$ as the bottom-right diagonal (again, including itself). We later remove these additional constraints via `delete!` and proceed normally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous configuration, we have the following conflicting positions:\n",
    "\n",
    "![Conflicts of the previous configuration are highlighted](img/8-queens04.png)\n",
    "\n",
    "We can test our fitness function `diag_constraints` on this board and see the number of conflicts in total:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "test = [1, 3, 5, 2, 6, 4, 7, 8]\n",
    "diag_constraints(test)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Going through each queen $q_i$ (with $i$ being the column number), we have the following conflicts:\n",
    "\n",
    "1. $q_1 = 2$\n",
    "2. $q_2 = 1$\n",
    "3. $q_3 = 0$\n",
    "4. $q_4 = 1$\n",
    "5. $q_5 = 1$\n",
    "6. $q_6 = 1$\n",
    "7. $q_7 = 2$\n",
    "8. $q_8 = 2$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evolutionary operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now need to choose our evolutionary operators, that is what we will use for **selection**, **crossover** and **mutation**.\n",
    "\n",
    "However, since we're dealing with permutations, we are restricted to use specific operators that do not violate our constraints.\n",
    "EvoLP contains some operators that can deal with permutation-based individuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Tournament parent selection with tournament size \\texttt{k}.\n",
       "\n"
      ],
      "text/markdown": [
       "Tournament parent selection with tournament size `k`.\n"
      ],
      "text/plain": [
       "  Tournament parent selection with tournament size \u001b[36mk\u001b[39m."
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc TournamentSelectionSteady"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Order 1 crossover (OX1) for permutation-based individuals.\n",
       "\n"
      ],
      "text/markdown": [
       "Order 1 crossover (OX1) for permutation-based individuals.\n"
      ],
      "text/plain": [
       "  Order 1 crossover (OX1) for permutation-based individuals."
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc OrderOneCrossover"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Swap mutation for permutation-based individuals.\n",
       "\n"
      ],
      "text/markdown": [
       "Swap mutation for permutation-based individuals.\n"
      ],
      "text/plain": [
       "  Swap mutation for permutation-based individuals."
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc SwapMutation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = TournamentSelectionSteady(5);\n",
    "C = OrderOneCrossover();\n",
    "M = SwapMutation();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logging statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `Logbook` to record statistics about our run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Logbook(LittleDict{AbstractString, Function, Vector{AbstractString}, Vector{Function}}(\"mean_eval\" => Statistics.mean, \"max_f\" => maximum, \"min_f\" => minimum, \"median_f\" => Statistics.median), NamedTuple{(:mean_eval, :max_f, :min_f, :median_f)}[])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "statnames = [\"mean_eval\", \"max_f\", \"min_f\", \"median_f\"]\n",
    "fns = [mean, maximum, minimum, median]\n",
    "thedict = LittleDict(statnames, fns)\n",
    "thelogger = Logbook(thedict)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we are ready to use all our building blocks to construct our own algorithm.\n",
    "In this case, we will use a _steady-state_ GA: instead of replacing the whole population, we will generate a fixed amount of candidate solutions and keep the best `n` individuals in the population.\n",
    "\n",
    "Let's do a summary then:\n",
    "\n",
    "- The representation is a **vector** of **permutations of integers** with values in the closed range $[1,8]$.\n",
    "- To **select** the **parents**, we use the `TournamentSelectionSteady` operator with a tournament size of 5.\n",
    "- To **recombine** the parents, we use the `OrderOneCrossover` operator.\n",
    "- To **mutate** a candidate solution, we use the `SwapMutation` operator.\n",
    "- To **select** the **survivors**, we _replace the worst_ individuals.\n",
    "- With a **population size** `pop_size` of 100.\n",
    "- With a **random initialisation** using the generator `rand_pop_int_perm`.\n",
    "- We will use a **crossover probability** of 100%.\n",
    "- And a **mutation probability** of 80%.\n",
    "\n",
    "We can then build our algorithm in a function, and use EvoLP's `Result` type for the return type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mySteadyGA (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function mySteadyGA(logbook, f, pop, k_max, S, C, M, mrate)\n",
    "    n = length(pop)\n",
    "    # Generation loop\n",
    "    runtime = @elapsed for _ in 1:k_max\n",
    "        fitnesses = f.(pop)\n",
    "        parents = select(S, fitnesses)  # this will return 2 parents\n",
    "        parents = vcat(parents, select(S, fitnesses))  # Extend the list with 2 more parents\n",
    "        offspring = [cross(C, pop[parents[1]], pop[parents[2]])]  # get first kid\n",
    "        offspring = vcat(offspring, [cross(C, pop[parents[3]], pop[parents[4]])])  # get 2nd\n",
    "        pop = vcat(pop, offspring) # add to population\n",
    "        \n",
    "        # Mutation loop\n",
    "        for i in eachindex(pop)\n",
    "            if rand() <= mrate\n",
    "                pop[i] = mutate(M, pop[i])\n",
    "            end\n",
    "        end\n",
    "        fitnesses = f.(pop)\n",
    "\n",
    "        # Log statistics\n",
    "        compute!(logbook, fitnesses)\n",
    "\n",
    "        # Find worst and remove\n",
    "        worst1 = argmax(fitnesses)\n",
    "        deleteat!(pop, worst1)\n",
    "        deleteat!(fitnesses, worst1)\n",
    "        \n",
    "        # And do the same again\n",
    "        worst2 = argmax(fitnesses)\n",
    "        deleteat!(pop, worst2)\n",
    "        deleteat!(fitnesses, worst2)\n",
    "    end\n",
    "\n",
    "    # Result reporting\n",
    "    best, best_i = findmin(f, pop)\n",
    "    n_evals = 2 * k_max * n + n\n",
    "    result = Result(best, pop[best_i], pop, k_max, n_evals, runtime)\n",
    "    return result\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To try our new algorithm, we just need to call its function with the appropriate arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "result  = mySteadyGA(thelogger, diag_constraints, population, 500, S, C, M, 0.8);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EvoLP provides convenient functions that we can use to get information about a result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimum(result) = 0\n",
      "optimizer(result) = [3, 5, 2, 8, 6, 4, 7, 1]\n",
      "f_calls(result) = 100100\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(mean_eval = 9.137254901960784, max_f = 22, min_f = 0, median_f = 8.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@show optimum(result)\n",
    "@show optimizer(result)\n",
    "@show f_calls(result)\n",
    "thelogger.records[end]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this is just a helper function to visualise our result as a chess board:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "drawboard (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function drawboard(x)\n",
    "    b = fill(\"◻\",8,8)\n",
    "    for i in 1:2:8\n",
    "        b[i,2:2:8] .= \"◼\"\n",
    "    end\n",
    "    for i in 2:2:8\n",
    "        b[i, 1:2:8] .= \"◼\"\n",
    "    end\n",
    "    for (i,j) in zip(x,1:8)\n",
    "        b[i, j] = \"♛\"\n",
    "    end\n",
    "    for i in 1:8\n",
    "        println(join(b[i,:]))\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "◻◼◻◼◻◼◻♛\n",
      "◼◻♛◻◼◻◼◻\n",
      "♛◼◻◼◻◼◻◼\n",
      "◼◻◼◻◼♛◼◻\n",
      "◻♛◻◼◻◼◻◼\n",
      "◼◻◼◻♛◻◼◻\n",
      "◻◼◻◼◻◼♛◼\n",
      "◼◻◼♛◼◻◼◻\n"
     ]
    }
   ],
   "source": [
    "drawboard(optimizer(result))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.2",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}