diff rDiff/src/locfit/Source/mexlf.c @ 0:0f80a5141704

version 0.3 uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rDiff/src/locfit/Source/mexlf.c	Thu Feb 14 23:38:36 2013 -0500
@@ -0,0 +1,110 @@
+/*=================================================================
+ * mexlocfit.c 
+ *
+ * starting a locfit interface.
+ *
+/* $Revision: 1.5 $ */
+#include "mex.h"
+#include "lfev.h"
+
+design des;
+lfit lf;
+int lf_error;
+
+extern void lfmxdata(), lfmxsp(), lfmxevs();
+
+void
+mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
+{
+    int  d, i, nvc[5], nvm, vc, nv, nc, mvc, lw[7];
+    double *y;
+    mxArray *iwkc, *temparray;
+    const char * fpcnames[] = {"evaluation_points","fitted_values","evaluation_vectors","fit_limits","family_link","kappa"};
+    const char * evsnames[] = {"cell","splitvar","lo","hi"};
+
+    mut_redirect(mexPrintf);
+    if (nrhs != 3) mexErrMsgTxt("mexlf requires 3 inputs.");
+    lf_error = 0;
+
+    lfit_alloc(&lf);
+    lfit_init(&lf);
+
+    lfmxdata(&lf.lfd,prhs[0]);
+    d = lf.lfd.d;
+
+    lfmxsp(&lf.sp,prhs[2],d);
+    lfmxevs(&lf,prhs[1]);
+
+    guessnv(&lf.evs,&lf.sp,&lf.mdl,lf.lfd.n,d,lw,nvc);
+    nvm = lf.fp.nvm = nvc[0];
+    vc = nvc[2];
+    if (ev(&lf.evs) != EPRES) lf.fp.xev = mxCalloc(nvm*d,sizeof(double));
+    lf.fp.lev = nvm*d;
+    lf.fp.wk = lf.fp.coef = mxCalloc(lw[1],sizeof(double));
+    lf.fp.lwk = lw[1];
+    lf.evs.iwk = mxCalloc(lw[2],sizeof(int));
+    lf.evs.liw = lw[2];
+    plhs[1] = mxCreateDoubleMatrix(1,lw[3],mxREAL);
+    lf.pc.wk = mxGetPr(plhs[1]);
+    lf.pc.lwk = lw[3];
+    lf.fp.kap = mxCalloc(lw[5],sizeof(double));
+/* should also allocate design here */
+    
+    if (!lf_error) startmodule(&lf,&des,NULL,NULL);
+
+/* now, store the results: 
+   plhs[0] stores informtion about fit points and evaluation structure.
+    it is now a matlab structure, not a cell   
+   fit_points.evaluation_points - matrix of evaluation points.
+   fit_points.fitted_values - matrix fitted values etc.
+   fit_points.evaluation_vectors - structure of 'integer' vectors {cell,splitvar,lo,hi}
+   fit_points.fit_limits - fit limit (matrix with 2 cols).
+   fit_points.family_link - [familt link] numeric vector.
+   fit_points.kappa - kap vector.
+*/
+    plhs[0] = mxCreateStructMatrix(1,1,6,fpcnames);
+    if ( plhs[0] == NULL ) {
+       printf("Problem with CreateStructMatrix for plhs[0]\n");fflush(stdout);
+    }
+
+    mxSetField(plhs[0],0,"evaluation_points",mxCreateDoubleMatrix(d,lf.fp.nv,mxREAL));
+    memcpy(mxGetPr(mxGetField(plhs[0],0,"evaluation_points")), lf.fp.xev, d*lf.fp.nv*sizeof(double));
+
+    mxSetField(plhs[0],0,"fitted_values",mxCreateDoubleMatrix(lf.fp.nv,lf.mdl.keepv,mxREAL));
+    for (i=0; i<lf.mdl.keepv; i++)
+      memcpy(&mxGetPr(mxGetField(plhs[0],0,"fitted_values"))[i*lf.fp.nv], &lf.fp.coef[i*nvm], lf.fp.nv*sizeof(double));
+    /* another bit to save here? -- split values, kdtree */
+   
+    temparray = mxCreateStructMatrix(1,1,4,evsnames);
+    if ( temparray == NULL ) {
+      printf("Problem with CreateStructMatrix for temparray\n");fflush(stdout);
+    }
+    mxSetField(plhs[0],0,"evaluation_vectors",temparray);
+    iwkc = mxGetField(plhs[0],0,"evaluation_vectors");
+    nv = lf.fp.nv;
+    nc = lf.evs.nce;
+    mvc = (nv>nc) ? nv : nc;
+    mxSetField(iwkc,0,"cell",mxCreateDoubleMatrix(vc,nc,mxREAL));  /* ce */
+    mxSetField(iwkc,0,"splitvar",mxCreateDoubleMatrix(1,mvc,mxREAL));  /* s  */
+    mxSetField(iwkc,0,"lo",mxCreateDoubleMatrix(1,mvc,mxREAL));  /* lo */
+    mxSetField(iwkc,0,"hi",mxCreateDoubleMatrix(1,mvc,mxREAL));  /* hi */
+    y = mxGetPr(mxGetField(iwkc,0,"cell"));
+    for (i=0; i<vc*nc; i++) y[i] = lf.evs.ce[i];
+    y = mxGetPr(mxGetField(iwkc,0,"splitvar"));
+    for (i=0; i<mvc; i++) y[i] = lf.evs.s[i];
+    y = mxGetPr(mxGetField(iwkc,0,"lo"));
+    for (i=0; i<mvc; i++) y[i] = lf.evs.lo[i];
+    y = mxGetPr(mxGetField(iwkc,0,"hi"));
+    for (i=0; i<mvc; i++) y[i] = lf.evs.hi[i];
+
+    mxSetField(plhs[0],0,"fit_limits",mxCreateDoubleMatrix(d,2,mxREAL));
+    memcpy(mxGetPr(mxGetField(plhs[0],0,"fit_limits")), lf.evs.fl, 2*d*sizeof(double));
+    
+    mxSetField(plhs[0],0,"family_link",mxCreateDoubleMatrix(1,2,mxREAL));
+    y = mxGetPr(mxGetField(plhs[0],0,"family_link"));
+    y[0] = fam(&lf.sp);
+    y[1] = link(&lf.sp);
+
+    mxSetField(plhs[0],0,"kappa",mxCreateDoubleMatrix(1,lf.mdl.keepc,mxREAL));
+    memcpy(mxGetPr(mxGetField(plhs[0],0,"kappa")),lf.fp.kap,lf.mdl.keepc*sizeof(double));
+}