ranvar.cc
上传用户:rrhhcc
上传日期:2015-12-11
资源大小:54129k
文件大小:15k
- /* -*- Mode:C++; c-basic-offset:8; tab-width:8; indent-tabs-mode:t -*- */
- /*
- * Copyright (c) Xerox Corporation 1997. All rights reserved.
- *
- * This program is free software; you can redistribute it and/or modify it
- * under the terms of the GNU General Public License as published by the
- * Free Software Foundation; either version 2 of the License, or (at your
- * option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License along
- * with this program; if not, write to the Free Software Foundation, Inc.,
- * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
- *
- * Linking this file statically or dynamically with other modules is making
- * a combined work based on this file. Thus, the terms and conditions of
- * the GNU General Public License cover the whole combination.
- *
- * In addition, as a special exception, the copyright holders of this file
- * give you permission to combine this file with free software programs or
- * libraries that are released under the GNU LGPL and with code included in
- * the standard release of ns-2 under the Apache 2.0 license or under
- * otherwise-compatible licenses with advertising requirements (or modified
- * versions of such code, with unchanged license). You may copy and
- * distribute such a system following the terms of the GNU GPL for this
- * file and the licenses of the other code concerned, provided that you
- * include the source code of that other code when and as the GNU GPL
- * requires distribution of source code.
- *
- * Note that people who make modified versions of this file are not
- * obligated to grant this special exception for their modified versions;
- * it is their choice whether to do so. The GNU General Public License
- * gives permission to release a modified version without this exception;
- * this exception also makes it possible to release a modified version
- * which carries forward this exception.
- */
- #ifndef lint
- static const char rcsid[] =
- "@(#) $Header: /cvsroot/nsnam/ns-2/tools/ranvar.cc,v 1.22 2008/02/01 21:39:43 tom_henderson Exp $ (Xerox)";
- #endif
- #include <stdio.h>
- #include "ranvar.h"
- RandomVariable::RandomVariable()
- {
- rng_ = RNG::defaultrng();
- }
- int RandomVariable::command(int argc, const char*const* argv)
- {
- Tcl& tcl = Tcl::instance();
- if (argc == 2) {
- if (strcmp(argv[1], "value") == 0) {
- tcl.resultf("%6e", value());
- return(TCL_OK);
- }
- }
- if (argc == 3) {
- if (strcmp(argv[1], "use-rng") == 0) {
- rng_ = (RNG*)TclObject::lookup(argv[2]);
- if (rng_ == 0) {
- tcl.resultf("no such RNG %s", argv[2]);
- return(TCL_ERROR);
- }
- return(TCL_OK);
- }
- }
- return(TclObject::command(argc, argv));
- }
- // Added by Debojyoti Dutta 12 October 2000
- // This allows us to seed a randomvariable with
- // our own RNG object. This command is called from
- // expoo.cc and pareto.cc
- int RandomVariable::seed(char *x){
-
- Tcl& tcl = Tcl::instance();
- rng_ = (RNG*)TclObject::lookup(x);
- if (rng_ == 0) {
- tcl.resultf("no such RNG %s", x);
- return(TCL_ERROR);
- }
- return(TCL_OK);
-
- }
- static class UniformRandomVariableClass : public TclClass {
- public:
- UniformRandomVariableClass() : TclClass("RandomVariable/Uniform"){}
- TclObject* create(int, const char*const*) {
- return(new UniformRandomVariable());
- }
- } class_uniformranvar;
- UniformRandomVariable::UniformRandomVariable()
- {
- bind("min_", &min_);
- bind("max_", &max_);
- }
- UniformRandomVariable::UniformRandomVariable(double min, double max)
- {
- min_ = min;
- max_ = max;
- }
- double UniformRandomVariable::value()
- {
- return(rng_->uniform(min_, max_));
- }
- static class ExponentialRandomVariableClass : public TclClass {
- public:
- ExponentialRandomVariableClass() : TclClass("RandomVariable/Exponential") {}
- TclObject* create(int, const char*const*) {
- return(new ExponentialRandomVariable());
- }
- } class_exponentialranvar;
- ExponentialRandomVariable::ExponentialRandomVariable()
- {
- bind("avg_", &avg_);
- }
- ExponentialRandomVariable::ExponentialRandomVariable(double avg)
- {
- avg_ = avg;
- }
- double ExponentialRandomVariable::value()
- {
- return(rng_->exponential(avg_));
- }
- /*
- Generates Erlang variables following:
-
- x^(k-1)*exp(-x/lambda)
- p(x; lambda, k) = ------------------------
- (k-1)!*lambda^k
- */
- static class ErlangRandomVariableClass : public TclClass {
- public:
- ErlangRandomVariableClass() : TclClass("RandomVariable/Erlang") {}
- TclObject* create(int, const char*const*) {
- return(new ErlangRandomVariable());
- }
- } class_erlangranvar;
-
- ErlangRandomVariable::ErlangRandomVariable()
- {
- bind("lambda_", &lambda_);
- bind("k_", &k_);
- }
- ErlangRandomVariable::ErlangRandomVariable(double lambda, int k)
- {
- lambda_ = lambda;
- k_ = k;
- }
- double ErlangRandomVariable::value()
- {
- ExponentialRandomVariable * expRV = new ExponentialRandomVariable(lambda_);
- double result=0;
- for (int i=0; i < k_; i++){
- result += expRV->value();
- }
- return result;
- }
- /*
- Generates Erlang variables following:
-
- x^(alpha-1)*exp(-x/beta)
- p(x; alpha, beta) = --------------------------
- Gamma(alpha)*beta^alpha
- */
- static class GammaRandomVariableClass : public TclClass {
- public:
- GammaRandomVariableClass() : TclClass("RandomVariable/Gamma") {}
- TclObject* create(int, const char*const*) {
- return(new GammaRandomVariable());
- }
- } class_gammaranvar;
- GammaRandomVariable::GammaRandomVariable()
- {
- bind("alpha_", &alpha_);
- bind("beta_", &beta_);
- }
- GammaRandomVariable::GammaRandomVariable(double alpha, double beta)
- {
- alpha_ = alpha;
- beta_ = beta;
- }
- double GammaRandomVariable::value()
- {
- // Proposed by Marsaglia in 2000:
- // G. Marsaglia, W. W. Tsang: A simple method for gereating Gamma variables
- // ACM Transactions on mathematical software, Vol. 26, No. 3, Sept. 2000
- if (alpha_ < 1) {
- double u = rng_->uniform(1.0);
- return GammaRandomVariable::GammaRandomVariable(1.0 + alpha_, beta_).value() * pow (u, 1.0 / alpha_);
- }
-
- double x, v, u;
- double d = alpha_ - 1.0 / 3.0;
- double c = (1.0 / 3.0) / sqrt (d);
- while (1) {
- do {
- x = rng_->normal(0.0, 1.0);
- v = 1.0 + c * x;
- } while (v <= 0);
- v = v * v * v;
- u = rng_->uniform(1.0);
- if (u < 1 - 0.0331 * x * x * x * x)
- break;
- if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
- break;
- }
- return beta_ * d * v;
- }
- static class ParetoRandomVariableClass : public TclClass {
- public:
- ParetoRandomVariableClass() : TclClass("RandomVariable/Pareto") {}
- TclObject* create(int, const char*const*) {
- return(new ParetoRandomVariable());
- }
- } class_paretoranvar;
- ParetoRandomVariable::ParetoRandomVariable()
- {
- bind("avg_", &avg_);
- bind("shape_", &shape_);
- }
- ParetoRandomVariable::ParetoRandomVariable(double avg, double shape)
- {
- avg_ = avg;
- shape_ = shape;
- }
- double ParetoRandomVariable::value()
- {
- /* yuck, user wants to specify shape and avg, but the
- * computation here is simpler if we know the 'scale'
- * parameter. right thing is to probably do away with
- * the use of 'bind' and provide an API such that we
- * can update the scale everytime the user updates shape
- * or avg.
- */
- return(rng_->pareto(avg_ * (shape_ -1)/shape_, shape_));
- }
- /* Pareto distribution of the second kind, aka. Lomax distribution */
- static class ParetoIIRandomVariableClass : public TclClass {
- public:
- ParetoIIRandomVariableClass() : TclClass("RandomVariable/ParetoII") {}
- TclObject* create(int, const char*const*) {
- return(new ParetoIIRandomVariable());
- }
- } class_paretoIIranvar;
- ParetoIIRandomVariable::ParetoIIRandomVariable()
- {
- bind("avg_", &avg_);
- bind("shape_", &shape_);
- }
- ParetoIIRandomVariable::ParetoIIRandomVariable(double avg, double shape)
- {
- avg_ = avg;
- shape_ = shape;
- }
- double ParetoIIRandomVariable::value()
- {
- return(rng_->paretoII(avg_ * (shape_ - 1), shape_));
- }
- static class NormalRandomVariableClass : public TclClass {
- public:
- NormalRandomVariableClass() : TclClass("RandomVariable/Normal") {}
- TclObject* create(int, const char*const*) {
- return(new NormalRandomVariable());
- }
- } class_normalranvar;
-
- NormalRandomVariable::NormalRandomVariable()
- {
- bind("avg_", &avg_);
- bind("std_", &std_);
- }
-
- double NormalRandomVariable::value()
- {
- return(rng_->normal(avg_, std_));
- }
- static class LogNormalRandomVariableClass : public TclClass {
- public:
- LogNormalRandomVariableClass() : TclClass("RandomVariable/LogNormal") {}
- TclObject* create(int, const char*const*) {
- return(new LogNormalRandomVariable());
- }
- } class_lognormalranvar;
-
- LogNormalRandomVariable::LogNormalRandomVariable()
- {
- bind("avg_", &avg_);
- bind("std_", &std_);
- }
-
- double LogNormalRandomVariable::value()
- {
- return(rng_->lognormal(avg_, std_));
- }
- static class ConstantRandomVariableClass : public TclClass {
- public:
- ConstantRandomVariableClass() : TclClass("RandomVariable/Constant"){}
- TclObject* create(int, const char*const*) {
- return(new ConstantRandomVariable());
- }
- } class_constantranvar;
- ConstantRandomVariable::ConstantRandomVariable()
- {
- bind("val_", &val_);
- }
- ConstantRandomVariable::ConstantRandomVariable(double val)
- {
- val_ = val;
- }
- double ConstantRandomVariable::value()
- {
- return(val_);
- }
- /* Hyperexponential distribution code adapted from code provided
- * by Ion Stoica.
- */
- static class HyperExponentialRandomVariableClass : public TclClass {
- public:
- HyperExponentialRandomVariableClass() :
- TclClass("RandomVariable/HyperExponential") {}
- TclObject* create(int, const char*const*) {
- return(new HyperExponentialRandomVariable());
- }
- } class_hyperexponentialranvar;
- HyperExponentialRandomVariable::HyperExponentialRandomVariable()
- {
- bind("avg_", &avg_);
- bind("cov_", &cov_);
- alpha_ = .95;
- }
- HyperExponentialRandomVariable::HyperExponentialRandomVariable(double avg, double cov)
- {
- alpha_ = .95;
- avg_ = avg;
- cov_ = cov;
- }
- double HyperExponentialRandomVariable::value()
- {
- double temp, res;
- double u = Random::uniform();
- temp = sqrt((cov_ * cov_ - 1.0)/(2.0 * alpha_ * (1.0 - alpha_)));
- if (u < alpha_)
- res = Random::exponential(avg_ - temp * (1.0 - alpha_) * avg_);
- else
- res = Random::exponential(avg_ + temp * (alpha_) * avg_);
- return(res);
- }
- static class WeibullRandomVariableClass : public TclClass {
- public:
- WeibullRandomVariableClass() : TclClass("RandomVariable/Weibull") {}
- TclObject* create(int argc, const char*const* argv) {
- if (argc == 6) {
- return (new WeibullRandomVariable ((double) atof(argv[4]),
- (double) atof(argv[5])));
- }
- else if (argc == 7) {
- // shape, scale, RNG
- RNG* rng = (RNG*)TclObject::lookup(argv[6]);
- return (new WeibullRandomVariable ((double) atof(argv[4]),
- (double) atof(argv[5]),
- rng));
- }
- else {
- return(new WeibullRandomVariable());
- }
- }
- } class_weibullranvar;
- WeibullRandomVariable::WeibullRandomVariable()
- {
- bind("shape_", &shape_);
- bind("scale_", &scale_);
- }
- WeibullRandomVariable::WeibullRandomVariable(double scale, double shape)
- {
- shape_ = shape;
- scale_ = scale;
- }
- WeibullRandomVariable::WeibullRandomVariable(double scale, double shape,
- RNG* rng)
- {
- shape_ = shape;
- scale_ = scale;
- rng_ = rng;
- }
- double WeibullRandomVariable::avg(void)
- {
- return 0;
- }
- double WeibullRandomVariable::value()
- {
- return(rng_->rweibull(scale_, shape_));
- }
-
- /*
- // Empirical Random Variable:
- // CDF input from file with the following column
- // 1. Possible values in a distrubutions
- // 2. Number of occurances for those values
- // 3. The CDF for those value
- // code provided by Giao Nguyen
- */
- static class EmpiricalRandomVariableClass : public TclClass {
- public:
- EmpiricalRandomVariableClass() : TclClass("RandomVariable/Empirical"){}
- TclObject* create(int, const char*const*) {
- return(new EmpiricalRandomVariable());
- }
- } class_empiricalranvar;
- EmpiricalRandomVariable::EmpiricalRandomVariable() : minCDF_(0), maxCDF_(1), maxEntry_(32), table_(0)
- {
- bind("minCDF_", &minCDF_);
- bind("maxCDF_", &maxCDF_);
- bind("interpolation_", &interpolation_);
- bind("maxEntry_", &maxEntry_);
- }
- int EmpiricalRandomVariable::command(int argc, const char*const* argv)
- {
- Tcl& tcl = Tcl::instance();
- if (argc == 3) {
- if (strcmp(argv[1], "loadCDF") == 0) {
- if (loadCDF(argv[2]) == 0) {
- tcl.resultf("%s loadCDF %s: invalid file",
- name(), argv[2]);
- return (TCL_ERROR);
- }
- return (TCL_OK);
- }
- }
- return RandomVariable::command(argc, argv);
- }
- int EmpiricalRandomVariable::loadCDF(const char* filename)
- {
- FILE* fp;
- char line[256];
- CDFentry* e;
- fp = fopen(filename, "r");
- if (fp == 0)
- return 0;
- if (table_ == 0)
- table_ = new CDFentry[maxEntry_];
- for (numEntry_=0; fgets(line, 256, fp); numEntry_++) {
- if (numEntry_ >= maxEntry_) { // resize the CDF table
- maxEntry_ *= 2;
- e = new CDFentry[maxEntry_];
- for (int i=numEntry_-1; i >= 0; i--)
- e[i] = table_[i];
- delete table_;
- table_ = e;
- }
- e = &table_[numEntry_];
- // Use * and l together raises a warning
- sscanf(line, "%lf %*f %lf", &e->val_, &e->cdf_);
- }
- fclose(fp);
- return numEntry_;
- }
- double EmpiricalRandomVariable::value()
- {
- if (numEntry_ <= 0)
- return 0;
- double u = rng_->uniform(minCDF_, maxCDF_);
- int mid = lookup(u);
- if (mid && interpolation_ && u < table_[mid].cdf_)
- return interpolate(u, table_[mid-1].cdf_, table_[mid-1].val_,
- table_[mid].cdf_, table_[mid].val_);
- return table_[mid].val_;
- }
- double EmpiricalRandomVariable::interpolate(double x, double x1, double y1, double x2, double y2)
- {
- double value = y1 + (x - x1) * (y2 - y1) / (x2 - x1);
- if (interpolation_ == INTER_INTEGRAL) // round up
- return ceil(value);
- return value;
- }
- int EmpiricalRandomVariable::lookup(double u)
- {
- // always return an index whose value is >= u
- int lo, hi, mid;
- if (u <= table_[0].cdf_)
- return 0;
- for (lo=1, hi=numEntry_-1; lo < hi; ) {
- mid = (lo + hi) / 2;
- if (u > table_[mid].cdf_)
- lo = mid + 1;
- else
- hi = mid;
- }
- return lo;
- }