+++ /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: Dec 10, 2005 */
-/*Changes:
- 12/10/05: It's multivariate now
- 12/16/05: Scaled <eps> and sigma(eps)
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <limits.h>
-#include <math.h>
-#include "routines/tsa.h"
-
-#define WID_STR "Determines the fraction of false nearest neighbors."
-
-char *outfile=NULL;
-char *infile=NULL;
-char stdo=1,dimset=0;
-char *column=NULL;
-unsigned long length=ULONG_MAX,exclude=0,theiler=0;
-unsigned int delay=1,maxdim=5,minemb=1;
-unsigned int comp=1,maxemb=5;
-unsigned int verbosity=0xff;
-double rt=2.0;
-double eps0=1.0e-5;
-double **series;
-double aveps,vareps;
-double varianz;
-
-#define BOX 1024
-int ibox=BOX-1;
-long **box,*list;
-unsigned int *vcomp,*vemb;
-unsigned long toolarge;
-
-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 ignore [default: 0]\n");
- fprintf(stderr,"\t-c columns to read [default: 1]\n");
- fprintf(stderr,"\t-m min. test embedding dimension [default: %u]\n",minemb);
- fprintf(stderr,"\t-M # of components,max. emb. dim. [default: %u,%u]\n",
- comp,maxemb);
- fprintf(stderr,"\t-d delay [default: 1]\n");
- fprintf(stderr,"\t-f escape factor [default: %.2lf]\n",rt);
- fprintf(stderr,"\t-t theiler window [default: 0]\n");
- fprintf(stderr,"\t-o output file [default: 'datafile'.fnn; without -o"
- " stdout]\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='+ information about the current state\n");
- fprintf(stderr,"\t-h show these options\n");
- exit(0);
-}
-
-void scan_options(int n,char **in)
-{
- char *out;
-
- if ((out=check_option(in,n,'l','u')) != NULL)
- sscanf(out,"%lu",&length);
- if ((out=check_option(in,n,'x','u')) != NULL)
- sscanf(out,"%lu",&exclude);
- if ((out=check_option(in,n,'c','s')) != NULL)
- column=out;
- if ((out=check_option(in,n,'m','u')) != NULL)
- sscanf(out,"%u",&minemb);
- if ((out=check_option(in,n,'M','2')) != NULL) {
- sscanf(out,"%u,%u",&comp,&maxemb);
- maxdim=comp*(maxemb+1);
- dimset=1;
- }
- if ((out=check_option(in,n,'d','u')) != NULL)
- sscanf(out,"%u",&delay);
- if ((out=check_option(in,n,'f','f')) != NULL)
- sscanf(out,"%lf",&rt);
- if ((out=check_option(in,n,'t','u')) != NULL)
- sscanf(out,"%lu",&theiler);
- if ((out=check_option(in,n,'V','u')) != NULL)
- sscanf(out,"%u",&verbosity);
- if ((out=check_option(in,n,'o','o')) != NULL) {
- stdo=0;
- if (strlen(out) > 0)
- outfile=out;
- }
-}
-
-void mmb(unsigned int hdim,unsigned int hemb,double eps)
-{
- unsigned long i;
- long x,y;
-
- for (x=0;x<BOX;x++)
- for (y=0;y<BOX;y++)
- box[x][y] = -1;
-
- for (i=0;i<length-(maxemb+1)*delay;i++) {
- x=(long)(series[0][i]/eps)&ibox;
- y=(long)(series[hdim][i+hemb]/eps)&ibox;
- list[i]=box[x][y];
- box[x][y]=i;
- }
-}
-
-char find_nearest(long n,unsigned int dim,double eps)
-{
- long x,y,x1,x2,y1,i,i1,ic,ie;
- long element,which= -1;
- double dx,maxdx,mindx=1.1,hfactor,factor;
-
- ic=vcomp[dim];
- ie=vemb[dim];
- x=(long)(series[0][n]/eps)&ibox;
- y=(long)(series[ic][n+ie]/eps)&ibox;
-
- for (x1=x-1;x1<=x+1;x1++) {
- x2=x1&ibox;
- for (y1=y-1;y1<=y+1;y1++) {
- element=box[x2][y1&ibox];
- while (element != -1) {
- if (labs(element-n) > theiler) {
- maxdx=fabs(series[0][n]-series[0][element]);
- for (i=1;i<=dim;i++) {
- ic=vcomp[i];
- i1=vemb[i];
- dx=fabs(series[ic][n+i1]-series[ic][element+i1]);
- if (dx > maxdx)
- maxdx=dx;
- }
- if ((maxdx < mindx) && (maxdx > 0.0)) {
- which=element;
- mindx=maxdx;
- }
- }
- element=list[element];
- }
- }
- }
-
- if ((which != -1) && (mindx <= eps) && (mindx <= varianz/rt)) {
- aveps += mindx;
- vareps += mindx*mindx;
- factor=0.0;
- for (i=1;i<=comp;i++) {
- ic=vcomp[dim+i];
- ie=vemb[dim+i];
- hfactor=fabs(series[ic][n+ie]-series[ic][which+ie])/mindx;
- if (hfactor > factor)
- factor=hfactor;
- }
- if (factor > rt)
- toolarge++;
- return 1;
- }
- return 0;
-}
-
-int main(int argc,char **argv)
-{
- char stdi=0;
- FILE *file=NULL;
- double min,inter=0.0,ind_inter,epsilon,av,ind_var;
- char *nearest,alldone;
- long i;
- unsigned int dim,emb;
- unsigned long donesofar;
-
- 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,NULL,verbosity);
- if (infile == NULL)
- stdi=1;
-
- if (outfile == NULL) {
- if (!stdi) {
- check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
- strcpy(outfile,infile);
- strcat(outfile,".fnn");
- }
- else {
- check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
- strcpy(outfile,"stdin.fnn");
- }
- }
- if (!stdo)
- test_outfile(outfile);
-
- if (column == NULL)
- series=(double**)get_multi_series(infile,&length,exclude,&comp,"",dimset,
- verbosity);
- else
- series=(double**)get_multi_series(infile,&length,exclude,&comp,column,
- dimset,verbosity);
-
- for (i=0;i<comp;i++) {
- rescale_data(series[i],length,&min,&ind_inter);
- variance(series[i],length,&av,&ind_var);
- if (i == 0) {
- varianz=ind_var;
- inter=ind_inter;
- }
- else {
- varianz=(varianz>ind_var)?ind_var:varianz;
- inter=(inter<ind_inter)?ind_inter:inter;
- }
- }
-
- check_alloc(list=(long*)malloc(sizeof(long)*length));
- check_alloc(nearest=(char*)malloc(length));
- check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
- for (i=0;i<BOX;i++)
- check_alloc(box[i]=(long*)malloc(sizeof(long)*BOX));
-
- if (!stdo) {
- file=fopen(outfile,"w");
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Opened %s for writing\n",outfile);
- }
- else {
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Writing to stdout\n");
- }
- check_alloc(vcomp=(unsigned int*)malloc(sizeof(int)*(maxdim)));
- check_alloc(vemb=(unsigned int*)malloc(sizeof(int)*(maxdim)));
- for (i=0;i<maxdim;i++) {
- if (comp == 1) {
- vcomp[i]=0;
- vemb[i]=i;
- }
- else {
- vcomp[i]=i%comp;
- vemb[i]=(i/comp)*delay;
- }
- }
- for (emb=minemb;emb<=maxemb;emb++) {
- dim=emb*comp-1;
- epsilon=eps0;
- toolarge=0;
- alldone=0;
- donesofar=0;
- aveps=0.0;
- vareps=0.0;
- for (i=0;i<length;i++)
- nearest[i]=0;
- if (verbosity&VER_USR1)
- fprintf(stderr,"Start for dimension=%u\n",dim+1);
- while (!alldone && (epsilon < 2.*varianz/rt)) {
- alldone=1;
- mmb(vcomp[dim],vemb[dim],epsilon);
- for (i=0;i<length-maxemb*delay;i++)
- if (!nearest[i]) {
- nearest[i]=find_nearest(i,dim,epsilon);
- alldone &= nearest[i];
- donesofar += (unsigned long)nearest[i];
- }
- if (verbosity&VER_USR1)
- fprintf(stderr,"Found %lu up to epsilon=%e\n",donesofar,epsilon*inter);
- epsilon*=sqrt(2.0);
- if (!donesofar)
- eps0=epsilon;
- }
- if (donesofar == 0) {
- fprintf(stderr,"Not enough points found!\n");
- exit(FALSE_NEAREST_NOT_ENOUGH_POINTS);
- }
- aveps *= (1./(double)donesofar);
- vareps *= (1./(double)donesofar);
- if (stdo) {
- fprintf(stdout,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
- aveps*inter,sqrt(vareps)*inter);
- fflush(stdout);
- }
- else {
- fprintf(file,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
- aveps*inter,sqrt(vareps)*inter);
- fflush(file);
- }
- }
- if (!stdo)
- fclose(file);
-
- if (infile != NULL)
- free(infile);
- if (outfile != NULL)
- free(outfile);
- free(series);
- free(list);
- free(nearest);
- for (i=0;i<BOX;i++)
- free(box[i]);
- free(box);
-
- return 0;
-}
git://source.jalview.org / jabaws.git / blobdiff