/*
  Copyright 2005 Allen B. Downey
 
    This file contains a simulator I used to study TCP performance.
    It is research-grade software, which is to say that is was
    designed primarily for my use, and is provided primarily as
    a reference implementation of the algorithms described in my
    papers.  The code should be readable, and demonstrably correct,
    but is not particularly user-friendly.  I may be able to answer
    questions about the algorithms implemented here, but in general
    I cannot provide support for anyone using this program.
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program; if not, see http://www.gnu.org/licenses/gpl.html
    or write to the Free Software Foundation, Inc., 51 Franklin St,
    Fifth Floor, Boston, MA  02110-1301  USA
*/

#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <assert.h>
#include "../../cdf/cdflib.h"

extern double round(double x);

#define RECV   1           // packet arrives at receiver
#define SENDER 2           // ACK arrives at sender
#define BOTTLE 3           // packet arrives at bottleneck
#define LAST   4           // packet arrives at the home stretch
#define OBS    5

typedef struct event {
  int type, seqno, tripdup;
  double time, size;
  struct conn *conn;
} Event;

typedef struct heap {
  Event **v;
  int next, size;
} Heap;

typedef struct dist {
  double *pdf;
  int n;
} Dist;

typedef struct conn {
  int id;
  FILE *tfp, *sfp, *cwfp, *distcwfp, *ffp, *queuefp, *seenfp;
  int done, last;
  int seqno, last_ack, slow_start;
  int start_seqno;
  double start_time, start_size;
  Heap *tripdups;
  double last_arr1, last_arr2;
  double total_sent, total_recv;
  double cw, ssthresh;
  double total_size;
  double total_obs;
  double pq_obs;
  double total_cw;
  double total_queue, total_pq;
  double queue;
  double mu, sigma;
  int num_drops;
  double t1, t2, t3, t4, rtt, b, bw;
  Vector *timing;
  struct conn *other;
} Conn;


Heap *h;
Conn *conn, *conn2, *conn3;
int started = 0;
int finished = 0;

double last_arr3 = 0.0;
double last_dep3 = 0.0;

double gtime = 0.0;

int no_delay = 1;
int done_obs = 0;
int first_delay_seqno = 0;

double mss = 1500.0;

int rseed = 17;
int drop1 = 300000000;
int kmin = 62;
int num_tripdup = 0;
int num_drop = 0;

int randdrop = 0;
int periodic = 1;
int verbose;

int ack_delay_seqno = -2;
int data1_delay_seqno = -2;
int data2_delay_seqno = -2;
double delay = 300.0;

Vector *starts, *cw1, *cw2, *recv1, *recv2;

void *mymalloc (int size) {
  void *result = malloc (size);
  assert (result != NULL);
  return result;
}

// DIST

Dist *make_dist (int n)
{
  Dist *dist = (Dist *) mymalloc (sizeof (Dist));
  dist->pdf = (double *) mymalloc (n * sizeof (double));
  dist->n = n;
  return dist;
}

void zero_dist (Dist *dist)
{
  int i;

  for (i=0; i<dist->n; i++) {
    dist->pdf[i] = 0.0;
  }
}

/* pdf_mean: use the pdf to calculate the mean */

double pdf_mean (Dist *dist)
{
  int i;
  double count = 0.0;
  double total = 0.0;

  for (i=0; i<dist->n; i++) {
    count += dist->pdf[i];
    total += dist->pdf[i] * i;
  }
  return total / count;
}


/* normalize_dist: using the pdf, normalize
   the distribution such that the pdf integrates to 1.0 */

void normalize_dist (Dist *dist)
{
  int i;
  double total = 0.0;

  for (i=0; i<dist->n; i++) {
    total += dist->pdf[i];
  }
  if (total == 0.0) return;
  for (i=0; i<dist->n; i++) {
    dist->pdf[i] /= total;
  }
}

/* normal_dist: generate a normal distribution with the given
   parameters and put the values into the distribution.  Note
   that the resulting dist is not normalized!!! */

void normal_dist (Dist *dist, double mean, double sd)
{
  int i;
  double xs;

  for (i=0; i<dist->n; i++) {
    xs = ((double)i - mean) / sd;
    dist->pdf[i] = exp (-0.5 * xs*xs);
  }
}

/* print_pdf: print the pdf array */

void print_pdf (FILE *fp, Dist *dist, char *name)
{
  int i;

  if (name != NULL) {
    fprintf (fp, "\n\"%s\n", name);
  }

  for (i=0; i<dist->n; i++) {
    fprintf (fp, "%d\t%5g\n", i, dist->pdf[i]);
  }
}

/* print_cdf2: print the pdf array as a cdf 
   (assumes that the pdf is normalized) */

void print_cdf2 (FILE *fp, Dist *dist, char *name)
{
  int i;
  double total = 0.0;

  if (name != NULL) {
    fprintf (fp, "\n\"%s\n", name);
  }

  for (i=0; i<dist->n; i++) {
    total += dist->pdf[i];
    fprintf (fp, "%d\t%5g\n", i, total);
  }
}


// EVENT

Event *make_event (int type, double time, Conn *conn) {
  Event *e = (Event *) mymalloc (sizeof (Event));
  e->type = type;
  e->seqno = 0;
  e->tripdup = 0;
  e->time = time;
  e->size = mss;
  e->conn = conn;
  return e;
}

Event *remake_event (Event *e, int type, double time) {
  e->type = type;
  e->time = time;
  return e;
}

void event_print (Event *e) {
  printf ("%.2lf\t%d\n", e->time, e->type);
}

Heap *make_heap (int size) {  
  Heap *h = (Heap *) mymalloc (sizeof (Heap));
  h->v = (Event **) mymalloc (size * sizeof (Event *));
  h->next = 1;
  h->size = size;
  return h;
}

int heap_is_empty (Heap *h) {
  return h->next == 1;
}

void heap_resize (Heap *h) {
  h->size *= 2;
  h->v = (Event **) realloc (h->v, h->size * sizeof (Event *));
}

Event *heap_cargo (Heap *h, int i) {
  if (i>0 && i < h->next) {
    return h->v[i];
  } else {
    return NULL;
  }
}

void heap_print (Heap *h) {
  int i;
  printf ("Heap contains:\n");
  for (i=1; i<h->next; i++) {
    event_print (heap_cargo (h, i));
  }
}

int heap_compare (Heap *h, int i, int j) {
  double t1, t2;
  Event *e1 = heap_cargo (h, i);
  Event *e2 = heap_cargo (h, j);
  if (e1 == NULL) {
    if (e2 == NULL) {
      return 0;
    } else {
      return -1;
    }
  }
  if (e2 == NULL) return 1;
  t1 = e1->time;
  t2 = e2->time;
  if (t1 < t2) return 1;
  if (t2 < t1) return -1;
  return 0;
}

void heap_swap (Heap *h, int i, int j) {
  Event *temp = h->v[i];
  h->v[i] = h->v[j];
  h->v[j] = temp;
}

void heap_add (Heap *h, Event *e) {
  int parent;
  int i = h->next;

  if (i == h->size) {
    heap_resize (h);
  }

  h->v[i] = e;
  h->next++;

  while (i>1) {
    parent = i/2;
    if (heap_compare (h, parent, i) > 0) break;
    heap_swap (h, i, parent);
    i = parent;
  }
}

void heap_reheapify (Heap *h, int i) {
  int left = i*2;
  int right = i*2+1;
  int winleft, winright;

  winleft = heap_compare (h, i, left);
  winright = heap_compare (h, i, right);
  if (winleft>=0 && winright>=0) return;

  if (heap_compare (h, left, right) > 0) {
    heap_swap (h, i, left);
    heap_reheapify (h, left);
  } else {
    heap_swap (h, i, right);
    heap_reheapify (h, right);
  }
}

Event *heap_peek (Heap *h) {
  if (heap_is_empty (h)) return NULL;
  return h->v[1];
}

Event *heap_remove (Heap *h) {
  Event *result;
  if (heap_is_empty (h)) return NULL;
  result = h->v[1];

  h->next--;
  h->v[1] = h->v[h->next];
  heap_reheapify (h, 1);
  return result;
}

void heap_mark_tripdup (Heap *h, int seqno, Conn *conn) {
  int i;
  Event *e, *t;

  // find the packet that will be the third dup ACK
  // and mark it as a tripdup
  for (i=1; i<h->next; i++) {
    e = heap_cargo (h, i);
    if (e->seqno == seqno) {
      e->tripdup = 1;
      return;
    }
  }

  // if we get here, then the packet that will be the
  // third dup doesn't exist yet.  We have to make a
  // pseudoevent and put it into the heap of tripdups,
  // which keeps track of future seqnos that will be
  // tripdups when they are created

  t = make_event (0, (double) seqno, NULL);
  heap_add (conn->tripdups, t);
  // printf ("adding seqno %d to the heap\n", seqno);
}

//-----------------------------------------------------------

FILE *open_file (char *name, int n)
{
  char filename[128];
  if (n == 1) {
    return fopen (name, "w");
  } else {
    sprintf (filename, "%s%d", name, n);
    return fopen (filename, "w");
  }
}

void calc_rtt (Conn *conn)
{
  conn->rtt = conn->t1 + conn->t2 + conn->t3 + conn->t4;    // rtt in ms
  conn->b = conn->rtt / conn->t3;                // bdp in packets
  conn->bw = 1 / conn->t3 * 1000.0;               // in packets / second
  //  printf ("rtt = %.1lf, b = %.1lf, bw = %.1lf\n",
  //  conn->rtt, conn->b, conn->bw);
}

void init_conn (Conn *conn) {
  static int count = 0;

  count++;
  conn->id = count;
  conn->tfp = open_file ("timing", count);
  conn->sfp = open_file ("send", count);
  conn->cwfp = open_file ("cw", count);
  conn->distcwfp = open_file ("cw.sim", count);
  conn->ffp = open_file ("inflight", count);
  conn->queuefp = open_file ("queue", count);
  conn->seenfp = open_file ("seen", count);
  conn->timing = make_vector (128);
}

Conn *make_conn ()
{
  Conn *conn = (Conn *) mymalloc (sizeof (Conn));
  double mult;

  conn->done = 0;
  conn->last = 0;
  conn->seqno = 1;
  conn->last_ack = 0;
  conn->last_arr1 = 0.0;
  conn->last_arr2 = 0.0;
  conn->start_seqno = 0;
  conn->start_time = 0.0;
  conn->start_size = 0.0;

  conn->slow_start = 1;
  conn->num_drops = 0;
  conn->tripdups = make_heap (128);
  conn->total_sent = 0.0;
  conn->total_recv = -1.0;
  conn->cw = 0.0;
  conn->queue = 0.0;
  conn->ssthresh = 150000.0;
  conn->total_size = 150000.0;
  conn->total_cw = 0.0;
  conn->total_queue = 0.0;
  conn->total_pq = 0.0;
  conn->pq_obs = 0.0;

  conn->mu = 0.0;
  conn->sigma = 0.0;

  conn->t1 = 100.0;
  conn->t2 = 30.0;
  conn->t3 = 20.0;
  conn->t4 = 50.0;

  mult = 1.0;
  conn->t1 *= mult;
  conn->t2 *= mult;
  conn->t3 *= mult;
  conn->t4 *= mult;


  calc_rtt (conn);
  init_conn (conn);
  return conn;
}

Conn *copy_conn (Conn *src)
{
  Conn *conn = (Conn *) mymalloc (sizeof (Conn));
  bcopy (src, conn, sizeof(Conn));
  init_conn (conn);
  return conn;
}
//---------------------------------------------------------

double rand_uniform ()
{
  return (double) random() / ((double) RAND_MAX + 1);
}

/* rand_exp: random number from an exponential distribution
   with mean mu */

double rand_exp (double mu)
{
  return -log(rand_uniform()) * mu;
}

/* rand_normal: this is gasdev from Numerical Recipes in C */

double rand_normal ()
{
  static int iset = 0;
  static double gset;
  double fac, r, v1, v2;

  // we produce two value each time, so if one is cached,
  // just return it; otherwise, produce two more

  if (iset == 0) {
    do {
      v1 = 2.0 * rand_uniform() - 1.0;
      v2 = 2.0 * rand_uniform() - 1.0;
      r = v1*v1 + v2*v2;
    } while (r >= 1.0);
    fac = sqrt (-2.0 * log (r)/r);
    gset = v1 * fac;
    iset = 1;
    return v2 * fac;
  } else {
    iset = 0;
    return gset;
  }
}

double rand_lognormal (double mu, double sigma, Conn *conn)
{
  double x, y;
  if (no_delay || conn->seqno < first_delay_seqno) {
    return 0.0;
  }

  x = rand_normal() * sigma + mu;
  y = exp (x);
  return y;
}

int random_drop (int k) {
  return (random() % k == 0);
}

//-------------------------------------------------------------

// AVG_THROUGH: find the average throughput during an interval
// with no drops, given initial congestion window (after halving)
// of c0 packets and k packets until the next drop.  The return 
// value is average throughput; the location referenced by t gets
// set to the elapsed time of the interval, in packet tts.

double avg_through (double c0, double k, double *t, Conn *conn) {
  double kca, cavg, tca, tsc;
  double b = conn->b;

  if (k == 0) {
    *t = 0.0;
    return 0.0;
  }

  if (c0 >= b) {
    *t = k;
    return conn->bw * (k-1) / k;
  }

  // kca is the number of packet received while in CA
  kca = (b*b - c0*c0) / 2.0;

  /* this special case is nice, but unnecessary
  if (kca >= k) {
    c1 = sqrt (2*k + c0*c0);
    cavg = 2.0/3.0 * (c1*c1*c1 - c0*c0*c0) / (c1*c1 - c0*c0);
    tca = b * kca / (cavg - 0.5);
    *t = tca;
    return conn->bw * (kca-1) / tca;
  }
  */

  // cavg is the average congestion window during CA
  cavg = 2.0/3.0 * (b*b*b - c0*c0*c0) / (b*b - c0*c0);

  // tca is the real time spent in CA (in packet tts)
  tca = b * kca / (cavg - 0.5);

  // tsc is the real time spent in SC
  tsc = k - kca;

  *t = tca + tsc;
  return conn->bw * (k-1) / (tca + tsc);  
}

void apply_weight_dist (Dist *dist, double w, double c0, double p) {
  int i, j;
  double d;

  // printf ("\nupdate %lf\n", c0);
  i = 0;
  while (1) {
    j = (int) ceil (sqrt (2*i + c0*c0));
    if (j >= dist->n) break;

    d = p * w;
    dist->pdf[j] += d;
    w -= d;
    // printf ("%d\t%lf\n", i, d);
    i++;
  }
}

void update_dist (Dist *d1, Dist *d2, double p) {
  int i;
  double w, c0;

  zero_dist (d2);
  for (i=0; i<d1->n; i++) {
    w = d1->pdf[i];
    c0 = i / 2.0;
    apply_weight_dist (d2, w, c0, p);
  }
}

void apply_weight_dist2 (Dist *dist, double w, double c0, double p, Conn *conn) {
  int i, j;
  double t, th, d;

  //printf ("\nupdate %lf\n", c0);
  i = 1;
  while (1) {
    th = avg_through (c0, i, &t, conn);
    j = (int) round (th);
    assert (j>=0 && j<dist->n);

    d = p * w;
    dist->pdf[j] += d * t;
    w -= d;
    if (d < 1e-5) break;

    // printf ("%d\t%d\t%lf\t%lf\n", i, j, t, d);
    i++;
  }
}

void update_dist2 (Dist *d1, Dist *d2, double p, Conn *conn) {
  int i;
  double w, c0;

  zero_dist (d2);
  for (i=0; i<d1->n; i++) {
    w = d1->pdf[i];
    c0 = i / 2.0;
    apply_weight_dist2 (d2, w, c0, p, conn);
  }
}

// CALC_CW_DIST: compute the distribution of congestion window
// given kmin average packets between drops

void calc_cw_dist (Conn *conn) {
  Dist *d1, *d2, *d3;
  FILE *fp;
  double p = 1.0 / kmin;
  double a, cavg, mean;
  int n;
  int i;

  // cavg is the mean cw given periodic drops
  // n = 4cavg = the upper bound of our approximate distribution
  a = sqrt (2.0/3.0 * kmin);
  cavg = 14.0/9.0 * a;
  n = (int) round (cavg * 4.0);

  //printf ("p = %lf\n", p);
  //printf ("a = %lf\n", a);
  //printf ("cavg = %lf\n", cavg);

  d1 = make_dist (n);
  d2 = make_dist (n);

  // d1 contains the initial estimate of the distribution of cw
  normal_dist (d1, cavg, cavg/3.0);
  normalize_dist (d1);

  fp = fopen ("dist.cw.pred", "w");

  // each iteration of the loop does two improvements:
  // the improvement of d1 goes into d2 and then
  // the improvement of d2 goes into d1
  // four iterations is _plenty_

  for (i=0; i<4; i++) {
    // print_pdf (fp, d2, "dist1");
    update_dist (d1, d2, p);
    normalize_dist (d2);

    // print_pdf (fp, d2, "dist1");
    update_dist (d2, d1, p);
    normalize_dist (d1);
  }
  print_cdf2 (fp, d1, NULL);
  fclose (fp);


  // d3 contains the distribution of throughputs, using
  // d1 as the distribution of cw at the beginning of each
  // interval

  d3 = make_dist (100);
  update_dist2 (d1, d3, p, conn);
  normalize_dist (d3);

  fp = fopen ("dist.through.pred", "w");
  print_cdf2 (fp, d3, NULL);
  fclose (fp);

  mean = pdf_mean (d3);
  printf ("through.analyis.random = %.3lf\n", mean);
}


// MAKE_DENOM: builds a vector that is the number of packets
// in queue as seen by the ith packet.  It is the denominator
// of the fraction when we divide observed inter-packet excess
// delay by the number of packet tts we think the delay reflects.

Vector *make_denom (int len, int icw, int b)
{
  int j, k, i = 0;
  int n = 1, m = 1;
  int cw;
  Vector *denom = make_vector (len);

  cw = icw;
  while (cw < b) {
    cw *= 2;
    m++;
  }
  cw /= 2;

  for (j=0; j<m; j++) {
    for (k=0; k<n; k++) {
      add_to_vector (denom, (double) k);
      i++;
      if (i == len) return denom;
    }
    n = n*2;
  }

  k = 2*cw-b;

  while (1) {
    add_to_vector (denom, (double) k);
    k++;
    i++;
    if (i == len) return denom;
  }
  return denom;
}

// PROCESS_TIMING: during slow start, we know which data packet
// is triggered by which ACK.  Therefore there are certain pairs
// of packets whose spacing is expected to be rtt + q tt
// where q is the number of packets we expect to be in queue
// when the packet arrives, and tt is transmission time on the
// bottleneck

/*
void process_timing (Conn *conn) {
  Vector *vector = conn->timing;
  int i;
  double diff;
  FILE *fp = fopen ("diff", "w");
  Vector *denom;

  denom = make_denom (vector->n, 2, 10);

  for (i=0; i<vector->n; i++) {
    if (denom->x[i] == 0.0) continue;
    if (2*i+1 >= vector->n) break;
    diff = vector->x[2*i+1] - vector->x[i];
    fprintf (fp, "%d\t%g\t%g\n", i, 
	     (diff - conn->rtt) / denom->x[i], denom->x[i]);
  }
  fclose (fp);
}
*/


void process_timing (Conn *conn) {
  Vector *timing = conn->timing;
  Vector *starts = make_vector (5);
  Vector *dt;
  double diff;
  int i;
  FILE *fp;

  fp = open_file ("diff", conn->id);

  add_to_vector (starts, timing->x[0]);
  for (i=1; i<timing->n; i++) {
    diff =  timing->x[i] -  timing->x[i-1];
    if (diff > conn->t3 * 5) {
      add_to_vector (starts, timing->x[i]);
    }
  }
  dt = difference_vector (starts);
  fprintf (fp, "%d\t%lf\n", 0, dt->x[0]);
  for (i=1; i<dt->n; i++) {
    if (dt->x[i] < dt->x[i-1]) break;
    fprintf (fp, "%d\t%lf\n", i, dt->x[i]);
  }
  fclose (fp);
}

// WRAP_IT_UP: print the results of various calculations at the
// end of the simulation

void wrap_it_up (Conn *conn) {
  int i;
  double a, b, kca, cavg, tca, tsc, th, q, t;
  double p, cw, queue, pq;
  Vector *dt, *dr1, *dr2;

  p = 1.0 / kmin;
  th = 1000.0 * (conn->total_recv - conn->start_size) / mss / 
    (gtime - conn->start_time);
  cw = conn->total_cw / conn->total_obs / mss;
  queue = conn->total_queue / conn->total_obs / mss;
  pq = conn->total_pq / conn->pq_obs;

  printf ("%d\t%.3f\t%f\t%.3f\t%.3f\t%.3f\t%.3f\n",
	  conn->id, th, 1.0/conn->rtt, conn->b, cw, queue, pq);

  /*
  printf ("\n");
  printf ("th_simulation = %.2lf\n", th);

  printf ("th_friendly = %.2lf\n", 1300 / conn->rtt / sqrt(p));
  printf ("bdp = %.0lf\n", conn->b);
  printf ("num obs = %.0lf\n", conn->total_obs);
  printf ("avg cw = %.2lf\n", conn->total_cw / conn->total_obs / mss);
  printf ("avg queue = %.2lf\n", conn->total_queue / conn->total_obs / mss);
  */

  // process_timing (conn);

  if (conn->id == 1) return;

  dt = difference_vector (starts);
  dr1 = difference_vector (recv1);
  dr2 = difference_vector (recv2);

  for (i=0; i<dt->n; i++) {
    /*
    printf ("%.1lf\t%.1lf\t%.1lf\t%.1lf\t%.1lf\t%.1lf\t%.3lf\n",
	    starts->x[i], dt->x[i],
	    cw1->x[i]/mss, cw2->x[i]/mss,
	    dr1->x[i]/mss, dr2->x[i]/mss,
	    dr1->x[i] / (dr1->x[i] + dr2->x[i]));
    if (i>0 && dt->x[i] < dt->x[i-1]) break;
    */
  }
  return;

  a = sqrt (2.0/3.0*kmin);
  b = conn->b;

  if (b > a) {
    kca = (b*b - a*a) / 2.0;
    cavg = 2.0/3.0 * (b*b*b - a*a*a) / (b*b - a*a);
    tca = b * kca / (cavg - 0.5);
    tsc = kmin - kca;
    th = conn->bw * (kmin-1) / (tsc + tca);
    q = 2*a - b;
  } else {
    kca = 0.0;
    cavg = 0.0;
    tca = 0.0;
    tsc = kmin - kca;
    th = conn->bw * (kmin-1) / (tsc + tca);
    q = 2*a - b;
  }
  
  calc_cw_dist (conn);
  // printf ("th_analysis = %.3lf\n", th);
  printf ("through.analysis.periodic = %.3lf\n", 
	  avg_through (a, (double)kmin, &t, conn));
  
  printf ("avg_delay = %.3lf\n", exp (conn->mu + conn->sigma*conn->sigma/2));
  printf ("a = %.3lf\n", a);
  printf ("b = %.3lf\n", b);
  printf ("kca = %.3lf\n", kca);
  printf ("cavg = %.3lf\n", cavg);
  printf ("tca = %.3lf\n", tca);
  printf ("real time = %.3lf\n", (tca + q) * conn->rtt / b);
  // printf ("avg throughput = %.2lf\n", 1000.0 * conn->total_size / mss / gtime);
  
  // assert (num_drop == num_tripdup);

}

//------------------------------------------------------------------------


// HANDLE_RECV: a data packet arrives at the receiver

void handle_recv (Event *e) {
  Conn *conn = e->conn;
  double arr;
  int n;
  double diff;
  Vector *timing = conn->timing;

  if (verbose) {
    printf ("%.3lf Received packet %d\n", gtime, e->seqno);
  }

  if (e->seqno == conn->start_seqno) {
    conn->start_time = gtime;
    conn->start_size = conn->total_recv;
    if (verbose) {
      printf ("%.3lf start time\n", gtime);
    }
  }

  if (conn->total_recv == -1.0) {
    conn->total_recv = 0.0;
  } else {
    conn->total_recv += e->size;
  }
  fprintf (conn->tfp, "%.3lf\t%.0lf\n", gtime, conn->total_recv/mss);
  add_to_vector (conn->timing, gtime);

  n = timing->n-1;
  if (n == 0) {
    diff = timing->x[n];
  } else {
    diff = timing->x[n] - timing->x[n-1];
  }
  if (conn->id==2 && diff > conn->t3 * 6) {
    add_to_vector (starts, gtime);
    add_to_vector (cw1, conn->other->cw);
    add_to_vector (cw2, conn->cw);
    add_to_vector (recv1, conn->other->total_recv);
    add_to_vector (recv2, conn->total_recv);
  }

  // if the connection was already done, just return
  if (conn->done) return;

  // if this connection is done, wrap it up
  if (conn->total_recv >= conn->total_size) {
    conn->done = 1;
    wrap_it_up (conn);
    done_obs = 1;

    // if this is the last connection, end the simulation
    finished++;
    if (finished == started) {
      exit (0);
    } else {
      return;
    }
  }
  
  // when will the ACK reach the sender?
  arr = gtime + conn->t1 + rand_lognormal (conn->mu, conn->sigma, conn);
  if (e->seqno == ack_delay_seqno) {
    arr += delay;
  }

  if (arr <= conn->last_arr2) {
    arr = conn->last_arr2 + .001;
  }
  conn->last_arr2 = arr;

  e->type = SENDER;
  e->time = arr;
  heap_add (h, e);
}

// SEND_PACKET: invoked by handle_sender whenever it is deemed that
// a packet will be transmitted

Event *send_packet (Conn *conn) {
  Event *e;
  Event *td;
  double arr;
  int suspect;

  fprintf (conn->sfp, "%.3lf\n", gtime);

  e = make_event (BOTTLE, gtime, conn);
  e->seqno = conn->seqno++;

  // does this packet get delayed?
  if (e->seqno == first_delay_seqno) {
    conn->start_time = gtime;
    conn->start_size = conn->total_recv;
    if (verbose) {
      printf ("%.3lf start time\n", gtime);
    }
  }

  // has this packet been marked to be a tripdup?
  td = heap_peek (conn->tripdups);
  if (td != NULL) {
    suspect = (int) td->time;
    if (e->seqno == suspect) {
      td = heap_remove (conn->tripdups);
      // printf ("removing seqno %d from the heap\n", e->seqno);
      e->tripdup = 1;
    }
  }

  // will this packet get dropped?
  if (e->seqno == drop1 ||
      (randdrop && (e->seqno > drop1 && random_drop(kmin))) ||
      (periodic && (e->seqno > drop1 && (e->seqno-drop1)%kmin == 0))) {

    e->size = 0.0;
  }
  conn->total_sent += e->size;

  // compute the arrival time at the next stop
  arr = gtime + conn->t2 + rand_lognormal (conn->mu, conn->sigma, conn);
  if (e->seqno == data1_delay_seqno) {
    arr += delay;
  }
  if (arr <= last_arr3) {
    arr = last_arr3 + .001;
  }
  last_arr3 = arr;
  e->time = arr;
  heap_add (h, e);

  if (verbose) {
    printf ("%.3lf Sending packet %d (total_sent = %lf)\n",
	    arr, e->seqno, conn->total_sent);
  }
  return e;
}


// HANDLE_SENDER: an ACK arrives at the sender 

void handle_sender (Event *e) {
  Conn *conn = e->conn;
  Event *e2;
  double inflight;
  int total = 0;

  // initialize the congestion window
  // or update last_ack
  if (e->seqno == 0) {
    conn->cw = mss;
  } else {
    conn->last_ack = e->seqno;
  }
  if (0) {
    printf ("%.3lf ack of %d, seqno is %d, last_ack = %d\n",
	    gtime, e->seqno, conn->seqno, conn->last_ack);
  }

  // if this packet was not dropped, update cw
  if (e->size > 0) {
    if (conn->slow_start) {
      conn->cw += mss;
    } else {
      conn->cw += mss * mss / conn->cw;
    }
    if (conn->cw > conn->ssthresh) conn->slow_start = 0;
  }

  // if the packet was dropped, don't increment size, and
  // go mark the subsequent packet that will be the tripdup
  if (e->size == 0) {
    num_drop++;
    heap_mark_tripdup (h, e->seqno+3, conn);
  }

  // if this packet is a tripdup, cut the congestion window
  if (e->tripdup) {
    if (++num_tripdup == 5) {
      conn->start_time = gtime;
      conn->start_size = conn->total_recv;
      if (verbose) {
	printf ("%.3lf start time\n", gtime);
      }
    }
    if (verbose) {
      printf ("%.3lf cw cut from %.2lf to %.2lf (avg = %.2lf, q = %.2lf)\n", 
	      gtime, conn->cw/mss, conn->cw/2/mss, conn->cw*7/9/mss, 
	      conn->queue/mss);
    }
    // fprintf (conn->distcwfp, "%.3lf\t%.3lf\n", gtime, conn->cw/mss);
    fprintf (conn->distcwfp, "%.3lf\n", conn->cw/mss);
    conn->cw /= 2;
    conn->slow_start = 0;
  }

  fprintf (conn->cwfp, "%.3lf\t%.3lf\n", gtime, conn->cw/mss);

  // how much data is in flight?
  inflight = mss * (conn->seqno - conn->last_ack - 1);
  assert (inflight >= 0);
  fprintf (conn->ffp, "%.3lf\t%.3lf\n", gtime, inflight/mss);

  free (e);

  // send as much data as possible
  while (inflight + mss <= conn->cw) {
    if (conn->total_sent >= conn->total_size) break;

    e2 = send_packet (conn);
    total++;
    inflight = mss * (conn->seqno - conn->last_ack - 1);
    assert (inflight >= 0);
  }

  fprintf (conn->ffp, "%.3lf\t%.3lf\n", gtime, inflight/mss);
  // printf ("%.1lf sent %d\n", gtime, total);
}


// HANDLE_BOTTLE: a packet arrives at the bottleneck

void handle_bottle (Event *e) {
  Conn *conn = e->conn;
  double arr, dep;

  // update the queue
  fprintf (conn->seenfp, "%.3lf\t%.3lf\n", gtime, conn->queue/mss);
  conn->queue += mss;
  fprintf (conn->queuefp, "%.3lf\t%.3lf\n", gtime, conn->queue/mss);

  // when does this packet depart (begin to transmit)?
  dep = gtime;
  if (dep <= last_dep3 + conn->t3) {
    dep = last_dep3 + conn->t3;
  }
  assert (dep > last_dep3);
  last_dep3 = dep;

  // when does this packet arrive at LAST?
  arr = dep + conn->t3;

  e = remake_event (e, LAST, arr);
  heap_add (h, e);
}


// HANDLE_LAST: a data packet arrives at the home stretch

void handle_last (Event *e) {
  Conn *conn = e->conn;
  double arr;

  // update queue
  conn->queue -= mss;
  fprintf (conn->queuefp, "%.3lf\t%.3lf\n", gtime, conn->queue/mss);

  // when will this packet get to the receiver?
  arr = gtime + conn->t4 + rand_lognormal (conn->mu, conn->sigma, conn);

  // does this packet get a special delay?
  if (e->seqno == data2_delay_seqno) {
    arr += delay;
  }

  // check FIFO requirement
  if (arr <= conn->last_arr1) {
    arr = conn->last_arr1 + .001;
  }
  conn->last_arr1 = arr;

  e->type = RECV;
  e->time = arr;
  heap_add (h, e);
}

// HANDLE_OBS: a data packet arrives at the home stretch

void handle_obs (Event *e) {
  Conn *conn = e->conn;
  Conn *c1, *c2;
  double pq;

  // printf ("%d obs %lf at %lf\n", conn->id, conn->total_obs, gtime);
  if (done_obs) {
    return;
  }

  conn->total_obs += 1.0;
  conn->total_cw += conn->cw;
  conn->total_queue += conn->queue;

  if (conn->id == 1) {
    c1 = conn;
    c2 = conn->other;
  } else {
    c2 = conn;
    c1 = conn->other;
  }
  if (c1->queue + c2->queue > 0.0) {
    pq = c1->queue / (c1->queue + c2->queue);
    conn->pq_obs += 1.0;
    conn->total_pq += pq;
  }

  //if (gtime > 20000.0) return;

  if (verbose) {
    printf ("%.3lf obs when cw = %.3lf, queue = %.3lf\n",
	    gtime, conn->cw, conn->queue);
  }

  e->time = gtime + rand_exp(250.0) + rand_exp(250.0);

  // the following are for equally-spaced queue observations
  // fprintf (conn->queuefp, "%.3lf\t%.3lf\n", gtime, conn->queue/mss);
  // e->time = gtime + 2.718;

  heap_add (h, e);
}

// HANDLE_EVENT: dispatch events to their handlers

void handle_event (Event *e) {

  gtime = e->time;

  switch (e->type) {
  case RECV:
    handle_recv (e);
    break;
  case SENDER:
    handle_sender (e);
    break;
  case BOTTLE:
    handle_bottle (e);
    break;
  case LAST:
    handle_last (e);
    break;
  case OBS:
    handle_obs (e);
    break;
  }
}

void start_conn (Conn *conn, double time)
{
  Event *e = make_event (RECV, time, conn);
  heap_add (h, e);
  started++;
}

void check_heap(Conn *conn)
{
  int i;
  Event *e;

  // if there are events for this connection in the queue, return
  for (i=1; i<h->next; i++) {
    e = heap_cargo (h, i);
    if (e->type != OBS && e->conn == conn) return;
  }

  // otherwise, simulate an event in 1 second
  // printf ("timeout!\n");
  e = make_event (SENDER, gtime+1000.0, conn);
  heap_add (h, e);
}

int main (int argc, char *argv[]) {
  Event *e;
  int c;
  double mult = 1.0;
  double q1, q2, ratio;
  double ssthresh_ratio = 1.0;
  double stagger;

  conn = make_conn ();
  opterr = 0;
  while ( (c = getopt (argc, argv, "i:l:h:x:f:m:n:g:a:o:p:q:d:s:t:k:rv")) != -1) {
    switch (c) {
    case 'd':
      drop1 = atoi (optarg);
      break;
    case 'f':
      first_delay_seqno = atoi (optarg);
      break;
    case 'k':
      kmin = atoi (optarg);
      break;
    case 't':
      conn->total_size = atoi (optarg) * mss;
      break;
    case 's':
      conn->ssthresh = atoi (optarg) * mss;
      break;
    case 'i':
      ssthresh_ratio = atof (optarg);
      break;
    case 'h':
      conn->start_seqno = atoi (optarg);
      break;
    case 'r':
      randdrop = 1;
      periodic = 0;
      break;
    case 'v':
      verbose = 1;
      break;
    case 'a':
      ack_delay_seqno = atoi (optarg);
      break;
    case 'o':
      data1_delay_seqno = atoi (optarg);
      break;
    case 'p':
      data2_delay_seqno = atoi (optarg);
      break;
    case 'q':
      delay = atof (optarg);
      break;
    case 'm':
      conn->mu = atof (optarg);
      no_delay = 0;
      break;
    case 'g':
      conn->sigma = atof (optarg);
      no_delay = 0;
      break;
    case 'l':
      mult = atof (optarg);
      break;
    case 'x':
      rseed = atoi (optarg);
      break;
    }
  }

  starts = make_vector (64);
  cw1 = make_vector (64);
  cw2 = make_vector (64);
  recv1 = make_vector (64);
  recv2 = make_vector (64);

  srandom(rseed);
  h = make_heap (128);

  conn2 = copy_conn (conn);
  conn3 = copy_conn (conn);
  conn->other = conn2;
  conn2->other = conn;

  mult = (conn->rtt * mult - conn->t3) / (conn->rtt - conn->t3);
  conn2->t1 *= mult;
  conn2->t2 *= mult;
  conn2->t4 *= mult;

  //conn2->total_size *= 2.0;

  conn2->ssthresh *= ssthresh_ratio;
  calc_rtt (conn2);
  stagger = 0.001;

  // uncomment this to make share4
  stagger = 1000.0;

  // Uncomment these to make the three connection figure
  conn2->ssthresh = 15 * mss;
  conn3->ssthresh = 5 * mss;

  start_conn (conn, 0.0);
  start_conn (conn2, stagger);
  start_conn (conn3, stagger*2);
  
  q1 = conn->ssthresh - conn->b*mss;
  q2 = conn2->ssthresh - conn2->b*mss;
  q1 = conn->ssthresh;
  q2 = conn2->ssthresh;
  ratio = q1 / (q1 + q2);
  // printf ("predicted throughput = %.2lf\n", conn->bw * ratio);

  /*
  e = make_event (OBS, 10000.0, conn);
  heap_add (h, e);
  e = make_event (OBS, 10000.0, conn2);
  heap_add (h, e);
  */

  while (1) {
    check_heap(conn);
    check_heap(conn2);
    e = heap_remove (h);
    // printf ("%d\n", e->type);
    handle_event (e);
  }

  return 0;
}