ranv2.cc

obs网络试验平台

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "ranv2.h"
#include "ranvar.h"
#include "rng.h"

static class GammaRandomVariableClass : public TclClass {
public:
	GammaRandomVariableClass() : TclClass("RandomVariable/Gamma") {}
	TclObject* create(int, const char*const*) {
		return(new GammaRandomVariable());
	}
} class_gammaranvar;

GammaRandomVariable::GammaRandomVariable()
{
	bind("avg", &avg_);
	bind("stdev", &stdev_);
}

GammaRandomVariable::GammaRandomVariable(double avg, double stdev)
{
        avg_ = avg;
        stdev_ = stdev;
}

double GammaRandomVariable::value()
{
     double alpha, beta;

     if (stdev_ <= 0.0)
        return (avg_);

     alpha = avg_*avg_/(stdev_*stdev_);
     beta = stdev_*stdev_/avg_; // We will ensure avg_ != 0 by forcing the user to enter a finite rate for
                                // the traffic and taking the inverse
     return( Gamma (beta, alpha) );
}

double GammaRandomVariable::erlang(double mean, int rep)
{
  double sum;
  int i;
  
  for (sum = 1, i = 1; i <= rep; i++)
    sum *= rng_->uniform (1.0);
  return( -mean * log(sum));
}

double GammaRandomVariable::Gamma(double beta, double alpha)
{
  double x, y, b, tmp;
  int a;
  
  x = y = 1;
  tmp = -1;
  if (alpha < 1) {
    while(x + y > 1) {
      x = pow( rng_->uniform(1.0), 1./alpha);
      y = pow( rng_->uniform(1.0), 1./(1. - alpha));
    }
    return(-x*log( rng_->uniform(1.0) )*beta/(x+y));
  } else if (alpha < 5) {
    while( rng_->uniform(1.0) > tmp) {
      a = (int) floor(alpha);
      b = alpha - a;
      x = erlang(alpha / a, a);
      tmp = pow((x/alpha),b)*exp(-b*(x/alpha-1.));
    }
    return(x * beta);
  } else {
    a = floor(alpha);
    if ( rng_->uniform(1.0) < (alpha - a))
      a++;
    x = erlang(beta, a);
    return(x);
  }
}

static class NegBinomRandomVariableClass : public TclClass {
public:
	NegBinomRandomVariableClass() : TclClass("RandomVariable/NegBinom") {}
	TclObject* create(int, const char*const*) {
		return(new NegBinomRandomVariable());
	}
} class_negbinomranvar;

NegBinomRandomVariable::NegBinomRandomVariable()
{
	bind("avg", &avg_);
	bind("sparm", &sparm_);
}

NegBinomRandomVariable::NegBinomRandomVariable(double avg, int sparm)
{
        avg_ = avg;
        sparm_ = sparm;
}

double NegBinomRandomVariable::value()
{
     double p;

     if (sparm_ <= 0)
        return (avg_);

     p = sparm_/(avg_ + sparm_ - 1.0);
     return( negbinom (sparm_, p) );
}

double NegBinomRandomVariable::geometric0 (double p)
{
     double tmp;

     if (p <= 0.0)
        return (0.0);
     tmp = log (rng_->uniform(1.0))/log(1.0 - p);
     return (floor(tmp));
}

double NegBinomRandomVariable::negbinom (int s, double p)
{
     int i;
     double x;

     for (i = 1, x= 0.0; i <= s; i++)
        x += geometric0 (p);

     x += 1.0;  // p was chosen to correspond to mean equal avg_ minus 1; add one back
                //   here.  This guarantees that the sample will always be positive.

     return (x);
}