/*
 * Code implementing and using the downey 1997 workload model
 *
 * Written by Allen Downey, SDSC
 *
 * This code is provided as is, with no warranty.
 * If you use it in your work, please include an appropriate
 * acknowledgement.
 */
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <assert.h>

#define NPROCS 400

#define NUM_PES 64          /* 64 processors */
#define MAX_PAR 6           /* maximum parallelism = 2^6 = 64 */
#define TWELVE_HOURS 43200

#define END_OF_TIME FLT_MAX

#define ARRIVAL 0
#define COMPLETION 1
#define DONE 2
#define PRUNED 3

#define UNSCHEDULED 0
#define QUEUED 1
#define RUNNING 2
#define COMPLETED 3

typedef struct event {
  short type, proc;
  float time;
} Event;

typedef struct state {
  short sum;
  short free, next_arrival;
  short stat[NPROCS];
  short pes[NPROCS];
  float start[NPROCS];
  float finish[NPROCS];
  float total_res_time;      /* this is used in find_opt for pruning,
				not for reporting stats */
} State;

/* A Heuristic is a function that takes a State and returns void */
typedef void (*Heuristic) (State *);

/* GLOBAL VARIABLES (and lots of 'em) */

float arrival[NPROCS];
float lifetime[NPROCS];
float A[NPROCS];
float sig[NPROCS];
float max_time[NPROCS];
float min_time[NPROCS];
int mark[NPROCS];
float efficiency[NPROCS];
int nprocs;

float time;

int strat;
int pat;
float param1, param2, param3;
float rho;

double total_time = 0.0;
double total_load = 0.0;
double total_real_load = 0.0;


/* WORKLOAD MODEL */

double factor[7][12] = { 1,1,1,1,1,1,0,0,0,0,0,0,
                         1.5, 1.5, 1.5, .5, .5, .5, 0,0,0,0,0,0,
			 .5, .5, .5, 1.5, 1.5, 1.5, 0,0,0,0,0,0,
			 1, 1, 1.5, 1.5, .5, .5, 0,0,0,0,0,0,
			 1, 1, .5, .5, 1.5, 1.5, 0, 0, 0, 0, 0, 0,
			 1.5, .5, 1.5, .5, 1.5, .5, 0,0,0,0,0,0,
			 .5, 1.5, .5, 1.5, .5, 1.5, 0,0,0,0,0,0 };


/* DRANDOM: choose a random floating-point number between zero and
   one, based on a 31-bit random integer */

double drandom ()
{
  return (double) (random() & 0x7fffffff) / (double) 0x7fffffff;
}

/* CHOOSE FROM EXPONENTIAL : choose a value from an exp distribution
   with the given parameter lambda */

double choose_from_exponential (double lambda)
{
  double x;

  do x = drandom (); while (x == 0.0);
  return -log (drandom ()) / lambda;
}

/* CHOOSE FROM LOG UNIFORM : low and high are the exponents of the
   range; i.e. low = 0 and high = 3 would have a range from 1 second
   to exp(3) seconds */

double choose_from_log_uniform (double low, double high)
{
  double x = drandom () * (high-low) + low;
  return exp(x);
}

double choose_from_multistage ()
{
  float x = drandom ();

  if (x < .9) {
    return choose_from_log_uniform (3.6, 11.1);
  } else {
    return choose_from_log_uniform (11.1, 14.5);
  }
}

double choose_lifetime ()
{
  if (param1 == 0.0 && param2 == 0.0) {
    return choose_from_multistage ();
  } else {
    return choose_from_log_uniform (param1, param2);
  }
}

double avg_lifetime ()
{
  if (param1 == 0.0 && param2 == 0.0) {
    return 64306.385766;
  } else {
    return (exp (param2) - exp(param1)) / (param2 - param1);
  }
}

/* CHOOSE PARALLELISM : log uniform distribution between 1 and 2^max */

double choose_parallelism (int max)
{
  double x = drandom() * max;
  return pow (2.0, x);
}

/* CHOOSE SIGMA : uniform distribution between 0 and 2 */

double choose_sigma ()
{
  return (drandom() * 2.0);
}

float calc_real_load (State *state)
{
  int i;
  float load = 0.0;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] != RUNNING) continue;
    load += state->pes[i] * efficiency[i];
  }
  return load / NUM_PES;
}


/* GENERATE STATISTICS */


void check_state (State *state, float t)
{
  static FILE *fp1 = NULL;
  static FILE *fp2 = NULL;
  static float last_time = 0.0;
  static float load = 0.0;
  static float real_load = 0.0;
  float dt;

  if (fp1 == NULL) {
    fp1 = fopen ("load_avg", "w");
    fp2 = fopen ("real_load_avg", "w");
  }

  /* stop monitoring the load after the arrival of the last job;
     that way, a strategy that does well will not be punished during
     the quiescent period if it runs out of work early */
  if (t > arrival[nprocs-1]) return;

  dt = t - last_time;
  total_time += dt;
  total_load += load * dt;
  total_real_load += real_load * dt;

  last_time = t;
  fprintf (fp1, "%lf\t%lf\n", t, load);
  fprintf (fp2, "%lf\t%lf\n", t, real_load);
  load = (float)(NUM_PES - state->free) / (float)NUM_PES;
  real_load = calc_real_load (state);
  fprintf (fp1, "%lf\t%lf\n", t, load);
  fprintf (fp2, "%lf\t%lf\n", t, real_load);
}

void print_sched (State *state)
{
  int i;
  float run_time, queue_time, slowdown, eff;
  FILE *fp = fopen ("sched_info", "w");

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == UNSCHEDULED || state->stat[i] == QUEUED) {
      state->pes[i] = 0;
      run_time = queue_time = -1.0;
      efficiency[i] = -1.0;
    }
    run_time = state->finish[i] - state->start[i];
    queue_time = state->start[i] - arrival[i];
    slowdown = (run_time + queue_time) / min_time[i];
    fprintf (fp, "%4d   %12.1lf   %12.1lf   %8.1lf   %5.2lf   ",
	     i, arrival[i], lifetime[i], A[i], sig[i]);
    fprintf (fp, "%4d   %12.1f   %12.1f   %12.1f   %12.3f   %8.4f\n",
	     state->pes[i], queue_time, run_time, min_time[i], slowdown,
	     efficiency[i]);
  }
}

void print_summary (State *state)
{
  int i;
  double hours, min;
  int num_procs = 0;
  double total_queue_time = 0.0;
  double total_res_time = 0.0;
  double total_slowdown = 0.0;
  int total_cluster_size = 0;

  /* summary statistics include all jobs */
  for (i=0; i<nprocs; i++) {
    assert (state->stat[i] == RUNNING || state->stat[i] == COMPLETED);
    num_procs++;
    total_cluster_size += state->pes[i];
    total_queue_time += state->start[i] - arrival[i];
    total_res_time += state->finish[i] - arrival[i];
    total_slowdown += (state->finish[i] - arrival[i]) / min_time[i];
  }
  hours = total_time/3600.0;
/*  printf ("%s\t%2d   %8.4lf   ", filename, strat, param1);  */
  printf (
"%8.2lf   %5d   %8.4lf   %8.4lf   %10.1lf   %10.1lf   %10.1lf   %6.2lf\n",
	  hours, num_procs, total_load/total_time, total_real_load/total_time,
	  total_queue_time/num_procs, total_res_time/num_procs,
	  total_slowdown/num_procs, (double)total_cluster_size/num_procs);
}



/* Sa: calculate the slowdown on n processors given average parallelism A
   and coefficient of variation (squared) cv2 */
 
float Sa (int n, float A, float cv2)
{
  if (cv2 <= 1.0) {
 
    /* low variance model */
    if (n <= A) {
      return A*n / (A + cv2/2 * (n-1));
    } else if (n < 2*A - 1) {
      return A*n / (cv2 * (A - 0.5) + n * (1 - cv2/2));
    } else {
      return A;
    }
  } else {
    /* high variance model */
    if (n < A*cv2 + A - cv2) {
      return n*A * (cv2+1) / (cv2 * (n+A-1) + A);
    } else {
      return A;
    }
  }
}
 
inline float run_time (int proc, int pes)
{
  float speedup = Sa (pes, A[proc], sig[proc]);
  return lifetime[proc] / speedup;
}

inline float min_total_res_time (State *state)
{
  int i;
  float total = state->total_res_time;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      total += time + min_time[i] - arrival[i];
    }
  }
  return total;
}


/* EVENTS */

void print_event (Event *event)
{
  switch (event->type) {
  case ARRIVAL:
    printf ("%d ARRIVES t=%g\n", event->proc, event->time);
    break;
  case COMPLETION:
    printf ("%d COMPLETES t=%g\n", event->proc, event->time);
    break;
  case DONE:
    printf ("DONE\n");
    break;
  }
}

/* GET_NEXT_EVENT: sets the fields of the event struct and returns
   the time of the next event.  If the event is an arrival, then
   the next_arrival field of state gets incremented. */

float get_next_event (State *state, Event *event)
{
  int i, proc;
  float next_arrival_time, min;

  if (state->next_arrival == nprocs) next_arrival_time = END_OF_TIME;
  else next_arrival_time = arrival[state->next_arrival];

  /* find next decedent */
  min = END_OF_TIME;
  for (i=0; i<nprocs; i++) {
    if (state->stat[i] != RUNNING) continue;
    if (state->finish[i] < min) {
      min = state->finish[i];
      proc = i;
    }
  }

  /* if the next death precedes the next arrival ... */
  if (min < next_arrival_time) {

    /* construct completion event */
    event->type = COMPLETION;
    event->time = min;
    event->proc = proc;
  } else {

    /* construct arrival event */
    if (next_arrival_time == END_OF_TIME) {
      event->type = DONE;
      event->time = END_OF_TIME;
    } else {
      event->type = ARRIVAL;
      event->time = next_arrival_time;
      event->proc = state->next_arrival;
      state->next_arrival++;
    }
  }
  return event->time;
}

void unget_event (State *state, Event *event)
{
  if (event->type == ARRIVAL) {
    state->next_arrival--;
  }
}

float process_events (State *state)
{
  float t;
  Event event[1];

  t = get_next_event (state, event);
  do {
#ifdef DEBUG
    print_event (event);
#endif

    switch (event->type) {
    case ARRIVAL:
      state->stat[event->proc] = QUEUED;
      break;
    case COMPLETION:
      state->stat[event->proc] = COMPLETED;
      state->free += state->pes[event->proc];
      check_state (state, t);
      break;
    case DONE:
      return event->time;
    default:
      assert (0);
    }
  } while (get_next_event (state, event) == t);
  unget_event (state, event);

  return t;
}

void start_job (State *state, int proc, int pes)
{
  float rt = run_time (proc, pes);

  state->free -= pes;
  state->stat[proc] = RUNNING;
  state->pes[proc] = pes;
  state->start[proc] = time;
  state->finish[proc] = time + rt;
  state->total_res_time += state->finish[proc] - arrival[proc];
  efficiency[proc] = lifetime[proc] / rt / pes;
  check_state (state, time);

#ifdef DEBUG
    printf ("%d STARTS t=%g, pes = %d, dur = %.1lf\n",
	    proc, time, pes, rt);
#endif
  /* if we have just scheduled the last job, we can stop */
  if (proc == nprocs-1) {
    print_summary (state);
    print_sched (state);
    exit (0);
  }
}


/* STRATEGIES: SMART FIFO */

int find_next_decedent (State *state)
{
  int i, proc;
  float min = END_OF_TIME;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] != RUNNING || mark[i] == 1) continue;
    if (state->finish[i] < min) {
      min = state->finish[i];
      proc = i;
    }
  }
  if (min == END_OF_TIME) return -1;
  mark[proc] = 1;
  return proc;
}

/* QUEUE TIME: calculates the actual queue time a job will experience */

float queue_time (State *state, int pes)
{
  int i, total;

  /* unmark all processes */
  for (i=0; i<nprocs; i++) {
    mark[i] = 0;
  }
  if (pes <= state->free) {
    return 0;
  }
  total = state->free;
  do {
    i = find_next_decedent (state);
    assert (i != -1);              
    total += state->pes[i];
  } while (total < pes);
  return state->finish[i] - time;
}

int best_decision (State *state, int proc)
{
  int i, best;
  float res_time, min = END_OF_TIME;

#ifdef DEBUG
  print_witnesses (state);
  printf ("%d\tfree\n\n", state->free);
#endif
  for (i=1; i<=NUM_PES; i++) {
#ifdef DEBUG
    printf ("%d\t%lf\t%lf\n", i, queue_time (state, i), run_time (proc, i));
#endif
    res_time = queue_time (state, i) + run_time (proc, i);
    if (res_time < min) {
      min = res_time;
      best = i;
    }
  }
#ifdef DEBUG
  printf ("\n");
#endif
  return best;
}

void smart_fifo (State *state)
{
  int i, pes;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = best_decision (state, i);
      if (pes <= state->free) {
	start_job (state, i, pes);
      } else {
	break;
      }
    }
  }
}

/* STRATEGIES: MAXIMUM POWER, GREEDY FIFO, SEVCIK */

float max_power_pes (int i)
{
  float ans;

  if (sig[i] <= 1.0) {
    ans = sig[i] * (A[i] - 0.5) / (1.0 - sig[i]/2.0);
    if (ans < A[i]) ans = A[i];
    return ans;
  } else {
    return (A[i] * (sig[i] + 1.0) - sig[i]) / sig[i];
  }
}

float max_useful_pes (int i)
{
  if (sig[i] == 0.0) {
    return A[i];
  } else if (sig[i] <= 1.0) {
    return 2*A[i] - 1.0;
  } else {
    return sig[i] * A[i] + A[i] - sig[i];
  }
}

void max_useful (State *state)
{
  int i, pes, max;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = state->free;
      max = (int) ceil (max_useful_pes (i));
      if (pes > max) pes = max;

      assert (pes >= 1 && pes <= NUM_PES);
      start_job (state, i, pes);
    }
    if (state->free == 0) break;
  }
}  

void sevcik (State *state)
{
  int i, pes, max;
  float h, s;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = state->free;

      h = rho * (A[i]-1);
      if (time > 12*3600) h = 0.0;
      s = sig[i]/2.0;
      if (s > 1) s = 1;
      max = (int) ceil (A[i] - s * h);
      if (pes > max) pes = max;

      assert (pes >= 1 && pes <= NUM_PES);
      start_job (state, i, pes);
    }
    if (state->free == 0) break;
  }
}  

void simp_sevcik (State *state)
{
  int i, pes, max;
  float h, s;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = state->free;

      h = rho * (A[i]-1);
      if (time > 12*3600) h = 0.0;

      /* all jobs as if sigma = 1.0 */
      s = 0.5;
      max = (int) ceil (A[i] - s * h);
      if (pes > max) pes = max;

      assert (pes >= 1 && pes <= NUM_PES);
      start_job (state, i, pes);
    }
    if (state->free == 0) break;
  }
}  

void time_sevcik (State *state)
{
  int i, pes, max, period;
  float h, s, k;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = state->free;

      period = time/7200;
      k = factor[pat][period];
      h = k * rho * (A[i]-1);

      /* all jobs as if sigma = 1.0 */
      s = 0.5;
      max = (int) ceil (A[i] - s * h);
      if (pes > max) pes = max;

      assert (pes >= 1 && pes <= NUM_PES);
      start_job (state, i, pes);
    }
    if (state->free == 0) break;
  }
}  

void mult_average (State *state)
{
  int i, pes, max;
  float h, s;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = state->free;

      max = (int) ceil (1.5*A[i]);
      if (pes > max) pes = max;

      assert (pes >= 1 && pes <= NUM_PES);
      start_job (state, i, pes);
    }
    if (state->free == 0) break;
  }
}

void average (State *state)
{
  int i, pes, max;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = state->free;

      max = (int) ceil (A[i]);
      if (pes > max) pes = max;

      assert (pes >= 1 && pes <= NUM_PES);
      start_job (state, i, pes);
    }
    if (state->free == 0) break;
  }
}  

void max_power (State *state)
{
  int i, pes, max;

  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      pes = state->free;

      max = (int) ceil (max_power_pes (i));
      if (pes > max) pes = max;

      assert (pes >= 1 && pes <= NUM_PES);
      start_job (state, i, pes);
    }
    if (state->free == 0) break;
  }
}  

/* ALLOCATE_EQUIPART : returns the number of jobs that were started */

int allocate_equipart (State *state)
{
  int i;
  int pes, max;
  int count = 0;

  if (state->free == 0) return 0;

  /* how many jobs are in queue ? */
  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {
      count++;
    }
  }
  if (count == 0) return 0;

  /* partition = min (equipartition, max_useful_pes) */

  pes = (int) ceil ((float) state->free / (float) count);
  for (i=0; i<nprocs; i++) {
    if (state->stat[i] == QUEUED) {

      max = (int) ceil (max_useful_pes (i));
      if (pes > max) pes = max;
      start_job (state, i, pes);
      break;
    }
  }
  assert (state->free >=0 && state->free <= NUM_PES);
  return 1;
}

void equipartition (State *state)
{
  /* keep calling allocate_equipart until there are no jobs in queue */
  while (allocate_equipart (state));
}  

/* FIND_HEURISTIC: from the given initial state, finds a schedule for
   the set of jobs, using the provided heuristic */

void find_heuristic (State *state, Heuristic heur)
{
  int i;
  Event event[1];

  while (1) {
    time = process_events (state);
    if (time == END_OF_TIME) {
      return;
    }

    if (state->free) {
      heur (state);
    }
  }
}

void generate_arrivals (double lambda)
{
  int i, period;
  double k;

  time = 0.0;
  for (i=0; i<NPROCS; i++) {
    arrival[i] = time;
    lifetime[i] = choose_lifetime ();
    A[i] = choose_parallelism (MAX_PAR);
    sig[i] = choose_sigma ();
    max_time[i] = run_time (i, 1);
    min_time[i] = run_time (i, NUM_PES);

    period = time/7200;
    k = factor[pat][period];
    time += choose_from_exponential (k * lambda);

    /* no arrivals after hour 18 */
    if (time > 18*3600) break;
  }
  if (i == NPROCS) {
    printf ("Not enough processes.\n");
    exit (1);
  }
  nprocs = i;
}

void read_arrivals (FILE *fp)
{
  int i;

  for (i=0; i<nprocs; i++) {
    if (fscanf (fp, "%f %f %f %f",
		&arrival[i], &lifetime[i], &A[i], &sig[i]) != 4) {
      printf ("Error: cannot read input file.\n");
      exit(1);
    }
/*    if (param1 == 0.0) sig[i] = 1.0;  */
    max_time[i] = run_time (i, 1);
    min_time[i] = run_time (i, NUM_PES);
  }
}

void print_arrivals ()
{
  int i;
  FILE *fp = fopen ("proc_info", "w");

  for (i=0; i<nprocs; i++) {
    fprintf (fp, "%10.3lf\t%10.2lf\t%10.2lf\t%10.2lf\n",
	     arrival[i], lifetime[i], A[i], sig[i]);
  }
}

#define NUM_HEURS 7
Heuristic heur[NUM_HEURS] = 
{ max_useful, equipartition, sevcik, simp_sevcik, average, max_power,
time_sevcik };

main (int argc, char *argv[])
{
  int i;
  State state[1];
  int seed;
  double lambda, mu;

  if (argc == 7) {
    rho = atof (argv[1]);      /* load */
    seed = atoi (argv[2]);     /* random # seed */
    strat = atoi (argv[3]);    /* which strategy */
    pat = atoi (argv[4]);      /* which arrival pattern */
    param1 = atof (argv[5]);   /* these parameters mean different things */
    param2 = atof (argv[6]);   /* depending on the other arguments */
  } else {
    printf ("Usage: sim load seed strategy tmax tmin.\n");
    exit (1);
  }

  assert (pat >=0 && pat <= 6);
  srandom (seed);

  /* mu is average lifetime */
  mu = avg_lifetime ();
  lambda = rho * NUM_PES / mu;
  
  generate_arrivals (lambda);
  print_arrivals ();

  /* initialize state */
  state->free = NUM_PES;
  state->next_arrival = 0;
  state->total_res_time = 0.0;
  for (i=0; i<nprocs; i++) {
    state->stat[i] = UNSCHEDULED;
  }
  time = 0.0;

  if (strat >= 0 && strat < NUM_HEURS) {
    find_heuristic (state, heur[strat]);
  } else {
    printf ("There is no strategy %d.\n", strat);
    exit(1);
  }

  print_summary (state);
  print_sched (state);

  return 0;
}