/* PURPOSE:
 * This program generates a synthetic workload (i.e. a stream of jobs)
 *    that satisfies a Hyper-Erlang distribution(s) of Common Order.
 *
 * The input parameters used in this program were derived from the observed
 *   workload on a 322-node partition of the Cornell Theory Center (CTC) SP2,
 *   over the period from June 25, 1996 to September 12, 1996.
 *
 * Detailed description of the models can be found in the paper titled:
 *            "Modeling of Workload in MPPs"
 * by Jann J. et al, published in the book ISBN 3-540-63574-2 by Springer Verlag,
 * 1997, Lecture Notes in Computer Science vol. 1291, pp. 95-116;
 * Job Scheduling Stategies for Parallel Processing.
 *
 * INPUT:
 * For each model:
 * 1. Number of jobs desired (njobs)
 * 2. Minimum & maximum number of processes per job (pmin & pmax)
 * 3. Inter-arrival time: lambda1, lambda2, order of the HErlang distrib, p1(=prob)
 * 4. Service/CPU time:   lambda1, lambda2, order of the HErlang distrib, p1(=prob)
 * Input to this sample program is hardcoded for simplicity, as this program is
 * meant for illustration purposes only.
 *
 * OUTPUT:
 * For each job in the workload:
 * 1. Number of processors for the job
 * 2. CPU time used by the job (=service time)
 * 3. Arrival time of the job
 *
 * For the whole workload of njobs :
 * 1. 1st, 2nd & 3rd moments of the inter-arrival-time        of the jobs
 * 2. 1st, 2nd & 3rd moments of the service-time (= CPU-time) of the jobs
 *
 * NOTES:
 * To keep this program simple, we have eliminated all checkings of the input.
 * If you use the program often,  we strongly recommend you to add code to
 *     catch error conditions.
 *
 * AUTHOR: Joefon Jann joefon@watson.ibm.com
 *         IBM T. J. Watson Research Center,
 *         P. O. Box 218, Yorktown Heights, NY 10598, USA
 * DATE:   11/1997
 */

#include <stdio.h>
#include <math.h>

#define min(a,b) (((a) < (b)) ? (a) : (b))
#define RANSTREAM 3
#define RANPOOLSIZE 10000

typedef struct
{       double p1,lambda1,lambda2,ur;
        int    order;
}  MODEL_DATA;

/* Global Variables: */
double ran_pool[RANSTREAM][RANPOOLSIZE];
int    ran_remain[RANSTREAM];

/* - - - - - Function prototypes: - - - - - - - - - */
double ur_to_HER_r( double ur,  MODEL_DATA model);

double H_Er_cdf(double x, void *model);

void solve_bisec( double xleft_0,
                  double xright_0,
                  double tol,
                  double *x,
                  double *y,
                  void *model_data,
                  double (*f)(double, void *),
                  int *rc);

double fn_to_solve(double x, void *);

int    fac(int k);

double str_rand(int str);

/* =================================================== */
main(int argc, char **argv)
{  MODEL_DATA inter_arrival_model, cpu_model;
   double inter_arrival_mu1, inter_arrival_mu2,
          inter_arrival_mu3;
   double cpu_mu1, cpu_mu2, cpu_mu3;
   double simtime, inter_arrival_time, cpu_time_used;
   double dtemp;

   int njobs,  pmin,pmax,p_needed;
   int ijob,   nmodels, imodel;
   int i,j;

   /* One needs to change the value of nmodels here to
    * reflect  the number of models (namely number of
    * lines from Tables 1 through 6 in the Results
    * section of this paper) concurrently used for
    * generating the job stream.
    */
   nmodels=1;

   for(imodel=0; imodel<nmodels; imodel++)
   {  simtime=0;

      /* The following INPUT section describes the parameters of model imodel,
       * and should be modified to correspond to your chosen model.
       */
      pmin = 65;
      pmax = 128;

      inter_arrival_model.lambda1 = 3.33e-05;
      inter_arrival_model.lambda2 = 6.07e-04;
      inter_arrival_model.order=2;
      inter_arrival_model.p1=5.03e-01;

      cpu_model.lambda1 = 1.19e-06;
      cpu_model.lambda2 = 3.99e-04;
      cpu_model.order=6;
      cpu_model.p1=2.07e-01;

      njobs=10000;
      /* End of the INPUT section - - - - - - - - - - - - - - - - - - */

      inter_arrival_mu1=0;
      inter_arrival_mu2=0;
      inter_arrival_mu3=0;

      cpu_mu1=0;
      cpu_mu2=0;
      cpu_mu3=0;

      ijob=0;

      /* Here we are using the number of jobs (njobs) as the termination criterion.
       * One may choose to use simtime (i.e. total length of time represented
       *                              in this workload) as a termination criterion.
       * In that case, use simtime in the following WHILE statement.
       */
      while( ijob < njobs)
      {  ijob++;

         dtemp=str_rand(0);
         inter_arrival_time =
                   ur_to_HER_r(dtemp,inter_arrival_model);
         simtime += inter_arrival_time;
         inter_arrival_mu1 +=inter_arrival_time;
         inter_arrival_mu2 +=(inter_arrival_time *
                              inter_arrival_time);
         inter_arrival_mu3 +=(inter_arrival_time *
                              inter_arrival_time *
                              inter_arrival_time);

         dtemp=str_rand(1);
         cpu_time_used = ur_to_HER_r(dtemp,cpu_model);
         cpu_mu1 +=  cpu_time_used;
         cpu_mu2 += (cpu_time_used  *
                     cpu_time_used);
         cpu_mu3 += (cpu_time_used  *
                     cpu_time_used  *
                     cpu_time_used);

         dtemp=str_rand(2);
         dtemp *= (pmax-pmin+1);
         j=dtemp;
         if( ((double) j) < dtemp ) { p_needed = j + pmin; }
         else                       { p_needed = j + pmin - 1; }

         /* One needs to replace the following section (enclosed by #if and #endif)
          * with appropriate statements to output the job stream information.
          * If multiple models are concurrently used, it might be desirable to
          * sort the final output by job arrival time.
          */
         #if 0
             if(ijob < 50 )    /* Just print some jobs to see the values */
             {  printf( "Job used %d CPUs and %e cummulative CPUtime",
                                     p_needed,   cpu_time_used);
                printf( " arrived at time %e \n", simtime);
             }
         #endif
      }  /* end of while loop */

      printf(" For Model %d :\n",imodel);
      inter_arrival_mu1 /= (double)ijob;
      inter_arrival_mu2 /= (double)ijob;
      inter_arrival_mu3 /= (double)ijob;
      printf("Inter-arrival Time Moments of the jobs are \n");
      printf( " %d  mu1 = %e  "
              "     mu2 = %e  "
              "     mu3 = %e \n",
             ijob, inter_arrival_mu1,
                   inter_arrival_mu2,
                   inter_arrival_mu3 );

      cpu_mu1 /= (double)ijob;
      cpu_mu2 /= (double)ijob;
      cpu_mu3 /= (double)ijob;
      printf("CPU Time Moments of the jobs are \n");
      printf(" %d mu1 = %e               "
                 "mu2 = %e "
                 "mu3 = %e\n",
             ijob, cpu_mu1, cpu_mu2, cpu_mu3 );

  } /* end of FOR loop over models */

}   /* end of main program */

/* ================================================== */
double ur_to_HER_r(double ur, MODEL_DATA m)
{  double xleft, xright;
   double x,y;
   int rc;

   m.ur = ur;
   xleft = 1.e-20;        /* a small +ve number */
   xright = -log(1.0-ur) / min(m.lambda1,m.lambda2);
   xright *= 100;

   solve_bisec(xleft, xright,
               (double) 1.e-12,
               &x, &y, &m,
               fn_to_solve, &rc);
   if(rc !=0 )
   {  printf("Unable to find the Random Number for
                              ur = %e \n",m.ur);
   }
   return(x);
}

/* =================================================== */
double fn_to_solve(double x, void *model_data)
{  MODEL_DATA *m;

   m = (MODEL_DATA *) model_data;
   return( H_Er_cdf(x, model_data)  -  m->ur );
}

/* =================================================== */
double H_Er_cdf(double x, void *model_data)
{  MODEL_DATA *p;
   double cdf, t1,t2, dtemp;
   int k;
   t1=0;
   t2=0;
   p = (MODEL_DATA *) model_data;

   for(k=0;  k< p->order;  k++)
   {  dtemp=fac(k);
      t1 += (pow(x*p->lambda1,k)/dtemp);
      t2 += (pow(x*p->lambda2,k)/dtemp);
   }
   cdf = 1.0 -    p->p1 *exp(-(p->lambda1*x))*t1
             - (1-p->p1)*exp(-(p->lambda2*x))*t2;
   return(cdf);
 }

/* =================================================== */
void solve_bisec(double xleft_0,
                 double xright_0,
                 double tol,
                 double *x,
                 double *y,
                 void *model_data,
                 double (*f)(double, void *),
                 int *rc)

{  double xleft, yleft,
          xright, yright,
          xnext,ynext;
   int    ic;

   xleft=xleft_0; xright=xright_0;
   yleft=(*f)(xleft, model_data);
   yright=(*f)(xright,model_data);
   *rc=1; ic=0;

   if( yleft*yright  >  0.0 ) { return; }
   /* No solution in this interval */

   while(ic < 100000 )
   {  ic++;
      xnext=(xleft+xright)/2.0;
      ynext=(*f)(xnext,model_data);
      if( ynext*yright  >  0.0 )
      {   xright=xnext;
          yright=ynext;
      }
      else
      {   xleft=xnext;
          yleft=ynext;
      }
      if( fabs(ynext)  <  tol )
      {   *x=xnext; *y=ynext; *rc=0; return;}
   }
}

/* =================================================== */
int fac(int k)
{  int i,j;
   j=1;

   if(k==0) return (j);
   j=1;
   for(i=1;i<=k;i++) j*=i;
   return(j);
}

/* =================================================== */
double str_rand(int strm)
{  int i;
   if(  ran_remain[strm] <=  0 )
   {  for(i=0;i<RANPOOLSIZE;i++)
          ran_pool[strm][i]=drand48();
          ran_remain[strm]=RANPOOLSIZE;
   }
   ran_remain[strm]--;
   return( ran_pool[strm] [RANPOOLSIZE-ran_remain[strm]]) ;
}