cuikenergy.c

Go to the documentation of this file.

#include "filename.h"
#include "box_list.h"
#include "error.h"
#include "defines.h"
#include "parameters.h"
#include "bioworld.h"


#include <stdio.h>
#include <math.h>

#define TRUE_COSTS 0

int main(int argc, char **arg)
{
  Tparameters parameters;
  TBioWorld bioWorld;
  Tfilename fparam,fsols,fcost;
  Tbox *box;
  FILE *fout;
  unsigned int n,k,nBoxes;
  double *point;
  Tlist boxList;
  Titerator i;
  double *cost,minCost,maxCost;
  double *conf;
  unsigned int nMin,nMax;
  double e;

  if (argc>2)
    { 
      /*Init parameters*/
      CreateFileName(NULL,arg[1],NULL,PARAM_EXT,&fparam);
      fprintf(stderr,"Reading parameters                : %s\n",GetFileFullName(&fparam));
      InitParametersFromFile(GetFileFullName(&fparam),&parameters);
      
      InitBioWorld(&parameters,arg[1],NO_UINT,&conf,&bioWorld);

      CreateFileName(NULL,arg[2],NULL,SOL_EXT,&fsols);
      fprintf(stderr,"Reading solution file             : %s\n",GetFileFullName(&fsols));
      ReadListOfBoxes(GetFileFullName(&fsols),&boxList);
      nBoxes=ListSize(&boxList);
      if (nBoxes==0)
        Error("Empty solution file in cuikenergy");

      /* Compute the energy for each conformation */
      InitIterator(&i,&boxList);
      First(&i);
      k=0;
      NEW(cost,nBoxes,double);
      while(!EndOfList(&i))
        {
          box=(Tbox *)GetCurrent(&i);
          if (k==0)
            {
              n=GetBoxNIntervals(box);
              NEW(point,n,double);
            }
          GetBoxCenter(NULL,point,box);
          e=BioWorldEnergy(&parameters,FALSE,point,(void *)&bioWorld);
          #if (TRUE_COSTS)
            cost[k]=e; 
          #else
            cost[k]=log(e);
          #endif

          fprintf(stderr,"%u -> %g (%g)\n",k,cost[k],e);

          if (k==0)
            {
              minCost=cost[k]; nMin=k;
              maxCost=cost[k]; nMax=k;
            }
          else
            {
              if (cost[k]>maxCost)
                { maxCost=cost[k]; nMax=k; }
              else
                {
                  if (cost[k]<minCost)
                    { minCost=cost[k]; nMin=k; }
                }
            }
          k++;
          Advance(&i);
        }
      fprintf(stderr,"Extreme energies                  : %g (%u) -- %g (%u)\n",minCost,nMin,maxCost,nMax);

      /* Normalize and write */
      CreateFileName(NULL,arg[2],NULL,COST_EXT,&fcost);
      fprintf(stderr,"Writing cost file                 : %s\n",GetFileFullName(&fcost));

      fout=fopen(GetFileFullName(&fcost),"w");
      if (!fout)
        Error("Output file for costs can not be opened");
      
      //if (maxCost>4*minCost)
      //        maxCost=4*minCost;
      for(k=0;k<nBoxes;k++)
        {
          //if (cost[k]>maxCost)
          //  fprintf(fout,"1e-5\n");
          //else
            {
              /* Low energy is good -> low values must be mapped to values
                 close to 1 to get green */
            
              #if (TRUE_COSTS)
                fprintf(fout,"%g\n",cost[k]);
              #else  
                e=((cost[k]-minCost)/(maxCost-minCost))*(1-1e-3)+1e-3;
                //e=pow(e,0.5);
                fprintf(fout,"%g\n",e);
              #endif
            }
        }
      fclose(fout);

      DeleteFileName(&fsols);
      DeleteFileName(&fcost);
      DeleteFileName(&fparam);
      
      free(conf);
      free(point);
      free(cost);
      DeleteListOfBoxes(&boxList);
      DeleteParameters(&parameters);
      DeleteBioWorld(&bioWorld);
    }
  else
    {
      fprintf(stderr,"  Wrong number of parameters.\n");
      fprintf(stderr,"  Use:\n");
      fprintf(stderr,"    cuikenergy <problem_name>.pdb  <sol_name>\n");
      fprintf(stderr,"   where:\n");
      fprintf(stderr,"      <problem_name> the input molecular information file (pdb).\n");
      fprintf(stderr,"      <sol_name> Configurations to evaluate.\n");
      fprintf(stderr,"  The extension (e.g., '.pdb') is required\n");
      fprintf(stderr,"  All the extensions managed by OpenBabel can be used\n");
      fprintf(stderr,"  The ouput file is <sol_name>.cost\n");
    }

  return(EXIT_SUCCESS);
}