Mercurial > repos > vipints > rdiff
diff rDiff/src/locfit/Source/mexpp.c @ 0:0f80a5141704
version 0.3 uploaded
author | vipints |
---|---|
date | Thu, 14 Feb 2013 23:38:36 -0500 (2013-02-15) |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rDiff/src/locfit/Source/mexpp.c Thu Feb 14 23:38:36 2013 -0500 @@ -0,0 +1,128 @@ +/*================================================================= + * mexpp.c + * + * starting a locfit interface. + * + * $Revision: 1.5 $ */ +#include "mex.h" +#include "lfev.h" + +extern void lfmxdata(), lfmxevs(); + +design des; +lfit lf; +int lf_error; + +void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[]) +{ int i, vc, nc, mvc, mg[MXDIM], where, m, wh, dir, nkap; + const mxArray *dat, *fpc, *sp, *pc, *evs, *iwkc; + double *y, *x, *fy, *fx[MXDIM], *kap, cv, *fh, *se, *lo, *hi, lev; + char band[16], wher[16], what[16], rest[16]; + + if (nrhs<6) mexErrMsgTxt("ppmex requires 5 inputs."); + lf_error = 0; + + mut_redirect(mexPrintf); + dat= mxGetField(prhs[1],0,"data"); + evs= mxGetField(prhs[1],0,"evaluation_structure"); + sp = mxGetField(prhs[1],0,"smoothing_parameters"); + fpc= mxGetField(prhs[1],0,"fit_points"); + pc = mxGetField(prhs[1],0,"parametric_component"); + mxGetString(prhs[2],band,16); + mxGetString(prhs[3],what,16); + mxGetString(prhs[4],rest,16); + dir = mxGetPr(prhs[5])[0]; + kap = mxGetPr(prhs[6]); + nkap = mxGetN(prhs[6]); + lev = mxGetPr(prhs[7])[0]; + lfmxdata(&lf.lfd,dat); + lfmxsp(&lf.sp,sp,lf.lfd.d); + lfmxevs(&lf,evs); /* this has initialized module */ + + if (rest[0]=='n') + wh = ppwhat(what); + else + wh = 64+restyp(rest); + + lf.fp.xev = mxGetPr(mxGetField(fpc,0,"evaluation_points")); + lf.lfd.d = lf.fp.d = mxGetM(mxGetField(fpc,0,"evaluation_points")); + lf.fp.nv = lf.fp.nvm = mxGetN(mxGetField(fpc,0,"evaluation_points")); + lf.fp.wk = mxGetPr(mxGetField(fpc,0,"fitted_values")); + lf.fp.h = NULL; + if (lf.mdl.alloc!=NULL) lf.mdl.alloc(&lf); + + iwkc = mxGetField(fpc,0,"evaluation_vectors"); + vc = mxGetM(mxGetField(iwkc,0,"cell")); + nc = mxGetN(mxGetField(iwkc,0,"cell")); + mvc = mxGetN(mxGetField(iwkc,0,"splitvar")); + lf.evs.iwk = mxCalloc(vc*nc+3*mvc,sizeof(int)); + lf.evs.nce = lf.evs.ncm = nc; + y = mxGetPr(mxGetField(iwkc,0,"cell")); + lf.evs.ce = lf.evs.iwk; + for (i=0; i<vc*nc; i++) lf.evs.ce[i] = y[i]; + y = mxGetPr(mxGetField(iwkc,0,"splitvar")); + lf.evs.s = &lf.evs.iwk[vc*nc]; + for (i=0; i<mvc; i++) lf.evs.s[i] = y[i]; + y = mxGetPr(mxGetField(iwkc,0,"lo")); + lf.evs.lo = &lf.evs.s[mvc]; + for (i=0; i<mvc; i++) lf.evs.lo[i] = y[i]; + y = mxGetPr(mxGetField(iwkc,0,"hi")); + lf.evs.hi = &lf.evs.lo[mvc]; + for (i=0; i<mvc; i++) lf.evs.hi[i] = y[i]; + lf.fp.hasd = deg(&lf.sp)>0; + + lf.fp.kap = mxGetPr(mxGetField(fpc,0,"kappa")); + + lf.pc.wk = mxGetPr(pc); + lf.pc.lwk = mxGetN(pc); + pcchk(&lf.pc,lf.fp.d,npar(&lf.sp),1); + haspc(&lf.pc) = !noparcomp(&lf.sp); + lf.pc.xtwx.st = JAC_EIGD; + + where = 0; + if (mxIsCell(prhs[0])) /* interpret as grid margins */ + { m = 1; + for (i=0; i<lf.fp.d; i++) + { fx[i] = mxGetPr(mxGetCell(prhs[0],i)); + mg[i] = mxGetM(mxGetCell(prhs[0],i)); + m *= mg[i]; + } + where = 2; + } + if (mxIsChar(prhs[0])) + { mxGetString(prhs[0],wher,16); + where = 3; /* data */ + m = mg[0] = lf.lfd.n; + if (wher[0] == 'f') /* or fit points */ + { where = 4; + m = mg[0] = lf.fp.nv; + } + } + if (where==0) /* interpret as numeric vector */ + { x = mxGetPr(prhs[0]); + m = mg[0] = mxGetM(prhs[0]); /* should check mxGetN == d */ + for (i=0; i<lf.lfd.d; i++) + fx[i] = &x[i*m]; + where=1; + } + + plhs[0] = mxCreateDoubleMatrix(m,1,mxREAL); /* fitted values */ + plhs[1] = mxCreateDoubleMatrix(m,1,mxREAL); /* std errors */ + fh = mxGetPr(plhs[0]); + se = mxGetPr(plhs[1]); + + preplot(&lf,fx,fh,se,band[0],mg,where,wh,dir); + + cv = 0.0; + if (band[0] != 'n') + { cv = critval(1-lev,kap,nkap,lf.lfd.d,TWO_SIDED,0.0,GAUSS); + plhs[2] = mxCreateDoubleMatrix(m,2,mxREAL); + lo = mxGetPr(plhs[2]); + hi = &lo[m]; + for (i=0; i<m; i++) + { lo[i] = fh[i]-cv*se[i]; + hi[i] = fh[i]+cv*se[i]; + } + } + else plhs[2] = mxCreateDoubleScalar(0.0); +}