# BROUTE: a benchmark suite for the implementation of standard vehicle routing
algorithms

Fabien Tricoire  Vienna University of Economics and Business, Welthandelsplatz
1, 1020 Wien, Austria
<EMAIL_ADDRESS>
###### Abstract

We introduce BROUTE, a benchmark suite for vehicle routing optimization
algorithms. We define a selection of algorithms traditionally used in vehicle
routing optimization. They capture essential features that are also relevant
in optimization algorithms for different application domains, like local
search move evaluation, memory allocation, dynamic programming, or insertion
and deletion from a list. Each algorithm is deterministic. We implement these
benchmark algorithms using a selection of programming languages and different
data structures. BROUTE is free, open-source, and can be used to inform early
decisions in projects that involve programming, such as which language to use.

##  1 Introduction

When it comes to implementing optimization algorithms and routing algorithms,
some important design decisions such as the choice of the programming language
must be made early. While most programming languages are Turing-complete, thus
allowing to run any deterministic algorithm, they vary in many aspects
including conciseness, ease of writing, ease of reading, ease of debugging and
runtime performance. Many of these aspects are subjective in nature, but
runtime performance can be measured objectively. We introduce a benchmark
suite for vehicle routing optimization algorithms called BROUTE, available at
https://github.com/fa-bien/broute, which focuses on CPU effort. CPU effort is
especially relevant when running the final version of a program, as it is a
direct measure of performance. However it also matters while developing and
debugging, as these tasks typically involve repeatedly running the program and
waiting for the runs to complete.

The goal of this suite is to inform early design decisions when programming
optimization algorithms, such as the choice of the programming language, or
how to represent a distance matrix. Readers may already know a few programming
languages and know both (i) which one they appreciate using the most and (ii)
which one offers the best performance, but they do not necessarily have a
quantified assessment of the trade-off between performance and convenience,
i.e. of the performance cost of convenience. Other readers may need to decide
which programming language their new student should learn and use for the
implementation of their research work. Finally, some readers may simply want
to know which language is the most efficient for which task. These decisions
can all be informed by using BROUTE.

###  1.1 Scope of this study

We want to compare the reasonably efficient implementation of various
algorithms used in combinatorial optimization. In particular, we look at five
algorithms that are commonly used in vehicle routing optimization. However
these algorithms capture features that are also common in other domains of
optimization, such as modifying subsequences, inserting and removing elements
in a sequence, or other operations on graphs.

Our goal is not to provide the most efficient implementation of each algorithm
for a set of languages, or even to use the most efficient algorithm to solve
the given problems, rather to use these algorithms as benchmarks, each
capturing certain features currently encountered in optimization. We put
ourselves in the situation of a programmer with limited experience having to
implement these algorithms; the performance of different programming languages
are compared in that context. for the time being we implement all benchmarks
in C++, Python, Java, JavaScript, Rust and Julia. We also assess the impact on
performance of some implementation decisions, e.g. using a flat matrix
representation.

###  1.2 What is measured

In all proposed benchmarks, the CPU time of the benchmarked algorithm is
measured. Every other task, such as reading input data or constructing an
initial solution for local search benchmarks, is not measured. However if an
auxilliary graph needs to be constructed every time the algorithm is run, then
we also measure the construction of said auxilliary graph. This is the case
for benchmarks _maxflow_ and _espprc_.

In general, compilation time is not measured. However, some implementations do
just-in-time (JIT) compilation, for example Julia, JavaScript, Pypy, Numba and
Java. In such cases, this JIT compilation time is included in the running time
of the algorithm. In general we measure _clock_ time, i.e. the CPU budget used
by the algorithm. In the case of JavaScript we measure _wall_ time as we are
not aware of any way to measure clock time in JavaScript.

In order to represent a tour in various heuristics, we use variable-size
vectors. Dynamic size containers are intuitive to use in algorithms for
vehicle routing, for instance to insert or delete a vertex from a tour. Such
containers can also be used to represent solutions in other domains, e.g.
scheduling or any application that requires the representation of sequences of
variable length. Inserting and deleting elements in a vector has linear worst-
case time complexity. However, for heuristics it can still be a good idea to
use vectors rather than, say, linked lists, since neighborhood exploration
typically requires to evaluate many moves then perform only one of these
moves. The complexity of performing a move is thus a secondary concern and
unlikely to be a performance bottleneck. Additionally, performing certain
moves can actually be more efficient using a vector than using a linked list,
e.g. inverting a subsequence in the context of 2-opt.

###  1.3 Limitations

Benchmarks are by definition limited in scope; any given benchmark measurement
is only valid within a certain environment. In our case such environmental
factors include the choice of compilers, compiler versions, operating system
and CPU family. Therefore all results are to be taken as indicators but not
absolute truths. The intent of this study is not to declare a language the
winner because it is consistently faster than others by a few percents, but
rather to assess what kind of performance losses can be expected by, for
example, using Python instead of C++. Another example is to assess what kind
of performance gains can be achieved by using a JIT compiler for Python.

Additionally, there are some implicit conditions for the good functioning of
interpreters and compilers. For instance Pypy is advertised to work better
when things are “kept simple”, which is generally good advice regardless of
the compiler. However, this typically requires discipline.

##  2 The benchmarks

We consider five different benchmark algorithms that perform tasks commonly
encountered in transport optimization. Together they capture key aspects of
computationally intensive tasks in heuristics and exact methods for transport
optimization, although they do not necessarily use the most efficient
algorithm to solve the problem they are tackling. Additionally, we believe
that they also capture such aspects for other fields of application of
operations research, such as e.g. scheduling. For the sake of simplicity, we
consider the symmetric travelling salesperson problem (TSP) as a base problem,
i.e. input data are in the form of a symmetric distance matrix while a
solution is a permutation. However some benchmark algorithms solve different
problems than the TSP based on these TSP data, as described below. We do not
provide abstract algorithms for each algorithms; instead, we refer to the
actual implementations, all available at https://github.com/fa-bien/broute. In
particular, the Python implementation is close to pseudo-code.

All five benchmarking algorithms described below are deterministic.
Additionally, the input provided to them is also deterministic, i.e. each
implementation receives the exact same inputs and performs the same operations
from these inputs. For that purpose, we generate various instance files, each
with a given size $n$, which represents the number of vertices in the graph,
and a number $p$ of permutations. A permutation is a solution to the TSP.
Vertex coordinates are generated randomly in $[0,100)$. Let $d_{ij}$ be the
Euclidean distance between $i$ and $j$, then the distance considered is
$c_{ij}=\lfloor 100d_{ij}\rfloor$. Since it is an integer number, it is
sufficient to use integer number representation to compute the cost of a TSP
solution; however certain benchmarks use floating-point number representation,
as explained below. Each instance file contains one $n\times n$ distance
matrix as well as $p$ randomly generated permutations, which are used as
starting _seeds_ for the algorithms (e.g. as starting solution for 2-opt).

Since the same instance files are given to each implementation, all
implementations perform the same operations and return the same result. In
order to control result integrity we use a mechanism similar to _checksum_.
The checksum calculation for applying a certain algorithm using a given
permutation as starting seed differs based on the algorithm and is explained
separately for each algorithm below. The checksum calculation for an instance
file is the sum of checksum values over all permutations in that instance
file.

Performances are measured per instance file, i.e. each time reported is for
running an algorithm $p$ times, using one distance matrix with the $p$
different seed permutations.

###  2.1 2-opt

The 2-opt heuristic was first described by Croes, (1958). One of the most
commonly used heuristics in vehicle routing, 2-opt improves a tour by
performing 2-exchanges. A 2-exchange consists in removing two edges from a
tour and reconnecting the tour with two other edges. It is equivalent to
inverting a sub-sequence of the tour, and can be performed in place using an
array or vector solution representation. If the distances are symmetric, then
each move can be evaluated in constant time.

The checksum for a given permutation seed is calculated as the number of
improvements found while applying first-improvement 2-opt until no improving
move exists, using this permutation as starting solution.

###  2.2 Or-opt

Or-opt is a heuristic that was first introduced by Or, (1976). It relies on
the exploration of a neighborhood that is a subset of all 3-exchanges. A
3-exchange consists in removing 3 edges from a tour and reconnecting the tour
in a different way. Of all these moves, Or-opt only considers those that shift
a sequence of 1, 2 or 3 vertices to a different position in the tour. Each
move can be evaluated in constant time.

The checksum for a given permutation seed is calculated as the number of
improvements found while applying first-improvement Or-opt until no improving
move exists, using this permutation as starting solution.

###  2.3 lns: Large neighbourhood search

Large neighborhood search (LNS) was originally introduced in conjunction with
constraint programming to solve the vehicle routing problem with time windows
(VRPTW) Shaw, (1998). Since it has been a very popular method for tackling a
wide variety of vehicle routing problems, see e.g. Pisinger and Ropke, (2010).

LNS iteratively (i) copies the _incumbent_ solution, (ii) destroys the copy by
removing some vertices from it, (iii) repairs the partial solution by re-
inserting previously removed customers in it and (iv) determines whether this
newly produced solution becomes the new incumbent solution or not. The method
relies, among other things, on randomness. This is an issue since we want each
benchmark to be deterministic. To remedy this issue, we design deterministic
schemes for both destroying and repairing stages. For each starting seed, we
apply a certain number of LNS iterations. At each iteration, the vertices at
even indices in the tour are removed. These removed vertices are stored in a
vector in the order in which they were removed, i.e. from smallest to largest
index. Then the cost of each possible insertion is calculated and the cheapest
insertion is performed. Because removed vertices are stored in a vector, the
algorithm’s behaviour in case of same-cost insertions is deterministic.

The checksum for a given permutation seed is calculated as the total cost of
insertions performed while applying 10 iterations of LNS using that
permutation as starting solution.

###  2.4 espprc and espprc-index: Dynamic programming for column generation

Many state-of-the-art exact methods for vehicle routing problems rely on
column generation (see e.g. Baldacci et al., (2012)). In such methods, a
significant amount of CPU effort is spent solving the _pricing subproblem_ ,
which differs depending on the problem at hand but is nonetheless often an
elementary shortest path problem with resource constraints (ESPPRC). While it
can be solved heuristically, in order to establish optimality it eventually
needs to be solved exactly as well. For this benchmark we implement a dynamic
programming algorithm similar to the one by Feillet et al., (2004) to solve
the ESPPRC, albeit without time windows. One difficulty in this problem stems
from the existence of negative-cost cycles in the graph.

In order to simulate an environment similar to the one encountered in column
generation for VRPs, we derive an ESPPRC instance from each TSP seed
permutation by using an auxilliary distance matrix $c^{\prime}$ from the
distance matrix $c$, as follows:

  1. 1.

For each vertex $j$, let $\pi_{j}$ be the cost of arc $(i,j)$, where $i$ is
the predecessor of $j$ in the seed permutation.

  2. 2.

For each arc $(i,j)$, set $c^{\prime}_{ij}$ to $c_{ij}-d_{j}$.

Additionally, we generate resource consumption as follows: we consider 6
resources numbered from 0 to 5, and any vertex $i$ consumes resource $r$ iff
the $r^{th}$ bit of $i$’s binary representation is 1 (considering that the
first bit is bit 0). That is, in C syntax, iff `i & (1 << r) > 0`. Capacity is
set to 1 for each resource. The algorithm then computes the elementary
shortest path with resource constraints on $c^{\prime}$ from 0 to 0 (1 to 1 in
Julia).

For this benchmark, floating-point numbers are used to represent cost, even
though the values are integer. This is to emulate the usual setting when
solving such problems in the context of column generation.

Each label needs (i) a reference to its predecessor, so that paths can be
reconstructed and (ii) a collection of references to its successors, in order
to recursively propagate dominance and later ignore labels that are successors
of an already proved-to-be-dominated label. Rust does not allow that, at least
not in the intuitive way where a label is encapsulated in a `struct`,
including references to predecessor and successors. This is because Rust is
focused on security, while this type of construct can lead to unsafe
operations if not used properly. This means that this benchmark cannot be
safely implemented with Rust, even though it is consistent and reliable in all
other implementations. In order to remediate this issue, we produce a second
implementation of the same algorithm, although in a less intuitive way: all
labels are stored in a collection, and the indices in this collection are used
in other collections to store all predecessors and successors. This second
implementation represents a different benchmark, which we call _espprc-index_
, whereas the base benchmark is called _espprc_. The version with indices is
more tedious to implement, but it may also be more efficient, so we also
implement it in languages other than Rust.

The checksum for a given permutation seed is the cost of the shortest
elementary path obtained when using that permutation seed to generate the
ESPPRC instance, truncated to its integer part.

###  2.5 maxflow: Maximum flow problem

Another popular framework for the exact solution of routing problems is
_branch and cut_. Branch and cut is typically applied to solve mixed-integer
linear programs; it is akin to branch and bound but only considers a subset of
all constraints explicitly. Other constraints are generated dynamically when
they are found to be violated while exploring the search tree. These
constraints are said to be _separated_. One early success story of branch and
cut is actually in application to the symmetric TSP (Padberg and Rinaldi,,
1987), by separating subtour elimination constraints. The separation procedure
involves finding a minimum capacity cut, which itself involves solving a
maximum flow problem. For this reason the maximum flow problem is important in
vehicle routing, and we dedicate one benchmark to solving it. The _maxflow_
benchmark implements the Edmonds-Karp algorithm Edmonds and Karp, (1972),
which is a specific version of the Ford-Fulkerson Ford and Fulkerson, (1956)
augmenting path method. The specificity resides in how the augmenting paths
are generated.

We derive a capacity graph $C$ from the distance matrix $c$ and a given TSP
solution as follows:

  1. 1.

For each vertex $j$, let $t_{j}$ be the cost of arc $(i,j)$, where $i$ is the
predecessor of $j$ in the TSP solution.

  2. 2.

For each arc $(i,j)$, set $C_{ij}$ to $\frac{c_{ij}}{1000}$ if $c_{ij}>t_{j}$,
otherwise set it to 0.

For any starting seed, the capacity graph is generated then the Edmonds-Karp
algorithm is applied with 0 as _source_ and every other node as _sink_. The
checksum for that permutation seed is the sum of all maximum flow values thus
obtained, truncated to its integer part.

The reason why capacity values are distance values divided by 1000 is that
otherwise, checksum values can exceed the largest possible integer
representation with 32 bits, which can be a problem with some languages. This
means that floating-point numbers are used, which is representative of the use
case of branch and cut for vehicle routing problems, where arc capacities
represent flow in the current solution to the routing problem.

##  3 Languages and implementations considered

We consider six different programming languages: C++, Python, Java, Julia,
Rust and Javascript. Any rule of selection is arbitrary in nature; however we
indicate below, for each of these languages, our perceived advantages of these
languages, which justify their selection. It is possible to add other
languages to this list in the future, and keep the benchmark collection alive
in an online form. For a given language several implementations are possible,
e.g. one with nested distance matrices and one with flat distance matrices. As
a convention, implementations are named in lower-case letters, for example
_python-flat-matrix_. We also use the simplest names for the best
implementations, i.e. _python_ uses a nested distance matrix because it
performs better than _python-flat-matrix_ , but _pypy_ uses a flat matrix
because it performs better than _pypy-nested-matrix_. Language versions and
options are summarised in Appendix A.

For each language, we may consider multiple _implementations_ , which differ
by which mechanisms or library they use. For instance with Python we can
compare an implementation using native lists against an implementation using
Numpy arrays. Implementations are detailed individually below for each
language. The code for all implementations can be found at
https://github.com/fa-bien/broute.

It is difficult to properly measure the popularity of a programming language.
There exist rankings, for example based on online searches TIOBE, (2020) or on
GitHub contributions GitHub, (2020). All these approaches have flaws. Basing
an index on online searches is subject to all the biases and exploits that
target online search. Measuring GitHub contributions is not clear, as there is
no trivial way to decide what should count among repositories, number of lines
written, number of commits, etc. Ultimately, language popularity for its own
sake does not need to be a crucial factor in the choice of a programming
language, especially for academic research. These aspects seem more important
to us:

  * •

The programs written in a given programming language should be easily run on
any computer. For that reason, we only look at languages that are free and
work under multiple operating systems. The ability to run under Linux is
especially important, since it is typically the operating system of high-
performance computing (HPC) clusters TOP500, (2021).

  * •

The language should be well documented, but also have online resources for
support. This in fact also depends on popularity.

###  3.1 C++

C++ is one of the most popular programming languages in the World, and has
been so for decades. Among its many features there is good runtime
performance. In fact the performance is good enough that whole operating
systems can be written in C++, for example Microsoft Windows. C++ has many
features and it is often possible to implement the same algorithm using
different subsets of the language; one criticism is that people and teams
evolve into using only their own subset of the language.

We develop two C++ implementations of every benchmark, using the C++98 and
C++14 standards, which we call _c++98_ and _c++14_. they both use a flat
distance matrix. With C++98 we use references to objects, while with C++14 we
use smart pointers such as `shared_ptr`. Comparing these two implementations
will allow us to evaluate the cost of smart pointers.

For each benchmark the reference CPU time is provided by the c++98
implementation. We also develop a variant of c++14 using a nested distance
matrix, and another using static arrays.

###  3.2 Python

Python is another of the most popular programming languages in the World,
especially for tasks of machine learning in recent years. The base Python
implementation is called _python_ and uses a nested distance matrix. There is
also an implementation using a flat distance matrix, called _python-flat-
matrix_ , and one using a nested distance matrix wrapped in function calls,
called _python-nested-matrix-function_.

One issue with Python is that much of the code is not compiled, rather
directly interpreted, which can lead to poor performance. One way to address
that issue is to rely on libraries that are implemented in C, however that is
not always possible. We consider alternative implementations in Python, all
bringing different solutions to this issue.

The simplest way to speed up the run time of a Python program is to run this
program with Pypy Pypy, (2021) instead of the standard CPython interpreter.
Pypy uses JIT compilation to speed things up. There is no need to modify the
Python program, which is a strong advantage compared to library-based
solutions or annotation-based solutions. Pypy’s website claims that it is on
average 4.2 times faster than CPython Pypy, (2021), however this is just an
average and actual numbers may vary a lot depending on the benchmark. We use
both the python and python-flat-matrix implementations with Pypy, resulting in
_pypy-nested-matrix_ and _pypy_ implementations, respectively.

We also implement the benchmarks using Numpy arrays instead of standard Python
lists. Numpy is one of the most popular external libraries for Python; it is
coded in C for runtime performance Harris et al., (2020). Numpy performs
especially well when using _array operations_ , for instance matrix
multiplication. However the benchmarks that we consider do not take advantage
of Numpy functions, only of its array data structure. Since Numpy provides
multi-dimensional arrays, we use a two-dimensional array to represent the
distancematrix. Additionally, we consider a flat matrix in implementation
_numpy-flat-matrix_.

Numpy itself does not provide any kind of JIT compilation, but another Python
package, Numba, performs JIT compilation when using Numpy arrays. Numba
requires code annotation, therefore we annotate the Numpy implementations to
take advantage of them with Numba.

###  3.3 Java

Java is a very popular programming language of the past two decades. Java code
is compiled to _bytecode_ , which is then run by a virtual machine (VM). The
bytecode compiled on one computer can be run by a VM on another computer.
Nonetheless, JIT compilation to native code happens dynamically at run time.
The syntax is based on C/C++ syntax. The base Java implementation, _java_ ,
uses a flat distance matrix and a `ArrayList<Integer>` object to represent a
tour. We also implement a version with nested distance matrix, _java-nested-
matrix_ , as well as a version with static arrays to represent a tour, _java-
static-arrays_.

###  3.4 Julia

Julia is a relatively recent programming language. It is a general-purpose
language but it is also particularly aimed at computationally intensive tasks
such as numerical analysis Bezanson et al., (2017). Version 1.0 was released
in 2018 and the language is well-documented and has a strong community.

Compilation happens in a JIT fashion. Like Numpy but unlike all other
programming languages considered, Julia supports multi-dimensional arrays,
therefore a native two-dimensional array can be used to represent a distance
matrix. We also develop an implementation using flat matrix representation in
Julia, called _julia-flat-matrix_ ; we will determine which version is better
through experiment.

One valuable feature of Julia is the fact that projects have an environment,
including package dependencies with specific versions. The language makes it
easy to _instantiate_ the project environment, i.e. recreate the same
environment in which the project was developed, with the same packages in the
same versions. This provides safety from backward-compatibility-breaking
library updates.

One notable issue with Julia is the inability to generate a compiled binary
executable file. In general, in order to run a Julia program, one must install
the Julia ecosystem. There are active projects to remedy this, but we are not
aware of any easy one-command solution.

###  3.5 Rust

Rust is also relatively recent, with version 1.0 released in 2015. It is a
general-purpose language with an emphasis on memory safety and on performance.
Notable software using Rust include Firefox, Dropbox and Cloudflare The Rust
Team, (2021).

Programs in Rust are pre-compiled, like in C++ for instance. Like Julia, the
language is well documented and has a strong community. Also like Julia,
projects have an environment, and it is easy to recreate the environment for a
given project, including dependencies and their specific version, which again
provides safety from backward-compatibility-breaking library updates. The
_rust_ implementation uses a flat matrix representation.

###  3.6 JavaScript

JavaScript is likely one of the most popular programming languages in the
World, as it is used on countless websites as well as mobile applications
based on web technology. JavaScript is not expected to perform as well as,
say, C++. However there are several reasons why it can be appealing to
implement routing algorithms in JavaScript:

  * •

Any program written in JavaScript can run on virtually any computer with a
modern web browser, including smartphones.

  * •

Integrating JavaScript with a web interface is especially easy.

  * •

JavaScript engines have received considerable attention from major companies
and been the subject of fierce competition in relation with web browsers
engines. As a result they have seen vast performance improvements over the
years and the trend is likely to continue.

We use Node.js, which allows to run JavaScript programs from the command line.
This means that there is in fact not even the need for a web browser. Node.js
currently uses Google Chrome v8’s JavaScript engine nodejs.org, (2021).

The _javascript_ implementation uses a flat matrix representation but there is
also an implementation using nested matrix representation, called _javascript-
nested-matrix_.

###  3.7 Considerations on how to implement a distance matrix

Perhaps the most straightforward way to implement a distance matrix in a
number of languages is to use nested arrays, i.e. each element of the main
array is an array representing a row of the distance matrix. This is typically
done in C++, Java, Python (although the structure is officially called a
_list_ and not an array), Javascript. C and C++ use pointers to achieve that
effect. All these languages use the same C syntax for looking up values in the
distance matrix: the distance between vertices $i$ and $j$ using distance
matrix $d$ is written d[i][j].

Another easy way to implement a distance matrix is to use what we call from
here on a _flat_ distance matrix representation, which is a single-dimensional
array containing all distance values. Assuming indices start at 0, the
distance from $i$ to $j$ in array $d$ can be coded as d[i*n+j], where $n$ is
the total number of vertices. This representation guarantees that the whole
distance matrix is stored in contiguous memory. Additionally, looking up in
arrays also takes time, so one lookup is better than two. The drawback is that
we have to pay a multiplication and an addition for each lookup, and that
lookups have to be wrapped in a function call for the sake of readability.
Such function calls can usually be inlined, i.e. the function content is
substituted to the function call, so there should be no runtime penalty from
using a function call. In general we expect a speedup from using a flat
representation. We implemented some of the benchmarks with both flat and
nested matrix representation in order to determine if there are significant
performance differences.

There are a few cases where the above considerations on using a flat
representation do not apply: Python does not allow inlining, while Numpy and
Julia provide multi-dimensional arrays. Nonetheless we also wrote a flat
matrix version of these implementations, in order to elicit any difference in
performance.

###  3.8 Other general cross-language considerations

Here we discuss considerations that are subjective in nature, as they relate
to individual experience implementing all benchmarks in all languages. The
author’s experience with all considered programming languages is
heterogeneous, which can cause bias in the perception of the ease of
implementation in a given language. To add to the subjectivity, each benchmark
is first implemented in Python then in other languages, with Rust typically
being the last one to be implemented. Therefore the languages implemented last
represent supposedly easier tasks, due to increased familiarity with the
algorithm. For all reasons mentioned above, measuring the time spent
implementing each benchmark in each language would be meaningless, therefore
we limit ourselves to subjective opinion here. Keeping this in mind, we see a
subset of perceived _easy_ languages, as in “easy to program in and get things
running with the expected outcome”, which consists of Python, Julia and
JavaScript. These easy languages typically allow to program the same
algorithms faster than the other ones with the same results, performance
considerations aside. On the other hand, we found the process of getting
benchmarks to work with Rust to be much more tedious than with any other
language, although this is certainly due to lack of experience. Nonetheless,
the JavaScript implementation was typically easier than the C++ or Java ones
despite having more experience with C++ and Java, so experience does not
explain everything. Ultimately, there is little doubt that any experienced
programmer in any of these languages would have no trouble implementing any of
these benchmarks. Nonetheless, this programmer with experience in mostly C++
and Python implemented every benchmark faster in Python, JavaScript or Julia
than in C++.

We also want to underline that no matter the language, it is possible to write
both efficient and inefficient code. Familiarity and experience are
beneficial; beginners in a language are more likely to write inefficient code.
This being said, not every language is equal when it comes to documenting how
to write efficient code. We note that the official Julia documentation has a
whole section dedicated to that topic Julia Documentation, (2021). There also
exists an automated tool to detect inefficient code patterns, based on the
official documentation JunoLab, (2020).

##  4 Implementation-specific notes

###  4.1 Python

CPython, Python’s official interpreter, does not allow to inline functions.
Since using a flat matrix requires to wrap matrix value lookup in a function,
using a flat distance matrix in Python is not beneficial. It would of course
be possible to inline by hand, i.e. formulate the correct 1-dimensional index
in the matrix for every lookup, but that would be impractical and defeat the
purposes and philosophy stated in Section 1. Running Python code that uses a
flat representation with a JIT compiler (e.g. Pypy, Numba) introduces
automatic inlining and is likely to remedy this.

On the ESPPRC benchmark, a label needs to store information about which
customers have already been visited in the partial path it represents. In
Python, preliminary testing reveals that using a set of integers for that
purpose is twice as fast as using a vector of booleans. Similar preliminary
testing in C++ reveals the opposite: using an array of booleans is faster.
This is the behaviour we would expect in general, the observation on Python is
the surprising one here. There are two reasons for this expectation: (i)
copying an array, represented in contiguous memory, can in general be
performed faster than copying a set structure and (ii) lookup in an array is
performed in $\mathcal{O}(1)$ while lookup in a set is performed in
$\mathcal{O}(log(n))$.

###  4.2 Julia

Julia has native multi-dimensional arrays, which may be more efficient than a
flat representation. We will determine which version to use based on
experiment. Since Julia’s arrays start with index 1, the code for flat
representation lookup of the distance between vertices $i$ and $j$ in matrix
$d$ is d[i*n-n+j]. It is worth noting that it involves one subtraction on top
of the addition and multiplication used for 0-indexed languages.

###  4.3 JavaScript

All benchmarks are repetitive in nature, as they successively perform similar
operations with different input data. In JavaScript, when our implementation
loads 40 instances and performs the same benchmarking operations on all 40
them, we observe a significant performance hit after a few instances (2-5
times slower). In order to remedy this, we re-start the program individually
for each instance. Advanced knowledge of this specific JavaScript engine might
allow to remove that performance hit, but this is clearly outside the scope of
this study. However it is worth noting that this can be a concern in general.

###  4.4 Numba

Numba is not yet feature-complete for Numpy. A direct consequence is that some
benchmarks cannot be implemented for Numba. For example `numpy.insert()` is
not implemented, but this function is needed for the $lns$ benchmark,
therefore the $lns$ benchmark cannot be implemented with Numba. Additionally,
we were not able to implement the ESPPRC benchmark in a satisfying manner. In
our experience, Numba error messages are sometimes unrelated to the issue
causing them, or too obscure to make sense, keeping in mind the premises of
this paper. It would likely be possible to implement a dynamic programming
algorithm for the espprc benchmark using Numba’s JIT compilation features,
however our conclusion is that it would not be possible to do so without
advanced knowledge of Numba.

Still, Numba is relevant enough to be benchmarked and its current performance
is an indication of what to expect in years to come, as the project matures.

###  4.5 Java

The ecosystem of Java implementations and versions can be hard to navigate.
There exist multiple compilers, multiple virtual machines, and there is a new
Java version every 6 months, sometimes breaking backwards compatibility. The
official reference implementation is OpenJDK since version 7 Oracle
Corporation, (2011). We use OpenJDK version 11, which is the current long-term
support version as of writing this article.

###  4.6 Rust

As mentioned in Section 2.4, Rust does not allow a given `struct` to have
multiple references to the same `struct` type. For this reason, a different
implementation of the espprc benchmark is needed. In general Rust puts
restrictions on what constitutes valid code, for the sake of safety. This is a
valuable feature in systems programming, which is the main application that
Rust was designed for. This also means that writing code that compiles and
runs is actually a more tedious task, but that once the code compiles, there
is also a higher chance that it will run according to expectations, like with
the Ada programming language for instance.

In the context of this study, these benefits did not appear relevant however:
each benchmark was programmed in all other languages as well, providing the
same results as with the Rust implementation, while no evidence of memory leak
was observed. As mentioned above, The Rust implementation was typically way
more time-consuming than other implementations, although one reason for that
is lack of experience. Nonetheless, in the context of this study, the
generally assumed benefits of using Rust are not apparent, with the notable
exception of good runtime performance.

##  5 Experiments

In this section, we explore some possibilities offered by BROUTE and perform
topic-specific comparisons. Unless indicated otherwise, all programs are run
concurrently on an desktop computer with an Intel Xeon E5-2650 v3, 20 cores, 3
GHz, running Linux. For every benchmark and every instance, CPU times are
reported as a ratio of a _reference time_ for that benchmark and instance. In
practice, for each benchmark and instance, the reference time is the time
required by the c++14 implementation for that benchmark and that instance.
Therefore the c++14 implementation always has a CPU ratio of 1.

A summary of all versions used for compilers and interpreters is provided in
Appendix A.

###  5.1 Benchmark data

Each instance file contains one $n\times n$ distance matrix as well as $p$
randomly generated permutations, which are used as starting _seeds_ for the
algorithms (e.g. as starting solution for 2-opt). For $n$ taking each value of
$\\{20,40,60\\}$ and fixing $p=1000$ we generate 40 instance files, for a
total of 120 instance files. Each instance file contains 1000 permutations,
therefore involves running each benchmark 1000 times. As mentioned in Section
2, distances are Euclidean distances multiplied by 100 and truncated to their
integer part. In the following, each data point represents the total CPU time
for running one benchmark implementation for each of the 1000 permutations of
one given instance file.

These data points are represented in two types of charts: box plots and
performance profiles. Box plots contain different boxes for different values
of $n$. Performance profiles aggregate all data points for a given
implementation into one profile, i.e. one curve. A performance profile is a
chart that indicates, for a given level of performance, the ratio of data
points produced with that performance level or a better one. In our case, the
performance level is the CPU effort required to produce a data point, as a
ratio of the CPU time taken by the fastest implementation over all
implementations compared in the same chart, for the same input file. For more
information, see Dolan and Moré, (2002). Box plots are produced with the
_ggplot2_ package Wickham, (2016) in R R Core Team, (2021). Performance
profiles are produced using a homemade script in Python, relying on matplotlib
Hunter, (2007).

###  5.2 A comparison of Python implementations

We first compare different implementations based on Python:

  * •

Implementation ’python’ uses CPython, the official Python interpreter

  * •

Implementation ’numpy’ uses Numpy, a Python library for numerical operations
that relies on code implemented in C. We also call it from CPython.

  * •

Implementation ’numba’ used Numba, a Python library that performs JIT
compilation. Numba is advertised to work well in conjunction with Numpy, i.e.
the functions from Numpy are compiled by Numba. Numba is under development. It
is also called from CPython.

  * •

Implementation ’pypy’ uses Pypy, a different Python interpreter that also
performs JIT compilation. The code is exactly the same as with implementation
’python’.

We first compare the performance of these four implementations on benchmarks
2-opt and Or-opt. The results for this comparison are summarised in Figure 1.

(a) 2-opt benchmark

(b) Or-opt benchmark

Figure 1: Python implementations: box plots and performance profiles

The box plots indicate that the performance of the numpy implementation is
always more than 400 times slower than the c++14 implementation. The
performance profiles show that it is also up to 600 times slower than the
fastest Python implementation in some cases. Overall, it is clear that the
Numpy implementation has the worst performance. This is easily explained:
every time a function from Numpy is called, it is wrapped in order to call the
C code for that function. At the same time, none of the benchmarks actually
use functions which would make this beneficial, i.e. functions that would
perform a lot of calculations written in C then return the result of this
calculation (e.g. matrix multiplication). Such operations are simply not
needed in these benchmarks. Therefore the cost of calling C code is paid
often, but the benefits are never reaped. Additionally, it is clear that the
base CPython implementation is also very slow, although not as much as the
Numpy one. However, the Numpy structures combined with JIT compilation from
Numba look promising.

We now compare Pypy and Numba only, on all benchmarks except for lns and
espprc, since they are not implemented with Numba as explained above. These
results are presented in Figure 2.

(a) 2-opt benchmark

(b) Or-opt benchmark

(c) maxflow benchmark

Figure 2: Pypy vs Numba: box plots and performance profiles

Pypy is up to almost twice as fast than Numba on the 2-opt benchmark, while
Numba is more than three times as fast in a majority of cases for the Or-opt
benchmark and up to twice as fast for the maxflow benchmark. We also note that
for the 2-opt benchmark, the performance gap reduces when instance size
increases, keeping in mind that the 2-opt benchmark is the fastest of the lot.
This may simply mean that Numba has a higher initial cost for JIT compilation.
Overall it seems that Numba is faster, except for very quick computations.
However we still see two reasons to use Pypy: it does not require any
modification to the Python code, and it is feature-complete, thus allowing to
implement all benchmarks. In fact the code for the CPython and Pypy
implementations of each benchmark is fully identical.

###  5.3 Impact of using a flat matrix

We now look at the impact of using a flat distance matrix representation,
across multiple languages. The expectation is that it should generally be at
least as fast, as long as function calls are inlined; however, in the case of
languages that have multi-dimensional arrays, it is harder to predict. We
first look at such a comparison for the c++14 implementation in Figure 3,
considering benchmarks 2-opt and Or-Opt. Benchmarks maxflow and espprc use
their own matrix structure, therefore we do not include them in this
comparison.

(a) 2-opt benchmark

(b) Or-opt benchmark

(c) lns benchmark

Figure 3: C++14: comparison of matrix representations

Although not depicted here, similar observations can be made with most
languages, with the notable exception of base Python (i.e. using CPython
without Numba). Since CPython does not do any kind of function inlining, the
performance cost of calling functions for each distance lookup outweighs the
performance benefit of using a flat, contiguous matrix representation.

We now compare a flat matrix representation in Julia against the “natural”
representation using Julia’s native multi-dimensional arrays. This comparison
is presented in Figure 4.

(a) 2-opt benchmark

(b) Or-opt benchmark

(c) lns benchmark

Figure 4: Julia: comparison of matrix representations

While it is not overwhelming, there can still be a benefit in using a flat
matrix representation.

Overall, it appears that it is never a significantly bad idea to use a flat
matrix representation, while it can bring benefits. Therefore, from this point
on, flat matrix representation is done by default.

###  5.4 Impact of using static arrays

In some cases, variable-size arrays are not necessary. This is the case for
example when implementing a 2-opt procedure for a TSP. In these cases, fixed-
size arrays are enough and their performance have to be at least as good as
the performance of variable-size arrays (like C++’s vectors). We run this
comparison for benchmarks 2-opt and Or-opt, which can be both implemented
easily using the fixed-size native arrays of C++, which are C native arrays,
and which work in a similar fashion in Java. We run this comparison for both
C++14 and Java, the results are summarized in Figure 5.

(a) C++14, 2-opt benchmark

(b) C++14, Or-opt benchmark

(c) Java, 2-opt benchmark

(d) Java, Or-opt benchmark

Figure 5: Performance of static arrays versus variable-size vectors.

Using static arrays in C++ can be up to almost twice as fast, while in Java it
can be up to five times as fast. One possible explanation why the gap is
larger in Java lies in the nature of the Java containers used: we use
`ArrayList` objects, which take as generic parameter the class of the objects
stored in the collection, e.g. `ArrayList<Integer>`. This parameter can only
be a class, not a native type, which matters in Java. This means that it is
not possible to have a collection of `int` variables, only a collection of
`Integer` objects, each wrapping an `int`.

###  5.5 A comparison of C++ implementations

We investigate the difference in performance between using C++14’s smart
pointers and older traditional C++ coding with references and/or C pointers.
For that purpose, we compare the C++14 and C++98 implementations of each
benchmark. Figure 6 summarizes this comparison.

(a) 2-opt benchmark

(b) Or-opt benchmark

(c) lns benchmark

(d) espprc benchmark

(e) maxflow benchmark

Figure 6: Performance of C++14 versus C++98.

The C++98 implementation is faster in all five benchmarks, although not always
by much. We conclude that if something can be implemented without using smart
pointers and without extra effort, then it is a good idea to do so.

###  5.6 General cross-language comparison

We now compare the performance of various programming languages when
performing the same tasks. We compare the following languages: C++, Java,
JavaScript, Julia, Python and Rust. For that purpose, we use what seems to be
the best choice for each language based on the above experiments, meaning that
we use the C++98 implementation for C++, and the Pypy interpreter when running
a Python program. There is a natural separation between compiled “fast”
languages (C++, Java, Julia, Rust) on one side and “interpreted” languages
(Python, JavaScript) on the other side. No language is exclusively
interpreted, as Pypy and Node use JIT compiling; this being said, there is a
performance gap between fast and interpreted languages, and in order to
produce more readable charts we compare them separately.

####  5.6.1 Cross-language comparison: “fast” languages

We compare the performance of fast languages in Figure 7.

(a) 2-opt benchmark

(b) Or-opt benchmark

(c) lns benchmark

(d) espprc benchmark

(e) maxflow benchmark

(f) espprc-index benchmark

Figure 7: Performance of “fast” languages.

There is no clear winner: each language is the fastest for at least one of the
benchmarks. There are however some relatively large factors in the execution
speed in some cases, up to a factor 8 in the case of Java for the Or-opt
benchmark. However in many cases the factor between two languages is never
more than three. This is reassuring for the researchers and students who want
to use something else than C++ and are worried about the performance hit.

####  5.6.2 Cross-language comparison: “interpreted” languages

We compare the performance of interpreted languages in Figure 8.

(a) 2-opt benchmark

(b) Or-opt benchmark

(c) lns benchmark

(d) espprc benchmark

(e) maxflow benchmark

Figure 8: Performance of “interpreted” languages.

Again there is no clear winner: both Python and JavaScript offer the best
performance for some of the benchmarks. More interestingly, these
implementations are never slower than ten times the effort required by the
baseline C++14 implementation (with the exception of one outlier for the
Python implementation of the Or-opt benchmark). In many cases, they are less
than five times slower even. This means that depending on the use case,
implementing routing optimization algorithms in Python or JavaScript can be a
viable (and convenient) option.

###  5.7 Cross-platform comparison

Many findings presented above can be reproduced on different platforms,
however in some cases different platforms produce different outcomes. We now
illustrate this with an example, looking at the comparison of Pypy and Numba
on the 2-opt benchmark. We run the 2-opt benchmark on the previously mentioned
Xeon CPU as well as on a Raspberry Pi 4 board with an ARMv7 CPU, 4 cores, 1.5
GHz. Pypy versions differ slightly: 7.0.0 on the ARMv7 CPU versus 7.3.3 on the
Xeon CPU. The Numba version is 0.42.0 on both platforms. Figure 9 summarises
and compares these runs. Ratios are computed independently for each platform,
e.g. the CPU time on ARMv7 is presented as a ratio of the CPU time taken by
the C++14 implementation on ARMv7.

(a) Computing cluster: Intel Xeon E5-2650 v3, 20 cores, 3 GHz.

(b) Raspberry Pi 4: ARMv7, 4 cores, 1.5 GHz.

Figure 9: Pypy versus Numba on 2-opt benchmark: different platforms.

The trends are opposite on different platforms. This may be due to several
reasons, for instance certain JIT compilers working better with certain CPUs,
or version differences for Pypy and/or Numba. Ultimately, we emphasize that it
would be unwise to generalise the observations made previously to different
settings. The goal of this article is not to provide final conclusions on what
tool to use, rather to introduce the BROUTE benchmark suite and encourage
using it when making early design decisions on implementation, using the
appropriate experimental platform to draw the appropriate conclusions.

##  6 Discussion

The question of what programming language to use when implementing
optimization algorithms is one that, by essence, does not have a final answer.
This being said, there is value in knowing how much performance is gained or
lost when using a certain language on a certain architecture, or how much can
be gained by using different data structures. The benchmark suite introduced
in this paper is free software and can be used for that purpose. We encourage
students and researchers in transport optimization to use it to support their
implementation decisions. Further experimental analyses, similar in flavor to
the ones conducted here, can also be good practice. For instance, comparing
the runtime performance of the same C++ code when compiled with different
compilers may provide a speedup for very little effort.

Additionally, some insight can already be taken straight out of the
experiments conducted above. Using languages like Julia or Rust does not incur
a significant runtime penalty when compared to more established languages like
C++, and can actually be faster sometimes. Interpreted languages like Python
or JavaScript nowadays have reasonably good runtime performance, usually less
than five times slower than pre-compiled languages, under the right conditions
(typically, using a JIT compiler). This means that the flexibility and ease of
use of these languages can be worth the performance hit, depending on the
application. Additionally, using JavaScript opens perspectives, as the
integration of optimization algorithms in web applications is easier than it
ever was, and these algorithms can now be run efficiently client-side.

This is only the initial version of BROUTE; contributions are encouraged,
concerning both languages and algorithms. For example a Matlab or C#
implementation of all benchmark algorithms would be welcome as long as they
can be run for free under the main operating systems. Implementing a new
benchmark in all currently supported languages would also be welcome, or even
re-implementing current benchmark algorithms in currently supported languages
but with different data structures.

## Acknowledgements

The author is grateful to Sebastian Leitner for his valuable support in
relation to the Rust implementation, and to Belma Turan for her valuable
comments.

## Declarations

### Conflicts of interest:

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

### Availability of data and material:

All data sets are licensed under the GNU General Public License version 3 and
available at https://github.com/fa-bien/broute.

### Code availability:

All code used in this work is licensed under the GNU General Public License
version 3 and available at https://github.com/fa-bien/broute.

## References

  * Baldacci et al., (2012) Baldacci, R., Mingozzi, A., and Roberti, R. (2012).  Recent exact algorithms for solving the vehicle routing problem under capacity and time window constraints.  European Journal of Operational Research, 218(1):1 – 6.
  * Bezanson et al., (2017) Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. (2017).  Julia: A fresh approach to numerical computing.  SIAM Review, 59(1):65–98.
  * Croes, (1958) Croes, G. (1958).  A method for solving traveling-salesman problems.  Operations Research, 6(6):791–908.
  * Dolan and Moré, (2002) Dolan, E. D. and Moré, J. J. (2002).  Benchmarking optimization software with performance profiles.  Mathematical Programming, 91(2):201–213.
  * Edmonds and Karp, (1972) Edmonds, J. and Karp, R. M. (1972).  Theoretical improvements in algorithmic efficiency for network flow problems.  Journal of the ACM, 19(2):248–264.
  * Feillet et al., (2004) Feillet, D., Dejax, P., Gendreau, M., and Gueguen, C. (2004).  An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems.  Networks, 44(3):216–229.
  * Ford and Fulkerson, (1956) Ford, L. R. and Fulkerson, D. R. (1956).  Maximal flow through a network.  Canadian Journal of Mathematics, 8:399–404.  Publisher: Cambridge University Press.
  * GitHub, (2020) GitHub (2020).  The 2020 state of the octoverse.  https://octoverse.github.com/.  [Online; accessed 17-Feb-2021].
  * Harris et al., (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del R’ıo, J. F., Wiebe, M., Peterson, P., G’erard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E. (2020).  Array programming with NumPy.  Nature, 585(7825):357–362.
  * Hunter, (2007) Hunter, J. D. (2007).  Matplotlib: A 2d graphics environment.  Computing in Science & Engineering, 9(3):90–95.
  * Julia Documentation, (2021) Julia Documentation (2021).  Performance tips.  https://docs.julialang.org/en/v1/manual/performance-tips/.  [Online; accessed 11-Feb-2021].
  * JunoLab, (2020) JunoLab (2020).  Traceur.  https://github.com/JunoLab/Traceur.jl.  [Online; accessed 11-Feb-2021].
  * nodejs.org, (2021) nodejs.org (2021).  Node.js.  https://nodejs.org.
  * Or, (1976) Or, I. (1976).  Traveling salesman-type combinatorial problems and their relation to the logistics of regional blood banking.  PhD thesis, Department of Industrial Engineering and Management Sciences, Northwestern University, Evanston, IL.
  * Oracle Corporation, (2011) Oracle Corporation (2011).  Moving to openjdk as the official java se 7 reference implementation.  https://blogs.oracle.com/java/moving-to-openjdk-as-the-official-java-se-7-reference-implementation.  [Online; accessed 10-Feb-2021].
  * Padberg and Rinaldi, (1987) Padberg, M. and Rinaldi, G. (1987).  Optimization of a 532-city symmetric traveling salesman problem by branch and cut.  Operations Research Letters, 6(1):1–7.
  * Pisinger and Ropke, (2010) Pisinger, D. and Ropke, S. (2010).  Large Neighborhood Search, pages 399–419.  Kluwer’s International Series.
  * Pypy, (2021) Pypy (2021).  Pypy.  https://pypy.org.  [Online; accessed 17-Feb-2021].
  * R Core Team, (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Shaw, (1998) Shaw, P. (1998).  Using constraint programming and local search methods to solve vehicle routing problems.  In Principles and Practice of Constraint Programming - CP98, pages 417–431. Springer-Verlag.
  * The Rust Team, (2021) The Rust Team (2021).  Rust.  https://rust-lang.org.  [Online; accessed 18-Feb-2021].
  * TIOBE, (2020) TIOBE (2020).  Tiobe programming community index definition.  https://www.tiobe.com/tiobe-index/programming-languages-definition/.  [Online; accessed 17-Feb-2021].
  * TOP500, (2021) TOP500 (2021).  Top500.  https://top500.org/.  [Online; accessed 17-Feb-2021].
  * Wickham, (2016) Wickham, H. (2016).  ggplot2: Elegant Graphics for Data Analysis.  Springer-Verlag New York.

##  Appendix A Summary of languages and versions

Language/tool | language or compiler version | options
---|---|---
C++14 | g++ 8.3.0 | -Wall -ansi -pedantic -O3 -std=c++14
C++98 | g++ 8.3.0 | -Wall -ansi -pedantic -O3 -std=c++98
CPython | 3.7.3 |
Pypy | 7.3.3-beta0 |
Numpy | 1.16.2 |
Numba | 0.42.0 |
Java | OpenJDK version 11 |
Julia | 1.6.0 | –check-bounds=no –inline=yes -O3 -t 1
Rust | 1.51 | cargo build –release
JavaScript | Node v14.15.4 |
Table 1: Overview of languages and versions.