--- /dev/null
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * TISEAN is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/*Author: Rainer Hegger. Last modified: Sep 3, 1999*/
+#include <math.h>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "routines/tsa.h"
+
+#define WID_STR "Estimates the maximal Lyapunov exponent using the Kantz\n\t\
+algorithm"
+
+#define BOX 128
+const unsigned int ibox=BOX-1;
+
+unsigned long length=ULONG_MAX;
+unsigned long exclude=0;
+unsigned long reference=ULONG_MAX;
+unsigned int maxdim=2;
+unsigned int mindim=2;
+unsigned int delay=1;
+unsigned int column=1;
+unsigned int epscount=5;
+unsigned int maxiter=50;
+unsigned int window=0;
+unsigned int verbosity=0xff;
+double epsmin=1.e-3,epsmax=1.e-2;
+char eps0set=0,eps1set=0;
+char *outfile=NULL;
+char *infile=NULL;
+
+double *series,**lyap;
+long box[BOX][BOX],*liste,**lfound,*found,**count;
+double max,min;
+
+void show_options(char *progname)
+{
+ what_i_do(progname,WID_STR);
+
+ fprintf(stderr," Usage: %s [options]\n",progname);
+ fprintf(stderr," Options:\n");
+ fprintf(stderr,"Everything not being a valid option will be "
+ "interpreted as a possible datafile.\nIf no datafile "
+ "is given stdin is read. Just - also means stdin\n");
+ fprintf(stderr,"\t-l # of data [default: whole file]\n");
+ fprintf(stderr,"\t-x # of lines to be ignored [default: 0]\n");
+ fprintf(stderr,"\t-c column to read [default: 1]\n");
+ fprintf(stderr,"\t-M maxdim [default: 2]\n");
+ fprintf(stderr,"\t-m mindim [default: 2]\n");
+ fprintf(stderr,"\t-d delay [default: 1]\n");
+ fprintf(stderr,"\t-r mineps [default: (data interval)/1000]\n");
+ fprintf(stderr,"\t-R maxeps [default: (data interval)/100]\n");
+ fprintf(stderr,"\t-# # of eps [default: 5]\n");
+ fprintf(stderr,"\t-n # of reference points [default: # of data]\n");
+ fprintf(stderr,"\t-s # of iterations [default: 50]\n");
+ fprintf(stderr,"\t-t time window [default: 0]\n");
+ fprintf(stderr,"\t-o outfile [default: 'datafile'.lyap]\n");
+ fprintf(stderr,"\t-V verbosity level [default: 3]\n\t\t"
+ "0='only panic messages'\n\t\t"
+ "1='+ input/output messages'\n\t\t"
+ "2='+ plus statistics'\n");
+ fprintf(stderr,"\t-h show these options\n");
+ exit(0);
+}
+
+void scan_options(int n,char **str)
+{
+ char *out;
+
+ if ((out=check_option(str,n,'l','u')) != NULL)
+ sscanf(out,"%lu",&length);
+ if ((out=check_option(str,n,'x','u')) != NULL)
+ sscanf(out,"%lu",&exclude);
+ if ((out=check_option(str,n,'c','u')) != NULL)
+ sscanf(out,"%u",&column);
+ if ((out=check_option(str,n,'M','u')) != NULL)
+ sscanf(out,"%u",&maxdim);
+ if ((out=check_option(str,n,'m','u')) != NULL)
+ sscanf(out,"%u",&mindim);
+ if ((out=check_option(str,n,'d','u')) != NULL)
+ sscanf(out,"%u",&delay);
+ if ((out=check_option(str,n,'r','f')) != NULL) {
+ eps0set=1;
+ sscanf(out,"%lf",&epsmin);
+ }
+ if ((out=check_option(str,n,'R','f')) != NULL) {
+ eps1set=1;
+ sscanf(out,"%lf",&epsmax);
+ }
+ if ((out=check_option(str,n,'#','u')) != NULL)
+ sscanf(out,"%u",&epscount);
+ if ((out=check_option(str,n,'n','u')) != NULL)
+ sscanf(out,"%lu",&reference);
+ if ((out=check_option(str,n,'s','u')) != NULL)
+ sscanf(out,"%u",&maxiter);
+ if ((out=check_option(str,n,'t','u')) != NULL)
+ sscanf(out,"%u",&window);
+ if ((out=check_option(str,n,'V','u')) != NULL)
+ sscanf(out,"%u",&verbosity);
+ if ((out=check_option(str,n,'o','o')) != NULL)
+ if (strlen(out) > 0)
+ outfile=out;
+}
+
+void put_in_boxes(double eps)
+{
+ unsigned long i;
+ long j,k;
+ static unsigned long blength;
+
+ blength=length-(maxdim-1)*delay-maxiter;
+
+ for (i=0;i<BOX;i++)
+ for (j=0;j<BOX;j++)
+ box[i][j]= -1;
+
+ for (i=0;i<blength;i++) {
+ j=(long)(series[i]/eps)&ibox;
+ k=(long)(series[i+delay]/eps)&ibox;
+ liste[i]=box[j][k];
+ box[j][k]=i;
+ }
+}
+
+void lfind_neighbors(long act,double eps)
+{
+ unsigned int hi,k,k1;
+ long i,j,i1,i2,j1,element;
+ static long lwindow;
+ double dx,eps2=sqr(eps);
+
+ lwindow=(long)window;
+ for (hi=0;hi<maxdim-1;hi++)
+ found[hi]=0;
+ i=(long)(series[act]/eps)&ibox;
+ j=(long)(series[act+delay]/eps)&ibox;
+ for (i1=i-1;i1<=i+1;i1++) {
+ i2=i1&ibox;
+ for (j1=j-1;j1<=j+1;j1++) {
+ element=box[i2][j1&ibox];
+ while (element != -1) {
+ if ((element < (act-lwindow)) || (element > (act+lwindow))) {
+ dx=sqr(series[act]-series[element]);
+ if (dx <= eps2) {
+ for (k=1;k<maxdim;k++) {
+ k1=k*delay;
+ dx += sqr(series[act+k1]-series[element+k1]);
+ if (dx <= eps2) {
+ k1=k-1;
+ lfound[k1][found[k1]]=element;
+ found[k1]++;
+ }
+ else
+ break;
+ }
+ }
+ }
+ element=liste[element];
+ }
+ }
+ }
+}
+
+void iterate_points(long act)
+{
+ double **lfactor;
+ double *dx;
+ unsigned int i,j,l,l1;
+ long k,element,**lcount;
+
+ check_alloc(lfactor=(double**)malloc(sizeof(double*)*(maxdim-1)));
+ check_alloc(lcount=(long**)malloc(sizeof(long*)*(maxdim-1)));
+ for (i=0;i<maxdim-1;i++) {
+ check_alloc(lfactor[i]=(double*)malloc(sizeof(double)*(maxiter+1)));
+ check_alloc(lcount[i]=(long*)malloc(sizeof(long)*(maxiter+1)));
+ }
+ check_alloc(dx=(double*)malloc(sizeof(double)*(maxiter+1)));
+
+ for (i=0;i<=maxiter;i++)
+ for (j=0;j<maxdim-1;j++) {
+ lfactor[j][i]=0.0;
+ lcount[j][i]=0;
+ }
+
+ for (j=mindim-2;j<maxdim-1;j++) {
+ for (k=0;k<found[j];k++) {
+ element=lfound[j][k];
+ for (i=0;i<=maxiter;i++)
+ dx[i]=sqr(series[act+i]-series[element+i]);
+ for (l=1;l<j+2;l++) {
+ l1=l*delay;
+ for (i=0;i<=maxiter;i++)
+ dx[i] += sqr(series[act+i+l1]-series[element+l1+i]);
+ }
+ for (i=0;i<=maxiter;i++)
+ if (dx[i] > 0.0){
+ lcount[j][i]++;
+ lfactor[j][i] += dx[i];
+ }
+ }
+ }
+ for (i=mindim-2;i<maxdim-1;i++)
+ for (j=0;j<=maxiter;j++)
+ if (lcount[i][j]) {
+ count[i][j]++;
+ lyap[i][j] += log(lfactor[i][j]/lcount[i][j])/2.0;
+ }
+
+ for (i=0;i<maxdim-1;i++){
+ free(lfactor[i]);
+ free(lcount[i]);
+ }
+ free(lcount);
+ free(lfactor);
+ free(dx);
+}
+
+int main(int argc,char **argv)
+{
+ char stdi=0;
+ double eps_fak;
+ double epsilon;
+ unsigned int i,j,l;
+ FILE *fout;
+
+ if (scan_help(argc,argv))
+ show_options(argv[0]);
+
+ scan_options(argc,argv);
+#ifndef OMIT_WHAT_I_DO
+ if (verbosity&VER_INPUT)
+ what_i_do(argv[0],WID_STR);
+#endif
+
+ infile=search_datafile(argc,argv,&column,verbosity);
+ if (infile == NULL)
+ stdi=1;
+
+ if (outfile == NULL) {
+ if (!stdi) {
+ check_alloc(outfile=(char*)calloc(strlen(infile)+6,1));
+ sprintf(outfile,"%s.lyap",infile);
+ }
+ else {
+ check_alloc(outfile=(char*)calloc(11,1));
+ sprintf(outfile,"stdin.lyap");
+ }
+ }
+ test_outfile(outfile);
+
+ series=get_series(infile,&length,exclude,column,verbosity);
+ rescale_data(series,length,&min,&max);
+
+ if (eps0set)
+ epsmin /= max;
+ if (eps1set)
+ epsmax /= max;
+
+ if (epsmin >= epsmax) {
+ epsmax=epsmin;
+ epscount=1;
+ }
+
+ if (reference > (length-maxiter-(maxdim-1)*delay))
+ reference=length-maxiter-(maxdim-1)*delay;
+ if ((maxiter+(maxdim-1)*delay) >= length) {
+ fprintf(stderr,"Too few points to handle these parameters!\n");
+ exit(LYAP_K__MAXITER_TOO_LARGE);
+ }
+
+ if (maxdim < 2)
+ maxdim=2;
+ if (mindim < 2)
+ mindim=2;
+ if (mindim > maxdim)
+ maxdim=mindim;
+
+ check_alloc(liste=(long*)malloc(sizeof(long)*(length)));
+ check_alloc(found=(long*)malloc(sizeof(long)*(maxdim-1)));
+ check_alloc(lfound=(long**)malloc(sizeof(long*)*(maxdim-1)));
+ for (i=0;i<maxdim-1;i++)
+ check_alloc(lfound[i]=(long*)malloc(sizeof(long)*(length)));
+ check_alloc(count=(long**)malloc(sizeof(long*)*(maxdim-1)));
+ for (i=0;i<maxdim-1;i++)
+ check_alloc(count[i]=(long*)malloc(sizeof(long)*(maxiter+1)));
+ check_alloc(lyap=(double**)malloc(sizeof(double*)*(maxdim-1)));
+ for (i=0;i<maxdim-1;i++)
+ check_alloc(lyap[i]=(double*)malloc(sizeof(double)*(maxiter+1)));
+
+ if (epscount == 1)
+ eps_fak=1.0;
+ else
+ eps_fak=pow(epsmax/epsmin,1.0/(double)(epscount-1));
+
+ fout=fopen(outfile,"w");
+ if (verbosity&VER_INPUT)
+ fprintf(stderr,"Opened %s for writing\n",outfile);
+ for (l=0;l<epscount;l++) {
+ epsilon=epsmin*pow(eps_fak,(double)l);
+ for (i=0;i<maxdim-1;i++)
+ for (j=0;j<=maxiter;j++) {
+ count[i][j]=0;
+ lyap[i][j]=0.0;
+ }
+ put_in_boxes(epsilon);
+ for (i=0;i<reference;i++) {
+ lfind_neighbors(i,epsilon);
+ iterate_points(i);
+ }
+ if (verbosity&VER_USR1)
+ fprintf(stderr,"epsilon= %e\n",epsilon*max);
+ for (i=mindim-2;i<maxdim-1;i++) {
+ fprintf(fout,"#epsilon= %e dim= %d\n",epsilon*max,i+2);
+ for (j=0;j<=maxiter;j++)
+ if (count[i][j])
+ fprintf(fout,"%d %e %ld\n",j,lyap[i][j]/count[i][j],count[i][j]);
+ fprintf(fout,"\n");
+ }
+ fflush(fout);
+ }
+ fclose(fout);
+ return 0;
+}
git://source.jalview.org / jabaws.git / blobdiff