//================================================== file = syntraf1a.c =====
//=  Program to generate synthetic packet traffic                           =
//=   - Bounded Pareto burst size and exponential interburst time           =
//=   - Modification of syntraf1.c to take target utilization as an input   =
//===========================================================================
//=  Notes: 1) Writes to a user specified output file                       =
//=             - File format is <interarrival_time packet_size>            =
//=         2) Bounded Pareto pdf is (k is lower bound, p is upper bound):  =
//=             - pdf is f(x) = ((a*k^a) / (1 - (k/p)^a))*x^(-a-1)          =
//=         3) See M. Crovella and M. Harchol-Balter, and C. Murta, "Task   =
//=            Assignment in a Distributed System: Improving Performance    =
//=            by Unbalancing Load," BUCS-TR-1997-018, October 1997.        =
//=-------------------------------------------------------------------------=
//= Example user input:                                                     =
//=                                                                         =
//=   ---------------------------------------------- syntraf1a.c -----      =
//=   -  Program to generate synthetic packet traffic                -      =
//=   -   - Bounded Pareto burst size, exponential interburst time   -      =
//=   ----------------------------------------------------------------      =
//=   Enter output file name =========================> output.dat          =
//=   Random number seed (greater than 0) ============> 1                   =
//=   Target utilizaiton (%) =========================> 1.0                 =
//=   Pareto alpha value =============================> 1.001               =
//=   Pareto k value (min burst size in bytes) =======> 1024.0              =
//=   Pareto p value (max burst size in bytes) =======> 1073741824.0        =
//=   Min packet length (bytes) ======================> 64                  =
//=   Max packet length (bytes) ======================> 1500                =
//=   Link data rate (bits per second) ===============> 1000000000.0        =
//=   Number of packets to generate ==================> 10000000.0          =
//=   ----------------------------------------------------------------      =
//=   ----------------------------------------------------------------      =
//=   -  Generating 10000000 packets...                                     =
//=   -   - Mean burst size     = 14111.821582 bytes                        =
//=   -   - Variance burst size = 1086346073616.395752 bytes                =
//=   -   - Mean intergap time  = 0.011177 sec                              =
//=   ----------------------------------------------------------------      =
//=   ----------------------------------------------------------------      =
//=   -  Summary statistics...                                              =
//=   -   - Total time of trace    = 11150.000292 sec                       =
//=   -   - Minimum size burst     =     1024 bytes                         =
//=   -   - Maximum size burst     = 1015481368 bytes                       =
//=   -   - Total number of bursts =   988343                               =
//=   -   - Total number of bytes  = 14333409182 bytes                      =
//=   -   - Utilization            = 1.028406 %                             =
//=   ----------------------------------------------------------------      =
//=-------------------------------------------------------------------------=
//= Example output file (first 10 entries of "output.dat" for above):       =
//=                                                                         =
//=   0.022679   1025                                                       =
//=   0.008724   1500                                                       =
//=   0.000012   1500                                                       =
//=   0.000009   1185                                                       =
//=   0.016988   1500                                                       =
//=   0.000006   690                                                        =
//=   0.004338   1075                                                       =
//=   0.000767   1500                                                       =
//=   0.000012   1500                                                       =
//=   0.000002   190                                                        =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 syntraf1a.c                                               =
//=-------------------------------------------------------------------------=
//=  Execute: syntraf1a                                                     =
//=-------------------------------------------------------------------------=
//=  Author: Ken Christensen                                                =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (10/11/05) - Genesis (from syntraf1.c)                    =
//===========================================================================
//----- Include files -------------------------------------------------------
#include <stdio.h>              // Needed for printf()
#include <stdlib.h>             // Needed for exit() and ato*()
#include <math.h>               // Needed for log() and pow()
#include <assert.h>             // Needed for assert() macro

//----- Function prototypes -------------------------------------------------
double bpareto(double a, double k, double p); // Returns a bounded Pareto rv
double expon(double mean_value);              // Returns an exponential rv
double rand_val(int seed);                    // Jain's RNG

//===== Main program ========================================================
void main(void)
{
  char     in_string[80];       // Input string
  FILE     *fp;                 // File pointer to output file
  double   target_util;         // Target utilization (%)
  double   a;                   // Bounded Pareto alpha value
  double   k;                   // Bounded Pareto k value (bytes)
  double   p;                   // Bounded Pareto p value (bytes)
  int      min_pktlen;          // Minimum packet length (bytes)
  int      max_pktlen;          // Maximum packet length (bytes
  double   data_rate;           // Link data rate (bit/sec)
  int      num_pkt;             // Number of packets to generate
  double   mean_burst_size;     // Computed Pareto mean burst size
  double   var_burst_size;      // Computed Pareto variance burst size
  double   mean_burst_time;     // Computed Pareto mean burst time
  double   intergap_time;       // Computed exponential burst gap time (sec)
  double   pareto_rv;           // Pareto random variable
  double   expon_rv;            // Exponential random variable
  double   time;                // Current time
  double   time_last;           // Last time
  double   min_burst;           // Minimum measured burst size
  double   max_burst;           // Maximum measured burst size
  double   sum_bytes;           // Sum of bytes
  double   sum_bursts;          // Sum of bursts
  int      burst_len;           // Number of packets in a burst
  int      last_pktlen;         // Length of last packet in a burst (bytes)
  int      pktlen;              // Packet length (bytes)
  int      i;                   // Loop counter

  // Output banner
  printf("---------------------------------------------- syntraf1a.c ----- \n");
  printf("-  Program to generate synthetic packet traffic                - \n");
  printf("-   - Bounded Pareto burst size, exponential burst gap time    - \n");
  printf("---------------------------------------------------------------- \n");

  // Prompt for output filename and then create/open the file
  printf("Enter output file name =========================> ");
  scanf("%s", in_string);
  fp = fopen(in_string, "w");
  if (fp == NULL)
  {
    printf("ERROR in creating output file (%s) \n", in_string);
    exit(1);
  }

  // Prompt for random number seed and then use it
  printf("Random number seed (greater than 0) ============> ");
  scanf("%s", in_string);
  rand_val(atoi(in_string));

  // Prompt for mean time between bursts
  printf("Target utilization (%) =========================> ");
  scanf("%s", in_string);
  target_util = atof(in_string);

  // Prompt for Pareto alpha value
  printf("Pareto alpha value =============================> ");
  scanf("%s", in_string);
  a = atof(in_string);

  // Prompt for Pareto k value
  printf("Pareto k value (min burst size in bytes) =======> ");
  scanf("%s", in_string);
  k = atof(in_string);

  // Prompt for Pareto p value
  printf("Pareto p value (max burst size in bytes) =======> ");
  scanf("%s", in_string);
  p = atof(in_string);

  // Prompt for minimum packet length (for last packet in a burst)
  printf("Min packet length (bytes) ======================> ");
  scanf("%s", in_string);
  min_pktlen = atoi(in_string);

  // Prompt for maximum packet length (used for a burst)
  printf("Max packet length (bytes) ======================> ");
  scanf("%s", in_string);
  max_pktlen = atoi(in_string);

  // Prompt for link data rate
  printf("Link data rate (bits per second) ===============> ");
  scanf("%s", in_string);
  data_rate = atof(in_string);

  // Prompt for number of packets to generate
  printf("Number of packets to generate ==================> ");
  scanf("%s", in_string);
  num_pkt = atoi(in_string);

  // Compute the needed mean burst gap time for the given target utilization
  mean_burst_size = (a*(pow(k,a))*((pow(k,(1-a))) - (pow(p,(1-a)))))/((a-1)*(1-pow((k/p),a)));
  var_burst_size = ((a*(pow(k,a))*((pow(k,(2-a))) - (pow(p,(2-a)))))/((a-2)*(1-(pow(k/p,a))))) - pow(mean_burst_size, 2.0);
  mean_burst_time = (8.0 * mean_burst_size) / data_rate;
  intergap_time = (mean_burst_time * (1.0 - (target_util / 100.0))) / (target_util / 100.0);

  //Output message and generate samples
  printf("---------------------------------------------------------------- \n");
  printf("-  Generating %d packets... \n", num_pkt);
  printf("-   - Mean burst size     = %f bytes \n", mean_burst_size);
  printf("-   - Variance burst size = %f bytes \n", var_burst_size);
  printf("-   - Mean intergap time  = %f sec   \n", intergap_time);
  printf("---------------------------------------------------------------- \n");

  // Generate the synthentic traffic
  time = time_last = sum_bytes = sum_bursts = 0.0;
  min_burst = p;
  max_burst = k;
  while(1)
  {
    // Pull random values
    pareto_rv = bpareto(a, k, p);
    expon_rv = expon(intergap_time);
    assert((pareto_rv >= k) && (pareto_rv <= p));
    assert(expon_rv >= 0.0);

    // Update minimum and maximum burst sizes
    if (pareto_rv < min_burst) min_burst = pareto_rv;
    if (pareto_rv > max_burst) max_burst = pareto_rv;
    assert(min_burst >= k);
    assert(max_burst <= p);

    // Determine the interburst time
    time = time + expon_rv;

    // Determine the number of packets in the burst
    burst_len = ceil(pareto_rv / max_pktlen);

    // Determine the size of the last packet in the burst
    if ((pareto_rv / max_pktlen) == burst_len)
      last_pktlen = max_pktlen;
    else
    {
      last_pktlen = ceil((double) max_pktlen *
        ((pareto_rv / max_pktlen) - floor(pareto_rv / max_pktlen)));
      if (last_pktlen < min_pktlen) last_pktlen = min_pktlen;
    }

    // Generate the burst
    for (i=0; i<burst_len; i++)
    {
      if (i< (burst_len - 1))
        pktlen = max_pktlen;
      else
        pktlen = last_pktlen;

      sum_bytes = sum_bytes + pktlen;
      time = time + ((8.0 * pktlen) / data_rate);
      fprintf(fp,"%f   %d \n", (time - time_last), pktlen);
      time_last = time;
      num_pkt--;
      if (num_pkt <= 0) break;
    }
    sum_bursts++;
    if (num_pkt <= 0) break;
  }

  // Output trace summary statistics and close the output file
  printf("---------------------------------------------------------------- \n");
  printf("-  Summary statistics... \n");
  printf("-   - Total time of trace    = %f sec      \n", time);
  printf("-   - Minimum size burst     = %8.0f bytes \n", min_burst);
  printf("-   - Maximum size burst     = %8.0f bytes \n", max_burst);
  printf("-   - Total number of bursts = %8.0f       \n", sum_bursts);
  printf("-   - Total number of bytes  = %8.0f bytes \n", sum_bytes);
  printf("-   - Utilization            = %f %%       \n",
    100.0 * ((8.0 * sum_bytes) / time) / data_rate);
  printf("---------------------------------------------------------------- \n");
  fclose(fp);
}

//===========================================================================
//=  Function to generate bounded Pareto distributed RVs using              =
//=    - Input:  a, k, and p                                                =
//=    - Output: Returns with bounded Pareto RV                             =
//===========================================================================
double bpareto(double a, double k, double p)
{
  double z;          // Uniform random number from 0 to 1
  double pb_rv;      // Computed bounded Pareto value to be returned

  // Pull a uniform random number (0.0 < z < 1.0)
  do
  {
    z = rand_val(0);
  }
  while ((z == 0.0) || (z == 1.0));

  // Generate the bounded Pareto rv using the inversion method
  pb_rv = -(z * pow(p,a) - z * pow(k, a) - pow(p,a)) / (pow(p, a) * pow(k,a));
  pb_rv = pow(pb_rv, (-1.0 / a));

  return(pb_rv);
}

//===========================================================================
//=  Function to generate exponentially distributed random variables        =
//=    - Input:  mean_value                                                 =
//=    - Output: Returns with exponentially distributed random variable     =
//===========================================================================
double expon(double mean_value)
{
  double z;          // Uniform random number (0 < z < 1)
  double exp_rv;     // Computed exponential value to be returned

  // Pull a uniform random number (0.0 < z < 1.0)
  do
  {
    z = rand_val(0);
  }
  while ((z == 0.0) || (z == 1.0));

  // Compute exponential random variable using inversion method
  exp_rv = -mean_value * log(z);

  return(exp_rv);
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - x_n = 7^5*x_(n-1)mod(2^31 - 1)                                    =
//=   - With x seeded to 1 the 10000th x value should be 1043618065       =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=     John Wiley & Sons, 1991. (Page 443, Figure 26.2)                  =
//=========================================================================
double rand_val(int seed)
{
  const long  a =      16807;  // Multiplier
  const long  m = 2147483647;  // Modulus
  const long  q =     127773;  // m div a
  const long  r =       2836;  // m mod a
  static long x;               // Random int value
  long        x_div_q;         // x divided by q
  long        x_mod_q;         // x modulo q
  long        x_new;           // New x value

  // Set the seed if argument is greater than zero and return zero
  if (seed > 0)
  {
    x = seed;
    return(0.0);
  }

  // RNG using integer arithmetic
  x_div_q = x / q;
  x_mod_q = x % q;
  x_new = (a * x_mod_q) - (r * x_div_q);
  if (x_new > 0)
    x = x_new;
  else
    x = x_new + m;

  // Return a random value between 0.0 and 1.0
  return((double) x / m);
}