ranvar.cc
上传用户:rrhhcc
上传日期:2015-12-11
资源大小:54129k
文件大小:15k
源码类别:

通讯编程

开发平台:

Visual C++

  1. /* -*- Mode:C++; c-basic-offset:8; tab-width:8; indent-tabs-mode:t -*- */
  2. /*
  3.  * Copyright (c) Xerox Corporation 1997. All rights reserved.
  4.  *
  5.  * This program is free software; you can redistribute it and/or modify it
  6.  * under the terms of the GNU General Public License as published by the
  7.  * Free Software Foundation; either version 2 of the License, or (at your
  8.  * option) any later version.
  9.  *
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  *
  15.  * You should have received a copy of the GNU General Public License along
  16.  * with this program; if not, write to the Free Software Foundation, Inc.,
  17.  * 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  18.  *
  19.  * Linking this file statically or dynamically with other modules is making
  20.  * a combined work based on this file.  Thus, the terms and conditions of
  21.  * the GNU General Public License cover the whole combination.
  22.  *
  23.  * In addition, as a special exception, the copyright holders of this file
  24.  * give you permission to combine this file with free software programs or
  25.  * libraries that are released under the GNU LGPL and with code included in
  26.  * the standard release of ns-2 under the Apache 2.0 license or under
  27.  * otherwise-compatible licenses with advertising requirements (or modified
  28.  * versions of such code, with unchanged license).  You may copy and
  29.  * distribute such a system following the terms of the GNU GPL for this
  30.  * file and the licenses of the other code concerned, provided that you
  31.  * include the source code of that other code when and as the GNU GPL
  32.  * requires distribution of source code.
  33.  *
  34.  * Note that people who make modified versions of this file are not
  35.  * obligated to grant this special exception for their modified versions;
  36.  * it is their choice whether to do so.  The GNU General Public License
  37.  * gives permission to release a modified version without this exception;
  38.  * this exception also makes it possible to release a modified version
  39.  * which carries forward this exception.
  40.  */
  41. #ifndef lint
  42. static const char rcsid[] =
  43.     "@(#) $Header: /cvsroot/nsnam/ns-2/tools/ranvar.cc,v 1.22 2008/02/01 21:39:43 tom_henderson Exp $ (Xerox)";
  44. #endif
  45. #include <stdio.h>
  46. #include "ranvar.h"
  47. RandomVariable::RandomVariable()
  48. {
  49. rng_ = RNG::defaultrng(); 
  50. }
  51. int RandomVariable::command(int argc, const char*const* argv)
  52. {
  53. Tcl& tcl = Tcl::instance();
  54. if (argc == 2) {
  55. if (strcmp(argv[1], "value") == 0) {
  56. tcl.resultf("%6e", value());
  57. return(TCL_OK);
  58. }
  59. }
  60. if (argc == 3) {
  61. if (strcmp(argv[1], "use-rng") == 0) {
  62. rng_ = (RNG*)TclObject::lookup(argv[2]);
  63. if (rng_ == 0) {
  64. tcl.resultf("no such RNG %s", argv[2]);
  65. return(TCL_ERROR);
  66. }
  67. return(TCL_OK);
  68. }
  69. }
  70. return(TclObject::command(argc, argv));
  71. }
  72. // Added by Debojyoti Dutta 12 October 2000
  73. // This allows us to seed a randomvariable with 
  74. // our own RNG object. This command is called from 
  75. // expoo.cc and pareto.cc 
  76. int  RandomVariable::seed(char *x){
  77.         
  78.         Tcl& tcl = Tcl::instance();
  79.                 rng_ = (RNG*)TclObject::lookup(x);
  80.                 if (rng_ == 0) {
  81.                         tcl.resultf("no such RNG %s", x);
  82.                         return(TCL_ERROR);
  83.                 }
  84.                 return(TCL_OK);
  85.  
  86. }
  87. static class UniformRandomVariableClass : public TclClass {
  88. public:
  89. UniformRandomVariableClass() : TclClass("RandomVariable/Uniform"){}
  90. TclObject* create(int, const char*const*) {
  91.  return(new UniformRandomVariable());
  92.  }
  93. } class_uniformranvar;
  94. UniformRandomVariable::UniformRandomVariable()
  95. {
  96. bind("min_", &min_);
  97. bind("max_", &max_); 
  98. }
  99. UniformRandomVariable::UniformRandomVariable(double min, double max)
  100. {
  101. min_ = min;
  102. max_ = max;
  103. }
  104. double UniformRandomVariable::value()
  105. {
  106. return(rng_->uniform(min_, max_));
  107. }
  108. static class ExponentialRandomVariableClass : public TclClass {
  109. public:
  110. ExponentialRandomVariableClass() : TclClass("RandomVariable/Exponential") {}
  111. TclObject* create(int, const char*const*) {
  112. return(new ExponentialRandomVariable());
  113. }
  114. } class_exponentialranvar;
  115. ExponentialRandomVariable::ExponentialRandomVariable()
  116. {
  117. bind("avg_", &avg_);
  118. }
  119. ExponentialRandomVariable::ExponentialRandomVariable(double avg)
  120. {
  121. avg_ = avg;
  122. }
  123. double ExponentialRandomVariable::value()
  124. {
  125. return(rng_->exponential(avg_));
  126. }
  127. /*
  128. Generates Erlang variables following:
  129.                    x^(k-1)*exp(-x/lambda)
  130. p(x; lambda, k) = ------------------------
  131.                        (k-1)!*lambda^k
  132. */
  133. static class ErlangRandomVariableClass : public TclClass {
  134. public:
  135. ErlangRandomVariableClass() : TclClass("RandomVariable/Erlang") {}
  136. TclObject* create(int, const char*const*) {
  137. return(new ErlangRandomVariable());
  138. }
  139. } class_erlangranvar;
  140.  
  141. ErlangRandomVariable::ErlangRandomVariable()
  142. {
  143. bind("lambda_", &lambda_);
  144. bind("k_", &k_);
  145. }
  146. ErlangRandomVariable::ErlangRandomVariable(double lambda, int k)
  147. {
  148. lambda_ = lambda;
  149. k_ = k;
  150. }
  151. double ErlangRandomVariable::value()
  152. {
  153. ExponentialRandomVariable * expRV = new ExponentialRandomVariable(lambda_);
  154. double result=0;
  155. for (int i=0; i < k_; i++){
  156. result += expRV->value();
  157. }
  158. return result;
  159. }
  160. /*
  161. Generates Erlang variables following:
  162.                      x^(alpha-1)*exp(-x/beta)
  163. p(x; alpha, beta) = --------------------------
  164.                          Gamma(alpha)*beta^alpha
  165. */
  166. static class GammaRandomVariableClass : public TclClass {
  167. public:
  168. GammaRandomVariableClass() : TclClass("RandomVariable/Gamma") {}
  169. TclObject* create(int, const char*const*) {
  170. return(new GammaRandomVariable());
  171. }
  172. } class_gammaranvar;
  173. GammaRandomVariable::GammaRandomVariable()
  174. {
  175. bind("alpha_", &alpha_);
  176. bind("beta_", &beta_);
  177. }
  178. GammaRandomVariable::GammaRandomVariable(double alpha, double beta)
  179. {
  180. alpha_ = alpha;
  181. beta_ = beta;
  182. }
  183. double GammaRandomVariable::value()
  184. {
  185. // Proposed by Marsaglia in 2000:
  186. // G. Marsaglia, W. W. Tsang: A simple method for gereating Gamma variables
  187. // ACM Transactions on mathematical software, Vol. 26, No. 3, Sept. 2000
  188. if (alpha_ < 1) {
  189. double u = rng_->uniform(1.0);
  190. return GammaRandomVariable::GammaRandomVariable(1.0 + alpha_, beta_).value() * pow (u, 1.0 / alpha_);
  191. }
  192. double x, v, u;
  193. double d = alpha_ - 1.0 / 3.0;
  194. double c = (1.0 / 3.0) / sqrt (d);
  195. while (1) {
  196. do {
  197. x = rng_->normal(0.0, 1.0);
  198. v = 1.0 + c * x;
  199. } while (v <= 0);
  200. v = v * v * v;
  201. u = rng_->uniform(1.0);
  202. if (u < 1 - 0.0331 * x * x * x * x)
  203. break;
  204. if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
  205. break;
  206. }
  207. return beta_ * d * v;
  208. }
  209. static class ParetoRandomVariableClass : public TclClass {
  210.  public:
  211. ParetoRandomVariableClass() : TclClass("RandomVariable/Pareto") {}
  212. TclObject* create(int, const char*const*) {
  213. return(new ParetoRandomVariable());
  214. }
  215. } class_paretoranvar;
  216. ParetoRandomVariable::ParetoRandomVariable()
  217. {
  218. bind("avg_", &avg_);
  219. bind("shape_", &shape_);
  220. }
  221. ParetoRandomVariable::ParetoRandomVariable(double avg, double shape)
  222. {
  223. avg_ = avg;
  224. shape_ = shape;
  225. }
  226. double ParetoRandomVariable::value()
  227. {
  228. /* yuck, user wants to specify shape and avg, but the
  229.  * computation here is simpler if we know the 'scale'
  230.  * parameter.  right thing is to probably do away with
  231.  * the use of 'bind' and provide an API such that we
  232.  * can update the scale everytime the user updates shape
  233.  * or avg.
  234.  */
  235. return(rng_->pareto(avg_ * (shape_ -1)/shape_, shape_));
  236. }
  237. /* Pareto distribution of the second kind, aka. Lomax distribution */
  238. static class ParetoIIRandomVariableClass : public TclClass {
  239.  public:
  240.         ParetoIIRandomVariableClass() : TclClass("RandomVariable/ParetoII") {}
  241.         TclObject* create(int, const char*const*) {
  242.                 return(new ParetoIIRandomVariable());
  243.         }
  244. } class_paretoIIranvar;
  245. ParetoIIRandomVariable::ParetoIIRandomVariable()
  246. {
  247.         bind("avg_", &avg_);
  248.         bind("shape_", &shape_);
  249. }
  250. ParetoIIRandomVariable::ParetoIIRandomVariable(double avg, double shape)
  251. {
  252.         avg_ = avg;
  253.         shape_ = shape;
  254. }
  255. double ParetoIIRandomVariable::value()
  256. {
  257.         return(rng_->paretoII(avg_ * (shape_ - 1), shape_));
  258. }
  259. static class NormalRandomVariableClass : public TclClass {
  260.  public:
  261.         NormalRandomVariableClass() : TclClass("RandomVariable/Normal") {}
  262.         TclObject* create(int, const char*const*) {
  263.                 return(new NormalRandomVariable());
  264.         }
  265. } class_normalranvar;
  266.  
  267. NormalRandomVariable::NormalRandomVariable()
  268.         bind("avg_", &avg_);
  269.         bind("std_", &std_);
  270. }
  271.  
  272. double NormalRandomVariable::value()
  273. {
  274.         return(rng_->normal(avg_, std_));
  275. }
  276. static class LogNormalRandomVariableClass : public TclClass {
  277.  public:
  278.         LogNormalRandomVariableClass() : TclClass("RandomVariable/LogNormal") {}
  279.         TclObject* create(int, const char*const*) {
  280.                 return(new LogNormalRandomVariable());
  281.         }
  282. } class_lognormalranvar;
  283.  
  284. LogNormalRandomVariable::LogNormalRandomVariable()
  285.         bind("avg_", &avg_);
  286.         bind("std_", &std_);
  287. }
  288.  
  289. double LogNormalRandomVariable::value()
  290. {
  291.         return(rng_->lognormal(avg_, std_));
  292. }
  293. static class ConstantRandomVariableClass : public TclClass {
  294.  public:
  295. ConstantRandomVariableClass() : TclClass("RandomVariable/Constant"){}
  296. TclObject* create(int, const char*const*) {
  297. return(new ConstantRandomVariable());
  298. }
  299. } class_constantranvar;
  300. ConstantRandomVariable::ConstantRandomVariable()
  301. {
  302. bind("val_", &val_);
  303. }
  304. ConstantRandomVariable::ConstantRandomVariable(double val)
  305. {
  306. val_ = val;
  307. }
  308. double ConstantRandomVariable::value()
  309. {
  310. return(val_);
  311. }
  312. /* Hyperexponential distribution code adapted from code provided
  313.  * by Ion Stoica.
  314.  */
  315. static class HyperExponentialRandomVariableClass : public TclClass {
  316. public:
  317. HyperExponentialRandomVariableClass() : 
  318. TclClass("RandomVariable/HyperExponential") {}
  319. TclObject* create(int, const char*const*) {
  320. return(new HyperExponentialRandomVariable());
  321. }
  322. } class_hyperexponentialranvar;
  323. HyperExponentialRandomVariable::HyperExponentialRandomVariable()
  324. {
  325. bind("avg_", &avg_);
  326. bind("cov_", &cov_);
  327. alpha_ = .95;
  328. }
  329. HyperExponentialRandomVariable::HyperExponentialRandomVariable(double avg, double cov)
  330. {
  331. alpha_ = .95;
  332. avg_ = avg;
  333. cov_ = cov;
  334. }
  335. double HyperExponentialRandomVariable::value()
  336. {
  337. double temp, res;
  338. double u = Random::uniform();
  339. temp = sqrt((cov_ * cov_ - 1.0)/(2.0 * alpha_ * (1.0 - alpha_)));
  340. if (u < alpha_)
  341. res = Random::exponential(avg_ - temp * (1.0 - alpha_) * avg_);
  342. else
  343. res = Random::exponential(avg_ + temp * (alpha_) * avg_);
  344. return(res);
  345. }
  346. static class WeibullRandomVariableClass : public TclClass {
  347.  public:
  348.         WeibullRandomVariableClass() : TclClass("RandomVariable/Weibull") {}
  349.         TclObject* create(int argc, const char*const* argv) {
  350.                 if (argc == 6) {
  351.                         return (new WeibullRandomVariable ((double) atof(argv[4]),
  352.                                                            (double) atof(argv[5])));
  353.                 }
  354.                 else if (argc == 7) {
  355.                         // shape, scale, RNG
  356.                         RNG* rng = (RNG*)TclObject::lookup(argv[6]);
  357.                         return (new WeibullRandomVariable ((double) atof(argv[4]),
  358.                                                            (double) atof(argv[5]),
  359.                                                            rng));
  360.                 }
  361.                 else {
  362.                         return(new WeibullRandomVariable());
  363.                 }
  364.         }                                           
  365. } class_weibullranvar;
  366. WeibullRandomVariable::WeibullRandomVariable()
  367. {        
  368. bind("shape_", &shape_);
  369.         bind("scale_", &scale_);
  370. }
  371. WeibullRandomVariable::WeibullRandomVariable(double scale, double shape)
  372. {
  373.         shape_ = shape;
  374.         scale_ = scale;
  375. }
  376. WeibullRandomVariable::WeibullRandomVariable(double scale, double shape, 
  377.      RNG* rng) 
  378. {
  379.         shape_ = shape;
  380.         scale_ = scale;
  381.         rng_ = rng;
  382. }
  383. double WeibullRandomVariable::avg(void)         
  384. {
  385.         return 0;
  386. }
  387. double WeibullRandomVariable::value()
  388. {
  389.         return(rng_->rweibull(scale_, shape_));
  390. }
  391.                                                                                
  392. /*
  393. // Empirical Random Variable:
  394. //  CDF input from file with the following column
  395. //   1.  Possible values in a distrubutions
  396. //   2.  Number of occurances for those values
  397. //   3.  The CDF for those value
  398. //  code provided by Giao Nguyen
  399. */
  400. static class EmpiricalRandomVariableClass : public TclClass {
  401. public:
  402. EmpiricalRandomVariableClass() : TclClass("RandomVariable/Empirical"){}
  403. TclObject* create(int, const char*const*) {
  404. return(new EmpiricalRandomVariable());
  405. }
  406. } class_empiricalranvar;
  407. EmpiricalRandomVariable::EmpiricalRandomVariable() : minCDF_(0), maxCDF_(1), maxEntry_(32), table_(0)
  408. {
  409. bind("minCDF_", &minCDF_);
  410. bind("maxCDF_", &maxCDF_);
  411. bind("interpolation_", &interpolation_);
  412. bind("maxEntry_", &maxEntry_);
  413. }
  414. int EmpiricalRandomVariable::command(int argc, const char*const* argv)
  415. {
  416. Tcl& tcl = Tcl::instance();
  417. if (argc == 3) {
  418. if (strcmp(argv[1], "loadCDF") == 0) {
  419. if (loadCDF(argv[2]) == 0) {
  420. tcl.resultf("%s loadCDF %s: invalid file",
  421.     name(), argv[2]);
  422. return (TCL_ERROR);
  423. }
  424. return (TCL_OK);
  425. }
  426. }
  427. return RandomVariable::command(argc, argv);
  428. }
  429. int EmpiricalRandomVariable::loadCDF(const char* filename)
  430. {
  431. FILE* fp;
  432. char line[256];
  433. CDFentry* e;
  434. fp = fopen(filename, "r");
  435. if (fp == 0) 
  436. return 0;
  437. if (table_ == 0)
  438. table_ = new CDFentry[maxEntry_];
  439. for (numEntry_=0;  fgets(line, 256, fp);  numEntry_++) {
  440. if (numEntry_ >= maxEntry_) { // resize the CDF table
  441. maxEntry_ *= 2;
  442. e = new CDFentry[maxEntry_];
  443. for (int i=numEntry_-1; i >= 0; i--)
  444. e[i] = table_[i];
  445. delete table_;
  446. table_ = e;
  447. }
  448. e = &table_[numEntry_];
  449. // Use * and l together raises a warning
  450. sscanf(line, "%lf %*f %lf", &e->val_, &e->cdf_);
  451. }
  452.         fclose(fp);
  453. return numEntry_;
  454. }
  455. double EmpiricalRandomVariable::value()
  456. {
  457. if (numEntry_ <= 0)
  458. return 0;
  459. double u = rng_->uniform(minCDF_, maxCDF_);
  460. int mid = lookup(u);
  461. if (mid && interpolation_ && u < table_[mid].cdf_)
  462. return interpolate(u, table_[mid-1].cdf_, table_[mid-1].val_,
  463.    table_[mid].cdf_, table_[mid].val_);
  464. return table_[mid].val_;
  465. }
  466. double EmpiricalRandomVariable::interpolate(double x, double x1, double y1, double x2, double y2)
  467. {
  468. double value = y1 + (x - x1) * (y2 - y1) / (x2 - x1);
  469. if (interpolation_ == INTER_INTEGRAL) // round up
  470. return ceil(value);
  471. return value;
  472. }
  473. int EmpiricalRandomVariable::lookup(double u)
  474. {
  475. // always return an index whose value is >= u
  476. int lo, hi, mid;
  477. if (u <= table_[0].cdf_)
  478. return 0;
  479. for (lo=1, hi=numEntry_-1;  lo < hi; ) {
  480. mid = (lo + hi) / 2;
  481. if (u > table_[mid].cdf_)
  482. lo = mid + 1;
  483. else
  484. hi = mid;
  485. }
  486. return lo;
  487. }