/home/peifer/diffit/createSpline_mex.cc

Go to the documentation of this file.
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 // Input
00028 #define DATA   prhs[0]
00029 #define ALPHA  prhs[1]
00030 #define MAXA   prhs[2]
00031 #define NOGCV  prhs[3]
00032 
00033 // Output
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   //mex input
00049   double *data_,*alpha_,*maxa_,*nogcv_;
00050   //mex output
00051   double *spline_,*alpha_out_;
00052    
00053   //initialise mexstuff
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   //writing output
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   //free memory
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 

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