ACO Implementation for Sequence Alignment
with Genetic Algorithms
Aaron Lee
lee42@fas.harvard.edu

Livia King
lking@fas.harvard.edu

December 1, 2002
Abstract
In this paper, we implement Ant Colony Optimization (ACO) for sequence alignment. ACO is a meta-heuristic recently
developed for nearest neighbor approximations in large, NP-hard search spaces. Here we use a genetic algorithm approach
to evolve the best parameters for an ACO designed to align two sequences. We then used the best parameters found to
interpolate approximate optimal parameters for a given string length within a range. The basis of our comparison is the
alignment given by the Needleman-Wunsch algorithm. We found that ACO can indeed be applied to sequence alignment.
While it is computationally expensive compared to other equivalent algorithms, it is a promising algorithm that can be readily
applied to a variety of other biological problems.

Motivation

found, the ants would start to wander back along the path
that they had traversed to the ant colony. Because the ant
who came back first with food would double the pheromone
on the trail, the next ant to leave the ant colony in search of
food who have a higher chance to take that path because ants
will tend to follow trails that have the highest pheromone
concentration [1].
In nature, it was found that this often led to suboptimal
solutions. The decay of the pheromone was occasionally not
fast enough to cause the ants to seek other solutions; hence,
there was a slight propensity to follow the path of the ant
that returned back to the ant colony first. In other words,
there was too much of a weight on the ant that reached the
first food source, because the first food source found was not
necessarily the closest food source globally.
In computers, this limitation was easily avoided by adding
a pheromone evaporation eﬀect where the ants’ pheromones
would decrease in concentration as time progressed, preventing the early convergence onto local minimums [5].

Ant colony optimization (ACO) is a novel algorithm recently developed for nearest neighbor approximations in
large, NP-hard search spaces [1]. The algorithm itself is more
accurately described as a meta-heuristic that must be applied
diﬀerently to diﬀerent problems. ACOs have been applied
successfully to both the Traveling Salesman Problem (TSP)
[2] and the Quadratic Assignment Problem (QAP) [3].
More recently, better methods of attaining optimal values for ACOs have been developed. In 1998, Bottee et al.
used genetic algorithms on top of ACOs to achieve locally
optimal values for ACOs very quickly [4].
As far as we have been able to tell, the ACO algorithm has
not yet been applied to biological problems. Most biological problems today face the same diﬃculty as other NP-hard
computational problems. The ACO algorithm seems particularly well-suited for these kinds of applications.
We propose using ant colony optimization with genetic
algorithms and applying the resulting ACO to two sequence
alignment.

Critical Review

Introduction

Dorigo, Caro, and Gambardella outlined in their 1999 paper
the basic characteristics of the ant algorithm and reported
on some of the uses for which it has been tested. They
list the ant algorithm as deriving many of its characteristics
from real ant colonies. For example, artificial ants, like real
ants, are part of a colony of cooperating individuals, and use
pheromones to facilitate indirect communication between

Ant Colony Optimization is a meta-heuristic that was developed after studying the behavior of ants and learning how
ants were able to solve constraint optimization problems in
nature. It was found that ants laid down pheromone trails
as they traversed the ground. As soon as a food source was
1

themselves and other ants. Although each of the individuals in a population is simple and not capable of deriving a
good solution to a problem, the population as a whole can
come up with a good solution, concept that is now known
as swarm intelligence. Both real ants and artificial ants make
decisions using local information.
Unlike real ants, however, artificial ants have some abilities that allow them to derive higher quality solutions than
their natural counterparts. Artificial ants can keep a memory of past actions. Their methods of pheromone depositing are also diﬀerent. They may deposit varying amounts of
pheromones based on the quality of the solution found and
might deposit pheromones after the generation of a complete solution rather than step by step. Other abilities like
the ability to look ahead several steps, or to backtrack, can
also be added to artificial ants.
The first problem ant algorithms were applied to was the
traveling salesman problem (TSP), a commonly researched
NP hard problem in which a theoretical salesman has to
find the shortest path connecting n cities. Three algorithms based on the ACO meta-heuristic were defined, and
the ”ant-cycle” algorithm (later called Ant System, or AS),
a variation in which the ants deposited pheromone after
they finished constructing their solution was found to be
the most eﬀective. The ant algorithm performed as well as
several other general purpose heuristics on relatively small
problems (involving 30 to 75 cities). However, the ant algorithm was not able to come up with the best known solution when it was applied to larger problems, although it was
able to find good solutions quickly. Later variations on the
AS algorithm, most notably ant colony systems (ACS), were
more successful and beat many other heuristics in terms of
solution quality and CPU time when run on standard TSP
problems.
ACO was next applied to the Quadratic Assignment
problem, which involves assigning n facilities to n locations
so as to minimize cost. In this case, the ACO algorithms
tested are the best available heuristic. ACO has also been
shown to be the best available heuristic by far in the sequential ordering problem. Other problems to which ACO has
been applied include Job shop scheduling, vehicle routing,
graph coloring, shortest common supersequence, connection oriented network routing, and connection-less network
routing.
From the results listed by Dorigo et. al, it can be seen that
ACO is a valid approach that has been successful for many
problems. However, there are some limitations to ACO’s
potential. First, because of the ant algorithm’s random nature, the success of ant algorithms very problem dependent
and it is hard to predict what variations on the algorithm
will be successful without actually implementing the algorithm and testing it. There are also limitations as to which
problems can be solved with ACO. The problem must lend

itself to being represented by ants traveling across a virtual
space. Each state available to the ant must also not have
many neighbor states, as too many neighbor states will minimize the chance of two ants appearing in the same state,
and render the pheromone communication worthless.

Design
ACO
We closely modeled our approach after Bottee and
Bonabeau’s algorithm for solving the traveling salesman
problem.
A sequence alignment between two sequences n and m
bases long is represented by a n x m Needleman-wunsch
scoring matrix. At the beginning of every generation, a set
number of ants start at the bottom right corner (which represents the end of an alignment) of the matrix, and traverse
the matrix until they reach the beginning of the alignment
(when the ant reaches either the top or left edge of the matrix). From any given position, an ant can move in one of
three directions: up, left, or diagonally up-left. As is the case
with the Needleman-Wunsch algorithm, a move in the former two directions represents a gap in one of the strands,
while a diagonal move is either a match or a mismatch.
At every step, the ants use a scoring algorithm to decide
what direction to move in. Every direction (up, left, or diagonal) is assigned a score based on the following formula:
Wφ

φi

· M Wm · R Wr

Where φi is the pheromones level for step i in a particular direction, M has a value of 2 if going in that direction
will result in a match and 1 if it results in a mismatch, R is
the regional weighting score that encourages ants to move
toward the diagonal of the matrix, and Wφ , Wm , and Wr
are the weights assigned to those values. The ant chooses a
number pi . If pi is greater than the parameter p, then the ant
simply moves in the direction with the highest score. If pi
is less than p, the ant chooses the direction using a weighted
probability based on the scores.
After each ant moves, it leaves a pheromone trail according to the formula.
φi+1 = (φi + φstep ) · ql
Where φi is the old pheromone level, phii+1 is the new
pheromone level, and ql is a number between 0 and 1 that
represents the local decay of pheromones.
As each ant moves, it also scores its path using a
Needleman-Wunsch like scoring scheme that adds points for
matches and subtracts points for mismatches and gaps. At
the end of each generation, the best path is reinforced according to the formula.
φi+1 = (φi + φstep · (Sg /Sh )
2

Variable
g
a
φ0
φstep
Wφ
Wm
Wr
ql
qg
p

Range
10 - 40
5 - 30
0-1
0-1
0 - 10
0 - 10
0 - 10
0-1
0-1
0-1

In calculating the output of the ACO given a set of parameters, we first showed empirically that the distribution
of scores followed a normal distribution. This allowed us to
easily throw out the data points that were outside x ± σ and
calculate a corrected mean (x corr ) during the evolution of the
parameters. This prevents random extreme values from affecting the mean and is a valid procedure after verifying that
the output from the ACO is normally distributed. We also
measured the total time to run all the trials of for each parameter (ttrials ). Using these two variables we can assign a
score to each parameter (Pi ) 1 :
Score(Pi ) =

Table 1: Variables optimized for in the genetic algorithm.

x corr (i)3
ttrials (i)

(1)

Using this scoring function, we were able to make a relative ranking of the current population of the highest scoring
parameters. We then kept the top 1% of the population
as “parents.” The next generation was then generated using
various parents and copy mutations (see Appendix to see the
actual implementation of the genetic algorithm).
Several runs for each string length was preformed, in
hopes of finding a set of parameters that came up often
which would represent a locally optimal solution. After the
best parameters for each of the strings were found, we can
interpolate and approximate the best set of parameters for
any string length up to a length of 100 characters.

where Sg is the score of the best path of this generation, and
Sh is the best score obtained so far overall. The entire matrix
also undergoes a global decay:
φi+1 = φi · qg
where qg is the global decay parameter.
The simulation keeps on looping through generations until the score for the best alignment converges, or until the
maximum number of generations is reached.

Genetic Algorithm

Results

The motivation for using a genetic algorithm is simple: all
the parameters above can be approximately determined using a genetic algorithm which evolves for the best set of parameters. Since hand-tweaking of ACO parameters can be a
futile exercise and often a waste of time, genetic algorithms
provide a convenient alternative albeit it may sometimes be
inaccurate. We used a very simple implementation of genetic algorithms to quickly find a set of parameters that were
locally optimal at discrete steps of search string length.
At the beginning of every generation, two strings are randomly generated. The first string is the template string
which is simply a random string filled with four possible
values (0, 1, 2, 3) at every given position. The second string
is generated by performing approximately 33% to 66% random mutations of the template string. This includes insertions of a random value, deletion at a given position, and a
point mutation of changing one value to another value at a
given point in the string. This generation process hopefully
produces two strings that are mostly similar but at the same
time contain significant diﬀerences in their sequences.
Using a fixed string length size for the two similar strings
that will be aligned, we evolved a best set of parameters.
Table 1 shows the variables and their possible range. We then
used a series of string lengths (10, 20, 30, 40, 50, 60, 70, 80,
90, 100) and attempted to evolve the best set of parameters
for each of these parameters.

Ant Colony Optimization Results
Using the same string generators for the genetic algorithm,
a string length of 50 was set. Then a random set of parameters were chosen. Afterwards, the ACO was run with these
parameters and the chosen strings for 100 repetitions. The
outputting scores were then gathered and compiled into a
histogram. This process was repeated several times. Figure
1 shows one such histogram. All of the outputs from the
ACOs followed a similar distribution to the one shown in
Figure 1.
Hence we were able to conclude that to a first order approximation, we can expect the ACO output to be normally
distributed. Given that the distribution is normal, we use
this knowledge to save time during the genetic algorithm
run and to make a good approximation of the mean score
given by a set of parameters as was shown above.

Genetic Algorithm Results
We explored the range of 10 - 100 characters for the string
length. Our genetic algorithm probed this string length
1 It was found that if x
corr is not more emphasized than time, then the
GA will start sacrificing score for faster and faster times. This is an undesirable result; raising to the third power provided better results.

3

range by attempting to find the best set of parameters in
10 character increments: we ran a separate genetic algorithm
for the string lengths 10, 20, 30, 40, 50, 60, 70, 80, 90,
and 100 characters. Each of these string lengths were run
multiple times with a population of 500 parameter sets, and
the genetic algorithm was stopped after the population had
reached an equilibrium for about 10 generations.
Figure 2 and Figure 3 show the trends for each of the variables.2 It is clear from these trends that most of the variables
follow a linear trend. The only exceptions are g and a which
both seem to require a lower limit and then proceed to grow
starting at 80 string length. (There is not enough data to
support what kind of growth the trend follows.) All other
variables increase linearly, decrease linearly, or remain constant.
Using these variables we were able to interpolate values
for any string length between 10 to 100. Using this system of
approximation, we were able to design an optimized form of
the ACO implementation that first calculated approximate
optimal values for the parameter set and then runs the ACO
on the input strings.

14
12
10

+

+

+|

+|

|

|

|

|

+

+

+

+

8

|

+

×|

+

×
|

6

|

×

×

×

×

×

×

×

×

4
2
0
0

20

40

60

80

100

String Length
Figure 2: This figure shows the trends for the variables in
the range of 0 to 15.
represents g. × represents a. +
|
represents Wφ . represents Wm . represents Wr .
provided a sanity check to make sure that the our algorithm
did indeed provide a satisfactory alignment.
In testing unalignable pairs of sequences, we choose sequences that were reversed of each other and no character
appears twice in these sequences. Type II alignments provide a test to see how the ACO deals with pairs of sequences
that it had not evolved to handle.
In order to test Type III, we generate a pair of sequences
in the similar manner as we did during the evolutions of the
parameters for genetic algorithms. Theoretically, it is these
sequences that should provide the most interesting alignments.
Table 2 shows an example alignment for each one of these
categories. There is not an objective measurement that
would provide a quantitative degree of similarity between
the alignments, hence we only provide a representative alignment of each of the categories instead.

Optimized ACO Results
Using the parameter set given above, we modified our ACO
implementation to use the parameter table and interpolate approximately optimal parameters for a given input
string (See Appendix A for the actual implementation). We
provide here some of the sample alignments outputted by
our implementation and compare them to the NeedlemanWunsch sequence alignment.
We divide the category of pairs of sequences into one of
three types: highly similar (Type I), unalignable (Type II),
and in between (Type III). This provides a qualitative basis
for comparison between the Needleman-Wunsch and our
ACO implementation.
To test highly similar pairs of sequences, we choose sequences that were nearly the same. Thus Type I alignments
2 Note: variables with similar ranges were grouped together to save
space.

Type
I

60
50
40
30

II
III

20
10
0

ACO

Needleman-Wunsch

abcdefgggghijklmnopq
|||||||
||||||||||
abcdefg---hijklmnopq
qpo-nmlkjihgfedcba
|
-abcdefghijklmnopq
-CACTTTTTCAGATCTATTG
|||||||||||| ||||
CTACTTTTTCAGATATATTC

abcdefgggghijklmnopq
|||||||
||||||||||
abcdefg---hijklmnopq
qponmlkjihgfedcba
|
abcdefghijklmnopq
C-ACTTTTTCAGATCTATTG
| |||||||||||| ||||
CTACTTTTTCAGATATATTC

Table 2: Sample alignments from both ACO and
Needleman-Wunsch. Type I are pairs of sequences that are
highly similar. Type II are pairs of sequences that are unalignable. Finally, Type III are pairs of sequences that are
somewhere in between.

45 47 49 51 53 55 57 59 61 63 65 67 69 71
Figure 1: A histogram of output of one ACO run with same
set of parameters for a 100 runs.
4

+

1.0

|

+

0.8

|

+

|

+

|

+

|

|

|

|

|

|

250

+
+

+
+

0.6

+

+

0.4

×

×

×

×

×
×

200
150

×
×

300

×

×

100
50

0.2

0

0
0

20

40

60

80

100

−50

String Length

++ + +
+++++ + +
+ +
++ ++ +++++++ ++++++
+++++++ + +
+++++++
+++ ++++++++++
+
+
+
+
+
+
+
+
+
+
+
+
+ +
+++++ +
+++++++++++++++++++++
++ +++ +++
++++++++++ +++
++++ ++++++++++ +
+
+
+
+
+
+
+
+ ++ +++++++++++++++++++++++++++++++++++++++
+
+ ++++++++
++++
+++++++++++++
++++++++
+++
+
+++++++++
++++
+++++++++
+++++
+++++++ +
++++++++++++++++
++++++++++
++++++
++++++++++++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ +++++++++++++++++++
++++++++++++++++++++
++++++
+++
++++
+++++++++++++++
+++
+++++++ ++ ++
++++
+++
+++++
+++++
++
++
+++++
+++++++++ ++++++++++
+++
+++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
++
+
+
+
+
+
++++ ++ + ++
+
+++++++++
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
++
+++++++++++++++++++ ++ ++ + + + ++
++
+
+++
+++++++++
+++
++++++ ++++++
+
+
+
+
+
+
+
+
+
+
+
+
+
+ +
++++
++++ + +
+
+++++++
+ ++
+++++
++
+++++++ + +++
++
++++++++ + + + ++ +++
++++++++++++
+
+
+
+
+
+
+
+
++ + + +
+
+

50

100

150

200

250

300

350

Figure 4: This is a scatter plot of scores from
both Needleman-Wunsch and ACO. The x-axis is the
Needleman-Wunsch scores, and the y-axis is the ACO
scores.

Figure 3: This figure shows the trends for the variables in
the range of 0 to 1. + represents φ0 . × represents φstep .
represents qg . represents qg . | represents p.

Some problems and sources of error were encountered
during the course of experimental runs of our implementation. We realize that there is an error in our genetic algorithm in that we should have optimized for high scoring consistency and low run time instead of optimizing for
high average score and low run time. This produced parameter sets for the ACO that were suboptimal for the problem
that we were addressing. Due to the time constraint of the
project, we were not able to fix this issue in time for the
deadline.
Also, the higher string lengths such as the 90 and 100
case were based on data that was not completely reliable in
that we did not have enough time to run the genetic algorithm multiple times with many generations due to the
highly computationally expensive nature of this algorithm.
This produced possibly suboptimal parameter sets for these
string lengths and hence cannot be as well trusted as the parameter sets generated for the lower part of the string length
range.
Finally, we realize that genetic algorithms have the danger
of producing suboptimal results, but we are confident in the
results we have obtained because multiple runs of the genetic
algorithm for a given string length gave results that were very
similar to each other, an indication that the parameter sets
generated were indeed near the global optimum.
A few improvements could be made to the ACO implementation for sequence alignment as well. The way that
it is currently implemented, there is a limitation in the accuracy of the outputted alignments. The ants crawling up
from the bottom right of the scoring matrix to the upper left
cause alignments to be very accurate at the beginning of the
traversal but not at the end of the alignment since the paths
of the ants will diverge towards the end. This limitation can
be easily overcome by adding an additional feature of having

Discussion
Analysis
At first glance, Table 2 provides some very promising results.
It shows that the ACO actually preforms at a level of reaching nearly the mathematically optimal solution. In Type I,
the alignments are exactly the same. For Type II, the ACO
alignment is suboptimal compared with the NeedlemanWunsch because there is an extra gap that is not necessary.
For Type III, the alignment is again nearly optimal: the main
diﬀerence comes at the beginning of the alignment where
the ACO could have made an extra match. It is interesting
to note that on several other alignments, the main diﬀerence between ACO and Needleman-Wunsch comes at the
beginning, a matter that we will address later.
Hence even from a completely qualitative standpoint, it
can be said that the ACO does provide satisfactory results.
While it can not be considered to be faster or more eﬃcient
than Needleman-Wunsch, it can, nevertheless, be applied to
sequence alignment.

Conclusion
As the results of the optimized ACO runs have shown, the
ACO algorithm is a powerful algorithm that can be readily
applied to sequence alignment. While it is true that this particular implementation of the ACO with respect to sequence
alignment is not algorithmically faster or more eﬃcient than
the Needleman-Wunsch dynamic programming equivalent,
the purpose of this paper was to show that the ACO is a
viable algorithm that can be used to solve computational biology problems perhaps more eﬃciently when searching an
NP hard search space.
5

non-scoring ants dropped in random places on the scoring
matrix. This way, areas around where the ants were dropped
will have better alignment paths for when the scoring ants
which started in the lower left corner reached the paths that
were reinforced by the non-scoring ants. If local alignment
is desired, then the random ants can be simply considered
to be scoring.
In conclusion, we found that the ACO algorithm is an
algorithm that can be used to solve biological problems.
While in this case, it is not necessarily better than the existing algorithm for two sequence alignment; we are confident that the ACO can be readily adapted for much more
complicated problems such as Multi-Sequence Alignment
(MSA). In this case, the search space is much larger and the
Needleman-Wunsch algorithm is out of the question; hence,
a good adaption of the ACO meta-heuristic might provide
a solution that rivals even the fastest MSA implementation
today. Additional applications of ACO to biological problems could be such problems as clustering or even an ACO
adapted algorithm for protein folding.

Contributions
Thanks to Qian Zhang for his suggestions on ways of improving our algorithm. To Charles Duan for helping us with
LATEX. To the Harvard NICE computing environment. To
our TF, Jon Radoﬀ for answering our never-ending questions. And to Professor Church for his instruction of this
class.

References
[1] Dorigo, M. and Gambardella, L.M. (1997). Ant Colony System: A Cooperative Learning Approach to the Traveling Salesman Problem. IEEE Trans. Evol. Comp. 1, 53-66.
[2] Lawler, E.L., Lenstra, J.K., Rinnooy-Kan, A.H.G., and
Shmoys, D.B. (eds) (1985). The Travelling Salesman Problem.
New York, NY: Wiley.
[3] Gambardella, L.M., Taillard, E.D. and Dorigo, M. (1998). Ant
Colonies for the QAP J. Operational Res. Soc. 1998. (in press).
[4] Bottee, H. and Bonabeau, E. (1998). Evolving Ant Colony Optimization Adv. Complex Systems 1998. 1, 149-159
[5] Colorni, A., Dorigo, M., Maniezzo, V. and Trubian, M. (1993).
Ant system for job-shop scheduling. Belg. J. Oper. Res., Stat.
and Comput. Sci. 34, 39-53.

6

Appendix
A ACO Implementation
1
2
3
4

#!/usr/local/bin/perl
use strict;
# Ant.pl
# Authors: Li-Wei King, Aaron Lee

5
6
7
8
9
10
11
12
13
14
15
16
17
18

###################################
# Default values just in case something goes terribly wrong
my $NUM_ANTS
= 10;
my $MAX_GEN
= 50;
my $PHER_STEP
= .1;
my $PHER_WEIGHT = .1;
my $MATCH_WEIGHT = .1;
my $REGION_WEIGHT = .00001;
my $PROB_PROB
= .5;
# Must be between 0 and 1. It is the probability
#
that the direction will be chosen based on probability.
my $INIT_PHER = .01;
# The amount of pheremone that an untraveled path has.
my $L_DECAY = .7;
# Local decay. Cannot be larger than 1.
my $G_DECAY = .7;
# Global decay. Cannot be larger than 1.

19
20
21
22
23

# Constants that won’t eventually be passed to the program.
use constant MATCH_BONUS => 5;
use constant GAP_PENALTY => -4;
use constant MISMATCH_PENALTY => -3;

24
25
26
27
28
29
30
31

############
# Useful abbreviations
#For encoding the direction an ant moves
use constant D => 0;
use constant L => 1;
use constant U => 2;
use constant N => -1;
# Undefined direction

32
33
34
35
36
37

# The index at which certain information is stored in the $ants array
use constant X_INDEX => 0;
use constant Y_INDEX => 1;
use constant DIR_INDEX => 2;
# Where the ant decides to move from that position
use constant SCORE_INDEX => 3; # The current score at that position

38
39
40
41
42
43
44
45

###########
# Table of GA optimized
my @GA_TABLE;
$GA_TABLE[10][0] = 10;
$GA_TABLE[40][0] = 10;
$GA_TABLE[70][0] = 10;
$GA_TABLE[100][0] = 15;

values
$GA_TABLE[20][0]
$GA_TABLE[50][0]
$GA_TABLE[80][0]

= 10; $GA_TABLE[30][0]
= 10; $GA_TABLE[60][0]
= 10; $GA_TABLE[90][0]

= 10;
= 10;
= 12;

46
47
48
49
50

$GA_TABLE[10][1]
$GA_TABLE[40][1]
$GA_TABLE[70][1]
$GA_TABLE[100][1]

=
=
=
=

5; $GA_TABLE[20][1]
5; $GA_TABLE[50][1]
5; $GA_TABLE[80][1]
10;

= 5; $GA_TABLE[30][1]
= 5; $GA_TABLE[60][1]
= 5; $GA_TABLE[90][1]

= 5;
= 5;
= 8;

$GA_TABLE[10][2]
$GA_TABLE[30][2]
$GA_TABLE[50][2]
$GA_TABLE[70][2]
$GA_TABLE[90][2]

=
=
=
=
=

.411191490;
.453581583;
.432201854;
.417636990;
.341437549;

$GA_TABLE[20][2]
$GA_TABLE[40][2]
$GA_TABLE[60][2]
$GA_TABLE[80][2]
$GA_TABLE[100][2]

=
=
=
=
=

.438349294;
.517059770;
.436950953;
.366514982;
.329430526;

$GA_TABLE[10][3]
$GA_TABLE[30][3]
$GA_TABLE[50][3]
$GA_TABLE[70][3]
$GA_TABLE[90][3]

=
=
=
=
=

9.434392207;
9.353249096;
9.290221874;
9.149204714;
9.332940631;

$GA_TABLE[20][3]
$GA_TABLE[40][3]
$GA_TABLE[60][3]
$GA_TABLE[80][3]
$GA_TABLE[100][3]

=
=
=
=
=

9.423857194;
9.284172996;
9.282690356;
9.064648465;
9.259328124;

$GA_TABLE[10][4]
$GA_TABLE[30][4]
$GA_TABLE[50][4]
$GA_TABLE[70][4]
$GA_TABLE[90][4]

=
=
=
=
=

6.109820365;
9.343075608;
10;
10;
10;

$GA_TABLE[20][4]
$GA_TABLE[40][4]
$GA_TABLE[60][4]
$GA_TABLE[80][4]
$GA_TABLE[100][4]

=
=
=
=
=

6.926580738;
9.244311660;
10;
10;
10;

$GA_TABLE[10][5]
$GA_TABLE[30][5]
$GA_TABLE[50][5]
$GA_TABLE[70][5]

=
=
=
=

3.909960135;
2.224402772;
2.142958734;
1.736982155;

$GA_TABLE[20][5]
$GA_TABLE[40][5]
$GA_TABLE[60][5]
$GA_TABLE[80][5]

=
=
=
=

2.350289525;
1.853908945;
1.915834968;
1.878705913;

51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

7

74

$GA_TABLE[90][5]

= 1.856936304; $GA_TABLE[100][5] = 1.862138526;

$GA_TABLE[10][6]
$GA_TABLE[30][6]
$GA_TABLE[50][6]
$GA_TABLE[70][6]
$GA_TABLE[90][6]

=
=
=
=
=

0.853763237;
0.827315180;
0.798933405;
0.730894471;
0.622498305;

$GA_TABLE[20][6]
$GA_TABLE[40][6]
$GA_TABLE[60][6]
$GA_TABLE[80][6]
$GA_TABLE[100][6]

=
=
=
=
=

0.830586273;
0.827352487;
0.728488635;
0.620157623;
0.628942392;

$GA_TABLE[10][7]
$GA_TABLE[30][7]
$GA_TABLE[50][7]
$GA_TABLE[70][7]
$GA_TABLE[90][7]

=
=
=
=
=

0.660878498;
0.606102419;
0.571805707;
0.579124731;
0.519128179;

$GA_TABLE[20][7]
$GA_TABLE[40][7]
$GA_TABLE[60][7]
$GA_TABLE[80][7]
$GA_TABLE[100][7]

=
=
=
=
=

0.635274652;
0.595236898;
0.577967010;
0.532124853;
0.515925041;

$GA_TABLE[10][8]
$GA_TABLE[30][8]
$GA_TABLE[50][8]
$GA_TABLE[70][8]
$GA_TABLE[90][8]

=
=
=
=
=

0.917907684;
1;
1;
1;
1;

$GA_TABLE[20][8]
$GA_TABLE[40][8]
$GA_TABLE[60][8]
$GA_TABLE[80][8]
$GA_TABLE[100][8]

=
=
=
=
=

1;
0.965550166;
1;
1;
1;

$GA_TABLE[10][9]
$GA_TABLE[30][9]
$GA_TABLE[50][9]
$GA_TABLE[70][9]
$GA_TABLE[90][9]

=
=
=
=
=

0.990544051;
1;
1;
1;
1;

$GA_TABLE[20][9]
$GA_TABLE[40][9]
$GA_TABLE[60][9]
$GA_TABLE[80][9]
$GA_TABLE[100][9]

=
=
=
=
=

1;
1;
1;
1;
1;

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

####################################
if ($#ARGV != 1) {
print "Usage: \n";
print "
ant.pl <SEG1> <SEG2>\n";
exit;
}

107
108
109

my $SEG1 = $ARGV[0];
my $SEG2 = $ARGV[1];

110
111
112
113
114
115
116

my $strlenavg = length($SEG1) + length($SEG2) / 2;
if ($strlenavg < 10 || $strlenavg > 100) {
print "The string lengths of SEG1 and SEG2 must have an average string
length betweeen 10 and 100!\n";
exit;
}

117
118

optimize_vars_run_aco($SEG1, $SEG2, $strlenavg);

119
120
121
122
123
124
125
126
127

sub optimize_vars_run_aco {
my $str1 = shift;
my $str2 = shift;
my $len = shift;
my @opt;
my $lowten = (int $len / 10) * 10;
my $highten = (int $len / 10) * 10 + 10;
my $wherex = $len - $lowten;

128

for (my $i = 0; $i <= 9; $i++) {
push(@opt, $wherex * ($GA_TABLE[$highten][$i]-$GA_TABLE[$lowten][$i])/10 + $GA_TABLE[$lowten][$i]);
}

129
130
131
132

run_aco($str1, $str2, $opt[0], $opt[1], $opt[2], $opt[3], $opt[4], $opt[5],
$opt[6], $opt[7], $opt[8], $opt[9]);

133
134
135

}

136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

sub run_aco {
# Arguments to the array
my $seq1 = shift;
my $seq2 = shift;
my $MAX_GEN = shift;
my $NUM_ANTS = shift;
my $PHER_STEP = shift;
my $PHER_WEIGHT = shift;
my $MATCH_WEIGHT = shift;
my $REGION_WEIGHT = shift;
my $INIT_PHER = shift;
my $L_DECAY = shift;
my $G_DECAY = shift;
my $prob_prob = shift;

151
152
153

#Other variables
my $high_score;

8

154
155
156
157
158

my
my
my
my
my

@best_path_x;
@best_path_y;
@best_path_dir;
$prev_score = 0;
# Used to figure out if there are repeats or not
$repeat_counter = 0;

159
160

print "Aligning: $seq1 with $seq2\n";

161
162
163
164

# Store the sequences in arrays
my @seqy = split(’’, $seq1);
my @seqx = split(’’, $seq2);

165
166
167
168
169
170
171
172
173
174
175

my @pher;
# Stores the pheremone values
# begin matrix initialization.
for (my $i = 0; $i <= $#seqy; $i++) {
for (my $j = 0; $j <= $#seqx; $j++) {
$pher[$i][$j][D] = $INIT_PHER;
$pher[$i][$j][L] = $INIT_PHER;
$pher[$i][$j][U] = $INIT_PHER;
}
}

176
177
178
179
180

# begin main generation loop
for (my $gen = 0; $gen < $MAX_GEN; $gen++) {
my %done = ();
# This hash stores a 1 if the ant is finished, 0 if it isn’t
my @ants = ();
# Stores the information for all the ants

181
182
183
184
185
186
187
188
189
190
191
192
193
194

# reinitialize ants. The ant array is organized thus:
# $ants[ant_id][move_number][information for that move]
for (my $ant = 0; $ant < $NUM_ANTS; $ant++) {
$ants[$ant][0][X_INDEX] = $#seqx;
$ants[$ant][0][Y_INDEX] = $#seqy;
$done{$ant} = 0;
if ($seqx[$#seqx] eq $seqy[$#seqy]) {
#$ants[$ant][0][SCORE_INDEX] = MATCH_BONUS;
$ants[$ant][0][SCORE_INDEX] = 0;
} else {
$ants[$ant][0][SCORE_INDEX] = 0;
}
}

195
196
197

my $ants_done = 0;
for (my $step = 0; $ants_done < $NUM_ANTS; $step++) {

198
199
200

# simulate the ants:
for (my $ant = 0; $ant < $NUM_ANTS; $ant++) {

201
202
203
204
205

#If the ant is done, then skip it.
if ($done{$ant} == 1) {
next;
}

206
207
208
209

# pick up where we left off:
my $ant_x = $ants[$ant][$step][X_INDEX];
my $ant_y = $ants[$ant][$step][Y_INDEX];

210
211
212
213
214
215

# get the direction that the ant should take, and update
# the direction index
my $dir = get_dir($ant_x, $ant_y,\@pher, $PHER_WEIGHT, $MATCH_WEIGHT,
$REGION_WEIGHT, $prob_prob, \@seqx, \@seqy);
$ants[$ant][$step][DIR_INDEX] = $dir;

216
217
218
219
220

# update pher matrix
$pher[$ant_y][$ant_x][$dir] += $PHER_STEP;
# Do local decay
$pher[$ant_y][$ant_x][$dir] *= $L_DECAY;

221
222
223
224
225
226
227
228
229
230
231
232
233

# move in the chosen direction and update the score.
# The score stored at each square includes
# match/mismatch score for this square as well as well
# as any gap penalty for the move that got the ant to
# that position .
if ($dir == D) {
$ants[$ant][$step + 1][X_INDEX] = $ant_x - 1;
$ants[$ant][$step + 1][Y_INDEX] = $ant_y - 1;
$ants[$ant][$step + 1][SCORE_INDEX] = $ants[$ant][$step][SCORE_INDEX];
if ($seqx[$ants[$ant][$step][X_INDEX]] eq $seqy[$ants[$ant][$step][Y_INDEX]]) {
$ants[$ant][$step + 1][SCORE_INDEX] += MATCH_BONUS;
} else {

9

$ants[$ant][$step + 1][SCORE_INDEX] += MISMATCH_PENALTY;

234

}

235
236

} elsif ($dir == U) {
$ants[$ant][$step
$ants[$ant][$step
$ants[$ant][$step
} elsif ($dir == L) {
$ants[$ant][$step
$ants[$ant][$step
$ants[$ant][$step
}

237
238
239
240
241
242
243
244
245

+ 1][X_INDEX] = $ant_x;
+ 1][Y_INDEX] = $ant_y - 1;
+ 1][SCORE_INDEX] = $ants[$ant][$step][SCORE_INDEX] + GAP_PENALTY;
+ 1][X_INDEX] = $ant_x - 1;
+ 1][Y_INDEX] = $ant_y;
+ 1][SCORE_INDEX] = $ants[$ant][$step][SCORE_INDEX] + GAP_PENALTY;

246

#Check if the ant is done and update score if needed
if ($ants[$ant][$step + 1][X_INDEX] == 0 || $ants[$ant][$step + 1][Y_INDEX] == 0) {
$ants[$ant][$step + 1][DIR_INDEX] = N;
$done{$ant} = 1;
if ($seqx[$ants[$ant][$step+1][X_INDEX]] eq $seqy[$ants[$ant][$step+1][Y_INDEX]]) {
$ants[$ant][$step + 1][SCORE_INDEX] += MATCH_BONUS;
} else {
$ants[$ant][$step + 1][SCORE_INDEX] += MISMATCH_PENALTY;
}

247
248
249
250
251
252
253
254
255
256

$ants_done++;

257

}

258

}

259
260

}

261
262
263
264
265
266
267
268

#Do a global update
#Find the best score
my $max_gen_score = $ants[0][$#{$ants[0]}][SCORE_INDEX];
my $max_score_ant = 0;
my @best_gen_path_x = ();
# The best path of the generation:
my @best_gen_path_y = ();
my @best_gen_path_dir = ();

x coordinates

269
270
271
272
273
274
275

for (my $a = 1; $a <= $#ants; $a++) {
if ($ants[$a][$#{$ants[$a]}][SCORE_INDEX] > $max_gen_score) {
$max_gen_score = $ants[$a][$#{$ants[$a]}][SCORE_INDEX];
$max_score_ant = $a;
}
}

276
277
278
279
280
281
282

#Get the best path
for(my $s = 0; $s <= $#{$ants[$max_score_ant]}; $s++) {
push (@best_gen_path_x, $ants[$max_score_ant][$s][X_INDEX]);
push (@best_gen_path_y, $ants[$max_score_ant][$s][Y_INDEX]);
push (@best_gen_path_dir, $ants[$max_score_ant][$s][DIR_INDEX]);
}

283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302

#Find out if the high score this time is the all time high.
if ($gen == 0) {
$high_score = $max_gen_score;
for (my $i = 0; $i <= $#best_gen_path_x; $i++) {
push (@best_path_x, $best_gen_path_x[$i]);
push (@best_path_y, $best_gen_path_y[$i]);
push (@best_path_dir, $best_gen_path_dir[$i]);
}
} elsif ($max_gen_score > $high_score) {
@best_path_x = ();
@best_path_y = ();
@best_path_dir = ();
$high_score = $max_gen_score;
for (my $i = 0; $i <= $#best_gen_path_x; $i++) {
push (@best_path_x, $best_gen_path_x[$i]);
push (@best_path_y, $best_gen_path_y[$i]);
push (@best_path_dir, $best_gen_path_dir[$i]);
}
}

303
304
305
306
307
308
309
310
311

#Reinforce the best path again based on how good it is. I’m
#using < rather than <= because the last direction is undefined.
if ($max_gen_score > 0 && $high_score > 0) {
for (my $i = 0; $i < $#best_gen_path_x; $i++) {
$pher[$best_gen_path_y[$i]][$best_gen_path_x[$i]][$best_gen_path_dir[$i]]
+= ($max_gen_score/$high_score) * $PHER_STEP;
}
}

312
313

#global_decay(\@pher, $G_DECAY);

10

314

#Check for convergence
if ($max_gen_score == $prev_score) {
$repeat_counter++;
} else {
$max_gen_score = $prev_score;
$repeat_counter = 0;
}

315
316
317
318
319
320
321
322

#If we get the same score 5 times, we’ll stop the calculation
if ($repeat_counter == 5) {
last;
}

323
324
325
326

}

327
328

my
my
my
my

329
330
331
332

@alignment = get_alignment(\@best_path_x, \@best_path_y, \@seqx, \@seqy);
$x_alignment = $alignment[0];
$dashes = $alignment[1];
$y_alignment = $alignment[2];

333

print
print
print
print
print

334
335
336
337
338

"The best alignment obtained was: \n";
"$x_alignment\n";
"$dashes\n";
"$y_alignment\n";
"The high score was: $high_score \n";

339
340

}

341
342
343
344
345
346
347

sub get_alignment {
#Arguments
my $best_x_ref = shift;
my $best_y_ref = shift;
my $seq_x_ref = shift;
my $seq_y_ref = shift;

# Stores the x coordinates of the best path

348
349
350
351
352
353
354

#Other variables
my $last_x;
my $last_y;
my $x_result = "";
my $y_result = "";
my $dashes = "";

# The x coordinate of the last step we traced

355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371

# Handle the cases where the alignment doesn’t start in the 0,0
# position
if($best_x_ref->[$#{@$best_x_ref}] != 0) {
for (my $i = 0; $i < $best_x_ref->[$#{@$best_x_ref}]; $i++) {
$x_result = $x_result.$seq_x_ref->[$i];
$y_result = $y_result."-";
$dashes = $dashes." ";
}
} elsif ($best_y_ref->[$#{@$best_y_ref}] != 0) {
for (my $i = 0; $i < $best_y_ref->[$#{@$best_y_ref}]; $i++) {
$x_result = $x_result."-";
$y_result = $y_result.$seq_y_ref->[$i];
$dashes = $dashes." ";
}
}

372
373
374
375
376
377
378

#First step. Start from the beginning of the alignment (the end of
# the best_x_ref and best_y_ref arrays
$last_x = $best_x_ref->[$#{@$best_x_ref}];
$last_y = $best_y_ref->[$#{@$best_y_ref}];
$x_result = $x_result.$seq_x_ref->[$last_x];
$y_result = $y_result.$seq_y_ref->[$last_y];

379
380
381
382
383
384

if ($seq_x_ref->[$last_x] eq $seq_y_ref->[$last_y]) {
$dashes = $dashes."|";
} else {
$dashes = $dashes." ";
}

385
386
387
388
389
390
391
392
393

for (my $i = $#{@$best_x_ref} - 1; $i >= 0; $i--) {
# Gap in the x if the x index has not changed. Append a gap to
# the sequence.
if ($best_x_ref->[$i] eq $last_x) {
$last_y = $best_y_ref->[$i];
$x_result = $x_result."-";
$y_result = $y_result."$seq_y_ref->[$best_y_ref->[$i]]";
$dashes = $dashes." ";

11

} elsif ($best_y_ref->[$i] eq $last_y) {
# Gap in y
$last_x = $best_x_ref->[$i];
$y_result = $y_result."-";
$x_result = $x_result."$seq_x_ref->[$best_x_ref->[$i]]";
$dashes = $dashes." ";
} else {
#Diagonal move
$last_x = $best_x_ref->[$i];
$last_y = $best_y_ref->[$i];
$y_result = $y_result."$seq_y_ref->[$best_y_ref->[$i]]";
$x_result = $x_result."$seq_x_ref->[$best_x_ref->[$i]]";
if ($seq_y_ref->[$best_y_ref->[$i]] eq $seq_x_ref->[$best_x_ref->[$i]]) {
$dashes = $dashes."|";
} else {
$dashes = $dashes." ";
}
}

394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411

}

412
413

return ($x_result, $dashes, $y_result);

414
415

}

416
417
418
419
420

sub global_decay {
#Arguments
my $pher_ref = shift;
my $g_decay = shift;

421

for (my $x = 0; $x <= $#{ $pher_ref }; $x++) {
for (my $y = 0; $y <= $#{ $pher_ref->[$x] }; $y++ ) {
for (my $d = 0; $d < 3; $d++) {
$pher_ref->[$x]->[$y]->[$d] *= $g_decay;
}
}
}

422
423
424
425
426
427
428
429

}

430
431
432
433
434
435
436
437
438
439
440
441

sub get_dir {
#Arguments
my $x = shift;
my $y = shift;
my $pher_ref = shift;
my $pher_weight = shift;
my $match_weight = shift;
my $region_weight = shift;
my $prob_prob = shift;
my $seqx = shift;
my $seqy = shift;

442
443
444
445

#Calculate the regional weighting for these coordinates.
#assumes that both sequences are the same length.
#TODO: MAKE SURE TO CHANGE REGIONAL WEIGHTINGS

This

446
447
448
449

my $reg_w_up;
my $reg_w_left;
my $reg_w_diag;

450
451
452
453
454
455
456
457
458
459
460
461
462
463

if ($x == $y) {
$reg_w_up = 1;
$reg_w_left = 1;
$reg_w_diag = 2;
} elsif ($x > $y) {
$reg_w_up = 1;
$reg_w_diag = 1.5;
$reg_w_left = 2;
} else {
$reg_w_up = 2;
$reg_w_left = 1;
$reg_w_diag = 1.5;
}

464
465
466
467
468

#Figure out whether going up, diagonal, or left will give us a match
my $match_up;
my $match_left;
my $match_diag;

469
470
471
472
473

if ($seqx->[$x] eq $seqy->[$y - 1]) {
$match_up = 2;
} else {
$match_up = 1;

12

}

474
475

if ($seqx->[$x - 1] eq $seqy->[$y - 1]) {
$match_diag = 2;
} else {
$match_diag = 1;
}

476
477
478
479
480
481

if ($seqx->[$x - 1] eq $seqy->[$y]) {
$match_left = 2;
} else {
$match_left = 1;
}

482
483
484
485
486
487

#Figure out the scores for each direction
my $up_score;
my $diag_score;
my $left_score;

488
489
490
491
492

$up_score = exp(log($pher_ref->[$y]->[$x]->[U]) * $pher_weight
+ log($match_up) * $match_weight
+ log($reg_w_up) * $region_weight);

493
494
495
496

$diag_score = exp(log($pher_ref->[$y]->[$x]->[D]) * $pher_weight
+ log($match_diag) * $match_weight
+ log($reg_w_diag) * $region_weight);

497
498
499
500

$left_score = exp(log($pher_ref->[$y]->[$x]->[L]) * $pher_weight
+ log($match_left) * $match_weight
+ log($reg_w_left) * $region_weight);

501
502
503
504

#Decide whether we’re gonna decide the direction based on the best
#score, or based on a weighted probability.
#
if (rand(1) > $prob_prob) {
#If we’re just gonna do it based on best score

505
506
507
508
509
510

my $high_score = max($up_score, $diag_score, $left_score);

511
512

if ($high_score == $up_score) {
return U;
} elsif ($high_score == $diag_score) {
return D;
} else {
return L;
}
} else {
#Determine the direction using probabilities

513
514
515
516
517
518
519
520
521
522

my $total_score = $up_score + $diag_score + $left_score;

523
524

my $rand_num = rand($total_score);

525
526

$total_score -= $up_score;
if ($rand_num > $total_score) {
return U;
}

527
528
529
530
531

$total_score -= $diag_score;
if ($rand_num > $total_score) {
return D;
} else {
return L;
}

532
533
534
535
536
537

}

538
539

}

540
541

sub max {

542

my $max = $_[0];

543
544

for(my $i = 1; $i <= $#_; $i++) {
if ($_[$i] > $max) {
$max = $_[$i];
}
}

545
546
547
548
549
550

return $max;

551
552

}

13

B Genetic Algorithm Implementation
1
2
3

#!/usr/local/bin/perl -w
use Time::HiRes;
use strict;

4
5

### BEGIN VARIABLES

6
7
8
9
10
11
12
13
14

# Constants
my $STR_LENGTH
my $POP_LIMIT
my $RECORD_FREQ
my $NUM_TEST_IND
my $KEEP_PARENTS
my $CROSSOVER
my $ZERO

=
=
=
=
=
=
=

20;
1000;
5;
7;
int ($POP_LIMIT * .1);
.2;
0.0000000001;

15
16
17
18
19
20

# Global variables
my @pop;
my $cur_gen
my $file_num
my $SCORE_INDEX = 10;

= 0;
= 0;

21
22

my @range = ();

23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

$range[0][0]
$range[0][1]
$range[1][0]
$range[1][1]
$range[2][0]
$range[2][1]
$range[3][0]
$range[3][1]
$range[4][0]
$range[4][1]
$range[5][0]
$range[5][1]
$range[6][0]
$range[6][1]
$range[7][0]
$range[7][1]
$range[8][0]
$range[8][1]
$range[9][0]
$range[9][1]

=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=

10;
40;
5;
30;
$ZERO;
1;
$ZERO;
10;
$ZERO;
10;
$ZERO;
5;
$ZERO;
1;
$ZERO;
1;
$ZERO;
1;
$ZERO;
1;

# $MAX_GEN
# $NUM_ANTS
# $PHER_STEP
# $PHER_WEIGHT
# $MATCH_WEIGHT
# $REGION_WEIGHT
# $INIT_PHER
# $L_DECAY
# $G_DECAY
# $PROB_PROB

44
45
46

### BEGIN MAIN PROGRAM

47
48
49

#randomize variables
randomize_vars();

50
51
52
53
54

while (1) {
#initialization
my $seq1 = getrandseq();
my $seq2 = mutate($seq1);

55

print
print
print
print

56
57
58
59

"Generation $cur_gen\n";
"-Testing sequences:\n";
"-Sequence 1 (".length($seq1)."): $seq1\n";
"-Sequence 2 (".length($seq2)."): $seq2\n";

60
61

test_pop($seq1, $seq2);

62
63

gen_next_pop();

64
65

if ($cur_gen == $RECORD_FREQ) {
print_pop();
$cur_gen = 0;
$file_num++;
}

66
67
68
69
70
71

$cur_gen++;

72
73

}

74
75

### BEGIN SUBROUTINES

76
77

sub getrandseq {

14

my $limit = $STR_LENGTH;
my $str = "";

78
79
80

for (my $i = 0; $i < $limit; $i++) {
my $aa = int rand(4);

81
82
83

$str .= $aa;

84

}

85
86

return $str;

87
88

}

89
90
91
92

sub mutate {
my $template = shift;
my $num_muts = int length($template) / (rand (2) + 1.5);

93

for (my $i = 0; $i < $num_muts; $i++) {
my $whichmut = int rand(3);
my $where = int rand(length($template));

94
95
96
97

if ($whichmut == 0) {
# Point mutation
substr($template,$where,1,int rand(4));
} elsif ($whichmut == 1) {
# Insertion
my $beg = substr($template,0,$where);
my $end = substr($template,$where);
my $num = int rand(4);
$template = "$beg$num$end";
} else {
# Deletion
my $beg = substr($template,0,$where);
my $end = substr($template,$where+1);
$template = "$beg$end";
}

98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

}

113
114

return $template;

115
116

}

117
118

sub randomize_vars {

119

for (my $i = 0; $i < $POP_LIMIT; $i++) {
for (my $j = 0; $j <= $#range; $j++) {
$pop[$i][$j] = rand($range[$j][1]-$range[$j][0]) + $range[$j][0];
}
$pop[$i][0] = int $pop[$i][0];
$pop[$i][1] = int $pop[$i][1];
}

120
121
122
123
124
125
126
127

}

128
129
130
131

sub test_pop {
my $seq1 = shift;
my $seq2 = shift;

132

my @winners;

133
134

print "-- Beginning population testing...\n";
for (my $i = 0; $i < $POP_LIMIT; $i++) {
my $score = score_individual($i, $seq1, $seq2);
$pop[$i][$SCORE_INDEX] = $score;
}

135
136
137
138
139
140

# Sort so that the winners are at the top and losers are at the bottom
@pop = sort {$b->[$SCORE_INDEX] <=> $a->[$SCORE_INDEX]} @pop;

141
142
143

print "-- DONE! \n";
print "-- Top three scores: ";
print "$pop[0][$SCORE_INDEX], $pop[1][$SCORE_INDEX], $pop[2][$SCORE_INDEX]\n";

144
145
146
147

}

148
149
150
151
152

sub score_individual {
my $ind = shift;
my $seq1 = shift;
my $seq2 = shift;

153
154

my @scores = ();

155
156
157

my $totalscore = 0;
my $totaltime = 0;

15

158

print "-- Testing string length: $STR_LENGTH\n";
print "-- Generation: ".($file_num*$RECORD_FREQ+$cur_gen)."\n";
print "-- Testing individual $ind \n";
print "----- Ants : $pop[$ind][1] \n";
print "----- Gens : $pop[$ind][0] \n";
for (my $i = 0; $i < $NUM_TEST_IND; $i++) {
my $score;
my $before = Time::HiRes::time();
$score = run_aco($seq1, $seq2, int $pop[$ind][0], int $pop[$ind][1],
$pop[$ind][2], $pop[$ind][3], $pop[$ind][4],
$pop[$ind][5], $pop[$ind][6], $pop[$ind][7],
$pop[$ind][8], $pop[$ind][9]);
$totaltime += Time::HiRes::time() - $before;
$totalscore += $score;
push(@scores, $score);
}

159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175

my $mean = $totalscore / $NUM_TEST_IND;
print "------- Mean:$mean\n";

176
177
178

my $stddev = 0;
for (my $i = 0; $i < $NUM_TEST_IND; $i++) {
my $term = ($mean - $scores[$i])*($mean - $scores[$i]);
$stddev += $term;
}
$stddev = sqrt($stddev / ($NUM_TEST_IND + 1));
print "------- Standard dev:$stddev\n";

179
180
181
182
183
184
185
186

my $newmean = 0;
my $datapoints = 0;
my $thrownout = 0;
for (my $i = 0; $i < $NUM_TEST_IND; $i++) {
if ($scores[$i] >= $mean - $stddev && $scores[$i] <= $mean + $stddev) {
$newmean += $scores[$i];
$datapoints++;
} else {
$thrownout++;
}
}
if ($datapoints == 0) {
$newmean = -10000;
} else {
$newmean /= $datapoints;
}

187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203

my $ga_score =
print "------print "------print "------print "------print "-------

204
205
206
207
208
209

($newmean * $newmean * $newmean) / $totaltime;
Good datapoints: $datapoints\n";
Thrown out data: $thrownout\n";
Corrected mean: $newmean\n";
Average time: ".$totaltime/$NUM_TEST_IND."\n";
Final GA score: ".$ga_score."\n";

210

return ($ga_score);

211
212

}

213
214
215
216
217

sub gen_next_pop {
# Populate the next generation
for (my $i = $KEEP_PARENTS; $i < $POP_LIMIT; $i++) {
my $cur_parent = int rand($KEEP_PARENTS);

218
219
220

for (my $j = 0; $j < $SCORE_INDEX; $j++) {
$pop[$i][$j] = $pop[$cur_parent][$j];

221
222
223
224
225
226

# Copy mutation
my $which_mut = rand(10);
if ($which_mut < 2) {
# Plus case
$pop[$i][$j] += .1 * rand(1) * ($range[$j][1] - $range[$j][0]);

227
228
229
230
231
232
233

if ($pop[$i][$j] > $range[$j][1]) {
$pop[$i][$j] = $range[$j][1];
}
} elsif ($which_mut < 4) {
# Minus case
$pop[$i][$j] -= .1 * rand(1) * ($range[$j][1] - $range[$j][0]);

234
235
236
237

if ($pop[$i][$j] < $range[$j][0]) {
$pop[$i][$j] = $range[$j][0];
}

16

} else {
# Do nothing case
}

238
239
240
241

# Branch migration
if (rand(1) < $CROSSOVER) {
$cur_parent = int rand($KEEP_PARENTS);
}

242
243
244
245

}

246

}

247
248

}

249
250
251

sub print_pop {
my $filename = "out-$file_num.txt";

252

open(FH, "> $filename");

253
254

printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
printf FH ("%13s
print FH "\n";

255
256
257
258
259
260
261
262
263
264
265
266

",
",
",
",
",
",
",
",
",
",
",

"MAX_GEN");
"NUM_ANTS");
"PHER_STEP");
"PHER_WEIGHT");
"MATCH_WEIGHT");
"REGION_WEIGHT");
"INIT_PHER");
"L_DECAY");
"G_DECAY");
"PROB_PROB");
"SCORE");

267
268

for (my $i = 0; $i < $POP_LIMIT; $i++) {
for (my $j = 0; $j <= $SCORE_INDEX; $j++){
printf FH ("%13.9f ", $pop[$i][$j]);
}
print FH "\n";
}
print FH "\n";
close FH;

269
270
271
272
273
274
275
276
277

}

278
279
280
281
282
283
284

####################################
# Put ACO code here
#
#
In the interest of space, the ACO code was not included here.
#
Please see Appendix A to see the implementation of the ACO.
#
The only change that was made was silencing all the output functions.

17

