/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 

Generated on Mon Jan 29 17:09:12 2007 for Diffit by  doxygen 1.4.6