00001 #include "mex.h"
00002 #include <iostream>
00003 #include <stdlib.h>
00004 #include <fstream>
00005 #include <math.h>
00006 #include <string>
00007 #include <getopt.h>
00008 #include "nr.h"
00009
00010 using namespace std;
00011
00012 #define TRUE 1
00013 #define FALSE 0
00014
00015
00016 typedef struct
00017 {
00018 double *t;
00019 double *y;
00020 double *sig;
00021 double maxLambda;
00022 long nKnots;
00023 int gcv;
00024 } glob;
00025
00026
00027
00028 #define DATA prhs[0]
00029 #define ALPHA prhs[1]
00030 #define MAXA prhs[2]
00031 #define NOGCV prhs[3]
00032
00033
00034
00035 #define SPLINE plhs[0]
00036 #define ALPHA_OUT plhs[1]
00037
00038 void mexFunction( int nlhs,
00039 mxArray *plhs[],
00040 int nrhs,
00041 const mxArray *prhs[]
00042 )
00043 {
00044
00045 glob globs;
00046 long i,j,n,col;
00047 double *g,*gam,alpha;
00048
00049 double *data_,*alpha_,*maxa_,*nogcv_;
00050
00051 double *spline_,*alpha_out_;
00052
00053
00054 data_=mxGetPr(DATA);
00055 alpha_=mxGetPr(ALPHA);
00056 maxa_=mxGetPr(MAXA);
00057 nogcv_=mxGetPr(NOGCV);
00058
00059 globs.nKnots=mxGetM(DATA);
00060 col=mxGetN(DATA);
00061
00062 n=globs.nKnots;
00063 g=dvector(1,n);
00064 gam=dvector(1,n);
00065
00066 globs.t=dvector(1,n);
00067 globs.y=dvector(1,n);
00068 globs.sig=dvector(1,n);
00069
00070 if(col==2)
00071 {
00072 for(i=1;i<=n;i++)
00073 {
00074 globs.t[i]=data_[i-1];
00075 globs.y[i]=data_[n+i-1];
00076 globs.sig[i]=1.;
00077 }
00078 }
00079 else if(col==3)
00080 {
00081 for(i=1;i<=n;i++)
00082 {
00083 globs.t[i]=data_[i-1];
00084 globs.y[i]=data_[n+i-1];
00085 globs.sig[i]=data_[2*n+i-1];
00086 }
00087 }
00088 else
00089 {
00090 mexWarnMsgTxt("Insufficient data structure. Exit.\n");
00091 return;
00092 }
00093
00094 alpha=fabs(alpha_[0]);
00095 globs.maxLambda=fabs(maxa_[0]);
00096 globs.gcv=(int)nogcv_[0];
00097
00098
00099 if(alpha > globs.maxLambda)
00100 alpha=globs.maxLambda;
00101
00102 if(globs.gcv==TRUE)
00103 splines_gcv(globs.t,globs.y,globs.sig,n,&alpha,g,gam,globs.maxLambda);
00104 else
00105 splines(globs.t,globs.y,globs.sig,n,alpha,g,gam);
00106
00107
00108
00109 SPLINE=mxCreateDoubleMatrix(n,3,mxREAL);
00110 spline_=mxGetPr(SPLINE);
00111 ALPHA_OUT=mxCreateDoubleMatrix(1,1,mxREAL);
00112 alpha_out_=mxGetPr(ALPHA_OUT);
00113 alpha_out_[0] = alpha;
00114
00115 spline_[0]=globs.t[1];
00116 spline_[n]=g[1];
00117 spline_[2*n]=0.;
00118 for(i=2;i<=n-1;i++)
00119 {
00120 spline_[i-1]=globs.t[i];
00121 spline_[n+i-1]=g[i];
00122 spline_[2*n+i-1]=gam[i-1];
00123 }
00124 spline_[n-1]=globs.t[n];
00125 spline_[2*n-1]=g[n];
00126 spline_[3*n-1]=0.;
00127
00128
00129 free_dvector(globs.y,1,n);
00130 free_dvector(globs.t,1,n);
00131 free_dvector(globs.sig,1,n);
00132 free_dvector(g,1,n);
00133 free_dvector(gam,1,n);
00134 }
00135