distributions.cc

题描述的是一个旅行商要到几个城市去

/* $Id: distributions.cc,v 1.3 2006-08-12 20:50:06 jonathan Exp $
 * Jonathan Ledlie, Harvard University.
 * Copyright 2005.  All rights reserved.
 */

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

/*
 * Parts of this are borrowed (with thanks) from:
 *
 * Perl Math::Random
 * Math::Random was put together by John Venier and Barry W. Brown
 * with help from SWIG.  For version 0.61, Geoffrey Rommel made various
 * cosmetic changes.
 *
 * Kenneth J. Christensen
 * Associate Professor
 * Department of Computer Science and Engineering
 * University of South Florida
 *
 */

float INVALID = -10099.99104;

double ParetoDistribution::next () {
  return (scale/(pow(randPct(),(1./shape))));
};

double ParetoDistribution::getMean() {
  if (shape > 1.)
    return (scale * shape) / (shape - 1.);
  else
    return 0;
}

double NormalDistribution::next () {
  double normal = stddev*snorm()+mean;
  return normal;
};

double NormalDistribution::getMean() {
  return mean;
}

double UniformDistribution::next () {
  return randPct();
};

double UniformDistribution::getMean() {
  return 0.;
}

int UniformIntegerDistribution::next () {
  return unifRand(0,maxElement);
};

double UniformIntegerDistribution::getMean() {
  return maxElement/2.;
}

double GnutellaBWDistribution::next () {
  double z = randPct();

  if (z < .2) {
    return 1.0672;
  } else if (z < .65) {
    return 10.672;
  } else if (z < .95) {
    return 106.72;
  } else if (z < .999) {
    return 1067.2;
  } else {
    return 10672;
  }
};

double GnutellaBWDistribution::getMean() {
  ASSERT (0);
  return 0.;
}

double computeZipfNormalization (double alpha, int elements) {
  double c = 0.;
  for (int i=1; i<=elements; i++) {
    c = c + (1.0 / pow((double) i, alpha));
  }
  c = 1.0 / c;
  return c;
}

ZipfDistribution::ZipfDistribution (double a, int n) : alpha(a), elements(n) {
  c = computeZipfNormalization (alpha, elements);

  // Map z to the value
  double sum_prob = 0.;
  for (int i=1; i<=elements; i++) {
    sum_prob = sum_prob + c / pow((double) i, alpha);
    zipfKeys.insert(pair<double,double>(sum_prob,randPct()));
  }
}

ZipfIntegerDistribution::ZipfIntegerDistribution (double a, int n) :
  alpha(a), elements(n) {
  c = computeZipfNormalization (alpha, elements);

  vector<int> elementList;
  for (int i=0; i<elements; i++) {
    elementList.push_back(i);
  }  
  random_shuffle (elementList.begin(),elementList.end());

  // Map z to the value
  double sum_prob = 0.;
  for (int i=1; i<=elements; i++) {
    sum_prob = sum_prob + c / pow((double) i, alpha);
    zipfKeys.insert(pair<double,int>(sum_prob,elementList[i-1]));
  }
}


double ZipfDistribution::next () {
  double z;                     // Uniform random number (0 < z < 1)
  z = randPct();

  map<double,double>::iterator p;
  p = zipfKeys.upper_bound (z);
  if (p == zipfKeys.end())
    p = zipfKeys.begin();

  return p->second;
};

int ZipfIntegerDistribution::next () {
  double z;                     // Uniform random number (0 < z < 1)
  z = randPct();

  map<double,int>::iterator p;
  p = zipfKeys.upper_bound (z);
  if (p == zipfKeys.end())
    p = zipfKeys.begin();

  return p->second;
};


double ZipfDistribution::getMean() {
  return alpha;
}

double ZipfIntegerDistribution::getMean() {
  return alpha;
}

double PoissonDistribution::next () {
  
  double u;
  while ((u = randPct()) == 0.);
  return -mean * log(u);
};

double PoissonDistribution::getMean() {
  return mean;
}

double ConstantDistribution::next () {
  return value;
};

double ConstantDistribution::getMean() {
  return value;
}

IntegerDistribution* DistributionFactory::newIntegerDistribution (char *str) {
  char distType;
  float para1, para2;
  getParameters (str, distType, para1, para2);
  //printf ("%c p1 %f p2 %f\n", distType, para1, para2);

  switch (distType) {
  case 'Z':
    if (para1 == INVALID || para2 == INVALID) return NULL;
    return new ZipfIntegerDistribution ((double)para1, (int)para2);
    break;
  case 'U':
    if (para1 == INVALID || para2 != INVALID) return NULL;
    return new UniformIntegerDistribution ((int)para1);
    break;
  default:
    return NULL;    
  }

  return NULL;
}



void getParameters (char *str, char &distType,
		    float &para1, float &para2) {
  para1 = INVALID;
  para2 = INVALID;
  int ret;

  if ( (sscanf (str, "%c/%f/%f", &distType, &para1, &para2) != 3) &&
       (sscanf (str, "%c/%f", &distType, &para1) != 2) &&
       (sscanf (str, "%c", &distType) != 1)) {
    return;
  }
}

Distribution* DistributionFactory::newDistribution (char *str) {
  char distType;
  float para1, para2;
  getParameters (str, distType, para1, para2);

  switch (distType) {
  case 'P':
    if (para1 == INVALID || para2 == INVALID) return NULL;
    return new ParetoDistribution (para1, para2);
    break;
  case 'N':
    if (para1 == INVALID || para2 == INVALID) return NULL;
    return new NormalDistribution (para1, para2);
    break;
  case 'Z':
    if (para1 == INVALID || para2 == INVALID) return NULL;
    return new ZipfDistribution ((double)para1, (int)para2);
    break;
  case 'U':
    if (para1 != INVALID || para2 != INVALID) return NULL;
    return new UniformDistribution ();
    break;
  case 'F':
    if (para1 == INVALID || para2 != INVALID) return NULL;
    return new PoissonDistribution (para1);
    break;
  case 'C':
    if (para1 == INVALID || para2 != INVALID) return NULL;
    return new ConstantDistribution (para1);
    break;
  case 'G':
    if (para1 != INVALID || para2 != INVALID) return NULL;
    return new GnutellaBWDistribution ();
    break;
  default:
    return NULL;    
  }

  return NULL;
}

double snorm(void)
/*
**********************************************************************
                                                                      
                                                                      
     (STANDARD-)  N O R M A L  DISTRIBUTION                           
                                                                      
                                                                      
**********************************************************************
**********************************************************************
                                                                      
     FOR DETAILS SEE:                                                 
                                                                      
               AHRENS, J.H. AND DIETER, U.                            
               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             
               SAMPLING FROM THE NORMAL DISTRIBUTION.                 
               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          
                                                                      
     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  
     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  
                                                                      
     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
     SUNIF.  The argument IR thus goes away.                          
                                                                      
**********************************************************************
     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
*/
{
static double a[32] = {
    0.0,              0.03917608550309, 0.07841241273311, 0.11776987457909,
    0.15731068461017, 0.19709908429430, 0.23720210932878, 0.27769043982157,
    0.31863936396437, 0.36012989178957, 0.40225006532172, 0.44509652498551,
    0.48877641111466, 0.53340970624127, 0.57913216225555, 0.62609901234641,
    0.67448975019607, 0.72451438349236, 0.77642176114792, 0.83051087820539,
    0.88714655901887, 0.94678175630104, 1.00999016924958, 1.07751556704027,
    1.15034938037600, 1.22985875921658, 1.31801089730353, 1.41779713799625,
    1.53412054435253, 1.67593972277344, 1.86273186742164, 2.15387469406144
};
static double d[31] = {
    0.0, 0.0, 0.0, 0.0,
    0.0,              0.26368432217502, 0.24250845238097, 0.22556744380930,
    0.21163416577204, 0.19992426749317, 0.18991075842246, 0.18122518100691,
    0.17360140038056, 0.16684190866667, 0.16079672918053, 0.15534971747692,
    0.15040938382813, 0.14590257684509, 0.14177003276856, 0.13796317369537,
    0.13444176150074, 0.13117215026483, 0.12812596512583, 0.12527909006226,
    0.12261088288608, 0.12010355965651, 0.11774170701949, 0.11551189226063,
    0.11340234879117, 0.11140272044119, 0.10950385201710
};
static double t[31] = {
    7.6738283767E-4,   2.30687039764E-3,  3.86061844387E-3,  5.43845406707E-3,
    7.05069876857E-3,  8.70839582019E-3,  1.042356984914E-2, 1.220953194966E-2,
    1.408124734637E-2, 1.605578804548E-2, 1.815290075142E-2, 2.039573175398E-2,
    2.281176732513E-2, 2.543407332319E-2, 2.830295595118E-2, 3.146822492920E-2,
    3.499233438388E-2, 3.895482964836E-2, 4.345878381672E-2, 4.864034918076E-2,
    5.468333844273E-2, 6.184222395816E-2, 7.047982761667E-2, 8.113194985866E-2,
    9.462443534514E-2, 0.11230007889456,  0.13649799954975,  0.17168856004707,
    0.22762405488269,  0.33049802776911,  0.58470309390507
};
static double h[31] = {
    3.920617164634E-2, 3.932704963665E-2, 3.950999486086E-2, 3.975702679515E-2,
    4.007092772490E-2, 4.045532602655E-2, 4.091480886081E-2, 4.145507115859E-2,
    4.208311051344E-2, 4.280748137995E-2, 4.363862733472E-2, 4.458931789605E-2,
    4.567522779560E-2, 4.691571371696E-2, 4.833486978119E-2, 4.996298427702E-2,
    5.183858644724E-2, 5.401138183398E-2, 5.654656186515E-2, 5.953130423884E-2,
    6.308488965373E-2, 6.737503494905E-2, 7.264543556657E-2, 7.926471414968E-2,
    8.781922325338E-2, 9.930398323927E-2, 0.11555994154118,  0.14043438342816,
    0.18361418337460,  0.27900163464163,  0.70104742502766
};
static long i;
static double snormboth,u,s,ustar,aa,w,y,tt;
//    u = ranf();
    u = randPct();
    s = 0.0;
    if(u > 0.5) s = 1.0;
    u += (u-s);
    u = 32.0*u;
    i = (long) (u);
    if(i == 32) i = 31;
    if(i == 0) goto S100;
/*
                                START CENTER
*/
    ustar = u-(double)i;
    aa = *(a+i-1);
S40:
    if(ustar <= *(t+i-1)) goto S60;
    w = (ustar-*(t+i-1))**(h+i-1);
S50:
/*
                                EXIT   (BOTH CASES)
*/
    y = aa+w;
    snormboth = y;
    if(s == 1.0) snormboth = -y;
    return snormboth;
S60:
/*
                                CENTER CONTINUED
*/
//    u = ranf();
    u = randPct();
    w = u*(*(a+i)-aa);
    tt = (0.5*w+aa)*w;
    goto S80;
S70:
    tt = u;
    //    ustar = ranf();
    ustar = randPct();
S80:
    if(ustar > tt) goto S50;
    //    u = ranf();
    u = randPct();
    if(ustar >= u) goto S70;
    //    ustar = ranf();
    ustar = randPct();
    goto S40;
S100:
/*
                                START TAIL
*/
    i = 6;
    aa = *(a+31);
    goto S120;
S110:
    aa += *(d+i-1);
    i += 1;
S120:
    u += u;
    if(u < 1.0) goto S110;
    u -= 1.0;
S140:
    w = u**(d+i-1);
    tt = (0.5*w+aa)*w;
    goto S160;
S150:
    tt = u;
S160:
    //    ustar = ranf();
    ustar = randPct();
    if(ustar > tt) goto S50;
    //    u = ranf();
    u = randPct();
    if(ustar >= u) goto S150;
    //    u = ranf();
    u = randPct();
    goto S140;
}