[Pyxida-cvs-commits] Pyxida/sim Makefile, NONE, 1.1 coord2nn.cc, NONE, 1.1 dist_driver.cc, NONE, 1.1 distributions.cc, NONE, 1.1 distributions.h, NONE, 1.1 energy.cc, NONE, 1.1 error.h, NONE, 1.1 genPerfectCoords.cc, NONE, 1.1 gossip.cc, NONE, 1.1 histogram.cc, NONE, 1.1 histogram.h, NONE, 1.1 math_util.cc, NONE, 1.1 math_util.h, NONE, 1.1 nn.cc, NONE, 1.1 node.cc, NONE, 1.1 node.h, NONE, 1.1 parameters.cc, NONE, 1.1 point.cc, NONE, 1.1 point.h, NONE, 1.1 rs.cc, NONE, 1.1 rs.h, NONE, 1.1 rs_alg.h, NONE, 1.1 rs_alg_greedy.cc, NONE, 1.1 rs_alg_simple.cc, NONE, 1.1 rs_alg_simple.h, NONE, 1.1 rs_aux.cc, NONE, 1.1 rs_stat.cc, NONE, 1.1 samples.cc, NONE, 1.1 samples.h, NONE, 1.1 sim.h, NONE, 1.1 sim_aux.cc, NONE, 1.1 triangle.cc, NONE, 1.1 triangleRPL.cc, NONE, 1.1 vivaldi.cc, NONE, 1.1 vs.cc, NONE, 1.1 zoneTest.cc, NONE, 1.1

[Jonathan L. <le...@us...>]

Update of /cvsroot/pyxida/Pyxida/sim
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv16219

Added Files:
	Makefile coord2nn.cc dist_driver.cc distributions.cc 
	distributions.h energy.cc error.h genPerfectCoords.cc 
	gossip.cc histogram.cc histogram.h math_util.cc math_util.h 
	nn.cc node.cc node.h parameters.cc point.cc point.h rs.cc rs.h 
	rs_alg.h rs_alg_greedy.cc rs_alg_simple.cc rs_alg_simple.h 
	rs_aux.cc rs_stat.cc samples.cc samples.h sim.h sim_aux.cc 
	triangle.cc triangleRPL.cc vivaldi.cc vs.cc zoneTest.cc 
Log Message:
vivaldi and routing simulator code import

--- NEW FILE: coord2nn.cc ---
/* $Id: coord2nn.cc,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2006.  All rights reserved.
 */

#include "node.h"

void printUsage () {
  cout << "coord2nn\n"
       << " -n nodeCount\n"
       << " -c coord file\n"
       << " -D destination ids less than or equal to this id\n"
       << " -x dimensions\n"
       << " -z use height\n"
       << " -k output k nearest neighbors\n"
       << " -e output distance to each neighbor\n";
  exit (-1);
}

int main (int argc, char **argv) {
  char c =  0;
  int MAX_DST = 0;
  char *coordFile = NULL;
  bool outputDistance = false;
  int numNeighbors = 1;
  while ((c = getopt(argc,argv,"x:n:c:zD:k:e")) >= 0) {
    switch (c) {
    case 'x':
      DIMENSIONS = atoi(optarg);
      break;
    case 'k':
      numNeighbors = atoi(optarg);
      break;
    case 'e':
      outputDistance = true;
      break;
    case 'D':
      MAX_DST = atoi(optarg);
      break;
    case 'n':
      nodeCount = atoi(optarg);
      break;
    case 'z':
      USE_HEIGHT = true;
      break;
    case 'c':
      coordFile = optarg;
      break;
    default:
      printUsage ();
      break;
    }
  }

  if (MAX_DST <= 0) {
    printUsage();
  }
  if (coordFile == NULL) {
    printUsage();
  }

  FILE *coordFP = NULL;
  if ((coordFP = fopen (coordFile, "r")) == NULL) {
    printf ("Cannot open file %s", coordFile);
    exit (-1);
  }

  map<int,Point*> node2coord;
  bool endOfFile = false;
  int nodeId;
  int amt;
  while (!endOfFile) {
    amt = fscanf (coordFP, "%d", &nodeId);
    if (amt <= 0) {
      endOfFile = true;
    } else {
      double vec[DIMENSIONS];
      double height;
      for (int d = 0; d < DIMENSIONS; d++) {
	float element;
	int amt = fscanf (coordFP, " %f", &element);
	ASSERT (amt > 0);
	vec[d] = element;
      }
      if (USE_HEIGHT) {
	amt = fscanf (coordFP, " h%f", &height);
	ASSERT (amt > 0);
      }
      amt = fscanf (coordFP, "\n");
      Point *p = new Point (vec, height);
      node2coord.insert(pair<int,Point*>(nodeId,p));
    }
  }

  for (map<int,Point*>::const_iterator dstIter = node2coord.begin();
       dstIter != node2coord.end(); dstIter++) {
    int dstId = dstIter->first;
    Point *dstPoint = dstIter->second;
    if (dstId > MAX_DST) {
      map<double,int> dist2id;
      for (map<int,Point*>::const_iterator srcIter = node2coord.begin();
       srcIter != node2coord.end(); srcIter++) {
	int srcId = srcIter->first;
	Point *srcPoint = srcIter->second;
	if (srcId <= MAX_DST) {
	  double distance = srcPoint->getEuclideanDistance(dstPoint);
	  dist2id.insert(pair<double,int>(distance,srcId));
	}
      }
      int nearCount = 0;
      printf ("%d", dstId);
      for (map<double,int>::const_iterator nearIter = dist2id.begin();
	   nearCount < numNeighbors &&
	   nearIter != dist2id.end(); nearIter++) {
	double nearestDistance = nearIter->first;
	int nearestId = nearIter->second;
	printf (" %d", nearestId);
	if (outputDistance) {
	  printf (" %f", nearestDistance);
	}
	nearCount++;
      }
      printf ("\n");
    }
  }
}
  
    

--- NEW FILE: dist_driver.cc ---
/* $Id: dist_driver.cc,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2005.  All rights reserved.
 */

#include <stdio.h>
#include <stdlib.h>
#include "error.h"
#include "distributions.h"

int main () {
  SRand(getpid());
  ZipfIntegerDistribution *z = new ZipfIntegerDistribution (2.0, 10);
  for (int i = 0; i<10; i++) {
    printf ("%d %d\n", i, z->next());
  }
}

--- NEW FILE: rs.cc ---
/* $Id: rs.cc,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 */

#include "distributions.h"
#include "rs.h"
#include "rs_alg.h"

int ringBase = 2;
Node *nodes;
int summaryFrequency = 1000;
FILE *runLogFP, *summaryFP, *traceFP;
FILE *hopFP;
bool useChurn = false;
bool useGossipSet = false;
bool recordDelayStretch = false;
bool useGreedyByZone = true;
bool usePlanarDistance = false;
bool useNodeWithNearestDistance = false;
void performNNcomparison ();

/*
 * The gist:
 * Take as input: a coordinates file, a latency matrix,
 *  and a routing algorithm.
 * The nodes get divided into two types: participants and destinations.
 * The simulator runs for a given number of rounds.
 * In every round (which is meant to represent a second)
 *   each node is either born or dies or stays alive or dead,
 *   a given number of queries take place
 *     queries originate from a random node.
 *     targets are either randomly or Zipf distributed
 *   each node can attempt to update its routing tables
 *    by sending one or more messages
 */

void printUsage (char* problem) {
  cout 
    << "Usage: rs\n"
    <<"\n"
    <<"\n           Coordinate parameters\n"
    <<"            -m routing strategy\n"
    <<"               g: simple greedy, 5 nearby neighbors, 0 long distance\n"
    <<"               s: Small World, 4 nearby neighbors, 1 random neighbor\n"
    <<"               t: Theta Graphs: angle=pi/16, levels=1\n"
    <<"               h: Hassan and Peleg's Scaled-Theta Graphs: angle=pi/16, levels=4\n"
    <<"               c: Ratnasamy's CAN: angle=pi/4, level=1\n"
    <<"               a: Abraham and Malkhi's Compact Routing\n"
    <<"            -s num sectors (-s and -l override parameters from -m)\n"
    <<"            -l num levels\n"
    <<"            -b base for ring radius (default "<<ringBase<<")\n"
    <<"            -o routing options\n"
    <<"               d: do routing greedily by distance, not by zone\n"
    <<"               p: use planar distance to determine zones, not euclidean\n"
    <<"               c: use nearest coordinate as nearest point, not nearest latency on routing path\n"
    <<"\n           Coordinate parameters\n"
    <<"            -c initial coordinates filename\n"
    <<"            -x coordinate dimensions\n"
    <<"            -z coordinates use height\n"
    <<"\n           Workload distribution parameters\n"
    <<"            -w queries per round\n"
    <<"            -K routing destinations (uniform or zipf), e.g. U/80/20\n"
    <<"               U/80/20 uniformly at random from \n"
    <<"\n           Node/Churn parameters\n"
    <<"            -r RTT filename (Euclidean distance btw nodes used if missing)\n"
    << "           -C use churn with this as mean lifetime, poisson distributed\n"
    <<"            -g gossip strategy (r and n are mutually exclusive)\n"
    <<"               r: random gossip\n"
    <<"               n: nearby neighbor gossip\n"
    <<"               j: seed tables by routing to our coord from random node\n"
    <<"               a: gossip everyone we know, not just tables\n"
    <<"            -n overlay node count\n"
    <<"            -t target nodes count\n"
    <<"               n+t = nodes in RTT file\n"
    <<"\n"
    <<"           Simulator parameters\n"
    <<"            -V types of output\n"
    <<"               r: log of each query (.log) \n"
    <<"               h: track hop distances (.hop) and delay stretch\n"
    <<"            -S random seed\n"
    <<"            -R rounds\n"
    <<"            -P start recording stats after this round (default: no delay)\n"
    <<"            -D debug\n"
    <<"            -F frequency of runtime summary output (once every "<<summaryFrequency<<" rounds is default)\n"
    <<"            -O output file prefix\n"
    <<"            -N compare coord to rtt nearest neighbours\n"
    <<"\n"
    <<"Distributions:\n"
    <<" Pareto       P/scale/shape\n"
    <<" Normal       N/mean/stddev\n"
    <<" Zipf         Z/alpha/num-elements\n"
    <<" Poisson      F/mean\n"
    <<" Uniform      U [0..1)\n"
    <<" Constant     C/value\n";

  if (problem) {
    cout << problem << endl;
  }
  exit (-1);
}

int main (int argc, char **argv) {
  Strategy strategy;
  NNAlg *nnAlg = NULL;
    
  char routingStrategyChar = ' ';
  char *churnDistStr = NULL, *keyDistStr = NULL, *outputPrefix = NULL,
    *rttFilename = NULL, *coordsFilename;
  IntegerDistribution *keyDist;
  Distribution *churnDist;
  DistributionFactory *distFactory = NULL;
  char c = 0;
  int seed = getpid();
  int rounds = 0;
  int meanLifetime = 0;
  bool gossipOnJoin = false;
#ifdef DEBUG
  debug = false;
#endif
  double queriesPerRound = -1.;
  bool finished = false;
  char *routingStrategyOptionsStr = NULL;
  char *gossipStrategyStr = NULL;
  char *outputTypeStr = NULL;
  int numSectors = 0;
  int numLevels = 0;
  bool perQueryOutput = false;
  int recordingStatsRound = 0;
  char *nn_alg_options = NULL;
  double lKatJoin = 0.;
  bool nnComparison = false;

  distFactory = new DistributionFactory ();
  
  while ( (c = getopt(argc,argv,"m:c:C:w:K:S:Df:O:n:u:r:x:t:F:R:o:g:s:l:V:b:P:zo:N")) >= 0) {
    switch (c) {
    case 'N':
      nnComparison = true;
      break;
    case 'g':
      gossipStrategyStr = optarg;
      break;
    case 'm':
      routingStrategyChar = *optarg;
      break;
    case 'V':
      outputTypeStr = optarg;
      break;
    case 'o':
      nn_alg_options = optarg;
      break;
    case 'b':
      ringBase = atoi(optarg);
      break;
    case 'P':
      recordingStatsRound = atoi(optarg);
      break;
    case 's':
      numSectors = atoi(optarg);
      break;
    case 'l':
      numLevels = atoi(optarg);
      break;
    case 'C':
      meanLifetime = atoi(optarg);
      useChurn = true;
      break;
    case 'K':
      keyDistStr = optarg;
      break;
    case 'w':
      queriesPerRound = atof (optarg);
      break;
    case 'r':
      rttFilename = optarg;
      break;
    case 'c':
      coordsFilename = optarg;
      break;
    case 'z':
      USE_HEIGHT = true;
      break;
    case 'x':
      DIMENSIONS = atoi (optarg);
      //printf ("set dim to %d\n", DIMENSIONS);
      break;
    case 'R':
      rounds = atoi (optarg);
      break;
    case 'n':
      nodeCount = atoi (optarg);
      break;
    case 't':
      dstNodeCount = atoi (optarg);
      break;
    case 'S':
      seed = atoi (optarg);
      break;
    case 'F':
      summaryFrequency = atoi (optarg);
      break;
#ifdef DEBUG
    case 'D':
      debug = true;
      break;
#endif
    case 'O':
      outputPrefix = optarg;
      break;
    default:
      printUsage ("Undefined option");
      break;
    }
  }

  systemNodeCount = nodeCount + dstNodeCount;
  SRand (seed);


  if (ringBase <= 1) 
    printUsage("Ring base must be greater than 1");
  if (recordingStatsRound < 0 || recordingStatsRound >= rounds) 
    printUsage("Invalid recording stats round");
  if (coordsFilename == NULL)
    printUsage("Error with coordinate filename");
  if (keyDistStr == NULL)
    printUsage("Error in key distribution");
  if (summaryFrequency <= 0)
    printUsage("f must be greater than 0");
  if (queriesPerRound < 0.)
    printUsage("Workload rate must be greater than or equal to 0.");
  if (outputPrefix == NULL)
    printUsage("Output file prefix is required.");
  if (DIMENSIONS <= 0)
    printUsage("Dimensions must be greater than 0");
  if (!useChurn && gossipStrategyStr != NULL) {
    printUsage("Cannot have a gossip strategy without churn");
  }
  if (gossipStrategyStr != NULL) {
    if ((strchr(gossipStrategyStr,'r') != NULL) &&
	(strchr(gossipStrategyStr,'n') != NULL)) {
      printUsage("Gossip strategies are mutually exclusive");
    }

    if (strchr(gossipStrategyStr,'r')) {
      gossipStrategy = Random;
    }
    if (strchr(gossipStrategyStr,'n')) {
      gossipStrategy = Local;
    }
    if (strchr(gossipStrategyStr,'j')) {
      gossipOnJoin = true;
    }
    if (strchr(gossipStrategyStr,'a')) {
      useGossipSet = true;
    }
  }

  if (useChurn && gossipStrategy == None) {
    printUsage("Must have a gossip strategy with churn");
  }

  if (outputTypeStr != NULL) {
    if (strchr(outputTypeStr,'r')) {
      perQueryOutput = true;
    }
    if (strchr(outputTypeStr,'h')) {
      recordDelayStretch = true;
    }
  }

  if (nn_alg_options != NULL) {
    if (strchr(nn_alg_options,'d')) {
      useGreedyByZone = false;
    }
    if (strchr(nn_alg_options,'p')) {
      usePlanarDistance = true;
    }
    if (strchr(nn_alg_options,'c')) {
      useNodeWithNearestDistance = true;
    }
  }

  /*
   * INITIALIZE SYSTEM
   */

  if (recordingStatsRound == 0) {
    recordingStats = true;
  }

  fprintf (stderr, "starting initialization\n");
  if (debug) printf ("pickDstNodes\n");
  pickDstNodes ();
  if (debug) printf ("readRTTmatrix\n");
  readRTTmatrix (rttFilename);
  if (debug) printf ("initNodes\n");
  initNodes (coordsFilename);
  if (debug) printf ("assignCoords\n");
  assignCoords(nodes);
  if (debug) printf ("initRankMap\n");
  initRankMap ();

  if (nnComparison) {
    performNNcomparison ();
    exit (0);
  } else {
    if (!verifyRTTmatrix ()) {
      fprintf (stderr, "RTT matrix is not complete");
      exit (-1);
    }
  }

  switch (routingStrategyChar) {
  case 'g':
    strategy = SimpleGreedy;
    printUsage ("not implemented");
    nnAlg = new NNSmallWorld(5, 0);
    break;
  case 's':
    strategy = SmallWorld;
    printUsage ("not implemented");
    nnAlg = new NNSmallWorld(4, 1);
    break;
  case 't':
    strategy = Theta;
    if (numLevels == 0) numLevels = 8;
    if (numSectors == 0) numSectors = 1;
    // This is what theta means
    ASSERT (numLevels == 1);
    nnAlg = new NNTheta(numSectors, numLevels, 1);
    break;
  case 'h':
    strategy = ScaledTheta;
    if (numLevels == 0) numLevels = 8;
    if (numSectors == 0) numSectors = 4;
    nnAlg = new NNTheta(numSectors, numLevels, ringBase);
    break;
  case 'c':
    strategy = CAN;
    printUsage ("not implemented");
    nnAlg = new NNTheta(4, 1, 1);
    break;
    /*
      case 'a':
      strategy = Compact;
      nnAlg = new NNCompact();
      break;
    */
  default:
    printUsage ("Unknown routing strategy");
    break;
  }


  // Initialise the routing tables according to the algorithm used
  if (!useChurn) {
    if (debug) printf ("initPerfectRoutingTables\n");
    bool wasDebug = debug;
    debug = false;
    for (int i = 0; i < nodeCount; i++) {
      nnAlg->initPerfectRoutingTables(nodes[i]);
    }
    debug = wasDebug;
    vector<int> rt;
    //nnAlg->printState (nodes[0], rt);
  }

  nodeUp = nodeCount;

  vector<int> allNodes;
  for (int i = 0; i < nodeCount; i++) {
    allNodes.push_back(i);
  }

  if (useChurn) {

    // stir and shake
    for (int i = 0; i < nodeCount; i++) {
      int targetLifetime = 10*meanLifetime;
      int myRound = 0;
      bool alive = true;
      while (myRound < targetLifetime) {
	alive = alive ? false : true;
	myRound += poisson(meanLifetime);
      }
      if (alive) {
	nodes[i].deathTime = myRound - targetLifetime;
	nodes[i].alive = true;
        if (debug) printf ("start alive %d DT %d\n", i, nodes[i].deathTime);
      } else {
	nodes[i].birthTime = myRound - targetLifetime;
	nodes[i].alive = false;
        if (debug) printf ("start dead %d BT %d\n", i, nodes[i].birthTime);
	nodeUp--;
      }
    }

    vector<int> aliveNodes;
    vector<int> aliveNodesShuffle;
    for (int i = 0; i < nodeCount; i++) {
      if (nodes[i].alive) {
	aliveNodes.push_back(i);
	aliveNodesShuffle.push_back(i);
      }
    }

    const int initialNodeKnowledge = 5;
    ASSERT (aliveNodes.size() >= initialNodeKnowledge);

    for (int i = 0; i < aliveNodes.size(); i++) {
      random_shuffle (aliveNodesShuffle.begin(),aliveNodesShuffle.end());
      for (int n = 0; n < initialNodeKnowledge; n++) {
	nnAlg->addNeighbour(nodes[aliveNodes[i]],nodes[aliveNodesShuffle[n]]);
      }

    }

  }

  Statistics *sumStats = new Statistics ();
  Statistics *interimStats = new Statistics ();

  keyDist = distFactory->newIntegerDistribution (keyDistStr);
  if (keyDist == NULL)
    printUsage ("Problem with key distribution");

  cRound = 0;

  // very detailed per query stats go here
  openOutputFile (runLogFP, outputPrefix, ".run");

  // other logging info goes here
  openOutputFile (traceFP, outputPrefix, ".trace");

  // hop distances
  if (recordDelayStretch) {
    openOutputFile (hopFP, outputPrefix, ".hop");
  }

  printf ("end of initialization\n");

  for (cRound = 0; cRound < rounds; ) {

    // Start recording stats
    if (!recordingStats && cRound == recordingStatsRound) {
      recordingStats = true;
      // Don't need to do this because we aren't adding to them below
      //sumStats->clear();
      fprintf (traceFP, "R %d Recording stats.\n", cRound);
    }

    DEBUG_P (("ROUND %d\n", cRound));

    /*
     * PER ROUND BIRTH AND DEATH
     */

    if (useChurn) {
      random_shuffle (allNodes.begin(),allNodes.end());

      for (int a = 0; a < nodeCount; a++) {
	int i = allNodes[a];
	if (nodes[i].birthTime == cRound) {
          ASSERT (!nodes[i].alive);
	  nodeUp++;
	  nodes[i].alive = true;
	  nodes[i].deathTime = cRound + poisson(meanLifetime) + 1;
	  nnAlg->clearRoutingTables(nodes[i]);
	  if (debug) 
	    cout << cRound << " node "<< i << " born, dT " << nodes[i].deathTime << endl;


	  // join procedure for bootstrapping routing tables
	  if (gossipOnJoin) {

	    // pick random alive node
	    int bootstrapNodeId = i;
	    while (bootstrapNodeId == i ||
		   !nodes[bootstrapNodeId].alive) {
	      bootstrapNodeId = unifRand (0, nodeCount);
	    }

	    // do query for self
	    int bootResultNodeId;
	    deque<int> bootNodeIds;
	    double bootDelayStretch = 0.;
            vector<double> bootHopDistance;

	    //debug = true;
            if (debug)
              printf ("R %d node %d starting query for self from node %d\n",
                      cRound, i, bootstrapNodeId);

	    //nnAlg->printState();

	    nnAlg->route(nodes[bootstrapNodeId], i, bootResultNodeId,
			 bootNodeIds, bootHopDistance, bootDelayStretch);
	    // route() pops this off
	    bootNodeIds.push_front(bootstrapNodeId);

	    // add everyone seen on the path to our RT
	    for (deque<int>::const_iterator bootIter = bootNodeIds.begin();
		 bootIter != bootNodeIds.end(); bootIter++) {
	      if (debug)
                printf ("R %d node %d learned about node %d from query for self\n",
                        cRound, i, *bootIter);
	      nnAlg->gossipPair(nodes[i],nodes[*bootIter]);

	      if (debug) {
		vector<int> rt;
		nnAlg->printState(nodes[i],rt);
	      }
	      //debug = false;
	    }
	  } else {
            // gossip once with somebody to bootstrap our RT
            nnAlg->gossip(nodes[i]);
          }

	  double gossipSize, angleTableSize, pctNearbyKnown;
	  nnAlg->state (nodes[i], gossipSize, angleTableSize, pctNearbyKnown);
	  if (debug)
	    printf ("join %d gS %.f aT %.f lK %.2f\n",
		    i, gossipSize, angleTableSize, pctNearbyKnown);
	  const double ALPHA = 0.05;
	  lKatJoin = (ALPHA*pctNearbyKnown)+((1.-ALPHA)*lKatJoin);


	} else if (nodes[i].deathTime == cRound) {
          ASSERT (nodes[i].alive);
	  nodeUp--;
	  nodes[i].alive = false;
	  nodes[i].birthTime = cRound + poisson(meanLifetime) + 1;
	  nnAlg->clearRoutingTables(nodes[i]);
	  if (debug)
	    cout << cRound << " node "<< i << " died, dT " << nodes[i].birthTime << endl;


	}
      }
    }

    if (debug) printf ("upFrac %.3f\n", nodeUp/(double)nodeCount);
    ASSERT (nodeUp > 1);

    /*
     * PER ROUND GOSSIP
     */


    if (useChurn) {
      DEBUG_P (("\nGOSSIP: %d\n", cRound));

      // shuffle order of nodes gossipping
      vector<int> aliveNodes;
      for (int i = 0; i < nodeCount; i++) {
	if (nodes[i].isAlive() && nodes[i].takeTurn()) {
          aliveNodes.push_back(i);
	}
      }

      random_shuffle (aliveNodes.begin(),aliveNodes.end());
      for (int i = 0; i < aliveNodes.size(); i++) {
	//if (aliveNodes[i] == 0) debug = true;
        nnAlg->gossip(nodes[aliveNodes[i]]);
	//debug = false;
      }

    }

    /*
     * PER ROUND QUERIES
     */

    if (recordingStats) {
    DEBUG_P (("QUERIES\n"));
    int queryCount = (int)(ceil(nodeUp * queriesPerRound));
    if (debug) {
      printf ("nodeUp %d: ", nodeUp);
      for (int i = 0; i < nodeCount; i++) {
        if (nodes[i].alive) {
          printf ("%d ",i);
        }
      }
      printf ("\n");
    }
    // figure out the source of this round's queries
    vector<int> srcNodeIds;
    srcNodeIds.reserve (queryCount);

    vector<int> dstNodeIds;
    dstNodeIds.reserve (queryCount);

    for (int i = 0; i < queryCount; i++) {
      int dstNodeId = -1;
      while (dstNodeId < 0 || !nodes[dstNodeId].alive)
	if (dstNodeCount == 0) {
	  dstNodeId = keyDist->next();
	} else {
	  dstNodeId = keyDist->next() + nodeCount;
	}

      int srcNodeId = dstNodeId;
      while (srcNodeId == dstNodeId || !nodes[srcNodeId].alive)
	srcNodeId = unifRand (0, nodeCount);

      // make num of target nodes zero + querysize = node count -K U/100
      // -n 100 -t 0 -S 1 -R 10 -O foo -K U/100 -m h
      ASSERT (srcNodeId >= 0 && srcNodeId < nodeCount);
      srcNodeIds.push_back (srcNodeId);
      dstNodeIds.push_back (dstNodeId);
    }

    ASSERT (dstNodeIds.size() == queryCount);
    ASSERT (srcNodeIds.size() == queryCount);

    int queryIndex = -1;
    int querySuccessCount = 0;

    for (int queryIndex = 0; queryIndex < queryCount; queryIndex++) {
      int dstNodeId = dstNodeIds[queryIndex];
      int srcNodeId = srcNodeIds[queryIndex];
      ASSERT (dstNodeId >= 0 && dstNodeId < systemNodeCount);
      ASSERT (srcNodeId >= 0 && srcNodeId < nodeCount);
      DEBUG_P (("Starting query from %d->%d\n", srcNodeId, dstNodeId));

      int resultNodeId;
      deque<int> hopNodeIds;
      double delayStretch = 0.;
      vector<double> hopDistance;
      //if (cRound == 500) debug = true;

      nnAlg->route(nodes[srcNodeId], dstNodeId, resultNodeId, hopNodeIds, hopDistance, delayStretch);
      //if (cRound == 500) exit (0);

      int nearestNodeId;
      int ncRank, latRank;
      double glRelLatPenalty, glAbsLatPenalty, ncRelLatPenalty, ncAbsLatPenalty;
      nnQueryOptimal (dstNodeId, resultNodeId,
		      nearestNodeId,
                      latRank,
		      glRelLatPenalty, glAbsLatPenalty,
		      ncRank, 
		      ncRelLatPenalty, ncAbsLatPenalty);

      // with euclidean distances we should always find the nearest
      // er..., maybe not

      if (ncRank != 0 && !useChurn && rttFilename == NULL) {
      //if (glRelLatPenalty > 26) {
        //(srcNodeId == 89 && dstNodeId == 91)){
        debug = true;
        nnAlg->route(nodes[srcNodeId], dstNodeId, resultNodeId, hopNodeIds, hopDistance, delayStretch);
        nnQueryOptimal (dstNodeId, resultNodeId, nearestNodeId, latRank, glRelLatPenalty, glAbsLatPenalty, ncRank, ncRelLatPenalty, ncAbsLatPenalty);
        debug = false;
        printf ("quitting because rank was %d src %d dst %d res %d", 
                ncRank, srcNodeId, dstNodeId, resultNodeId);
        exit (0);
      }

      interimStats->nnQuery (ncRank, glRelLatPenalty,
			     glAbsLatPenalty, ncRelLatPenalty,
			     ncAbsLatPenalty,
			     hopNodeIds.size(), hopDistance, delayStretch);
      if (perQueryOutput) {
      fprintf (runLogFP, "src %d dst %d res %d near %d rank %d grlp %.2f galp %.2f nrlp %.2f nalp %.2f hS %d delay %.3f\n", 
	       srcNodeId, dstNodeId, resultNodeId, nearestNodeId, ncRank, 
	       glRelLatPenalty, glAbsLatPenalty, 
	       ncRelLatPenalty, ncAbsLatPenalty, hopNodeIds.size(), delayStretch);
      }
      if (debug) {
	printf ("src %d dst %d res %d near %d rank %d grlp %.2f galp %.2f nrlp %.2f nalp %.2f hS %d delay %.3f\n", 
	       srcNodeId, dstNodeId, resultNodeId, nearestNodeId, ncRank, 
	       glRelLatPenalty, glAbsLatPenalty, 
	       ncRelLatPenalty, ncAbsLatPenalty, hopNodeIds.size(), delayStretch);
      }


    }
    }


    /*
     * RECORD END OF ROUND
     */

    interimStats->nodesUp (nodeUp);
    cRound++;

    if (cRound % summaryFrequency == 0) {
      nnAlg->printState ();
      double gossipSize, angleTableSize, pctNearbyKnown;
      nnAlg->state(gossipSize, angleTableSize, pctNearbyKnown);
      interimStats->gossipState(gossipSize, angleTableSize, pctNearbyKnown);
      if (recordingStats)
	sumStats->add (interimStats);
      fprintf (runLogFP, "R %d ", cRound);
      fprintf (stdout, "R %d ", cRound);
      fprintf (stdout, "R %d lKatJoin %.2f\n", cRound, lKatJoin);
      fprintf (runLogFP, "R %d lKatJoin %.2f\n", cRound, lKatJoin);
      interimStats->print (runLogFP);
      interimStats->print (stdout);
      interimStats->clear ();
    }

  } // END OF ROUNDS

  /*
   * PRINT STATS
   */

  double gossipSize, angleTableSize, pctNearbyKnown;
  nnAlg->state(gossipSize, angleTableSize, pctNearbyKnown);
  interimStats->gossipState(gossipSize, angleTableSize, pctNearbyKnown);

  sumStats->add (interimStats);
  openOutputFile (summaryFP, outputPrefix, ".sum");

  // print maxNeighbors
  /*
    fprintf (summaryFP, "dst_dist %s q_per_r %d rtts %s dim %d node_count %d dst_count %d seed %d prefix %s\n",
    keyDistStr, queriesPerRound, rttFilename,
    DIMENSIONS, nodeCount, dstNodeCount, seed, outputPrefix);
  */
	   
  sumStats->print (summaryFP);
  sumStats->print (stdout);
  if (recordDelayStretch)
    sumStats->printHops (hopFP);

  // destroy nodes and route tables

  printf ("Shutting down.\n");

  delete nnAlg;
  deleteNodes ();
  delete keyDist;
  delete interimStats;
  delete sumStats;
  delete distFactory;
}

void nnQueryOptimal (int targetNodeId, int resultNodeId, 
		     int &nearestNodeId,
                     int &latRank,
		     double &glRelLatPenalty, 
		     double &glAbsLatPenalty,
		     int &ncRank,
		     double &ncRelLatPenalty, 
		     double &ncAbsLatPenalty) {

  ASSERT (resultNodeId >= 0 && resultNodeId < nodeCount);

  map<int,vector<int> >::const_iterator targetLatIter =
    targetId2rankLat2nodeId.find (targetNodeId);
  //ASSERT (targetIter != targetId2rankLat2nodeId.end());
  //if (debug) printf ("nnQueryOptimal target %d\n", targetNodeId);
  if (targetLatIter == targetId2rankLat2nodeId.end()) {
    ASSERT (targetLatIter != targetId2rankLat2nodeId.end());
  }

  map<int,vector<int> >::const_iterator targetDistIter =
    targetId2rankDist2nodeId.find (targetNodeId);
  if (targetDistIter == targetId2rankDist2nodeId.end()) {
    ASSERT (targetDistIter != targetId2rankDist2nodeId.end());
  }


  vector<int> rankLat2nodeId = targetLatIter->second;
  vector<int> rankDist2nodeId = targetDistIter->second;

  latRank = -1;
  for (int i = 0; i < rankLat2nodeId.size(); i++) {
    if (rankLat2nodeId[i] == resultNodeId) {
      latRank = i;
      break;
    }
  }
  if (latRank < 0) {
    printf ("target %d res %d size %d\n", targetNodeId,
	    resultNodeId, rankLat2nodeId.size());
    ASSERT (latRank >= 0);
  }


  ncRank = -1;
  for (int i = 0; i < rankDist2nodeId.size(); i++) {
    if (rankDist2nodeId[i] == resultNodeId) {
      ncRank = i;
      break;
    }
  }
  if (ncRank < 0) {
    printf ("target %d res %d size %d\n", targetNodeId,
	    resultNodeId, rankDist2nodeId.size());
    ASSERT (ncRank >= 0);
  }


  // find the nearest alive node to the target
  int nearestLatNodeId = -1;
  for (int i = 0; i < rankLat2nodeId.size(); i++) {
    if (nodes[rankLat2nodeId[i]].alive) {
      nearestLatNodeId = rankLat2nodeId[i];
      break;
    } else {
      if (debug) printf ("skipping node w rankLat %d\n", i);
    }
  }
  ASSERT (nearestLatNodeId >= 0 && nearestLatNodeId < nodeCount);

  int nearestDistNodeId = -1;
  for (int i = 0; i < rankDist2nodeId.size(); i++) {
    if (nodes[rankDist2nodeId[i]].alive) {
      nearestDistNodeId = rankDist2nodeId[i];
      break;
    } else {
      if (debug) printf ("skipping node w rankDist %d\n", i);
    }
  }
  ASSERT (nearestDistNodeId >= 0 && nearestDistNodeId < nodeCount);

  double resultLat = getLatency(targetNodeId,resultNodeId);
  double nearestLatLat = getLatency(targetNodeId,nearestLatNodeId);
  // dstNodeCount == 0 we are not doing NN
  ASSERT ((dstNodeCount == 0 && resultLat >= 0) || (resultLat > 0));
  ASSERT ((dstNodeCount == 0 && nearestLatLat >= 0) || (nearestLatLat > 0));
  if (nearestLatLat > resultLat) {
    ASSERT (nearestLatLat <= resultLat);
  }

  double nearestDistLat = getLatency(targetNodeId,nearestDistNodeId);
  ASSERT ((dstNodeCount == 0 && nearestDistLat >= 0) || (nearestDistLat > 0));
  if (nearestLatLat > nearestDistLat) {
    ASSERT (nearestLatLat <= nearestDistLat);
  }

  

  glAbsLatPenalty = resultLat - nearestLatLat;
  glRelLatPenalty = 0;
  if (nearestLatLat > 0)
    glRelLatPenalty = glAbsLatPenalty / nearestLatLat;

  ncAbsLatPenalty = resultLat - nearestDistLat;
  ncRelLatPenalty = 0;
  if (nearestDistLat > 0)
    ncRelLatPenalty = ncAbsLatPenalty / nearestDistLat;

  if (debug) {
    cout 
      << "target "<< targetNodeId
      << " nearestLat "<< nearestLatNodeId
         << " " << nodes[nearestLatNodeId].vec
      << " distance=" << getDistance (nodes[nearestLatNodeId],nodes[targetNodeId])
         << " latency= " << nearestLatLat
      << " nearestDist "<< nearestDistNodeId
         << " " << nodes[nearestDistNodeId].vec
      << " distance=" << getDistance (nodes[nearestDistNodeId],nodes[targetNodeId])
         << " latency= " << nearestDistLat
      << " result "<< resultNodeId
         << " " << nodes[resultNodeId].vec
         << " distance=" << getDistance(nodes[resultNodeId],nodes[targetNodeId])
         << " latency= " << resultLat << endl;
  }
}

void performNNcomparison () {
  
  // for each node, walk through the list of coordinates in order of distance
  // and output:
  // - absolute latency penalty from the top 1,3,5,7,15 coordinates

  if (dstNodeCount != 0) {
    fprintf (stderr, "use -T 0 for full matrix\n");
    exit (-1);
  }

  vector<int> range;
  for (int p = 0; p <= 4; p++) {
    int rMax = (int)(pow(2,p));
    range.push_back(rMax); 
  }

  for (map<int,vector<int> >::const_iterator targetDistIter =
	 targetId2rankDist2nodeId.begin ();
       targetDistIter != targetId2rankDist2nodeId.end (); targetDistIter++) {

    // who are we considering
    int targetNodeId = targetDistIter->first;

    // what is the latency of his nearest latency neighbour
    map<int,vector<int> >::const_iterator targetLatIter =
      targetId2rankLat2nodeId.find (targetNodeId);
    ASSERT (targetLatIter != targetId2rankLat2nodeId.end());
    vector<int> rankLat2nodeId = targetLatIter->second;
    if (rankLat2nodeId[0] != targetNodeId) {
      printf ("target %d rank0 %d rank1 %d rank2 %df\n",
	      targetNodeId, rankLat2nodeId[0],
	      rankLat2nodeId[1], rankLat2nodeId[2]);
      ASSERT (rankLat2nodeId[0] == targetNodeId);
    }
    int nearestLatNodeId = rankLat2nodeId[1];
    ASSERT (nearestLatNodeId >= 0 && nearestLatNodeId < nodeCount);
    double nearestLatency = getLatency(nearestLatNodeId,targetNodeId);
    if (nearestLatency <= 0) {
      printf ("target %d nearest %d nearestLat %.3f\n",
	      targetNodeId, nearestLatNodeId, nearestLatency);
      ASSERT (nearestLatency > 0);
    }

    // Repeatedly walk the list of nodes that are sorted
    // by coordinate distance to the target.
    // Pick the one whose latency to the target is lowest
    // and compare that with the true lowest latency.
    vector<int> rankDist2nodeId = targetDistIter->second;

    printf ("%d", targetNodeId);

    for (int r = 0; r < range.size(); r++) {
      double minLatency = -1;
      int rStep = 0;
      // skip the first one, its the guy itself
      vector<int>::const_iterator nodeIter = rankDist2nodeId.begin();
      nodeIter++;
      for (; rStep < range[r] && nodeIter != rankDist2nodeId.end(); nodeIter++) {
	int nodeId = *nodeIter;

	ASSERT (nodeId != targetNodeId);

	if (nodeId < 0 || nodeId >= nodeCount) {
	  printf ("r %d rStep %d target %d nodeId %d\n",
		  r, rStep, targetNodeId, nodeId);
	  ASSERT (nodeId >= 0 && nodeId < nodeCount);
	}
	ASSERT (nodeId >= 0 && nodeId < nodeCount);
	
	if (debug)
	  printf ("r %d rStep %d nodeId %d dist %.2f lat %.2f target %d\n", r, rStep, nodeId,
		nodes[nodeId].vec->getEuclideanDistance(nodes[targetNodeId].vec),
		  getLatency (nodeId, targetNodeId), targetNodeId);

	//if (nodeId == targetNodeId) {
	//exit (-1);
	  ASSERT (nodeId != targetNodeId);
	//}
	double nodeLatency = getLatency (nodeId, targetNodeId);
	ASSERT (nodeLatency != 0);

	// It can happen that we have a close distance to a guy
	// but his distance isn't in the latency matrix

	if (nodeLatency > 0) {

	  if (minLatency == -1 || nodeLatency < minLatency) {
	    minLatency = nodeLatency;
	  }
	  rStep++;
	}
      }
      ASSERT (minLatency > 0);
      ASSERT (minLatency >= nearestLatency);
      double absLatencyPenalty = minLatency - nearestLatency;
      printf (" r%d %.2f", range[r], absLatencyPenalty);
    }
    printf ("\n");

  }
}

--- NEW FILE: error.h ---
/* $Id: error.h,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2005.  All rights reserved.
 */

#ifndef ERROR_H
#define ERROR_H

#include <stdio.h>
#include <stdlib.h>

// Comment out to remove all debugging paths
#define DEBUG 1

#ifdef DEBUG
extern bool debug;
#define DEBUG_P(c) if (debug) printf c;
#define IF_DEBUG(c)   if (debug) c;
#else
#define DEBUG_P(c) 
#define IF_DEBUG(c)
#endif

#define ASSERT(c) { \
  if (!(c)) { \
    fprintf (stderr, \
	     "ASSERTION FAILURE: %s:%d\n", __FILE__, __LINE__); \
    fprintf (stderr, "\tfailed condition: %s\n", #c); \
    exit (-1); \
  } \
}

//#define ASSERT(c) 

#endif

--- NEW FILE: math_util.h ---
/* $Id: math_util.h,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2005.  All rights reserved.
 */

#ifndef MATH_UTIL_H
#define MATH_UTIL_H

//#ifndef LINUX
//#include <sunmath.h>
//#endif

#include <vector>
#include <algorithm>
#include <math.h>
#include <iostream>
using namespace::std;

#define Long long long

#define Rand() drand48()
#define SRand(c) srand48(c)

#define ABS(x) ((x > 0) ? x : -1*x)


double randPct ();
int poisson (int mean);
double poisson (double mean);
int unifRand (int a, int b);
vector<int> unifRands (int count, int a, int b);
double mean (vector<double> v);
double mean (vector<int> v);
double stddev (vector<double> v, double mean);
double percentile (vector<double> v, double percentile);
int percentile (vector<int> v, double percentile);
double square (double x);
double log16 (double x);
double log8 (double x);

#endif


--- NEW FILE: histogram.h ---
/* $Id: histogram.h,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2006.  All rights reserved.
 */

#ifndef HISTOGRAM_H
#define HISTOGRAM_H

#include <math.h>
#include <iostream>
#include <vector>
#include <algorithm>
#include <map>
#include "error.h"
using namespace::std;

class Histogram {
  friend ostream& operator << (ostream& os, Histogram* h);
 private:
  int size;
  double sum;
  // TODO some kind of sorted list would be better
  vector<double> values;
  map<double,int> value2count;
  void add (double value, int count, bool doCoerce);
  double findBucket (double value);

 public:
  const double bucketWidth;
  Histogram (double bw) : bucketWidth(bw) { size = 0; sum = 0;}
  void add (Histogram *);
  void add (double value, int count) { add (value, count, true); }
  void add (double value) { add (value, 1); }
  void add (int value) { add ((double)value,1); }
  void add (int value, int count) { add ((double)value,count); }
  void clear () { value2count.clear(); size = 0; values.clear(); }
  double mean ();
  double percentile (double pct);
};

#endif

--- NEW FILE: Makefile ---
# $Id: Makefile,v 1.1 2007/03/15 12:20:49 ledlie Exp $

#PROGRAMS =	triangle
#PROGRAMS =	triangleRPL
#PROGRAMS =	gossip
#PROGRAMS =	coord2nn vivaldi
PROGRAMS =	rs
#PROGRAMS =	nn
#PROGRAMS =	dist_driver
#PROGRAMS =	triangle genPerfectCoords rs vivaldi
#PROGRAMS =	genPerfectCoords
#PROGRAMS =	vivaldiChurn
#PROGRAMS =	histogram
#PROGRAMS = zone
#PROGRAMS =	vivaldi

LD =            ld

OPT =           -O2
#OPT =           -O0


#CXX =           CC
#CXX =           c++
CXX =           g++
#CXX =           purify CC
#CXX =           quantify CC
#CXX =           cc

#SYMBOLS =       
SYMBOLS =       -g
#SYMBOLS =       -Wall -g3 -O0
#SYMBOLS =       -Wall
#GPROF = 	-pg

#LDFLAGS =	-static -lm

LDFLAGS =	-lm 
#LDFLAGS =	-lsunmath -lm
#LDFLAGS =	-lm -lstlport -static
#LDFLAGS =	-lm -lstlport -lpthread
#LDFLAGS =	-lm -lstlport -lpthread -static
#LDFLAGS =	-lm -lstlport_gcc

#INC =		-I/deas/software/include
#INC =		-I/usr/include/stlport
#LIBS =		-L/deas/software/lib
#LIBS =		-L/usr/lib
#INC =		-I/usr/include/stlport

CXXFLAGS =      $(SYMBOLS) $(GPROF) $(OPT) $(STLPORT_FLAGS) $(INC) $(GPROF) $(LIBS)
CFLAGS =      $(SYMBOLS) $(GPROF) $(OPT) $(INC) $(GPROF) $(LIBS)
MAKEFILE =      Makefile


GOSSIP_OBJ = gossip.o
NN_OBJ = math_util.o point.o node.o samples.o nn.o
TRIANGLE_OBJ = triangle.o
TRIANGLE_RPL_OBJ = triangleRPL.o
SIM_OBJ = math_util.o point.o node.o samples.o sim_aux.o histogram.o distributions.o
HISTOGRAM_OBJ = histogram.o
VIVALDI_OBJ = $(SIM_OBJ) vivaldi.o vs.o
RS_OBJ = $(SIM_OBJ) rs_stat.o rs_aux.o rs_alg_greedy.o rs.o
DIST_DRIVER_OBJ = math_util.o distributions.o dist_driver.o
ZONE_OBJ = point.o math_util.o zoneTest.o
GENPERFECTCOORDS_OBJ = math_util.o sim_aux.o point.o genPerfectCoords.o
COORD2NN_OBJ = math_util.o sim_aux.o point.o node.o samples.o coord2nn.o
VIVALDI_SRC = error.h distributions.h histogram.h math_util.h node.h point.h samples.h sim.h vivaldi.h math_util.cc point.cc node.cc samples.cc sim_aux.cc histogram.cc distributions.cc vivaldi.cc vs.cc Makefile

all:            $(PROGRAMS)

vivaldiChurn:	$(VIVALDI_OBJ)
	$(CXX) -o $@ $(VIVALDI_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

vivaldi:	$(VIVALDI_OBJ)
	$(CXX) -o $@ $(VIVALDI_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

nn:	$(NN_OBJ)
	$(CXX) -o $@ $(NN_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

gossip:	$(GOSSIP_OBJ)
	$(CXX) -o $@ $(GOSSIP_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

# keep rs obj around for other simulations that are currently using it
rs:	$(RS_OBJ)
	$(CXX) -o $@ $(RS_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

histogram:	$(HISTOGRAM_OBJ)
	$(CXX) -o $@ $(HISTOGRAM_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

dist_driver:	$(DIST_DRIVER_OBJ)
	$(CXX) -o $@ $(DIST_DRIVER_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

genPerfectCoords:	$(GENPERFECTCOORDS_OBJ)
	$(CXX) -o $@ $(GENPERFECTCOORDS_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

coord2nn:	$(COORD2NN_OBJ)
	$(CXX) -o $@ $(COORD2NN_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

zone:	$(ZONE_OBJ)
	$(CXX) -o $@ $(ZONE_OBJ) $(LDFLAGS) $(LIBS) $(GPROF)

# GEODATA->/local/data/nc/geo

runVS: vivaldi
	./vivaldi -n 2500 -h 0 -g sSAnR -k 0 -m e -t 64 -w 64  -l /wild/geo/meridian/meridian.rtt -r 5000000 -I 10000 -s 2 -O /tmp/4dh -x 4 -z


runRS: all
#	./rs -c $(GEODATA)/perfect/x2.n10.coord -x 2 -w 1 -r $(GEODATA)/perfect/x2.n10.rtt -n 8 -t 2 -S 1 -R 10 -O foo -K U/2 -D -m t
#	./rs -c x2/x2.n10.coord -x 2 -w 1 -r x2/x2.n10.rtt -n 8 -t 2 -S 1 -R 10 -O foo -K U/2 -m g
#	./rs -c /wild/geo/perfect/x2.n100.coord -x 2 -w 1 -n 90 -t 10 -S 1 -R 1000000 -O /tmp/p -K U/10 -m h -l 3 -s 8 -C 100000 -g jn -V r -b 8
#	./rs -c /wild/geo/perfect/x2.n100.coord -x 2 -w 1 -n 90 -t 10 -S 1 -R 10000 -O /tmp/p -K U/10 -m h -l 3 -s 8 -b 8 -C 500 -g jn -P 5000
#	./rs -c /wild/geo/perfect/x2.n2048.coord -x 2 -w 1 -n 2048 -t 0 -S 23326 -R 1 -O /tmp/p -K U/2048 -m h -l 4 -s 13
#	./rs -c x2.n100.v3.coord -x 2 -w 100 -n 90 -t 10 -S 1 -R 1 -O /tmp/p -K U/10 -m h -s 4 -l 10 -b 8
#	./rs -c /wild/geo/mit/MIT.n100.x2.wH.appcoord -x 2 -w 100 -n 90 -t 10 -S 1 -R 1 -O /tmp/p -K U/10 -m h -s 4 -l 10 -b 8 -z -r /wild/geo/mit/MIT.n100.rtt
#	./rs -c x3.n1024.coord -x 3 -w 10 -n 900 -t 124 -S 1 -R 1 -O /tmp/p -K U/124 -m h -s 4 -l 10 -b 8
	./rs -c x2.n1024.coord -x 2 -w 10 -n 900 -t 124 -S 1 -R 1 -O /tmp/p -K U/124 -m h -s 4 -l 10 -b 8


tarball:
	tar cf vivaldi-sim.tar $(VIVALDI_SRC)
	gzip -f vivaldi-sim.tar

.cc.o:
	$(CXX) $(CXXFLAGS) -c $< -o $@

.c.o:
	$(C) $(CFLAGS) -c $< -o $@

clean:
	rm -f core* *~ *.o $(PROGRAMS) *.o

--- NEW FILE: energy.cc ---
#include "vivaldi.h"
#include "node.h"

bool debug = false;

void addToWindow (Point *p, deque<Point*> &initialDist, 
		  deque<Point*> &currentDist, int windowSize) {

    if (A.size() < windowSize) {
      A.push_back (p);
    } else {
      B.push_back (p);
      if (B.size() > windowSize) {
	cout << "need to drop" << endl;
	B.pop_front ();
      }
    }

    if (A.size() >= windowSize &&
	B.size() >= windowSize) {
      double dist = compareDistributions (A, B);
      cout << "dist = " << dist << endl;
    }

}

double findEnergyDistance
(const deque<Point*> &A, const deque<Point*> &B) {
  double sum = 0., abSum = 0., aSum = 0., bSum = 0.;

  if (debug) {
    cout << "A size " << A.size() 
	 << " B size " << B.size() << endl;

    cout << "A: " << endl;
    for (int i = 0; i < A.size(); i++) {
      cout << i << ": " << *A[i] << endl;
    }

    cout << "B: " << endl;
    for (int i = 0; i < B.size(); i++) {
      cout << i << ": " << *B[i] << endl;
    }
  }

  for (int i = 0; i < A.size(); i++) {
    for (int j = 0; j < B.size(); j++) {
      abSum += A[i]->getEuclideanDistance (B[j]);
    }
  }
  abSum *= 2./(A.size()*B.size());
  if (debug)
    cout << "abSum = " << abSum << endl;


  for (int i = 0; i < A.size(); i++) {
    for (int j = 0; j < A.size(); j++) {
      aSum += A[i]->getEuclideanDistance (A[j]);
    }
  }
  aSum *= 1./(A.size()*A.size());
  if (debug)
    cout << "aSum = " << aSum << endl;

  for (int i = 0; i < B.size(); i++) {
    for (int j = 0; j < B.size(); j++) {
      bSum += B[i]->getEuclideanDistance (B[j]);
    }
  }
  bSum *= 1./(B.size()*B.size());
  if (debug)
    cout << "bSum = " << bSum << endl;

  sum = abSum - aSum - bSum;
  if (debug)
    cout << "psum = " << sum << endl;

  sum *= (A.size()*B.size())/(A.size()+B.size());
  return sum;
}

char *usage = \
				"energy\n" \
				" -f sample file (3-d)\n" \
				" -d debug\n" \
				" -w window size\n";

int main (int argc, char **argv) {
  char c =  0;
  FILE *fp;
  char *sampleFile = NULL;
  int windowSize = 0;

  deque<Point*> A, B;

  while ((c = getopt(argc,argv,"w:f:d")) >= 0) {
    switch (c) {
    case 'f':
      sampleFile = optarg;
      break;
    case 'w':
      windowSize = atoi(optarg);
      break;
    case 'd':
      debug = true;
      break;
    default:
      printf (usage);
      exit (-1);
      break;
    }
  }

  if (sampleFile == NULL || windowSize <= 0) {
    printf (usage);
    exit (-1);
  }    

  if ((fp = fopen (sampleFile, "r")) == NULL) {
    printf ("Cannot open file %s", sampleFile);
    exit (-1);
  }

  if (DIMENSIONS != 3) {
    printf ("Bad number of dimenions\n");
    exit (-1);
  }

  float x, y, z;
  double *v = new double[DIMENSIONS];
  while (fscanf (fp, "%f %f %f\n", &x, &y, &z) > 0) {
    //while (fscanf (fp, "%f %f\n", &x, &y, &z) > 0) {

    //printf ("x %f y %f\n", x, y);
    printf ("x %f y %f z %f\n", x, y, z);
    v[0] = x;
    v[1] = y;
    v[2] = z;

    Point *p = new Point (v);
    cout << "new point " << *p << endl;

    addToWindow (p, initialDistSmall, currentDistSmall, smallWindowSize);

    if (A.size() < windowSize) {
      A.push_back (p);
    } else {
      B.push_back (p);
      if (B.size() > windowSize) {
	cout << "need to drop" << endl;
	B.pop_front ();
      }
    }

    if (A.size() >= windowSize &&
	B.size() >= windowSize) {
      double dist = compareDistributions (A, B);
      cout << "dist = " << dist << endl;
    }

  }

  fclose (fp);
  return 0;
}

--- NEW FILE: samples.cc ---
/* $Id: samples.cc,v 1.1 2007/03/15 12:20:50 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2005.  All rights reserved.
 */

#include "sim.h"
#include "samples.h"

int PING_HISTORY_COUNT = 4;
double PING_SAMPLE_PERCENTILE = .25;
int MIN_HISTORY_COUNT = 0;

Samples::Samples() {
  ewma = 0;
  jitterCount = 0;
  weightedError = 0;
  vec = new Point ();
  appVector = new Point ();
}

Samples::~Samples() {
  delete vec;
  delete appVector;
}

void Samples::addSample (double sample, int myId, int yourId, 
                           Point *yourCoord, double yourWeightedError,
                           Point *yourAppCoord, int timestamp) {
  stamp = timestamp;

  if (PING_HISTORY_COUNT == -1) {
    ewma = PING_SAMPLE_PERCENTILE * sample +
      (1.-PING_SAMPLE_PERCENTILE) + ewma;
    sample = ewma;
  }

  weightedError = yourWeightedError;
  vec->assign (yourCoord);
  appVector->assign (yourAppCoord);

  // toss the sample from the front
  if (samples.size() > PING_HISTORY_COUNT)
    samples.pop_front ();
  // add this guy to the back
  samples.push_back (sample);
}

double Samples::getSample () {

  if (samples.size() > MIN_HISTORY_COUNT) {
    // sort em
    deque<double> sortedSamples (samples);
    sort (sortedSamples.begin(), sortedSamples.end());

    int percentile = (int)(samples.size() * PING_SAMPLE_PERCENTILE);
    double sample_at_percentile = sortedSamples[percentile];

    return sample_at_percentile;
  } else {
    return -1;
  }
}

void Samples::print () {
  printf ("sample size %d\n", samples.size());
  for (int i = 0; i < samples.size(); i++) {
    printf ("sample %d %f\n", i, samples[i]);
  }
}

ostream& operator << (ostream& os, Samples *s) {
  os << "[sC " << s->vec << " aC " << s->appVector 
     << " wE " << s->weightedError << " s ";
  for (int i = 0; i < s->samples.size(); i++) {
    os << s->samples[i] << " ";
  } 
  os << "]";
  return os;
}

--- NEW FILE: genPerfectCoords.cc ---
/* $Id: genPerfectCoords.cc,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 */

#include "sim.h"
#include "node.h"

void printUsage (char* problem) {
  cout 
    << "Usage: genPerfectCoords\n"
    <<"            -x dimensions\n"
    <<"            -s random seed\n"
    <<"            -n node count\n"
    <<"            -S side length\n"
    <<"            -O output file prefix\n"
    <<"            -g output format\n"
    <<"               m: matrix\n"
    <<"               c: coordinates\n"
    <<"               r: rtt file\n"
    <<"\n";

  if (problem) {
    cout << problem << endl;
  }
  exit (-1);
}


int main (int argc, char **argv) {
  char *outputPrefix = NULL;
  char *outputFormat = NULL;
  char c = 0;
  int seed = getpid();
  int nodeCount = 0;
  double sideLength = 10.;

  double sum = 0.;

  USE_HEIGHT = false;

  while ( (c = getopt(argc,argv,"s:S:x:O:n:g:")) >= 0) {
    switch (c) {
    case 'g':
      outputFormat = optarg;
      break;
    case 'n':
      nodeCount = atoi (optarg);
      break;
    case 'x':
      DIMENSIONS = atoi (optarg);
      break;
    case 's':
      seed = atoi (optarg);
      break;
    case 'S':
      sideLength = atof (optarg);
      break;
    case 'O':
      outputPrefix = optarg;
      break;
    default:
      printUsage ("Undefined option");
      break;
    }
  }
  
  //printf ("seed %d\n", seed);

  SRand (seed);
  srandom (seed);
  if (outputPrefix == NULL)
    printUsage("Output file prefix is required.");
  if (DIMENSIONS <= 0)
    printUsage("Dimensions must be greater than 0");
  if (outputFormat == NULL)
    printUsage("Output format required.");

  FILE *coordOutputFP, *rttOutputFP, *matrixOutputFP;

  Point **points = new Point* [nodeCount];
  for (int i = 0; i < nodeCount; i++) {
    points[i] = Point::getRandomPoint (sideLength);
  }

  if (strchr(outputFormat,'c') != NULL) {
    openOutputFile (coordOutputFP, outputPrefix, ".coord");
    for (int i = 0; i < nodeCount; i++) {  
      fprintf (coordOutputFP, "%d ", i);
      points[i]->print (coordOutputFP);
      fprintf (coordOutputFP, "\n");
    }

    closeOutputFile (coordOutputFP);
  }

  int pointCount = 0;
  double avg = 0.;
  if (strchr(outputFormat,'r') != NULL) {
    openOutputFile (rttOutputFP, outputPrefix, ".rtt");

    for (int src = 0; src < nodeCount; src++) {  
      for (int dst = 0; dst < nodeCount; dst++) {  
	if (src != dst) {
	  double dist = points[src]->getEuclideanDistance(points[dst]);
	  //dist = dist * dist;
	  sum += dist;
	  pointCount++;
	  ASSERT (dist > 0);
	  fprintf (rttOutputFP, "%d %d %f\n", src, dst, dist);
	}
      }
    }
    closeOutputFile (rttOutputFP);
    if (pointCount > 0)
      avg = sum / (double)pointCount;
    fprintf (stderr, "avg dist between points= %f\n", avg);
  }

  if (strchr(outputFormat,'m') != NULL) {
    openOutputFile (matrixOutputFP, outputPrefix, ".matrix");
    for (int src = 0; src < nodeCount; src++) {  
      for (int dst = 0; dst < nodeCount; dst++) {  
	double dist = points[src]->getEuclideanDistance(points[dst]);
	ASSERT (dist >= 0.);
	if (dst > 0) fprintf (matrixOutputFP," ");
	fprintf (matrixOutputFP, "%f", dist);
      }
      fprintf (matrixOutputFP,"\n");
    }
    closeOutputFile (matrixOutputFP);
  }


}

--- NEW FILE: point.h ---
/* $Id: point.h,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2006.  All rights reserved.
 */

#ifndef POINT_H
#define POINT_H

class Point {
  friend ostream& operator << (ostream& os, const Point *p);
public:
  double *v;
  double height;
  int stamp;

  Point ();
  Point (int stamp);
  Point (const Point *v);
  Point (const double *vec, double h);
  ~Point ();
  void clear ();
  void scale (double s);
  double length ();
  double lengthWithoutHeight ();
  void bump ();
  void minibump ();
  void add (const Point *adjustment);
  bool isInInitialState ();
  void checkHeight ();
  double getAbsoluteDistance (const Point *b) const;
  double getEuclideanDistance (const Point *b) const;
  double getPlanarDistance (const Point *b) const;
  void assign (const Point *vec);
  Point* getDirection (const Point *p) const;
  Point* getForce (const Point *p);
  bool operator<(const Point *rhs) const;
  bool operator<(int rhs);
  static Point* getRandomPoint (double sideLength);
  static Point* getRandomUnitVector ();
  void print (FILE *fp) const;
};



#endif

--- NEW FILE: rs_alg_simple.cc ---
/* $Id: rs_alg_simple.cc,v 1.1 2007/03/15 12:20:50 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2005.  All rights reserved.
 */

#include "sim.h"
#include "node.h"
#include "rs.h"
#include "rs_alg_simple.h"

NNSimple::NNSimple(int nodeCount, int maxNeigbours) : maxNeighbours(maxNeighbours) { 
	neighbours = new set<int>[nodeCount]; 
}

void NNSimple::initPerfectRoutingTables(Node& node) {
  // would be nice if the number of neighbors known in each case was the same
  // and were log(n)
  // or 4 x dim

  if (options && strchr(options, 'n') != NULL) {
    // nearest-only

  } else {
    // random
    
  for (int i = 0; i < maxNeighbors; i++) {
    int yourId = node.id;
    // isNeighbor(myId,yourId)) {
    while (yourId == node.id || rtts[node.id][yourId] <= 0.
	   || contains(neighbours[node.id], yourId)) {
      yourId = unifRand (0, nodeCount);
    }
    neighbours[node.id].insert (yourId);
  }
  }
  
}

void NNSimple::printState(Node& node) {
	
}

void NNSimple::route(Node& srcNode, int targetNodeId, int &resultNodeId, deque<int> &hopNodeIds, double &delay) {

  // TODO problem with this is we are assuming we know the coord
  // for the targetNodeId, i.e. that it has been computed via
  // SBON's spring-relaxation or via proxy coords.
  // So, if we add coord computation and update to the sim
  // then these coords of destination nodes 
  // are going to need to evolve along with it

  // TODO look at SBON code

  // add simple greedy search first

  // two choices??

  // Note: assumes current knowledge of all of our neighbors coordinates
  // This is fine if they aren't changing

  Point *target = nodes[targetNodeId].vec;

  // could create nnQuery object like in Java
  set<int> waiting;
  set<int> completed;

  // add ourselves to kick off the query
  waiting.insert (srcNode.id);
  if (debug) {
    cout << "starting query for target "<< target << endl;
  }

  while (waiting.size() > 0) {
    set<int>::iterator waitingIter = waiting.begin();
    int currentHop = *waitingIter;
    hopNodeIds.push_back (currentHop);
    waiting.erase (waitingIter);

    set<int> suggestedNodeIds;
    if (debug) {
      cout << "using node "<< currentHop << endl;
    }
    findNNQueryHops(nodes[currentHop], target, *neighbours, suggestedNodeIds);
    completed.insert(currentHop);

    for (set<int>::const_iterator suggestedIter = suggestedNodeIds.begin(); 
         suggestedIter != suggestedNodeIds.end(); suggestedIter++) {
      int suggestedNodeId = *suggestedIter;
      if (!contains(waiting,suggestedNodeId) &&
	        !contains(completed,suggestedNodeId)) {
	    	waiting.insert(suggestedNodeId);
      }
    }
  }

  hopNodeIds.pop_front ();

  map<double,int> dist2id;
  orderByDistance (target, completed, dist2id);

  map<double,int>::const_iterator dist2idIter = dist2id.begin();
  resultNodeId = dist2idIter->second;

}

void nnQueryOptimal (int targetNodeId, int resultNodeId, 
		     int &nearestNodeId,
		     int &rankOfResult,
		     double &relLatPenalty, double &absLatPenalty) {

  map<int,vector<int> >::const_iterator targetIter =
    targetId2rank2nodeId.find (targetNodeId);
  ASSERT (targetIter != targetId2rank2nodeId.end());

  vector<int> rank2nodeId = targetIter->second;

  bool found = false;
  for (int i = 0; i < rank2nodeId.size(); i++) {
    //printf ("i %d n %d res %d\n", i, rank2nodeId[i], resultNodeId);
    if (rank2nodeId[i] == resultNodeId) {
      rankOfResult = i;
      found = true;
      break;
    }
  }
  ASSERT (found);

  //vector<int>::const_iterator nearestIter = nodeId2rank.begin();
  //nearestNodeId = nearestIter->second;
  nearestNodeId = rank2nodeId[0];
  ASSERT (nearestNodeId >= 0 && nearestNodeId < nodeCount);

  double resultLat = rtts[targetNodeId][resultNodeId];
  double nearestLat = rtts[targetNodeId][nearestNodeId];
  ASSERT (resultLat > 0);
  ASSERT (nearestLat > 0);
  ASSERT (nearestLat <= resultLat);

  absLatPenalty = resultLat - nearestLat;
  relLatPenalty = absLatPenalty / nearestLat;

}

// returns number of messages used in maintenance
int Node::takeTurn () {
  return 0;
  //DEBUG_P (("R %d NODE %d taking turn\n", cRound, id));
  /*
  if (noChurn) {

  } else {
    // maintain routing tables
    // ping neighbors
    ASSERT (false);

  }
  */
}

void Node::assignCoords () {
    map<int,Point*>::const_iterator iter = nodeId2precoord.find(id);
    ASSERT (iter != nodeId2precoord.end());
    vec->assign (iter->second);
}

void Node::birth () {
  ASSERT (!alive);
  nodeUp++;
  alive = true;
  birthTime = cRound;
  deathTime = -1;
  DEBUG_P (("R %d NODE %d starting birth\n", cRound, id));

  //if (noChurn) {
    // set our own pre-computed coord
  	assignCoords ();

		// TODO FIXME: initialise the routing tables
    //nnInit();
    assert (0);
    /*
  } else {

    // learn about five nodes that are up
    // setNextRepairRoutesTime ();
    ASSERT (false);
  }
    */
}

void Node::death () {
  ASSERT (alive);
  DEBUG_P (("R %d NODE %d death\n", cRound, id));

  nodeUp--;
  alive = false;
  deathTime = cRound;
}

bool Node::isAlive () {
  return alive;
}

/*
void Node::setNextRepairRoutesTime () {
  nextRepairRoutesTime = cRound + (int)(poisson ((double)repairRoutesActionInterval));
}
*/

--- NEW FILE: rs_alg_greedy.cc ---
/* $Id: rs_alg_greedy.cc,v 1.1 2007/03/15 12:20:49 ledlie Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2005.  All rights reserved.
 */

#include "sim.h"
#include "node.h"
#include "rs.h"
#include "rs_alg.h"

GossipStrategy gossipStrategy = None;

// This method initialises the routing tables
// by offering to add all nodes to each node's routing table
void NNGreedy::initPerfectRoutingTables(Node& node) {
  DEBUG_P(("initPerfectRoutingTable, nodeCount=%d\n", nodeCount));


  // necessary for initial zone
[...1319 lines suppressed...]
  ASSERT (angle >= -M_PI && angle <= M_PI);

  for (double ray = -M_PI; ray <= M_PI; ray += sectorIncrement) {
    if (angle >= ray && angle < ray+sectorIncrement) {
      break;
    }
    sectorIndex++;
  }

  // happens if angle is PI
  if (sectorIndex == numSectors) sectorIndex--;

  //if (sectorIndex < 0 || sectorIndex >= numSectors)
  ASSERT (sectorIndex >= 0 && sectorIndex < numSectors);
  delete p;

  //if (debug) printf ("me %d S %d\n", me.id, sectorIndex);
  return sectorIndex;
}


--- NEW FILE: rs_stat.cc ---
/* $Id: rs_stat.cc,v 1.1 2007/03/15 12:20:50 ledlie Exp $
 */

#include "rs.h"

void Statistics::clear () {
  startTime = cRound;

  upNodeCount = 0;
  birthCount = 0;
  deathCount = 0;

  queryCount = 0;
  //queryMsgCount = 0;
  maintMsgCount = 0;

  queryRank.clear();
  queryGlRelLatPenalty.clear();
  queryGlAbsLatPenalty.clear();
  queryNcRelLatPenalty.clear();
  queryNcAbsLatPenalty.clear();
  queryHopCount.clear();
  queryDelay.clear();

  for (int i =0; i < MAX_HOPS; i++) {
    hopCount[i] = 0;
    for (int j =0; j < MAX_HOPS; j++) {
      hopDistance[i][j] = 0.;
    }
  }

  gossipSize = 0.;
  angleTableSize = 0.;
  pctNearbyKnown = 0.;
  gossipStatCount = 0;

}

void Statistics::gossipState(double gS, double aTS, double pNK) {
  gossipStatCount++;
  gossipSize += gS;
  angleTableSize += aTS;
  pctNearbyKnown += pNK;
}

void Statistics::add (Statistics *s) {
  upNodeCount += s->upNodeCount;
  birthCount += s->birthCount;
  deathCount += s->deathCount;

  queryCount += s->queryCount;
  maintMsgCount += s->maintMsgCount;

  gossipStatCount += s->gossipStatCount;
  gossipSize += s->gossipSize;
  angleTableSize += s->angleTableSize;
  pctNearbyKnown += s->pctNearbyKnown;

  for (int i = 0; i < s->queryCount; i++) {
    queryRank.push_back (s->queryRank[i]);
    queryGlRelLatPenalty.push_back (s->queryGlRelLatPenalty[i]);
    queryGlAbsLatPenalty.push_back (s->queryGlAbsLatPenalty[i]);
    queryNcRelLatPenalty.push_back (s->queryNcRelLatPenalty[i]);
    queryNcAbsLatPenalty.push_back (s->queryNcAbsLatPenalty[i]);
    queryHopCount.push_back (s->queryHopCount[i]);
    if (recordDelayStretch)
      queryDelay.push_back (s->queryDelay[i]);
  }

  if (recordDelayStretch) {
  for (int i =0; i < MAX_HOPS; i++) {
    hopCount[i] += s->hopCount[i];
    for (int j =0; j < MAX_HOPS; j++) {
      hopDistance[i][j] += s->hopDistance[i][j];
    }
  }
  }

}

void Statistics::nodesUp (int count) {
  upNodeCount = count;
}

void Statistics::printHops (FILE *fp) {

  int medianHopCount = percentile (queryHopCount,.5);

  // gives the pct covered * median distance
  // i.e. pct of each hop such that total is medianHopCount
  for (int i = 0; i < medianHopCount; i++) {
    fprintf (fp, "%d %.4f\n", i, hopDistance[medianHopCount][i]/(double)hopCount[medianHopCount]*(double)medianHopCount);
  }

}

void Statistics::print (FILE *fp) {
  double diff = cRound - startTime;
  ASSERT (diff >= 0.);
  if (diff < 1.) diff = 1.;

  fprintf (fp, "up %d queries %d ", upNodeCount, queryCount);

  printf ("size %d\n", queryRank.size());

  fprintf (fp, "rk %.2f rk50 %4d rk95 %4d ",
	   mean (queryRank), percentile (queryRank, .5), percentile (queryRank,.95));
  fprintf (fp, "grlp %.2f grlp50 %.2f grlp80 %.2f grlp95 %.2f ",
	   mean (queryGlRelLatPenalty), percentile (queryGlRelLatPenalty,.5), percentile (queryGlRelLatPenalty,.8), 
	   percentile (queryGlRelLatPenalty,.95));
  fprintf (fp, "galp %.2f galp50 %.2f galp80 %.2f galp95 %.2f ",
	   mean (queryGlAbsLatPenalty), percentile (queryGlAbsLatPenalty,.5), percentile (queryGlAbsLatPenalty,.8), 
	   percentile (queryGlAbsLatPenalty,.95));
  fprintf (fp, "nrlp %.2f nrlp50 %.2f nrlp80 %.2f nrlp95 %.2f ",
	   mean (queryNcRelLatPenalty), percentile (queryNcRelLatPenalty,.5), percentile (queryNcRelLatPenalty,.8), 
	   percentile (queryNcRelLatPenalty,.95));
  fprintf (fp, "nalp %.2f nalp50 %.2f nalp80 %.2f nalp95 %.2f ",
	   mean (queryNcAbsLatPenalty), percentile (queryNcAbsLatPenalty,.5), percentile (queryNcAbsLatPenalty,.8), 
	   percentile (queryNcAbsLatPenalty,.95));
  fprintf (fp, "hc %.2f hc50 %4d hc95 %4d ",
	   mean (queryHopCount),
	   percentile (queryHopCount,.5), percentile (queryHopCount, .95));
  fprintf (fp, "de %.2f de50 %.2f de95 %.2f ",
	   mean (queryDelay),
	   percentile (queryDelay,.5), percentile (queryDelay,.95));
  fprintf (fp, "g %.2f at %.2f lk %.4f",
           gossipSize/gossipStatCount, angleTableSize/gossipStatCount, pctNearbyKnown/gossipStatCount);
  fprintf (fp, "\n");

}


Statistics::Statistics () {
  /*
  queryRank = new Histogram (1.);
  queryGlRelLatPenalty = new Hist...
 
[truncated message content]