#include <stdio.h> 
#include <stdlib.h> 
#include <iostream> 
#include <math.h> 
 
int dim=108931;    // dimension of database
int dim_query=7392;// dimension of query set

float **matrix;    // e-values of matrix edges  for protein i
int **matrixind;   // indices of matrix edges for protein i
int *size;         // number of edges for protein i
char *line[100000]; 
 
float sigma=100;   // sigma, controls size of local weights   
float alpha=0.95;  // alpha parameter, controls amount of "clustering" 
int max_loops=20;  // number of iterations of algorithm 
  
 
void load_matrix(char *fname)    // load the whole matrix 
{ 
  FILE *f; 
  int index,amount,i,j,k,ind; 
  char c; 
  char buffer[20]; 
  float eval;int flag=0; 
   
  f=fopen(fname,"r");  
  if (f==NULL) {printf("ERROR: File %s not found!\n",fname); exit(0);} 

  j=fscanf(f,"%d",&dim);
  printf("Proteins in database: %d\n",dim);  
  matrix=new float*[dim+1];
  matrixind=new int*[dim+1];
  size=new int[dim+1];
  size[0]=0; 

  for(k=1;k<dim;k++) 
  { 
    j=fscanf(f,"%d %d",&index,&amount); 
       if ((index%500)==0) 
        printf("loading %2.1f%% complete (%dth entry)... \n",100.0*(((float)index)/((float)dim)),index); 
    
    float *row= new float[amount]; 
    int *rowind= new int[amount]; 
     
  for(i=0;i<amount;i++) 
    { fscanf(f,"%d",&ind); rowind[i]=ind;} 
 
  for(i=0;i<amount;i++) 
    {fscanf(f,"%g",&eval); row[i]=eval;} 
 
   matrix[index]=row; matrixind[index]=rowind; size[index]=amount; 
  } 
  fclose(f); 
} 
 
 
FILE *fquery; 
int query_index=0,query_size; 
float *query_row=NULL; 
int *query_rowind=NULL; 
 
 
void load_query_matrix(char *fname,int start)   
 // START to load the detailed query matrix 
{ 
  int index,amount,i,j,k,ind; 
  char c;float eval;int flag=0; 
 
  fquery=fopen(fname,"r");  
  if (fquery==NULL) {printf("ERROR: File %s not found!\n",fname); exit(0);} 

  j=fscanf(fquery,"%d",&dim_query);
  printf("Number of query proteins: %d\n",dim_query);  
  
  for(j=0;j<start-1;j++) // SKIP PAST THE START 
    { 
      fscanf(fquery,"%d %d",&index,&amount); 
      for(i=0;i<amount;i++) 
       { fscanf(fquery,"%d",&ind);} 
      for(i=0;i<amount;i++) 
       {fscanf(fquery,"%g",&eval);} 
    } 
} 
 
void next_query_matrix() // get next query 
{ 
  static int first_time=1; 
  int i,j,k,ind; 
  char c;float eval;int flag=0; 
  float sum=0,v; 
 
    j=fscanf(fquery,"%d %d",&query_index,&query_size); 

    //printf("->%d %d\n",query_index,query_size);
 
    if (!first_time) 
      {  
	delete query_row; delete query_rowind;  
      } else first_time=0; 
    query_row= new float[query_size]; 
    query_rowind= new int[query_size]; 
     
  for(i=0;i<query_size;i++) 
    { fscanf(fquery,"%d",&ind); query_rowind[i]=ind;} 
 
  for(i=0;i<query_size;i++) 
    {fscanf(fquery,"%g",&eval); query_row[i]=eval;} 
 
  // now normalize this row 
   
  for(i=0;i<query_size;i++)  // calc exp, add up 
    {v=query_row[i]; v=exp(-v/sigma); 
     query_row[i]=v; sum+=v; 
    }  
  if(sum>1e-20)              // otherwise we get nans 
    for(i=0;i<query_size;i++)  // so,  normalize 
      query_row[i]/=sum; 
   
} 
 
 
 
void normalize()  // normalize data to be a transition matrix 
{ 
  char c; 
  printf("[normalizing]\n"); 
  int i,j; float sum,v; 
  for(j=1;j<=dim;j++)  // do each row in turn ...  
    { 
      sum=0; 
 
      for(i=0;i<size[j];i++)  // calc exp, add up 
	{v=matrix[j][i]; v=exp(-v/sigma); 
	matrix[j][i]=v; sum+=v; 
	} 
 
      if(sum>1e-20)              // otherwise we get nans 
	for(i=0;i<size[j];i++)  // so,  normalize 
	  matrix[j][i]/=sum; 
    } 
  //  printf("         [...done.]\n"); 
   
} 
 
 
 
 
void ranking(int start, int end, char *fnam) 
{   
  float *y,*yold,*yorig;
  y=new float[dim+1];
  yold=new float[dim+1];
  yorig=new float[dim+1];
  
  int *row_ind; float *row_eval; char c; 
  float tmp2,tmp; 
  int i,j,k,tt,loops;  
  char fname[1000];  
     
   if (end>dim_query) end=dim_query;  
     
    sprintf(fname,"rm %s_%d-%d",fnam,start,end); system(fname); 
    // delete pre-existing file, we will concatenate each query on to this 
 
for(k=start;k<=end;k++) // all  examples - make matrix!!!! 
{ 
  next_query_matrix(); // load next query in query file 
  // printf("%d %d %d \n",query_index,query_size,k);   
 
  for(i=0;i<=dim+2;i++)   // setup labels 
    {y[i]=0.0; yorig[i]=0;}     
  y[k]=1.0;               // query point  

  row_ind=query_rowind; row_eval=query_row; // take query from query file 
 
  for(i=0;i<query_size;i++)  // calculate yorig - distance of each pt to query 
    yorig[row_ind[i]]=row_eval[i];      
   
   //******* start of loop ******** 
 
  for(loops=0;loops<max_loops;loops++)  // iterations of propagation 
    {   
      for(i=0;i<=dim+1;i++) yold[i]=y[i];  
 
      for(j=1;j<=dim;j++)  
	{ 
	  if(j==k) continue;           // don't update query 
          row_ind=matrixind[j]; row_eval=matrix[j]; 
	  tmp=0;  
 
	  for(i=0;i<size[j];i++)  
	     { 
	       tmp+=row_eval[i]*yold[row_ind[i]]; // weighted comb of neighbors 
	    }  
	  y[j]=(alpha*tmp)+yorig[j];  // new y is old y +  
	                             // new weighted combination of neighbors 
	}  
      //printf("query %d iteration %d : y[2]=%g \n",k-1,loops,y[2]); 
    } 
 
  sprintf(fname,"tmp.rankprop%d-%d",start,end); 
  FILE *f=fopen(fname,"w"); 
  fprintf(f,"%d ",k); 
  for(int i=1;i<=dim;i++) 
    fprintf(f,"%g ", y[i] ); 
  fprintf(f,"\n"); 
  fclose(f); 
  if(k>start)    
    sprintf(fname,"cat tmp.rankprop%d-%d >> %s_%d-%d",start,end,fnam,start,end); 
  else 
    sprintf(fname,"cat tmp.rankprop%d-%d > %s_%d-%d",start,end,fnam,start,end); 
 
  //printf("[%s]\n",fname); 
  printf("query %d ranked - %2.1f%% complete\n",k,100.0*(((float)k+1-start)/((float)end-start+1)));

  system(fname); 
} 

sprintf(fname,"rm tmp.rankprop%d-%d",start,end); // clean up temp file
system(fname);
} 
 

 
int main(int argc, char *argv[]) 
{ 
 if (argc<6) 
   { 
     printf("Usage: rankprop <matrix_file> <query_file> <output_file> <start> <end> [-s <value>] [-i <value>] [-a <value>]\n\n"); 
     printf("Performs the propagation ranking algorithm on the given queries in <query_file>.\n");
     printf("<start> and <end> are indices into the query_matrix in order to run a subset of the queries.\n");
     printf("For each query, propagation is performed on the network given in <matrix_file>.\n");
     printf("The ranking of the examples in the network is given as output by <output_file> \n(the ranking scores of each item in the matrix file).\n"); 
     printf("-a,-s and -i options set the parameters slpha, sigma and the iterations of the algorithm.");
     printf("\n  FILE FORMATS:\n");
     printf("  Both <matrix_file> and <query_file> are of the following repeated format:\n"); 
     printf("  <number of entries>\n");
     printf("  for each entry:\n");
     printf("    <index> <number of edges> <indices of edges> <e-values of edges>.\n"); 
     printf("  (Indices should start from 1.)\n");
     exit(0); 
   } 
  
/*
Example matrix file:
2          <--- number of proteins
1 2        <--- protein 1 has 2 edges
1 2        <--- edge indices are to proteins 1 & 2
0.1 100    <--- with the following weights
2 2        <--- protein 2 also has 2 edges
1 2        <--- edge indices are also to proteins 1 & 2 (fully connected)
100 0.1    <--- with the following weights (small means highly similar)
*/


  for(int i=6;i<argc;i+=2) // get optional parameters
  {
    if (argv[i][1]=='a') alpha= atoi(argv[i+1]);
    if (argv[i][1]=='s') sigma= atoi(argv[i+1]);
    if (argv[i][1]=='i') max_loops= atoi(argv[i+1]);
  }	

  printf("Running RankProp.\n");
  printf("Parameters: alpha=%g sigma=%g iterations=%d.\n",alpha,sigma,max_loops);

  int start,end; 
  start=atoi(argv[4]); // can start at different positions  
  end=atoi(argv[5]);  // for parallel implementation 

  load_matrix(argv[1]); 
  load_query_matrix(argv[2],start); // load queries from a different file 
  printf("Only run queries %d-%d\n",start,end); 
                   
  normalize();            // do rbf, normalization    
 
  ranking(start,end,argv[3]);       // for each query, do ranking  
                                   // output part of matrix to file   
  
}