#include <iostream>
#include <sstream>
#include <fstream>

#include "TStopwatch.h"

#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TCanvas.h"
#include "TString.h"
#include "TStyle.h"
#include "TChain.h"
#include "TH2.h"
#include "TH1.h"
#include "TF1.h"
#include "TTree.h"
#include "TKey.h"
#include "TMath.h"
#include "TRandom.h"
#include "TRandom3.h"
#include "Riostream.h"
#include "TVirtualFitter.h"
#include <stdio.h>
#include <math.h>

#include "Math/GoFTest.h"


void callcenter (const int Nclerks=10, double Tarrival=25, double MaxWait=1200, double Tave=5, double Tave_sigma=1, int maxtime=50000, int Ntoys=100) {

  // good combo: 5, 35, 1200, 5, 0.5, 50000, 100
  // another one: 8, 60, 1200, 5.5, 1, 100000, 100

  const double Nc = Nclerks; 
  //double maxtime = 1000000; // seconds of simulation

  int * ClientsPerClerk = new int[Nc];
  int * BusySince = new int[Nc];
  int * BusyUntil = new int[Nc];
  bool * Free = new bool[Nc];
  for (int ic=0; ic<Nclerks; ic++) {
    Free[ic]=true;
    BusyUntil[ic]=0;
    BusySince[ic]=0;
  }
  int NClientsInLine=0;
  int * ArrivalTime = new int[100000];
  int * ServedTime = new int[100000];
  int * est_time = new int[100000];
  int * est_time_unc = new int[100000];
  int TotClients=0;
  int ServedClients=0;
  int DroppedClients=0;
  int NowFree=Nclerks;
  //int Ntoys=1000;

  TH1D * WaitingTime = new TH1D ("WaitingTime", "Waiting time", 100, 0., MaxWait);
  TH1D * DeltaWaiting = new TH1D ("DeltaWaiting", "Delta waiting time", 100, -MaxWait, MaxWait);
  TH1D * CallTime = new TH1D ("CallTime", "Call duration", 100, 0., 1000.);
  TH1D * Est_time = new TH1D ("Est_time", "Estimated waiting time", 100, 0., 2000.);
  TH1D * Est_time_unc = new TH1D ("Est_time_unc", "Estimated uncertainty on waiting time", 100, 0., 400.);
  TH2D * TimeVsUnc = new TH2D ("TimeVsUnc", "Estimated time vs uncertainty", 100, 0., 1200., 100, 0., 400.);
  gStyle->SetOptStat(111111);

  // Time cycle
  for (int time=0; time<maxtime; time++) {
   
    if (time/10000 == 0.0001*time) {
      cout << endl;
      cout << "Time = " << time << "; Clients in line = " << NClientsInLine << "; Served clients = " << ServedClients << endl;
    }
    if (TotClients==100000) continue; // break cycle when we reach max

    // arrives a client
    if (gRandom->Uniform()<1/Tarrival) {
      NClientsInLine++;
      ArrivalTime[TotClients]=time;
      //Served[TotClients]=false;

      // The system estimates the waiting time for this client
      est_time[TotClients]=0;
      est_time_unc[TotClients]=0;
      if (NowFree==0) {
	double ave_est_time = 0;
	double ave_est_time2= 0;
	for (int toy=0; toy<Ntoys; toy++) {
	  int est_time_min=100000;
	  for (int ic=0; ic<Nclerks; ic++) {
	    int DurationThisCall = (int)exp(gRandom->Gaus(Tave,Tave_sigma));
	    do {
	      DurationThisCall=(int)exp(gRandom->Gaus(Tave,Tave_sigma));
	    } while (DurationThisCall<time-BusySince[ic]);
	    int est_time_this = DurationThisCall-(time-BusySince[ic]);
	    // add the other clients in line he must serve
	    for (int iclientsperclerk=0; iclientsperclerk<(NClientsInLine-1)/Nclerks; iclientsperclerk++) {
	      est_time_this+=(int)exp(gRandom->Gaus(Tave,Tave_sigma));
	    }
	    if (est_time_this<est_time_min) est_time_min=est_time_this;
	  }	
	  ave_est_time = ave_est_time + (est_time_min)/(double)Ntoys;
	  ave_est_time2= ave_est_time2+ pow(est_time_min,2)/(double)Ntoys;
	} // end toy
	est_time[TotClients]= ave_est_time;
	est_time_unc[TotClients]= sqrt(ave_est_time2-ave_est_time*ave_est_time);
	Est_time->Fill(ave_est_time);
	Est_time_unc->Fill(est_time_unc[TotClients]);
	TimeVsUnc->Fill(ave_est_time,est_time_unc[TotClients]);
      }

      TotClients++;
    } // if a new client is arriving

    // tired clients drop off 
    for (int i=ServedClients; i<TotClients; i++) {
      if (time-ArrivalTime[i]>MaxWait) {
	NClientsInLine--;
	ServedClients++;
	DroppedClients++;
      }
    }

    // management of clerks
    NowFree=0;
    for (int ic=0; ic<Nclerks; ic++) {
      // When a clerk is done we free her
      if (time>BusyUntil[ic]) {
	Free[ic]=true;
      }
      // free clerks take new client
      if (Free[ic]) {
	if (NClientsInLine>0) {
	  Free[ic]=false;
	  NClientsInLine--;
	  ServedTime[ServedClients]=time;
	  WaitingTime->Fill((double)(ServedTime[ServedClients]-ArrivalTime[ServedClients]));
	  if (est_time[ServedClients]>0) 
	    DeltaWaiting->Fill((double)(ServedTime[ServedClients]-ArrivalTime[ServedClients]-est_time[ServedClients]));
	  ServedClients++;
	  // when will this clerk be free again ?
	  double in_line = exp(gRandom->Gaus(Tave,Tave_sigma));
	  BusySince[ic]=time;
	  BusyUntil[ic]=time+(int)in_line;
	  CallTime->Fill(in_line);
	}
      }
      // see how many free clerks there are
      if (Free[ic]) NowFree++;
    } // end ic cycle

  }

  TCanvas * C = new TCanvas ("C", "Waiting time", 500, 500);
  C->Divide(2,3);
  C->cd(1);
  C->GetPad(1)->SetLogy();
  WaitingTime->Draw();
  C->cd(2);
  CallTime->Draw();
  C->cd(3);
  Est_time->Draw();
  C->cd(4);
  Est_time_unc->Draw();
  C->cd(5);
  DeltaWaiting->Draw();
  C->cd(6);
  TimeVsUnc->Draw();
  // End of program
  // --------------
  cout << "Goodbye." << endl;

}