Message Boards Message Boards

Ant Colony - Speeding up simulation & Extra Functionality

Posted 11 years ago
Following up on Bernat's posting here http://community.wolfram.com/groups/-/m/t/136023?p_p_auth=2Af2XNsw, I attempted my first cut on simulating an ant colony hunting for food.
The code runs too slow, it processes 3.5 cycles (ticks) per second for 25 ants.

Would be very grateful for ideas on how to:
  1. Improve performance of the code.
  2. Mathematica Functions/Ideas for adding diffusion to the pheromone field (represented by a matrix) without adding excessive performance cost (played with converting the matrix to an image > blurring > converting it back to a matrix) which adds and excessive time penalty.
Hope the code is clear enough. As written it takes 517 seconds in my PC.
  (*Define the world size*)
  ymax = 200; xmax = 400;
  (*Locate the colony somwhere in the world*)
  nest = {RandomInteger[{1, xmax}], RandomInteger[{1, ymax}]};
  
  maxAnts = 25;
  ticks = 1;
  maxSteps = 2000;
  maxSpeed = 6;
 maxAngle = 13 \[Pi]/16;
 maxFoodPiles = 8;
 smellRange = 12;
 pheromoneDeposit = 1;
 foodMultiplier = 1000;
 pheromoneMultiplier = 10;
 fieldCaptureInterval = 10;
 discriminationFactor = 16;
 noFood = 0;
 withFood = 1;
 antSymbol = \[Infinity];
 evaporationRate = ConstantArray[(1 - 0.005), {xmax, ymax}];
 pheromoneField = ConstantArray[0., {xmax, ymax}];
 pheromoneField[[First@nest, Last@nest]] =
   100000; (*Added a lot of pheromone at the nest to 'help' the ants \
 get back to the nest*)
 foodField = ConstantArray[0., {xmax, ymax}];
 (*Scatter the food in the world*)
 foodPiles =
   Transpose[{RandomInteger[{1, xmax}, maxFoodPiles],
     RandomInteger[{1, ymax}, maxFoodPiles]}];
 (foodField[[#[[1]], #[[2]]]] = 10) & /@ foodPiles;
 (*We'll record the history of the fields in these variables*)
 pheromoneHistory = {pheromoneField};
 foodHistory = {foodField};
 (*Functions*)
 (*Find out all positions that the ant will step as it moves \
 forward/Keep ants within the world*)
 calculateLoc[location_, direction_, steps_] :=
  Module[{loc, result = {}},
   Do[loc = Flatten[
      location + Round[{temp*Cos[direction], temp*Sin[direction]}]];
    loc = Min[#] & /@ Transpose[{loc, {xmax, ymax}}];
    loc = Max[#] & /@ Transpose[{loc, {1, 1}}];
    result = Append[result, loc], {temp, steps}];
   result]
 (*Return the distribution probability of the direction that the ant \
 will head to*)
 (*It will add the smell of pheromones + food if no food is being \
 carried *)
 probDist[location_, direction_, hasFood_, range_, angle_, steps_] :=
  Module[{dist, total, bins, vals},
   dist = Flatten[
     Table[{direction + y ,
       Flatten[location +
         Round[{x*Cos[direction + y ],
           x*Sin[direction + y ]}]]}, {y, -angle/2, angle/2,
       angle/steps}, {x, 1, range, 1}], 1];
   dist = DeleteCases[dist,
     Alternatives[{_, {a_ /; a < 1, _}}, {_, {a_ /;
         a > xmax, _}}, {_, {_, b_ /; b < 1}}, {_, {_,
        b_ /; b > ymax}}]];
   (dist[[#,
        2]] = ((1 - hasFood) foodMultiplier foodField[[
           First@dist[[#, 2]], Last@dist[[#, 2]]]] +
         pheromoneField[[First@dist[[#, 2]], Last@dist[[#, 2]]]])) & /@
     Range[Length[dist]]; bins = Union[dist[[All, 1]]];
   vals = Cases[dist, {#, y_} -> y] & /@ bins;
   vals = (Total[#] & /@ vals);
   total = Total[vals];
   If[total == 0,
    UniformDistribution[{direction - angle/2, direction + angle/2}],
    EmpiricalDistribution[Rule[vals, bins]]]]
 (*Routine Determines where the ant should head next*)
 (*Also updates the pheromone field*)
 update[agent[id_, symbol_, location_, direction_, hasFood_]] :=
  Module[{temp, dist, dir, speed, loc, path, status, foundDestination,
    multiplier},
   dir = RandomVariate[
     probDist[location, direction, hasFood, smellRange, maxAngle,
      discriminationFactor]];
   speed = RandomVariate[UniformDistribution[{1, maxSpeed}], 1];
   loc = calculateLoc[location, dir, speed];
   status = hasFood;
   If[status == noFood,
       multiplier = 1;
       foundDestination =
        Select[loc, foodField[[#[[1]], #[[2]]]] > 0 &];
       If[Length@foundDestination != 0, status = withFood;
        dir = dir + \[Pi]; multiplier = pheromoneMultiplier;
        foundDestination = First@foundDestination;
        foodField[[foundDestination[[1]], foundDestination[[2]]]] =
         foodField[[foundDestination[[1]], foundDestination[[2]]]] -
          1;],
       multiplier = pheromoneMultiplier;
       foundDestination = Cases[loc, nest];
       If[Length@foundDestination != 0, status = noFood;
        dir = dir + \[Pi]; multiplier = 1;]
       ]
      (pheromoneField[[#[[1]], #[[2]]]] =
       pheromoneField[[#[[1]], #[[2]]]] +
        multiplier pheromoneDeposit) & /@ loc;
 
  agent[id, symbol, Last@loc, dir, status]]

update[agentList_List] :=
Module[{set = update[#] & /@ agentList},
  pheromoneField = evaporationRate pheromoneField;
  If[Mod[ticks, fieldCaptureInterval] == 0,
   pheromoneHistory = Append[pheromoneHistory, pheromoneField];
   foodHistory = Append[foodHistory, foodField];]; ticks++; set]

(*Routines to display the ants in the animate section. Regretably \
can't be used due to memory leak*)
display[agent[id_, symbol_, location_, direction_, hasFood_]] :=
Text[Rotate[symbol, direction], location]
display[list_] := Module[{lst = display[#] & /@ list}, Graphics[lst]]

(*initialize ants*)
(*Initialize the colony - in which direction will the ants start \
walking out*)
ants = Table[
   agent[id, antSymbol, nest,
    RandomVariate[UniformDistribution[{0, 2 \[Pi]}]], noFood], {id,
    maxAnts}];

antHistory = NestList[update, ants, maxSteps]; // AbsoluteTiming

(*Did the ants eat the food? Check status of the food piles and how \
the pheromone field changed with time*)
foodField[[#[[1]], #[[2]]]] & /@ foodPiles
Manipulate[
Show[ColorNegate@ImageRotate@Image@pheromoneHistory[[i]],
  Graphics[{Red, PointSize[Large], Point[nest], PointSize[Medium],
    Blue, Point[foodPiles]}]], {i, 1, Length@pheromoneHistory, 1}]
POSTED BY: Diego Zviovich
3 Replies
There are a number of things I would change in the code before using Compile.
DeleteCases, Union, Append, Total, and the statistics functions
Instead of Map[f, list] I would use Table  or even replace
a = f /@ list

with
Do[a[[i]]=f[[i]],{i,1,Length[list]}]

Instead of DeleteCases, Append, etc., a common trick is to have a dummy value like -1 stand for an unused entry.  Then your array can have a fixed length, and deleting a case would be changing its value to -1.

It's not pretty, but to use Compile it helps to be in the spirit of writing lower level code.

It is possible that shorter Mathematica code would speed things up too.  Maybe find out which function is taking the most time?
POSTED BY: Todd Rowland
Maybe it could be processed by Compile?
POSTED BY: Daniel Lichtblau
Posted 11 years ago
Hi Daniel, thanks for taking the time to write.
cant compile the functions
for the calculateloc Iget the following message Compile::cpts: "The result after evaluating Insert[result,loc,-1] should be a tensor. Nontensor lists are not supported at present; evaluation will proceed with the uncompiled function."
for probDist I get "Compile does not support patterns: \
DeleteCases[dist,{_,{a_/;a<1,_}}|{_,{a_/;a>xmax,_}}|{_,{_,b_/;b<1}}|{_\
,{_,b_/;b>ymax}}] cannot be compiled; evaluation will use the \
uncompiled function."
POSTED BY: Diego Zviovich
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract