interpoint.c

代码用于估计关联维数。包括G-P算法(corrint.m)

/*=================================================================
 *
 * INTERPOINT.C	
 *              .MEX file corresponding to INTERPOINT.M
 *              returns the number of interpoint distance of the 
 *              embedded time series in each bin
 *
 * The calling syntax is:
 *
 *		np = interpoint(y,de,tau,bins,nref,nt)
 *
 *
 * This is a MEX-file for MATLAB.  
 *=================================================================*/
/* $Revision: 1.5 $ */
#include <math.h>
#include "mex.h"
#include <stdlib.h>
#include <math.h>


void ipb(double	*data, /* the embedded data (d-by-n) */
	 int ndata, /* n-points */
	 int *de, /* embedding dimensions  */
	 int nde, /* number of embedding dimensions */
	 int tau, /* embedding lag */
	 double *bins, /* the bins */
	 int nbins, /* the number of bins */
	 long int nref, /* number of reference vectors to choose */
	 int nt, /* the Thieler band (forbidden points) */
	 double *binpop) /* the population in each bin (to be computed) */
     /*actually calculate the interpoint binning*/
{
  long int i,idum;
  int xi,yi,k,p,q; /* overdoes on indices */
  int nextde,maxde,maxn;
  double dist,val;
  bool toosmall,useall;

  /*seed the RNG using the standard one*/
  idum=rand();

  /* determnie the maximum embedding window */
  maxde=*(de+nde-1);
  maxn=ndata-(maxde-1)*tau;
  
  /*initialise random number generator */
  srand((unsigned)time(NULL));

  /*use all the data or random reference vectors? */
  if (nref==0) {
    useall=1;
    if (nt>0) {
      nref=(maxn-2*nt+2);
      nref=nref*(maxn-2*nt+1);;
      nref += 2*(nt-1)*(maxn-nt+1)-(nt-1)*nt;
    } else {
      nref=maxn;
      nref=nref*nref;
    }
    xi=0;
    yi=0+nt-1;
    /*  nref=maxn*(maxn-nt); */
  } else {
    useall=0;
  }

  for (i=0; i<nref; i++) {

    /* get xi and yi */
    if (useall) {
      /* use all data and therefore just increment xi and yi */
      yi++;
      if (yi==maxn) {
	yi=0;
	xi++;
      }
      /* ensure that the points aren't too close  */     
      while (abs(xi-yi)<nt) {
	yi++;
	if (yi==maxn) {
	  yi=0;
	  xi++;
	}
      }
    } else {
      /* choose xi and yi randomly between 0 and n */
      /* this RNG is supposed to have a period of O(10^9) */
	/* xi and yi are in [1,maxn] */
      idum=1664525*idum+1013904223;
      xi=(int)((idum/2147483648.0+1)*maxn/2.0); 
      idum=1664525*idum+1013904223;
      yi=(int)((idum/2147483648.0+1)*maxn/2.0); 

      /* check that the points aren't too close  */     
      while (abs(xi-yi)<nt) {
	idum=1664525*idum+1013904223;
      	yi=(int)((idum/2147483648.0+1)*maxn/2.0); 
	}

      /*xi=(int)(rand()/(RAND_MAX+1.0)*maxn)+1;*/ 
      /* I know that the C/C++ rand functions are crap, but it'll do */
      /* yi=(int)(rand()/(RAND_MAX+1.0)*maxn)+1;*/
      /* find points that aren't too close  */     
      /* while (abs(xi-yi)<nt) {*/
      /*	yi=(int)(rand()/(RAND_MAX+1.0)*maxn)+1; */
      /*} */

    } /* if (useall) */
      
    /* compute the square of the distance for these points  */
    dist=0;
    nextde=*(de);
    p=0;
    q=0;
    k=0;
    for (k=0; k<maxde; k++) {
      val = *(data+xi+k*tau)-*(data+yi+k*tau);
      dist += val*val;
      if (k==nextde-1) { /*is this one of the de we should check */
	
	/* find the bin this distance fits in  */
	if ( (*(bins+p))>dist) {
	  (*(binpop+p+q*(nbins+1)))++; /*increment the relevant bin */
	} else {
	  toosmall=1;
	  while (toosmall && p<nbins) {
	    if ( (*(bins+p))>dist) {
	      toosmall=0;
	      (*(binpop+p+q*(nbins+1)))++; /*increment the relevant bin */
	    }
	    p++;
	  } /* while */
	  p--;
	}
	/*increment nextde*/
	q++;
	nextde=*(de+q);
	
      } /* if (k==nextde) */
	
    } /* for (k= ...) */
    
  }/* for i */
  
}


void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray *prhs[] )
     /* the MATLAB mex wrapper function */     
{ 
    double *data,*bins,*thearg,*binpop;
    int mrows,ncols,i,j;
    int d,n,nbins,tau,nt,nde;
    int *de;
    long int nref;
    bool warnThem;
    
    /* Check for proper number of arguments */
    
    if (nrhs > 6) { 
	mexErrMsgTxt("Too many input arguments."); 
    }
    if (nrhs < 4) {
        mexErrMsgTxt("Insufficient input arguments.");
    }
    if (nrhs < 5) {
      /* set nref=0 */
      nref=0;
    } 
    if (nrhs < 6) {
      /* set nt=0 */
      nt=0;
    } 
    if (nlhs > 1) {
	mexErrMsgTxt("Too many output arguments."); 
    } 
    
    /*Assign a pointer to the input matrix*/
    data = mxGetPr(prhs[0]);

    /* The first input argument (y) */     
    mrows = mxGetM(prhs[0]); 
    ncols = mxGetN(prhs[0]);
    /* check that y is a vector */
    if ((mrows!=1) && (ncols!=1)) {
      mexWarnMsgTxt("First input should be a vector");
    }
    n=mrows*ncols; /* number of data points */
    
    /* The second input argument (de) */
    thearg =mxGetPr(prhs[1]);
    mrows = mxGetM(prhs[1]);
    ncols = mxGetN(prhs[1]);
    /* check that all bins ascending */
    warnThem=0;
    nde=mrows*ncols;
    for (i=0; i<(nde-1); i++) 
      if ( *(thearg+i+1)<*(thearg+i)) 
	warnThem=1;
    if (warnThem) {
      mexErrMsgTxt("Error in interpoint: Second input not ascending");
    }
    /* alloc memory and make a copy of de (squaring as we go)*/
    /* and,  ... check that it is an integer */
    warnThem=0;
    de = (int *) calloc(nde,sizeof(int));
    for (i=0; i<nde; i++) {
      *(de+i) = (int)(*(thearg+i));
      if (floor(*(thearg+i))!=*(thearg+i)) {
	warnThem=1;
	}
    }
    if (warnThem) {
      mexWarnMsgTxt("Second input being rounded to integer values");
    }

    /* The third input argument (tau) */
    thearg = mxGetPr(prhs[2]);
    mrows = mxGetM(prhs[2]);
    ncols = mxGetN(prhs[2]);
    /* check that x is a scalar */
    if ((mrows!=1) || (ncols!=1)) {
      mexWarnMsgTxt("Third input should be a scalar");
    }
    /* check that it is an integer */
    if (floor(*(thearg))!=*(thearg)) {
      mexWarnMsgTxt("Third input being rounded to an integer value");
    }
    tau = (int)(*(thearg));
    
    /* The fourth input arg (bins) */
    thearg = mxGetPr(prhs[3]);
    mrows = mxGetM(prhs[3]);
    ncols = mxGetN(prhs[3]);
    /* check that x is a vector */
    if ((mrows!=1) && (ncols!=1)) {
      mexWarnMsgTxt("Fourth input should be a vector");
    }
    nbins = mrows*ncols;
    /* check that all bins ascending */
    warnThem=0;
    for (i=0; i<(nbins-1); i++) 
      if ( *(thearg+i+1)<*(thearg+i)) 
	warnThem=1;
    if (warnThem) {
      mexErrMsgTxt("Error in interpoint: Second input not ascending");
    }
    /*alloc memory and make a copy of bins (squaring as we go)*/
    bins = (double *) calloc(nbins,sizeof(double));
    for (i=0; i<nbins; i++) 
      *(bins+i) = (*(thearg+i))*(*(thearg+i));

    /* The fifth input argument (nref) */
    if (nrhs>4){
      thearg =mxGetPr(prhs[4]);
      mrows = mxGetM(prhs[4]);
      ncols = mxGetN(prhs[4]);
      /* check that x is a scalar */
      if ((mrows!=1) || (ncols!=1)) {
	mexWarnMsgTxt("Fifth input should be a scalar");
      }
      /* check that it is an integer */
      warnThem=0;
      if (floor(*(thearg))!=*(thearg)) {
	mexWarnMsgTxt("Fifth input being rounded to an integer value");
	}
      nref=(long int)(*(thearg));
    }

    /* The sixth input argument (nt) */
    if (nrhs==6){
      thearg =mxGetPr(prhs[5]);
      mrows = mxGetM(prhs[5]);
      ncols = mxGetN(prhs[5]);
      /* check that x is a scalar */
      if ((mrows!=1) || (ncols!=1)) {
	mexWarnMsgTxt("Sixth input should be a scalar");
      }
      /* check that it is an integer */
      warnThem=0;
      if (floor(*(thearg))!=*(thearg)) {
	mexWarnMsgTxt("Fourth input being rounded to an integer value");
	}
      nt=(int)(*(thearg));
    }    

    /* Create a matrix for the return argument */ 
    plhs[0] = mxCreateDoubleMatrix(nbins+1,nde, mxREAL); 

    /* Assign pointer to the ouput matrix */ 
    binpop = mxGetPr(plhs[0]);
        
    /* Do the actual computations in a subroutine */
    ipb(data,n,de,nde,tau,bins,nbins,nref,nt,binpop);

    return;
    
}