tdkdtree.c

Go to the documentation of this file.

#include "tdkdtree.h"
#include "queue.h"
#include <stdio.h>
#include <stdlib.h>

// STATIC FUNCTIONS DECLARATION.
static void globaltreebuild(); 
static infokdnode* infokdnodebuild(); 
static kdtree* localtreebuild(); 
static kdnode* tdbuild(); 
static kdnode* tdfindcutvalues(); 
static int arrange(); 
static void rnn(); 
static double distsqrd2(); 
static void destroytdkdtreenode(); 
static double** y = NULL;

static int n = 0;

static int k = 0;

static int* perm = NULL;

static kdnode* root = NULL;

static int cutoff = 15;

static int bigbucket = 0;

kdtree* newtdkdtree
(
   double** param_y     ,
   int      param_n     ,
   int      param_k     ,
   int      p           ,
   int      q
)
{
   kdtree*  tree  =  NULL  ;
   int      i     =  0     ;

   tree = (kdtree*)malloc(sizeof(kdtree));
   tree->y = y = param_y;
   tree->n = n = param_n;
   tree->k = k = param_k;
   perm = (int*)malloc(sizeof(int)*n);
   for( i = 0 ; i < n ; i++ ) perm[i] = i;
   tree->perm = perm;
   if( p == 1 ) // Monoprocessing.
      tree->root = root = tdbuild(0,n-1);
   else // Multiprocessing.
   {
      globaltreebuild(p);
      tree = localtreebuild(tree,p,q);
   }

   return tree;   
}

static void globaltreebuild
(
   int      p
)
{
   qnode*   node     =  NULL  ;
   qnode*   lonode   =  NULL  ;
   qnode*   hinode   =  NULL  ;
   int      i        =  0     ;

   qbuild();
   node = (qnode*)malloc(sizeof(qnode));
   node->info = infokdnodebuild(0,n-1);
   qinsert(node);
   for( i = 2 ; i <= p ; i++ )
   {
      node = qretrieve();
      lonode = (qnode*)malloc(sizeof(qnode));
      lonode->info = infokdnodebuild(node->info->lopt_l,node->info->hipt_l);
      hinode = (qnode*)malloc(sizeof(qnode));
      hinode->info = infokdnodebuild(node->info->lopt_h,node->info->hipt_h);
      qinsert(lonode);
      qinsert(hinode);
      qdelete();
   }

   return;
}

static infokdnode* infokdnodebuild
(
   int   l  ,
   int   u
)
{
   infokdnode* info  =  NULL  ;
   kdnode*     p     =  NULL  ;
   int         m     =  0     ;

   info = (infokdnode*)malloc(sizeof(infokdnode));
   p = (kdnode*)malloc(sizeof(kdnode));
   p = tdfindcutvalues(l,u,p);
   info->cutdim = p->cutdim;
   info->cutval = p->cutval;
   m = arrange(l,u,info->cutval,info->cutdim);
   info->lopt_l = l;
   info->hipt_l = m;
   info->lopt_h = m+1;
   info->hipt_h = u;

   return info;
}

static kdtree* localtreebuild
(
   kdtree*  tree  ,
   int      p     ,
   int      q
)
{
   int      id       =  0     ;
   qnode*   node     =  NULL  ;
   int      i        =  0     ;
   int      j        =  0     ;
   int*     newperm  =  NULL  ;

   for( id = 0 ; id < p ; id++ )
   {
      node = qretrieve();
      if( id == q )
      {
         tree->n = n = node->info->hipt_h - node->info->lopt_l + 1;
         newperm = (int*)malloc(sizeof(int)*n);
         j = 0;
         for( i = node->info->lopt_l ; i <= node->info->hipt_h ; i++ )
         {
            newperm[j] = perm[i];
            j++;
         }
         if( perm ) free(perm);
         tree->perm = perm = newperm;
         tree->root = root = tdbuild(0,n-1);
         break;
      }
      qdelete();
   }
   qdestroy();

   return tree;
}

static kdnode* tdbuild
(
   int   l  ,
   int   u
)
{
   kdnode*  p  =  NULL  ;
   int      m  =  0     ;
   
   p = (kdnode*)malloc(sizeof(kdnode));
   if( u - l + 1 <= cutoff )
   {
      p->bucket = 1;
      p->lopt = l;
      p->hipt = u;
      p->loson = p->hison = NULL;
   }
   else
   {
      p->bucket = 0;
      p = tdfindcutvalues(l,u,p);
      m = arrange(l,u,p->cutval,p->cutdim);
      if( bigbucket )
      {
         p->bucket = 1;
         p->lopt = l;
         p->hipt = u;
         p->loson = p->hison = NULL;
         
         return p;
      }
      p->loson = tdbuild(l,m);
      p->hison = tdbuild(m+1,u);
   }
   
   return p;
}

static kdnode* tdfindcutvalues
(
   int      l  ,
   int      u  ,
   kdnode*  p
)
{
   int      i           =  0     ;
   int      j           =  0     ;
   double   sum         =  0.0   ;
   double   mean        =  0.0   ;
   double   temp        =  0.0   ;
   double   deviation   =  0.0   ;
   double   maxspread   =  -HUGE ;

   if( u - l + 1 > 1000 )
      u = l + 1000 - 1; 
   for( i = 0 ; i < k ; i++)
   {
      sum = 0.0;
      for( j = l ; j <= u ; j++)
         sum += y[i][perm[j]];
      mean = sum / (u - l + 1);
      sum = 0.0;
      for( j = l ; j <= u ; j++)
      {
         temp = y[i][perm[j]] - mean;
         if( temp < 0.0 )
            temp *= -1;
         sum += temp;
      }
      deviation = sum / (u - l + 1);
      if(deviation > maxspread )
      {
         p->cutdim = i;
         p->cutval = mean;
         maxspread = deviation; 
      }
   }
   
   return p;
}

static int arrange
(
   int      l        ,
   int      u        ,
   double   cutval   ,
   int      cutdim
)
{
   int      n     =  0     ;
   int      n_l   =  0     ;
   int      n_r   =  0     ;
   int      i     =  0     ;
   double   value =  0.0   ;
   int      m     =  0     ;

   n = u - l + 1;
   do
   {
      i = perm[l + n_l];
      value = y[cutdim][i];
      if (value <= cutval )
         n_l++;
      else
      {
         perm[l + n_l] = perm[u - n_r];
         perm[u - n_r] = i;
         n_r++;
      } 
   }while( n_l + n_r != n );
   m = l + n_l - 1;
   bigbucket = 0;
   if( n > cutoff && ( n_l == 0 || n_r == 0 ) )
      bigbucket = 1;

   return m;
}

static int nntarget = 0;
static int nnfound = 0;
static double nndist2 = 0.0;

int tdnn
(
   int      j              ,
   double*  param_nndist2
)
{
   nntarget = j;
   nndist2 = HUGE;
   rnn(root);
   *param_nndist2 = nndist2;  

   return nnfound;
}

static void rnn
(
   kdnode*  p
)
{
   int      i        =  0     ;
   double   thisdist =  0.0   ;  
   double   diff     =  0.0   ;
   
   if( p->bucket )
   {
      for( i = p->lopt ; i <= p->hipt ; i++ )
      {
         if( perm[i] != nntarget )
         {
            thisdist = distsqrd2(perm[i],nntarget);
            if( thisdist < nndist2 )
            {
               nndist2 = thisdist;
               nnfound = perm[i];
            }
         }
      }
   }
   else
   {
      diff = y[p->cutdim][nntarget] - p->cutval;   
      if( diff < 0.0 )
      {
         rnn(p->loson);
         if( nndist2 > diff * diff )
            rnn(p->hison);
      }
      else
      {
         rnn(p->hison);
         if( nndist2 > diff * diff )
            rnn(p->loson);
      }
   }
}

static double distsqrd2
(
   int   i  ,
   int   j
)
{
   int      d     =  0     ;
   double   diff  =  0.0   ;
   double   dist2 =  0.0   ;
   
   for( d = 0 ; d < k ; d++ )
   {
      diff = y[d][i] - y[d][j];
      dist2 += diff * diff;
   }
   
   return dist2;
}

void destroytdkdtree
(
   kdtree*  tree
)
{
   if( tree )
   {
      if( tree->root )
         destroytdkdtreenode(tree->root);
      free(tree);
   }
   if( perm ) free(perm);
   y = NULL;
   n = 0;
   k = 0;
   perm = NULL;
   root = NULL;
}

static void destroytdkdtreenode
(
   kdnode*  node
)
{
   if( node )
   {
      if( node->loson ) destroytdkdtreenode(node->loson);
      if( node->hison ) destroytdkdtreenode(node->hison);
      free(node);
   }
}

// END OF FILE.