#include <stdio.h>#include <math.h>#include <stdlib.h>#include <string.h>#include "ranv2.h"#include "ranvar.h"#include "rng.h"#include "debug.h"const double Y_ALPHA = 3.10;  /* 0.999 quantile of normal distribution */const double HPD = 1.2564708; /* number of poles (and zeros) per decade for flicker filter ( = 1/log10(2.5^2) ); agrees with Barnes, Jarvis, and Greenhall */void init_flick (int n, double h, double *a, double *b, double gamma);void flicker (int N, double *y, double *y_old, double *a, double *b);/* Implement a Self-Similar traffic source.  The batch arrival * time process is self-similar with Hurst Exponenet H and * specified mean arrival rate.  The times * between batches are Gamma Distributed (and correlated, in * accordance with H).  The batch size process is self-similar, * with Hurst Exponenet (not necessarily the same as the Hurst * Exponent for the arrival times).  The batch size distribution * is negative binomial. */class SelfSimilar_Traffic : public TrafficGenerator { public:	SelfSimilar_Traffic();	~SelfSimilar_Traffic();	virtual double next_interval(int&);	double next_batch_size();        virtual void timeout();	double avg_rate_;     /* mean arrival rate of traffic batchs */        double sdt_;  /* std deviation of time between batches */        double avg_batch_size_; /* mean batch size */        int sb_;  /* neg binomial distribution variance parameter (number of convolved geometric distributions) */        double Ht_;  /* Hurst Exponent for batch arrival process */        double Hb_;  /* Hurst Exponent for batch size */        double starttime_; /* time (s) that this self-similar traffic source is turned on */        double stoptime_; /* time (s) that this self-similar traffic source is turned off */ protected:        void init();        GammaRandomVariable interval_; /* batch inter-arrival time (sec) */        NegBinomRandomVariable batch_size_; /* average batch size (bytes)*/ private:        int Nt;  /* order of flicker filters */        int nt;  /* number of decades the flicker filters will simulate */        /* Note:  both the batch arrival and batch size flicker filters have the same order, because this is based on the total number of batches produced */        /* arrival process parameters and states */        double ht;  /* number of poles per decade for arrival process flicker filter (may be a fraction) */        double *at;  /* pointer to array of poles for arrival process flicker filter */        double *bt;  /* pointer to array of zeros for arrival process flicker filter */        double gamt;  /* power spectral density noise exponent for arrival process flicker filter */        double gamt0;  /* equals gamt for 0 <= gamt < 2.0; else gamt - 2.0; for gamt > 2, we run the flicker filter for gamt0 and then integrate */        double *xt;  /* pointer to current state vector array plus white noise input in element zero (array size = Nt+1)*/        double *xt_old;  /* pointer to state vector array plus white noise input in element zero at previous iteration (array size = Nt+1)*/        double yt;  /* current flicker noise value for arrival time process */        double yt_old;  /* previous iteration flicker noise value for arrival time process */        /* batch size process parameters and states */        double hb;  /* number of poles per decade for batch size process flicker filter (may be a fraction) */        double *ab;  /* pointer to array of poles for batch size process flicker filter */        double *bb;  /* pointer to array of zeros for batch size process flicker filter */        double gamb;  /* power spectral density noise exponent for batch size process flicker filter */        double gamb0;  /* equals gamb for 0 <= gamb < 2.0; else gamb - 2.0; for gamb > 2, we run the flicker filter for gamb0 and then integrate */        double *xb;  /* pointer to current state vector array plus white noise input in element zero (array size = Nt+1)*/        double *xb_old;  /* pointer to state vector array plus white noise input in element zero at previous iteration (array size = Nt+1)*/        double yb;  /* current flicker noise value for batch size process */        double yb_old;  /* previous iteration flicker noise value for batch size process */};static class SelfSimilarTrafficClass : public TclClass { public:	SelfSimilarTrafficClass() : TclClass("Application/Traffic/SelfSimilar") {}	TclObject* create(int, const char*const*) {		return (new SelfSimilar_Traffic());	}} class_selfsimilar_traffic;SelfSimilar_Traffic::SelfSimilar_Traffic() : interval_(0.0, 0.0), batch_size_(0.0, 0){     double alphat, betat;     double tot_time;     int nd0;     double msr;     char s[250];     bind_bw("rate", &avg_rate_);     bind_bw("std_dev_inter_batch_time", &sdt_);     bind("batchsize", &avg_batch_size_);     bind("sb", &sb_);     bind("Ht", &Ht_);     bind("Hb", &Hb_);     bind("starttime", &starttime_);     bind("stoptime", &stoptime_);     msr = 1.0/(avg_rate_ * sdt_);     if (msr > 25.0)     {        sprintf (s, "Error -- ratio of mean to std dev for self-similar process exceeds 25; it would be better to use a CBR traffic model");        Debug::debug(s);        sprintf (s, "avg_rate = %f   mean interarrival time = %f   std dev interarrival time = %f   ratio = %f", avg_rate_, 1.0/avg_rate_, sdt_, msr);        Debug::debug(s);        exit(1);     }/* * Calculate number of decades that flicker filter must be valid * over * as 0.999 quantile of number of cells that will arrive in * simulation * time ; assume the distribution is normal (asymptotic) and * interarrival times are independent (renewal process).  These * are * approximations, as the interarrival times for a flicker process * are correlated; but note that here we are only determining the * filter order.  Also calculate number of poles and zeros (filter * stages). */     tot_time = starttime_ - stoptime_;     if (sdt_ > 0.0)     {        alphat = 1.0 / (avg_rate_ * avg_rate_ * sdt_ * sdt_);        betat = sdt_ * sdt_ * avg_rate_;        nd0 = log10 ( ceil (Y_ALPHA *                      sqrt (sdt_*sdt_*avg_rate_*avg_rate_*                           avg_rate_*tot_time) +                              avg_rate_*tot_time));        if (nd0 < 8) /* set minimum flicker filter order to 8; this                   * is what was used in previous simulations of                   * flicker noise in clocks */           nd0 = 8;        nt = nd0;     }     else     {        alphat = 0.0;        betat = 0.0;        nt = nd0 = 8;     }     Nt = ceil (nd0*HPD);/* * Allocate arrays for parameters and states for both flicker * filters */     at = new double[Nt];     bt = new double[Nt];     xt = new double[Nt+1];     xt_old = new double[Nt+1];     ab = new double[Nt];     bb = new double[Nt];     xb = new double[Nt+1];     xb_old = new double[Nt+1];/* * Initialize fractional flicker filter states and outputs to zero */     memset (xt, 0, (Nt+1)*sizeof(double) );     memset (xt_old, 0, (Nt+1)*sizeof(double) );     memset (xb, 0, (Nt+1)*sizeof(double) );     memset (xb_old, 0, (Nt+1)*sizeof(double) );     yt = yt_old = yb = yb_old = 0.0;}SelfSimilar_Traffic::~SelfSimilar_Traffic(){/* * Free arrays for parameters and states for both flicker * filters */     delete [] at;     delete [] bt;     delete [] xt;     delete [] xt_old;     delete [] ab;     delete [] bb;     delete [] xb;     delete [] xb_old;}void SelfSimilar_Traffic::init(){/* * Calculate PSD exponent for arrival and batch size processes, * from Hurst Exponents */     gamt = 2*Ht_ + 1;     if (gamt < 2.0)        gamt0 = gamt;     else        gamt0 = gamt - 2.0;     gamb = 2*Hb_ + 1;     if (gamb < 2.0)        gamb0 = gamb;     else        gamb0 = gamb - 2.0;/* * Initialize mean and std dev for Gamma and Neg Binomial random * variables */     interval_.setavg( 1.0/avg_rate_);     interval_.setstdev( sdt_);     batch_size_.setavg( avg_batch_size_);     batch_size_.setsparm( sb_);     /* * Invoke init_flick to initialized flicker filter parameters */     init_flick (nt, HPD, at, bt, gamt0);     init_flick (nt, HPD, ab, bb, gamb0);/* * Set packet type */     if (agent_)        if (agent_->get_pkttype() != PT_TCP &&                             agent_->get_pkttype() != PT_TFRC)           agent_->set_pkttype(PT_SELFSIM);}double SelfSimilar_Traffic::next_interval(int& size){     double val;     int icount0 = 1;     char s[250];     if (sdt_ > 0.0)     {        xt[0] = interval_.value() - interval_.avg();        if (gamt0 > 0.0)           flicker (Nt, xt, xt_old, at, bt);        else //white noise case (note that gamt0 < 0 is not allowed)           xt[Nt] = xt[0];        if (gamt >= 2.0) /* Integrate the noise process for                 * PSD exponent of 2 or greater (i.e., integrate                 * Gaussian noise motion to get fractional                 * Brownian motion) */           yt = yt_old + xt[Nt];        else           yt = xt[Nt];     }     else        yt = yt_old = 0.0;// Note that val is the time we add to the current time to get// the next time.  Since the filter output is the offset from// the nominal (average) time, we must add the average to get//  the next nominal time, and then subtract the previous// offset and add the current offset     val = yt - yt_old + interval_.avg();// If val <= 0, print error message; batch interarrival time must be > 0//          Then exit.     if (val <= 0.0)     {        sprintf (s, "Error - self-similar traffic model produced non-positive batch interarrival time.  Std deviation is likely too large relative to the mean interarrival time.");        Debug::debug(s);        sprintf (s, "Set stdev to <= 0.1 times mean.");        Debug::debug(s);        sprintf (s, "Ht = %f     rate = %f     mean interarriv time = %f    std dev interarr time = %f",Ht_, avg_rate_, interval_.avg(), interval_.stdev());        Debug::debug(s);        exit (1);     }     yt_old = yt;     return (val);}double SelfSimilar_Traffic::next_batch_size(){     double val;     if (sb_ > 0)     {        xb[0] = batch_size_.value() - batch_size_.avg();        flicker (Nt, xb, xb_old, ab, bb);        if (gamb >= 2.0) /* Integrate the noise process for                 * PSD exponent of 2 or greater (i.e., integrate                 * Gaussian noise motion to get fractional                 * Brownian motion) */           yb = yb_old + xb[Nt];        else           yb = xb[Nt];     }     else        yb = yb_old = 0.0;// Note that here we simply return batch size; there is no// need to subtract the previous offset as was done for// arrival times.     val = yb + batch_size_.avg();     yb_old = yb;     if (val < 1.0) // for batch size process, we enforce a hard lower limit                    // of 1.0; no need to report if < 1.        val = 1.0;     return (val);}void SelfSimilar_Traffic::timeout(){     if (! running_)        return;     /* send a packet */     size_ = next_batch_size();     agent_->sendmsg(size_);     /* figure out when to send the next one */     nextPkttime_ = next_interval(size_);     double tnextt = nextPkttime_;     /* schedule it */     if (nextPkttime_ > 0)     timer_.resched(nextPkttime_);}