#include <stdio.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <string.h>
#include <stdbool.h>
#include <time.h>

#define frand() ((float)rand() / (RAND_MAX + 1.0))

#define NUMBYTESOFDATAPAIR 8
#define NUMBITSOFNEURON 32

//#define FOLDERPATH "C:/Users/Joachim Reif/Workbench/Git Repositories/fpga-webinterface/textfiles/output/"

typedef struct event_tt
{

   int tag;
   float t;

} event_t;

typedef struct ev_tt
{
   int size;
   float t[1024];
   struct ev_tt *next;

} ev_t;

typedef struct ex_tt
{
   int tag;
   float rate;
   int *bins; /* used for correlation */
   int size;
   float sumXi; /* sum of pulses        */
   float nXi;   /* number of timesteps  */
   int nEvents;
   int nBins;

   struct ev_tt *theList;
   struct ev_tt *topList;
   struct ex_tt *next;
} ex_t;

typedef struct sys_tt
{

   int modulus;
   ex_t **hTable;

   float h;
   int size;
   int first_id;
   int first_index;
   int last_id;
   int last_index;
   int n;
   char *label;
   char *variable;
   float tLim;
   int precision;
   int exponent;
   int num_time_steps;

   float BW_PEARSON_CORR;
   float BW_PULSE_RATE;
   float BW_CV_ISI;
   float BIN_PEARSON;

   int PR_SET;
   bool compare;
   char *stdFolder;

} sys_t;

char *cpyStr(char *in)
{

   int size;
   char *b;

   for (size = 0; in[size]; size++)
      ;

   b = (char *)malloc(sizeof(char) * (size + 1));

   for (size = 0; in[size]; size++)
   {
      b[size] = in[size];
   }
   b[size] = 0;

   return b;
}

int hashIt(int dataIn, sys_t *sim)
{

   int idx;
   int modulus = sim->modulus;

   idx = dataIn % modulus;

   return idx;
}

ex_t *newEntry(int tag, float t)
{

   ex_t *entry = (ex_t *)malloc(sizeof(ex_t));
   ev_t *theList = (ev_t *)malloc(sizeof(ev_t));

   theList->next = NULL;
   theList->size = 1;
   theList->t[0] = t;

   entry->next = NULL;
   entry->theList = entry->topList = theList;
   entry->tag = tag;
   entry->bins = NULL;

   return entry;
}

void insertValueInList(float* list,int sizeList, float value)
{
	float swapper = 0.0;
	for (int i = 0; i < sizeList+1; ++i)
	{
		if (list[i] == value)
		{
			return ;
		}
		else if (list[i] >value)
		{
			swapper = list[i];
			list[i] = value;
			value = swapper;
		}
	}
	list[sizeList + 1] = value;
	return;
}

void addEntry(float t, ex_t *entry)
{

   ev_t *topList = entry->topList;
   ev_t *nextList;
   int s;

   if (topList->size == 1024)
   {
      entry->topList = topList->next = nextList = (ev_t *)malloc(sizeof(ev_t));

      nextList->next = NULL;
      nextList->size = 1;
      nextList->t[0] = t;
   }
   else
   {
      s = topList->size;
      topList->t[s++] = t;
      topList->size = s;

	  //insertValueInList(topList->t, s, t);
	  //topList->size += 1;
   }
}

void insertEvent(event_t *aEvent, sys_t *sim)
{

   int idx;

   int tag = aEvent->tag;
   float t = aEvent->t;
   ex_t **hTable = sim->hTable;
   ex_t *entry;

   idx = hashIt(tag, sim);

   entry = hTable[idx];

   if (entry == NULL)
   {
      /* new entry */
      hTable[idx] = newEntry(tag, t);
   }
   else
   {

      while (entry)
      {
         if (entry->tag == tag)
         {
            addEntry(t, entry);
            break;
         }
         if (entry->next == NULL)
         {
            entry->next = newEntry(tag, t);
            break;
         }
      }
   }
}

void initHash(sys_t *sim)
{

   int modulus = 99991;
   int i;
   ex_t **hTable;

   hTable = (ex_t **)malloc(sizeof(ex_t *) * modulus);

   for (i = 0; i < modulus; hTable[i++] = NULL)
      ;

   sim->modulus = modulus;
   sim->hTable = hTable;
}

void assignSize(char *scanLine, sys_t *sim)
{
   int size;
   sscanf(scanLine, "%d", &size);
   sim->size = size;
}

void assignDt(char *scanLine, sys_t *sim)
{
   float h;
   sscanf(scanLine, "%f", &h);
   sim->h = h;
}

void assignN(char *scanLine, sys_t *sim)
{
   int n;
   sscanf(scanLine, "%d", &n);
   sim->n = n;
}

void assignFix(char *scanLine, sys_t *sim)
{
   int first_index;
   sscanf(scanLine, "%d", &first_index);
   sim->first_index = first_index;
}

void assignFid(char *scanLine, sys_t *sim)
{
   int first_id;
   sscanf(scanLine, "%d", &first_id);
   sim->first_id = first_id;
}

void assignLix(char *scanLine, sys_t *sim)
{
   int last_index;
   sscanf(scanLine, "%d", &last_index);
   sim->last_index = last_index;
}

void assignLid(char *scanLine, sys_t *sim)
{
   int last_id;
   sscanf(scanLine, "%d", &last_id);
   sim->last_id = last_id;
}

void assignLBL(char *scanLine, sys_t *sim)
{
   sim->label = cpyStr(scanLine);
}

void assignVar(char *scanLine, sys_t *sim)
{
   sim->variable = cpyStr(scanLine);
}

void assignPrc(char *scanLine, sys_t *sim)
{
   int precision;
   sscanf(scanLine, "%d", &precision);
   sim->precision = precision;
}

void assignExp(char *scanLine, sys_t *sim)
{
   int exponent;
   sscanf(scanLine, "%d", &exponent);
   sim->exponent = exponent;
}

void assignNumTimeSteps(char *scanLine, sys_t *sim)
{
   int num_time_steps;
   sscanf(scanLine, "%d", &num_time_steps);
   sim->num_time_steps = num_time_steps;
}


size_t read_int32(FILE *ptr, int *result)
{
   return fread(result, 4, 1, ptr);
}

size_t read_utf8(FILE *ptr, int length, char *result)
{
   result[length] = 0;
   return fread(result, 1, length, ptr);
}

void printSummary(sys_t* sim)
{
   printf("----------------------------------\n");
   printf("-----    SUMMARY HEADER      -----\n");
   printf("----------------------------------\n");
   if (sim->label)
      printf(" label        = '%s'\n", sim->label);
   else
      printf(" label        = <NULL>\n");
   printf(" popSize      = %d [neurons]\n", sim->size);
   if (sim->variable)
      printf(" variable     = %s\n", sim->variable);
   else
      printf(" variable     = <NULL>\n");
   printf(" h            = %2.4f [ms]\n", sim->h);
   printf(" first_id     = %d\n", sim->first_id);
   printf(" first_index  = %d\n", sim->first_index);
   printf(" last_id      = %d\n", sim->last_id);
   printf(" last_index   = %d\n", sim->last_index);
   printf(" n            = %d [spike events]\n", sim->n);
   printf(" precision    = %d [bits mantissa]\n", sim->precision);
   printf(" exponent     = %d [bits exponent]\n", sim->exponent);
   printf(" numTimeSteps = %d \n", sim->num_time_steps);
   printf("----------------------------------\n");
   printf("-----   SUMMARY PARAMETER    -----\n");
   printf("----------------------------------\n");
   printf(" tLim         = %2.2f [ms]\n", sim->tLim);
   printf(" BW_PEARSON_CR= %2.9f [1]\n", sim->BW_PEARSON_CORR);
   printf(" BW_PULSE_RATE= %2.9f [1/s]\n", sim->BW_PULSE_RATE);
   printf(" BW_CV_ISI    = %2.9f [1]\n", sim->BW_CV_ISI);
   printf(" PR_SET       = %d [neurons per population]\n", sim->PR_SET);
   printf(" BIN_PEARSON  = %2.2f [ms]\n", sim->BIN_PEARSON);
   printf("----------------------------------\n");
   printf("----------------------------------\n");
}


int countTags(sys_t *sim)
{

   char toolStr[] = "countTags";

   int i;
   int j;
   int idx;
   int c = 0;
   int cnt = 0;
   ex_t *entry;
   ev_t *theList;
   float tLim = sim->tLim;
   float tEnd;
   int size;
   int nEmpty;

   nEmpty = 0;
   for (i = 0; i < sim->modulus; i++)
   {
      entry = sim->hTable[i];
      while (entry)
      {
         /* ------------------------ */
         /* spikes of a neuron       */

         size = 0;
         theList = entry->theList;
         while (theList)
         {
            for (j = 0; j < theList->size; j++)
            {
               tEnd = theList->t[j];
               if (tEnd > tLim)
               {
                  size++;
               }
            }
            theList = theList->next;
         }

         /* -------------------------------------- */
         /* number of spikes                       */

         entry->size = size;

         if (size != 0)
         {
            entry->bins = (int *)malloc(sizeof(int) * size);
            entry->size = size;         /* number of spikes */
            entry->sumXi = (float)size; /* dto.             */
            float simTime = sim->num_time_steps * sim->h;
            entry->nXi = (simTime - tLim) / sim->BIN_PEARSON; /* number of bins   */

            c = 0;
            theList = entry->theList;
            while (theList)
            {
               for (j = 0; j < theList->size; j++)
               {
                  if (theList->t[j] > tLim)
                  {
                     /* which bin is it */
                     //printf("%d: %f\n",entry->tag,  theList->t[j]);
                     idx = (int)((theList->t[j] - tLim) / sim->BIN_PEARSON);

					 if (entry->tag == 4 && idx == 0)
					 {
						 printf("time %f\n", theList->t[j]);
					 }
                     entry->bins[c++] = idx;
                  }
               }
               theList = theList->next;
            }

            if (c != size)
            {
               printf("(%s) BUG!  c=%d  size=%d\n", toolStr, c, size);
               exit(1);
            }

            for (j = 1; j < size; j++)
            {
               if (entry->bins[j - 1] >= entry->bins[j])
               {
                 printf("Bug! neuron %d   %d(%d): %d %d\n", entry->tag, j, size, entry->bins[j - 1], entry->bins[j]);
                 exit(1);
               }
            }
         }
         else
         {
            nEmpty++;
         }

         entry = entry->next;
         cnt++;
      }
   }

   printf("(%s) active neurons  : %d\n", toolStr, cnt);
   printf("            inactive neurons: %d\n", sim->size - cnt);
   printf("            spike-less      : %d (have spikes, but outside the valid window t > %3.1f [ms])\n", nEmpty, tLim);

   if(tLim + 0.1 >= sim->num_time_steps * sim->h){
      printf("(ERROR) tLim = %3.1f ms is bigger too close to simulation time = %3.1f. Adjust tLim maybe!\n", tLim, sim->num_time_steps * sim->h);
      exit(0); 
   }

   if(cnt - nEmpty < sim->PR_SET){
      printf("(ERROR) Number of active neurons is smaller than PR_SET!");
      exit(0);
   }

   return cnt;
}

static int cmpFlt(const void *p1, const void *p2)
{

   float *p1_ = (float *)p1;
   float *p2_ = (float *)p2;

   if ((*p1_) > (*p2_))
      return 1;
   if ((*p1_) < (*p2_))
      return -1;
   return 0;
}

float gKernel(float tArg, float n, float h)
{

   float f0 = 1.0 / (n * h * sqrt(2.0 * 3.1415926));

   return (f0 * exp(-0.5 * pow(tArg / h, 2.0)));
}

void updateBins(float t, float n, float h, float b, float *bins, int nBins, int iMin)
{

   /* bin size is b */

   int i;
   float tMax = t + 4.0 * h;
   float tMin = t - 4.0 * h;
   float tArg;
   float tCurr;

   /* get the mean index */

   i = -iMin + (tMin) / b;

   tCurr = tMin;
   for (; tCurr <= tMax;)
   {
      tArg = tCurr - t;

      if ((i >= 0) && (i < (nBins - iMin)))
         bins[i] += gKernel(tArg, n, h);

      i++;
      tCurr += b;
   }
}

void friedmanDraconis(char *fr, float *allT, int nSamples, float h, float RATE_MAX, float RATE_MIN, FILE *fd)
{

   char toolStr[] = "friedmanDiaconis";
   int i;
   int i25;
   int i75;
   float IQR;
   float b;
   float *bins;
   float bSum;
   float tLast;
   int nBins;
   int iMin;
   float n = (float)nSamples;

   qsort(allT, nSamples, sizeof(float), cmpFlt);

   i25 = nSamples >> 2;
   i75 = i25 + (nSamples >> 1);

   IQR = 2.0 * (allT[i75] - allT[i25]) / (pow(n, 1.0 / 3.0));

   if (IQR == 0.0)
   {
      printf("(%s:%s) WARNING. setting bin width to bandwidth\n", toolStr, fr);
      IQR = h;
   }

   b = IQR;

   nBins = (int)(RATE_MAX / b);
   iMin = (int)(RATE_MIN / b); /* could be negative */
   
   //printf("nBins: %f, iMin: %f\n", RATE_MAX, RATE_MIN);

   if ((nBins - iMin) > 1024.0)
   {
      b = (RATE_MAX - RATE_MIN) / 1024.0;
      nBins = (int)(RATE_MAX / b);
      iMin = (int)(RATE_MIN / b); /* could be negative */
   }

   bins = (float *)malloc(sizeof(float) * (nBins - iMin));

   for (i = 0; i < (nBins - iMin); bins[i++] = 0.0)
      ;

   /* ------------------------------------------ */
   /* sample the samples using kernels           */

   for (i = 0; i < nSamples; i++)
   {
      updateBins(allT[i], n, h, b, bins, nBins, iMin);
   }

   /* -----------------------------------------  */
   /* output to disc                             */

   if (fd)
   {
      bSum = 0.0;

      for (i = 1; i < (nBins - iMin); i++)
      {
         /* integrate trapezoid */
         bSum += 0.5 * ((bins[i] + bins[i - 1]) * b);
      }

      tLast = RATE_MIN;
      for (i = 1; i < (nBins - iMin); i++)
      {
         float val = bins[i] / bSum;
		 //printf("%2.5f %2.5f\n", tLast, val);
         fprintf(fd, "%2.5f %2.5f\n", tLast, val);
         tLast += b;
      }
   }

   free(bins);
}

void getRates(sys_t *sim, int nActive, int nInactive, FILE *fd)
{

   char toolStr[] = "RATES";
   int i;
   int j;
   int idx;
   int c = 0;
   float tLast;
   float RATE;
   float RATE_MAX = 0.0;
   float RATE_MIN = 1000.0;
   float tLim = sim->tLim;
   ex_t *entry;
   ev_t *theList;
   float *allT;

   int nSamples = nActive + nInactive;

   allT = (float *)malloc(sizeof(float) * nSamples);

   for (i = 0; i < nSamples; allT[i++] = 0.0)
      ;

   idx = 0;
   RATE = 0.0;
   for (i = 0; i < nInactive; i++)
   {
      allT[idx++] = RATE;
   }

   for (i = 0; i < sim->modulus; i++)
   {
      entry = sim->hTable[i];
      c = 0;
      while (entry)
      {
         theList = entry->theList;
         tLast = 0.0;
         while (theList)
         {
            for (j = 0; j < theList->size; j++)
            {
               tLast = theList->t[j];
               if (tLast > tLim)
               {
                  c++;
               }
            }
            theList = theList->next;
         }

         /* everything in [ms] */

         RATE = 1000.0 * ((float)c) / (tLast - tLim);
		 //printf("Rate: %f\n", RATE);
		 if(RATE_MAX < RATE)
			 RATE_MAX = RATE;
		 if(RATE_MIN > RATE)
			 RATE_MIN = RATE;

         allT[idx++] = RATE;

         entry = entry->next;
      }
   }

   if (nSamples != idx)
   {
      printf("(%s) BUG! sample error: nSamples:%d != %d:idx\n", toolStr, nSamples, idx);
      exit(1);
   }
   else
   {
      /* */
   }

   friedmanDraconis(toolStr, allT, nSamples, sim->BW_PULSE_RATE, RATE_MAX, RATE_MIN, fd);

   free(allT);
}

void getCV(sys_t *sim, int nActive, int nInactive, FILE *fd)
{

   char toolStr[] = "CV-ISI";
   int i;
   int j;
   int idx;
   float t;
   float CV;
   float ISI;
   float ISInn;
   float ISImu;
   float ISIsig;
   float RATE_MIN = 0.0;
   float RATE_MAX = 2.0;
   float tLim = sim->tLim;
   float tLast;
   ex_t *entry;
   ev_t *theList;
   float *allT;

   int nSamples = nActive + nInactive;

   allT = (float *)malloc(sizeof(float) * nSamples);

   for (i = 0; i < nSamples; allT[i++] = 0.0)
      ;

   idx = 0.0;
   CV = 0.0;
   for (i = 0; i < nInactive; i++)
   {
      allT[idx++] = CV;
   }

   for (i = 0; i < sim->modulus; i++)
   {
      entry = sim->hTable[i];
      ISImu = ISIsig = ISInn = 0.0;
      while (entry)
      {
         theList = entry->theList;
         int listCount = 0;
         tLast = -1.0;
         while (theList)
         {
            if (listCount == 0 && (theList->size <= 2))
            {
               break;
            }
            for (j = 0; j < theList->size; j++)
            {

               t = theList->t[j];
               ISI = t - tLast;
               if ((t > tLim) && (tLast >= 0.0))
               {
                  ISInn += 1.0;
                  ISImu = ISImu + (pow(ISI, 1.0) - ISImu) / ISInn;
                  ISIsig = ISIsig + (pow(ISI, 2.0) - ISIsig) / ISInn;
               }
               tLast = t;
            }
            theList = theList->next;
            ++listCount;
         }

         ISIsig = (ISIsig - pow(ISImu, 2.0));

         if (ISIsig >= 0.0)
            ISIsig = sqrt(ISIsig);
         else
            ISIsig = 0.0;

         if (ISImu < -0.0000001 || ISImu > 0.0000001)
         {
            CV = ISIsig / ISImu;
         }
         else
         {
            CV = 0;
         }

         allT[idx++] = CV;

         entry = entry->next;
      }
   }

   if (nSamples != idx)
   {
      printf("(%s) BUG! sample error: nSamples:%d != %d:idx\n", toolStr, nSamples, idx);
      exit(1);
   }
   else
   {
      /* */
   }

   friedmanDraconis(toolStr, allT, nSamples, sim->BW_CV_ISI, RATE_MAX, RATE_MIN, fd);

   free(allT);
}

int correlateBins(ex_t *A, ex_t *B)
{

   int idxA, idxB;

   int aA, aB;
   int c;

   idxA = idxB = 0;

   aA = A->bins[idxA++];
   aB = B->bins[idxB++];
   c = 0;

   while ((idxA <= A->size) && (idxB <= B->size))
   {

      if (aA > aB)
      {
         aB = B->bins[idxB++];
      }
      else if (aA < aB)
      {
         aA = A->bins[idxA++];
      }
      else
      {
         c++;
         aA = A->bins[idxA++];
         aB = B->bins[idxB++];
      }
   }

   return (c);
}

void getPC(sys_t *sim, int nActive, int nInactive, FILE *fd)
{

   /* pearson correlation coefficient */

   char toolStr[] = "pearsonCorrelation";
   ex_t *entry;
   ex_t *entryA;
   ex_t *entryB;

   ex_t **allE = (ex_t **)malloc(sizeof(ex_t *) * nActive);
   int i;
   int j;
   int i1, i2;
   int c;

   float cAA;
   float cBB;
   float cAB;
   float *theF = (float *)malloc(sizeof(float) * sim->PR_SET * sim->PR_SET);
   int nSamples;


   j = 0;
   for (i = 0; i < sim->modulus; i++)
   {
      entry = sim->hTable[i];

      while (entry)
      {
         if (entry->size != 0)
         {
            allE[j++] = entry;
         }
         entry = entry->next;
      }
   }

   if (j > nActive)
   {
      printf("(getPC) BUG!\n");
      exit(1);
   }
   /* shuffle */

   nActive = j;

   //printf("(%s) ... shuffle %d pairs\n",toolStr,nActive);

   for (j = 0; j < nActive; j++)
   {

      do
      {
         i1 = rand() % nActive;
         i2 = rand() % nActive;
      } while (i1 == i2);

      /* swap */
      entry = allE[i1];
      allE[i1] = allE[i2];
      allE[i2] = entry;
   }
   
   
   nSamples = 0;
   for (i = 0; i < sim->PR_SET; i++)
   {
      entryA = allE[i];
      cAA = sqrt(entryA->sumXi * (1.0 - ((entryA->sumXi) / entryA->nXi)));

      for (j = 0; j < sim->PR_SET; j++)
      {
         if (i != j)
         {
            entryB = allE[j];
            cBB = sqrt(entryB->sumXi * (1.0 - ((entryB->sumXi) / entryB->nXi)));

            c = correlateBins(entryA, entryB);

            cAB = (((float)c) - entryB->sumXi * ((entryA->sumXi) / entryA->nXi)) / (cAA * cBB);

			//printf("cab: %f\n", cAB);

            theF[nSamples++] = cAB;
         }
      }
   }

   printf("(%s) %d nSamples processed\n",toolStr,nSamples);

   friedmanDraconis(toolStr, theF, nSamples, sim->BW_PEARSON_CORR, 0.15, -0.05, fd);
   
   free(allE);
   free(theF);
   
}

void freeEvents(ev_t *theList)
{

   ev_t *dList;

   if (theList)
   {
      dList = theList;
      theList = theList->next;
      free(dList);
   }
}

void freeHash(int modulus, ex_t **hTable)
{

   int i;
   ex_t *entry, *nextEntry;
   ;

   if (hTable)
   {
      for (i = 0; i < modulus; i++)
      {
         entry = hTable[i];
         while (entry)
         {
            freeEvents(entry->theList);
            nextEntry = entry;
            entry = entry->next;
            if (nextEntry->bins)
               free(nextEntry->bins);
            else
            {
               //printf("no bins for %d\n",nextEntry->tag);
            }
            free(nextEntry);
         }
      }
      free(hTable);
   }
}

void freeSim(sys_t *sim)
{

   if (sim)
   {

      freeHash(sim->modulus, sim->hTable);
      if (sim->label)
         free(sim->label);
      if (sim->variable)
         free(sim->variable);
      free(sim);
   }
}

void showError(void)
{
   printf("(showError) wrong arument/option\n");
   exit(1);
}

float roundData(float number, int fractionalDigits){
	float num = number * (pow(10, fractionalDigits));
	if(num >= 0){
		num = num + 0.5;
	} else{
		num = num - 0.5;
	}
	int tempNumber = (int) num;
	number = (float)tempNumber/(pow(10,fractionalDigits));
	return number;
}

size_t read_bytesX(FILE *ptr, unsigned char* result, int numberOfBytes)
{
	return fread(result, numberOfBytes, 1, ptr);
}

void readBinaryData(char* filename, int numBitsOfNeuron, int numBytesOfDataPair, sys_t* sim)
{	
	
	int numTimeSteps = sim->num_time_steps;
	int maxNeuronId = sim->last_index;
	
	FILE *fd = fopen(filename, "rb");
	int zeroCounter = 0;
	
	if(fd)
	{
		char* data = (char*)calloc(numBytesOfDataPair, sizeof(char));
		char* bits = (char*)malloc(numBytesOfDataPair*8*sizeof(char));
		while(!feof(fd))
		{
			if(read_bytesX(fd, data, numBytesOfDataPair))
			{
				int neuronID = 0;
				int timeStep = 0;
					
				for(int i=0; i<numBytesOfDataPair; ++i)
				{
					for(int b=0; b<8; ++b)
					{
						bits[8*i+b] = (data[i] >> b) & 1;
					}
				}
					
				for(int i=0; i<numBitsOfNeuron; ++i)
				{
					neuronID += bits[i] << i; 
				}
				for(int i=numBitsOfNeuron; i<numBytesOfDataPair*8; ++i)
				{
					timeStep += bits[i] << (i-numBitsOfNeuron);
				}
				
				float tt = (float)timeStep/10;
				
				event_t aEvent;
				aEvent.t = tt;
				aEvent.tag = neuronID;
				
				insertEvent(&aEvent, sim);
				
				if(neuronID > maxNeuronId)
					maxNeuronId = neuronID;
				
				if(timeStep > numTimeSteps)
					numTimeSteps = timeStep;
			}
		}
		
		//printf("Counter: %d\n", zeroCounter);
		free(data);
		free(bits);
		fclose(fd);
		
		sim->size = maxNeuronId+1;
	
		printf("MaxNeuronId %d, NumTimeSteps %d\n", maxNeuronId, numTimeSteps);
	
		sim->first_index = 0;
		sim->last_index = maxNeuronId;
		
		sim->first_id = 0;
		sim->last_id = maxNeuronId;
		
		sim->num_time_steps = numTimeSteps;
	}
	else
	{
		fprintf(stderr, "The File %s do not exist.\n", filename);
		exit(EXIT_FAILURE);
	}
	
}

int lenghtOfArray(const char* arr)
{
	int size = 0;
	while(*arr)
	{
		++size;
		++arr;
	}
	return size;
}


void freeStringArray(char** str, int numElements)
{
	free(str);
}

int create(int numFiles, char* files, char* outDir)
{
	printf("Start calculation\n");
	
	if(numFiles < 1){
		fprintf(stderr, "ERROR: No File selected\n");
		exit(EXIT_FAILURE);
	}
	//printf("outDir: %s\n", outDir);
	
	
	clock_t start, end;
	double cpu_time_used;
	
	start = clock();
	//char* outDir = ".";
	float tLim = 0;
	float bw_pearson = 0.0002;
	//float bw_pearson = 0.002;
	float bw_pulse_rate = 0.3;
	float bw_cv_isi = 0.04;
	float bin_pearson = 2.0;
	float pr_set = 100;
	
	float timeStep = 0.1;
	int precision = 52;
	int exponent = 11;
	
	sys_t *sim;
	
	int nActive;
	int nInactive;
	int nNeurons;
	int isReference = 0;
	int n;
   
	FILE *fd;
	char fdName[4096];
	FILE *fdBin;
	FILE *fdCV;
	FILE *fdPC;
	
	sim = (sys_t*)malloc(sizeof(sys_t));
	
	sim->h = 0.100; /*  [ms]  */
	sim->n = -1;
	sim->size = -1;
	sim->first_index = -1;
	sim->first_id = -1;
	sim->last_index = -1;
	sim->last_id = -1;
	sim->num_time_steps = -1;
	sim->variable = NULL;
	sim->label = NULL;
	sim->tLim = tLim; /*  [ms];  */
	sim->precision = precision;
	sim->exponent = exponent;
	
	sim->BW_PEARSON_CORR = bin_pearson;
	sim->BW_PULSE_RATE = bw_pulse_rate;
	sim->BW_CV_ISI = bw_cv_isi;
	sim->BIN_PEARSON = bin_pearson;
	sim->PR_SET = pr_set;
	
	sim->size = 0;
	sim->first_index = 0;
	sim->last_index = 0;
	sim->first_id = 0;
	sim->last_id = 0;
	sim->num_time_steps = 0;
	
	
	initHash(sim);
	
	printf("Reading Data\n");

	readBinaryData(files, NUMBITSOFNEURON, NUMBYTESOFDATAPAIR, sim);

	printf("Free memory\n");
	
	if (sim->num_time_steps == -1)
	{
		fprintf(stderr, "ERROR: num_time_steps must either be set in header or as command line parameter!");
		exit(EXIT_FAILURE);
	}
	
	if (sim->size > 0)
	{
		nNeurons = sim->size;
	}
	else
	{
		fprintf(stderr, "(main) ERROR: unable to determine number of neurons!\n");
		exit(EXIT_FAILURE);
	}
	
	//sim->PR_SET = ((int)(nNeurons*0.3)>200)?200:(int)(nNeurons*0.3);
	nActive = countTags(sim);
	nInactive = nNeurons - nActive;
	
	if (outDir == NULL)
	{
		fprintf(stderr, "(main) ERROR: No output directory is selected!\n");
		exit(EXIT_FAILURE);
	}
	else
	{
		sprintf(fdName, "%s/pr_output.txt", outDir);
	}
	
	fdBin = fopen(fdName, "wb+");
	if (fdBin == NULL)
	{
		fprintf(stderr, "ERROR: unable to create %s for write!\n", fdName);
		exit(EXIT_FAILURE);
	}
	
	printf("Start calculating Puls Rate\n");
	getRates(sim, nActive, nInactive, fdBin);
	
	printf("(main) %s created\n", fdName);
	fclose(fdBin);
	
	sprintf(fdName, "%s/cv_output.txt", outDir);
	fdCV = fopen(fdName, "wb+");
	if(fdCV == NULL)
	{
		fprintf(stderr, "ERROR: unable to create %s for write!\n", fdName);
		exit(EXIT_FAILURE);
	}
	
	printf("Start calculating Coefficient of Variation\n");
	getCV(sim, nActive, nInactive, fdCV);
	printf("(main) %s created\n", fdName);
	fclose(fdCV);
	
	sprintf(fdName, "%s/pc_output.txt", outDir);
	fdPC = fopen(fdName, "wb+");
	if(fdPC == NULL)
	{
		fprintf(stderr, "ERROR: unable to create %s for write!\n", fdName);
		exit(EXIT_FAILURE);
	}
	
	printf("Start calculating Pearson Correlation\n");
	getPC(sim, nActive, nInactive, fdPC);
	("(main) %s created\n", fdName);
	fclose(fdPC);
	
	freeSim(sim);
	
	end=clock();
	cpu_time_used = ((double)(end - start)) /CLOCKS_PER_SEC;
	printf("Time needed in total: %f\n", cpu_time_used);
	
	return EXIT_SUCCESS;
}
