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
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
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
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
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
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
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
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
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
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