//  Simulacao de Monte Carlo do modelo SIS (linguagem C)

// VOCE DEVE RENOMEAR O ARQUIVO ANTES DE USAR, MUDE PARA: sis_campo_medio.c


// INSTRUCOES PARA USUARIOS DE LINUX:

// COMPILACAO: gcc sis_campo_medio.c -o prog -lm

// ISSO VAI GERAR UM ARQUIVO EXECUTAVEL COM O NOME: prog 

// DEPOIS BASTA FAZER ./prog



// Suscetiveis (S): representados por +1

// Infectados (I): representados por -1



#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define N 10000 // tamanho da populacao

#define tmax 10000 // tempo total da simulacao

#define density 0.98 // densidade inicial de Suscetiveis

#define alpha 0.1 // probabilidade da transicao I-->S

#define lambda 0.15 // probabilidade da transicao S-->I


#define nsamples 10  // numero de simulacoes independentes ("ensembles")


// gerador de numero aleatorios:
#define rfactor (long int)302875106592253 // numero "magico"
#define rmaxint (double)18446744073709551616.0 //maior inteiro p/ 64 bits


void main()
{
  
  int state[N],state_temp[N];
  int n1,cont,cont2,realizations,k1;
  int S_cont,I_cont,node1,node2,count;
  int i,j,l,t,t1,t_discar,sample;
  double factor,S_ave,I_ave;
  double Prob,num,num2,num3,denom,numer,prob;
  double S[tmax],I[tmax];
  FILE *fp,*fq;
  char out_name[100],out_name2[100];

  unsigned r0,e;

  long int r;

  
  // abrindo o arquivo de dados:
  sprintf(out_name,"S_I_S_x_time_mean_field_N%d_lambda%4.3lf_alpha%3.2lf.dat",N,lambda,alpha); 
  fp=fopen(out_name,"w");


  r0=3412341; //semente do gerador de numeros aleatorios

  r=r0;
  for(k1=0;k1<100;k1++) r*=rfactor; 
  

  for(t=0;t<tmax;t++)
    {
      S[t]=0.0;
      I[t]=0.0;
    }

  S_ave=0.0;
  I_ave=0.0;

  
  for(sample=0;sample<nsamples;sample++) // laco de simulacoes diferentes
    {
      
      // zerando os contadores:
      
      S_cont=0;
      I_cont=0;
      
      
      for(i=0;i<N;i++)
	{
	  // configuracao inicial de individuos:
	  
	  r*=rfactor;
	  num=((unsigned long)r)/rmaxint;	
	  
	  if(num<0.98)  // Suscetiveis em t=0
	    {
	      state[i]=1;
	      state_temp[i]=1;
	      S_cont++;
	    }
	  else // Infectados em t=0
	    {
	      state[i]=-1;
	      state_temp[i]=-1;
	      I_cont++;
	    }
	}
      
      t1=0;
      
      S[t1]+=((double) S_cont);
      
      I[t1]+=((double) I_cont);



      //################################################

      //                 Dinamica

      //################################################
      
      
      
      for(t=1;t<tmax;t++) // laco de tempo
	{	  
	  for(i=0;i<N;i++) // percorrendo a populacao em ordem
	    {

	      if(state[i]==-1) // Se o agente "i" esta' Infectado
		{
		  r*=rfactor;
		  num=((unsigned long)r)/rmaxint;
		  
		  if(num<alpha) 
		    {
		      state_temp[i]=1; // "i" se torna Suscetivel com
		                       // probabilidade "alpha"
		    }	
		}
	      
	      else if(state[i]==1) // Se o agente "i" esta' Susceptivel
		{
		  do
		    {
		      // escolhemos um "vizinho" de "i" aleatoriamente
		      r*=rfactor;
		      num2=((unsigned long)r)/rmaxint;
		      
		      node2=((int)N*num2);
		      
		    }while(node2==i); // garantir que o vizinho nao e' "i"

		  if(state[node2]==-1) // Se o vizinho esta' Infectado
		    {
		      r*=rfactor;
		      num3=((unsigned long)r)/rmaxint;
		      
		      if(num3<lambda) 
			{
			  state_temp[i]=-1; // "i" se torna Infectado com
			                    // probabilida "lambda"
			}		 
		    } 
		}
		  
	    } // fim do laco para os individuos


	  // atualizacao paralela dos estados apos 1 passo de tempo:
	  for(j=0;j<N;j++) state[j]=state_temp[j];

	  // zerando os contadores das populacoes S e I
	  S_cont=0;
	  I_cont=0;
	  
	  // contando os agentes nos estados S e I
	  for(i=0;i<N;i++)
	    {
	      if(state[i]==1) S_cont++;
	      else if(state[i]==-1) I_cont++;
	    }


	  // acumulando o numero de agentes em cada estado no instante t:
	  S[t]+=((double) S_cont);
	  
	  I[t]+=((double) I_cont);
	  

	  	  
	} // fim do laco de tempo

            
      
    } //fim do laco de simulacoes diferentes

  

  
  // imprimindo os dados de S x t e I x t no arquivo:
  // 1a coluna: t
  // 2a coluna: densidade de agentes S
  // 3a coluna: densidade de agentes I
  
  for(t=0;t<tmax;t++)
    {
      S[t]/=N*nsamples;
      I[t]/=N*nsamples;
      fprintf(fp,"%d\t%12.10lf\t%12.10lf\n",t,S[t],I[t]);
    }

  fclose(fp);

  
}
