+++ /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: Nov 30, 2000 */
-/*Changes:
- 12/11/05: Going multivariate
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include <limits.h>
-#include "routines/tsa.h"
-
-#define WID_STR "Performs simple noise reduction."
-
-#define BOX (unsigned int)512
-
-unsigned long length=ULONG_MAX,exclude=0;
-unsigned int comp=1,embed=5,delay=1,iterations=1,alldim;
-unsigned int verbosity=0x3;
-char *column=NULL;
-double eps=1.0e-3,epsvar;
-
-char *outfile=NULL,epsset=0,stdo=1,epsvarset=0;
-char *infile=NULL;
-double **series,**corr,*interval,*min,*hcor;
-long **box,*list,**nf;
-unsigned int **indexes;
-char dimset=0;
-
-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 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 no. of comp.,embedding dim. [default: %u,%u]\n",
- comp,embed);
- fprintf(stderr,"\t-d delay [default: 1]\n");
- fprintf(stderr,"\t-i iterations [default: 1]\n");
- fprintf(stderr,"\t-r neighborhoud size [default: (interval of data)/1000]\n");
- fprintf(stderr,"\t-v neighborhoud size (in units of the std. dev. of the "
- "data \n\t\t(overwrites -r) [default: not set]\n");
- fprintf(stderr,"\t-o output file name [Default: 'datafile'.laz.n,"
- "\n\t\twhere n is the number of the last iteration,"
- "\n\t\twithout -o the last iteration is written to 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='+ write output of all iterations to files'\n\t\t"
- "4='+ write the number of neighbors found for each point\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','c')) != NULL) {
- column=out;
- dimset=1;
- }
- if ((out=check_option(in,n,'m','2')) != NULL)
- sscanf(out,"%u,%u",&comp,&embed);
- if ((out=check_option(in,n,'d','u')) != NULL)
- sscanf(out,"%u",&delay);
- if ((out=check_option(in,n,'i','u')) != NULL)
- sscanf(out,"%u",&iterations);
- if ((out=check_option(in,n,'r','f')) != NULL) {
- epsset=1;
- sscanf(out,"%lf",&eps);
- }
- if ((out=check_option(in,n,'V','u')) != NULL)
- sscanf(out,"%u",&verbosity);
- if ((out=check_option(in,n,'v','f')) != NULL) {
- epsvarset=1;
- sscanf(out,"%lf",&epsvar);
- }
- if ((out=check_option(in,n,'o','o')) != NULL) {
- stdo=0;
- if (strlen(out) > 0)
- outfile=out;
- }
-}
-
-unsigned int correct(unsigned long n)
-{
- int i,i1,i2,j,j1,k;
- int ibox=BOX-1;
- unsigned int hdel,hcomp;
- double epsinv,dx;
- long element,nfound=0;
-
- epsinv=1./eps;
-
- for (i=0;i<alldim;i++)
- hcor[i]=0.0;
-
- i=(int)(series[0][n]*epsinv)&ibox;
- j=(int)(series[comp-1][n-(embed-1)*delay]*epsinv)&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) {
- for (k=0;k<alldim;k++) {
- hcomp=indexes[0][k];
- hdel=indexes[1][k];
- dx=fabs(series[hcomp][n-hdel]-series[hcomp][element-hdel]);
- if (dx > eps)
- break;
- }
- if (k == alldim) {
- nfound++;
- for (k=0;k<alldim;k++) {
- hcomp=indexes[0][k];
- hdel=indexes[1][k];
- hcor[k] += series[hcomp][element-hdel];
- }
- }
- element=list[element];
- }
- }
- }
- for (k=0;k<alldim;k++) {
- hcomp=indexes[0][k];
- hdel=indexes[1][k];
- corr[hcomp][n-hdel] += hcor[k]/nfound;
- nf[hcomp][n-hdel]++;
- }
-
- return nfound;
-}
-
-int main(int argc,char **argv)
-{
- char *ofname;
- char stdi=0;
- int iter;
- unsigned int *nmf;
- unsigned long n,i;
- double dav,dvar,maxinterval,maxdvar;
- FILE *file=NULL;
-
- 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));
- check_alloc(ofname=(char*)calloc(strlen(infile)+9,(size_t)1));
- sprintf(outfile,"%s.laz",infile);
- }
- else {
- check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
- check_alloc(ofname=(char*)calloc((size_t)14,(size_t)1));
- sprintf(outfile,"stdin.laz");
- }
- }
- else
- check_alloc(ofname=(char*)calloc(strlen(outfile)+10,(size_t)1));
-
-
- 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);
-
- check_alloc(interval=(double*)malloc(sizeof(double)*comp));
- check_alloc(min=(double*)malloc(sizeof(double)*comp));
-
- maxinterval=maxdvar=0.0;
- for (i=0;i<comp;i++) {
- rescale_data(series[i],length,&min[i],&interval[i]);
- if (interval[i] > maxinterval) maxinterval=interval[i];
- variance(series[i],length,&dav,&dvar);
- if (dvar > maxdvar) maxdvar=dvar;
- }
- alldim=comp*embed;
-
- check_alloc(nmf=(unsigned int*)malloc(sizeof(int)*length));
- check_alloc(list=(long*)malloc(sizeof(long)*length));
- check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
- for (n=0;n<BOX;n++)
- check_alloc(box[n]=(long*)malloc(sizeof(long)*BOX));
-
- check_alloc(nf=(long**)malloc(sizeof(long*)*comp));
- check_alloc(corr=(double**)malloc(sizeof(double*)*comp));
- for (i=0;i<comp;i++) {
- check_alloc(nf[i]=(long*)malloc(sizeof(long)*length));
- check_alloc(corr[i]=(double*)malloc(sizeof(double)*length));
- }
-
- indexes=make_multi_index(comp,embed,delay);
-
- if (epsset)
- eps/=maxinterval;
- else
- eps=1.0/1000.;
-
- if (epsvarset)
- eps=epsvar*maxdvar;
-
- for (iter=1;iter<=iterations;iter++) {
- make_multi_box2(series,box,list,length,BOX,comp,embed,delay,eps);
- for (n=0;n<length;n++) {
- for (i=0;i<comp;i++) {
- corr[i][n]=0.0;
- nf[i][n]=0;
- }
- nmf[n]=1;
- }
-
- check_alloc(hcor=(double*)malloc(sizeof(double)*alldim));
- for (n=(embed-1)*delay;n<length;n++)
- nmf[n]=correct(n);
- free(hcor);
-
- for (n=0;n<length;n++)
- for (i=0;i<comp;i++)
- if (nf[i][n])
- series[i][n]=corr[i][n]/nf[i][n];
-
- if ((verbosity&VER_USR1) && (iter < iterations)) {
- sprintf(ofname,"%s.%d",outfile,iter);
- test_outfile(ofname);
- file=fopen(ofname,"w");
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Opened %s for writing\n",ofname);
- if (stdo && (iter == iterations)) {
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Writing to stdout\n");
- }
- for (n=0;n<length;n++) {
- if (stdo && (iter == iterations)) {
- if (verbosity&VER_USR2) {
- for (i=0;i<comp;i++)
- fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
- fprintf(stdout,"%u\n",nmf[n]);
- }
- else {
- fprintf(stdout,"%e",series[0][n]*interval[0]+min[0]);
- for (i=1;i<comp;i++)
- fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
- fprintf(stdout,"\n");
- }
- }
- if (verbosity&VER_USR2) {
- for (i=0;i<comp;i++)
- fprintf(file,"%e ",series[i][n]*interval[i]+min[i]);
- fprintf(file,"%u\n",nmf[n]);
- }
- else {
- fprintf(file,"%e",series[0][n]*interval[0]+min[0]);
- for (i=1;i<comp;i++)
- fprintf(file," %e",series[i][n]*interval[i]+min[i]);
- fprintf(file,"\n");
- }
- }
- fclose(file);
- }
- if (iter == iterations) {
- if (!stdo || (verbosity&VER_USR1)) {
- sprintf(ofname,"%s.%d",outfile,iter);
- test_outfile(ofname);
- file=fopen(ofname,"w");
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Opened %s for writing\n",ofname);
- if (stdo && (iter == iterations)) {
- if (verbosity&VER_INPUT)
- fprintf(stderr,"Writing to stdout\n");
- }
- }
- for (n=0;n<length;n++) {
- if (stdo) {
- if (verbosity&VER_USR2) {
- for (i=0;i<comp;i++)
- fprintf(stdout,"%e ",series[i][n]*interval[i]+min[i]);
- fprintf(stdout,"%u\n",nmf[n]);
- }
- else {
- fprintf(stdout,"%e",series[0][n]*interval[0]+min[0]);
- for (i=1;i<comp;i++)
- fprintf(stdout," %e",series[i][n]*interval[i]+min[i]);
- fprintf(stdout,"\n");
- }
- }
- if (!stdo || (verbosity&VER_USR1)) {
- if (verbosity&VER_USR2) {
- for (i=0;i<comp;i++)
- fprintf(file,"%e ",series[i][n]*interval[i]+min[i]);
- fprintf(file,"%u\n",nmf[n]);
- }
- else {
- fprintf(file,"%e",series[0][n]*interval[0]+min[0]);
- for (i=1;i<comp;i++)
- fprintf(file," %e",series[i][n]*interval[i]+min[i]);
- fprintf(file,"\n");
- }
- }
- }
- if (!stdo || (verbosity&VER_USR1))
- fclose(file);
- }
- }
-
- /*cleaning up */
- for (i=0;i<comp;i++) {
- free(series[i]);
- free(nf[i]);
- free(corr[i]);
- }
- free(series);
- free(nf);
- free(corr);
-
- for (i=0;i<2;i++)
- free(indexes[i]);
- free(indexes);
-
- free(list);
- free(nmf);
- free(interval);
- free(min);
-
- for (i=0;i<BOX;i++)
- free(box[i]);
- free(box);
-
- if (outfile != NULL)
- free(outfile);
- if (ofname != NULL)
- free(ofname);
- if (infile != NULL)
- free(infile);
- if (column != NULL)
- free(column);
- /* end cleaning up */
-
- return 0;
-}
git://source.jalview.org / jabaws.git / blobdiff