From 3b5d512b7fcccbba52732d8b20238f513ecb4e99 Mon Sep 17 00:00:00 2001
From: Joao Carreira <jfmcarreira@gmail.com>
Date: Tue, 30 Apr 2013 10:57:41 +0100
Subject: [PATCH] Initialization of the GIT repository for the sqngen app

---
 ErrorPattern_Generator_Script |  33 ++
 Makefile                      |  29 ++
 sqnanalyzer.c                 | 147 +++++++
 sqngen.c                      | 783 ++++++++++++++++++++++++++++++++++
 4 files changed, 992 insertions(+)
 create mode 100755 ErrorPattern_Generator_Script
 create mode 100644 Makefile
 create mode 100644 sqnanalyzer.c
 create mode 100644 sqngen.c

diff --git a/ErrorPattern_Generator_Script b/ErrorPattern_Generator_Script
new file mode 100755
index 0000000..bb7c4ba
--- /dev/null
+++ b/ErrorPattern_Generator_Script
@@ -0,0 +1,33 @@
+#!/bin/bash
+
+#OUTPUT_DIR=ErrorPattern_Ballet_QP20_PCKT3778_N5_P1.1.7_B1.1.5
+
+SEQNAME=BookArrivalHEVC
+NLOOP=15
+NPACKETS=100
+PLOSS_INI=50
+PLOSS_INC=10
+PLOSS_FIN=50
+BLENGTH_INI=1
+BLENGTH_INC=1
+BLENGTH_FIN=1
+
+OUTPUT_DIR="ErrorPattern_"$SEQNAME"_PCKT"$NPACKETS_A"_N"$NLOOP"_P"$PLOSS_INI"-"$PLOSS_INC"-"$PLOSS_FIN"_B"$BLENGTH_INI"-"$BLENGTH_INC"-"$BLENGTH_FIN
+
+if [ $PLOSS_INC = 0 ]; then
+  PLOSS_INC=1
+fi
+
+if [ $PLOSS_INI = 0 ]; then
+  PLOSS_INI=$PLOSS_INI+$PLOSS_INC
+fi
+
+mkdir $OUTPUT_DIR
+
+for (( PLOSS=$PLOSS_INI; PLOSS<=$PLOSS_FIN; PLOSS=PLOSS+$PLOSS_INC ))  
+do
+	for (( BLENGTH=$BLENGTH_INI; BLENGTH<=$BLENGTH_FIN; BLENGTH=BLENGTH+$BLENGTH_INC ))  
+	do 
+		./SQNGen $NLOOP $NPACKETS 0 loss_GI $PLOSS $BLENGTH -o $OUTPUT_DIR"/"$SEQNAME
+	done
+done
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..698a6f6
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,29 @@
+###
+###     Makefile for GilberModel
+###
+###          
+###
+
+
+OUTPUT_DIR := Bin
+PROJECTS := sqnanalyzer sqngen
+
+export DBG
+export STC 
+export M32
+export OPENMP
+
+.HONY: default all $(PROJECTS) #distclean clean tags dependP
+
+default: all
+
+all: $(PROJECTS)
+
+sqngen:
+	g++ -Wall sqngen.c -lgsl -lgslcblas -o SQNGen
+sqnanalyzer:
+	g++ -Wall sqnanalyzer.c -o SQNAnalyzer
+
+clean:
+	rm SQNGen
+	rm SQNAnalyzer
\ No newline at end of file
diff --git a/sqnanalyzer.c b/sqnanalyzer.c
new file mode 100644
index 0000000..961a723
--- /dev/null
+++ b/sqnanalyzer.c
@@ -0,0 +1,147 @@
+// to compile g++ -Wall sqnanalyzer.c -o sqnanalyzer
+//
+//syntax:
+//	sqngenanalyzer
+//
+//This program is free software; you can redistribute it and/or
+//modify it under the terms of the GNU General Public License
+//s published by the Free Software Foundation; either version
+//2 of the License, or (at your option) any later version.
+//
+ 
+//Author:	Fabio Ludovici <fabio.ludovici at yahoo.it>
+//Version 1.0 final
+//July 2009
+// 	
+//This program was used to test the loss generation algorithm implemented in Netem2
+//available at http://netgroup.uniroma2.it/twiki/bin/view.cgi/Main/NetEm2
+// 
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdbool.h>
+#include <math.h>
+
+#define UNO 4294967295
+
+#include <iostream>
+#include <fstream>
+using namespace std;
+
+void blind_external_check(int length, int gmin_max, int sequence[], int sample, float *ploss, float *mean_length_of_burst, float *density, float *pisolated_loss, float *mean_length_of_good_burst) {
+	bool this_is_burst;
+	int num_of_bursts;
+	int length_of_all_bursts;
+	int num_of_good_bursts;
+	int length_of_all_good_bursts;
+	int num_of_drops;
+	int num_of_burst_drops;
+	int num_of_isolated_drops;
+	for(int gmin=2;gmin<gmin_max*2;gmin=gmin*2) {
+		this_is_burst=false;
+		num_of_bursts=0;
+		num_of_good_bursts=0;
+		length_of_all_bursts=0;
+		length_of_all_good_bursts=0;
+		num_of_drops=0;
+		num_of_burst_drops=0;
+		num_of_isolated_drops=0;
+		for (int k=0;k<length;k++) {
+			if (this_is_burst==true) {	//BURST ON this_is_burst=true
+				if ((sequence[k-1]==1) & (sequence[k]==0)) 
+					num_of_good_bursts++;
+				if (sequence[k]==0) {
+					length_of_all_bursts++;
+					length_of_all_good_bursts++;
+				} else if (sequence[k]==1) {
+					num_of_drops++;
+					num_of_burst_drops++;
+					length_of_all_bursts++;		
+				}
+				if ((k<(length-gmin)) & (sequence[k]==1) & (sequence[k+1]==0)) {
+					int zeros=1;
+					for(int l=1;l<gmin;l++) {
+						if (sequence[k+1+l]==0) zeros++;
+					}
+					if (zeros==gmin) 
+						this_is_burst=false;
+				}
+			} else {			//BURST OFF this_is_burst=false
+				if (k==(length-1)) {
+					if(sequence[k]==1) {
+						num_of_drops++;
+						num_of_isolated_drops++;
+					}
+				} else if ((sequence[k]==1) & (sequence[k+1]==1)) {
+					this_is_burst=true;
+					num_of_drops++;
+					num_of_bursts++;
+					num_of_burst_drops++;
+					length_of_all_bursts++;
+				} else if ((sequence[k]==1) & (sequence[k+1]==0)) {
+					num_of_drops++;
+					num_of_isolated_drops++;
+				}	
+			}
+		}
+		printf("\nFor gmin=%d the generated sequence has the following characteristics (from sequence analysis):\n",gmin);
+		int i=log2(gmin/2)+1;
+		int j=log2(gmin_max/2)+1;
+		ploss[sample*j+i-1]=100*((float)num_of_drops)/((float)length);
+		pisolated_loss[sample*j+i-1]=100*((float)num_of_isolated_drops)/((float)length);
+		if (num_of_bursts>0) 
+			mean_length_of_burst[sample*j+i-1]=(float) length_of_all_bursts/num_of_bursts;
+		else 
+			mean_length_of_burst[sample*j+i-1]=0;
+		if (num_of_good_bursts>0) 
+			mean_length_of_good_burst[sample*j+i-1]=(float) length_of_all_good_bursts/num_of_good_bursts;
+		else 
+			mean_length_of_good_burst[sample*j+i-1]=0;
+		density[sample*j+i-1]=100*(((float)num_of_burst_drops)/((float)length_of_all_bursts));
+		printf("P_loss %.3f%% Mean_length_of_burst %.3f Density %.3f%% P_isol %.3f%% Mean_length_of_good_burst %.3f\n", ploss[sample*j+i-1], mean_length_of_burst[sample*j+i-1], density[sample*j+i-1], pisolated_loss[sample*j+i-1], mean_length_of_good_burst[sample*j+i-1]);
+		fstream data;
+		data.open("data.dat", ios::out | ios::app);
+		data<<gmin<<";"<<fixed<<ploss[sample*j+i-1]<<";"<<mean_length_of_burst[sample*j+i-1]<<";"<<density[sample*j+i-1]<<";"<<pisolated_loss[sample*j+i-1]<<";"<<mean_length_of_good_burst[sample*j+i-1]<<";";
+		data.close();				 
+	}
+}
+
+int main (int argc, char *argv[]) {
+	if (argc<2) {
+		printf("Error: too few arguments\n");
+		printf("Syntax is: sqnanalyzer filename.dat\n");
+		return 0;
+	}
+	int num_of_samples=1;
+	int gmin_max=128;
+	char * buffer;
+	int length;
+	ifstream file (argv[1], ios::in|ios::binary|ios::ate);
+	length=file.tellg();
+	file.seekg (0, ios::beg);
+	buffer = new char [length];
+	file.read (buffer, length);
+	file.close();
+	int sequence[length];
+	printf("\nThe sequence is: ");
+	for(register int i=0; i<(length); i++) {
+		if(buffer[i]==';') 
+			break;
+		else 
+			sequence[i]=buffer[i];
+	if (sequence[i]==48) 
+		sequence[i]=0;
+	else if (sequence[i]==49) 
+		sequence[i]=1;
+	printf("%d", sequence[i]); 
+	}
+	float ploss[num_of_samples*gmin_max];
+	float pisolated_loss[num_of_samples*gmin_max];
+	float mean_length_of_burst[num_of_samples*gmin_max];
+	float density[num_of_samples*gmin_max];
+	float mean_length_of_good_burst[num_of_samples*gmin_max];
+	int sample=1; 
+	printf("\nLength is %d: ", length);
+	blind_external_check (length, gmin_max, sequence, sample, ploss, mean_length_of_burst, density, pisolated_loss, mean_length_of_good_burst);
+}
diff --git a/sqngen.c b/sqngen.c
new file mode 100644
index 0000000..7767b2f
--- /dev/null
+++ b/sqngen.c
@@ -0,0 +1,783 @@
+//to compile use the command g++ -Wall sqngen.c -lgsl -lgslcblas -o sqngen
+//syntax:
+//	sqngen num_of_samples length gmin_max loss_4state p13 p31 [p32 p23 [p14]]
+//OR	sqngen num_of_samples length gmin_max loss_GI [ploss [burst_length [density [pisol [good_burst_length]]]]
+//OR 	sqngen num_of_samples length gmin_max loss_bern_4s p
+//OR 	sqngen num_of_samples length gmin_max loss_gilb_4s p r [1-h]
+//OR 	sqngen num_of_samples length gmin_max loss_gilb_ell_4s p r [1-h [1-k]]
+//OR 	sqngen num_of_samples length gmin_max loss_gilb p r [1-h]
+//OR 	sqngen num_of_samples length gmin_max loss_gilb_ell p r [1-h [1-k]]
+//p13 p31 p32 p23 p14 ploss density pisol p r 1-k 1-h are percentage
+//num_of_samples length gmin_max burst_length good_burst_length are integer
+//
+
+//This program is free software; you can redistribute it and/or
+//modify it under the terms of the GNU General Public License
+//s published by the Free Software Foundation; either version
+//2 of the License, or (at your option) any later version.
+//
+
+//Author:	Fabio Ludovici <fabio.ludovici at yahoo.it>
+//Version 2.4
+//January 2010
+// 	
+//This program was used to test the loss generation algorithms implemented in Netem2
+//available at http://netgroup.uniroma2.it/twiki/bin/view.cgi/Main/NetEm2
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+#include <stdbool.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+#define ONE 4294967295u
+
+#include <iostream>
+using std::fixed;
+#include <fstream>
+using namespace std;
+
+double p31;
+double p13;
+double p23=1;
+double p32=0;
+double p14=0;
+
+double p;
+double r;
+double h;
+double k;
+
+double burst_length;
+double good_burst_length;
+double rho;
+double rho_min;
+double pisol;
+double ploss;
+
+void arg_error() {
+	printf("Error: too few arguments\n");
+	printf("Syntax is: sqngen num_of_samples length gmin_max loss_4state p13 p31 [p32 p23 p14]\n");
+	printf("or 	   sqngen num_of_samples length gmin_max loss_GI ploss [burst_length [density [pisol [good_burst_length]]]]\n");
+	printf("or 	   sqngen num_of_samples length gmin_max loss_bern_4s p\n");
+	printf("or 	   sqngen num_of_samples length gmin_max loss_gilb_4s p r [1-h]\n");
+	printf("or 	   sqngen num_of_samples length gmin_max loss_gilb_ell_4s p r [1-h [1-k]]\n");
+	printf("or 	   sqngen num_of_samples length gmin_max loss_gilb p r [1-h]\n");
+	printf("or 	   sqngen num_of_samples length gmin_max loss_gilb_ell p r [1-h [1-k]]\n");
+	printf("p13 p31 p32 p23 p14 ploss density pisol p r 1-k 1-h are percentage\n");
+	printf("num_of_samples length gmin_man burst_length good_burst_length are integer\n");
+}
+
+void sqn_gen(gsl_rng * q, FILE *out_file, char* sequence, int length, double p13,  double p31, double p32, double p23, double p14) {
+	double GI_random;
+	int state=1;
+	int tags=0;
+	//if(length>0) printf("The sequence is:\n");
+	gsl_rng_set(q, time(NULL));
+	for (int i=0;i<(length+tags);i++) {
+		GI_random=gsl_rng_get(q);
+		GI_random=GI_random/ONE;
+		switch (state) {
+
+		case 1:
+			sequence[i]=0;
+			//printf("0");
+			fprintf(out_file, "1\n");
+			if((0<GI_random)&&(GI_random<p13)){
+				state=3;
+				sequence[i+1]='S';
+				//printf("S");
+				tags++;
+				i++;
+				break;
+			} else if((p13<GI_random)&&(GI_random<(1-p14))){
+				state=1;
+				break;
+			} else if(((1-p14)<GI_random) && (GI_random<1)){
+				state=4;
+				break;
+			}
+
+		case 2:
+			sequence[i]=0;
+			//printf("0");
+			fprintf(out_file, "1\n");
+			if((0<GI_random)&&(GI_random<p23)){
+				state=3;
+				break;
+			} else {
+				state=2;
+				break;
+			}
+		case 3:
+			sequence[i]=1;
+			//printf("1");
+			fprintf(out_file, "0\n");
+			if((0<GI_random)&&(GI_random<p32)){
+				state=2;
+				break;
+			} else if((p32<GI_random)&&(GI_random<(p31+p32))){
+				state=1;
+				sequence[i+1]='E';	
+				//printf("E");
+				i++;
+				tags++;
+				break;
+			} else if((((p31+p32))<GI_random) && (GI_random<1)){
+				state=3;
+				break;
+			}
+		case 4:
+			sequence[i]=1;
+			//printf("1");
+			fprintf(out_file, "0\n");
+			state=1;
+			break;
+		}
+	}
+}
+
+void gilb_ell_gen(gsl_rng * q, char* sequence, int length, double p,  double r, double k, double h) {
+	double gilb_ell_random;
+	int state=1;
+	int tags=0;
+	if(length>0) printf("The sequence is:\n");
+	gsl_rng_set(q, time(NULL));
+	for (int i=0;i<(length+tags);i++) {
+		switch (state) {
+		case 1:	
+			gilb_ell_random=gsl_rng_get(q);
+			gilb_ell_random=gilb_ell_random/ONE;
+			if(gilb_ell_random<(1-k)) {
+				sequence[i]=1;
+				printf("1");
+			} else {
+				sequence[i]=0;
+				printf("0");
+				}
+			gilb_ell_random=gsl_rng_get(q);
+			gilb_ell_random=gilb_ell_random/ONE;
+			if(gilb_ell_random<p) {
+				sequence[i+1]='S';
+				printf("S");
+				i++;
+				tags++;
+				state=2;
+			}
+			break;
+		case 2:	
+			gilb_ell_random=gsl_rng_get(q);
+			gilb_ell_random=gilb_ell_random/ONE;
+			if(gilb_ell_random<(1-h)) {
+				sequence[i]=1;
+				printf("1");
+			} else {
+				sequence[i]=0;
+				printf("0");
+			}
+			gilb_ell_random=gsl_rng_get(q);
+			gilb_ell_random=gilb_ell_random/ONE;
+			if(gilb_ell_random<r) {
+				sequence[i+1]='E';
+				printf("E");
+				i++;
+				tags++;	
+				state=1;
+			}
+			break;
+		}
+	};
+}
+
+void internal_check(int length, char sequence[], int sample, double *ploss_bis, double *mean_length_of_burst_bis, double *density_bis, double *pisolated_loss_bis, double *mean_length_of_good_burst_bis) {
+	bool gilb_ell_burst=false;
+	bool real_burst=false;
+	bool good_burst=false;
+	int num_of_bursts_bis=0;
+	int length_of_all_bursts_bis=0;
+	int num_of_good_bursts_bis=0;
+	int length_of_all_good_bursts_bis=0;
+	int num_of_drops_bis=0;
+	int num_of_burst_drops_bis=0;
+	int num_of_isolated_drops_bis=0;
+	int initial_zeros=0;
+	int final_zeros=0;
+	int ones=0;
+	int tags=0;
+	for (int k=0;k<(length+tags);k++) {
+	if (sequence[k]=='S') 	{
+		gilb_ell_burst=true;
+		num_of_bursts_bis++;
+		tags++;
+	}
+	if (sequence[k]=='E') {
+		gilb_ell_burst=false;
+		real_burst=false;
+		good_burst=false;
+		tags++;
+		if (ones==0) {
+			num_of_bursts_bis--;
+				}
+		length_of_all_bursts_bis=length_of_all_bursts_bis-initial_zeros-final_zeros;
+		length_of_all_good_bursts_bis=length_of_all_good_bursts_bis-initial_zeros-final_zeros;
+		ones=0;
+		initial_zeros=0;
+		final_zeros=0;
+	}
+	if ((gilb_ell_burst==true) & (ones==0) & (sequence[k]==1)) {
+		real_burst=true;
+	}
+	if ((gilb_ell_burst==true) & ((sequence[k]==1) || (sequence[k]==0))) {
+		length_of_all_bursts_bis++;
+		if (sequence[k]==0) {
+			if (real_burst==false)
+				initial_zeros++;
+			else
+				final_zeros++;
+			if (real_burst==true) {
+				length_of_all_good_bursts_bis++;
+				if (good_burst==false) {
+					good_burst=true;
+					num_of_good_bursts_bis++;
+				}
+			}
+		}
+	if (sequence[k]==1) {
+		if (real_burst==false) real_burst=true;
+		if (good_burst==true) good_burst=false;
+		ones++;
+		final_zeros=0;
+		num_of_drops_bis++;
+		num_of_burst_drops_bis++;
+	}
+	}
+	if ((gilb_ell_burst==false) & (sequence[k]==1)) {
+		num_of_drops_bis++;
+		num_of_isolated_drops_bis++;
+	}
+	}
+	printf("\nThe generated sequence has the following characteristics (from the internal check):\n");
+	ploss_bis[sample]=100*((double)num_of_drops_bis)/((double)length);
+	pisolated_loss_bis[sample]=100*((double)num_of_isolated_drops_bis)/((double)length);
+	if (num_of_bursts_bis>0) mean_length_of_burst_bis[sample]=(double) length_of_all_bursts_bis/num_of_bursts_bis;
+	else mean_length_of_burst_bis[sample]=0;
+	if (num_of_good_bursts_bis>0) mean_length_of_good_burst_bis[sample]=(double) length_of_all_good_bursts_bis/num_of_good_bursts_bis;
+	else mean_length_of_good_burst_bis[sample]=0;
+	density_bis[sample]=100*(double) num_of_burst_drops_bis/length_of_all_bursts_bis;
+	printf("P_loss %.3f%% Mean_length_of_burst %.3f Density %.3f%% P_isol %.3f%% Mean_length_of_good_burst %.3f\n", ploss_bis[sample], mean_length_of_burst_bis[sample], density_bis[sample], pisolated_loss_bis[sample], mean_length_of_good_burst_bis[sample]);
+	fstream data;
+	data.open("data.dat", ios::out | ios::app);
+	data<<fixed<<ploss_bis[sample]<<";"<<mean_length_of_burst_bis[sample]<<";"<<density_bis[sample]<<";"<<pisolated_loss_bis[sample]<<";"<<mean_length_of_good_burst_bis[sample]<<";";
+	data.close();
+}
+
+void external_blind_check(int length, int gmin_max, char sequence[], int sample, double *ploss_ter, double *mean_length_of_burst_ter, double *density_ter, double *pisolated_loss_ter, double *mean_length_of_good_burst_ter) {
+	bool this_is_burst=false;
+	int num_of_bursts_ter=0;
+	int length_of_all_bursts_ter=0;
+	int num_of_good_bursts_ter=0;
+	int length_of_all_good_bursts_ter=0;
+	int num_of_drops_ter=0;
+	int num_of_burst_drops_ter=0;
+	int num_of_isolated_drops_ter=0;
+	int tags=0;  
+	for (int k=0;k<(length+tags);k++) {
+	if((sequence[k]=='S') || (sequence[k]=='E')) {
+		sequence[k]=sequence[k+1];
+		tags++;
+	}
+	}
+	for(int gmin=2;gmin<gmin_max*2;gmin=gmin*2) {
+		this_is_burst=false;
+		num_of_bursts_ter=0;
+		num_of_good_bursts_ter=0;
+		length_of_all_bursts_ter=0;
+		length_of_all_good_bursts_ter=0;
+		num_of_drops_ter=0;
+		num_of_burst_drops_ter=0;
+		num_of_isolated_drops_ter=0;
+		for (int k=0;k<length;k++) {
+			if (this_is_burst==true) {	/* BURST ON this_is_burst=true */
+				if ((sequence[k-1]==1) & (sequence[k]==0)) num_of_good_bursts_ter++;
+				if (sequence[k]==0)	{
+					length_of_all_bursts_ter++;
+					length_of_all_good_bursts_ter++;
+				} else if (sequence[k]==1) {
+					num_of_drops_ter++;
+					num_of_burst_drops_ter++;
+					length_of_all_bursts_ter++;
+				}      
+				if ((k<(length-gmin)) & (sequence[k]==1) & (sequence[k+1]==0)) {
+				int zeros=1;
+				for(int l=1;l<gmin;l++)	{
+					if (sequence[k+1+l]==0) zeros++;
+				}
+				if (zeros==gmin) this_is_burst=false;
+				}
+				} else {	/* BURST OFF this_is_burst=false */
+					if (k==(length-1)) {
+						if(sequence[k]==1) {
+						num_of_drops_ter++;
+							num_of_isolated_drops_ter++;
+						}
+					} else if ((sequence[k]==1) & (sequence[k+1]==1)) {
+						this_is_burst=true;
+						num_of_drops_ter++;
+						num_of_bursts_ter++;
+						num_of_burst_drops_ter++;
+						length_of_all_bursts_ter++;
+					} else if ((sequence[k]==1) & (sequence[k+1]==0)) {
+						num_of_drops_ter++;
+						num_of_isolated_drops_ter++;
+					}
+				}
+		}
+	printf("\nFor gmin=%d the generated sequence has the following characteristics (from sequence analysis):\n",gmin);
+	int i=log2(gmin/2)+1;
+	int j=log2(gmin_max/2)+1;
+	ploss_ter[sample*j+i-1]=100*((double)num_of_drops_ter)/((double)length);
+	pisolated_loss_ter[sample*j+i-1]=100*((double)num_of_isolated_drops_ter)/((double)length);
+	if (num_of_bursts_ter>0) mean_length_of_burst_ter[sample*j+i-1]=(double) length_of_all_bursts_ter/num_of_bursts_ter;
+	else mean_length_of_burst_ter[sample*j+i-1]=0;
+	if (num_of_good_bursts_ter>0) mean_length_of_good_burst_ter[sample*j+i-1]=(double) length_of_all_good_bursts_ter/num_of_good_bursts_ter;
+	else mean_length_of_good_burst_ter[sample*j+i-1]=0;
+	density_ter[sample*j+i-1]=100*(((double)num_of_burst_drops_ter)/((double)length_of_all_bursts_ter));
+	printf("P_loss %.3f%% Mean_length_of_burst %.3f Density %.3f%% P_isol %.3f%% Mean_length_of_good_burst %.3f\n", ploss_ter[sample*j+i-1], mean_length_of_burst_ter[sample*j+i-1], density_ter[sample*j+i-1], pisolated_loss_ter[sample*j+i-1], mean_length_of_good_burst_ter[sample*j+i-1]);
+	fstream data;
+	data.open("data.dat", ios::out | ios::app);
+	data<<gmin<<";"<<fixed<<ploss_ter[sample*j+i-1]<<";"<<mean_length_of_burst_ter[sample*j+i-1]<<";"<<density_ter[sample*j+i-1]<<";"<<pisolated_loss_ter[sample*j+i-1]<<";"<<mean_length_of_good_burst_ter[sample*j+i-1]<<";";
+	data.close();
+	}
+}
+						
+void statistics(char *filename, int num_of_samples, double ploss[], double mean_length_of_burst[], double density[], double pisolated_loss[], double mean_length_of_good_burst[], int gmin, int gmin_max) {
+	double ploss_mean=0;
+	double mean_length_of_burst_mean=0;
+	double density_mean=0;
+	double pisolated_loss_mean=0;
+	double mean_length_of_good_burst_mean=0;
+	double ploss_variance=0;
+	double mean_length_of_burst_variance=0;
+	double density_variance=0;
+	double pisolated_loss_variance=0;
+	double mean_length_of_good_burst_variance=0;
+	double ploss_stddev=0;
+	double mean_length_of_burst_stddev=0;
+	double density_stddev=0;
+	double pisolated_loss_stddev=0;
+	double mean_length_of_good_burst_stddev=0;
+	double ploss_cx=0;
+	double mean_length_of_burst_cx=0;
+	double density_cx=0;
+	double pisolated_loss_cx=0;
+	double mean_length_of_good_burst_cx=0;
+	int i=log2(gmin/2)+1;
+	int j=log2(gmin_max/2)+1;
+	for(int sample=0;sample<num_of_samples;sample++) {
+		ploss_mean=ploss_mean+ploss[sample*j+i-1];
+		mean_length_of_burst_mean=mean_length_of_burst_mean+mean_length_of_burst[sample*j+i-1];
+		density_mean=density_mean+density[sample*j+i-1];
+		pisolated_loss_mean=pisolated_loss_mean+pisolated_loss[sample*j+i-1];
+		mean_length_of_good_burst_mean=mean_length_of_good_burst_mean+mean_length_of_good_burst[sample*j+i-1];
+	}
+	printf("\nThese are the mean values:\n");
+	ploss_mean=ploss_mean/num_of_samples;
+	mean_length_of_burst_mean=mean_length_of_burst_mean/num_of_samples;
+	density_mean=density_mean/num_of_samples;
+	pisolated_loss_mean=pisolated_loss_mean/num_of_samples;
+	mean_length_of_good_burst_mean=mean_length_of_good_burst_mean/num_of_samples;
+	printf("ploss %.3f%% mean_length_of_burst %.3f density %.3f%% pisol %.3f%% mean_length_of_good_burst %.3f\n", ploss_mean, mean_length_of_burst_mean, density_mean, pisolated_loss_mean, mean_length_of_good_burst_mean); 
+	for(int sample=0;sample<num_of_samples;sample++) { 
+		ploss_variance=ploss_variance+(ploss[sample*j+i-1]-ploss_mean)*(ploss[sample*j+i-1]-ploss_mean);
+		mean_length_of_burst_variance=mean_length_of_burst_variance+(mean_length_of_burst[sample*j+i-1]-mean_length_of_burst_mean)*(mean_length_of_burst[sample*j+i-1]-mean_length_of_burst_mean);
+		density_variance=density_variance+(density[sample*j+i-1]-density_mean)*(density[sample*j+i-1]-density_mean);
+		pisolated_loss_variance=pisolated_loss_variance+(pisolated_loss[sample*j+i-1]-pisolated_loss_mean)*(pisolated_loss[sample*j+i-1]-pisolated_loss_mean);
+		mean_length_of_good_burst_variance=mean_length_of_good_burst_variance+(mean_length_of_good_burst[sample*j+i-1]-mean_length_of_good_burst_mean)*(mean_length_of_good_burst[sample*j+i-1]-mean_length_of_good_burst_mean);
+	}
+	printf("\nThese are the variance values:\n");
+	ploss_variance=ploss_variance/(num_of_samples-1);
+	mean_length_of_burst_variance=mean_length_of_burst_variance/(num_of_samples-1);
+	density_variance=density_variance/(num_of_samples-1);
+	pisolated_loss_variance=pisolated_loss_variance/(num_of_samples-1);
+	mean_length_of_good_burst_variance=mean_length_of_good_burst_variance/(num_of_samples-1);  
+	printf("ploss %.3f mean_length_of_burst %.3f density %.3f pisol %.3f mean_length_of_good_burst %.3f\n", ploss_variance, mean_length_of_burst_variance, density_variance, pisolated_loss_variance, mean_length_of_good_burst_variance); 
+	printf("\nThese are the standard deviation values:\n"); 
+	ploss_stddev=sqrt(ploss_variance);
+	mean_length_of_burst_stddev=sqrt(mean_length_of_burst_variance);
+	density_stddev=sqrt(density_variance);
+	pisolated_loss_stddev=sqrt(pisolated_loss_variance);
+	mean_length_of_good_burst_stddev=sqrt(mean_length_of_good_burst_variance);
+	printf("ploss %.3f%% mean_length_of_burst %.3f density %.3f%% pisol %.3f%% mean_length_of_good_burst %.3f\n", ploss_stddev, mean_length_of_burst_stddev, density_stddev, pisolated_loss_stddev, mean_length_of_good_burst_stddev);
+	printf("\nThese are the variation coefficients:\n"); 
+	ploss_cx=ploss_stddev/ploss_mean;
+	mean_length_of_burst_cx=mean_length_of_burst_stddev/mean_length_of_burst_mean;
+	density_cx=density_stddev/density_mean;
+	if(pisolated_loss_mean>0) 
+		pisolated_loss_cx=pisolated_loss_stddev/pisolated_loss_mean;
+	else 
+		pisolated_loss_cx=0;
+	mean_length_of_good_burst_cx=mean_length_of_good_burst_stddev/mean_length_of_good_burst_mean;
+	printf("ploss %.3f mean_length_on_burst %.3f density %.3f pisol %.3f mean_length_of_good_burst %.3f\n", ploss_cx, mean_length_of_burst_cx, density_cx, pisolated_loss_cx, mean_length_of_good_burst_cx);   
+	fstream data;
+	data.open(filename, ios::out | ios::app);
+	data <<"Average Prob. Loss; Average Burst Length; Average Density; Average Prob. Loss Isol.; Average Good Burst Length;\n";	data<<fixed<<ploss_mean<<";"<<mean_length_of_burst_mean<<";"<<density_mean<<";"<<pisolated_loss_mean<<";"<<mean_length_of_good_burst_mean<<";\n";
+	
+	data <<"Variance Prob. Loss; Variance Burst Length; Variance Density; Variance Prob. Loss Isol.; Variance Good Burst length;\n";	data<<ploss_variance<<";"<<mean_length_of_burst_variance<<";"<<density_variance<<";"<<pisolated_loss_variance<<";"<<mean_length_of_good_burst_variance<<";\n";
+	
+	data <<"Standard Dev. Prob. Loss; Standard Dev. Burst Length; Standard Dev. Density; Standard Dev. Prob. Loss Isol.; Standard Dev. Good Burst length;\n";
+data<<ploss_stddev<<";"<<mean_length_of_burst_stddev<<";"<<density_stddev<<";"<<pisolated_loss_stddev<<";"<<mean_length_of_good_burst_stddev<<";\n";
+	data <<"Variation Prob. Loss; Variation Burst Length; Variation Density; Variation Prob. Loss Isol.; Variation Good Burst length;\n";
+	data<<ploss_cx<<";"<<mean_length_of_burst_cx<<";"<<density_cx<<";"<<pisolated_loss_cx<<";"<<mean_length_of_good_burst_cx<<";\n";
+	data.close(); 
+}
+
+int main (int argc, char *argv[]) {
+	int num_of_samples;
+	int length;
+	int gmin_max;
+	fstream data;
+	char *filename_common;
+	fprintf(stdout,"********************** SQNGEN (Edited by Joao Carreira) ***********************\n");
+	if (argc<7) {
+		arg_error();
+		return 0;
+	}
+	for (int i=0;i<argc;i++)
+	{
+		if (strcmp (argv[i],"-o") == 0)
+		{
+			filename_common = argv[i+1];
+			argc = argc - 2;
+			break;
+		}
+	}
+	if (argv[1]) num_of_samples=atoi(argv[1]);
+	if (argv[2]) length=atoi(argv[2]);
+	gmin_max=atoi(argv[3]);
+	if (strcmp (argv[4],"loss_4state") == 0) {
+		if (argc<7)	{
+			arg_error();
+			return 0;
+		}			
+		if (argc>7 && (atof(argv[6])+atof(argv[7])>100)) {printf("\nError: p31+p32>100\n"); return 0;}
+		if (argc>9 && (atof(argv[5])+atof(argv[9])>100)) {printf("\nError: p13+p14>100\n"); return 0;}
+		p13=atof(argv[5])/100;
+		p31=atof(argv[6])/100;
+		if (argc>7) 
+			p32=atof(argv[7])/100;
+		else 
+			p32=0;
+		if (argc>8) 
+			p23=atof(argv[8])/100;
+		else 
+			p23=1;
+		if (argc>9) 
+			p14=atof(argv[9])/100;
+		else 
+			p14=0;
+	} else if (strcmp (argv[4],"loss_GI") == 0) { 
+		if (argc<6) {
+			arg_error();
+			return 0;
+		}
+		/* Change here to get a 0.1% resolution */
+		//ploss=(double) atof(argv[5])/100; 
+		ploss=(double) atof(argv[5])/1000; 
+		if (ploss==0) {
+			printf("\nError: ploss=0\n"); 
+			return 0;
+		}
+		if (argc>8 && (atof(argv[8])>atof(argv[5]))) {
+			printf("\nError: pisol>ploss\n"); 
+			return 0;
+		}
+		if (argc>9 && (atof(argv[9])>atof(argv[6]))) {
+			printf("\nError: good_burst_length>burst_length\n"); 
+			return 0;
+		}
+		if (argc>6) 
+			burst_length=atof(argv[6]);
+		else 
+			burst_length=1/(1-ploss);
+		if (argc>7) 
+			rho=atof(argv[7])/100;
+		else 
+			rho=1;
+		if (argc>8) 
+			pisol=atof(argv[8])/100;
+		if (argc>9) 
+			good_burst_length=atoi(argv[9]);
+		else 
+			good_burst_length=0;
+		p13=(pisol-ploss)/(burst_length*(pisol*(rho-1)-rho+ploss));
+		p31=1/(burst_length*rho);
+		if (rho<1) {	
+		if (good_burst_length>0) {
+				p23=(1/good_burst_length);
+				p32=(1-rho)/(rho*good_burst_length);
+			}
+			else { 
+				p32=(1+rho*rho*burst_length-rho-burst_length*rho)/(rho-rho*burst_length);
+				p23=(burst_length*rho-1)/(burst_length-1);
+			}
+		} else {
+			p32=0;
+			p23=1;
+		}
+		if (pisol>0) 
+			p14=(double) pisol/(1-pisol);
+		else 
+			p14=0;
+	} else if (strcmp (argv[4],"loss_gilb_4s") == 0) {
+		if (argc<7) {
+			arg_error();
+			return 0;
+		}
+		p=(atof(argv[5]))/100;
+		r=(atof(argv[6]))/100;
+		if(argv[7]) 
+			h=1-((atof(argv[7]))/100);
+		else 
+			h=0;
+		ploss=(p-h*p)/(p+r);
+		burst_length=(1-h+r*r*h)/(r*(1-h+r*h));
+		if(argv[7]) 
+			rho=((1-h+r*h)*(1-h+r*h)) / (1-h+r*r*h);
+		else 
+			rho=1;	
+		p13=(-ploss)/(burst_length*(-rho+ploss));
+		p31=1/(burst_length*rho);
+		if (rho<1) {
+        		p32=(1+rho*rho*burst_length-rho-burst_length*rho)/(rho-rho*burst_length);
+			p23=(burst_length*rho-1)/(burst_length-1);
+		} else {
+			p32=0;
+			p23=1;
+		}
+		p14=0;
+	} else if (strcmp (argv[4],"loss_gilb_ell_4s") == 0) {
+		if (argc<7) {
+			arg_error();
+			return 0;
+		}
+		p=(atof(argv[5]))/100;
+		r=(atof(argv[6]))/100;
+		if(argv[7]) 
+			h=1-((atof(argv[7]))/100);
+		else 
+			h=0;
+		if(argv[8]) 
+			k=1-((atof(argv[8]))/100);
+		else 
+			k=1;
+		ploss=(r-k*r+p-h*p)/(p+r);
+		burst_length=(1-h+r*r*h)/(r*(1-h+r*h));
+		if(argv[7]) 
+			rho=((1-h+r*h)*(1-h+r*h)) / (1-h+r*r*h);
+		else 
+			rho=1;
+		if(argv[8]) 
+			pisol=(r-k*r)/(p+r);
+		else 
+			pisol=0;
+		p13=(pisol-ploss)/(burst_length*(pisol*(rho-1)-rho+ploss));
+		p31=1/(burst_length*rho);
+		if (argv[7]) {	
+        		p32=(1+rho*rho*burst_length-rho-burst_length*rho) / (rho-rho*burst_length);
+			p23=(burst_length*rho-1) / (burst_length-1);
+		} else 
+			p32=0;
+		if (pisol>0) 
+			p14=pisol/(1-pisol);
+		else 
+			p14=0;
+	} else if (strcmp (argv[4],"loss_bern") == 0) {
+		p=atof(argv[5])/100;
+		r=1-p;
+		h=0;
+		k=1;
+	} else if (strcmp (argv[4],"loss_gilb") == 0) {
+		if (argc<7) {
+			arg_error();
+			return 0;
+		}			
+		p=atof(argv[5])/100;
+		r=atof(argv[6])/100;
+		if(argc>7) {
+			h=1-((atof(argv[7]))/100);
+			k=0;
+		} else {
+	    		h=0;
+	    		k=1;
+	     	}
+	} else if (strcmp (argv[4],"loss_gilb_ell") == 0) {
+		if (argc<7) {
+			arg_error();
+			return 0;
+		}
+		p=atof(argv[5])/100;
+		r=atof(argv[6])/100;
+		if(argc>7){
+			h=1-((atof(argv[7]))/100);
+			k=0;
+		} else {	
+	    		h=0;
+	    		k=1;
+	      	}
+		if(argc>8) {
+			h=1-((atof(argv[7]))/100);
+			k=1-((atof(argv[8]))/100);
+		}
+		else {
+	    		k=1;
+	      	}	
+	}
+	char sequence[2*length];
+	double ploss_bis[num_of_samples];
+	double pisolated_loss_bis[num_of_samples];
+	double mean_length_of_burst_bis[num_of_samples];
+	double density_bis[num_of_samples];
+	double mean_length_of_good_burst_bis[num_of_samples];// good burst length
+	double ploss_ter[num_of_samples*gmin_max];
+	double pisolated_loss_ter[num_of_samples*gmin_max];
+	double mean_length_of_burst_ter[num_of_samples*gmin_max];
+	double density_ter[num_of_samples*gmin_max];
+	double mean_length_of_good_burst_ter[num_of_samples];	// good burst length
+	gsl_rng * q=gsl_rng_alloc (gsl_rng_taus2);
+	data.open("data.dat", ios::out);
+	data << "Sample; Num of Samples; Legth; Model; PLoss; BurstLegth";
+	data.close(); 
+	for (int sample=0;sample<num_of_samples;sample++) {  
+		FILE *out_file = NULL;
+		char filename[100];
+		sprintf(filename,"%s_P%d_B%d_%d_ERRORPATTERN.dat", filename_common, (int)(ploss*1000), (int)burst_length, sample);
+		out_file = fopen (filename, "w");
+		if (out_file == NULL)
+		{
+			printf("Error opening the error pattern file\n");
+			return 0;
+		}
+		fprintf(out_file, "nPackets: %d\nVector:\n",length);
+		data.open("data.dat", ios::out | ios::app);
+		data <<"\n"<<(sample+1)<<";";
+		data <<num_of_samples<<";";
+		data <<length<<";";
+		data <<gmin_max<<";";
+		data <<argv[4]<<";";		//chosen model
+		data <<atof(argv[5])<<";";	//saving input values
+		if(argv[6]) data <<atof(argv[6])<<";"; else data <<"N/A"<<";";
+		if((argv[6]) && (argv[7])) data <<atof(argv[7])<<";"; else data <<"100"<<";";
+		if((argv[6]) && (argv[7]) && (argv[8])) data <<atof(argv[8])<<";"; else data <<"0"<<";";
+		if((argv[6]) && (argv[7]) && (argv[8]) && (argv[9])) data <<atof(argv[9])<<";"; else data <<"N/A"<<";";
+		data.close(); 
+		//sleep (1); 
+		if ((strcmp (argv[4],"loss_bern") == 0) || (strcmp (argv[4],"loss_gilb") == 0) || (strcmp (argv[4],"loss_gilb_ell") == 0)) 
+			gilb_ell_gen(q, sequence, length, p, r, k, h);
+		else 
+			sqn_gen(q, out_file, sequence, length, p13, p31, p32, p23, p14);
+		fclose (out_file);
+		printf ("The selected generation option is %s with the following input parameters:\n", argv[4]);
+		if (strcmp (argv[4],"loss_4state") == 0) 
+			printf("p13=%.3f%% \np31=%.3f%% \np32=%.3f%% \np23=%.3f%% \np14=%.3f%%\n", (double)100*p13/1, (double)100*p31/1, (double)100*p32/1, (double)100*p23/1, (double)100*p14/1);
+		else if ((strcmp (argv[4],"loss_GI") == 0) & (good_burst_length==0))	
+			printf("P_loss %.3f%% Mean_length_of_burst %.3f%% Density %.3f%% P_isol %.3f%% Mean_length_of_good_burst N/A\n", 100*ploss, burst_length, 100*rho, 100*pisol);
+		else if ((strcmp (argv[4],"loss_GI") == 0) & (good_burst_length==0))	
+			printf("P_loss %.3f%% Mean_length_of_burst %.3f%% Density %.3f%% P_isol %.3f%% Mean_length_of_good_burst N/A\n", 100*ploss, burst_length, 100*rho, 100*pisol);
+		else if ((strcmp (argv[4],"loss_GI") == 0) & (good_burst_length!=0))	
+			printf("P_loss %.3f%% Mean_length_of_burst %.3f%% Density %.3f%% P_isol %.3f%% Mean_length_of_good_burst %.3f%%\n", 100*ploss, burst_length, 100*rho, 100*pisol, good_burst_length);
+		else if (strcmp (argv[4],"loss_gilb_4s") == 0) 
+			printf("p=%.3f%% r=%.3f%% 1-h=%.3f%%\n", 100*p, 100*r, 100*(1-h)/1);
+		else if (strcmp (argv[4],"loss_gilb_ell_4s") == 0) 
+			printf("p=%.3f%% r=%.3f%% 1-h=%.3f%% 1-k=%.3f%%\n", 100*p, 100*r, 100*(1-h), 100*(1-k)); 	    
+		else if (strcmp (argv[4],"loss_bern") == 0)	
+			printf("p=%.3f%%\n", (double)100*p/1);
+		else if (strcmp (argv[4],"loss_gilb") == 0)	
+			printf("p=%.3f%% r=%.3f%% 1-h=%.3f%%\n", (double)100*p/1, (double)100*r/1, (double) 100*(1-h)/1);
+		else if (strcmp (argv[4],"loss_gilb_ell") == 0)	
+			printf("p=%.3f%% r=%.3f%% 1-h=%.3f%% 1-k=%.3f%%\n", (double)100*p, (double)100*r, (double) 100*(1-h), (double) 100*(1-k)); 
+		if ((strcmp (argv[4],"loss_GI") == 0) || (strcmp (argv[4],"loss_bern_4s") == 0) || (strcmp (argv[4],"loss_gilb_4s") == 0) || (strcmp (argv[4],"loss_gilb_ell_4s") == 0)) {
+			printf ("The correspondent transition probabilities are:\n");
+			printf("p13=%.3f%% \np31=%.3f%% \np32=%.3f%% \np23=%.3f%% \np14=%.3f%%\n", (double)100*p13/1, (double)100*p31/1, (double)100*p32/1, (double)100*p23/1, (double)100*p14/1);
+		}
+		if ((strcmp (argv[4],"loss_GI") != 0) & (strcmp (argv[4],"loss_gilb") != 0) & (strcmp (argv[4],"loss_gilb_ell") != 0)) {
+			printf ("The correspondent GI (General and Intuitive) parameters are:");
+			printf ("P_loss=%.3f%% Mean_length_of_burst=%.3f Density=%.3f%% P_isol=%.3f%% Mean_length_of_good_burst %.3f \n", (double) 100*(p13*p23+p14*p23*p31) / (p13*p23+p23*p31+p14*p23*p31+p13*p32), (double) (p32+p23) / (p23*p31), (double) 100*(p13*p23) / (p13*p23+p13*p32), (double) 100*(p14*p23*p31) / (p14*p23*p31+p23*p31), (double) 1/p23 );
+		}	
+		if(length>0) {
+			internal_check(length, sequence, sample, ploss_bis, mean_length_of_burst_bis, density_bis, pisolated_loss_bis, mean_length_of_good_burst_bis );
+			if (gmin_max>0) 
+				external_blind_check (length, gmin_max, sequence, sample, ploss_ter, mean_length_of_burst_ter, density_ter, pisolated_loss_ter, mean_length_of_good_burst_ter);
+		}
+	 }
+	if (num_of_samples>1) {
+		char filename[100];
+		sprintf(filename,"%s_P%d_B%d_EPStatistics.csv", filename_common, atoi(argv[5]), atoi(argv[6]) );
+		printf("INTERNAL CHECK STATISTICS\n");
+		data.open(filename, ios::out);
+		data << "Internal Check Statistics\n";
+		data << "N Samples; Length of vector; Model;";
+		if (strcmp (argv[4],"loss_GI") == 0){
+		  data << "PLoss; Burst Legth;\n";
+		}
+		data <<num_of_samples<<";";
+		data <<length<<";";
+		data <<argv[4]<<";";		//chosen model
+		data <<fixed<<atof(argv[5])<<";";	//saving internal check values
+		if(argv[6]) 
+			data <<atof(argv[6])<<";"; 
+		else 
+			data <<"N/A"<<";";
+		if((argv[6]) && (argv[7])) 
+			data <<atof(argv[7])<<";"; 	
+		else 
+			data <<"100"<<";";
+		if((argv[6]) && (argv[7]) && (argv[8])) 
+			data <<atof(argv[8])<<";"; 
+		else 
+			data <<"0"<<";";
+		if((argv[6]) && (argv[7]) && (argv[8]) && (argv[9])) 
+			data <<atof(argv[9])<<";"; 
+		else 
+			data <<"N/A"<<";";
+		data <<"\n";
+		data.close();
+		statistics(filename, num_of_samples, ploss_bis, mean_length_of_burst_bis, density_bis, pisolated_loss_bis, mean_length_of_good_burst_bis, 2, 2);
+	}
+	if ((gmin_max>0) && (num_of_samples>1))	{
+		char filename[100];
+		sprintf(filename,"%s_P%d_B%d_EPStatistics.csv", filename_common, atoi(argv[5]), atoi(argv[6]) );
+		if (gmin_max>0) 
+			printf("\nEXTERNAL CHECK STATISTICS\n");
+		for (int gmin=2;gmin<gmin_max*2;gmin=gmin*2) {
+			printf("\nFor gmin=%d",gmin);
+			data.open("statistics.dat", ios::out | ios::app);
+			data <<"external"<<";"<<gmin<<";";
+			data <<num_of_samples<<";";
+			data <<length<<";";
+			data <<argv[4]<<";";	//chosen model
+			data <<fixed<<atof(argv[5])<<";";
+			if(argv[6]) 
+				data <<atof(argv[6])<<";"; 
+				else data <<"N/A"<<";";
+			if((argv[6]) && (argv[7]))
+				data <<atof(argv[7])<<";"; 
+			else 
+				data <<"N/A"<<";";
+			if((argv[6]) && (argv[7]) && (argv[8])) 
+				data <<atof(argv[8])<<";"; 
+			else 
+				data <<"N/A"<<";";
+			if((argv[6]) && (argv[7]) && (argv[8]) && (argv[9])) 
+				data <<atof(argv[9])<<";"; 
+			else 
+				data <<"N/A"<<";";
+			data.close();
+			statistics(filename, num_of_samples, ploss_ter, mean_length_of_burst_ter, density_ter, pisolated_loss_ter, mean_length_of_good_burst_ter, gmin, gmin_max);	
+		}
+	}
+	fprintf(stdout,"********************************* Exit SQNGEN *********************************\n");
+}