testRigidCLL.cpp

Go to the documentation of this file.

#include "RigidCLL.h"

#include <iostream>
#include <algorithm> 
#include <stdlib.h>    
#include <time.h> 
#include <math.h>
#include <unistd.h>
#include <sys/resource.h>

#define RAND_DOUBLE(s) (((double)rand()*2*s)/(double)RAND_MAX-s)

#define N_EXPERIMENTS 50

#define REPETITIONS 10000

double getTime();

void RandomTranslation(double sx,double sy,double sz,RigidTransform *t);

void PrintInteractions(const string &label,list<AtomPair> *interactions);

bool CmpInteractions(list<AtomPair> *interactions1,list<AtomPair> *interactions2);



/********************************************************************************/
/********************************************************************************/
/********************************************************************************/

double getTime()
{
  struct rusage usage;

  getrusage(RUSAGE_SELF,&usage);

  /* Accumulate the user and system time and convert to miliseconds */
  return((double)(usage.ru_utime.tv_sec +usage.ru_stime.tv_sec )*1000.0+
         (double)(usage.ru_utime.tv_usec+usage.ru_stime.tv_usec)/1000.0);
}

void RandomTranslation(double sx,double sy,double sz,RigidTransform *t)
{
  (*t)[0][0]=1;(*t)[0][1]=0;(*t)[0][2]=0;(*t)[0][3]=RAND_DOUBLE(sx);
  (*t)[1][0]=0;(*t)[1][1]=1;(*t)[1][2]=0;(*t)[1][3]=RAND_DOUBLE(sy);
  (*t)[2][0]=0;(*t)[2][1]=0;(*t)[2][2]=1;(*t)[2][3]=RAND_DOUBLE(sz);
  (*t)[3][0]=0;(*t)[3][1]=0;(*t)[3][2]=0;(*t)[3][3]=1;
}

void PrintInteractions(const string &label,list<AtomPair> *interactions)
{
  unsigned int i,ni;
  list<AtomPair>::iterator it;

  cout << "  " << label << "\n";
  interactions->sort();
  ni=interactions->size();
  cout << "    Number of interactions: " << ni << "\n    Interactions:";
  it=interactions->begin();
  if (ni<20)
    {
      for(i=0;i<ni;i++,++it)
        cout << (*it);
    }
  else
    {
      for(i=0;i<ni;i++,++it)
        {
          if (i<5)
            cout << (*it);
          else
            {
              if (i==5)
                cout << "   .....  ";
              else
                {
                  if (i>ni-6)
                    cout << (*it);
                }
            }
        }
    }
  cout << "\n\n";
}

bool CmpInteractions(list<AtomPair> *interactions1,list<AtomPair> *interactions2)
{
  unsigned int n1,n2,i;
  bool equal;
  list<AtomPair>::iterator it1,it2;

  interactions1->sort();
  interactions2->sort();

  n1=interactions1->size();
  n2=interactions2->size();

  equal=(n1==n2);
  it1=interactions1->begin();
  it2=interactions2->begin();
  for(i=0;((equal)&&(i<n1));i++,++it1,++it2)
    equal=((*it1)==(*it2));

  return(equal);
}

int main(int argc, char **arg)
{
  if (argc>1)
    {      
      RigidCLL molecule(arg[1]); /* reads the rigid information from a file */
      unsigned int i,nr;
      double sx,sy,sz;
      RigidTransform *t;
      bool clashRCLL,clashCLL,clashBF;
      list<AtomPair> interactionsRCLL;
      list<AtomPair> interactionsCLL;
      list<AtomPair> interactionsBF;
      double cutoff;
      unsigned int ne,rep;
      double start,end;
      double time1,time2;
      double ts;
      unsigned int rs;
      
      ts=(double)sysconf(_SC_CLK_TCK); /*tics per second*/
      
      /* Init the random number generator */
      rs=(unsigned int)time(NULL);
      srand(rs);

      /* Set up RigidCLL for steric clash */
      molecule.SetCutOff(0.0);

      /* Get information about the size of the problem and the enviroment */
      nr=molecule.GetNumRigids();
  
      /* Initialize the transform. */
      t=new RigidTransform[nr];
      for(i=0;i<nr;i++)
        RandomTranslation(0,0,0,&(t[i])); 
      
      /* initialize the collision system. Wee need this to have the
         bounding box available */
      molecule.WarmUp(t);

      /* Get information about the bounding box to adjust
         the size of the random perturbations */
      sx=molecule.GetBBSize(0);
      sy=molecule.GetBBSize(1);
      sz=molecule.GetBBSize(2);

      /* First we detect steric clashes */

      /* Random motion of all rigids but the 'fixed'. */
      for(i=1;i<nr;i++)
        RandomTranslation(0.05*sx,0.05*sy,0.05*sz,&(t[i])); 

      /* Now determine if there is a collision */
      /* First using the RigidCLL algorihtm */
      cout << "************************************************\n";
      cout << "Steric clash test\n\n";
      cout << "  RigidCLL   : ";
      clashRCLL=molecule.StericClash(t);
      if (clashRCLL)
        cout << "  Steric clash detected between the rigids\n";
      else
        cout << "  No steric clash detected between the rigids\n";

     /* now using the CLL approach */
      cout << "  CLL        : ";
      clashCLL=molecule.CLLStericClash(t);
      if (clashCLL)
        cout << "  Steric clash detected between the rigids\n";
      else
        cout << "  No steric clash detected between the rigids\n";

      /* and now the brute force approach to verify the results */
      cout << "  Brute force: ";
      clashBF=molecule.BruteForceStericClash(t);
      if (clashBF)
        cout << "  Steric clash detected between the rigids\n\n";
      else
        cout << "  No steric clash detected between the rigids\n\n";

      /* Verify the output */
      if (clashRCLL!=clashBF)
        cout << "  ERROR: discrepancy between RigidCLL and Brute Force \n";
      if (clashCLL!=clashBF)
        cout << "  ERROR: discrepancy between CLL and Brute Force \n";

      /* Repeat the execution many times with different perturbations
         to get average execution time */

      cout << "  Computing average execution time: " ;
      cout.flush();
      time1=0.0;
      time2=0.0;
      for(ne=0;ne<N_EXPERIMENTS;ne++)
        {
          for(i=1;i<nr;i++)
            RandomTranslation(0.05*sx,0.05*sy,0.05*sz,&(t[i])); 
          
          /* ensure a clean grid and remove the time to generate
             the grid out of the query time */
          molecule.WarmUp(t);

          start=getTime();
          for(rep=0;rep<REPETITIONS;rep++)
            molecule.StericClash(t);
          end=getTime();
          time1+=(end-start);

          /* ensure a clean grid and remove the time to generate
             the grid out of the query time */
          molecule.WarmUp(t);

          start=getTime();
          for(rep=0;rep<REPETITIONS;rep++)
            molecule.CLLStericClash(t);
          end=getTime();
          time2+=(end-start);

          cout << ".";cout.flush();
        }
      cout << "\n";

      time1/=(double)(N_EXPERIMENTS*REPETITIONS);
      cout<<"    Execution time RigidCLL (ms): " << time1 << "\n";

      time2/=(double)(N_EXPERIMENTS*REPETITIONS);
      cout<<"    Execution time CLL      (ms): " << time2 << "\n";

      cout<<"  Time ratio RigidCLL/CLL       : " << time1/time2*100 << " %\n";

      /* Now detect interactions within a given cutoff */

      cout << "\n************************************************\n";

      /* Set up RigidCLL for interaction detection */
      /* get the user provide cutOff or use the default one */
      if (argc>2)
        cutoff=atof(arg[2]);
      else
        cutoff=2.0;
      /* avoid degenerate cutoff values */
      if (cutoff<=0)
        cutoff=2;

      /* Set the new cutoff. This will re-initialize the problem */
      molecule.SetCutOff(cutoff);

      cout << "Interaction test with cutoff=" << molecule.GetCutOff() << "\n\n";
      
      /* Set up a new random perturbation for the interaction test */
      for(i=1;i<nr;i++)
        RandomTranslation(0.05*sx,0.05*sy,0.05*sz,&(t[i]));
      
      /* remove the time to create the grid out of the query */
      molecule.WarmUp(t);

      /* First detect the interactions using the RigidCLL algorithm */
      molecule.Interactions(t,&interactionsRCLL);
      PrintInteractions("RigidCLL",&interactionsRCLL);

     /* now using the original CLL algorithm */
      molecule.CLLInteractions(t,&interactionsCLL);
      PrintInteractions("CLL (without considering rigid clusters at all)",&interactionsCLL);

      /* and finally the brute force approach to verify the results */
      molecule.BruteForceInteractions(t,&interactionsBF);
      PrintInteractions("Brute Force (taking into account the rigid clusters)",&interactionsBF);

      /* Verify the output assuming that Brute Force is correct */
      cout << "  Verifying the outputs:\n";
      if (CmpInteractions(&interactionsRCLL,&interactionsBF))
        cout << "    The output of RigidCLL is correct (equal to Brute Force)\n\n";
      else
        cout << "    ERROR: discrepancy between RigidCLL and Brute Force\n\n";

      interactionsRCLL.clear();
      interactionsCLL.clear();
      interactionsBF.clear();

      /* Repeat the execution many times with different perturbations
         to get average execution time */
      cout << "  Computing average execution time: " ;
      cout.flush();

      time1=0.0;
      time2=0.0;
      for(ne=0;ne<N_EXPERIMENTS;ne++)
        {
          for(i=1;i<nr;i++)
            RandomTranslation(0.05*sx,0.05*sy,0.05*sz,&(t[i]));

          /* ensure a clean grid and remove the time to generate
             the grid out of the query time */
          molecule.WarmUp(t);
        
          for(rep=0;rep<REPETITIONS;rep++)
            {
              start=getTime();
              molecule.Interactions(t,&interactionsRCLL);
              end=getTime();
              time1+=(end-start);
              interactionsRCLL.clear();
            }

          /* ensure a clean grid and remove the time to generate
             the grid out of the query time */
          molecule.WarmUp(t); 

          for(rep=0;rep<REPETITIONS;rep++)
            {
              start=getTime();
              molecule.CLLInteractions(t,&interactionsCLL); 
              end=getTime();
              time2+=(end-start);
              interactionsCLL.clear();
            }
          
          cout << ".";cout.flush();
        }
      cout << "\n";

      time1/=(double)(N_EXPERIMENTS*REPETITIONS);
      cout<<"    Execution time RigidCLL (ms): " << time1 << "\n";

      time2/=(double)(N_EXPERIMENTS*REPETITIONS);
      cout<<"    Execution time CLL      (ms): " << time2 << "\n";

      cout<<"    Time ratio                  : " << time1/time2*100 << " %\n";

      cout << "\n************************************************\n\n";

      /* Release the allocated memory */
      delete [] t;

      /* and exit */
      return(EXIT_SUCCESS);
    }
  else
    {
      cerr << "Wrong number of parameters\n";
      cerr << "Use\n";
      cerr << "    testRigidCLL filename [cutoff]\n";
      cerr << "where 'filename' is the file with the atom/rigid information\n";
      cerr << "      'cutoff' is the theshold to use when detecting interacting atoms\n";
      cerr << "               This parameter is optional. If not give, 2 is used.\n";
      return(EXIT_FAILURE);
    }
}