/*

DIETA1.C

A PC-compatible program to calculate two indices of intrapopulation diet variation
based on network theory
 
*/

#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <fcntl.h>
#include <sys\types.h>
#include <sys\stat.h>
#include <io.h>
#include <math.h>

float moduloS(float x);
float moduloB(float x);

#define SIM                   1
#define NAO                   0
#define MAX_DIMENSAO          2000
#define MAX_ARQUIVOS          10000
#define MAX_SIMULACAO         10000

#define CR                    0x0d // '\r'
#define LF                    0x0a // '\n'
#define TAB                   0x09 // '\t'
#define ESPACO                0x20 // ' '
#define VIRGULA               0x2c // ','
#define PONTO                 0x2e // '.'
#define PONTO_VIRGULA         0x3b // ';'

float Matriz[MAX_DIMENSAO][MAX_DIMENSAO];
float R[MAX_DIMENSAO][MAX_DIMENSAO];
float P[MAX_DIMENSAO][MAX_DIMENSAO];
float POriginal[MAX_DIMENSAO][MAX_DIMENSAO];
float w[MAX_DIMENSAO][MAX_DIMENSAO];
float Cw[MAX_DIMENSAO];
float Ei[MAX_DIMENSAO];
float EiVarE[MAX_DIMENSAO];
float Somai[MAX_DIMENSAO];
float Somaj[MAX_DIMENSAO];
unsigned char B[MAX_DIMENSAO][MAX_DIMENSAO];


float WmedioSimulacao[MAX_SIMULACAO];
float ESimulacao[MAX_SIMULACAO];
float CWSimulacao[MAX_SIMULACAO];
float CwSSimulacao[MAX_SIMULACAO];

unsigned int ProbIni[MAX_DIMENSAO];
unsigned int ProbFim[MAX_DIMENSAO];
unsigned int Probabilidadej[MAX_DIMENSAO];

float ValorLido, temp, temp2, W, Wmedio, Wmax, Wmij, Wmjk, Wmki, Wraiz3, E, Ecalculado, F, CW, CwS, SomaMatriz, SomaR, Si, carga;
unsigned int PE, PCwS;
float EObservado, CwSObservado, Emedio, varE, mediaCwS;

char NomesLinhas[MAX_DIMENSAO][15];


int  linha, linmax, linmaxOriginal, numread, flag, indice, indice_arquivos, dimensao, TipoDadoMatriz;
int  imprimeBoot, imprimeCw, imprimeE, rodaSaramaki; 
int  coluna, colmax, colmaxOriginal, virgula, expo, simulacao, formula, replicas, weightFactor, Ki, rodada, fase2;

int linhaRetirada, calculandoVarE;

char buf_aux[50];

char lixo[12];
FILE *stream, *arqin;

char nomeA[MAX_PATH]="MATRIZ.TXT";
char nome[MAX_PATH]="matriz.TXT";
char nomeX[MAX_PATH]="Indices[matriz].TXT";
char nomeV[MAX_PATH]="VarE[matriz].TXT";

struct tm *Temp;
time_t tempo;

unsigned int base_probabilidade, ProbabilidadeSorteada;
float alatorio_float;
unsigned int aleatorio;

static int i, j, k;
char fileName[MAX_PATH];
char *file = &fileName[0];


//////////////////////////////////////////////////////////////////////////////

BOOL getFiles(char *filename) // Parm shall point to an array of MAX_PATH chars
{
WIN32_FIND_DATA buffer;
LPWIN32_FIND_DATA fileInfo = &buffer;
static HANDLE hdl = 0;

   if (filename == NULL)
      return FALSE;

again:

   if (hdl == 0) // First call in the interaction
   { if ((hdl = FindFirstFile(filename,fileInfo)) == INVALID_HANDLE_VALUE)
	    return (BOOL)hdl = 0; // Answer FALSE and prepare to next interaction
   }
   else // Subsequent call
    if (!FindNextFile(hdl,fileInfo)) // Get next file
 	      return (BOOL)hdl = 0; // Answer FALSE and prepare to next interaction

   if (strncmp("P[", fileInfo->cFileName, 2) == 0) goto again;

   strcpy(filename,fileInfo->cFileName); // Copy file name
   return TRUE;
}

//////////////////////////////////////////////////////////////////////////////

void main()
{
static int rando,passatempo,anterior;
static unsigned char vetor[4] = {'|', '/', '-', '\\'}; 

   /* Seed the random-number generator with current time so that
    * the numbers will be different every time we run.
    */
   srand((unsigned)time(NULL));

   indice_arquivos=0;
   lixo[0]=0x30;

   printf("_______________DIETA1_____________\n\n");
   printf("Date:    07/22/2009\n");
   printf("Authors: Paulo Guimaraes\n");
   printf("         Paulo R. Guimaraes Junior\n");
   printf("         Marcio S. Araujo\n");
   printf("__________________________________\n\n");

   printf("Please enter the file containing the data (e.g. datafile.txt)\n");
   printf("- ");
   scanf("%s", fileName);

tenteDeNovo1:
   printf("\nWhat type of data is it?\n");
   printf("(1)Proportions, (2)Integers, (3)Non-integer values - ");
   scanf("%d", &TipoDadoMatriz);
   if ((TipoDadoMatriz > 3) || (TipoDadoMatriz == 0))
   {
      printf("Wrong value, try again...\n");
      goto tenteDeNovo1;
   }

   if (TipoDadoMatriz == 2)
   {
      printf("\nWould you like to run a Monte Carlo bootstrap? (y/n) - ");
      scanf("%s", &lixo[0]);
      if ((lixo[0] == 'Y') || (lixo[0] == 'y'))
      {
         simulacao=SIM;

         printf("\nHow many replicates? (1 to 10000) - ");
         scanf("%d", &replicas);
         if ((replicas == 0) || (replicas > 10000))
         {
            replicas = 1;
         }
tenteDeNovo2:
         printf("\nHow would you like to calculate the population's diet proportions?\n");
         printf("(1) Numerical sum, (2) Average proportion - ");
         scanf("%d", &formula);
         if ((formula == 0) || (formula > 2))
         {
            printf("Wrong value, try again...\n");
            goto tenteDeNovo2;
         }
      }
      else
      {
         simulacao=NAO;
      }
   }

   printf("\nWhat is the 'weight factor'(1 to 100)? - ");
   scanf("%d", &weightFactor);
   if ((weightFactor == 0) || (weightFactor > 100))
   {
      weightFactor = 100;
   }

   lixo[0] = 0x30;

   imprimeBoot=NAO;
   imprimeCw=NAO;
   imprimeE=NAO;

   if (simulacao == SIM)
   {
#if 0
tenteDeNovo3:

   lixo[0] = 0x30;
      printf("\nWhat clustering coefficient you want to use? Type (s) for Saramaki's and (b) for Barrat's - ");
      scanf("%s", &lixo[0]);
      if ((lixo[0] == 'S') || (lixo[0] == 's'))
      {
         rodaSaramaki=SIM;
      }
      else if ((lixo[0] == 'B') || (lixo[0] == 'b'))
      {
         rodaSaramaki=NAO;
      }
      else 
      {
         printf("Wrong option, try again...\n");
         goto tenteDeNovo3;
      }
#endif
//      rodaSaramaki=SIM;

      printf("\nWould you like to print the file Bootstrap.txt (y/n)? - ");
      scanf("%s", &lixo[0]);
      if ((lixo[0] == 'Y') || (lixo[0] == 'y'))
      {
         imprimeBoot=SIM;
      }
   }

   rodaSaramaki=SIM;

   lixo[0] = 0x30;


#if 0 // retirar quando quiser imprimir
   printf("\nWould you like to print the file Cw[Matrix].TXT (y/n)? - ");
   scanf("%s", &lixo[0]);
   if ((lixo[0] == 'Y') || (lixo[0] == 'y'))
   {
      imprimeCw=SIM;
   }

   lixo[0] = 0x30;

   printf("\nWould you like to print the file E[Matriz].TXT (y/n)? - ");
   scanf("%s", &lixo[0]);
   if ((lixo[0] == 'Y') || (lixo[0] == 'y'))
   {
      imprimeE=SIM;
   }
#endif

Processa_Arquivo:

   rodada=0;
   fase2=NAO;
   mediaCwS=0;

   if (getFiles(file) == 0)
   {
      if (indice_arquivos == 0)
      {
         printf("\nFile %s not found or end of file\n", file);
         scanf("%s", &lixo);
      }
      return; // fim de arquivos
   }

   for (i=0; i<MAX_PATH; i++)
   {
      nomeA[i] = fileName[i];
   }
   
   if ((arqin=fopen(nomeA, "rt")) != NULL)
   {
		indice=linha=coluna=linmax=colmax=0;
      virgula=NAO;
      flag=1;
      do
      {
         numread = fread( buf_aux, sizeof( char ), 1, arqin );
         if (numread == 1)
         {
            if (flag == 1)
            {
               if ((buf_aux[0] == ESPACO) || (buf_aux[0] == TAB))
               {
                  NomesLinhas[linha][indice] = '\0';
                  flag=0;
                  ValorLido=0;
                  virgula=NAO;
                  indice=1;
                  continue;
               }
               else
               {
                  NomesLinhas[linha][indice] = buf_aux[0];
                  ++indice;
                  continue;
               }
            }
         	switch(buf_aux[0])
         	{
               case 0x09: // tab
               case ' ':
                  Matriz[linha][coluna] = ValorLido;
                  if (++coluna > 2000)
                  {
                     printf("\nFile has more than 2,000 columns! Enter any letter to abort\n");
                     scanf("%d", &lixo[0]);
                     exit(0);
                  }
                  ValorLido=0;
                  indice=1;
                  virgula=NAO;
                  break;
               case '.':
               case ',':
                  virgula=SIM;
                  break;
					case 0x30: // 0
					case 0x31: // 1
               case 0x32: // 2
               case 0x33: // 3
					case 0x34: // 4
					case 0x35: // 5
               case 0x36: // 6
               case 0x37: // 7
					case 0x38: // 8
					case 0x39: // 9
                  if (virgula)
                  {
                     expo=1;
                     temp = (float)(buf_aux[0]-0x30);
                     for (i=0; i<indice; i++) expo = expo*10;
   						ValorLido += (float)(temp/(float)(expo));
                     ++indice;
                  }
                  else
                  {
                     ValorLido = (ValorLido*10) + (float)(buf_aux[0]-0x30);
                  }
            		break;
            	case 0x0d: // CR
               case 0x0a: // LF
                  Matriz[linha][coluna] = ValorLido;
                  if (++coluna > 2000)
                  {
                     printf("\nFile has more than 2,000 columns! Enter any letter to abort\n");
                     scanf("%d", &lixo[0]);
                     exit(0);
                  }
                  if (++linha > 2000)
                  {
                     printf("\nFile has more than 2,000 rows! Enter any letter to abort\n");
                     scanf("%d", &lixo[0]);
                     exit(0);
                  }
                  if (colmax == 0)
                  {
                     colmax=coluna;
                  }
            		coluna=0;
                  ValorLido=0;
                  virgula=NAO;
                  linmax=linha;
                  // indica que vai ler outro nome de linha
                  flag=1;
                  indice=0;
                  break;
            	default:
            		break;
            }
         }
      }while(numread);
      fclose(arqin);
      if (coluna != 0) ++linmax; // não tem linha em branco no fim...
   }
   else
   {
      printf("\nFile %s not found\n", nomeA);
      scanf("%s", &lixo);
      return;
   }

   for(i=0; i<linmax; i++)
   {
      for(j=0; j<colmax; j++)
      {
         R[i][j] = Matriz[i][j];
      }
   }

   linmaxOriginal = linmax;
   colmaxOriginal = colmax;
   calculandoVarE = NAO;

RepeteCalculoFase2: // Loop para Tipo 2 (Contagens) com Simulação

   for(i=0; i<linmax; i++)
   {
      for(j=0; j<colmax; j++)
      {
         P[i][j] = R[i][j];
      }
   }

   tempo = time(NULL);
   Temp =  localtime(&tempo);
//   if (passatempo != Temp->tm_sec)
   {
      passatempo = Temp->tm_sec;
      ++rando;
      printf("%c\r", vetor[rando%4]);                   
   }

   // Passo 1. Transformando em proporções
   // Se Tipo 2 (Contagens) ou 3 (Numeros Decimais) calcula a matriz proporções
   if (TipoDadoMatriz != 1)
   {
      for (i=0; i<linmax; i++)
      {
         temp=0;
         for (j=0; j<colmax; j++)
         {
            temp += R[i][j];
         }
         for (j=0; j<colmax; j++)
         {
            if (temp == 0) temp=1;
            P[i][j] = R[i][j]/temp;
         }
      }
   }

   if (fase2 == NAO)
   {
      for (i=0; i<linmax; i++)
      {
         for (j=0; j<colmax; j++)
         {
            POriginal[i][j] = P[i][j];
         }
      }
   }
   
   // Passo 2. Calculando PS
   for (i=0; i<linmax; i++)
   {
      for (k=0; k<linmax; k++)
      {
         if (k == i) continue;
         temp=0;
         for (j=0; j<colmax; j++)
         {
            if (P[i][j] > P[k][j])
               temp += (P[i][j] - P[k][j]);
            else
               temp += (P[k][j] - P[i][j]);

         }
         w[i][k] = (float)1.0 - (float)0.5*temp;
      }
   }

   for (i=0; i<linmax; i++)
   {
      w[i][i]=0; // zera a diagonal
   }

   // Passo 3. Calculando a especialização média da população
   W = 0;
   for (i=0; i<linmax; i++)
   {
      Ei[i]=0;
      for (j=0; j<linmax; j++)
      {
         W += w[i][j];
         Ei[i] += w[i][j];
      }
      if (linmax > 1)
         Ei[i] = (float)(Ei[i]/(float)(linmax-1));
      else
         Ei[i] = 0;
   }
   if (linmax > 1)
      Wmedio = (float)(W / (float)(linmax*(linmax-1)));
   else
      Wmedio = 0;

   E = (float)(1.0 - Wmedio);
   
   if (fase2 == NAO)
   {
      EObservado = E;
      PE = 0;
      if (calculandoVarE == SIM)
      {
         EiVarE[linhaRetirada] = E;
         goto retiraProximaLinha;
      }
      else
      {
         Ecalculado = E;
      }
   }
   else
   {
      if (E >= EObservado) ++PE;
   }

   // Passo 4. Calculando o coeficiente de agregação da população
   // 2/(Ki(Ki-1)*Sj,k(Wmij*Wmjk*Wmkj)**1/3, Wmij = wij/Max(Wnm)
   Wmax=0;
   for (i=0; i<linmax; i++)
   {
      for (j=0; j<linmax; j++)
      {
         if (w[i][j] > Wmax) Wmax = w[i][j];
      }
   }
   
   CW = 0;
   if (rodaSaramaki == SIM)
   {
      for (i=0; i<linmax; i++)
      {
         Ki=0;
         for (j=0; j<linmax; j++)
         {
            if (w[i][j] > 0) ++Ki;
         }
         F=0;
         for (j=0; j<linmax; j++)
         {
            for (k=0; k<linmax; k++)
            {
               if ((i != j) && (i != k) && (j != k))
               {
                  if ((w[i][j] > 0) && (w[i][k] > 0) && (w[j][k] > 0))
                  {
                     Wmij = (float)(w[i][j]/Wmax);
                     Wmjk = (float)(w[j][k]/Wmax);
                     Wmki = (float)(w[k][i]/Wmax);
                     Wraiz3 = (float)(Wmij*Wmjk*Wmki);
                     Wraiz3 = (float)pow((float)Wraiz3, (float)(float)1.0/(float)3.0);
                     F += Wraiz3;
                  }
               }
            }
         }
         temp = (float)(Ki - (float)1.0);
         if (temp == 0) temp = 1.0;
         Cw[i] = (float)((float)((float)1.0/(float)(Ki*temp)) * F); 
         CW += Cw[i];
      }
   }
   else
   {
      for (i=0; i<linmax; i++)
      {
         Ki=0;
         Si=0;
         for (j=0; j<linmax; j++)
         {         
            if (w[i][j] > 0) ++Ki;
            if (w[i][j] > 0) Si=(float)(Si + w[i][j]);
         }
         F=0;
         for (j=0; j<linmax; j++)
         {
            for (k=0; k<linmax; k++)
            {            
               if ((i != j) && (i != k) && (j != k))
               {
                  if ((w[i][j] > 0) && (w[i][k] > 0) && (w[j][k] > 0))
                  {
                     carga = (float)((float)(w[i][j] + w[i][k])/(float)2);
                     F = (float)(F + carga);
                  }
               }
            }
         }
         temp = (float)(Ki - (float)1.0);
         if (temp == 0) temp = 1.0;
         temp2=(float)((float)1.0/(float)(Si*temp));
         Cw[i] = (float)(temp2 * F);
         CW += Cw[i];
      }
   }
   CW = CW/linmax;

   CwS = 0;
   if ((CW + Wmedio) != 0)
   {
      CwS = (CW - Wmedio)/(CW + Wmedio);
   }

   if (fase2 == NAO)
   {
      CwSObservado = CwS;
      PCwS = 0;
   }

   if (fase2 == SIM)
   {
      WmedioSimulacao[rodada-1] = Wmedio;
      ESimulacao[rodada-1] = E;
      CWSimulacao[rodada-1] = CW;
      CwSSimulacao[rodada-1] = CwS;
      goto RepeteSimulacao;
   }

   if (imprimeCw)
   {
      nome[0] = 'C';
      nome[1] = 'w';
      nome[2] = '[';
      for (i=0; i<(MAX_PATH-1); i++)
      {
         nome[3+i] = nomeA[i];
         if (nomeA[i] == '.')
            break;
      }
      nome[3+i] = ']';
      nome[4+i] = '.';
      nome[5+i] = 't';
      nome[6+i] = 'x';
      nome[7+i] = 't';
      nome[8+i] = '\0';
      stream=fopen(nome,"at");
      for (i=0; i<linmax; i++)
      {
         fprintf(stream, "Cw[%d]=%06.4f\n", i, Cw[i]);
      }
      fclose(stream);   
   }

   if (imprimeE)
   {
      nome[0] = 'E';
      nome[1] = '[';
      for (i=0; i<(MAX_PATH-1); i++)
      {
         nome[2+i] = nomeA[i];
         if (nomeA[i] == '.')
            break;
      }
      nome[2+i] = ']';
      nome[3+i] = '.';
      nome[4+i] = 't';
      nome[5+i] = 'x';
      nome[6+i] = 't';
      nome[7+i] = '\0';
      stream=fopen(nome,"at");
      for (i=0; i<linmax; i++)
      {
         fprintf(stream, "E[%d]=%06.4f\n", i, Ei[i]);
      }
      fclose(stream);   
   }

   // Passo 5. Calculando a matriz binária
   for (i=0; i<linmax; i++)
   {
      for (j=0; j<linmax; j++)
      {
         if (w[i][j] >= Wmedio)
         {
            B[i][j] = 1;
         }
         else
         {
            B[i][j] = 0;
         }
      }
   }

   // Passo 6. Imprimir Indices[matriz].txt
   for (i=0; i<(MAX_PATH-1); i++)
   {
      nomeX[i+8] = nomeA[i];
      if (nomeA[i] == '.')
         break;
   }
   j=i+8;
   nomeX[j++] = ']';
   nomeX[j++] = '.';
   nomeX[j++] = 't';
   nomeX[j++] = 'x';
   nomeX[j++] = 't';
   nomeX[j++]   = '\0';

   stream=fopen(nomeX,"at");
// fprintf(stream, "<W>=%06.4f E=%06.4f Cw=%06.4f\n",Wmedio, E, CW);  
   // foi mudado para imprimir somente E e CwS
   fprintf(stream, "E=%06.4f CwS=%06.4f\n", E, CwS);
   fclose(stream);   

   for (i=0; i<MAX_PATH; i++)
   {
      nome[i] = nomeA[i];
      if (nomeA[i] == '.') break;
   }
   nome[++i] = 'm';
   nome[++i] = 'a';
   nome[++i] = 't';
   nome[++i] = '\0';

   lixo[0] = '"';
   lixo[1] = '\0';
   stream=fopen(nome,"at");
   fprintf(stream, "%s%d\n","*Vertices       ", linmax);
   for (i=0; i<linmax; i++)
   {
      fprintf(stream, "       %d %s%s%s\n", (i+1), lixo, NomesLinhas[i], lixo);
   }
   fprintf(stream, "%s\n","*Matrix");
   
   for (i=0; i<linmax; i++)
   {
      for (j=0; j<linmax; j++)
      {
         fprintf(stream, "          %7.03f", (float)w[i][j]*weightFactor);
      }
      fprintf(stream, "\n");
   }   
   fclose(stream);   
   

   // Imprime B[Matriz].mat
   nome[0] = 'B';
   nome[1] = '[';
   for (i=0; i<(MAX_PATH-1); i++)
   {
      nome[2+i] = nomeA[i];
      if (nomeA[i] == '.') break;
   }
   nome[2+i] = ']';
   nome[3+i] = '.';
   nome[4+i] = 'm';
   nome[5+i] = 'a';
   nome[6+i] = 't';
   nome[7+i] = '\0';

   lixo[0] = '"';
   lixo[1] = '\0';
   stream=fopen(nome,"at");
   fprintf(stream, "%s%d\n","*Vertices       ", linmax);
   for (i=0; i<linmax; i++)
   {
      fprintf(stream, "       %d %s%s%s\n", (i+1), lixo, NomesLinhas[i], lixo);
   }
   fprintf(stream, "%s\n","*Matrix");
   
   for (i=0; i<linmax; i++)
   {
      for (j=0; j<linmax; j++)
      {
         fprintf(stream, "          %d", B[i][j]);
      }
      fprintf(stream, "\n");
   }   
   fclose(stream);   


   // Passo 7. Calcular var(E)
   if (fase2 == NAO)
   {
      linhaRetirada = 0;
      calculandoVarE = SIM;
      while (linhaRetirada < linmaxOriginal)
      {
         for(i=0; i<linhaRetirada; i++)
         {
            for(j=0; j<colmaxOriginal; j++)
            {
               R[i][j] = Matriz[i][j];
            }
         }

         for(i=linhaRetirada+1; i<linmaxOriginal; i++)
         {
            for(j=0; j<colmaxOriginal; j++)
            {
               R[i-1][j] = Matriz[i][j];
            }
         }

         // Calcula Ei
         linmax = linmaxOriginal - 1;

         goto RepeteCalculoFase2;

retiraProximaLinha:
         ++linhaRetirada;
      }
      
      calculandoVarE = NAO;
      linmax = linmaxOriginal;

      Emedio = 0;
      for (i=0; i<linmax; i++)
      {
         Emedio += EiVarE[i];
      }
      Emedio = Emedio/(float)linmax;
      Emedio = Emedio * Emedio;

      for (i=0; i<linmax; i++)
      {
         EiVarE[i] = EiVarE[i] * EiVarE[i];
      }

      varE = 0;
      for (i=0; i<linmax; i++)
      {
         varE += EiVarE[i];
      }
      varE = varE  - ((float)linmax * Emedio);
      
      if (linmax > 1)
      {
         varE = varE / ((float)(linmax*(linmax - 1)));
      }
      else
      {
         varE = 0;
      }

      // imprime varE VarE[matriz].txt
      for (i=0; i<(MAX_PATH-1); i++)
      {
         nomeV[i+5] = nomeA[i];
         if (nomeA[i] == '.')
            break;
      }
      j=i+5;
      nomeV[j++] = ']';
      nomeV[j++] = '.';
      nomeV[j++] = 't';
      nomeV[j++] = 'x';
      nomeV[j++] = 't';
      nomeV[j++] = '\0';
      stream=fopen(nomeV,"at");
      fprintf(stream, "Var(E)=%012.10f\n", varE);
      fclose(stream);   

   }

RepeteSimulacao:   
   // Passo 8. Null model analyses (Apenas para a opção 2 do "tipo de dado")
   if ((TipoDadoMatriz == 2) && (simulacao))
   {
      // Zera a matriz R[i][j]
      for (i=0; i<linmax; i++)
      {
         for (j=0; j<colmax; j++)
         {
            R[i][j]=0;
         }
      }

      // Soma linhas
      SomaMatriz=0;
      for(i=0; i<linmax; i++)
      {
         Somai[i]=0;
         for (j=0; j<colmax; j++)
         {
            Somai[i] += Matriz[i][j];
         }
         SomaMatriz += Somai[i];
      }

      base_probabilidade=0;
      // Verifica o tipo de fórmula
      if (formula == 1)
      {
         // Numérica
         // Soma colunas usando a Matriz Original
         for (j=0; j<colmax; j++)
         {
            Somaj[j]=0;
            for(i=0; i<linmax; i++)
            {
               Somaj[j] += Matriz[i][j];
            }
            Probabilidadej[j] = (unsigned int)((Somaj[j]/SomaMatriz)*32000);
            base_probabilidade += Probabilidadej[j];
         }
      }
      else
      {
         // Proporções Médias
         // Soma colunas usando a Matriz de Proporções
         for (j=0; j<colmax; j++)
         {
            Somaj[j]=0;
            for(i=0; i<linmax; i++)
            {
               Somaj[j] += POriginal[i][j];
            }
            Probabilidadej[j] = (unsigned int)((Somaj[j]/linmax)*32000);
            base_probabilidade += Probabilidadej[j];
         }
      }

      ProbIni[0] = 0;
      ProbFim[0] = Probabilidadej[0];
      for (i=1; i<colmax; i++)
      {
         ProbIni[i] = ProbFim[i-1] + 1;
         ProbFim[i] = ProbIni[i] + Probabilidadej[i];
      }

      for (i=0; i<linmax; i++)
      {
SorteiaMaisUma:
         aleatorio = rand(); /* numeros ate RAND_MAX (0x7fff) */
         alatorio_float = (float)(aleatorio)/(float)RAND_MAX;
         ProbabilidadeSorteada = (unsigned int)(alatorio_float*base_probabilidade);
         for (j=0; j<colmax; j++)
         {
            if ((ProbabilidadeSorteada >= ProbIni[j]) && (ProbabilidadeSorteada <= ProbFim[j]))
            {
               R[i][j] += 1;
               break;
            }
         }

         SomaR=0;
         for (j=0; j<colmax; j++)
         {
            SomaR += R[i][j];
         }
         if (SomaR < Somai[i]) goto SorteiaMaisUma;
      }

      fase2=SIM;
      if (++rodada <= replicas) goto RepeteCalculoFase2;

      mediaCwS=0;
      for (i=0; i<replicas; i++)
      {
         mediaCwS += CwSSimulacao[i];
      }
      mediaCwS = (float)(mediaCwS/(float)replicas);

      for (i=0; i<replicas; i++)
      {
         if (rodaSaramaki == SIM)
         {
            // Usar nova fórmula de modulo
            //if (moduloS(CwSSimulacao[i]) >= moduloS(CwSObservado)) ++PCwS;
            if (moduloB(CwSSimulacao[i]) >= moduloS(CwSObservado)) ++PCwS;
         }
         else
         {
            if (moduloB(CwSSimulacao[i]) >= moduloB(CwSObservado)) ++PCwS;
         }
      }

      nome[0] = 'b';
      nome[1] = 'o';
      nome[2] = 'o';
      nome[3] = 't';
      nome[4] = '[';

      for (i=0; i<(MAX_PATH-5); i++)
      {
         nome[5+i] = nomeA[i];
         if (nomeA[i] == '.')
            break;
      }
      
      nome[5+i] = ']';
      nome[6+i] = '.';
      nome[7+i] = 't';
      nome[8+i] = 'x';
      nome[9+i] = 't';
      nome[10+i] = '\0';

      if (imprimeBoot)
      {
         stream=fopen(nome,"at");
         for (i=0; i<replicas; i++)
         {
            if (i == 0) fprintf(stream, "<O>\tE\tCw\tCws\n");
            fprintf(stream, "%06.4f\t%06.4f\t%06.4f\t%06.4f\n",WmedioSimulacao[i], ESimulacao[i], CWSimulacao[i], CwSSimulacao[i]);
         }
         fclose(stream);   
      }

      nome[0] = 'P';
      nome[1] = '[';
      for (i=0; i<(MAX_PATH-1); i++)
      {
         nome[2+i] = nomeA[i];
         if (nomeA[i] == '.')
            break;
      }
      nome[2+i] = ']';
      nome[3+i] = '.';
      nome[4+i] = 't';
      nome[5+i] = 'x';
      nome[6+i] = 't';
      nome[7+i] = '\0';
      stream=fopen(nome,"at");
      if (replicas == 0) replicas=1;
      fprintf(stream, "P(E)=%06.4f P(CwS)=%06.4f\n", 
         (float)((float)PE/(float)replicas), 
         (float)((float)PCwS/(float)replicas));
      fclose(stream);   
     
   }

   ++indice_arquivos;
   goto Processa_Arquivo;


}   


float moduloS(float x)
{
   if(x < 0) return (-x);
   else return (x);
}

float moduloB(float x)
{
   if(x < mediaCwS) return (mediaCwS-x);
   else return (x-mediaCwS);
}