/* a_ucsf.c 
 *
 * This file consists of functions useful for 
 * identifying surface atoms by the UCSF approach
 * (UCSF = University of California at San Francisco) 
 * -----------------------------------------------------
 */

#include "main.h"

/* ------------------------------------------------------------------------- */
/* main function of this file 
 *  - uses the UCSF algorithm 
 *  - works on the input given as a parameter                                */
void identify_surf_at_ucsfa(SEZNAM *p_sf, char soub[]) {
/* ------------------------------------------------------------------------- */
     
/* definitions and declarations */
  SEZNAM *p_sa=NULL; /* auxiliary variable */
  float *min=NULL, *max=NULL; /* include information about min and max coordinate axis x, y, z */
  int i, j, k, pom; /* auxiliary variables */
  float prvniX, posledniX, prvniY, posledniY, prvniZ, posledniZ; /* auxiliary variables */
  int pocetBunekX=0, pocetBunekY=0, pocetBunekZ=0; /* auxiliary variables */
  int finX=0, finY=0, finZ=0; /* auxiliary variables */
  int ***M; /* matrix used by the UCSF algorithm */
  

/* determining min and max x-, y- and z-coordinates */
  if ((min = (float *) malloc (3*sizeof(float))) == NULL) {
      printf ("There is no more free space on heap.\nThe program will be terminated.\n");
      exit(1);
  }
  if ((max = (float *) malloc (3*sizeof(float))) == NULL) {
      printf ("There is no more free space on heap.\nThe program will be terminated.\n");
      exit(1);
  }
  
/* init */
    min[0] = pow(2.0,sizeof(float)*8-1)-1; /* maximum */
    min[1] = pow(2.0,sizeof(float)*8-1)-1;
    min[2] = pow(2.0,sizeof(float)*8-1)-1;
    max[0] = 0.0;
    max[1] = 0.0;
    max[2] = 0.0;
    
  /* finding the borders for M */
  p_sa = p_sf;
  if (p_sa == NULL) {
    printf("An internal error appeared during reading data about molecule.\n");
    printf("The program will be terminated.\n");
    exit(1);
  }

  while (p_sa != NULL) {
    if (p_sa->x < min[0]) {min[0] = p_sa->x;}  
    if (p_sa->y < min[1]) {min[1] = p_sa->y;}  
    if (p_sa->z < min[2]) {min[2] = p_sa->z;}  
    if (p_sa->x > max[0]) {max[0] = p_sa->x;}  
    if (p_sa->y > max[1]) {max[1] = p_sa->y;}  
    if (p_sa->z > max[2]) {max[2] = p_sa->z;}  
    p_sa = p_sa->dalsi_seznam;
  }
 
  /* init for definition of M */
    /* for axis X: finding the lower bound of M 
       (and reserve for the rest of needed cubes) */
  posledniX = prvniX = min[0]-(VAZBA-0.6)-VAZBA;
  i = 0;
  while (posledniX < max[0]) { posledniX += VAZBA; i++;}
  pocetBunekX = i+1; /* 1 extra is for a cube of right atom neighbour */

  posledniY = prvniY = min[1]-(VAZBA-0.6)-VAZBA;
  i = 0; 
  while (posledniY < max[1]) { posledniY += VAZBA; i++;}
  pocetBunekY= i+1;

  posledniZ = prvniZ = min[2]-(VAZBA-0.6)-VAZBA;
  i = 0;
  while (posledniZ < max[1]) { posledniZ += VAZBA; i++;}  
  pocetBunekZ = i+1;
  
  /* definition of matrix M */
     /* memory allocation for axis X */
  if ((M = (int ***) malloc((pocetBunekX) * sizeof(int**))) == NULL) {
      printf ("There is no more free space on heap.\nThe program will be terminated.\n");
      exit(1);            
  }
     /* memory allocation for axis Y */
  for (i=0; i<pocetBunekX; i++) {
     if ((M[i] = (int **) malloc((pocetBunekY) * sizeof(int*))) == NULL) {
         printf ("There is no more free space on heap.\nThe program will be terminated.\n");
         exit(1);             
     }
  }
     /* memory allocation for axis Z */
  for (i=0; i< pocetBunekX; i++) {
      for (j=0; j< pocetBunekY; j++) {
         if ((M[i][j] = (int *) malloc((pocetBunekZ) * sizeof(int))) == NULL) {
             printf ("There is no more free space on heap.\nThe program will be terminated.\n");
             exit(1);             
         }
      }
  }
  /* init of M */
  for (i=0; i<pocetBunekX; i++) {
    for (j=0; j<pocetBunekY; j++) {
       for (k=0; k<pocetBunekZ; k++) {
           M[i][j][k]=0;
       }
    }
  }

 /* --- main functinality -> the core of the algorithm --- */
 
 p_sa = p_sf;
 if (p_sa == NULL) {
    printf("An internal error appeared during reading data about molecule.\n");
    printf("The program will be terminated.\n");
    exit(1);
  }

 while (p_sa != NULL) {
     i=0;
     while (p_sa->x >= (i*VAZBA+prvniX)) {i++;} /* looks for X-coordinate */
     finX = i;
     i=0;
     while (p_sa->y >= (i*VAZBA+prvniY)) {i++;} /* looks for Y-coordinate */
     finY = i;
     i=0;
     while (p_sa->z >= (i*VAZBA+prvniZ)) {i++;} /* looks for Z-coordinate */
     finZ = i;

     /* cube M[finX-1][finY-1][finZ-1] contains the atom */

     for (i=finX-2; i<finX+1; i++) {
         for (j=finY-2; j<finY+1; j++) {
             for (k=finZ-2; k<finZ+1; k++) {
                 M[i][j][k]++; /* increasing grid */
             }
         }
     }
     /* remembering grid coordinates related to the atom in structure SEZNAM */
     p_sa->u_x = finX;
     p_sa->u_y = finY;
     p_sa->u_z = finZ;
     
     p_sa = p_sa -> dalsi_seznam;
 }

/* systematic trying to find blank matrix points (for determining the surface points)*/
  p_sa = p_sf;
  if (p_sa == NULL) {
    printf("An internal error appear during reading data about molecule.\n");
    printf("The program will be terminated.\n");
    exit(1);
  }

  while (p_sa != NULL) {
     pom = 0;
     /* cutting off the bounds of grid M which means that we have an empty grid point */
     if ((p_sa->u_x-3 < 0)||(p_sa->u_y-3 < 0)||(p_sa->u_z-3 < 0)) {
        p_sa->ucsf = 1; pom = 1; 
     }
     if ((p_sa->u_x+1 >= pocetBunekX)||(p_sa->u_y+1 >= pocetBunekY)||(p_sa->u_z+1 >= pocetBunekZ)) {
        p_sa->ucsf = 1; pom = 1;
     }
     if (pom == 0) {
       for (i=p_sa->u_x-3; i<p_sa->u_x+2; i++) {
          for (j=p_sa->u_y-3; j<p_sa->u_y+2; j++) {
              for (k=p_sa->u_z-3; k<p_sa->u_z+2; k++) {
                  if (M[i][j][k] ==0) {p_sa->ucsf = 1; pom=1; break;}
              }
              if (pom==1) {break;}
          }
          if (pom==1) {break;}
        }
     
     }
     p_sa = p_sa->dalsi_seznam;
  }
    
     
 /* deallocation */
  free ((void *) min);
  min = NULL;
  free ((void *) max);
  max = NULL;
  
  free ((void ***) M);
  M = NULL;
}

/* ------------------------------------------------------------------------- */
/* writes the results of the UCSF algorithm to the PBD file                  */
/* supposes that in != out (checked in higher function)                      */
void vypis_povrchove_atomy_ucsfa_soubor(SEZNAM *p_sf, char in[], char out[], int q){
/* ------------------------------------------------------------------------- */
  FILE *fw, *fr;
  SEZNAM *p_sa=NULL;
  char radek[SLOUPCE], outRadek[SLOUPCE], klicove_slovo[6];
  int uroven, i, counter=0;
  char *p_c;

  fr = fopen(in,"r");
  if (fr == NULL) {
    printf("The program was unable to open input file.\n");
    printf("It causes that executing of this program must be terminated.\n");
    exit(1);
  }
  fw = fopen(out,"w");
  if (fw == NULL) {
    printf("The program was unable to open output file.\n");
    printf("It causes that executing of this program must be terminated.\n");
    exit(1);
  }
  p_sa = p_sf;
  if (p_sa == NULL) {
    printf("An internal error appeared during reading data about molecule.\n");
    printf("The program will be terminated.\n");
    exit(1);
  }

  while (fgets(radek, SLOUPCE, fr) != NULL) { 
      strncpy(klicove_slovo,radek,6);
      klicove_slovo[6]='\0';
      p_c = strtok(klicove_slovo," ");
      p_c[strlen(p_c)]='\0';

      for (i=0; i<SLOUPCE+1; i++) {outRadek[i]=radek[i];}
      uroven=0;
      if (strcmp(p_c,"HETATM")==0) {uroven = 1;}
      if (strcmp(p_c,"ATOM")==0) {uroven = 2;} 
      switch (uroven) {
         case 0: /* all except for HETATM and ATOM */
	   fprintf(fw,"%s",outRadek);
           break;
 
         case 1: /* HETATM */
         case 2: /* ATOM */
	   if (p_sa->ucsf == 1) {
	     outRadek[17]='S';
	     outRadek[18]='U';
	     outRadek[19]='R';
	     counter++;
	   } else {
	     outRadek[17]='I';
	     outRadek[18]='N';
	     outRadek[19]='N';
	   }
	   fprintf(fw,"%s",outRadek);
	   p_sa = p_sa -> dalsi_seznam;
	   break;

         default : 
     	   /* opened file on this line does not begin with any of known keywords */
           fprintf(stderr,"\nFile %s does not have needed structure.\n",
                   in);
	   fclose(fr);
           exit(1);
      }

  }
  if (!q) {
     printf("\nResults:\n%d surface atoms identified.\n",counter);
     printf("Final PDB file: %s\n",out);
  }
  fclose(fr);
  fclose(fw);
}
/* ------------------------------------------------------------------------- */
/* writes the results of the UCSF approach to the PBD file */
/* supposes that in != out (checked in higher function)*/
void vypis_povrchove_atomy_ucsfa_struktura(SEZNAM *p_sf){
/* ------------------------------------------------------------------------- */
  SEZNAM *p_sa=NULL;

  p_sa = p_sf;
  if (p_sa == NULL) {
    printf("An internal error appeared during reading data about molecule.\n");
    printf("The program will be terminated.\n");
    exit(1);
  }

  while (p_sa != NULL) { 
	   if (p_sa->ucsf == 1) { /* SUR */ p_sa->sur_atom = 1; }
                       else { /* INN */ p_sa->sur_atom = 0; } 
	   p_sa = p_sa -> dalsi_seznam;
  }
}


/* ------------------------------------------------------------------------- */
/* writes results received by UCSF
 * algorithm for every single atom */
void vypis_ucsf_atomu(SEZNAM *p_sf) {
/* ------------------------------------------------------------------------- */
  SEZNAM *p_sa=NULL;

  p_sa = p_sf;
  while (p_sa != NULL) {
    printf("%d %s: UCSF = %d\n",p_sa->cislo,p_sa->ch_znacka,p_sa->ucsf);
    p_sa = p_sa->dalsi_seznam;
  }
}
