从时间序列中作庞卡来截面源程序.cpp

计算互信息和作庞卡来截面的C程序

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

poin.c -- started 8-26-96 by Eric R. Weeks

email: [e-mail]
web: http://www.physics.emory.edu/~weeks/software/poincare.html
home page: http://www.physics.emory.edu/~weeks/

This program is public domain, but please leave my name &
contact information intact.  I make no guarantees about bugs,
but if you find a problem with this program and contact me,
I may fix it.  Currently there is no guarantee this program
will work for anything but 3-D.

----------------------------------------

Why reinvent the wheel?  There's already a program called
"ponsec" on my computer to do this same thing.  Well, the man
page for ponsec does not match the current version of ponsec,
and I can't quite figure out how to make it work.  So, I'm
writing my own.

v01:  8-26-96: started from shell.c
      3-30-97: wow, discovered I had already started writing this program
               last summer.  Let's see what it does.  Added -o option
			added -abcxyzd options
v02:  4-11-97: add in ability to read in 3-D data (-M)
      4-11-97: add in uvec, vvec (-V); added -O
v03:  4-15-97: modify command line options (replace -abc with -n,
               replace -xyz with -x).  Add -T
      6-16-97: modify -P output slightly

 *******************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.14159265358979323846264338328
#define EE 2.71828182845904523536
#define MAX 200000

/* limit of 10-dimensional vectors */
#define MAXDIM 10

char *prgname;

main(argc,argv)
int argc;
char **argv;
{
	int c,numpts,delay,embed;
	double x,y,z,w;
	FILE *fp;
	extern int optind;
	extern char *optarg;
	float data[MAX][MAXDIM];
	int t,r,n;
	float coord[MAXDIM],lastcoord[MAXDIM];
	int side,lastside,timeseriesonly,psdraw;
	float vx,vy,vz;
	float qvec[MAXDIM],proj,proj1,proj2;
	float pvec[MAXDIM],nvec[MAXDIM];
	int direction,count1=0,count2=0,indim,usevecs,oneside=0;
	float uvec[MAXDIM],vvec[MAXDIM];
	float strobeperiod=-1.0,tt=0.0;

	delay = 1;
	embed = 3;
	indim = 1;
	timeseriesonly = 0;
	psdraw = 0;
	usevecs = 0;
	pvec[0] = 0.0; pvec[1] = 0.0; pvec[2] = 0.0;
	nvec[0] = 0.0; nvec[1] = 0.0; nvec[2] = 1.0;
	direction = 0;
	while ((c = getopt(argc, argv, "ht:m:oPnxd:M:VOT:")) != EOF)
	switch (c)  {
	case 'h': fprintf(stderr,"Usage: %s [options] [<] [file]\n\n",argv[0]);
fprintf(stderr," -h   : this help message\n");
fprintf(stderr,"INPUT:\n");
fprintf(stderr," -t # : delay in samples [%d]\n",delay);
fprintf(stderr," -m # : embedding dimension [%d]\n",embed);
fprintf(stderr," -M # : input data is this dimension -");
	fprintf(stderr," don't embed [%d]\n",indim);
fprintf(stderr,"POINCARE INFORMATION:\n");
fprintf(stderr," -n # # # : normal to poincare plane:");
	fprintf(stderr," [(%.1f, %.1f, %.1f)]\n",nvec[0],nvec[1],nvec[2]);
fprintf(stderr," -x # # # : point plane passes through:");
	fprintf(stderr," [(%.1f, %.1f, %.1f)]\n",pvec[0],pvec[1],pvec[2]);
fprintf(stderr," -T # : print out points with this strobing period");
	fprintf(stderr," [default:plane]\n");
fprintf(stderr," -d # : direction (0=norm, 1=-norm, 2=both) [%d]\n",direction);
fprintf(stderr," -V   : output coordinates in Poincare plane");
	fprintf(stderr," [default: 3D pts]\n");
fprintf(stderr,"OUTPUT:\n");
fprintf(stderr," -o   : output time series only (with delay coords)\n");
fprintf(stderr," -O   : print data points on one side of plane\n");
fprintf(stderr," -P   : make output psdraw compatible:\n");
fprintf(stderr,"          poin -P data | psdraw -C -c 0.1 > pic.ps\n");
fprintf(stderr,"\nInput data must be in one column only unless -M used!\n");
		   exit(1);
		   break;
	case 't': delay = atoi(optarg);
			break;
	case 'm': embed = atoi(optarg);
			break;
	case 'o': timeseriesonly = 1;
			break;
	case 'P': psdraw = 1;
			break;
	case 'x': pvec[0] = atof(argv[optind]);
			optind++;
			pvec[1] = atof(argv[optind]);
			optind++;
			pvec[2] = atof(argv[optind]);
			optind++;
			break;
	case 'n': nvec[0] = atof(argv[optind]);
			optind++;
			nvec[1] = atof(argv[optind]);
			optind++;
			nvec[2] = atof(argv[optind]);
			optind++;
			break;
	case 'd': direction = atoi(optarg);
			if ((direction<0) || (direction>2))  {
				fprintf(stderr,"error: -d option can be 0, 1, or 2\n");
				fprintf(stderr,"only.  Setting to default (0)\n");
				direction = 0;
			}
			break;
	case 'M': indim = atoi(optarg);
			delay = 0;
			embed = indim;
			break;
	case 'V': usevecs = 1;
			break;
	case 'O': oneside = 1;
			break;
	case 'T': strobeperiod = atof(optarg);
			break;
        }
	if (embed>3)  {
		fprintf(stderr,"Currently program not designed for\n");
		fprintf(stderr,"> 3 dimensions.  Trying anyway...\n");
	}
	if (embed>MAXDIM)  {
		fprintf(stderr,"Sorry, dimension must be less than\n");
		fprintf(stderr,"%d.  Change MAXDIM in program to fix.\n",MAXDIM);
		exit(1);
	}

	argc -= (optind-1) ; argv += (optind-1) ;
	fp = (argc > 1) ? fopen(*++argv, "r") : stdin;

	numpts = 0;
	switch (indim)  {
	case 1:
		while ((fscanf(fp,"%lf",&x)) == 1)  {
			data[numpts][0] = x;
			++numpts;
		}
		break;
	case 2:
		while ((fscanf(fp,"%lf %lf",&x,&y)) == 2)  {
			data[numpts][0] = x;
			data[numpts][1] = y;
			++numpts;
		}
		break;
	case 3:
		while ((fscanf(fp,"%lf %lf %lf",&x,&y,&z)) == 3)  {
			data[numpts][0] = x;
			data[numpts][1] = y;
			data[numpts][2] = z;
			++numpts;
		}
		break;
	case 4:
		while ((fscanf(fp,"%lf %lf %lf %lf",&x,&y,&z,&w)) == 4)  {
			data[numpts][0] = x;
			data[numpts][1] = y;
			data[numpts][2] = z;
			data[numpts][3] = w;
			++numpts;
		}
		break;
	default:
		fprintf(stderr,"I don't know how to read in more than 4\n");
		fprintf(stderr,"dimensions.  Exiting...\n");
		exit(1);
		break;
	}

	if (numpts<=delay*embed)  {
		fprintf(stderr,"ERROR: not enough points for delay coords\n");
		fprintf(stderr,"ERROR: numpts = %d\n",numpts);
		fprintf(stderr,"ERROR: delay  = %d\n",delay );
		exit(1);
	}

	if (timeseriesonly==1)  {
		if (indim>1)  {
			fprintf(stderr,"not quite sure what to do: -M option used\n");
			fprintf(stderr,"with -o option.  Taking my best guess.\n");
		}
		for (t=0;t<numpts-delay*embed;t++)  {
			for (r=0;r<embed;r++)
				for (n=0;n<indim;n++)
					printf("%f ",data[t+r*delay][n]);
			printf("\n");
		}
		exit(0);
	}

	if (usevecs==1)  {
		/* find uvec, vvec:  unit vectors, perpendicular to nvec */
		x = sqrt(nvec[0]*nvec[0]+nvec[1]*nvec[1]+nvec[2]*nvec[2]);
		if (x==0.0)  {
			fprintf(stderr,"error! normal vector has zero length\n");
			exit(1);
		}
		for (r=0;r<3;r++)  nvec[r] /= x;
		if (fabs(nvec[0]) > 0.01)  {
			uvec[0] = -nvec[1]/nvec[0];
			uvec[1] = 1.0;
			uvec[2] = 0.0;
		} else if (fabs(nvec[1]) > 0.01) {
			uvec[0] = 1.0;
			uvec[1] = -nvec[0]/nvec[1];
			uvec[2] = 0.0;
		} else {
			uvec[0] = 1.0;
			uvec[1] = 0.0;
			uvec[2] = -nvec[0]/nvec[2];
		}
		x = sqrt(uvec[0]*uvec[0]+uvec[1]*uvec[1]+uvec[2]*uvec[2]);
		for (r=0;r<3;r++)  uvec[r] /= x;
		vvec[0] = nvec[1]*uvec[2] - nvec[2]*uvec[1];
		vvec[1] = nvec[2]*uvec[0] - nvec[0]*uvec[2];
		vvec[2] = nvec[0]*uvec[1] - nvec[1]*uvec[0];
		fprintf(stderr,"nvec: %f %f %f\n",nvec[0],nvec[1],nvec[2]);
		fprintf(stderr,"uvec: %f %f %f\n",uvec[0],uvec[1],uvec[2]);
		fprintf(stderr,"vvec: %f %f %f\n",vvec[0],vvec[1],vvec[2]);
	}


	for (t=0;t<numpts-delay*embed;t++)  {
		if (indim==1)
			for (r=0;r<embed;r++)  coord[r] = data[t+r*delay][0];
		else
			for (r=0;r<indim;r++)  coord[r] = data[t][r];
		/* now coord[] contains the embed-dimensional data pt */

		if (strobeperiod<0.0)  {
			proj = 0.0;
			for (r=0;r<embed;r++)
				proj += nvec[r]*(coord[r] - pvec[r]);
			if (proj<0.0)  side = 1; else side = 2;

			if ((oneside==1)&&(side==2-direction))  {
				for (r=0;r<embed;r++)  printf("%f ",coord[r]);
				printf("\n");
			} else {
				if ((t>0)&&(  ((side>lastside)&&(direction != 1)) ||
						((side<lastside)&&(direction != 0)) ))  {
					/* crossing occurred */
					if ((side>lastside)&&(direction != 1))
						count1++;
					if ((side<lastside)&&(direction != 0))
						count2++;
					proj2 = 0.0;
					for (r=0;r<embed;r++)
						proj2 += (coord[r]-lastcoord[r])*nvec[r];
					proj /= proj2;
					if (usevecs == 0)  {
						for (r=0;r<embed;r++)  {
						  qvec[r]=proj*lastcoord[r]+(1.0-proj)*coord[r];
						  printf("%f ",qvec[r]);
						}
					} else {
						x = 0.0; y = 0.0;
						for (r=0;r<embed;r++)  {
						  qvec[r]=proj*lastcoord[r]+(1.0-proj)*coord[r];
						  x += qvec[r]*uvec[r];
						  y += qvec[r]*vvec[r];
						}
						printf("%f %f ",x,y);
					}
					if (psdraw==0)
						printf("\n");
					else
						if (side>lastside)
							printf("1 0 0.2\n");
						else
							printf("0 0.2 1\n");
				}
			}

			lastside = side;
		} else {
			/* strobeperiod > 0.0 */
			tt+=1.0;
			if (tt>strobeperiod)  {
				count1++;
				count2++;
				tt -= strobeperiod;
				/* now tt is between 0 and 1 */
				for (r=0;r<embed;r++)
					qvec[r] = coord[r]*tt+lastcoord[r]*(1.0-tt);
				if (usevecs == 0)  {
					for (r=0;r<embed;r++)  printf("%f ",qvec[r]);
				} else {
					x = 0.0; y = 0.0;
					for (r=0;r<embed;r++)  {
					  x += qvec[r]*uvec[r];
					  y += qvec[r]*vvec[r];
					}
					printf("%f %f ",x,y);
				}
				printf("\n");
			}
		}
		for (r=0;r<embed;r++)  lastcoord[r] = coord[r];
	}
	if (direction != 1)
fprintf(stderr,"%d crossing in direction of normal (red)\n",count1);
	if (direction != 0)
fprintf(stderr,"%d crossing in opposite direction of normal (blue)\n",count2);

}