- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# /home/peifer/diffit/createSpline.cc

Go to the documentation of this file.

    00001 #include <iostream>
    00002 #include <stdlib.h>
    00003 #include <fstream>
    00004 #include <math.h>
    00005 #include<string>
    00006 #include <getopt.h>
    00007 #include "nr.h"
    00008
    00009 using namespace std;
    00010
    00011 #define TRUE 1
    00012 #define FALSE 0
    00013
    00014 char usage[]="\n\t    --- createSpline : Creating spline data for diffit ---\n \
    00015  createSpline, v 1.0 25/Jul./2005  Univ. of Freiburg, Author: Martin Peifer\n\n \
    00016       SYNOPSIS \n \t createSpline  [Options]  <File>\n\n \
    00017       DESCRIPTION\n \
    00018       \t -h -? -help \t help\n \
    00019       \t -alpha <val> \t smoothing parameter\n \
    00020       \t -nogcv \t switches cross validation off\n \
    00021       \t -max  <val> \t upper bound for linesearch\n \
    00022       \t -o     <file>\t output file\n";
    00023
    00024 typedef struct
    00025 {
    00026   double *t;
    00027   double *y;
    00028   double *sig;
    00029   double maxLambda;
    00030   double alpha;
    00031   long nKnots;
    00032   int gcv;
    00033   char *outFile;
    00034   char *inFile;
    00035   int outset;
    00036 } glob;
    00037
    00038
    00039 void parseTheOpts(int argc, char *argv[],glob *globs)
    00040 {
    00041   int longindex,opt;
    00042   long i,j,columns=0,colvor=0;
    00043   double dum;
    00044   char cdum;
    00045
    00046   ifstream in,inp;
    00047
    00048   //option definition
    00049   static struct option longopts[]={
    00050     {"help", 0, 0,  'h'},
    00051     {"help", 0, 0,  '?'},
    00052     {"alpha"  , 1, 0,   1 },
    00053     {"nogcv"  , 0, 0,   2 },
    00054     {"max"   , 1, 0,   3 },
    00055     {"o"      , 1, 0,   4 },
    00056     {0, 0, 0, 0}
    00057   };
    00058
    00059   //initialise struct glob
    00060   globs->alpha=1;
    00061   globs->gcv=TRUE;
    00062   globs->maxLambda=1e10;
    00063   globs->nKnots=0;
    00064   globs->outset=FALSE;
    00065
    00066   while ((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
    00067     {
    00068       switch(opt)
    00069         {
    00070         case 'h':
    00071           cout << usage << endl;
    00072           exit(1);
    00073           break;
    00074         case '?':
    00075           cerr << usage << endl;
    00076           exit(1);
    00077           break;
    00078         }
    00079     }
    00080   //getFileName
    00081   if((argc-optind)!=1)
    00082     {
    00083       cerr << "No/Too many datafile(s) specified.  \n\n";
    00084       exit(1);
    00085     }
    00086   else
    00087     {
    00088       globs->inFile = new char[strlen(argv[optind])+1];
    00089       strcpy(globs->inFile,argv[optind]);
    00090     }
    00091
    00092 #ifdef DEBUG
    00093   cout << "Debug ->  Input filename : " << globs->inFile << endl;
    00094 #endif
    00095   optind=0;
    00096
    00097   //determine data structure and read data
    00098
    00099   in.open(globs->inFile);
    00100   if(in==NULL)
    00101     {
    00102       cerr << "Input file does not exist.\n";
    00103       exit(1);
    00104     }
    00105
    00106   i=0;
    00107   while(!in.eof())
    00108     {
    00109       cdum=in.peek();
    00110       while((cdum == ' ' || cdum == '\t' || cdum == '\n'))
    00111         {
    00112           if(cdum=='\n')
    00113             {
    00114               if(i==0)
    00115                 colvor=columns;
    00116               i=1;
    00117               if((columns != 2 && columns != 3)|| colvor != columns)
    00118                 {
    00119                   cerr << "Insufficient data format in line " <<  globs->nKnots+1 << ".\n";
    00120                   exit(1);
    00121                 }
    00122               globs->nKnots++;
    00123               colvor=columns;
    00124               columns=0;
    00125             }
    00126           cdum=in.get();
    00127           cdum=in.peek();
    00128         }
    00129       columns++;
    00130       in >> dum;
    00131     }
    00132   columns=colvor;
    00133 #ifdef DEBUG
    00134   cout << "Debug -> Columns & Lines : "<< columns << "  " << globs->nKnots << endl;
    00135 #endif
    00136   // allocate memory
    00137   in.close();
    00138
    00139   globs->t=dvector(1,globs->nKnots);
    00140   globs->y=dvector(1,globs->nKnots);
    00141   globs->sig=dvector(1,globs->nKnots);
    00142
    00143   //finally, read the data
    00144
    00145   inp.open(globs->inFile);
    00146   for(i=1;i<=globs->nKnots;i++)
    00147     {
    00148       if(columns==2)
    00149         {
    00150           inp >> globs->t[i] >> globs->y[i];
    00151           globs->sig[i]=1.;
    00152         }
    00153       else
    00154         inp>> globs->t[i] >> globs->y[i] >> globs->sig[i];
    00155 #ifdef DEBUG
    00156       cout << globs->t[i] << "\t" <<  globs->y[i] << "\t" <<  globs->sig[i] << endl;
    00157 #endif
    00158     }
    00159   inp.close();
    00160
    00161   //Parse the remaining parameters
    00162
    00163   while ((opt=getopt_long_only(argc,argv,"h?",longopts,&longindex)) != -1)
    00164     {
    00165       switch(opt)
    00166         {
    00167         case 'h':
    00168           cout << usage << endl;
    00169           exit(1);
    00170           break;
    00171         case '?':
    00172           cerr << usage << endl;
    00173           exit(1);
    00174           break;
    00175         case 1:
    00176           globs->alpha=fabs(atof(optarg));
    00177           break;
    00178         case 2:
    00179           globs->gcv=FALSE;
    00180           break;
    00181         case 3:
    00182           globs->maxLambda=fabs(atof(optarg));
    00183           break;
    00184         case 4:
    00185           globs->outFile=new char[strlen(optarg)+1];
    00186           strcpy(globs->outFile,optarg);
    00187           globs->outset=TRUE;
    00188           break;
    00189         default:
    00190           cerr << endl;
    00191           cerr << "Parsing parameters produced errors.\n\n";
    00192           exit(1);
    00193         }
    00194     }
    00195
    00196   if(globs->outset==FALSE)
    00197     globs->outFile="spline.dat";
    00198 }
    00199
    00200
    00201
    00202 main(int argc,char *argv[])
    00203 {
    00204
    00205   glob globs;
    00206   long i,n;
    00207   double *g,*gam,alpha;
    00208
    00209   ofstream out;
    00210
    00211   parseTheOpts(argc,argv,&globs);
    00212
    00213   n=globs.nKnots;
    00214   g=dvector(1,n);
    00215   gam=dvector(1,n);
    00216   out.open(globs.outFile);
    00217   alpha=globs.alpha;
    00218
    00219   if(alpha > globs.maxLambda)
    00220     alpha=globs.maxLambda;
    00221
    00222   if(globs.gcv==TRUE)
    00223     splines_gcv(globs.t,globs.y,globs.sig,n,&alpha,g,gam,globs.maxLambda);
    00224   else
    00225     splines(globs.t,globs.y,globs.sig,n,alpha,g,gam);
    00226
    00227   //writing output
    00228
    00229   out << n <<endl;
    00230   out << globs.t[1] << "\t" << g[1] << "\t" << "0" << endl;
    00231   for(i=2;i<=n-1;i++)
    00232     out << globs.t[i] << "\t" << g[i] << "\t" << gam[i-1] << endl;
    00233   out << globs.t[n] << "\t" << g[n] << "\t" << "0" << endl;
    00234
    00235   out.close();
    00236   //free memory
    00237   delete globs.inFile;
    00238   if(globs.outset==TRUE)
    00239     delete globs.outFile;
    00240   free_dvector(globs.y,1,n);
    00241   free_dvector(globs.t,1,n);
    00242   free_dvector(globs.sig,1,n);
    00243   free_dvector(g,1,n);
    00244   free_dvector(gam,1,n);
    00245 }
    00246