/*
+-----------------------------------------------------------------------+
| BooleNet computes attractors in the synchronous Boolean network model |
| v1.1 - April 2005                                                     |
|                                                                       |
| Maxim Teslenko                                                        |
| maximt@kth.se                                                         |
| KTH, Stockholm, Sweden                                                |
+-----------------------------------------------------------------------+
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <assert.h>
#include <cudd.h>

#include <unistd.h>

#define MAXIMUM_NUMBER_OF_RELEVANT_NODES_COUNTED 100

#define PRINTF(X,Y)
//#define PRINTF(X,Y) printf(X,Y) 

#define PRINT(X) printf(#X":%d\n",X) 
//#define PRINT(X)

#define ASSERT_ON

#ifdef ASSERT_ON
#define ASSERT(x)           assert(x)  
#else
#define ASSERT(x)
#endif 

#define MAX_NETWORK_SIZE    10000
#define DEAD_BIT            ((unsigned)0x80000000)   /* dead edge */
#define KILL_EDGE(x)        (x = x | DEAD_BIT)       /* to mark edge dead */
#define IS_DEAD(x)          (((x) & DEAD_BIT) != 0 ) /* to check whether it is dead */
#define AI_BDD_NODES        500000

typedef unsigned Edge;           

typedef struct vertex {
  char                  f;                          /* associated Boolean funcion*/
  Edge                  left;                       /* left child                */
  Edge                  right;                      /* right child               */
  unsigned              out_count;                  /* # fanouts                 */
  Edge                 *out;                        /* fanout structure          */
  char                  type;                       /* 0 if redundant, 1 else    */
  char                  component;                  /* connected component */
  DdNode               *var;
  DdNode               *var_next;
} Vertex;


  /* Reading from file informaition*/
  int number_of_var,redunt_number_of_var,number_print_var;
  int *var2red_var;
  int *redunt_var2var;

#define MAX_LENGTH 10000
#define MAX_NUMBER_OF_INPUTS 20
   /*---------------------------------*/


//static Vertex         *vertex_array    = NULL;
//static Edge           *red_array       = NULL;
//static unsigned        number_vertices = 0;
static unsigned        size_nonred_array;

static int *var_index_array12_23,*var_index_array123_312,*var_index_array13_12;
static DdNode *cube_quantfy1,*cube_quantfy2,*cube_quantfy3;
static DdNode **var_array;

static DdNode **dd_array1;     /*Array to store bdd for every node of RBN*/
static DdNode **dd_array2;     /*Array to store bdd for every node of RBN*/
static DdNode **dd_array_step; /*Array to store one step bdd for every node of RBN*/

static void            RBN_init_global_var(DdManager *bdd_manage);
static void            RBN_quit_global_var(DdManager *bdd_manager);

DdNode*                ReadNET(char *fileName,DdManager *bdd_manager);
DdNode*                Next_State(DdManager *bdd_manager,DdNode* current_state);
/* --------------------------------------------------------------------------- */
int 
main(int argc, char *argv[]) 
{

  DdManager *bdd_manager;
  unsigned  number_of_atractors;
  DdNode   *covered_states;
  char      *temp_string;


    /* building BDD for the RBN */
    bdd_manager = Cudd_Init(0, 0, 1024, CUDD_CACHE_SLOTS*4, 0);
    //Cudd_DisableGarbageCollection(bdd_manager);
    //Cudd_AutodynDisable(bdd_manager);
    
    FILE *fp_bdd_out;
    if ((fp_bdd_out = fopen("temp.tmp", "w+")) == NULL) {
      fprintf(stderr, "Couldn't open file: temp.tmp\n");
      exit (1);
    }

    Cudd_SetStdout(bdd_manager, fp_bdd_out);

    number_of_atractors=0;
    //RBN_init_global_var(bdd_manager);

    /*----------------- Partitioned T^2n generaition --------------*/
    DdNode *bdd_node0 ,*bdd_node_atractor_states;
    DdNode *current_state,*next_state,*tmp;

    bdd_node0=ReadNET(argv[1],bdd_manager);

    //printf("T^2n contains %d node(s)\n", Cudd_DagSize(bdd_node0));
    bdd_node_atractor_states=Cudd_bddExistAbstract( bdd_manager , bdd_node0 , cube_quantfy1 );
    Cudd_Ref(bdd_node_atractor_states);
    Cudd_RecursiveDeref(bdd_manager, bdd_node0);

    //printf("Atractor bdd Contains %d node(s)\n\n", Cudd_DagSize(bdd_node_atractor_states));

    /*------- Move atractor states to original variable ----------*/
    tmp=Cudd_bddPermute(bdd_manager,bdd_node_atractor_states,var_index_array123_312);
    Cudd_Ref(tmp);
    Cudd_RecursiveDeref(bdd_manager, bdd_node_atractor_states );
    bdd_node_atractor_states=tmp;

    /*----------- Going through all states of all atractor  --------------*/
    covered_states=Cudd_Not(Cudd_ReadOne(bdd_manager));
    Cudd_Ref(covered_states);
    int atractor_counter=0,atractor_counter_size=0;
    temp_string = (char *)realloc(NULL, sizeof(char)*(number_of_var+2) );

    while(1){
      atractor_counter++;
      current_state=Cudd_bddPickOneMinterm(bdd_manager,bdd_node_atractor_states,var_array,size_nonred_array);
      Cudd_Ref(current_state);
      covered_states=current_state;
      Cudd_Ref(covered_states);

      atractor_counter_size=0;

      /*----------- Going through all states of one atractor  --------------*/
      while(1){
	atractor_counter_size++;
	Cudd_PrintMinterm(bdd_manager,current_state);
	fseek(fp_bdd_out,SEEK_SET,0);
	fgets(temp_string, number_print_var+2 , fp_bdd_out);
	fseek(fp_bdd_out,SEEK_SET,0);
	puts(&temp_string[1]);

	next_state=Next_State(bdd_manager,current_state);
	Cudd_Ref(next_state);
	Cudd_RecursiveDeref( bdd_manager , current_state );
	current_state=next_state;
	tmp = Cudd_bddOr( bdd_manager , current_state, covered_states );
	if( tmp == covered_states){
	  Cudd_RecursiveDeref( bdd_manager , current_state );
	  break;
	}
	Cudd_Ref(tmp);
	Cudd_RecursiveDeref(bdd_manager, covered_states);
	covered_states=tmp;
      }
      printf("Attractor %d is of length %d\n",atractor_counter,atractor_counter_size);
      //PRINT(atractor_counter_size);
      puts(" ");
      /*----------- End of going through all states of one atractor  --------------*/

      tmp = Cudd_bddXor(bdd_manager,bdd_node_atractor_states,covered_states);
      Cudd_Ref(tmp);
      Cudd_RecursiveDeref(bdd_manager, covered_states);
      Cudd_RecursiveDeref(bdd_manager, bdd_node_atractor_states);
      if( tmp == Cudd_Not(Cudd_ReadOne(bdd_manager)) ){
	Cudd_RecursiveDeref(bdd_manager, tmp );
	break;
      }
      bdd_node_atractor_states=tmp;
    }
    //PRINT(atractor_counter);
    printf("The number of attractors is %d\n",atractor_counter);


    RBN_quit_global_var(bdd_manager);

    fclose(fp_bdd_out);
    free(dd_array1);
    free(dd_array2);
    free(dd_array_step);
    free(temp_string);
 
    return 0;
}


/* NOTE: variable count start from node 1*/
static void
RBN_init_global_var(DdManager *bdd_manager)
{
  int i;
  DdNode *temp_node12,*temp_node23;

  var_index_array12_23=(int *)realloc(NULL,sizeof(int) *(3*size_nonred_array+1));

  var_index_array123_312=(int *)realloc(NULL,sizeof(int) *(3*size_nonred_array+1));

  var_index_array13_12=(int *)realloc(NULL,sizeof(int) *(3*size_nonred_array+1));

  var_array=(DdNode **)realloc(NULL,sizeof(DdNode *) *(size_nonred_array+1));

  for(i=1;i<=(size_nonred_array<<1);i++){
    var_index_array12_23[i]=i+size_nonred_array;
  }

  for(i=1;i<=size_nonred_array;i++){
    var_index_array123_312[i]=i+(size_nonred_array<<1);
  }
  for(i=1+size_nonred_array;i<=(size_nonred_array<<1);i++){
    var_index_array123_312[i+size_nonred_array]=i;
    var_index_array123_312[i]=i-size_nonred_array;
  }

  for(i=1;i<=size_nonred_array;i++){
    var_index_array13_12[i]=i;
  }
  for(i=1+size_nonred_array;i<=(size_nonred_array<<1);i++){
    var_index_array13_12[i+size_nonred_array]=i;
  }

  cube_quantfy1=Cudd_bddIthVar(bdd_manager,1);
  Cudd_Ref(cube_quantfy1);//One ref needed for variable
  Cudd_Ref(cube_quantfy1);//One ref will be deref later when cube_quantfy1 is calculated
  var_array[0]=cube_quantfy1;
  var_array[size_nonred_array]=cube_quantfy1;
  for(i=2;i<=size_nonred_array;i++){
    temp_node12=Cudd_bddIthVar(bdd_manager,i);
    Cudd_Ref(temp_node12);
    var_array[i-1]=temp_node12;
    temp_node23=Cudd_bddAnd(bdd_manager,temp_node12 ,cube_quantfy1);
    Cudd_Ref(temp_node23);
    Cudd_RecursiveDeref(bdd_manager,cube_quantfy1);
    cube_quantfy1=temp_node23;
  }

  cube_quantfy3=Cudd_bddIthVar(bdd_manager,(size_nonred_array<<1)+1);
  Cudd_Ref(cube_quantfy3);
  Cudd_Ref(cube_quantfy3);
  for(i=(size_nonred_array<<1)+2;i<=3*size_nonred_array;i++){
    temp_node12=Cudd_bddIthVar(bdd_manager,i);
    Cudd_Ref(temp_node12);
    temp_node23=Cudd_bddAnd(bdd_manager,temp_node12 ,cube_quantfy3);
    Cudd_Ref(temp_node23);
    Cudd_RecursiveDeref(bdd_manager,cube_quantfy3);
    cube_quantfy3=temp_node23;
  }

  cube_quantfy2=Cudd_bddIthVar(bdd_manager,size_nonred_array+1);
  Cudd_Ref(cube_quantfy2);
  Cudd_Ref(cube_quantfy2);
  for(i=size_nonred_array+2;i<=(size_nonred_array<<1);i++){
    temp_node12=Cudd_bddIthVar(bdd_manager,i);
    Cudd_Ref(temp_node12);
    temp_node23=Cudd_bddAnd(bdd_manager,temp_node12 ,cube_quantfy2);
    Cudd_Ref(temp_node23);
    Cudd_RecursiveDeref(bdd_manager,cube_quantfy2);
    cube_quantfy2=temp_node23;
  }
}

static void
RBN_quit_global_var(DdManager *bdd_manager)
{
  DdNode *temp_node;
  int i;

  Cudd_RecursiveDeref(bdd_manager,cube_quantfy1);
  Cudd_RecursiveDeref(bdd_manager,cube_quantfy2);
  Cudd_RecursiveDeref(bdd_manager,cube_quantfy3);

  for(i=1;i<=(3*size_nonred_array);i++){
    temp_node=Cudd_bddIthVar(bdd_manager,i);
    Cudd_RecursiveDeref(bdd_manager,temp_node); 
  }

  free(var_index_array13_12);
  free(var_index_array12_23);
  free(var_index_array123_312);
  free(var_array);
}



/* Read in NET file and returns 2^n transition function*/
DdNode* 
ReadNET(char *fileName,DdManager *bdd_manager)
{

  char temp[MAX_LENGTH];
  char slask[MAX_LENGTH];
  int i, j , k , tmp_1;

  int input_array[MAX_NUMBER_OF_INPUTS];
  int number_of_inputs,current_number_of_inputs=0;
  int current_var;
  FILE *fp;

  DdNode* tmr,*tmr2,*tmr3;
  DdNode** tmpp;

  

  /* Open NET file for reading */
  if ((fp = fopen(fileName, "r")) == NULL) {
    fprintf(stderr, "Couldn't open NET file: %s\n", fileName);
    exit (1);
  }

  /* Reading number of non-redundent var*/
  while(!feof(fp)) {
    fgets(temp,MAX_LENGTH, fp);
    if (strncmp(temp,".v",2)==0) {
      /* Find size of the cubes */
      sscanf(temp, "%s %i", slask, &number_of_var);
      break;
    }
  }
  var2red_var=(int*) calloc(number_of_var+1, sizeof(int));

  size_nonred_array=number_of_var;

  dd_array1=(DdNode **)calloc((size_nonred_array+1),sizeof(DdNode *));
  dd_array2=(DdNode **)malloc(sizeof(DdNode *) * (size_nonred_array+1));
  dd_array_step=(DdNode **)malloc(sizeof(DdNode *) * (size_nonred_array+1));


  redunt_number_of_var=number_of_var;
  redunt_var2var=(int*) calloc(redunt_number_of_var+1, sizeof(int));


  number_print_var=number_of_var;

  int count_nodes;
  for(count_nodes=1; count_nodes<=number_of_var; count_nodes++ ){
    redunt_var2var[count_nodes]=count_nodes;
    var2red_var[count_nodes]=count_nodes;
  }

 
  if(feof(fp)) {
    fprintf(stderr, "Wrong format of input file. End of file is reached but no node function is red" );
    exit(1);
  }

  
  RBN_init_global_var(bdd_manager);

  /*-------------------- Filling information about nodes  -------------------*/

  while(!feof(fp)){
    fgets(temp,MAX_LENGTH, fp);
  next_vertex:
    if (strncmp(temp,".n",2)==0) {
      /* Find size of the cubes */

      i=3;
      sscanf(&temp[i],"%d",&current_var);

      if( current_var < 0 || current_var >redunt_number_of_var){
	fprintf(stderr, "Wrong format of input file.The varible %d in string:\n %s exceeds number of declared variables.\n", current_var,temp);
	exit (1);
      }
      

      i++;
      /* go to first space */
      while(temp[i]!=' '){
	if(temp[i]=='\n'){
	  fprintf(stderr, "Wrong format of input file.Wrong format of a string: %s\n", temp );
	  exit (1);
	}
	i++;
      }
      i++;
      /* go to the end of sequense of spaces */
      while(temp[i]==' '){
	i++;
      }
      if(temp[i]=='\n'){
	fprintf(stderr, "Wrong format of input file.Wrong format of the string: %s\n", temp );
	exit (1);
      }
      sscanf(&temp[i],"%d",&number_of_inputs);

      if(number_of_inputs > MAX_NUMBER_OF_INPUTS){
	fprintf(stderr, "Wrong format of input file.Number of inputs exceeds allowed number of inputs %s\n", temp  );
	exit (1);
      }
      

      while(1){
	i++;
	while(temp[i]!=' '){
	  if(temp[i]=='\n'){
	    goto end_loops;
	  }
	  i++;
	}
	i++;
	/* go to the end of sequense of spaces */
	while(temp[i]==' '){
	  i++;
	}
	if(temp[i]=='\n'){
	  goto end_loops;
	}
	sscanf(&temp[i],"%d",&tmp_1);

	if( tmp_1 < 0 || tmp_1 >redunt_number_of_var){
	  fprintf(stderr,"Wrong format of input file.The varible %d in string:\n %s exceeds number of declared variables.\n", tmp_1,temp);
	  exit (1);
	}
	
	if( redunt_var2var[tmp_1] == 0){
	  fprintf(stderr, "Wrong format of input file.One of the input was not declared in list of inputs %s\n", temp  );
	  exit (1);
	}

	input_array[current_number_of_inputs++]= redunt_var2var[tmp_1];

      }
    end_loops:


      if(current_number_of_inputs!=number_of_inputs){
	fprintf(stderr, "Wrong format of input file.Declared number of inputs is not equal to actual number of input %s\n", temp  );
	exit (1);
      }


      /*---------------  Reading in the function of the node ---------------*/
      current_number_of_inputs=0;
      
      DdNode* tmp_cube;
      DdNode* tmp_cube2; 
      DdNode* tmp_all=Cudd_Not(Cudd_ReadOne(bdd_manager));
      Cudd_Ref( tmp_all );


      while(!feof(fp)){
	fgets(temp,MAX_LENGTH, fp);
	if( temp[0]=='0' || temp[0]=='1' || temp[0]=='-'){
	  
	  tmp_cube=Cudd_ReadOne(bdd_manager);
	  Cudd_Ref(tmp_cube);

	  for(i=0;i<number_of_inputs;i++){
	    if(temp[i]=='-')
	      continue;
	    if(temp[i]=='1'){
	      tmp_cube2=Cudd_bddAnd(bdd_manager,tmp_cube ,Cudd_bddIthVar(bdd_manager,input_array[i]));
	      Cudd_Ref( tmp_cube2 );
	      Cudd_RecursiveDeref( bdd_manager , tmp_cube );
	      tmp_cube = tmp_cube2;

	    }else{
	      if(temp[i]!='0'){
		fprintf(stderr, "Wrong format of input file.Wrong line %s\n", temp  );
		exit (1);	
	      }

	      tmp_cube2=Cudd_bddAnd(bdd_manager,tmp_cube ,Cudd_Not(Cudd_bddIthVar(bdd_manager,input_array[i])));
	      Cudd_Ref( tmp_cube2 );
	      Cudd_RecursiveDeref( bdd_manager , tmp_cube );
	      tmp_cube = tmp_cube2;
	      
	    }
	  }
	  
	  tmp_cube2=Cudd_bddOr(bdd_manager,tmp_cube ,tmp_all);
	  Cudd_Ref( tmp_cube2 );
	  Cudd_RecursiveDeref( bdd_manager , tmp_cube );
	  Cudd_RecursiveDeref( bdd_manager , tmp_all );
	  tmp_all=tmp_cube2;

	}else{

	  if( redunt_var2var[current_var] == 0){
	    fprintf(stderr, "Wrong format of input file.The node:%d was not declared \n", current_var  );
	    exit (1);
	  }
	  if( dd_array1[redunt_var2var[current_var]] != NULL){
	    fprintf(stderr, "Wrong format of input file.The node:%d was defined multiple times \n", current_var  );
	    exit (1);
	  }
	  dd_array1[redunt_var2var[current_var]]=tmp_all;
	  dd_array_step[redunt_var2var[current_var]]=tmp_all;
	  goto next_vertex;
	}	
      }

      if( redunt_var2var[current_var] == 0){
	fprintf(stderr, "Wrong format of input file.The node:%d was not declared \n", current_var  );
	exit (1);
      }
      if( dd_array1[redunt_var2var[current_var]] != NULL){
	fprintf(stderr, "Wrong format of input file.The node:%d was defined multiple times \n", current_var  );
	exit (1);
      }

      dd_array1[redunt_var2var[current_var]]=tmp_all;
      dd_array_step[redunt_var2var[current_var]]=tmp_all;

    }/* -------- end processing input informaition of a vertex ------------*/
    
  }
/*end processing the file*/

  for(i = 1; i <=size_nonred_array; i++) {
    if(dd_array1[i] == NULL){
      fprintf(stderr, "Wrong format of input file.The node:%d was not associated with function \n", var2red_var[i]  );
      exit (1);
    }
  }

  /*------------------- Generating T^2n function in partitioned way -------------- */
  for(i = 1; i <=size_nonred_array; i++) {
    tmr=Cudd_bddPermute(bdd_manager,dd_array1[i],var_index_array12_23);
    Cudd_Ref(tmr);
    for(j = 1; j <=size_nonred_array; j++){
      tmr2=Cudd_bddCompose(bdd_manager,tmr,dd_array1[j],j+size_nonred_array);
      Cudd_Ref(tmr2);
      Cudd_RecursiveDeref(bdd_manager,tmr);
      tmr=tmr2;
    }
    dd_array2[i]=tmr;
  }

  for(k = 0; k <size_nonred_array; k++){
    for(i = 1; i <=size_nonred_array; i++) {
      tmr=Cudd_bddPermute(bdd_manager,dd_array2[i],var_index_array12_23);
      Cudd_Ref(tmr);
      for(j = 1; j <=size_nonred_array; j++){
	tmr2=Cudd_bddCompose(bdd_manager,tmr,dd_array2[j],j+size_nonred_array);
	Cudd_Ref(tmr2);
	Cudd_RecursiveDeref(bdd_manager,tmr);
	tmr=tmr2;
      }
      Cudd_RecursiveDeref(bdd_manager,dd_array1[i]);
      dd_array1[i]=tmr;
    }
    tmpp=dd_array1;
    dd_array1=dd_array2;
    dd_array2=tmpp;
  }

  tmr = Cudd_ReadOne(bdd_manager);
  Cudd_Ref(tmr);
  for(i = 1; i <=size_nonred_array; i++) {
    tmr2 = Cudd_bddXnor(bdd_manager,dd_array2[i] , Cudd_bddIthVar(bdd_manager,i+size_nonred_array));
    Cudd_Ref(tmr2);
    tmr3=Cudd_bddAnd(bdd_manager,tmr, tmr2);
    Cudd_Ref(tmr3);
    Cudd_RecursiveDeref(bdd_manager, tmr);
    Cudd_RecursiveDeref(bdd_manager, tmr2);
    tmr=tmr3;
  }


  return tmr;
}

/* The functon returns pointer to bdd of next state computed for  "current state", both current state and next state depend on 1 set of variables */
DdNode*
Next_State(DdManager *bdd_manager,DdNode* current_state){

  DdNode* result, *tmp;
  int i;

  result=Cudd_ReadOne(bdd_manager);
  Cudd_Ref(result);

  for(i = 1; i <=size_nonred_array; i++) {
    tmp=Cudd_bddAnd( bdd_manager , current_state , dd_array_step[i]);
    if( tmp == current_state ){
      /* i-th var is one in next state*/
      tmp=Cudd_bddAnd(bdd_manager, result , Cudd_bddIthVar(bdd_manager,i) );
    }else{
      ASSERT( tmp == Cudd_Not( Cudd_ReadOne(bdd_manager) ) );
      /* i-th var is zero in next state*/
      tmp=Cudd_bddAnd(bdd_manager, result , Cudd_Not(Cudd_bddIthVar(bdd_manager,i)) );
    }
    Cudd_Ref( tmp );
    Cudd_RecursiveDeref( bdd_manager , result );
    result=tmp;
  }
  
  return result;

}