This tutorial will showcase how to use some of the building blocks provided by LP to solve a combinatorial problem.

In this example, we will solve the 8-queens puzzle. 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.

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:

## Implementing the solution

We can solve this problem using a **Genetic Algorithm (GA)** that deals with the constraints, using a steady-state approach and logging statistics on every iteration.

We will use the following modules:

`using Statisticsusing EvoLPusing OrderedCollections`

### Dealing with constraints

To implement the solution, we need to handle the constraints in some way. For this problem, we can divide the constraints in *vertical*, *horizontal* and *diagonal* clashes between the queens.

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:

For the phenotype above, the genotype representation would look like this:

`x = [1, 3, 5, 2, 6, 4, 7, 8]`

The queen in the first column is in row 1. The queen in the second column is in row 3, and so on.

EvoLP contains a convenient `permutation_vector_pop`

:

permutation_vector_pop

`julia> @doc permutation_vector_pop`

We can now use the generator to initialise our population:

`pop_size = 100population = permutation_vector_pop(pop_size, 8, 1:8)first(population, 3)`

`3-element Vector{Vector{Int64}}: [3, 7, 4, 8, 2, 1, 5, 6] [5, 4, 7, 6, 2, 3, 1, 8] [8, 6, 5, 1, 2, 3, 4, 7]`

To deal with the *diagonal* constraints, we can use the **fitness function**.

The penalty of a queen is the number of queens she can check. The penalty of a board configuration would then be the sum of all penalties of all queens, and this is what we want to **minimise**. So let's build our fitness function step by step.

Assume a queen $q$ is in a position $(i, j)$. Then, we can define the diagonal neighbourhood as the following:

- Top-left: $(i-1, j-1)$
- Top-right: $(i-1, j+1)$
- Bottom-left: $(i+1, j-1)$
- Bottom-right: $(i+1, j+1)$

We can then use this information to iterate in all directions and check how many queens are there in the diagonals.

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.

`function diag_constraints(x) # rows are values in x # columns are indices from 1:8 fitness = [] for q in 1:8 tl = collect(zip(x[q]:-1:1, q:-1:1)) tr = collect(zip(x[q]:-1:1, q:1:8)) bl = collect(zip(x[q]:1:8, q:-1:1)) br = collect(zip(x[q]:1:8, q:1:8)) constraints = Set(vcat(tl, tr, bl, br)) delete!(constraints, (x[q], q)) q_fit = sum([(i, j) in constraints ? 1 : 0 for (i, j) in zip(x, 1:8)]) push!(fitness, q_fit) end return sum(fitness)end`

To handle the corners and not specify "emtpy" diagonals, we consider the position $(i,j)$ of a queen to count as a "clash" itself. 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.

Using the same configuration as before, we have the following conflicting positions:

We can test our fitness function `diag_constraints`

on this board and see the number of conflicts in total:

`test = [1, 3, 5, 2, 6, 4, 7, 8]diag_constraints(test)`

`10`

Going through each queen $q_i$ (with $i$ being the column number), we have the following number of conflicts: $q_1 = 2$, $q_2 = 1$, $q_3 = 0$, $q_4 = 1$, $q_5 = 1$, $q_6 = 1$, $q_7 = 2$, $q_8 = 2$

### Evolutionary operators

We now need to choose our evolutionary operators: what we will use for **selection**, **crossover** and **mutation**.

However, since we're dealing with permutations, we are restricted to use specific operators that do not end up destroying feasible solutions and therefore violate our constraints. EvoLP contains some operators that can deal with permutation-based individuals:

`julia> @doc TournamentSelector`

`Tournament parent selection with tournament size ``k``.`

`julia> @doc OX1Recombinator`

`Order 1 crossover (OX1) for permutation-based individuals.`

`julia> @doc SwapMutator`

`Swap mutation for permutation-based individuals.`

We can now *instantiate* them and continue.

`S = TournamentSelector(5);C = OX1Recombinator();M = SwapMutator();`

### Logging statistics

We can use the Logbook to record statistics about our run:

`statnames = ["mean_eval", "max_f", "min_f", "median_f"]fns = [mean, maximum, minimum, median]thedict = LittleDict(statnames, fns)thelogger = Logbook(thedict)`

`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)}[])`

### Constructing our own algorithm

And now we are ready to use all our building blocks to construct our own algorithm. 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 each iteration.

Let's do a summary then:

- The representation is a
**vector**of**permutations of integers**with values in the closed range $[1,8]$. - To
**select**the**parents**, we use the TournamentSelector operator with a tournament size of 5. - To
**recombine**the parents, we use the OX1Recombinator operator. - To
**mutate**a candidate solution, we use the SwapMutator operator. - To
**select**the**survivors**, we*replace the worst*individuals. - With a
**population size**`pop_size`

of 100. - With a
**random initialisation**using the generator permutation_vector_pop. - We will use a
**crossover probability**of 100%. - And a
**mutation probability**of 80%.

We can then build our algorithm in a function, and use EvoLP's Result type for the return:

`function mySteadyGA(logbook, f, pop, k_max, S, C, M, mrate) n = length(pop) # Generation loop for _ in 1:k_max fitnesses = f.(pop) parents = select(S, fitnesses) # this will return 2 parents parents = vcat(parents, select(S, fitnesses)) # Extend the list with 2 more parents offspring = [cross(C, pop[parents[1]], pop[parents[2]])] # get first kid offspring = vcat(offspring, [cross(C, pop[parents[3]], pop[parents[4]])]) # get 2nd pop = vcat(pop, offspring) # add to population # Mutation loop for i in eachindex(pop) if rand() <= mrate pop[i] = mutate(M, pop[i]) end end fitnesses = f.(pop) # Log statistics compute!(logbook, fitnesses) # Find worst and remove worst1 = argmax(fitnesses) deleteat!(pop, worst1) deleteat!(fitnesses, worst1) # And do the same again worst2 = argmax(fitnesses) deleteat!(pop, worst2) deleteat!(fitnesses, worst2) end # Result reporting best, best_i = findmin(f, pop) n_evals = 2 * k_max * n + n result = Result(best, pop[best_i], pop, k_max, n_evals) return resultend`

To try our new algorithm, we just need to call its function with the appropriate arguments:

`result = mySteadyGA(thelogger, diag_constraints, population, 500, S, C, M, 0.8);`

EvoLP provides convenient functions that we can use to get information about a result:

`@show optimum(result)@show optimizer(result)@show f_calls(result)thelogger.records[end]`

`optimum(result) = 0optimizer(result) = Any[5, 1, 8, 6, 3, 7, 2, 4]f_calls(result) = 100100(mean_eval = 9.392156862745098, max_f = 20, min_f = 0, median_f = 8.0)`

And this is just a helper function to visualise our result as a chess board:

`function drawboard(x) b = fill("◻",8,8) for i in 1:2:8 b[i,2:2:8] .= "◼" end for i in 2:2:8 b[i, 1:2:8] .= "◼" end for (i,j) in zip(x,1:8) b[i, j] = "♛" end for i in 1:8 println(join(b[i,:])) endend`

`drawboard(optimizer(result))`

`◻♛◻◼◻◼◻◼◼◻◼◻◼◻♛◻◻◼◻◼♛◼◻◼◼◻◼◻◼◻◼♛♛◼◻◼◻◼◻◼◼◻◼♛◼◻◼◻◻◼◻◼◻♛◻◼◼◻♛◻◼◻◼◻`