/* -*-	Mode:C++; c-basic-offset:8; tab-width:8; indent-tabs-mode:t -*- */
/*
 * Copyright (c) 2001 University of Mannheim, Praktische Informatik IV
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *	This product includes software developed by Praktische Informatik IV,
 *      University of Mannheim.
 * 4. Neither the name of Praktische Informatik IV nor of University of 
 *    Mannheim may be used to endorse or promote products derived from this
 *    software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY U. MANNHEIM AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL U. MANNHEIM OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <math.h>
 
#include "pcc-sink.h"
#include "formula.h"

#define PRINT_JCW 1
#define PRINT_JPD_PROB 0
#define PRINT_JPD 0
#define PRINT_JPD_VERBOSE 0
#define PRINT_NUMPACK_JPD 0

static class PccSinkClass : public TclClass {
public:
  	PccSinkClass() : TclClass("Agent/PCCSink") {}
  	TclObject* create(int, const char*const*) {
     		return (new PccSinkAgent());
  	}
} //--		(now - last_report_sent >= rtt_/NumFeedback_)) {
class_pccSink; 


PccSinkAgent::PccSinkAgent() : Agent(PT_PCC_ACK), nack_timer_(this),
			       reexp_timer_(this), off_timer_(this)
{
    bind("packetSize_", &size_);  //-- not used
	bind("InitHistorySize_", &hsz);
	bind("NumFeedback_", &NumFeedback_);
	bind ("printLoss_", &printLoss_);

	bind("df_", &df_);
	bind("bval_", &bval_);

	// for WALI ONLY
	bind ("NumSamples_", &numsamples);
	bind ("discount_", &discount);
	bind ("smooth_", &smooth_);

	bind ("totalReceived_", &total_received_);

	//PCC-specific
	bind ("protLossEvents_", &protLossEvents);
	bind ("protRTTs_", &protRTTs);
	bind ("protMaxTime_", &protMaxTime);
	bind ("Toff_", &Toff);
	bind ("Texp_", &Texp);
//--	bind ("rateDiscount_", &rateDiscount);
	bind ("adjustRate_", &adjust_rate);

	bind ("runtime_", &runtime_);

	bind ("pcc_recv_nr_", &pcc_recv_nr_);
	bind ("always_on_", &always_on_);

	rate_ = 0;
	rtt_ =  0; 
//--        sqrtrtt_=1; // is 1 a good default for sqrtrtt? (probably better than 0)
	reported_rtt_ =  0; 
	last_timestamp_ = 0;
	last_arrival_ = 0;
	last_report_sent=0;

	maxseq = -1;
	lastloss = 0;
	lastloss_round_id = -1 ;
        round_id = 0;
	UrgentFlag = 0;

	rtvec_ = NULL;
	tsvec_ = NULL;
	lossvec_ = NULL;

	// used by WALI and EWMA
	last_sample = 0;

	// used only for WALI 
	sample = NULL ; 
	weights = NULL ;
	mult = NULL ;
	sample_count = 1 ;
	mult_factor_ = 1.0;
	init_WALI_flag = 0;

	//PCC
	wsize = (int)(Toff/Texp); //set size of pwindow to the number of Texp in Toff
	if (Toff/Texp>wsize) wsize++; //if integer division isn't exact, add one
	pwindow = (double *)malloc(sizeof(double)*wsize);
	p_ts_window = (double *)malloc(sizeof(double)*wsize);
	init_pcc_values();

	noAckCount = 0;
	flowID = (int)(Random::uniform() * 100);
}

/*
 * Receive new data packet.  If appropriate, generate a new report.
 */
void PccSinkAgent::recv(Packet *pkt, Handler *)
{
	hdr_pcc *pcch = hdr_pcc::access(pkt); 
	double now = Scheduler::instance().clock();
	//double p = -1;

	/* update the round trip time */
	if (pcch->timestamp_echo > 0) {  //only update rtt if there is a valid timestamp echo
	  if (PRINT_JPD_VERBOSE) printf("((R%i)) updating rtt at %f with echo %f, offset %f\n",
					flowID, now, pcch->timestamp_echo, pcch->timestamp_offset);
	  update_rtt (pcch->timestamp_echo + pcch->timestamp_offset, now);
	}

	if (on && protect) {  //eventuell reicht hier protect, da protect=>on
	  if (protected_start <= 0) {  //first packet of protected time
	    protected_start = now;
	  } else {
	    if (pcch->round_id > round_id)
	      protected_rounds++;
	    protected_time = now - protected_start; //we compare protected_time to protMaxTime below
	  }
	}

	//if we receive packets while off, accept them for 2 RTTs (to give the sender time
	// to react to our off-packet), ignore (->discard) all packets after that
	//discarding packets is meant as an incentive to senders to adhere to the protocol
	//rcv_while_off is necessary since rtt_ can change over time so that without it,
	// one packet might be discarded and one with a higher sequence number received
	// (because of an increase in rtt_ between the two packets)
	//an alternative to rcv_while_off would be making 2*rtt_ fix (e.g. as rtt_at_off),
	// however, this might result in packets from a sender adhering to the protocol
	// being dropped because of a sharp increase in rtt_
	//this is of course still possible because of the smoothing of rtt_, hence the
	// use of 2*rtt_ instead of rtt_
	if (on || (((now - exp_failed_time) <= 2*rtt_) && rcv_while_off)) {
	  if (maxseq < 0) init_vectors(pcch, now);

	  /* for the time being, we will ignore out of order and duplicate 
	     packets etc. -> (JPD) round_id has been secured against packet reordering */
          int seqno = pcch->seqno ;
	  if (seqno > maxseq) {
	    total_received_++;
	    UrgentFlag = 0;
	    //only update round_id if the new one is larger
	    //this ensures that round_id increases monotonically, thus we will not
	    // count extra rounds if packets are reordered
	    if (pcch->round_id > round_id) round_id = pcch->round_id;
	    psize_=pcch->psize;
	    sendrate=pcch->sendrate;
	    last_arrival_=now;
	    last_timestamp_=pcch->timestamp;
	    rtvec_[seqno%hsz]=now;
	    tsvec_[seqno%hsz]=last_timestamp_;	
	    lossvec_[seqno%hsz] = RCVD;
	    int i = maxseq+1 ;
	    if (i < seqno) {
	      double delta = (tsvec_[seqno%hsz]-tsvec_[maxseq%hsz])/(seqno-maxseq) ; 
	      double tstamp = tsvec_[maxseq%hsz]+delta ;
	      while(i < seqno) {
		rtvec_[i%hsz]=now;	
		tsvec_[i%hsz]=tstamp;	
		if ((tsvec_[i%hsz]-lastloss > rtt_) && 
		    (pcch->round_id > lastloss_round_id)) {
		  // Lost packets are marked as "LOSS_EVENT"
		  // at most once per RTT.
		  lossvec_[i%hsz] = LOSS_EVENT;
		  if (PRINT_JPD) printf("((R%i)) New loss event!\n", flowID);
		  if (PRINT_JPD_PROB && now>=200 && now<=runtime_)
		    printf("lossevent\n");
		  UrgentFlag = 1 ; //(JPD) immediately send feedback on new loss event
		  lastloss = tstamp;
		  if (pcch->round_id > lastloss_round_id)
		    lastloss_round_id = pcch->round_id;
		  if (on && protect) {  //(see above)
		    protected_losses++;
		  }
		}
		else {
		  // This lost packet is marked "LOSS" (instead of "LOSS_EVENT")
		  // because it does not begin a loss event.
		  lossvec_[i%hsz] = LOSS; 
		}
		i++;
		tstamp = tstamp+delta;
	      }
	    }
//--		int x = maxseq; ??
	    maxseq = pcch->seqno ;
//--	    if ((rtt_ > SMALLFLOAT) && 
//--		(now - last_report_sent >= rtt_/NumFeedback_)) {
	    if (now - last_report_sent >= rtt_/NumFeedback_
		|| rtt_ >= 2*reported_rtt_ || rtt_ <= 0.5*reported_rtt_) {
	      //we need to send a feedback packet before we have a valid rtt
	      // so that we get a valid echo from which to calculate that rtt
	      UrgentFlag = 1;
	      if (rtt_ <= 0) rtt_ = 1;  //set rtt_ to some decent default value
	      reported_rtt_ = rtt_;
	    }
	    if (on && protect) {
	      if ((protected_losses >= protLossEvents && protected_rounds >= protRTTs)
		  || protected_time > protMaxTime) {
		if (PRINT_JPD_VERBOSE) printf("((R%i)) protected_losses: %i (%i), protected_rounds: %i ", flowID, protected_losses, protLossEvents, protected_rounds);
		//conduct the experiment
		experiment();
		protect = 0;
	      }
	    }
	    if (UrgentFlag) {
	      UrgentFlag = 0 ;
	      nextpkt();
	    }
	  }
	} else if (rcv_while_off && ((now - exp_failed_time) > 2*rtt_)) {
	  rcv_while_off = 0;
	  discard_after = maxseq;
	}
	Packet::free(pkt);
}

/*
 * Conduct the experiment to see if the flow may remain on.
 */
void PccSinkAgent::experiment()
{
	double now = Scheduler::instance().clock();

	double prob;  //probability to stay on
	double loss;
	double tcprate;
	double rand;
	int failed = 0;
	double temp_toff = Toff;

	loss = est_loss();
	//avoid division by 0 in case protected time ended with no losses because of protMaxTime
	if (loss > 0) { 
	  tcprate = p_to_b(loss, rtt_, 4*rtt_, psize_, bval_);
	  if (adjust_rate) {
	    //use model from paper (Ramesh/Rhee: "Issues...") to compute estimate for
	    // loss event rate based on (overly aggressive) tcprate as sendrate
	    loss = (1 - pow((1-lossrate),(int)(rtt_*tcprate/psize_))) / (rtt_*tcprate/psize_);
	    //next, recompute tcprate with this loss event rate estimate
	    tcprate = p_to_b(loss, rtt_, 4*rtt_, psize_, bval_);
	  }
	}
	else tcprate = MAXRATE;

    tcprate = 62500; // JCW

//--	est_rate();  //update the sendrate estimation --> done by sender
	if (PRINT_JPD) printf("((R%i)) sendrate: %f; tcprate: %f\n", flowID, sendrate, tcprate);
	if (PRINT_JPD_PROB && now>=200 && now<=runtime_ && tcprate<MAXRATE) {
	  printf("rates %f\n", tcprate);
	  printf("losseventrate %f\n", loss);
	}

	if (protect) {  //first experiment after protected time - set to 0 in recv
	  if (PRINT_JPD) printf("((R%i)) protected time ends after %f; ", flowID, protected_time);
	  if (PRINT_JPD_PROB && now>=200 && now<=runtime_) printf("\nprottime %f\n", protected_time);
	  prob = ((Toff+protected_time)*tcprate - protected_time*sendrate)/(Toff*sendrate);
	}
	// this usually works, but isn't in compliance with theory: non-ideal formula needs to be
	// used for all calculations up to the end of the first off time! Also, for the send
	// rate below the fraction bar (only!), the adjusted (effective) rate needs to be used

	else {  //ongoing evaluation
	  //if we have made other experiments within the last Toff, we adjust
	  // the sendrate to a fictive value representing this flow's expected
	  // average sending rate (probabilistically), by multiplying it with the
	  // prob values ( <= 1 ) that occurred in that time
	  double tmprate = sendrate;
	  for (int i=0; i<wsize && pwindow[i]>0 && p_ts_window[i]>0
//--		 && (p_ts_window[i]-now)<=(Toff/Texp); i++) {
		 && (now-p_ts_window[i])<=Toff; i++) {
	    tmprate*=pwindow[i];
	  }
	  prob = tcprate/tmprate;
	}
	if (PRINT_JPD) printf("prob = %f; ", prob);

    if (always_on_) // JCW
        prob = 1;

	if (prob < 0) {  //turn flow off, adjust (temporaray) Toff
	  //note: prob < 0 can only happen directly after protected time (otherwise, prob is
	  // calculated as tcprate/sendrate, neither of which can be negative)
	  failed = 1;
	  //instead of using the parameter Toff, schedule off_timer such that
	  // with this value, prob would have been exactly 0 (not <0)
	  //we don't adjust the actual Toff-parameter (yet); this may be an option
	  // though to prevent having p<0 too often ***
	  temp_toff = protected_time * (sendrate-tcprate) / tcprate;
	  if (PRINT_JPD_PROB && now>=200 && now<=runtime_) printf ("\nsmallerzero\n");
	} else if (prob < 1) {  // 0 <= prob < 1 => conduct the random experiment
	  rand = Random::uniform();
	  if (PRINT_JPD) printf("rand = %f; ", rand);
	  if (PRINT_JPD_PROB && now>=200 && now<=runtime_) printf ("\navgrand %f %f\n", prob, rand);
	  if (prob <= rand) {  //we fail -> turn flow off (leave temp_toff at Toff)
	    failed = 1;
	  } else {
	    if (protect) {
	      //after using the adjusted value for the experiment, reset prob to tcprate/sendrate
	      //thus, if all conditions are equal at the next experiment, prob will be exactly one
	      // (and not >1 as it would be with the adjusted value); more importantly, if
	      // conditions get slightly worse, prob will be <1
	      prob = tcprate/sendrate;
	    }
	    //update probabilities window
	    shift_array(pwindow, wsize, prob);
	    shift_array(p_ts_window, wsize, now);
	    if (PRINT_JPD) {
	      printf("pwindow = [");
	      for (int i = 0; i<wsize && pwindow[i]>0; i++) {
		printf("%f,", pwindow[i]);
	      }
	      printf("]\n");
	    }
	  }
	} else {  //prob >= 1 -> enter a 1 into the array to shift it in a neutral way
	         // *** could be removed for better performance, since we ckeck timestamp
	        //  when multiplying with array values
	  shift_array(pwindow, wsize, 1);
	  shift_array(p_ts_window, wsize, now);
	  if (PRINT_JPD_PROB && now>=200 && now<=runtime_) printf ("\ngreaterone\n");
	}

    if (PRINT_JCW) {
	  printf("%i %f %f %f %f %f %i %f\n", pcc_recv_nr_, now, tcprate, loss, rtt_, prob, failed, tcprate/sendrate);
    }

	if (failed) {
	  exp_failed_time = now;
	  on = 0;
	  UrgentFlag = 1;
	  rcv_while_off = 1;
	  off_timer_.resched(temp_toff); //schedule off_timer
	  return;
	}
	else {
	  //reschedule timer for next experiment (application parameter Texp)
	  reexp_timer_.resched(Texp);
	  return;
	}  
}

/*
 * Schedule sending this report, and set timer for the next one.
 */
void PccSinkAgent::nextpkt()
{
	sendpkt();

	double now = Scheduler::instance().clock();

	//reschedule feedback if the flow is on and for 4 RTTs afterwards
	//(send several NACKs in case some get lost)
	if (on || ((now - exp_failed_time) <= 4*rtt_)) {
	  /* schedule next report 1.5*rtt/NumFeedback later */
	  if (rtt_ > 0.0 && NumFeedback_ > 0) 
//--	    nack_timer_.resched(1.5*rtt_/NumFeedback_); //*** warum 1.5?
	    nack_timer_.resched(rtt_/NumFeedback_);
	} else if (PRINT_JPD) printf("((R%i)) Not rescheduling NACKs at %f\n", flowID, now);
}

/*
 * Create report message, and send it.
 */
void PccSinkAgent::sendpkt()
{
	double now = Scheduler::instance().clock();

//--	if (PRINT_JPD_VERBOSE) {
//--	  /* update estimate of fair rate */
//--	  double loss = est_loss();
//--	  printf("((R%i)) loss = %f;", flowID, loss);
//--	  rate_ = p_to_b(loss, rtt_, 4 * rtt_, psize_, bval_);
//--	  printf(" friendly rate = %f; ", rate_);
//--	}
	//check if the sender is currently allowed to send
	if (on) rate_ = MAXRATE;
	else rate_ = 0;
	if (PRINT_JPD) printf("rate = %f\n", rate_);

	/*don't send an ACK unless we've received new data*/
	/*if we're sending slower than one packet per RTT, don't need*/
	/*multiple responses per data packet.*/
	//do send an ACK every 10th time to ensure that at a high loss rate
	// the ACK still gets through at some time

	noAckCount = (noAckCount++) % 10;
	if (last_arrival_ >= last_report_sent || noAckCount == 9) {

		Packet* pkt = allocpkt();
		if (pkt == NULL) {
			printf ("error allocating packet\n");
			abort(); 
		}
	
		hdr_pcc_ack *pcc_ackh = hdr_pcc_ack::access(pkt);
	
		pcc_ackh->timestamp=now;
		pcc_ackh->rate = rate_;
		pcc_ackh->rtt = rtt_;
		if (PRINT_JPD_PROB && now>=200 && now<=runtime_ && (int)(now*100)%100==0)
		  printf("rtt %f\n", rtt_);
		//		if (discard_after >= 0)
		pcc_ackh->discard_after_seqno = discard_after;
		last_report_sent = now; 
		send(pkt, 0);
		if (PRINT_JPD) printf(" packet sent!\n");
		if (PRINT_NUMPACK_JPD) printf("Received %i packets, announced rate: %f\n", total_received_, rate_);
	} else if (PRINT_JPD_VERBOSE) printf("((R%i)) Not sending NACK because no new data has arrived\n", flowID);
}

void PccSinkAgent::update_rtt (double tao, double now) 
{
	/* fine grained RTT estimate for use in the equation */
	if (rtt_ > 0) {
		rtt_ = df_*rtt_ + ((1-df_)*(now - tao));
//--		sqrtrtt_ = df_*sqrtrtt_ + ((1-df_)*sqrt(now - tao));
	} else {
		rtt_ = now - tao;
//--		sqrtrtt_ = sqrt(now - tao);
	}
}

//--/*
//-- * compute smoothed estimate of the sender's rate
//-- */
//--void PccSinkAgent::est_rate ()
//--{
//--        double last = rtvec_[maxseq%hsz]; 
//--        int rcvd = 0;
//--        int i = maxseq;
//--	double rate_sample = 0;
//--	double interval = 0;
//--
//--	if (sendrate <= 0) {  //first rate estimation after protected time
//--	  //average sending rate over all of protected time (minus one rtt,
//--	  // to be on the safe side) to get an initial EWMA value
//--	  interval = protected_time - rtt_;
//--	} else {  //average rate over the last 2 RTTs
//--	  interval = 2*rtt_;
//--	}
//--	//count packets received or lost during the interval
//--	while (i > 0 && (last - rtvec_[i%hsz] < interval)) {
//--	  if (lossvec_[i%hsz] != UNKNOWN) {
//--	    rcvd++;
//--	    i--;
//--	  }
//--	}
//--
//--        if ((last - rtvec_[i%hsz]) > 0) {
//--	  rate_sample = rcvd*psize_/(last - rtvec_[i%hsz]);
//--        } else {
//--	  rate_sample = 0;
//--        }
//--
//--	if (sendrate <= 0)
//--	  sendrate = rate_sample;
//--	else
//--	  sendrate = rateDiscount * sendrate + (1-rateDiscount) * rate_sample;
//--	return;
//--}

int PccSinkAgent::command(int argc, const char*const* argv) 
{
	if (argc==2) {
		if (strcmp(argv[1],"stop")==0) {
			stop();
			return TCL_OK;
		}
	} else if (argc == 3) {
		if (strcmp(argv[1], "weights") == 0) {
			/* 
			 * weights is a string of numbers, seperated by + signs
			 * the first number is the total number of weights.
			 * the rest of them are the actual weights
			 * this overrides the defaults
			 */
			char *w ;
			w = (char *)calloc(strlen(argv[2])+1, sizeof(char)) ;
			if (w == NULL) {
				printf ("error allocating w\n");
				abort();
			}
			strcpy(w, (char *)argv[2]);
			numsamples = atoi(strtok(w,"+"));
			sample = (int *)malloc((numsamples+1)*sizeof(int));
			weights = (double *)malloc((numsamples+1)*sizeof(double));
			mult = (double *)malloc((numsamples+1)*sizeof(double));
			fflush(stdout);
			if (sample && weights) {
				int count = 0 ;
				while (count < numsamples) {
					sample[count] = 0;
					mult[count] = 1;
					char *w;
					w = strtok(NULL, "+");
					if (w == NULL)
						break ; 
					else {
						weights[count++] = atof(w);
					}	
				}
				if (count < numsamples) {
					printf ("error in weights string %s\n", argv[2]);
					abort();
				}
				sample[count] = 0;
				weights[count] = 0;
				mult[count] = 1;
				free(w);
				return (TCL_OK);
			}
			else {
				printf ("error allocating memory for smaple and weights:2\n");
				abort();
			}
		}
	}
	return (Agent::command(argc, argv));
}

void PccSinkAgent::stop()
{
  double now = Scheduler::instance().clock();
  if (PRINT_JPD_VERBOSE) printf("Receiver stopped at %f\n", now);
  nack_timer_.force_cancel();
  reexp_timer_.force_cancel();
  off_timer_.force_cancel();
  on = 0;               //make sure no more NACKs
  exp_failed_time = 0; // are scheduled
}

void PccNackTimer::expire(Event *) {
  if (PRINT_JPD_VERBOSE) printf("((R%i)) NackTimer expired.\n", a_->flowID);
  a_->nextpkt();
}

void PccReExpTimer::expire(Event *) {
  if (PRINT_JPD_VERBOSE) printf("((R%i)) ReExpTimer expired.\n", a_->flowID);
  a_->experiment();
}

void PccOffTimer::expire(Event *) {
  if (PRINT_JPD_VERBOSE) printf("((R%i)) OffTimer expired.\n", a_->flowID);
  a_->init_pcc_values(); //reset values to on and protected (etc.)
  a_->nextpkt();        //restart sending acks so the sender can resume
}

void PccSinkAgent::print_loss(int sample, double ave_interval)
{
	double now = Scheduler::instance().clock();
	double drops = 1/ave_interval;
	printf ("time: %7.5f sample 0: %5d loss_rate: %7.5f \n", 
		now, sample, drops);
	//printf ("time: %7.5f send_rate: %7.5f\n", now, sendrate);
	//printf ("time: %7.5f maxseq: %d\n", now, maxseq);
}

void PccSinkAgent::print_loss_all(int *sample) 
{
	double now = Scheduler::instance().clock();
	printf ("%f: sample 0: %5d 1: %5d 2: %5d 3: %5d 4: %5d\n", 
		now, sample[0], sample[1], sample[2], sample[3], sample[4]); 
}

/// WALI Code
double PccSinkAgent::est_loss() 
{
	double now = Scheduler::instance().clock();
	int i;
	double ave_interval1, ave_interval2; 
	int ds ; 

	if (!init_WALI_flag) {
		init_WALI () ;
	}
	// sample[i] counts the number of packets since the i-th loss event
	// sample[0] contains the most recent sample.
	for (i = last_sample; i <= maxseq ; i ++) {
		sample[0]++;
		if (lossvec_[i%hsz] == LOSS_EVENT) {
		        //  new loss event
			// double now = Scheduler::instance().clock();
			sample_count ++;
			shift_array (sample, numsamples+1, 0); 
			multiply_array(mult, numsamples+1, mult_factor_);
			shift_array (mult, numsamples+1, 1.0); 
			mult_factor_ = 1.0;
		}
	}

	// calculate loss rate over the same number of packets and with analogue weighting
	int packets_in_sample_array = 0;
	for (int i=0; i<numsamples+1; i++) {
	  packets_in_sample_array += sample[i];
	}
	double numlosses = 0;
	//!!potential error source if maxseq wraps around
	// implicit assumption: maxseq >= packets_in_sample_array (holds for "reasonably"
	// long experiments)
	int middle = maxseq - (packets_in_sample_array / 2);
	for (int i=maxseq; i > (maxseq-packets_in_sample_array); --i) {
	  if (lossvec_[i%hsz] == LOSS || lossvec_[i%hsz] == LOSS_EVENT) {
	    if (i>middle) numlosses += 1.0;
	    else numlosses += (double)(i - maxseq + packets_in_sample_array) / (double)((packets_in_sample_array+1)/2 + 1);
	  }
	}
	lossrate = numlosses / (packets_in_sample_array * 0.75);
	if (PRINT_JPD_PROB && now>=200 && now<=runtime_) printf ("lossrate %f\n", lossrate);
	//--return lossrate; //for experiments using loss rate instead of loss event rate

	if (PRINT_JPD_VERBOSE) {
	  printf("((R%i)) sample[", flowID);
	  for (int i=0; i<numsamples+1; i++) {
	    if (sample[i] > 0) printf("%i,", sample[i]);
	    else {
	      break;
	    }
	  }
	  printf("]\n");
	}
	last_sample = maxseq+1 ; 

	if (sample_count>numsamples+1)
		// The array of loss intervals is full.
		ds=numsamples+1;
    	else
		ds=sample_count;

	if (sample_count == 1) 
		// no losses yet
		return 0; 

	/* do we need to discount weights? */
	if (sample_count > 1 && discount && sample[0] > 0) {
		double ave = weighted_average(1, ds, 1.0, mult, weights, sample);
		int factor = 2;
		double ratio = (factor*ave)/sample[0];
		double min_ratio = 0.5;
		if ( ratio < 1.0) {
			// the most recent loss interval is very large
			mult_factor_ = ratio;
			if (mult_factor_ < min_ratio) 
				mult_factor_ = min_ratio;
		}
	}
	// Calculations including the most recent loss interval.
	ave_interval1 = weighted_average(0, ds, mult_factor_, mult, weights, sample);
	// The most recent loss interval does not end in a loss event.
	// Include the most recent interval in the calculations
	// only if this increases the estimated loss interval.
	ave_interval2 = weighted_average(1, ds, mult_factor_, mult, weights, sample);
	if (ave_interval2 > ave_interval1)
		ave_interval1 = ave_interval2;
	if (ave_interval1 > 0) { 
		if (printLoss_ > 0) {
			print_loss(sample[0], ave_interval1);
			print_loss_all(sample);
		}
		return 1/ave_interval1; 
	} else {
                double avg = 0;
                for (int i = 0 ; i <= numsamples; i++) {
                        fprintf(stderr, "%i ", sample[i]);
                        avg += sample[i];
                }
                fprintf(stderr, " --- avg %f %i\n", avg / (numsamples + 1), numsamples);
                exit(2);
        }
}


// Calculate the weighted average.
double PccSinkAgent::weighted_average(int start, int end, double factor, double *m, double *w, int *sample)
{
	int i; 
	double wsum = 0;
	double answer = 0;
	if (smooth_ == 1 && start == 0) {
		if (end == numsamples+1) {
			// the array is full, but we don't want to uses
			//  the last loss interval in the array
			end = end-1;
		} 
		// effectively shift the weight arrays 
		for (i = start ; i < end; i++) 
			if (i==0)
				wsum += m[i]*w[i+1];
			else 
				wsum += factor*m[i]*w[i+1];
		for (i = start ; i < end; i++)  
			if (i==0)
			 	answer += m[i]*w[i+1]*sample[i]/wsum;
			else 
				answer += factor*m[i]*w[i+1]*sample[i]/wsum;
	        return answer;
	} else {
		for (i = start ; i < end; i++) 
			if (i==0)
				wsum += m[i]*w[i];
			else 
				wsum += factor*m[i]*w[i];
		for (i = start ; i < end; i++)  
			if (i==0)
			 	answer += m[i]*w[i]*sample[i]/wsum;
			else 
				answer += factor*m[i]*w[i]*sample[i]/wsum;
	        return answer;
	}
}

// Shift array a[] up, starting with a[sz-2] -> a[sz-1].
void PccSinkAgent::shift_array(int *a, int sz, int defval) 
{
	int i ;
	for (i = sz-2 ; i >= 0 ; i--) {
		a[i+1] = a[i] ;
	}
	a[0] = defval;
}
void PccSinkAgent::shift_array(double *a, int sz, double defval) {
	int i ;
	for (i = sz-2 ; i >= 0 ; i--) {
		a[i+1] = a[i] ;
	}
	a[0] = defval;
}

// Multiply array by value, starting with array index 1.
// Array index 0 of the unshifted array contains the most recent interval.
void PccSinkAgent::multiply_array(double *a, int sz, double multiplier) {
	int i ;
	for (i = 1; i <= sz-1; i++) {
		double old = a[i];
		a[i] = old * multiplier ;
	}
}

void PccSinkAgent::init_WALI () {
	int i;
	if (numsamples < 0)
		numsamples = DEFAULT_NUMSAMPLES ;	
	if (smooth_ == 1) {
		numsamples = numsamples + 1;
	}
	sample = (int *)malloc((numsamples+1)*sizeof(int));
	weights = (double *)malloc((numsamples+1)*sizeof(double));
	mult = (double *)malloc((numsamples+1)*sizeof(double));
	for (i = 0 ; i < numsamples+1 ; i ++) {
		sample[i] = 0 ; 
	}
	if (smooth_ == 1) {
		int mid = int(numsamples/2);
		for (i = 0; i < mid; i ++) { //--effizienter hier bis <= mid
			weights[i] = 1.0;
		}
		for (i = mid; i <= numsamples; i ++){
			weights[i] = 1.0 - (i-mid)/(mid + 1.0);
		}
	} else {
		int mid = int(numsamples/2);
		for (i = 0; i < mid; i ++) {
			weights[i] = 1.0;
		}
		for (i = mid; i <= numsamples; i ++){
			weights[i] = 1.0 - (i+1-mid)/(mid + 1.0);
		}
	}
	for (i = 0; i < numsamples+1; i ++) {
		mult[i] = 1.0 ; 
	}
	init_WALI_flag = 1;  /* initialization done */
}

/*
 * initialize loss/return-time/timestamp-vectors
 */
void PccSinkAgent::init_vectors(hdr_pcc *pcch, double now)
{
  maxseq = pcch->seqno - 1 ;
  rtvec_=(double *)malloc(sizeof(double)*hsz);
  tsvec_=(double *)malloc(sizeof(double)*hsz);
  lossvec_=(char *)malloc(sizeof(double)*hsz);
  if (rtvec_ && lossvec_) {
    int i;
    for (i = 0; i < hsz ; i ++) {
      lossvec_[i] = UNKNOWN;
      rtvec_[i] = -1; 
      tsvec_[i] = -1; 
    }
    for (i = 0; i <= maxseq ; i++) {
      lossvec_[i] = NOLOSS ; 
      rtvec_[i] = now ;
      tsvec_[i] = last_timestamp_ ;
    }
  }
  else {
    printf ("error allocating memory for packet buffers\n");
    abort (); 
  }
  return;
}

/*
 * initialize values used specifically for PCC
 * this is done at startup as well as every time a sender resumes after an off period
 */
void PccSinkAgent::init_pcc_values() {
  on = 1;
  protect = 1;
  protected_time = 0;
  protected_start = 0;
  protected_losses = 0;
  protected_rounds = 0;
  exp_failed_time = LARGE_DOUBLE; //so that (now - exp_failed_time) < 0 < 3*rtt_
  for (int i=0; i<wsize; i++) pwindow[i] = -1;
  for (int i=0; i<wsize; i++) p_ts_window[i] = -1;
  discard_after = -1;
  sendrate = 0;
  rcv_while_off = 0;
  //off_nack_count = 0;
  //also reset last_report_sent so sendpkt will actually send a packet
  last_report_sent = 0;
  return;
}