#include <stdio.h> 
#include <stdlib.h> 
#include <iostream> 
#include <math.h> 
 
float *matrix[110000];  
int *matrixind[110000];  
int size[110000];  
char line[100000]; 
int dim=108931; 
int dim_scop=7392; 
 
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()    // load the whole matrix 
{ 
  FILE *f; 
  int index,amount,i,j,k,ind; 
  char c; 
  char buffer[20]; 
  float eval;int flag=0; 
  size[0]=0; 
 
  f=fopen("rankprop_matrix.txt","r");  
 if (f==NULL) {printf("ERROR: File rankprop_matrix.txt not found!\n"); exit(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 *fscop; 
int scop_index=0,scop_size; 
float *scop_row=NULL; 
int *scop_rowind=NULL; 
 
 
void load_scop_matrix(int start)    // START to load the detailed scop matrix 
{ 
  int index,amount,i,j,k,ind; 
  char c;float eval;int flag=0; 
 
    fscop=fopen("rankprop_scop_matrix.txt","r");  
 if (fscop==NULL) {printf("ERROR: File rankprop_matrix.txt not found!\n"); exit(0);} 
 
 //printf("skipstart --> %d \n",start); 
  for(j=0;j<start;j++) // SKIP PAST THE START 
    { 
      //      printf("skip --> %d %d\n",j,start); 
      fscanf(fscop,"%d %d",&index,&amount); 
      for(i=0;i<amount;i++) 
       { fscanf(fscop,"%d",&ind);} 
      for(i=0;i<amount;i++) 
       {fscanf(fscop,"%g",&eval);} 
    } 
  printf("loaded query matrix, starting at query %d\n",start); 
} 
 
void next_scop_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(fscop,"%d %d",&scop_index,&scop_size); 
 
    if (!first_time) 
      {  
	delete scop_row; delete scop_rowind;  
      } else first_time=0; 
    scop_row= new float[scop_size]; 
    scop_rowind= new int[scop_size]; 
     
  for(i=0;i<scop_size;i++) 
    { fscanf(fscop,"%d",&ind); scop_rowind[i]=ind;} 
 
  for(i=0;i<scop_size;i++) 
    {fscanf(fscop,"%g",&eval); scop_row[i]=eval;} 
 
  // now normalize this row 
   
  for(i=0;i<scop_size;i++)  // calc exp, add up 
    {v=scop_row[i]; v=exp(-v/sigma); 
     scop_row[i]=v; sum+=v; 
    }  
  if(sum>1e-20)              // otherwise we get nans 
    for(i=0;i<scop_size;i++)  // so,  normalize 
      scop_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) 
{   
  float y[110000],yold[110000],yorig[110000];  
  int query_size,*row_ind; float *row_eval; char c; 
  float tmp2,tmp; 
  int i,j,k,tt,loops;  
   char fname[1000];  
 
   if (end>dim_scop) end=dim_scop;  
     
    sprintf(fname,"rm rankprop_results_%d.txt",start); system(fname); 
    // delete pre-existing file, we will concatenate each query on to this 
 
for(k=start+1;k<=end+1;k++) // all  examples - make matrix!!!! 
{ 
  next_scop_matrix(); // load next query in scop file 
   printf("%d %d %d \n",scop_index,scop_size,k);   
 
  for(i=0;i<=dim+2;i++)   // setup labels 
    {y[i]=0.0; yorig[i]=0;}     
  y[k]=1.0;               // query point  
 
  query_size=scop_size;  
  row_ind=scop_rowind; row_eval=scop_row; // take query from scop 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",start); 
  FILE *f=fopen(fname,"w"); 
  fprintf(f,"%d ",k-1); 
  for(int i=1;i<=dim_scop;i++) 
    fprintf(f,"%g ", y[i] ); 
  fprintf(f,"\n"); 
  fclose(f); 
  if(k>start)    
    sprintf(fname,"cat tmp.rankprop%d >> rankprop_results_%d.txt",start,start); 
  else 
    sprintf(fname,"cat tmp.rankprop%d > rankprop_results_%d.txt",start,start); 
 
  printf("[%s]\n",fname); 
  system(fname); 
} 
} 
  
 
 
 
 
 
int main(int argc, char *argv[]) 
{ 
 if (argc<3) 
   { 
     printf("Usage: rankprop <start_index> <end_index>\n\n"); 
     printf(" Loads the file rankprop_matrix.txt, which gives the PSI-BLAST distance measure as input.\n"); 
     printf(" Then, the RankProp algorithm is executed.\n The real-valued rank outputs for SCOP are saved in a file rankprop_results_<start_index>.\n"); 
     exit(0); 
   } 
 
  int start,end; 
  start=atoi(argv[1]); // can start at different positions  
                        // for parallel implementation 
 
  end=atoi(argv[2]); // can start at different positions  
                       // for parallel implementation 
 
  printf("start on %dth query\n",start); 
  printf("end on %dth query\n",end); 
 
  load_scop_matrix(start); // load queries from a different file 
  load_matrix(); 
 
                    
  normalize();      // do rbf, normalization    
 
  ranking(start,end);     // for each pt in scop, do ranking  
                          // output part of matrix to file   
 
  
 
}