+++ /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 5, 2004 */
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <limits.h>
-#include <ctype.h>
-#include "routines/tsa.h"
-
-#define WID_STR "Tests for nonstationarity by means of the average\n\t\
-forecast error for a zeroth order fit"
-
-
-#ifndef _MATH_H
-#include <math.h>
-#endif
-
-/*number of boxes for the neighbor search algorithm*/
-#define NMAX 128
-
-unsigned int nmax=(NMAX-1);
-long **box,*list;
-unsigned long *found;
-double *series,*series1,*series2;
-double interval,min,epsilon;
-
-char epsset=0,causalset=0;
-char *infile=NULL;
-char *outfile=NULL,stdo=1,centerset=0;
-char *firstwindow,*secondwindow,**window;
-unsigned int COLUMN=1,pieces;
-unsigned int verbosity=0xff;
-int DIM=3,DELAY=1,MINN=30,STEP=1;
-int firstoffset= -1,secondoffset= -1;
-double EPS0=1.e-3,EPSF=1.2;
-unsigned long LENGTH=ULONG_MAX,exclude=0,center,causal;
-
-void show_options(char *progname)
-{
- what_i_do(progname,WID_STR);
- fprintf(stderr," Usage: %s -# [other 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 to use [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 embedding dimension [default: 3]\n");
- fprintf(stderr,"\t-d delay [default: 1]\n");
- fprintf(stderr,"\t-# # of pieces [no default]\n");
- fprintf(stderr,"\t-1 which pieces for the first window "
- "[default: 1-pieces]\n");
- fprintf(stderr,"\t-2 which pieces for the second window "
- "[default: 1-pieces]\n");
- fprintf(stderr,"\t-n # of reference points in the window [default: all]\n");
- fprintf(stderr,"\t-k minimal number of neighbors for the fit "
- "[default: 30]\n");
- fprintf(stderr,"\t-r neighborhoud size to start with "
- "[default: (data interval)/1000]\n");
- fprintf(stderr,"\t-f factor to increase size [default: 1.2]\n");
- fprintf(stderr,"\t-s steps to forecast [default: 1]\n");
- fprintf(stderr,"\t-C width of causality window [default: steps]\n");
- fprintf(stderr,"\t-o output file [default: 'datafile.nsz',"
- " without -o: stdout]\n");
- fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
- "0='only panic messages'\n\t\t"
- "1='+ input/output messages'\n");
- fprintf(stderr,"\t-h show these options\n");
- fprintf(stderr,"\n\t The -# option has to be set\n");
- exit(0);
-}
-
-void parse_minus(char *str,char *array,char *wopt)
-{
- int cm=0,i,strl,n1,n2;
-
- strl=strlen(str);
- for (i=0;i<strl;i++)
- if (str[i] == '-')
- cm++;
- if (cm > 1) {
- fprintf(stderr,"Invalid string for the %s option! "
- "Please consult the help-page\n",wopt);
- exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
- }
- if (cm == 0) {
- sscanf(str,"%d",&n1);
- n1--;
- if (n1 < 0) {
- fprintf(stderr,"Numbers in %s option must be larger than 0!\n",wopt);
- exit(NSTAT_Z__NOT_UNSIGNED_FOR_OPTION);
- }
- if (n1 >= pieces) {
- fprintf(stderr,"Numbers in %s option must be smaller than %u!\n",wopt,
- pieces);
- exit(NSTAT_Z__TOO_LARGE_FOR_OPTION);
- }
- array[n1]=1;
- }
- else {
- sscanf(str,"%d-%d",&n1,&n2);
- n1--;
- n2--;
- if ((n1 < 0) || (n2 < 0)) {
- fprintf(stderr,"Numbers in %s option must be larger than 0!\n",wopt);
- exit(NSTAT_Z__NOT_UNSIGNED_FOR_OPTION);
- }
- if ((n1 >= pieces) || (n2 >= pieces)) {
- fprintf(stderr,"Numbers in %s option must be smaller than %u!\n",wopt,
- pieces+1);
- exit(NSTAT_Z__TOO_LARGE_FOR_OPTION);
- }
- if (n2 < n1) {
- i=n1;
- n1=n2;
- n2=i;
- }
- for (i=n1;i<=n2;i++)
- array[i]=1;
- }
-}
-
-void parse_comma(char *str,char *array,char *wopt)
-{
- unsigned int strl,i,cp=1,which,iwhich;
- char **hstr;
-
- strl=strlen(str);
- for (i=0;i<strl;i++)
- if (str[i] == ',')
- cp++;
-
- if (cp == 1) {
- parse_minus(str,array,wopt);
- return ;
- }
-
- check_alloc(hstr=(char**)malloc(sizeof(char*)*cp));
- for (i=0;i<cp;i++)
- check_alloc(hstr[i]=(char*)calloc(strl,1));
-
- which=iwhich=0;
- for (i=0;i<strl;i++) {
- if (str[i] != ',')
- hstr[which][iwhich++]=str[i];
- else {
- which++;
- iwhich=0;
- }
- }
- for (i=0;i<cp;i++) {
- if (hstr[i][0] == '\0') {
- fprintf(stderr,"Invalid string for the %s option! "
- "Please consult the help-page\n",wopt);
- exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
- }
- if (!isdigit(hstr[i][strlen(hstr[i])-1])) {
- fprintf(stderr,"Invalid string for the %s option! "
- "Please consult the help-page\n",wopt);
- exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
- }
- parse_minus(hstr[i],array,wopt);
- }
- for (i=0;i<cp;i++)
- free(hstr[i]);
- free(hstr);
-}
-
-void parse_out(char *str,char *array,char *which)
-{
- unsigned int i;
- char test;
-
- for (i=0;i<pieces;i++)
- array[i]=0;
-
- for (i=0;i<strlen(str);i++) {
- test= (str[i] == '-') || (str[i] == ',') || isdigit(str[i]);
- if (!test) {
- fprintf(stderr,"Invalid string for the %s option! "
- "Please consult the help-page\n",which);
- exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
- }
- }
- if (!isdigit(str[strlen(str)-1])) {
- fprintf(stderr,"Invalid string for the %s option! "
- "Please consult the help-page\n",which);
- exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
- }
- parse_comma(str,array,which);
-}
-
-void parse_offset(char *str,int *iwhich,char *array,char *which)
-{
- int i,strl;
-
- if (str[0] != '+')
- return;
- strl=strlen(str);
- for (i=1;i<strl;i++)
- if (!isdigit(str[i])) {
- fprintf(stderr,"Invalid string for the %s option! "
- "Please consult the help-page\n",which);
- exit(NSTAT_Z__INVALID_STRING_FOR_OPTION);
- }
- sscanf(str,"+%d",iwhich);
- for (i=0;i<pieces;i++)
- array[i]=0;
-}
-
-void scan_options(int n,char **in)
-{
- unsigned int i;
- char *out,piecesset=0;
-
- 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','u')) != NULL)
- sscanf(out,"%u",&COLUMN);
- if ((out=check_option(in,n,'m','u')) != NULL)
- sscanf(out,"%u",&DIM);
- if ((out=check_option(in,n,'d','u')) != NULL)
- sscanf(out,"%u",&DELAY);
- if ((out=check_option(in,n,'V','u')) != NULL)
- sscanf(out,"%u",&verbosity);
- if ((out=check_option(in,n,'#','u')) != NULL) {
- sscanf(out,"%u",&pieces);
- if (pieces < 1)
- pieces=1;
- piecesset=1;
- check_alloc(firstwindow=(char*)malloc(pieces));
- check_alloc(secondwindow=(char*)malloc(pieces));
- for (i=0;i<pieces;i++)
- firstwindow[i]=secondwindow[i]=1;
- check_alloc(window=(char**)malloc(sizeof(char*)*pieces));
- for (i=0;i<pieces;i++)
- check_alloc(window[i]=(char*)malloc(pieces));
- }
- if (!piecesset) {
- fprintf(stderr,"\tThe -# option wasn't set. Please add it!\n");
- exit(NSTAT_Z__OPTION_NOT_SET);
- }
- if ((out=check_option(in,n,'1','s')) != NULL) {
- parse_offset(out,&firstoffset,firstwindow,"-1");
- if (firstoffset == -1)
- parse_out(out,firstwindow,"-1");
- }
- if ((out=check_option(in,n,'2','s')) != NULL) {
- parse_offset(out,&secondoffset,secondwindow,"-2");
- if (secondoffset == -1)
- parse_out(out,secondwindow,"-2");
- }
- if ((out=check_option(in,n,'n','u')) != NULL) {
- sscanf(out,"%lu",¢er);
- centerset=1;
- }
- if ((out=check_option(in,n,'k','u')) != NULL)
- sscanf(out,"%u",&MINN);
- if ((out=check_option(in,n,'r','f')) != NULL) {
- epsset=1;
- sscanf(out,"%lf",&EPS0);
- }
- if ((out=check_option(in,n,'f','f')) != NULL)
- sscanf(out,"%lf",&EPSF);
- if ((out=check_option(in,n,'s','u')) != NULL)
- sscanf(out,"%u",&STEP);
- if ((out=check_option(in,n,'C','u')) != NULL) {
- sscanf(out,"%lu",&causal);
- causalset=1;
- }
- if ((out=check_option(in,n,'o','o')) != NULL) {
- stdo=0;
- if (strlen(out) > 0)
- outfile=out;
- }
-}
-
-double make_fit(long act,unsigned long number)
-{
- double casted=0.0,*help;
- int i;
-
- help=series1+STEP;
- for (i=0;i<number;i++) {
- casted += help[found[i]];
- }
- casted /= number;
-
- return sqr(casted-series2[act+STEP]);
-}
-
-int main(int argc,char **argv)
-{
- char stdi=0;
- char alldone,*done,sdone;
- long i,first,second,pstart;
- unsigned long *hfound;
- unsigned long actfound;
- unsigned long clength;
- double *rms,av,error;
- FILE *file=NULL;
-
- if (scan_help(argc,argv))
- show_options(argv[0]);
-
- scan_options(argc,argv);
-
- if (!causalset)
- causal=STEP;
-
-#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)+5,(size_t)1));
- sprintf(outfile,"%s.nsz",infile);
- }
- else {
- check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
- sprintf(outfile,"stdin.nsz");
- }
- }
- if (!stdo)
- test_outfile(outfile);
-
- series=(double*)get_series(infile,&LENGTH,exclude,COLUMN,verbosity);
-
- rescale_data(series,LENGTH,&min,&interval);
-
- check_alloc(list=(long*)malloc(sizeof(long)*LENGTH));
- check_alloc(found=(unsigned long*)malloc(sizeof(long)*LENGTH));
- check_alloc(hfound=(unsigned long*)malloc(sizeof(long)*LENGTH));
- check_alloc(done=(char*)malloc(sizeof(char)*LENGTH));
- check_alloc(box=(long**)malloc(sizeof(long*)*NMAX));
-
- for (i=0;i<NMAX;i++)
- check_alloc(box[i]=(long*)malloc(sizeof(long)*NMAX));
-
- if (epsset)
- EPS0 /= interval;
-
- clength=(LENGTH-(DIM-1)*DELAY)/pieces;
- if ((clength-(DIM-1)*DELAY-STEP) < MINN) {
- fprintf(stderr,"You chose too many pieces and will never find enough"
- " neighbors!\n");
- exit(NSTAT_Z__TOO_MANY_PIECES);
- }
- check_alloc(rms=(double*)malloc(sizeof(double)*pieces));
- for (i=0;i<pieces;i++) {
- series1=series+i*clength;
- variance(series1,clength,&av,&rms[i]);
- }
-
- pstart=(DIM-1)*DELAY;
- if (!centerset)
- center=clength-STEP;
- else
- center=(center < (clength-STEP-pstart)) ? center : clength-STEP-pstart;
-
- if (stdo) {
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Writing to stdout\n");
- }
- else {
- file=fopen(outfile,"w");
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Opened %s for writing\n",outfile);
- }
- for (first=0;first<pieces;first++)
- for (second=0;second<pieces;second++)
- window[first][second]=firstwindow[first]&&secondwindow[second];
- if (firstoffset != -1) {
- for (second=0;second<pieces;second++)
- for (first=second-firstoffset;first<=second+firstoffset;first++)
- if ((first >= 0) && (first < pieces))
- window[first][second]=secondwindow[second];
- }
- if (secondoffset != -1) {
- for (first=0;first<pieces;first++)
- for (second=first-secondoffset;second<=first+secondoffset;second++)
- if ((second >= 0) && (second < pieces))
- window[first][second]=firstwindow[first];
- }
-
- free(firstwindow);
- free(secondwindow);
-
- for (first=0;first<pieces;first++) {
- sdone=0;
- for (second=0;second<pieces;second++) {
- if (window[first][second]) {
- sdone=1;
- series1=series+first*clength;
- series2=series+second*clength;
- for (i=0;i<LENGTH;i++)
- done[i]=0;
- alldone=0;
- epsilon=EPS0/EPSF;
- error=0.0;
- while (!alldone) {
- alldone=1;
- epsilon*=EPSF;
- make_box(series1,box,list,clength-STEP,NMAX,(unsigned int)DIM,
- (unsigned int)DELAY,epsilon);
- for (i=pstart;i<pstart+center;i++)
- if (!done[i]) {
- actfound=find_neighbors(series1,box,list,series2+i,clength,NMAX,
- (unsigned int)DIM,(unsigned int)DELAY,
- epsilon,hfound);
- actfound=exclude_interval(actfound,i-causal+1,
- i+causal+pstart-1,hfound,found);
- if (actfound >= MINN) {
- error += make_fit(i,actfound);
- done[i]=1;
- }
- alldone &= done[i];
- }
- }
- if (stdo)
- fprintf(stdout,"%ld %ld %e\n",first+1,second+1,
- sqrt(error/center)/rms[second]);
- else {
- fprintf(file,"%ld %ld %e\n",first+1,second+1,
- sqrt(error/center)/rms[second]);
- fflush(file);
- }
- }
- }
- if (sdone) {
- if (stdo)
- fprintf(stdout,"\n");
- else
- fprintf(file,"\n");
- }
- }
-
- if (!stdo)
- fclose(file);
-
- if (outfile != NULL)
- free(outfile);
- free(list);
- free(found);
- free(hfound);
- free(done);
- for (i=0;i<NMAX;i++)
- free(box[i]);
- free(box);
- for (i=0;i<pieces;i++)
- free(window[i]);
- free(window);
- free(rms);
- free(series);
-
- return 0;
-}
-
git://source.jalview.org / jabaws.git / blobdiff