//          - Source Treminal Network Unrealiability -
 //      Permutation Monte Carlo over a Static Network Model
 //  Leslie Murray - FCEIA, Universidad Nacional de Rosario, Argentina
 //                      - July, 2023 -

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

#define U genrand_real1()
#define semilla(s) init_genrand(s)


#define expo(lmbd) (-log(U)/lmbd)
// Generador de numeros aleatorios con distribucion exponencial de tasa "lmbd"


typedef struct link{
  double rlb;       // link reliability
  double lambda;    // rate to sample the link from exponential ditribution
  double tau;       // time when the link goes from failed to operative
  int op;           // 1 if the link is operative and 0 if it remains failed
  int node1;        // structure link models connection between node1 and node2
  int node2;        // (in both directions) i.e., from 1 to 2 and from 2 to 1
}lnk, *pt_lnk;

int s;              // network source node
int t;              // network target node
int numnodes;       // network number of nodes
int numlinks;       // network number of links

int **list;         // matrix to hold the adjacency lists

lnk *E;             // array to hold the set of all links such that: link from
                    // x to y and from y to x share the same bucket. Links are
                    // numbered from 0 to numlinks-1

                    // Nodes are numbered from 1 to numnodes (just like in the
                    // text files where the network is loaded from). In the
                    // data structures, node 0 is reserved for some algorithms
                    // tricks

void reset();       // sets tau=0 for all the links (among other things)

int next();         // returns the number of the next link to become operative

int *visited;       // array of nodes such that 1 indicates that BFS has
                    // visited the node and 0 that it hasn't

int connected;      // 1 if there's a path of operative links connecting s and
                    // t and 0 otherwise

double *lim;        // array to set limits to sort from discrete distribution

double previous_tau;// the latest time a node became operative

void Initialize(char* filename);
                    // takes the text file with the network info as input,
                    // builds up the rest of the structures and initializes all
                    // the variables

int Phi();          // returns 1 if, for the set of links that are already
                    // operative, the network is operative, and 0 otherwise

void DFS(int node); // performs Depth First Search from "node" and, as soon as
                    // node "t" is visited sets connected=1 and stops

double tot;         // the sum of the rates (lambdas) of all the nodes that are
                    // still failed

void accept(int n); // to be executed when the time just sorted (tau of link n)
                    // is accepted

void set_tau(int l);
                    // Sort the time to become operative to link l

double *D;          // array to hold the rates of the exponential random
                    // variables of the times between commutations

int m;              // index to save the number of occupied position of D


void set_perm();    // sets a links permutation and updates array D


double prob_perm(); // Estimates the probability that the permutation
                    // built by set_perm() puts the network operative
                    // in time t>1


int seed;           // seed for the random numbers generator

/********************************************************************
* The executable "e" runs as:  ./e <File> <Size> <Seed>, where:
* File: is the text file with the network (graph) info
* Size: is the number of trajectories launched from 0
* Seed: is the seed for the random numbers generator
********************************************************************/
int main(int argc, char *argv[]){
  int i,j,defined,size;
  double S,S2,X,V,Q,t;
  clock_t  ti;

  if(argc<4){
    printf("\n Some input data is missing! Run as:\n");
    printf("\n ./executable <archivo> <N_mc> <semilla>\n\n");
    exit(1);
  }
  if(argc>4){
    printf("\n You've entered more data than necessary! Run as:\n");
    printf("\n ./executable <archivo> <N_mc> <semilla>\n\n");
    exit(1);
  }
  if((size=atoi(argv[2]))<1){
    printf("\n The number of replications can not be less than 1! Run as:\n");
    printf("\n ./executable <archivo> <N_mc> <semilla>\n\n");
    exit(1);
  }

  if((seed=atoi(argv[3]))<0){
    printf("\n The seed can not be negative! Run as:\n");
    printf("\n ./executable <archivo> <N_mc> <semilla>\n\n");
    exit(1);
  }

  Initialize(argv[1]);


  //------------- THE PERMUTATION MONTE CARLO ALGORITHM -------------
  S=0.0;
  S2=0.0;

  ti = clock();

  for(i=0;i<size;i++){
    set_perm();
    Q=prob_perm();

    S=S+Q;
    S2=S2+Q*Q;

  }

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


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

  //------------------ PRINT OF OUTPUT AND RESULTS ------------------
  printf("\n      PERMUTATION MONTE CARLO     seed=%d",seed);
  printf("\n      Network: %s   Replications: %d   ExecTime=%f",argv[1],size,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)",sqrt(V)/Q,100*sqrt(V)/Q);
  printf("\n      --------------------------------------------------------------\n\n");
  //-----------------------------------------------------------------

  return 0;

}


/********************************************************************
* Initialize reads from the file that holds the network data all the
* information necessary to build up and allocate the data structures
* It also initializes the random numbers generator.
*********************************************************************/
void Initialize(char* filename){

  FILE* fp;
  int i,j,k,l,nodex,nodey,num,skip;
  double reliability;

  // Open the file with the network info and scan (at the top of it) for:
  // - source and target nodes
  // - number of nodes
  // - number of links
  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);

/******************************************************************************/
//ESTA MODIFICACION SE DEBE A QUE ESTE PROGRAMA FUE ORIGINALMENTE CONCEBIDO PARA
//ARCHIVOS DE REDES EN LAS QUE LOS NODOS SE NUMERAN DESDE 1 HASTA numnodes, PERO
//AQUI SE ESTAN USANDO ARCHIVOS CON LOS NODOS NUMERADOS DESDE 0 HASTA numnodes-1
  s++;
  t++;
/******************************************************************************/

  // Allocate array E and load it from the network file
  if((E=(lnk*)calloc(numlinks,sizeof(lnk)))==NULL){
    printf("\n Fail attempting to allocate memory...\n");
    exit(1);
  }

  // Clear array E before load it
  for(i=0;i<numlinks;i++){
    E[i].rlb=0.0;
    E[i].lambda=0.0;
    E[i].tau=0.0;
    E[i].op=0;
    E[i].node1=0;
    E[i].node2=0;
  }

  // Now load array E from the file
  l=0;
  for(i=1;i<=numnodes;i++){
    fscanf(fp, "%d%d", &nodex, &num);

/******************************************************************************/
//ESTA MODIFICACION SE DEBE A QUE ESTE PROGRAMA FUE ORIGINALMENTE CONCEBIDO PARA
//ARCHIVOS DE REDES EN LAS QUE LOS NODOS SE NUMERAN DESDE 1 HASTA numnodes, PERO
//AQUI SE ESTAN USANDO ARCHIVOS CON LOS NODOS NUMERADOS DESDE 0 HASTA numnodes-1
    nodex++;
/******************************************************************************/

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

/******************************************************************************/
//ESTA MODIFICACION SE DEBE A QUE ESTE PROGRAMA FUE ORIGINALMENTE CONCEBIDO PARA
//ARCHIVOS DE REDES EN LAS QUE LOS NODOS SE NUMERAN DESDE 1 HASTA numnodes, PERO
//AQUI SE ESTAN USANDO ARCHIVOS CON LOS NODOS NUMERADOS DESDE 0 HASTA numnodes-1
      nodey++;
/******************************************************************************/

      skip=0;
      for(k=0;k<numlinks;k++){
        if(E[k].node1==nodey && E[k].node2==nodex){
          skip=1;
          break;
        }
      }
      if(!skip){
        E[l].node1=nodex;
        E[l].node2=nodey;
        E[l].rlb=reliability;
        E[l].lambda=-log(1.0-reliability);
        l++;
      }
    }
  }

  if(numlinks!=l){
    printf("\n                    WARNING! \n");
    printf("\n%d links have been loaded and the\n",l);
    printf("\n  network is supposed to have %d \n\n",numlinks);
    exit(1);
  }

  fclose(fp);

  // Select a seed value and initialize the pseudorandom numbers generator
  semilla(seed);


  // Allocate matrix lists of size (numnodes+1)*(numnodes+1)
  // so as not to use "0" for the nodes
  if((list=(int**)calloc(numnodes+1,sizeof(int*)))==NULL){
    printf("\n Fail attempting to allocate memory...\n");
    exit(1);
  }
  for(i=0;i<=numnodes;i++){
    if((list[i]=(int*)calloc(numnodes+1,sizeof(int)))==NULL){
      printf("\n Fail attempting to allocate memory...\n");
      exit(1);
    }
  }

  // Initialize matrix lists, filling it with "0"
  for(i=1;i<=numnodes;i++)
    for(j=1;j<=numnodes;j++){
      list[i][j]=0;
    }


  // Allocate the visited array and set 0 for all nodes (not visited)
  if((visited=(int*)calloc(numnodes+1,sizeof(int)))==NULL){
    printf("\n Fail attempting to allocate memory...\n");
    exit(1);
  }

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

  // Allocate the lim array to aid sorting from a discrete distribution
  if((lim=(double*)calloc(numlinks,sizeof(double)))==NULL){
    printf("\n Fail attempting to allocate memory...\n");
    exit(1);
  }

  tot=0.0;
  for(i=0;i<numlinks;i++)
    tot=tot+E[i].lambda;

  lim[numlinks-1]=tot;
  for(i=numlinks-2;i>=0;i--)
   lim[i]=lim[i+1]-E[i+1].lambda;

  // Allocate array D
  if((D=(double*)calloc(numlinks,sizeof(double)))==NULL){
    printf("\n Fail attempting to allocate memory...\n");
    exit(1);
  }

  for(i=0;i<numlinks;i++)
    D[i]=0.0;

  m=0;

  connected=0;
  previous_tau=0.0;

  return;
}


/********************************************************************
* returns 1 if, for the set of links that are already operative, the
* network is operative, and 0 otherwise
*********************************************************************/
int Phi(){

  int i,j,k,n1,n2;

  // Clear the prior Adjcency Lists
  for(i=1;i<=numnodes;i++){
    k=1;
    while(list[i][k]>0){
      list[i][k++]=0;
    }
  }

  // Build up the adjacency lists for all the sampled links
  for(i=0;i<numlinks;i++){
    if(E[i].op){
      n1=E[i].node1;
      n2=E[i].node2;
      k=1;
      while(list[n1][k]>0){
        k++;
      }
      list[n1][k]=n2;
      k=1;
      while(list[n2][k]>0){
        k++;
      }
      list[n2][k]=n1;
    }
  }

  connected=0;

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

  DFS(s);

  return connected;
}


/********************************************************************
* makes a Depth First Search from node "s" and:
* - set connected=1 if node "t" is reached
* - set connected=0 if node "t" is not reached
********************************************************************/
void DFS(int node){

  int k;

  if(connected)
    return;

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

  visited[node]=1;

  k=1;
  while(list[node][k]>0){
    if(!visited[list[node][k]])
      DFS(list[node][k]);
    k++;
  }

  return;
}


/********************************************************************
* clears tau and op for matrix E, and also previous_tau
********************************************************************/
void reset(){

  int i;

  tot=0.0;
  for(i=0;i<numlinks;i++)
    tot=tot+E[i].lambda;

  lim[numlinks-1]=tot;
  for(i=numlinks-2;i>=0;i--)
   lim[i]=lim[i+1]-E[i+1].lambda;

  m=0;
  for(i=0;i<numlinks;i++){
    E[i].tau=0.0;
    E[i].op=0;

    D[i]=0.0;
  }

  previous_tau=0.0;

  return;

}


/********************************************************************
* It randomly selects a link to become operative and sets the
* exponentially commutation time to it
********************************************************************/
int next(){

  int i,j,k;
  double rnd;

  rnd=U*tot;

//  k=0;
//  while(lim[k]<=rnd)
//    k++;

  for(k=0;k<numlinks;k++)
    if(lim[k]>rnd)
      break;

  if(k>=numlinks)              // there are no more links left to set tau to
    return -1;

  E[k].tau=previous_tau+expo(tot);

  return k;

}


/********************************************************************
* It accepts the link sampled by next() and updates some values
********************************************************************/
void accept(int n){

  int i;

  E[n].op=1;

  previous_tau=E[n].tau;

  tot=tot-E[n].lambda;

  lim[n]=-1.0;

  for(i=n+1;i<numlinks;i++){
    if(E[i].op==0)
      lim[i]=lim[i]-E[n].lambda;
  }

  return;
}


/********************************************************************
* Sets the corresponding commutation time to link "l"
********************************************************************/
void set_tau(int l){

  E[l].tau=previous_tau+expo(tot);

  return;

}


/********************************************************************
* Sets a links permutation and updates array D
********************************************************************/
void set_perm(){
  int i,j;

  reset();

  for(i=0;i<numlinks;i++){
    j=next();
    D[m]=tot;
    m++;
    accept(j);
    if(Phi())
      break;
  }

}



/********************************************************************
* Estimates the probability that the permutation built by set_perm()
* puts the the network operative in time t>1
********************************************************************/
double prob_perm(){
  int i,j,k;
  double prod,prod_sum,sum;
//for(i=0;i<m;i++)printf("D[%d]=%f\n",i,D[i]);
  prod=1.0;
  for(i=0;i<m;i++){
    prod=prod*D[i];
  }

  sum=0.0;
  for(j=0;j<m;j++){

    prod_sum=1.0;
    for(k=0;k<m;k++){
      if(k!=j)
        prod_sum=prod_sum*(D[k]-D[j]);
    }
    prod_sum=D[j]*prod_sum;

    sum=sum+exp(-D[j])/prod_sum;

  }
  sum=prod*sum;
//printf("%d   %f\n",m-1,sum);
  return sum;

}