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>
#ifdef _MPNN
  #include <DNN/ANN_C.h>
#endif

//#define FIXED_SET

#ifdef FIXED_SET
double _point[38][8]={
  {-1.262280,2.039516,-1.241662,-1.121883,1.464420,-2.198507,1.457702,-1.109374 },
  {-1.237220,1.936589,-1.034433,-1.362698,1.509152,-2.025385,1.470437,-1.237948 },
  {-1.013849,1.950820,-1.295149,-1.183568,1.455763,-2.055501,1.515250,-1.356307 },
  {-1.482196,2.007648,-0.996732,-1.283119,1.534123,-2.183768,1.385654,-0.987308 },
  {-1.500418,2.097868,-1.220372,-1.033085,1.478405,-2.344629,1.351986,-0.864703 },
  {-1.058694,2.038497,-1.494067,-0.945908,1.372108,-2.234390,1.524408,-1.210050 },
  {-0.798914,1.929518,-1.574707,-0.976722,1.373544,-2.116169,1.558028,-1.453671 },
  {-0.945818,2.037145,-1.752772,-0.721324,1.236064,-2.330143,1.579138,-1.212565 },
  {-1.271433,2.116328,-1.528064,-0.812196,1.331857,-2.393527,1.459837,-0.966149 },
  {-1.741267,2.010276,-0.887950,-1.279929,1.584056,-2.287866,1.257626,-0.759013 },
  {-1.501615,1.884796,-0.762412,-1.528952,1.543581,-2.035084,1.412669,-1.072706 },
  {-1.264555,1.811661,-0.803575,-1.606490,1.510936,-1.867359,1.477011,-1.320378 },
  {-1.003100,1.834422,-1.071949,-1.441231,1.496033,-1.862884,1.501577,-1.483956 },
  {-1.645768,2.140992,-1.346015,-0.834454,1.431947,-2.503390,1.247494,-0.654193 },
  {-0.772193,1.827176,-1.340938,-1.248031,1.463509,-1.909529,1.518583,-1.600614 },
  {-1.170367,2.110745,-1.780438,-0.590770,1.168126,-2.476101,1.524563,-0.960380 },
  {-1.443594,2.150651,-1.626567,-0.630068,1.254252,-2.545875,1.376031,-0.740363 },
  {-0.772232,1.705545,-1.111414,-1.513798,1.491066,-1.702485,1.482258,-1.724730 },
  {-1.019783,1.707103,-0.855377,-1.681259,1.484966,-1.688372,1.489860,-1.578979 },
  {-0.534333,1.802012,-1.741788,-0.900786,1.350070,-2.073238,1.555067,-1.650844 },
  {-1.760679,1.875779,-0.639120,-1.535263,1.580386,-2.148980,1.310228,-0.841421 },
  {-1.495381,1.733306,-0.544721,-1.760710,1.500579,-1.881177,1.453162,-1.170318 },
  {-1.810300,2.087329,-1.079026,-1.044453,1.559473,-2.440484,1.162465,-0.600172 },
  {-1.984806,1.987968,-0.813940,-1.252543,1.641118,-2.384327,1.089403,-0.526641 },
  {-1.697993,2.147238,-1.551797,-0.593581,1.306569,-2.632096,1.192044,-0.503621 },
  {-1.880391,2.122402,-1.298290,-0.781258,1.482256,-2.588584,1.056138,-0.421648 },
  {-0.928724,2.033960,-1.978375,-0.483489,1.057093,-2.442046,1.627235,-1.132694 },
  {-0.681670,1.919328,-1.901170,-0.672762,1.202390,-2.265498,1.618211,-1.428134 },
  {-1.256384,1.658986,-0.596002,-1.835289,1.462431,-1.699324,1.492230,-1.425470 },
  {-0.543469,1.673458,-1.361725,-1.324151,1.484989,-1.760384,1.469352,-1.832039 },
  {-1.363260,2.127833,-1.865133,-0.407814,1.059721,-2.613543,1.444583,-0.719095 },
  {-2.039651,2.055739,-1.027143,-1.002035,1.620787,-2.517956,0.964684,-0.367795 },
  {-2.007873,1.841635,-0.545726,-1.522485,1.626393,-2.256107,1.166803,-0.605743 },
  {-1.780364,1.712648,-0.399743,-1.762945,1.524117,-2.026988,1.367265,-0.902002 },
  {-1.028464,1.558295,-0.645619,-1.909291,1.426565,-1.515046,1.478929,-1.673280 },
  {-0.795581,1.580700,-0.896500,-1.754151,1.464109,-1.517622,1.451276,-1.816512 },
  {-1.643063,1.578628,-0.275135,-1.936803,1.435267,-1.851053,1.461556,-1.088357 },
  {-0.538870,1.554160,-1.151812,-1.573927,1.493638,-1.553494,1.414819,-1.951598 }
};

double _q[8]={-0.419985,1.755759,-1.993156,-0.661874,1.224372,-2.184798,1.589697,-1.651615 };
#endif

#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++;
                  expansion=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);

  #ifdef FIXED_SET
    tpg=TOPOLOGY_R;
    nPoints=38;
    dim=8;
    searchR=0.4*2;
  #endif

  /* 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-(POINTS_ON_SPHERE?0.1:0);
      u=+1+(POINTS_ON_SPHERE?0.1:0);
      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
        #ifdef FIXED_SET
          for(j=0;j<dim;j++)
            point[i][j]=_point[i][j];
        #else
          RandomPointInRectangle(point[i],&ambient);
        #endif
      #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,topology,pLeaf,expansion,1,&(point[i]),&i);
      else
        AddPoint2KDtree(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))); 

  //PrintKDtree(stdout,1,Mkdt);

  #if (PLOTKDT)
    if  (dim==3)
      PlotKDtree("kdtree.gcl",POINTS_ON_SPHERE,Mkdt);
  #endif

  /* Define the query point */
  printf("Defining the query point\n\n");
  NEW(q,dim,double);
  #ifdef FIXED_SET
    for(j=0;j<dim;j++)
      q[j]=_q[j];
  #else
    RandomPointInRectangle(q,&ambient);
  #endif
  
  /* 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(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(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(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(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);
}