Guide / Tutorial

Installation and tests

BrkgaMpIpr can be installed using the Julia package manager. From the Julia REPL, type ] to enter the Pkg REPL mode and run

pkg> add BrkgaMpIpr

BrkgaMpIpr also provides a thorough unit testing that aims to harden and make the code ready for production environments. From Pkg REPL, just run

pkg> test BrkgaMpIpr

Although BrkgaMpIpr should work fine on Julia >= 1.2, some tests can fail. This issue occurs because BrkgaMpIpr uses the JLD package to save the population and results. JLD uses the HDF5 package, which produces slightly different binaries of different Julia versions. Although the tests may fail in those cases, BrkgaMpIpr is functional for regular usage. In the table below, you can see the testing fails due to JDL binary incompatibility.

Julia version	Windows	Linux	Mac OS X
1.2
1.3
1.4

TL;DR

The best way to keep it short is to look in the examples folder on the git repo. From main_minimal.jl, which solves the Travelling Salesman Problem (TSP). This is a trimmed copy:

using BrkgaMpIpr

include("tsp_instance.jl")
include("tsp_decoder.jl")

if length(ARGS) < 4
    println("Usage: julia main_minimal.jl <seed> <config-file> " *
            "<num-generations> <tsp-instance-file>")
    exit(1)
end

seed = parse(Int64, ARGS[1])
configuration_file = ARGS[2]
num_generations = parse(Int64, ARGS[3])
instance_file = ARGS[4]

instance = TSP_Instance(instance_file)

brkga_data, control_params = build_brkga(
    instance, tsp_decode!, MINIMIZE, seed, instance.num_nodes,
    configuration_file
)

initialize!(brkga_data)

evolve!(brkga_data, num_generations)

best_cost = get_best_fitness(brkga_data)
@show best_cost

You can identify the following basic steps:

1. Create a data structure inherited from AbstractInstance to hold your input data. This object is passed to the decoder function (example tsp_instance.jl);

2. Implement a decoder function. This function translates a chromosome (array of numbers in the interval [0,1]) to a solution for your problem. The decoder must return the solution value or cost to be used as fitness by BRKGA (example tsp_decoder.jl);

3. Load the instance and other relevant data;

4. Use build_brkga to create a BrkgaData that represents the internal state of the BRKGA-MP-IPR algorithm;

5. Use initialize! to init the BRKGA state;

6. Call evolve! to optimize;

7. Call get_best_fitness and/or get_best_chromosome to retrieve the best solution.

main_minimal.jl provides a very minimal example to understand the necessary steps to use the BRKGA-MP-IPR framework. However, main_complete.jl provides a full-featured code, handy for scientific use, such as experimentation and paper writing. This code allows fine-grained control of the optimization, shows several features of BRKGA-MP-IPR such as the resets, chromosome injection, and others. It also logs all optimization steps, creating outputs easy to be parsed. You should use this code for serious business and experimentation.

Getting started

BrkgaMpIpr is pretty simple, and you must provide one required data structure representing the problem instance, and one required decoder function to translate chromosomes to solutions.

Before you go further, please take a look at the examples folder in the git repo. We will use parts of that code in this guide. There, we solve the classical Traveling Salesman Problem. Given a set of cities and the distances between them (full weighted undirect graph), one must find a minimum-cost tour among all cities, such that each city is visited only once (i.e., find a Hamiltonian cycle of minimum cost). These are the files:

• tsp_instance.jl: contains the input data structures and helper functions;

• tsp_decoder.jl: contains the decoder function for TSP;

• greedy_tour.jl: simple heuristic that computes a greedy tour;

• config.conf: example of parameter settings;

• main_minimal.jl: minimal code useful to understand and test the framework. You should start here! Please take a look on this file before continue this tutorial;

• main_complete.jl: full-featured code, handy for scientific use, such as experimentation and paper writing. This code allows fine-grained control of the optimization, shows several features of BRKGA-MP-IPR such as the path-relinking calls, resets, chromosome injection, and others. It also logs all optimization steps, creating outputs easy to be parsed. You should use this code for serious business and experimentation;

• instances: folder containing some TSP instances for testing.

When you call main_minimal.jl or main_complete.jl without arguments, they show the usage. For example, assuming you are using a terminal:

$ julia main_minimal.jl
Usage: julia main_minimal.jl <seed> <config-file> <num-generations> <tsp-instance-file>

$ julia main_complete.jl
Usage:
  main_complete.jl -c <config_file> -s <seed> -r <stop_rule> -a <stop_arg> -t <max_time> -i <instance_file> [--no_evolution]
  main_complete.jl (-h | --help)

So, this is a possible output whe calling main_minimal.jl:

$ julia main_minimal.jl 27000001 config.conf 100 instances/brazil58.dat
Reading data...
Building BRKGA data and initializing...
Evolving 100 generations...
best_cost = 37552.0

For main_complete.jl, the output is more verbose, since we want to capture as much information as possible to do some statistical analysis. The output should be something close to this:

$ julia main_complete.jl -c config.conf -s 2700001 -r Generations -a 100 -t 60 -i instances/brazil58.dat
------------------------------------------------------
> Experiment started at 2019-02-13T18:40:11.789
> Instance: instances/brazil58.dat
> Configuration: config.conf
> Algorithm Parameters:
>  - population_size 2000
>  - elite_percentage 0.3
>  - mutants_percentage 0.15
>  - num_elite_parents 2
>  - total_parents 3
>  - bias_type LOGINVERSE
>  - num_independent_populations 3
>  - pr_number_pairs 0
>  - pr_minimum_distance 0.15
>  - pr_type PERMUTATION
>  - pr_selection BESTSOLUTION
>  - alpha_block_size 1.0
>  - pr_percentage 1.0
>  - exchange_interval 200
>  - num_exchange_indivuduals 2
>  - reset_interval 600
> Seed: 2700001
> Stop rule: GENERATIONS
> Stop argument: 100
> Maximum time (s): 60.0
> Number of parallel threads for decoding: 1
------------------------------------------------------

[18:40:11.87] Reading TSP data...
Number of nodes: 58

[18:40:11.906] Generating initial tour...
Initial cost: 30774.0

[18:40:11.909] Building BRKGA data...
New population size: 580

[18:40:12.092] Initializing BRKGA data...

[18:40:12.247] Warming up...

[18:40:12.771] Evolving...
* Iteration | Cost | CurrentTime
* 1 | 30774 | 0.03
* 34 | 30751 | 0.83
* 35 | 30507 | 0.85
* 36 | 30088 | 0.87
* 38 | 30023 | 0.93
* 39 | 29882 | 0.95
* 40 | 29665 | 0.97
* 41 | 29131 | 1.00
* 57 | 28221 | 1.38
* 66 | 28211 | 1.59
* 83 | 28200 | 2.01
* 86 | 28129 | 2.08
* 91 | 28118 | 2.19
[18:40:15.171] End of optimization

Total number of iterations: 100
Last update iteration: 91
Total optimization time: 2.40
Last update time: 2.19
Large number of iterations between improvements: 33
Total path relink time: 0.00
Total path relink calls: 0
Number of homogenities: 0
Improvements in the elite set: 0
Best individual improvements: 0

% Best tour cost: 28118
% Best tour: 22 8 1 30 13 40 25 9 32 20 53 50 4 18 44 24 58 5 27 43 12 57 23 54 55 2 41 35 10 52 51 47 49 3 48 39 29 36 17 26 19 6 28 14 37 34 56 46 15 45 33 21 11 16 42 38 31 7

Instance,Seed,NumNodes,TotalIterations,TotalTime,TotalPRTime,PRCalls,NumHomogenities,NumPRImprovElite,NumPrImprovBest,LargeOffset,LastUpdateIteration,LastUpdateTime,Cost
brazil58.dat,2700001,58,100,2.40,0.00,0,0,0,0,33,91,2.19,28118

I hope by now you got your system set up and running. Let's see the essential details on how to use the BrkgaMpIpr.

First things first: basic data structures and decoder function

All problem information must be encapsulated in a struct inherited from AbstractInstance. AbstractInstance is an empty-abstract struct required in the signature of the decoder function (described further). For example, assume we want to solve the Traveling Salesman Problem. One possible instance struct could be:

struct TSP_Instance <: AbstractInstance
    num_nodes::Int64
    distances::Array{Float64}
end

So, note that we have the number of nodes/citiesnum_nodes, and the distance matrix distances. We may need some additional code to load the instance and to deal with the distance matrix (omitted here, see tsp_instance.jl).

The second and most important requirement is the decoder function. The mandatory signature of the decoder is

decode!(chromosome::Array{Float64, 1},
        problem_instance::AbstractInstance,
        rewrite::Bool)::Float64

chromosome is a vector of numbers in the interval [0, 1] to be decoded. problem_instance is the data structure containing information about the problem. Such data is used by the decoder to build a solution. rewrite is an argument that indicates if the decoder should rewrite the chromosome, in case of local search / local improvements be performed during the decoder process. This flag is critical if you intend to use the Implicit Path Relink (details on path_relink!). The decoder must return a Float64 that is used as the fitness to rank the chromosomes. In general, fitness is the cost/value of the solution, but you may want to use it to penalize solutions that violate the problem constraints, for example.

In our TSP example, we have a very simple decoder that generates a permutation of nodes, and compute the cost of the cycle from that permutation. Note the used of function distance that returns the distance between two nodes and it is defined on tsp_instance.jl. Also note that, although we do not rewrite the chromosome in this example, we need rewrite in the function signature.

function tsp_decode!(chromosome::Array{Float64}, instance::TSP_Instance,
                     rewrite::Bool)::Float64

    permutation = Array{Tuple{Float64, Int64}}(undef, instance.num_nodes)
    for (index, key) in enumerate(chromosome)
        permutation[index] = (key, index)
    end

    sort!(permutation)

    cost = distance(instance, permutation[1][2], permutation[end][2])
    for i in 1:(instance.num_nodes - 1)
        cost += distance(instance, permutation[i][2], permutation[i + 1][2])
    end

    return cost
end

With the instance data and the decoder ready, we can build the BRKGA data structures and perform the optimization.

Building BRKGA-MP-IPR data structures

BrkgaMpIpr framework revolves over a single data structure called BrkgaData that represents the internal state of the BRKGA-MP-IPR algorithm. Since this structure has no constructor, you must build it using one of the Building functions. There are two build_brkga methods:

• load the parameters from a file:

function build_brkga(
    problem_instance::AbstractInstance,
    decode_function!::Function,
    opt_sense::Sense,
    seed::Int64,
    chromosome_size::Int64,
    configuration_file::String,
    evolutionary_mechanism_on::Bool = true
)::Tuple{BrkgaData, ExternalControlParams}

• load the parameters from a hand-made parameter object:

function build_brkga(
    problem_instance::AbstractInstance,
    decode_function!::Function,
    opt_sense::Sense,
    seed::Int64,
    chromosome_size::Int64,
    brkga_params::BrkgaParams,
    evolutionary_mechanism_on::Bool = true,
)::BrkgaData

Both methods require a problem_instance to be used in the decode_function! as commented before.

You also must indicate whether you are minimizing or maximizing through optimization Sense.

A good seed also must be provided for the (pseudo) random number generator (according to this paper). BrkgaMpIpr uses the Mersenne Twister [1] [2].

The chromosome_size also must be given. It indicates the length of each chromosome in the population. In general, this size depends on the instance and how the decoder works.

Another common argument is evolutionary_mechanism_on which is enabled by default. When disabled, no evolution is performed. The algorithm only decodes the chromosomes and ranks them. On each generation, all population is replaced excluding the best chromosome. This flag helps on implementations of simple multi-start algorithms.

As said before, the difference between the two methods is from where the algorithm's hyper-parameters come from. In the first version, the algorithm reads the BRKGA, IPR, and extra parameters from a simple text file that looks like this (see config.conf for detailed example):

population_size 2000
elite_percentage 0.30
mutants_percentage 0.15
num_elite_parents 2
total_parents 3
bias_type LOGINVERSE
num_independent_populations 3
pr_number_pairs 0
pr_minimum_distance 0.15
pr_type PERMUTATION
pr_selection BESTSOLUTION
alpha_block_size 1.0
pr_percentage 1.0
exchange_interval 200
num_exchange_indivuduals 2
reset_interval 600

When reading such file, the algorithm ignores all blank lines, and lines starting with #. During the building process, the building method creates a BrkgaParams object and a ExternalControlParams object. BrkgaParams contains all hyper-parameters regarding BRKGA and IPR methods and is stored in the brand-new BrkgaData returned by the method. ExternalControlParams are parameters that can be used outside the BRKGA-MP-IPR to control several aspects of the optimization. For instance, interval to apply path relink, reset the population, perform population migration, among others. Although their presence is required on the config file, they are not mandatory to the BRKGA-MP-IPR itself.

In the second method, we assume we already have a BrkgaParams object, and we just pass it directly to the function. Note that such param data is deep-copied inside BrkgaData.

Let's take a look in the example from main_minimal.jl:

seed = parse(Int64, ARGS[1])
configuration_file = ARGS[2]
num_generations = parse(Int64, ARGS[3])
instance_file = ARGS[4]

instance = TSP_Instance(instance_file)

brkga_data, control_params = build_brkga(
    instance, tsp_decode!, MINIMIZE, seed, instance.num_nodes,
    configuration_file
)

This code gets some arguments from the command line and loads a TSP instance. After that, it builds brkga_data. Note that we pass the instance data and the tsp_decode! function. Since we are looking for cycles of minimum cost, we ask for the algorithm MINIMIZE. The starting seed is also given. Since tsp_decode! considers each chromosome key as a node/city, the length of the chromosome must be the number of nodes, i.e., instance.num_nodes. Finally, we also pass the configuration file.

Let's take a look in a more elaborated example (main_complete.jl):

brkga_params, control_params = load_configuration(configuration_file)
...
brkga_params.population_size = min(brkga_params.population_size,
                                   10 * instance.num_nodes)
...
brkga_data = build_brkga(instance, tsp_decode!, MINIMIZE, seed,
                         instance.num_nodes, brkga_params, perform_evolution)

Here, we first load the configuration file using the helper function load_configuration. Then, we modify the population size to be the minimum between the original parameter and 10x the number of nodes (making the population size proportional to the instance size is a very common strategy used in genetic algorithms). Then, we call build_brkga using the param data instead of the configuration file (there is an additional parameter perform_evolution to tune the evolution either on or off, not necessary at this moment).

Now, we have a BrkgaData which will be used in all other functions during the optimization. Note that we can build several BrkgaData objects using different parameters, decoders, or instance data. These structures can be evolved in parallel and mixed-and-matched at your will. Each one holds a self-contained BRKGA state including populations, fitness information, and a state of the random number generator.

Initialization and Warm-start solutions

Before starting the optimization, we need to initialize BrkgaData using initialize! function. This procedure initializes the populations and others data structures of the BRKGA. If an initial population (warm start) is supplied, the initialization method completes the remaining individuals, if they do not exist. This method also performs the initial decoding of the chromosomes. Therefore, depending on the decoder implementation, this can take a while, and you may want to time such procedure. The syntax is pretty straightforward:

initialize!(brkga_data)

Warm-start solutions

One good strategy is to bootstrap the main optimization algorithm with good solutions from fast heuristics [1, 2, 3] or even from relaxations of integer linear programming models [4].

To do it, you must set these initial solutions before call initialize!. Since BRKGA-MP-IPR does not know the problem structure, you must encode the warm-start solution as chromosomes (vectors in the interval [0, 1]). In other words, you must do the inverse process that decode! does. For instance, this is a piece of code from main_complete.jl showing this process:

initial_cost, initial_tour = greedy_tour(instance)
...
keys = sort(rand(instance.num_nodes))
initial_chromosome = zeros(instance.num_nodes)
for i in 1:instance.num_nodes
    initial_chromosome[initial_tour[i]] = keys[i]
end
...
set_initial_population!(brkga_data, [initial_chromosome])
initialize!(brkga_data)

Here, we create one incumbent solution using the greedy heuristic greedy_tour() found in greedy_tour.jl. It gives us initial_tour which is a sequence of nodes to be visited. In the next four lines, we encode initial_tour. First, we create a vector of sorted random keys. Note that this is the same order that tsp_decode! uses. We then create the initial_chromosome, and fill it up with keys according to the nodes' order in initial_tour. Finally, we use set_initial_population! to assign the incumbent to the initial population. Note that initial_chromosome in between braces because set_initial_population! takes a vector of chromosomes. See its signature:

set_initial_population!(brkga_data::BrkgaData,
                        chromosomes::Array{Array{Float64, 1}, 1})

Indeed, you can have as much warm-start solutions as you like, limited to the size of the population. Just remember:

Optimization time: evolving the population

Once all data is set up, it is time to evolve the population and perform other operations like path-relinking, shaking, migration, and others. The call is pretty simple:

evolve!(brkga_data::BrkgaData, num_generations::Int64 = 1)

evolve! evolves all populations for num_generations.

For example, in main_minimal.jl, we just evolve the population for a given number of generations directly and then extract the best solution cost.

evolve!(brkga_data, num_generations)
best_cost = get_best_fitness(brkga_data)

On main_complete.jl, we have fine-grained control on the optimization. There, we have a main loop that evolves the population one generation at a time and performs several operations as to hold the best solution, to check whether it is time for path relink, population reset, among others. The advantage of that code is that we can track all optimization details.

Accessing solutions/chromosomes

Since Julia does not offer encapsulation mechanisms to keep data private within data structures, you can access all chromosomes, fitness, and other data members directly from BrkgaData. However, we do not recommend that, unless you are sure what you are doing. So, BrkgaMpIpr offers some helper functions.

Usually, we want to access the best chromosome after some iterations. You can use the companion functions:

get_best_fitness(brkga_data::BrkgaData)::Float64

get_best_chromosome(brkga_data::BrkgaData)::Array{Float64, 1}

get_best_fitness returns the value/fitness of the best chromosome across all populations.

get_best_chromosome returns a copy of the best chromosome across all populations. You may want to extract an actual solution from such chromosome, i.e., to apply a decoding function that returns the actual solution instead only its value.

You may also want to get a copy of specific chromosome for a given population using get_chromosome.

get_chromosome(brkga_data::BrkgaData,
               population_index::Int64,
               position::Int64)::Array{Float64, 1}

For example, you can get the 3rd best chromosome from the 2nd population using

third_best = get_chromosome(brkga_data, 2, 3)

Now, suppose you get such chromosome or chromosomes and apply a quick local search procedure on them. It may be useful to reinsert such new solutions in the BRKGA population for the next evolutionary cycles. You can do that using inject_chromosome!.

inject_chromosome!(brkga_data::BrkgaData,
                   chromosome::Array{Float64, 1},
                   population_index::Int64,
                   position::Int64,
                   fitness::Float64 = Inf)

Note that the chromosome is put in a specific position of a given population. If you do not provide the fitness, inject_chromosome! will decode the injected chromosome. For example, the following code injects a random chromosome keys into the population #1 in the last position (population_size), i.e., it will replace the worst solution:

keys = sort(rand(instance.num_nodes))
inject_chromosome!(brkga_data, keys, 1, brkga_data.params.population_size)

Implicit Path Relink

The Implicit Path Relinking (IPR) is a nice addition to the standard BRKGA framework, and it provides an excellent way to create hybrid heuristics and push the optimization further. The good thing about IPR is that you do not need to worry about the path relink implementation, which can be long and tedious if done by hand or customized per problem.

BrkgaMpIpr provides a friendly interface to use IPR directly from the BRKGA population, and you only must provide a few functions and arguments to have a Path Relink algorithm ready to go. This is the main signature of path_relink!

path_relink!(brkga_data::BrkgaData,
             pr_type::PathRelinkingType,
             pr_selection::PathRelinkingSelection,
             compute_distance::Function,
             affect_solution::Function,
             number_pairs::Int64,
             minimum_distance::Float64,
             block_size::Int64 = 1,
             max_time::Float64 = 0.0,
             percentage::Float64 = 1.0
)::PathRelinkingResult

The first argument is the BrkgaData as usual. The 2nd argument defines the type of implicit path relink to be performed (PathRelinkingType). The DIRECT path relink exchanges the keys of two chromosomes directly, and it is usually more suitable to or threshold representations, i.e., where the key values are used to some kind of discretization, such as " if x < 0.5, then 0, otherwise 1." The PERMUTATION path relink switches the order of a key according to its position in the other chromosome. Usually, this kind of path relink is more suitable to permutation representations, where the chromosome induces an order or permutation. For example, chromosome [0.4, 0.7, 0.1] may induce the increasing order (3, 1, 2). More details about threshold and permutation representations in this paper.

PathRelinkingSelection defines how the algorithm picks the chromosomes for relinking. BESTSOLUTION selects, in the order, the best solution of each population. RANDOMELITE chooses uniformly random solutions from the elite sets.

The next argument is a function to compute the distance between two chromosomes such signature must be

compute_distance(vector1::Array{Float64, 1},
                 vector2::Array{Float64, 1})::Float64

If the value returned by compute_distance() is greater than or equal to minimum_distance, the algorithm will perform the path relink between the two chromosomes. Otherwise, it will look for another pair of chromosomes. The algorithm will try number_pairs chromosomes before gives up. In the presence of multiple populations, the path relinking is performed between elite chromosomes from different populations, in a circular fashion. For example, suppose we have 3 populations. The framework performs 3 path relinkings: the first between individuals from populations 1 and 2, the second between populations 2 and 3, and the third between populations 3 and 1. In the case of just one population, both base and guiding individuals are sampled from the elite set of that population.

Note that in traditional path relink algorithms, compute_distance() depends on the problem structure. On IPR, you can use a generic distance function, or provide one that incorporates more knowledge about the problem. BrkgaMpIpr provides a function to compute the (modified) Hamming distance for threshold representations (hamming_distance), and a function that computes the Kendall Tau distance distance for permutation representations (kendall_tau_distance). Again, details about threshold and permutation representations in this paper.

As a simple example, suppose you are using a threshold representation where each chromosome key can represent one of 3 different values (a ternary threshold representation). So, one possible way to compute the distance between two chromosomes can be:

function value(key::Float64)::Float64
    return key < 0.33 ? 0.0 : (key < 0.66 ? 1.0 : 2.0)
end

function compute_distance(vector1::Array{Float64, 1},
                          vector2::Array{Float64, 1})::Float64
    total = 0.0
    for i in 1:length(vector1)
        total += abs(value(vector1[i]) - value(vector2[i]))
    end
    return total
end

To avoid changes that do not lead to new solutions, we must verify if such key exchanges affect the solution. For that, we must pass a function with the signature:

affect_solution(block1::SubArray{Float64, 1},
                block2::SubArray{Float64, 1})::Bool

affect_solution two gets partial chromosomes/block of genes block1 and block2 and checks whether changing the keys from block1 to block2 affects the solution. For instance, suppose that the alleles/keys are used as threshold such that values > 0.5 activate a feature. Suppose we have block1 = [0.3, 0.4, 0.1] and block2 = [0.4, 0.1, 0.2]. Since all values are below 0.5, changing the keys from block1 to block2 do not change the solution, and therefore, we can drop such change (and subsequently decoding). The blocks can hold only one key/allele, sequential key blocks, or even the whole chromosome. Note that affect_solution is crucial to the IPR performance since this function helps to avoid exploring regions already surveyed. Also, note that affect_solution can incorporate some problem knowledge.

(x, y) -> true

block_size defines the number of keys / size of the chromosome block to be exchanged during the direct path relink. This parameter is also critical for IPR performance since it avoids too many exchanges during the path building. Usually, we can compute this number based on the size of the chromosome by some factor (alpha_block_size in the configuration file), chosen by you. Again, details here.

The last two parameters are stopping criteria. The algorithm stops either when max_time seconds is reached or percentage% of the path is built.

Let's see the example on main_complete.jl. Remember, since we are solving the TSP, we want to use the permutation-based IPR, and the Kendall Tau distance functions.

result = path_relink!(
    brkga_data,
    brkga_params.pr_type,
    brkga_params.pr_selection,
    kendall_tau_distance,
    affect_solution_kendall_tau,
    brkga_params.pr_number_pairs,
    brkga_params.pr_minimum_distance,
    1, #block_size doesn't not matter for permutation.
    maximum_time - (time() - start_time),
    brkga_params.pr_percentage
)

Note that most parameters come from BrkgaParams. The maximum IPR time is set to the remaining time for optimization (global maximum_time minus the elapsed time time() - start_time.

path_relink! returns a PathRelinkingResult object which defines the status of the IPR optimization. These status are described on PathRelinkingResult.

If the found solution is the best solution found so far, IPR replaces the worst solution by it. Otherwise, IPR computes the distance between the found solution and all other solutions in the elite set, and replaces the worst solution by it if and only if the found solution is, at least, minimum_distance from all them.

Important notes about IPR

The API will call decode!() function always with writeback = false. The reason is that if the decoder rewrites the chromosome, the path between solutions is lost and inadvertent results may come up. Note that at the end of the path relinking, the method calls the decoder with writeback = true in the best chromosome found to guarantee that this chromosome is re-written to reflect the best solution found.

BrkgaMpIpr path_relink! implementation is multi-threaded. Instead of to build and decode each chromosome one at a time, the method builds a list of candidates, altering the alleles/keys according to the guide solution, and then decode all candidates in parallel. Note that $O(chromosome\_size^2 ~/~ block\_size)$ additional memory is necessary to build the candidates, which can be costly if the chromosome_size is very large.

Shaking and Resetting

Sometimes, BRKGA gets stuck, converging to local maxima/minima, for several iterations. When such a situation happens, it is a good idea to perturb the population, or even restart from a new one completely new. BrkgaMpIpr offers shake! function, an improved variation of the original version proposed in this paper.

shake!(brkga_data::BrkgaData,
       intensity::Int64,
       shaking_type::ShakingType,
       population_index::Int64 = Inf64)

shake! function gets an intensity parameter that measures how many times the perturbation is applied on the elite set for a given population_index (if not given, all populations are shaken). This method offers two generic/implicit ShakingTypes. With CHANGE, direct modifications are done in the keys/alleles. This kind of shaking is recommended when the chromosome uses direct or threshold representations. SWAP exchanges keys/alleles inducing new permutations. For representational definitions, please read this paper. For instance, the following code shakes all populations using 10 swap moves.

shake!(brkga_data, 10, SWAP)

Sometimes, even shaking the populations does not help to escape from local maxima/minima. So, we need a drastic measure, restarting from scratch the role population. This can be easily accomplished with reset!.

reset!(brkga_data)

Multi-population and migration

Multi-population or island model was introduced in genetic algorithms in this paper. The idea is to evolve parallel and independent populations and, once a while, exchange individuals among these populations. In several scenarios, this approach is very beneficial for optimization.

BrkgaMpIpr is implemented using such island idea from the core. If you read the guide until here, you may notice that several methods take into account multiple populations. To use multiple populations, you must set BrkgaParams.num_independent_populations with 2 ou more populations, and build BrkgaData from such parameters.

The immigration process is implemented by

exchange_elite!(brkga_data::BrkgaData, num_immigrants::Int64)

exchange_elite! copies num_immigrants from one population to another, replacing the worst num_immigrants individuals from the recipient population. Note that the migration is done for all pairs of populations. For instance, the following code exchanges 3 best individuals from each population:

exchange_elite!(brkga_data, 3)

Simulating the standard BRKGA

Sometimes, it is a good idea to test how the standard BRKGA algorithm performs for a problem. You can use BrkgaMpIpr framework to quickly implement and test a standard BRKGA.

First, you must guarantee that, during the crossover, the algorithm chooses only one elite individual and only one non-elite individual. This is easily accomplished setting num_elite_parents = 1 and total_parents = 2. Then, you must set up a bias function that ranks the elite and no-elite individual according to the original BRKGA bias parameter $\rho$ (rho).

You can use set_bias_custom_function! for that task. The given function receives the index of the chromosome and returns a ranking for it. Such ranking is used in the roulette method to choose the individual from which each allele comes to build the new chromosome. Since we have one two individuals for crossover in the standard BRKGA, the bias function must return the probability to one or other individual. In the following code, we do that with a simple if...else lambda function.

# create brkga_params by hand or reading from a file,
# then set the following by hand.
brkga_params.num_elite_parents = 1
brkga_params.total_parents = 2

rho = 0.75
set_bias_custom_function!(brkga_data, x -> x == 1 ? rho : 1.0 - rho)
initialize!(brkga_data)

Here, we first set the num_elite_parents = 1 and total_parents = 2 as explained before. Following, we set a variable rho = 0.75. This is the $\rho$ from standard BRKGA, and you may set it as you wish. Then, we set the bias function as a very simple lambda function:

x -> x == 1 ? rho : 1.0 - rho

So, if the index of the chromosome is 1 (elite individual), it gets a 0.75 rank/probability. If the index is 2 (non-elite individual), the chromosome gets 0.25 rank/probability.

Reading and writing parameters

Although we can build a BrkgaData by set up a BrkgaParams object manually, the easiest way to do so is to read such parameters from a configuration file. For this, we can use load_configuration that reads a simple plain text file and returns a tuple of BrkgaParams and ExternalControlParams objects. For instance,

brkga_params, control_params = load_configuration("tuned_conf.txt")

The configuration file must be plain text such that contains pairs of parameter name and value. This file must list all fields from BrkgaParams and ExternalControlParams, even though you do not use each one. In examples folder we have config.conf that looks like this:

population_size 2000
elite_percentage 0.30
mutants_percentage 0.15
num_elite_parents 2
total_parents 3
bias_type LOGINVERSE
num_independent_populations 3
pr_number_pairs 0
pr_minimum_distance 0.15
pr_type PERMUTATION
pr_selection BESTSOLUTION
alpha_block_size 1.0
pr_percentage 1.0
exchange_interval 200
num_exchange_indivuduals 2
reset_interval 600

It does not matter whether we use lower or upper cases. Blank lines and lines starting with # are ignored. The order of the parameters should not matter either. And, finally, this file should be readble for both C++, Julia, and Python framework versions.

In some cases, you define some of the parameters at the running time, and you may want to save them for debug or posterior use. To do so, you can use write_configuration, call upon a BrkgaParams object or BrkgaData object. For example,

write_configuration("my_new_parameters.conf", brkga_params, external_params)
# or
write_configuration("crazy_parameters.txt", brkga_data)

(Probable Valuable) Tips

Algorithm warmup

When using Julia code, it is an excellent idea to dry-run all functions you may use and, mainly, the ones you want to time. The reason is that Julia uses lazy evaluation when live-compiling the code, i.e., it compiles as it goes. Another advantage is the memory location effects of our data (principle of locality), that can be brought closer to the processor (L2/L3 caches) during the running. Obliviously, this depends on how you implement and use your data structures.

In main_complete.jl, we have the following piece of code to warmup mainly the decoder and other functions. Note that we just deep-copy brkga_data, and then, we may lose the principle of locality.

bogus_data = deepcopy(brkga_data)
evolve!(bogus_data, 2)
path_relink!(bogus_data, brkga_params.pr_type, brkga_params.pr_selection,
             (x, y) -> 1.0, (x, y) -> true, 0, 0.5, 1, 10.0, 1.0)
get_best_fitness(brkga_data)
get_best_chromosome(brkga_data)
bogus_data = nothing

Complex decoders and timing

Some problems require complex decoders while for others, the decoder contains local search procedures, that can be time-consuming. In general, the decoding is the most time-expensive component of a BRKGA algorithm, and it may skew some stopping criteria based on running time. Therefore, if your decoder is time-consuming, it is a good idea to implement a timer or chronometer kind of thing inside the decoder.

Testing for stopping time uses several CPU cycles, and you need to be careful when/where to test it, otherwise, you spend all the optimization time doing system calls to the clock.

IMHO, the most effective way to do it is to test time at the very end of the decoding. If the current time is larger than the maximum time allowed, simple return Inf or -Inf according to your optimization direction. In this way, we make the solution invalid since it violates the maximum time allowed. The BRKGA framework takes care of the rest.

Multi-threading

Since Moore's law is not holding its status anymore, we, simple mortals, must appeal to the wonders of multi-threading. This paradigm can be tricky to code, and Amdahl's law plays against us. Several genetic algorithms, and in particular, BRKGA, can use parallel solution evaluation (or decoding), which makes the use of multi-threading relatively straightforward. BrkgaMpIpr is not different, and it uses Julia multi-threading capabilities to do so.

First, as commented several times in this guide, the decoder must be THREAD-SAFE. So, each thread must have its own read/write data structures and may share other read-only data. The simplest way to do it is to create those structures inside the decoder (like most people do). But be aware, this strategy slows down the algorithm significantly depending on the size and format of the structures, and I do not recommend it.

IMHO, the best way to do that is to preallocate the data structure per thread (using Threads.nthreads()), and pass them to the decoder through the problem instance. Then, inside the decoder, you can use Threads.threadid() and recover the memory you want to use.

Let's see a simple example considering the TSP example. tsp_decode! uses a single array to create the permutation of nodes. Let's pre-allocate its memory per thread. So, in TSP_Instance, we pre-allocate copies of such array, one for each thread:

using Base.Threads

# Declare as a type to make the code shorter and more readable.
PermutationArray = Array{Tuple{Float64, Int64}, 1}

struct TSP_Instance <: AbstractInstance
    num_nodes::Int64
    distances::Array{Float64}

    # Permutations arrays per thread, to be pre-allocated.
    permutation_per_thread::Array{PermutationArray, 1}

    function TSP_Instance(filename::String)
        #... Code for loading here

        # Allocate the main array to create references for other arrays per thread.
        permutation_per_thread = Array{PermutationArray, 1}(undef, nthreads())

        # Pre-allocate the permutation arrays, one for each thread.
        for i in 1:nthreads()
            permutation_per_thread[i] = PermutationArray(undef, num_nodes)
        end

        new(num_nodes, distances, permutation_per_thread)
    end
end

Then, in tsp_decode!, we simply refer to the array according to the local thread ID:

function tsp_decode!(chromosome::Array{Float64}, instance::TSP_Instance,
                     rewrite::Bool)::Float64
    permutation = instance.permutation_per_thread[threadid()]
    #...
    #...
end

Note that to pre-allocate decoding structures inside the object holding the instance is not the most elegant and decoupled code we can write. However, to decouple the decoding data from the instance data requires that we pass another data object to the decoder. To do this explicitly, we may get an embroidered API. We could do it implicitly, by creating a Singleton object to hold the decoding data. However, this also reduces (a lot) the clarity and objectivity of the code. In C++ code, this is much easier accomplished by creating a Decoder object that can hold data members as much as methods. Therefore, when creating the Decoder object, we can pre-allocate all data structures we need.