test-cuik-kdtree.c

Go to the documentation of this file.

#include "cuik-kdtree.h"
#include "random.h"
#include "definitions.h"
#include "vector.h"

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/resource.h>
#include <DNN/ANN_C.h>

#define POINTS_ON_SPHERE 0

#define N_POINTS 50000

#define DIM 8

#define TOPOLOGY TOPOLOGY_S

#define POINTS_IN_LEAF 25

#define SAMPLING_EXPANSION 1.0

#define R 2.0

#define REP 1000


#define N_OPTIONS 9

int CmpUInt(const void *a,const void *b);

double getTime();

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

int CmpUInt(const void *a,const void *b)
{
  return(*((unsigned int*)a)-*((unsigned int*)b));
}

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);
}

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

int main(int argc, char **arg)
{
  char *option[N_OPTIONS]={"-h","-n","-d","-t","-l","-r","-m","-s","-e"};

  unsigned int nPoints,dim,tpg,pLeaf,rep;
  double searchR,expansion;

  /* The two kd-trees to compare */
  TKDtree *Mkdt;
  #ifdef _MPNN
    TKDTree *Ykdt;
  #endif

  /* Several auxiliary variables */
  boolean found;
  unsigned int i,j;
  double **point,*q;
  Trectangle ambient;
  unsigned int *topology;
  unsigned int id;
  double l,u,d,r2;
  unsigned int *nn,n,m;
  double *dst,dMin;
  unsigned int idMin;

  /* Random seed */
  unsigned int ri;
  
  /* To measure time */
  time_t t;
  double tStart,tEnd;

  /* Set the parameters to the default values */
  nPoints=N_POINTS;
  dim=DIM;
  tpg=TOPOLOGY;
  pLeaf=POINTS_IN_LEAF;
  expansion=SAMPLING_EXPANSION;
  searchR=R;
  rep=REP;

  /* Process the input parameters, if any */
  i=1;
  while (i<argc)
    {
      found=FALSE;
      j=0;
      while((!found)&&(j<N_OPTIONS))
        {
          if (strcmp(arg[i],option[j])==0)
            {
              switch(j)
                {
                case 0: /* -h */
                  fprintf(stderr,"cuik-kdtree test program\n");
                  fprintf(stderr,"Usage: \n");
                  fprintf(stderr,"  test-cuik-kdtree <option> <ooption> <option> ...\n");
                  fprintf(stderr,"Options: \n");
                  fprintf(stderr,"  -n <value>   Number of random  points to generate.\n");
                  fprintf(stderr,"                 Default value: %u\n",N_POINTS);
                  fprintf(stderr,"  -d <value>   Dimension of the  points to generate.\n");
                  fprintf(stderr,"                 Default value: %u\n",DIM);
                  fprintf(stderr,"  -t <R|S>     Topology (R or S) of the  points to generate.\n");
                  fprintf(stderr,"               The same topology is used for all the dimensions.\n");
                  fprintf(stderr,"                 Default value: %s\n",(TOPOLOGY==TOPOLOGY_R?"R":"S"));
                  fprintf(stderr,"  -l <value>   Maximum number of points in each leaf.\n");
                  fprintf(stderr,"                 Default value: %u\n",POINTS_IN_LEAF);
                  fprintf(stderr,"  -r <value>   Search radius (for the in-ball queries).\n");
                  fprintf(stderr,"                 Default value: %f\n",R);
                  fprintf(stderr,"  -m <value>   Number of reptetions to gather execution times.\n");
                  fprintf(stderr,"                 Default value: %u\n",REP);
                  fprintf(stderr,"  -e <value>   Sampling expansion (not used so far in the test).\n");
                  fprintf(stderr,"                 Default value: %f\n",SAMPLING_EXPANSION);
                  fprintf(stderr,"  -s <value>   Random seed.\n");
                  fprintf(stderr,"                 Default value: arbitrarily fixed\n");
                  fprintf(stderr,"  -h           This help.\n"); 
                  fprintf(stderr,"Examples: \n");
                  fprintf(stderr,"  test-cuik-kdtree\n"); 
                  fprintf(stderr,"  test-cuik-kdtree -t R -r 0.75\n"); 
                  fprintf(stderr,"  test-cuik-kdtree -n 100000\n"); 
                  exit(EXIT_SUCCESS);
                  break;

                case 1: /* -n */
                  found=TRUE;
                  i++;
                  nPoints=(unsigned int)atoi(arg[i]);
                  if (nPoints<2)
                    {
                      fprintf(stderr,"Too smalll number of points (-n).\n");
                      exit(EXIT_FAILURE);
                    }
                  break;

                case 2: /* -d */
                  found=TRUE;
                  i++;
                  dim=(unsigned int)atoi(arg[i]);
                  if (dim<1)
                    {
                      fprintf(stderr,"Too smalll dimension (-d).\n");
                      exit(EXIT_FAILURE);
                    }
                  break;

                case 3: /* -t */
                  found=TRUE;
                  i++;
                  
                  if (arg[i][0]=='R')
                    tpg=TOPOLOGY_R;
                  else
                    {
                      if (arg[i][0]=='S')
                        tpg=TOPOLOGY_S;
                      else
                        {
                          fprintf(stderr,"Unknown topology (-t).\n");
                          exit(EXIT_FAILURE);
                        }
                    }
                  break;

                case 4: /* -l */
                  found=TRUE;
                  i++;
                  pLeaf=(unsigned int)atoi(arg[i]);
                  if (pLeaf<1)
                    {
                      fprintf(stderr,"Too smalll number of points in each leaf (-l).\n");
                      exit(EXIT_FAILURE);
                    }
                  break;

                case 5: /* -r */
                  found=TRUE;
                  i++;
                  searchR=(double)atof(arg[i]);
                  if (searchR<=0)
                    {
                      fprintf(stderr,"The search radius must be positive (-r).\n");
                      exit(EXIT_FAILURE);
                    }
                  break;

                case 6: /* -m */
                  found=TRUE;
                  i++;
                  rep=(unsigned int)atoi(arg[i]);
                  if (rep<1)
                    {
                      fprintf(stderr,"The number of repetitions must be positive (-m).\n");
                      exit(EXIT_FAILURE);
                    }
                  break;

                case 7: /* -s */
                  found=TRUE;
                  i++;
                  ri=(unsigned int)atof(arg[i]);
                  break;

                case 8: /* -e */
                  found=TRUE;
                  i++;
                  ri=(unsigned int)atof(arg[i]);
                  break;
                }
            }
          j++;
        }

      if (!found)
        {
          fprintf(stderr,"Unknown command line option %s.\nUse -h for help.\n",arg[i]);
          exit(EXIT_FAILURE);
        }
      i++;
    }

  /* Init the random number generator  */
  t=time(NULL);  /* Get the time at which input files have been read */
  ri=(unsigned int)t;
  //ri=1386189846;
  randomSet(ri);

  /* Print information about */
  printf("\n\n************************************************************\n");
  printf("Parameters\n");
  printf("  Random seed   : %u\n",ri);
  printf("  Topology      : %s\n",(tpg==TOPOLOGY_R?"R":"S"));
  printf("  No. points    : %u\n",nPoints);
  printf("  Dimension     : %u\n",dim);
  printf("  Points in leaf: %u\n",pLeaf);
  printf("  Search radius : %f\n",searchR);
  if (rep>1)
    printf("  Repetitions   : %u\n",rep);
  printf("\n");

  /* define the ambient space */
  InitRectangle(dim,NULL,NULL,&ambient);
  
  if (tpg==TOPOLOGY_R)
    {
      l=-1;
      u=1;
      topology=NULL;
    }
  else
    {
      l=-M_PI;
      u=+M_PI;
      NEW(topology,dim,unsigned int);
      for(i=0;i<dim;i++)
        topology[i]=tpg; /* == TOPOLOGY_S */
    }
   for(i=0;i<dim;i++)
     {
       SetRectangleLowerLimit(i,l,&ambient);
       SetRectangleUpperLimit(i,u,&ambient);
     }

  /* Define the points: uniform distribution on a dim-sphere */
  printf("Generating %u random points\n\n",nPoints);
  NEW(point,nPoints,double*);
  for(i=0;i<nPoints;i++)
    {
      NEW(point[i],dim,double);
      #if (POINTS_ON_SPHERE)
        for(j=0;j<dim;j++)
          point[i][j]=randomNormal();
        VectorNormalize(dim,point[i]);
      #else
        RandomPointInRectangle(point[i],&ambient);
      #endif
    }

  /* Build the two kd-trees */
  printf("Building kd-tree\n");
  #ifdef _MPNN
    tStart=getTime();
    if (topology==NULL)
      Ykdt=CreateKDTree(dim);
    else
      Ykdt=CreateKDTreeT(dim,(int *)topology);
    for(i=0;i<nPoints;i++)
      AddPoint2KDTree(point[i],Ykdt);
    tEnd=getTime();
    printf("  MPNN        : %7.3f ms\n",tEnd-tStart);
    printf("  MPNN x point: %7.3f ms\n",(tEnd-tStart)/(double)nPoints);
  #endif

  tStart=getTime();
  for(i=0;i<nPoints;i++)
    {
      if (i==0)
        Mkdt=InitKDTree(&ambient,pLeaf,expansion,1,&(point[i]),&i);
      else
        AddPoint2KDtree(topology,point[i],i,&Mkdt);
    }
  tEnd=getTime();
  printf("  Cuik        : %7.3f ms\n",tEnd-tStart);
  printf("  Cuik x point: %7.3f ms\n\n",(tEnd-tStart)/(double)nPoints);
  //printf("    Height    : %u (%u)\n\n",Mkdt->height,(unsigned int)ceil(1+log2((double)nPoints/(double)pLeaf))); 

  /* Define the query point */
  printf("Defining the query point\n\n");
  NEW(q,dim,double);
  RandomPointInRectangle(q,&ambient);

  /* Query the trees */
  #ifdef _MPNN
    printf("Searching NN\n");
    printf("  MPNN kd-tree      : ");
    QueryKDTree(q,(int *)&id,&d,Ykdt);
    printf("%u\n",id);
  #endif

  printf("  Cuik kd-tree      : ");
  NearestNeighbour(topology,q,Mkdt);  
  printf("%u\n",id);

  printf("  Brute-force search: ");  
  dMin=INF;
  for(i=0;i<nPoints;i++)
    {
      d=VectorSquaredDistanceTopologyMin(dMin,dim,topology,point[i],q);
      if (d<dMin)
        {
          idMin=i;
          dMin=d;
        }
    }
  printf("%u\n\n",idMin);

  printf("Searching in-ball\n");
  #ifdef _MPNN
    printf("  MPNN kd-tree      : ");
    QueryKDTreeR(q,searchR,(int *)(&n),(int **)(&nn),&dst,Ykdt);
    if (n>0)
      {
        qsort(nn,n,sizeof(unsigned int),CmpUInt);
        printf("[%u] ",n);
        if (n>10)
          {
            for(i=0;i<5;i++)
              printf("%u ",nn[i]);
            printf("... ");
            for(i=n-5;i<n;i++)
              printf("%u ",nn[i]);
          }
        else
          {
            for(i=0;i<n;i++)
              printf("%u ",nn[i]);
          }
        if (n>0)
          {
            free(nn);
            free(dst);
          }
      }
    else
      printf("no point");
    printf("\n");
  #endif

  printf("  Cuik kd-tree      : ");
  NeighboursInBall(topology,q,searchR,&n,&nn,Mkdt);
  if (n>0)
    {
      qsort(nn,n,sizeof(unsigned int),CmpUInt);
      printf("[%u] ",n);
      if (n>10)
        {
          for(i=0;i<5;i++)
            printf("%u ",nn[i]);
          printf("... ");
          for(i=n-5;i<n;i++)
            printf("%u ",nn[i]);
        }
      else
        {
          for(i=0;i<n;i++)
            printf("%u ",nn[i]);
        }
      if (n>0)
        free(nn);
    }
  else
    printf("no point");
  printf("\n");
      
  printf("  Brute-force search: ");  
  r2=searchR*searchR;
  n=0;
  m=100;
  NEW(nn,m,unsigned int);
  for(i=0;i<nPoints;i++)
    {
      if (VectorSquaredDistanceTopologyMin(r2,dim,topology,point[i],q)<r2)
        {
          if (n==m)
            MEM_DUP(nn,m,unsigned int);
          nn[n]=i;
          n++;
        }
    }
  if (n>0)
    {
      printf("[%u] ",n);
      if (n>10)
        {
          for(i=0;i<5;i++)
            printf("%u ",nn[i]);
          printf("... ");
          for(i=n-5;i<n;i++)
            printf("%u ",nn[i]);
        }
      else
        {
          for(i=0;i<n;i++)
            printf("%u ",nn[i]);
        }
    }
  else
    printf("no point");
  free(nn);
  printf("\n\n");

  /* If we have to collect statistics about the execution time.... */
  if (rep>1)
    {
      printf("Average time per NN query\n");
 
      #ifdef _MPNN
        tStart=getTime();
        for(j=0;j<rep;j++)
          {
            RandomPointInRectangle(q,&ambient);
            QueryKDTree(q,(int *)&id,&d,Ykdt);
          }
        tEnd=getTime();
        printf("  MPNN       : %f ms\n",(tEnd-tStart)/(double)rep);
      #endif

      tStart=getTime();
      for(j=0;j<rep;j++)
        {
          RandomPointInRectangle(q,&ambient);
          id=NearestNeighbour(topology,q,Mkdt);
        }
      tEnd=getTime();
      printf("  Cuik       : %f ms\n",(tEnd-tStart)/(double)rep);

      tStart=getTime();
      for(j=0;j<rep;j++)
        {
          dMin=INF;
          for(i=0;i<nPoints;i++)
            {
              d=VectorSquaredDistanceTopologyMin(dMin,dim,topology,point[i],q);
              if (d<dMin)
                {
                  idMin=i;
                  dMin=d;
                }
            }
        }
      tEnd=getTime();
      printf("  Brute-force: %f ms\n\n",(tEnd-tStart)/(double)rep);

      printf("Average time per in-ball query\n");

      #ifdef _MPNN
        tStart=getTime();
        for(j=0;j<rep;j++)
          {
            RandomPointInRectangle(q,&ambient);
            QueryKDTreeR(q,searchR,(int *)(&n),(int **)(&nn),&dst,Ykdt);
            if (n>0)
              {
                free(nn);
                free(dst);
              }
          }
        tEnd=getTime();
        printf("  MPNN       : %f ms\n",(tEnd-tStart)/(double)rep);
      #endif

      tStart=getTime();
      for(j=0;j<rep;j++)
        {
          RandomPointInRectangle(q,&ambient);
          NeighboursInBall(topology,q,searchR,&n,&nn,Mkdt);
          if (n>0)
            free(nn);
        }
      tEnd=getTime();
      printf("  Cuik       : %f ms\n",(tEnd-tStart)/(double)rep);

      tStart=getTime();
      for(j=0;j<rep;j++)
        {
          RandomPointInRectangle(q,&ambient);
          n=0;
          m=100;
          NEW(nn,m,unsigned int);
          for(i=0;i<nPoints;i++)
            {
              if (VectorSquaredDistanceTopologyMin(r2,dim,topology,point[i],q)<r2)
                {
                  if (n==m)
                    MEM_DUP(nn,m,unsigned int);
                  nn[n]=i;
                  n++;
                }
            }
          free(nn);
        }
      tEnd=getTime();
      printf("  Brute force: %f ms\n\n",(tEnd-tStart)/(double)rep);
    }

  /* Release the information */
  free(q);

  for(i=0;i<nPoints;i++)
    free(point[i]);
  free(point);

  if (topology!=NULL)
    free(topology);
  DeleteRectangle(&ambient);

  DeleteKDtree(Mkdt);
  #ifdef _MPNN
    DeleteKDTree(Ykdt);
  #endif

  return(EXIT_SUCCESS);
}