vivaldinode.c

P2P模拟器

      if (_c[i] > _tsize/2) _c[i] = -(_tsize - _c[i]);
      if (_c[i] < -_tsize/2) _c[i] = _tsize + _c[i];
      assert (fabs(_c[i]) < _tsize/2);
    }
  } else {
    assert (_model_type == -1);
  }

  if (usinght && _c._ht <= 1000) // 1000 is 1ms
    _c._ht = 1000; 

}

void
VivaldiNode::initial_triangulation (vector<Sample> isamps)
{
  double prev_f = RAND_MAX;
  long iterations = 0;

  while (fabs(_pred_err - prev_f) > 0.0001
	 && iterations++ < 10000) {
    Coord f = net_force (_c, isamps);
    _c = _c + (f * 0.05); // XXX what timestep?
    prev_f = _pred_err;
    update_error (isamps);
    if (ip () == 1) cerr << ip () << " " << _pred_err << "\n";
  }

  cerr << ip () << ": " << _pred_err << "\n";
  return;
}

VivaldiNode::Coord
VivaldiNode::real_coords ()
{
  Coord ret;

  Topology *t = Network::Instance ()->gettopology ();
  Euclidean *e = dynamic_cast<Euclidean *>(t);
  if (e) {
    pair<int, int> c = e->getcoords (ip ());
    ret.init2d (c.first, c.second);
  }
  return ret;
}

/*
 * Coordinate manipulation code.
 */

ostream&
operator<< (ostream &s, VivaldiNode::Coord &c) 
{
  for (uint i = 0; i < c._v.size(); i++){
    if (i)
      s << ",";
    s << c._v[i];
  }
  if (usinght)
    s << ",ht=" << c._ht;
  return s;
}


VivaldiNode::Coord
operator-(VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  VivaldiNode::Coord c;
  assert (a._v.size () == b._v.size ());
  for (unsigned int i = 0; i < a._v.size (); i++) 
    c._v.push_back (a._v[i] - b._v[i]);
  if (usinght)
    c._ht = a._ht+b._ht;
  return c;
}

VivaldiNode::Coord
operator+(VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  VivaldiNode::Coord c;
  assert (a._v.size () == b._v.size ());
  for (unsigned int i = 0; i < a._v.size (); i++) 
    c._v.push_back(a._v[i] + b._v[i]);
  if (usinght)
    c._ht = a._ht+b._ht;
  return c;
}

VivaldiNode::Coord
operator/(VivaldiNode::Coord c, double x)
{
  for (unsigned int i = 0; i < c._v.size (); i++) 
    c._v[i] /= x;
  if (usinght)
    c._ht /= x;
  return c;
}

VivaldiNode::Coord
operator*(VivaldiNode::Coord c, double x)
{
  for (unsigned int i = 0; i < c._v.size (); i++) 
    c._v[i] *= x;
  if (usinght)
    c._ht *= x;
  return c;
}

double
operator* (VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  assert (a.dim () == b.dim ());
  double ret = 0.0;
  for (int i = 0; i < a.dim (); i++)
    ret += a[i]*b[i];
  return ret;
}

double
length(VivaldiNode::Coord c)
{
  double l = 0.0;
  for (unsigned int i = 0; i < c._v.size (); i++) 
    l += c._v[i]*c._v[i];
  l = sqrt(l);
  if (usinght)
    l += c._ht;
  return l;
}

double
plane_length(VivaldiNode::Coord c)
{
  double l = 0.0;
  for (unsigned int i = 0; i < c._v.size (); i++) 
    l += c._v[i]*c._v[i];
  l = sqrt(l);
  return l;
}


VivaldiNode::Coord 
cross (VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  assert (a.dim () == b.dim ());

  VivaldiNode::Coord ret(a.dim ());
  ret[0] = a[1]*b[2] - a[2]*b[1];
  ret[1] = a[2]*b[0] - a[0]*b[2];
  ret[2] = a[0]*b[1] - a[1]*b[0];
  return ret;
}

VivaldiNode::Coord
rotate_arb (VivaldiNode::Coord axis, VivaldiNode::Coord vec,
	    double angle)
{
  //unitize axis
  VivaldiNode::Coord a = axis / (length (axis));
  
  //calculate the rotation matrix
  // do the multiply
  assert (a.dim () == vec.dim ());
  VivaldiNode::Coord ret (vec.dim());
  
  //helpful intermediates
  double cos_a = cos (angle);
  double sin_a = sin (angle);
  //let's just write out all nine components of the rotation
  // matrix. I know that this is slow. Hopefully -O2 will eat this up.
  // This is from Strang, Intro. To Linear Algebra, p. 371
  double M_11 = cos_a + (1 - cos_a)*a[0]*a[0];
  double M_12 = (1 - cos_a)*a[0]*a[1] - sin_a*a[2];
  double M_13 = (1 - cos_a)*a[0]*a[2] + sin_a*a[1];
  double M_21 = (1 - cos_a)*a[0]*a[1] + sin_a*a[2];
  double M_22 = cos_a + (1 - cos_a)*a[1]*a[1];
  double M_23 = (1 - cos_a)*a[1]*a[2] - sin_a*a[0];
  double M_31 = (1 - cos_a)*a[0]*a[2] - sin_a*a[1];
  double M_32 = (1 - cos_a)*a[1]*a[2] + sin_a*a[0];
  double M_33 = cos_a + (1 - cos_a)*a[2]*a[2];

  //now let's write out the result of the matrix
  // multiplication in terms of the variables above

  ret[0] = M_11*vec[0] + M_12*vec[1] + M_13*vec[2];
  ret[1] = M_21*vec[0] + M_22*vec[1] + M_23*vec[2];
  ret[2] = M_31*vec[0] + M_32*vec[1] + M_33*vec[2];
  
  return ret;

}

double
flatearth_dist(VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  double d = 0.0;
  assert (a._v.size () == b._v.size ());
  for (unsigned int i = 0; i < a._v.size (); i++) 
    d += (a._v[i] - b._v[i])*(a._v[i] - b._v[i]);
  d = sqrt(d);
  if (usinght)
    d += a._ht + b._ht;
  return d;
}

double 
VivaldiNode::toroidal_dist (VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  assert (a.dim () == b.dim());
  double d_sum = 0.0;
  for (int i = 0; i < a.dim (); i++)
    {
      double x1, x2;
      if (a[i] > b[i]) { x1 = a[i]; x2 = b[i]; }
      else { x1 = b[i]; x2 = a[i]; }
      double d1 = x1 - x2;
      double d2 = (_tsize - x2) + x1;
      double d;
      if (d1 > d2) d = d1; 
      else d = d2;
      d_sum += d*d;
    }
  return sqrt (d_sum);
}


// return the angle between to points on a sphere in radians
double 
VivaldiNode::spherical_dist_arc (VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  //coords are in cartesian space, vectors to points on a sphere
  // of radius RADIUS
  assert (a._v.size () == 3);

  //vectors are constructed to be RADIUS long, so 
  // we don't calculate the norm
  double cos_angle = (a * b)/(_radius*_radius); 
  
  return acos (cos_angle);
}

double 
VivaldiNode::spherical_dist (VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  double arc_dist = spherical_dist_arc (a, b);
  return arc_dist * _radius;
}


double
VivaldiNode::dist(VivaldiNode::Coord a, VivaldiNode::Coord b)
{
  if (_model_type == MODEL_SPHERE)
    return spherical_dist (a, b);
  else if (_model_type == MODEL_EUCLIDEAN) 
    return flatearth_dist (a, b);
  else 
    return toroidal_dist (a,b);
    
}

/*
 * Use generic simplex code to determine the best placement.
 */

static inline double
square(double x)
{
  return x*x;
}

static double
simplex_error(double *x, int dim, void *v)
{
  VivaldiNode *node;
  
  node = (VivaldiNode*)v;
  if(usinght)
    dim--;
  VivaldiNode::Coord c(0);
  for(int i=0; i<dim; i++){
    c._v.push_back(x[i]);
  }
  if(usinght)
    c._ht = x[dim];
    
  double tot=0;
  for(uint i=0; i<node->_samples.size(); i++){
    double estimate = node->dist(c, node->_samples[i]._c);
    double actual = node->_samples[i]._latency;
    tot += square((actual-estimate)/actual);
  }
  return tot;
}

static VivaldiNode::Coord
run_perfect_spring(VivaldiNode *node)
{
  VivaldiNode::Coord c = node->_c;
  VivaldiNode::Coord f;

  do{
    f = node->net_force(c, node->_samples);
    c = c+f;
  }while(length(f) > .01);
  return c;
}

static VivaldiNode::Coord
run_simplex(VivaldiNode *node)
{
  int dim = node->_samples[0]._c._v.size();
  int dimh = dim + (usinght!=0);

  Simplex *plex = allocsimplex(dimh, simplex_error, (void*)node, .01);
  double *best = 0;
  double ans;
  for(int i=0; i<1000; i++)
    if(!stepsimplex(plex, &ans, &best))
      break;
  VivaldiNode::Coord c(0);
  for(int i=0; i<dim; i++){
    c._v.push_back(best[i]);
  }
  if(usinght)
    c._ht = best[dim];
  free(plex);
  return c;
}