/*******************************************************************************
    OpenAirInterface
    Copyright(c) 1999 - 2014 Eurecom

    OpenAirInterface 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 3 of the License, or
    (at your option) any later version.


    OpenAirInterface 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 OpenAirInterface.The full GNU General Public License is
    included in this distribution in the file called "COPYING". If not,
    see <http://www.gnu.org/licenses/>.

  Contact Information
  OpenAirInterface Admin: openair_admin@eurecom.fr
  OpenAirInterface Tech : openair_tech@eurecom.fr
  OpenAirInterface Dev  : openair4g-devel@eurecom.fr

  Address      : Eurecom, Compus SophiaTech 450, route des chappes, 06451 Biot, France.

*******************************************************************************/


/*! \file otg_models.c
* \brief function containing the OTG TX traffic generation functions 
* \author M. Laner and navid nikaein
* \date 2013
* \version 0.1
* \company Eurecom
* \email: navid.nikaein@eurecom.fr
* \note
* \warning
*/

#include "UTIL/MATH/oml.h"
#include "otg_models.h"
#include "UTIL/LOG/log.h"



/*all the TARMA stuff:
	this is based on the modelling framework provided in 
 	"M. Laner, P. Svoboda and M. Rupp, Modeling Randomness in Network Traffic, SIGMETRICS/Performance'12, London, UK, 2012". 

*/

/* this function must be called for each random process  and each time it shall produce a new random sample, 
 which is the return value 
 input:		- inputSamples: an array of length TARMA_NUM_PROCESSES with iid gaussian samples 
 			  which must be the same for xcorrelated processes
 			- proc: pointer to tarmaProcess_t
 output:	- the next value in the random process
*/
double tarmaCalculateSample( double inputSamples[], tarmaProcess_t *proc){  
  /* struct with all definitions of the random process and its history */
  int cnt;
  double X, Y, Ypower, Z;

  /* calculate the current sample X' by weighting the input samples */
  X=0;
  for( cnt=0; cnt<TARMA_NUM_PROCESSES ; cnt++){
	X += proc->inputWeight[cnt]*inputSamples[cnt];
  }

  /* update the MA history: shift everything back for one lag */
  for( cnt=TARMA_NUM_MA_MAX-1; cnt>0; cnt--){
	proc->maHist[cnt] = proc->maHist[cnt-1];
  }
  proc->maHist[0]=X;  /* update the MA history with the current filter input value */

  /* update the AR history: shift everything back for one lag */
  for( cnt=TARMA_NUM_AR_MAX-1; cnt>0; cnt--){
	proc->arHist[cnt] = proc->arHist[cnt-1];
  }
  proc->arHist[0]=0;  /* caution: arHist[0] must not be used for the filter processing! */

  /* filter processing */
  Y=0;
  for( cnt=0; cnt<TARMA_NUM_MA_MAX; cnt++){  /* moving average part, start at lag 0 = current X */
	Y += proc->maWeight[cnt] * proc->maHist[cnt];
  }
  for( cnt=1; cnt<TARMA_NUM_AR_MAX; cnt++){  /* auto regressive part, start at lag 1 = former Y */
	Y += -proc->arWeight[cnt] * proc->arHist[cnt];
  }
  proc->arHist[0]=Y;  /* update the AR history with the current filter output value */

  /* polynomial transformation */
  Z=0;
  Ypower=1;
  for(cnt=0; cnt<TARMA_NUM_POLY_MAX; cnt++){
	Z += proc->polyWeight[cnt]*Ypower;
	Ypower *= Y;
  }

  LOG_D(OTG,"TARMA_DEBUG: tarmaCalculateSamples called: Z=%f\n",Z);
  //tarmaPrintProc(proc);
  return Z;
}

/*
 this function must be called each time before size and itd are generated for the respective packet.
	it updates the input samples stored in the structure stream,
 	ps: calling it more often than for each packet is also fine 
 input: 	- the stream to be updated
 output:  
 */
void tarmaUpdateInputSample (tarmaStream_t *stream){
  int cnt;
  LOG_T(OTG,"TARMA_DEBUG: tarmaUpdateInputSample(%p)\n", stream);
  if(stream){
	for(cnt=0; cnt<TARMA_NUM_PROCESSES; cnt++){
	  stream->tarma_input_samples[cnt]=gaussian_dist(10000,1)-10000;
    LOG_D(OTG,"TARMA_DEBUG:   %f\n",stream->tarma_input_samples[cnt]);
	}
  }
}

/*
 this function initializes the given stream, it allocates the respective memory if necessary (stream=0)
 	and sets all the respective values to zero.
 input: 	- steam to be initialized (or zero if memory shall be allocated)
 output:	- the same stream
 */
tarmaStream_t *tarmaInitStream(tarmaStream_t *stream){
  if(stream==0){
	stream=(tarmaStream_t*) malloc(sizeof(tarmaStream_t));
  }
  tarmaProcess_t *proc;
  int cntvar, cntp, cntma, cntar, cntpy;
  for(cntvar=0; cntvar<2; cntvar++){
	if(cntvar==0){
	  proc=&(stream->tarma_idt);
	}else{
	  proc=&(stream->tarma_size);
	}
	for(cntp=0; cntp<TARMA_NUM_PROCESSES; cntp++){
	  proc->inputWeight[cntp]=0;
	}
	for(cntma=0; cntma<TARMA_NUM_MA_MAX; cntma++){
	  proc->maWeight[cntma]=0;
	  proc->maHist[cntma]=0;
	}
	for(cntar=0; cntar<TARMA_NUM_AR_MAX; cntar++){
	  proc->arWeight[cntar]=0;
	  proc->arHist[cntar]=0;
	}
	for(cntpy=0; cntpy<TARMA_NUM_POLY_MAX; cntpy++){
	  proc->polyWeight[cntpy]=0;
	}
  }
  LOG_D(OTG,"TARMA_DEBUG: tarmaInitStream(%p) called\n", stream);
  return stream;
}

/* 
 this function initializes the input stream according to the openarena traffic modelling
 input:		- stream to be initialized
 output: 	
 */
void tarmaSetupOpenarenaDownlink(tarmaStream_t *stream){
  stream->tarma_size.inputWeight[0]=1;
  stream->tarma_size.maWeight[0]=0.6;
  stream->tarma_size.maWeight[1]=-1.04;
  stream->tarma_size.maWeight[2]=0.44;
  stream->tarma_size.arWeight[0]=1;
  stream->tarma_size.arWeight[1]=-1.971;
  stream->tarma_size.arWeight[2]=0.971;
  stream->tarma_size.polyWeight[0]=157;
  stream->tarma_size.polyWeight[1]=67.2;
  stream->tarma_size.polyWeight[2]=-11.2;
  stream->tarma_size.polyWeight[3]=-1.9;
  stream->tarma_size.polyWeight[4]=18.6;
  stream->tarma_size.polyWeight[5]=1.8;
  //tarmaPrintStreamInit(stream);
}

/* 
 this function prints the actual history of a random process
 input:		- random process
 output: 	
 */
void tarmaPrintProc(tarmaProcess_t *proc){
  char prefix[]="OTG TARMA DEBUG: ";
  int cntma, cntar;

  printf("%s tarmaPrintProc(%p) called\n", prefix, proc);
  printf("%s     ma history:\n",prefix);
  for(cntma=0; cntma<TARMA_NUM_MA_MAX; cntma++){
	printf("%s       ma[%d]: %f\n",prefix,cntma,proc->maHist[cntma]);
  }
  printf("%s     ar history:\n",prefix);
  for(cntar=0; cntar<TARMA_NUM_AR_MAX; cntar++){
	printf("%s       ar[%d]: %f\n",prefix,cntar,proc->arHist[cntar]);
  }
}

/* 
 this function prints the fixed (weigths) parameters of the input stream
 input:		- stream to be printed
 output: 	
 */
void tarmaPrintStreamInit(tarmaStream_t *stream){
  char prefix[]="OTG TARMA DEBUG: ";
  tarmaProcess_t *proc;
  int cntvar, cntp, cntma, cntar, cntpy;

  printf("%s tarmaPrintStreamInit(%p) called\n", prefix, stream);
  for(cntvar=0; cntvar<2; cntvar++){
	if(cntvar==0){
	  proc=&(stream->tarma_idt);
	  printf("%s     variable: idt\n",prefix);
	}else{
	  proc=&(stream->tarma_size);
	  printf("%s     variable: size\n",prefix);
	}
	printf("%s       input processes\n",prefix);
	for(cntp=0; cntp<TARMA_NUM_PROCESSES; cntp++){
	  printf("%s         w_pr[%d]: %f\n",prefix,cntp,proc->inputWeight[cntp]);
	}
	printf("%s       moving average weights\n",prefix);
	for(cntma=0; cntma<TARMA_NUM_MA_MAX; cntma++){
	  printf("%s         w_ma[%d]: %f\n",prefix,cntma,proc->maWeight[cntma]);
	}
	printf("%s       auto-regressive weights\n",prefix);
	for(cntar=0; cntar<TARMA_NUM_AR_MAX; cntar++){
	  printf("%s         w_ar[%d]: %f\n",prefix,cntar,proc->arWeight[cntar]);
	}
	printf("%s       polynomial factors\n",prefix);
	for(cntpy=0; cntpy<TARMA_NUM_POLY_MAX; cntpy++){
	  printf("%s         w_py[%d]: %f\n",prefix,cntpy,proc->polyWeight[cntpy]);
	}
  }
}


/*
 everything about video modeling which is an extended version of tarma modeling
 */

/* 
 this function calculates the packet size (frame size) for each video frame,
 	the idt is always fixed to the framerate (i.e., between 20 and 50 fps, default=25)
 input:		- vidoe stream
 output: 	- frame size
 */
double tarmaCalculateVideoSample(tarmaVideo_t *video){
  double size=0;
  double inputsamples[TARMA_NUM_PROCESSES];
  tarmaProcess_t *proc;
  int frameidx;
  int cntp, cntpy;

  if(video){
	proc=&(video->tarma_size);
	frameidx=video->tarmaVideoGopStructure[video->tarmaVideoFrameNumber];
  LOG_D(OTG,"TARMA_DEBUG: tarmaCalculateVideoSample(%p) called\n", video);
  LOG_D(OTG,"TARMA_DEBUG:     frameidx=%d\n",frameidx);
	if(frameidx>=0 && frameidx<=TARMA_NUM_FRAME_TYPES){
	  for(cntpy=0; cntpy<TARMA_NUM_POLY_MAX; cntpy++){
		proc->polyWeight[cntpy]=video->polyWeightFrame[frameidx][cntpy];
	  }
	  for(cntp=0; cntp<TARMA_NUM_PROCESSES; cntp++){
		inputsamples[cntp]=gaussian_dist(10000,1)-10000;
	  }
	  size=tarmaCalculateSample(inputsamples,proc);
	}
	/*get ready for the next frame*/
	video->tarmaVideoFrameNumber++;
	if(video->tarmaVideoFrameNumber>=TARMA_NUM_GOP_MAX){
	  video->tarmaVideoFrameNumber=0;
	}
	if(video->tarmaVideoGopStructure[video->tarmaVideoFrameNumber]<0 || 
	   video->tarmaVideoGopStructure[video->tarmaVideoFrameNumber]>TARMA_NUM_FRAME_TYPES){
		 video->tarmaVideoFrameNumber=0;
	   }
  }
  return size;
}

/* 
 this function initializes a video stream (and allocates memory if needed)
 input:		- video stream (or zero if memory should be allocated)
 output: 	- the same video stream
 */
tarmaVideo_t *tarmaInitVideo(tarmaVideo_t *video){
  tarmaProcess_t *proc;
  int cntp, cntma, cntar, cntpy, cntgop, cnttype;

  LOG_D(OTG,"TARMA_DEBUG: tarmaInitVideo(%p) called\n", video);
  if(video==0){
	video=(tarmaVideo_t*) malloc(sizeof(tarmaVideo_t));
  }
  proc=&(video->tarma_size);
  for(cntp=0; cntp<TARMA_NUM_PROCESSES; cntp++){
	proc->inputWeight[cntp]=0;
  }
  for(cntma=0; cntma<TARMA_NUM_MA_MAX; cntma++){
	proc->maWeight[cntma]=0;
	proc->maHist[cntma]=0;
  }
  for(cntar=0; cntar<TARMA_NUM_AR_MAX; cntar++){
	proc->arWeight[cntar]=0;
	proc->arHist[cntar]=0;
  }
  for(cntpy=0; cntpy<TARMA_NUM_POLY_MAX; cntpy++){
	proc->polyWeight[cntpy]=0;
  }
  video->tarmaVideoFrameNumber=0;
  for(cntgop=0; cntgop<TARMA_NUM_GOP_MAX; cntgop++){
	video->tarmaVideoGopStructure[cntgop]=-1;
  }
  for(cnttype=0; cnttype<TARMA_NUM_FRAME_TYPES; cnttype++){
	for(cntpy=0; cntpy<TARMA_NUM_POLY_MAX; cntpy++){
	  video->polyWeightFrame[cnttype][cntpy]=0;
	}
  }
  return video;
}

/* 
 this function initializes a video stream according to a tarma model
 input:		- video to be initialized
 			- compression rate: 1...dvd quality
 output: 	
 */
void tarmaSetupVideoGop12(tarmaVideo_t *video, double compression){
  if(video){
	video->tarma_size.inputWeight[0]=1;
	video->tarma_size.maWeight[0]=0.47;
	video->tarma_size.maWeight[1]=-0.829;
	video->tarma_size.maWeight[2]=0.358;
	video->tarma_size.arWeight[0]=1;
	video->tarma_size.arWeight[1]=-1.984;
	video->tarma_size.arWeight[2]=0.984;
	video->tarmaVideoGopStructure[0]=0; /*i frame*/
	video->tarmaVideoGopStructure[1]=2; /*b frame*/
	video->tarmaVideoGopStructure[2]=2;
	video->tarmaVideoGopStructure[3]=1; /*p frame*/
	video->tarmaVideoGopStructure[4]=2;
	video->tarmaVideoGopStructure[5]=2;
	video->tarmaVideoGopStructure[6]=1;
	video->tarmaVideoGopStructure[7]=2;
	video->tarmaVideoGopStructure[8]=2;
	video->tarmaVideoGopStructure[9]=1;
	video->tarmaVideoGopStructure[10]=2;
	video->tarmaVideoGopStructure[11]=2;
	if(compression<1){compression=1;}
	video->polyWeightFrame[0][0]=55400/compression; /*i frame*/
	video->polyWeightFrame[0][1]=32300/compression;
	video->polyWeightFrame[0][2]=10400/compression;
	video->polyWeightFrame[0][3]=-100/compression;
	video->polyWeightFrame[0][4]=-800/compression;
	video->polyWeightFrame[0][5]=74/compression;
	video->polyWeightFrame[1][0]=20900/compression; /*p frame*/
	video->polyWeightFrame[1][1]=19100/compression;
	video->polyWeightFrame[1][2]=7150/compression;
	video->polyWeightFrame[1][3]=-740/compression;
	video->polyWeightFrame[1][4]=-130/compression;
	video->polyWeightFrame[1][5]=20/compression;
	video->polyWeightFrame[2][0]=11700/compression; /*b frame*/
	video->polyWeightFrame[2][1]=10300/compression;
	video->polyWeightFrame[2][2]=4320/compression;
	video->polyWeightFrame[2][3]=700/compression;
	video->polyWeightFrame[2][4]=-90/compression;
	video->polyWeightFrame[2][5]=-2/compression;
	tarmaPrintVideoInit(video);
  }
}

/* 
 this function prints the fixed (weigths) parameters of the input video stream
 input:		- video stream to be printed
 output: 	
 */
void tarmaPrintVideoInit(tarmaVideo_t *video){
  char prefix[]="OTG TARMA DEBUG: ";
  tarmaProcess_t *proc;
  int cntp, cntma, cntar, cntpy, cntgop, cnttype;

  printf("%s tarmaPrintVideoInit(%p) called\n", prefix, video);
  proc=&(video->tarma_size);
  printf("%s      input process weights\n",prefix);
  for(cntp=0; cntp<TARMA_NUM_PROCESSES; cntp++){
	 printf("%s        w[%d]=%f\n",prefix,cntp,proc->inputWeight[cntp]);
  }
  printf("%s      ma weights\n",prefix);
  for(cntma=0; cntma<TARMA_NUM_MA_MAX; cntma++){
	printf("%s        ma[%d]=%f\n",prefix,cntma,proc->maWeight[cntma]);
  }
  printf("%s      ar weights\n",prefix);
  for(cntar=0; cntar<TARMA_NUM_AR_MAX; cntar++){
	printf("%s        ar[%d]=%f\n",prefix,cntar,proc->arWeight[cntar]);
  }
  printf("%s      polynomial weights\n",prefix);
  for(cntpy=0; cntpy<TARMA_NUM_POLY_MAX; cntpy++){
	printf("%s        p[%d]=%f\n",prefix,cntpy,proc->polyWeight[cntpy]);
  }
  printf("%s      video frame number=%d\n",prefix,video->tarmaVideoFrameNumber);
  printf("%s      gop structure\n",prefix);
  for(cntgop=0; cntgop<TARMA_NUM_GOP_MAX; cntgop++){
	printf("%s        gop[%d]=%d\n",prefix,cntgop,video->tarmaVideoGopStructure[cntgop]);
  }
  printf("%s      frame poly weights\n",prefix); 
  for(cntpy=0; cntpy<TARMA_NUM_POLY_MAX; cntpy++){
	printf("%s      [%d] ",prefix,cntpy);
	for(cnttype=0; cnttype<TARMA_NUM_FRAME_TYPES; cnttype++){
	  printf("     %05.3f",video->polyWeightFrame[cnttype][cntpy]);
	}
	printf("\n");
  }
}





/* 
 everything about background traffic modeling according to 
 "M. Laner, P. Svoboda, S. Schwarz, M. Rupp, Users in Cells: a Data Traffic Analysis, WCNC'12, Paris, France, 2012".
 */

/* 
 this function creates a random realization of the mean data rate of a background-session 
 input:		
 output: 	- random realization of R_s (B/s)
 */
double backgroundRateRnd(){
  double rate;
  rate=pow(10,lognormal_dist(1.3525, 0.1954))/8; /*Byte/s*/;
  if(rate>BACKGROUND_RATE_MAX){ 
	rate=BACKGROUND_RATE_MAX;
  }else if(rate<0){
	rate=0;
  }
  return rate;
}

/* 
 this function creates a random realization of the session duration of a background-session
 input:		
 output: 	- random realization of D_s (ms)
 */
int backgroundSessionDurationRnd(){
  int duration;
  duration = ceil(pow(10,exponential_dist(1/0.3591))*1000); /*ms*/
  if(duration<0){
	duration=0;
  }else if(duration>1<<30){
	duration=1<<30;
  }
  return duration;
}

/* 
 this function initializes a background stream and allocates memory if needed
 input:		- the bg-stream (or zero if memory shall be allocated)
 			- the mean number of users (sessions) to be expected for this stream
 output: 	- the same bg-stream
 */
backgroundStream_t *backgroundStreamInit(backgroundStream_t *stream, double lambda_n){
  int cnts, numactivenow;
  if(stream==0){
	stream=(backgroundStream_t*) malloc(sizeof(backgroundStream_t));
  }

  stream->meanNumSessions=lambda_n;
  stream->lastUpdateTime=0;
  numactivenow=poisson_dist(lambda_n); /*how many sessions are active (i.e., started before simulation)*/
  for(cnts=0; cnts<BACKGROUND_NUM_ACTIVE_MAX; cnts++){
	if(cnts<numactivenow){
	  /*these sessions are already active for unknown time and have therefore a random phaseshift in ds*/
	  stream->activeSessions[cnts].meanSessionRate=backgroundRateRnd();
	  stream->activeSessions[cnts].endTime=ceil((double)backgroundSessionDurationRnd()*uniform_rng()); /*ms*/
	}else{
	  stream->activeSessions[cnts].meanSessionRate=0;
	  stream->activeSessions[cnts].endTime=-1;
	}
  }

  LOG_D(OTG,"BACKGROUND_USERS DEBUG: backgroundStreamInit(%p) called\n", stream);
  backgroundPrintStream (stream);
  return stream;
}

/* 
 this function updates a given background stream
 input:		- bg-stream
			- current time
 output: 	
 */
void backgroundUpdateStream(backgroundStream_t *stream, int ctime){
  int numNewSessions, cnts, period;

  LOG_D(OTG,"BACKGROUND DEBUG: backgroundUpdateStream(stream*=%p,ctime=%d) called\n", stream, ctime);
  if(stream){
	period=ctime-stream->lastUpdateTime;
	numNewSessions=poisson_dist(stream->meanNumSessions/5710*period);
	for(cnts=0; cnts<BACKGROUND_NUM_ACTIVE_MAX; cnts++){
	  if(stream->activeSessions[cnts].endTime<ctime){
		if(numNewSessions>0){
		  stream->activeSessions[cnts].meanSessionRate=backgroundRateRnd();
		  stream->activeSessions[cnts].endTime=
			ctime-period*uniform_rng()+backgroundSessionDurationRnd(); /*ms*/
		  numNewSessions--;
		}else{
		  stream->activeSessions[cnts].meanSessionRate=0;
		  stream->activeSessions[cnts].endTime=-1;
		}
	  }
	}
	stream->lastUpdateTime=ctime;
	backgroundPrintStream (stream);
  }
}

/* 
 this function generates the packet size for a given bg-stream and a given idt,
 	the distribution of the idt can thereby be arbitrary (const idt=20ms recommended)
 input:		- bg-stream
	   		- current time
	   		- idt of the current packet
 output: 	- packet size of the current packet
 */
double backgroundCalculateSize(backgroundStream_t *stream, int ctime, int idt){
  int cnts, cntact=0;
  double size=0;
  double mrate=0;

  backgroundUpdateStream(stream, ctime);
  LOG_D(OTG,"BACKGROUND DEBUG: backgroundCalculateSize(stream*=%p,idt=%d,ctime=%d) called\n", stream, idt, ctime);
  if(stream){
	for(cnts=0; cnts<BACKGROUND_NUM_ACTIVE_MAX; cnts++){
	  if(stream->activeSessions[cnts].endTime>ctime){
		mrate+=stream->activeSessions[cnts].meanSessionRate;
		cntact++;
	  }
	}
	size=mrate*idt/1000;
	LOG_D(OTG,"BACKGROUND DEBUG:     cntact=%02d, idt=%05d, agg_mrate=%05.1f, size=%04.1f\n", cntact, idt, mrate,size);
  }
  return size;
}

/* 
 this function prints the current status of the bg-stream
 input:		- bg-stream
 output: 	
 */
void backgroundPrintStream(backgroundStream_t *stream){
  int cnts;
  
  LOG_D(OTG,"BACKGROUND DEBUG: backgroundPrintStream(%p)\n", stream);
  if(stream){
    LOG_D(OTG,"BACKGROUND DEBUG:     meanNumSessions(lambda_n)=%f\n",stream->meanNumSessions);
	for(cnts=0; cnts<BACKGROUND_NUM_ACTIVE_MAX; cnts++){
	  LOG_D(OTG,"BACKGROUND DEBUG:       session[%d] -> mrate=%06.3f, etime=%05d\n",
	        cnts, stream->activeSessions[cnts].meanSessionRate, stream->activeSessions[cnts].endTime);
	}
  }
}