pkdfnn.10.10.01.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <math.h>
#include "tsa.h"
#include "tdkdtree.h"
#include "mpi.h"

// STATIC FUNCTIONS DECLARATION.
static void show_options();
static void scan_options();
static double standarddeviation();
static void makedelayvectors();
static FILE* prgerr = NULL;
static FILE* prgout = NULL;
static unsigned long N = ULONG_MAX;
static int med = 1;
static int m = 5;
static int tau = 1;
static double R_tol = 10.0; 
static unsigned long exclude  = 0;
static unsigned int column = 1;
static char* outfile = NULL;
static char* infile = NULL;
static unsigned int verbosity = 0xff;
static double* x = NULL;

int main
(
   int      argc  ,
   char**   argv
)
{
   double      R_A   =  0.0   ; /* attractor size */
   double      A_tol =  0.0   ; /* attractor size tolerance */
   int         n     =  0     ; /* number of delay vectors */
   double**    y     =  NULL  ; /* delay vectors */
   int         d     =  0     ; /* dimensions iterator */
   int         i     =  0     ; /* auxiliar variable */
   kdtree*     tree  =  NULL  ; /* the kd-tree */
   int         j     =  0     ; /* auxiliar variable */
   double      R_i   =  0.0   ; /* normalized distance */
   double      r_i   =  0.0   ; /* n-th coordinates distance */
   int*        fnn   =  NULL  ; /* false nearest neighbors found */
   int*        gfnn  =  NULL  ; /* global false nearest neighbors found */
   int         p     =  0     ; /* number of MPI processes */
   int         q     =  0     ; /* identifier of local process */
   double      t_s   =  0.0   ;
   double      t_f   =  0.0   ;

   MPI_Init(&argc,&argv);
   MPI_Comm_rank(MPI_COMM_WORLD,&q);
   MPI_Comm_size(MPI_COMM_WORLD,&p);
   prgerr = stderr;
   prgout = stdout;
   if( scan_help(argc,argv) )
   {
      if( q == 0 )
         show_options(argv[0]);
      MPI_Finalize();
      
      return 0;      
   }
   scan_options(argc,argv);
   infile = search_datafile(argc,argv,&column,verbosity);
   if( !infile )
   {
      if( q == 0 )
         fprintf(prgerr,"Input file not found. Exiting!\n");
      MPI_Finalize();
      
      return 0;
   }
   if( q == 0 && outfile )
   {
      prgout = fopen(outfile,"w");
      fprintf(prgerr,"Opened %s for writing\n",outfile);
   }
   x = get_series(infile,&N,exclude,column,verbosity);
// BEGIN FNN METHOD.
   if( q == 0 ) t_s = MPI_Wtime();
   R_A = standarddeviation(x,N);
   A_tol = 2.0;
   check_alloc(fnn = (int*)malloc(sizeof(int)*(m+1)));
   if( q == 0 )
   {
      check_alloc(gfnn = (int*)malloc(sizeof(int)*(m+1)));
      for( d = med ; d <= m ; d++ )
         gfnn[d] = 0;
   }
   for( d = med ; d <= m ; d++ )
   {
      fnn[d] = 0;
      n = N - d*tau;
      check_alloc( y = (double**)malloc(sizeof(double*)*d) );
      for( i = 0 ; i < d ; i++ )
         check_alloc( y[i] = (double*)malloc(sizeof(double)*n) );
      makedelayvectors(n,d,y);
      tree = newtdkdtree(y,n,d,p,q);
      // Top-down searching.
      for( i = 0 ; i < tree->n ; i++ )
      {
         j = tdnn(tree->perm[i],&R_i);
         r_i = fabs(x[tree->perm[i]+d*tau]-x[j+d*tau]);
         R_i = sqrt(R_i);
         if( r_i / R_i > R_tol || r_i / R_A > A_tol )
            fnn[d]++;
      }
      // Free local memory.
      if( tree ) destroytdkdtree(tree);
      if( y )
      {
         for( i = 0 ; i < d ; i++ )
            free(y[i]);
         free(y);
      }
   }
   MPI_Reduce(fnn,gfnn,m+1,MPI_INT,
         MPI_SUM,0,MPI_COMM_WORLD);
   if( q == 0 )
   {
      for( d = med ; d <= m ; d++ )
         fprintf(prgout,"%12d %10.3f %12d %12d\n",
            d,(double)gfnn[d]/(N - d*tau),gfnn[d],N - d*tau);

      t_f = MPI_Wtime();
      fprintf(prgerr,"RUNTIME %-.3lf\n",t_f-t_s);
   }
// END FNN METHOD.
   // Free global memory.
   if( q == 0 && gfnn ) free(gfnn);
   if( fnn ) free(fnn);
   if( q == 0 && outfile )
      fclose(prgout);
   if( outfile ) free(outfile);
   if( infile ) free(infile);
   if( x ) free(x);
   MPI_Finalize();

   return 0;
}

static double standarddeviation
(
   double*  x  ,
   int      N
)
{
   int      i           = 0   ;
   double   media       = 0.0 ;
   double   deviation   = 0.0 ;
   double   temp        = 0.0 ;
   double   R_A         = 0.0 ;

   media = 0.0;
   for ( i = 0 ; i < N ; i++ )
      media += x[i];
   media /= N;
   deviation = 0.0;
   for ( i = 0 ; i < N ; i++ )
   {
      temp = x[i] - media;
      deviation += temp * temp;
   }
   R_A = sqrt( deviation / N );
   
   return R_A; 
}

static void makedelayvectors
(
   int      n     ,
   int      d     ,
   double** y  
)
{
   int   i  =  0  ;
   int   j  =  0  ;
   
   for( i = 0 ; i < n ; i++ )
      for( j = 0 ; j < d ; j++ )
         y[j][i] = x[ i + j * tau ];

   return;
}

static void scan_options
(
   int      argc  ,
   char**   argv
)
{
   char* arg   =  NULL  ; // argument

   if( ( arg = check_option(argv,argc,'l','u') ) != NULL )
   {
      sscanf(arg,"%lu",&N);
      free(arg);
   }
   if( ( arg = check_option(argv,argc,'x','u') ) != NULL )
   {
      sscanf(arg,"%lu",&exclude);
      free(arg);
   }
   if( ( arg = check_option(argv,argc,'c','u') ) != NULL )
   {
      sscanf(arg,"%u",&column);
      free(arg);
   }
   if( ( arg = check_option(argv,argc,'m','u') ) != NULL )
   {
      sscanf(arg,"%d",&med);
      free(arg);
   }
   if( ( arg = check_option(argv,argc,'M','u') ) != NULL )
   {
      sscanf(arg,"%d",&m);
      free(arg);
   }
   if( ( arg = check_option(argv,argc,'d','u') ) != NULL )
   {
      sscanf(arg,"%d",&tau);
      free(arg);
   }
   if( ( arg = check_option(argv,argc,'f','f') ) != NULL )
   {
      sscanf(arg,"%lf",&R_tol);
      free(arg);
   }
   if( ( arg = check_option(argv,argc,'o','o') ) != NULL )
   {
      if( strlen(arg) > 0 )
         outfile = arg;
   }     
}


#define WID_STR "Computing minimal embedding dimension using a\n" \
   "        false nearest neighbors method by Kennel et al. (1991).\n" \
   "        This program is a kd-tree data structure approach\n" \
   "        based on a top-down search algorithm."

static void show_options
(
   char* progname
)
{
   what_i_do(progname,WID_STR);
   fprintf(prgerr,"Usage: %s [options]\n",progname);
   fprintf(prgerr,"Options:\n");
   fprintf(prgerr,"\t-l # of data [default: whole file]\n");
   fprintf(prgerr,"\t-x # of lines to be ignored [default: 0]\n");
   fprintf(prgerr,"\t-c column to read [default: 1]\n");
   fprintf(stderr,"\t-m minimal embedding dimension [default: 1]\n");
   fprintf(stderr,"\t-M maximal embedding dimension [default: 5]\n");
   fprintf(prgerr,"\t-d delay [default: 1]\n");
   fprintf(stderr,"\t-f tolerance factor [default: 10.0]\n");
   fprintf(prgerr,"\t-o outfile [default: standard output]\n");
   fprintf(prgerr,"\t-h show these options\n");
   fprintf(prgerr,"\n");
   fprintf(prgerr,
      "Everything not being a valid option will be interpreted as a possible\n"
      "datafile. If no datafile is given, the program %s is halted.\n"
      "The file reading interface code was copied from TISEAN 3.0.0 (C)\n"
      "by R. Hegger, H. Kantz, T. Schreiber (1998-2007)\n\n",
      progname);  
}

---