The Bee Colony in action

22 May 2011

In our previous installments, we laid the groundwork of our Bee Colony Algorithm implementation. Today, it’s time to put the bees to work, searching for an acceptable solution to the Traveling Salesman problem.

We will approach the search as a Sequence: starting from an initial hive and solution, we will unfold it, updating the state of the hive and the current best solution at each step. Let’s start with the hive initialization. Starting from an initial route, we need to create a pre-defined number of each Bee type, and provide them with an initial destination:

There is probably a more elegant way to do this, but this is good enough: we use a simple List comprehension to generate a list on the fly, yielding the appropriate number of each type of bees, and assigning them a shuffled version of the starting route.

Next, we need a function to update our hive, and the current best solution available:

The function takes 2 arguments: a Tuple of the current state of affairs (the hive and the best solution), and a random number generator. First, we use List.map, to apply the Search function we defined earlier to each bee, returning a new List of Tuples, containing each bee and an option containing the result of its search (a new solution, or None). We then extract the new solutions from that list, using List.choose to retrieve the second element of the Tuples when they are not None, and use a List.fold to find the new best solution from that list, if it is better than the current best solution. Finally, we update the Hive, by having each Bee perform its Waggle dance and share its new information via List.map, and applying the Activate function to promote some of the inactive bees to Active status – and return a Tuple of the updated hive and updated best solution.

We are now almost done – we can now define a Solve function, which will call Initialize to create an initial hive, and unfold an infinite Sequence of solutions:

It’s now time to see this in action. Let’s add a file to our project, and create a small console application Program.fs. First, we need a good test case: let’s generate a list of cities named A to Z, located on a circle:

In about three minutes, starting from a path of length 3450, we found a solution of length 946, out of a possible best of 627. The list has long stretches of correctly reverse-sorted cities (it got D to R properly ordered), with a few misplaced cities. Not too bad!

Before going further, here is the complete code I wrote so far, in its current state. First, the Solver.fs file, which contains the algorithm logic:

moduleSolveropenSystemletprobaFalsePositive=0.1// proba to incorrectly pick a worse solutionletprobaFalseNegative=0.1// proba to miss an improved solutionlettripsLimit=100// number of trips without improvements a bee can makeletprobaConvince=0.8// proba to convince a bee to target a better solutionletSwapIndexPairs(rng:Random)list=seq{foriin(List.lengthlist-1)..-1..1doyield(i,rng.Next(i+1))}letSwapIndexMapindex(moveIndex,toIndex)=ifindex=moveIndexthentoIndexelifindex=toIndexthenmoveIndexelseindexletSwapindexPairlist=letlength=List.lengthlistList.permute(funindex->SwapIndexMapindexindexPair)listletShuffle(rng:Random)list=letlength=List.lengthlistletindexPairs=SwapIndexPairsrnglistSeq.scan(funcurrentListindexPair->SwapindexPaircurrentList)listindexPairs|>Seq.nth(length-1)letSwapWithNextIndexMapindexswapIndexlength=letflipWith=(swapIndex+1)%lengthSwapIndexMapindex(swapIndex,flipWith)letSwapWithNextswapIndexlist=letlength=List.lengthlistList.permute(funindex->SwapWithNextIndexMapindexswapIndexlength)listletSwapRandomNeighborslist=letrandom=newRandom()letindex=random.Next(0,List.lengthlist)SwapWithNextindexlisttypeCity={Name:string;X:float;Y:float;}letDistance(city1,city2)=((city1.X-city2.X)**2.0+(city1.Y-city2.Y)**2.0)**0.5letCircuitCostlist=seq{foriin0..(List.lengthlist-1)->list.[i]yieldList.headlist}|>Seq.pairwise|>Seq.mapDistance|>Seq.sumtypeSolution={Route:List<City>;Cost:float}letEvaluate(route:List<City>)={Route=route;Cost=CircuitCostroute}typeBee=|ScoutofSolution|ActiveofSolution*int|InactiveofSolutionletSearchbee(rng:Random)=matchbeewith|Scoutsolution->letnewSolution=Evaluate(Shufflerngsolution.Route)ifnewSolution.Cost<solution.Costthen(Scout(newSolution),Some(newSolution))else(bee,None)|Active(solution,visits)->letnewSolution=Evaluate(SwapRandomNeighborssolution.Route)letproba=rng.NextDouble()ifnewSolution.Cost<solution.Costthenifproba<probaFalseNegativethen(Active(solution,(visits+1)),None)else(Active(newSolution,0),Some(newSolution))elseifproba<probaFalsePositivethen(Active(newSolution,0),Some(newSolution))else(Active(solution,(visits+1)),None)|Inactivesolution->(bee,None)letSolutionbee=matchbeewith|Scout(solution)->solution|Inactive(solution)->solution|Active(solution,trips)->solutionlet(|RequiresUpdate|)bee=matchbeewith|Scout(solution)->false|Inactive(solution)->true|Active(solution,trips)->trips>tripsLimitletWaggle(solutions:List<Solution>)(bee:Bee)(rng:Random)=matchbeewith|RequiresUpdatetrue->letcurrentSolution=SolutionbeeletnewSolution=List.fold(funaccelement->ifelement.Cost<acc.Cost&&rng.NextDouble()<probaConvincethenelementelseacc)currentSolutionsolutionsInactive(newSolution)|_->beeletPromotebee=matchbeewith|Inactive(solution)->Active(solution,0)|_->beeletActivaternginactivesbees=letinactiveIndexes=List.mapi(funib->matchbwith|Inactive(solution)->Some(i)|_->None)bees|>List.choose(fune->e)|>Shufflerngletpromoted=(List.lengthinactiveIndexes)-inactivesletpromotedIndexes=Seq.takepromotedinactiveIndexes|>Seq.toListbees|>List.mapi(funib->ifList.exists(fune->e=i)promotedIndexesthenPromotebelseb)letCountInactiveshive=hive|>List.choose(funb->matchbwith|Inactive(solution)->Some(b)|_->None)|>List.lengthletInitializenScoutsnActivesnInactivescities(rng:Random)=[foriin1..nScoutsdoletsolution=Evaluate(Shufflerngcities)yieldScout(solution)foriin1..nActivesdoletsolution=Evaluate(Shufflerngcities)yieldActive(solution,0)foriin1..nActivesdoletsolution=Evaluate(Shufflerngcities)yieldInactive(solution)]letUpdate(hive,currentBest:Solution)rng=letsearchResult=List.map(funb->Searchbrng)hiveletnewSolutions=List.choose(fune->snde)searchResultletnewBest=List.fold(funbestsolution->ifbest.Cost<solution.Costthenbestelsesolution)currentBestnewSolutionsletinactives=CountInactiveshiveletupdatedHive=searchResult|>List.map(funb->WagglenewSolutions(fstb)rng)|>Activaternginactives(updatedHive,newBest)letSolvescoutsactivesinactivesroute=letrng=newRandom()lethive=InitializescoutsactivesinactivesrouterngletinitialBest=List.map(funb->Solutionb)hive|>List.minBy(funs->s.Cost)Seq.unfold(funh->Some(h,(Updatehrng)))(hive,initialBest)

At that stage, we have a running algorithm, which seems to be doing what we expect, in about 200 lines of code. I am sure that code could be improved, and I would love to hear comments or suggestions to make it better!

My next objective will be to parallelize that code, to make it hopefully faster. Most of the code should be suitable for this, because we are operating on immutable structures, except for one issue, the random number generator, which uses the non-thread-safe Random() class. Stay tuned!