//         Anticonfiabilidad sobre Modelo Fuente-Terminal (s-t)
//            F-Monte Carlo sobre Modelo Estatico de Redes
//   Leslie Muray - FCEIA, Universidad Nacional de Rosario, Argentina
//                         - Julio, 2023 -

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include "mt.c"


#define G genrand_int32()
#define semilla(s) init_genrand(s)


int s;
// s: nodo fuente u origen.

int t;
// t: nodo terminal o destino.

int numnodes;
// numnodes: numero de nodos.

int numlinks;
// numlinks: numero de enlaces.

int **adj;
// Es la matriz de adyacencia de dimension [numnodes][numnodes]. Un "0" en la
// posicion [i][j] significa que no hay enlace entre los nodos "i" y "j",
// mientras que un "1" en la posicion [i][j] significa que hay enlace entre los
// nodos "i" y "j"

double **rlb;
// Lo mismo que "adj" solo que en la posición de cada "1" en la matriz "adj", la
// matriz "rlb" tiene el valor de confiabilidad del enlace correspondiente.

int **smp;
// "smp" lleva un registro actualizado de los enlaces sorteados. Es una mtatriz
// de dimension [numnodes][numnodes]. Se inicia con ceros y cada vez que se
// sortea un enlace entre los nodos (i,j), se ponene en "1" los pares (i,j) y
// (j,i).

int **adj_f;
// Lo mismo que "adj" pero luego de una configuracion generada por Conf(k).

void clean();
// Completa con ceros "smp" y "adj_f" y los deja listos para una nueva falla.

int *visited;
// Es un arreglo de dimension [numnodes]. Si visited[n]=1 significa que el nodo
// "n" ha sido visitado por DFS mientras que visited[n]=0 indica que el nodo "n"
// no ha sido visitado por DFS.

int connected;
// connected=1 indica que los nodos "s" y "t" estan conectados, connected=0
// indica que "s" y "t" no estan conectados.

void Initialize(char* filename);
// filename es el nombre del archivo que contiene la topologia de la red.
// Initialize es la función que asigna memoria para todas las estructuras de
// dimensión variable e inicializa todas la variables.

int Phi();
// Es la funcion de Estructura. Devuelve 1 si los nodos "s" y "t" estan
// conectados y 0 en caso contrario.

void DFS(int node);
// Hace un Depth First Search a partir de "nodo". Se inicia con nodo=s y al
// final fija el valor de "connected" segun el recorrido haya alcanzado el nodo
// "t", o no.

void Conf(int k);
// Genera una configuracion fijando, con distribucion uniforme, k enlacces
// en 0 y los restantes en 1.

int *link;
// Arreglo de dimension [numlinks] para ser usado en Conf()

int *pos;
// Arreglo de dimension [numlinks] para ser usado en Conf()

double mass(int m, int j);
// Calcula la masa de probabilidad del conjunto de configuraciones con j enlaces
// sanos, siendo m el total de enlaces.

double pot(int j);
// Calcula la q elevado a la potencia j, con el mismo mecanismo que trabaja la
// función mass(), por lo que se evita el calculo directo y se hace la potencia 
// de forma iterativa, sin riesgos de overflows ni inestabilidades



double beta(int k);
// es el coeficiente beta, cuando hay k links en cero y numlink-k links en uno

int size;
// es el numero de pruebas de monte carlo (global) para evaluar la varianza

int Num;
// es el numero de pruebas de monte carlo por etapa, para las estimaciones de
// cada uno de los beta

double p;
// la confiabilidad de cada enlace (es para hacer mas simple las cuentas, ya que
// la red se carga de un archivo que permite diferentes valores por link)

double q;
// la anti-confiabilidad de cada enlace (es para hacer mas simple las cuentas,
// ya que la red se carga de un archivo que permite diferentes valores por link)

int seed; 
// Guarda el valor de la semilla del generador de numeros aleatorios.


/*******************************************************************************
* El ejecutable "e", resultante de la compilacion, se ejecuta como
* ./e <archivo> <pruebas> <semilla>, donde:
*   <archivo> nombre del archivo con los datos de la red
*   <pruebas> numero de prunas de Monte Carlo
*   <semilla> semilla del generador de numeros aleatorios
/******************************************************************************/
int main(int argc, char *argv[]){
  int i,j,k,l;
  double X,V,S,Q,t;
  clock_t  ti;

  if(argc<5){
    printf("\n ¡Faltan argumentos! Ejecutar como:\n");
    printf("\n ./ejecutable <archivo> <N_mc> <N_etapa> <semilla>\n\n");
    exit(1);
  } 
  if(argc>5){
    printf("\n ¡Se han insertado argumentos de mas! Ejecutar como:\n\n");
    printf("\n ./ejecutable <archivo> <N_mc> <N_etapa> <semilla>\n\n");
    exit(1);
  }
  if((size=atoi(argv[2]))<1){
    printf("\n ¡El numero de pruebas no puede ser < 1! Ejecutar como:\n");
    printf("\n ./ejecutable <archivo> <N_mc> <N_etapa> <semilla>\n\n");
    exit(1);
  }
  if((Num=atoi(argv[3]))<1){
    printf("\n ¡El numero de pruebas no puede ser < 1! Ejecutar como:\n");
    printf("\n ./ejecutable <archivo> <N_mc> <N_etapa> <semilla>\n\n");
    exit(1);
  }
  if((seed=atoi(argv[4]))<0){
    printf("\n ¡La semilla no puede ser negativa! Ejecutar como:\n");
    printf("\n ./ejecutable <archivo> <N_mc> <N_etapa> <semilla>\n\n");
    exit(1);
  }

  Initialize(argv[1]);



  //-------------------- EL ALGORITMO DE MONTE CARLO CRUDO ---------------------

  ti = clock();

  S=0;
  X=0.0;
  V=0.0;

  for(i=0;i<size;i++){

    X=0.0;
    for(j=1;j<numlinks;j++){
      X=X+beta(j)*mass(numlinks,j);
    }
      X=X+pot(numlinks);


    S+=X;
    V+=X*X;
  }

  Q=(double)S/size;
  V=(V/size-Q*Q)/(size-1);

  t=(double)(clock()-ti)/CLOCKS_PER_SEC;

  //---------------------------- SALIDA Y RESULTADOS ---------------------------
 
 printf("\n    F-MONTE CARLO     seed=%d",seed);
  printf("\n    Network: %s   Replications: %d   Trials p/stage: %d   ExecTime=%f",argv[1],size,Num,t);
  printf("\n\n     Q=%1.16f = %1.2e  (Unreliability)",Q,Q);
  printf("\n     V=%1.16f = %1.2e  (Variance)",V,V);
  printf("\n    SD=%1.16f = %1.2e  (Standard Deviation)",sqrt(V),sqrt(V));
  printf("\n    RE=%1.16f = %1.2f%%     (Relative Error)\n",sqrt(V)/Q,
  100*sqrt(V)/Q);
  printf("    ---------------------------------------------------------------");
  printf("\n\n");
  printf("\t\t V*t = %1.16f = %1.2e \n\n",V*t,V*t);


  free(adj);
  free(adj_f);
  free(rlb);
  free(smp);
  free(visited);
  free(link);
  free(pos);

  return 0;

}

/******************************************************************************/

void Initialize(char* filename){

  FILE* fp;
  int i,j,node1,node2,num;
  double reliability;


  // Se abre el archivo con la informacion sobre la red y se lo recorre para
  // leer los datos uno tas otro
  if((fp=fopen(filename,"r"))==NULL){
   printf("\n Fail attempting to open a disk file...\n");
    exit(1);
  }
  fscanf(fp,"%d",&s);
  fscanf(fp,"%d",&t);
  fscanf(fp,"%d",&numnodes);
  fscanf(fp,"%d",&numlinks);


  // Se asigna la matriz adj de [numnodes][numnodes]
  adj=(int**)calloc(numnodes,sizeof(int*));
  for(i=0;i<numnodes;i++)
    adj[i]=(int*)calloc(numnodes,sizeof(int));

  // Se inicializa a cero la matriz adj
  for(i=0;i<numnodes;i++)
    for(j=0;j<numnodes;j++){
      adj[i][j]=0;
    }


  // Se asigna la matriz adj_f de [numnodes][numnodes]
  adj_f=(int**)calloc(numnodes,sizeof(int*));
  for(i=0;i<numnodes;i++)
    adj_f[i]=(int*)calloc(numnodes,sizeof(int));

  // Se inicializa a cero la matriz adj_f
  for(i=0;i<numnodes;i++)
    for(j=0;j<numnodes;j++){
      adj_f[i][j]=0;
    }


  // Se asigna la matriz rlb de [numnodes][numnodes]
  rlb=(double**)calloc(numnodes,sizeof(double*));
  for(i=0;i<numnodes;i++)
    rlb[i]=(double*)calloc(numnodes,sizeof(double));

  // Se inicializa a cero la matriz rlb
  for(i=0;i<numnodes;i++)
    for(j=0;j<numnodes;j++){
      rlb[i][j]=0.0;
    }


  // Se asigna la matriz smp de [numnodes][numnodes]
  smp=(int**)calloc(numnodes,sizeof(int*));
  for(i=0;i<numnodes;i++)
    smp[i]=(int*)calloc(numnodes,sizeof(int));

  // Se inicializa a cero la matriz smp
  for(i=0;i<numnodes;i++)
    for(j=0;j<numnodes;j++){
      smp[i][j]=0;
    }


  // Allocate the vector links of size [numlinks]
  if((link=(int*)calloc(numlinks,sizeof(int)))==NULL){
    printf("\n Fail attempting to allocate memory...\n");
    exit(1);
  }

  // Allocate the vector links of size [numlinks]
  if((pos=(int*)calloc(numlinks,sizeof(int)))==NULL){
    printf("\n Fail attempting to allocate memory...\n");
    exit(1);
  }




  // Se recorre cada una de las lineas siguientes del archivo con la informacion
  // sobre la red y se van colocando los valores leidos en las variables
  // correspondientes
  for(i=0;i<numnodes;i++){
    fscanf(fp, "%d%d", &node1, &num);

    for(j=0;j<num;j++){
      fscanf(fp, "%d%lf", &node2, &reliability);

      adj[node1][node2]=1;
      rlb[node1][node2]=reliability;

      adj[node2][node1]=1;
      rlb[node2][node1]=reliability;

      p=reliability;
      q=1-p;
    }
  }

  fclose(fp);

  // Se inicializa el generador de numeros aleatorios con la semilla elegida
  semilla(seed);


  // Se asigna el arreglo visited
  visited=(int*)calloc(numnodes,sizeof(int));

  // Se inicializa el arreglo visited con ceros
  for(i=0;i<numnodes;i++)
    visited[i]=0;

  connected=0;

  return;

}

/******************************************************************************/

void Conf(int k){

  int i,j;

  clean();


  // a esta altura todos los enlaces estan en 1. aqui defino y guardo en pos[],
  // el numero de los k enlaces que voy a pasar a 0  (clean() los puso en 1).

  //for(i=0;i<numlinks;i++) link[i] = i;

  for(i=0; i<k;i++) {
    j = G % (numlinks-i);
    pos[i] = link[j];
    link[j] = link[numlinks-i-1];
  }


  for(i=0; i<numlinks;i++) {
    link[i]=1;
  }
  for(i=0; i<k; i++){
    link[pos[i]]=0;
  }


int x;
x=0;
  for(i=0;i<numnodes;i++){

    for(j=i+1;j<numnodes;j++){

      if(adj[i][j]){

        adj_f[i][j]=link[x];
        adj_f[j][i]=link[x];


        x++;
      }
    }
  }



  return;

}

/******************************************************************************/

int Phi(){

  int k;

  // Se resetean "connected" y "visited" antes de ejecutar DFS
  connected=0;

  for(k=0;k<numnodes;k++)
    visited[k]=0;


  DFS(s);

  return connected;
}

/******************************************************************************/

void DFS(int node){

  int i,j,k;

  if(connected)
    return;

  if(node == t){
    connected=1;
    return;
  }

  visited[node]=1;

  for(i=0;i<numnodes;i++){
    if(!visited[i] && adj_f[node][i]){
      DFS(i);
    }
  }

  return;
}

/******************************************************************************/

void clean(){

  int i,j,x;


  x=0;
  for(i=0;i<numnodes;i++){
    for(j=i+1;j<numnodes;j++){
      adj_f[i][j]=0;
      adj_f[j][i]=0;

      if(adj[i][j]){
        link[x]=x;
        x++;
      }
    }
  }


  return;

}

/******************************************************************************/

double beta(int k){

int i,S;
double Z;

   S=0;
   Z=0.0;

  for(i=0;i<Num;i++){
    Conf(k);
    Z=(1-Phi());
    S+=Z;
  }

  return (double)S/Num;

}


/******************************************************************************/

double mass(int m, int j){

  int i;
  double a;
  long double ce = 1.0;


  a=exp((j*log(q))/(m-j));


  for(i=m;i>=(j+1);i--){
    ce = ce * (i*p*a/(i-j));
  }


  return ce;

}


/******************************************************************************/

double pot(int j){

  int i;
  long double ce = 1.0;

  for(i=0;i<j;i++){
    ce = ce * q;
  }

  return ce;

}