Mercurial > repos > vipints > rdiff

comparison rDiff/src/locfit/Source/mexpp.c @ 0:0f80a5141704

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
version 0.3 uploaded
author vipints
date Thu, 14 Feb 2013 23:38:36 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0f80a5141704
1 /*=================================================================
2 * mexpp.c
3 *
4 * starting a locfit interface.
5 *
6 * $Revision: 1.5 $ */
7 #include "mex.h"
8 #include "lfev.h"
9
10 extern void lfmxdata(), lfmxevs();
11
12 design des;
13 lfit lf;
14 int lf_error;
15
16 void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
17 { int i, vc, nc, mvc, mg[MXDIM], where, m, wh, dir, nkap;
18 const mxArray *dat, *fpc, *sp, *pc, *evs, *iwkc;
19 double *y, *x, *fy, *fx[MXDIM], *kap, cv, *fh, *se, *lo, *hi, lev;
20 char band[16], wher[16], what[16], rest[16];
21
22 if (nrhs<6) mexErrMsgTxt("ppmex requires 5 inputs.");
23 lf_error = 0;
24
25 mut_redirect(mexPrintf);
26 dat= mxGetField(prhs[1],0,"data");
27 evs= mxGetField(prhs[1],0,"evaluation_structure");
28 sp = mxGetField(prhs[1],0,"smoothing_parameters");
29 fpc= mxGetField(prhs[1],0,"fit_points");
30 pc = mxGetField(prhs[1],0,"parametric_component");
31 mxGetString(prhs[2],band,16);
32 mxGetString(prhs[3],what,16);
33 mxGetString(prhs[4],rest,16);
34 dir = mxGetPr(prhs[5])[0];
35 kap = mxGetPr(prhs[6]);
36 nkap = mxGetN(prhs[6]);
37 lev = mxGetPr(prhs[7])[0];
38 lfmxdata(&lf.lfd,dat);
39 lfmxsp(&lf.sp,sp,lf.lfd.d);
40 lfmxevs(&lf,evs); /* this has initialized module */
41
42 if (rest[0]=='n')
43 wh = ppwhat(what);
44 else
45 wh = 64+restyp(rest);
46
47 lf.fp.xev = mxGetPr(mxGetField(fpc,0,"evaluation_points"));
48 lf.lfd.d = lf.fp.d = mxGetM(mxGetField(fpc,0,"evaluation_points"));
49 lf.fp.nv = lf.fp.nvm = mxGetN(mxGetField(fpc,0,"evaluation_points"));
50 lf.fp.wk = mxGetPr(mxGetField(fpc,0,"fitted_values"));
51 lf.fp.h = NULL;
52 if (lf.mdl.alloc!=NULL) lf.mdl.alloc(&lf);
53
54 iwkc = mxGetField(fpc,0,"evaluation_vectors");
55 vc = mxGetM(mxGetField(iwkc,0,"cell"));
56 nc = mxGetN(mxGetField(iwkc,0,"cell"));
57 mvc = mxGetN(mxGetField(iwkc,0,"splitvar"));
58 lf.evs.iwk = mxCalloc(vc*nc+3*mvc,sizeof(int));
59 lf.evs.nce = lf.evs.ncm = nc;
60 y = mxGetPr(mxGetField(iwkc,0,"cell"));
61 lf.evs.ce = lf.evs.iwk;
62 for (i=0; i<vc*nc; i++) lf.evs.ce[i] = y[i];
63 y = mxGetPr(mxGetField(iwkc,0,"splitvar"));
64 lf.evs.s = &lf.evs.iwk[vc*nc];
65 for (i=0; i<mvc; i++) lf.evs.s[i] = y[i];
66 y = mxGetPr(mxGetField(iwkc,0,"lo"));
67 lf.evs.lo = &lf.evs.s[mvc];
68 for (i=0; i<mvc; i++) lf.evs.lo[i] = y[i];
69 y = mxGetPr(mxGetField(iwkc,0,"hi"));
70 lf.evs.hi = &lf.evs.lo[mvc];
71 for (i=0; i<mvc; i++) lf.evs.hi[i] = y[i];
72 lf.fp.hasd = deg(&lf.sp)>0;
73
74 lf.fp.kap = mxGetPr(mxGetField(fpc,0,"kappa"));
75
76 lf.pc.wk = mxGetPr(pc);
77 lf.pc.lwk = mxGetN(pc);
78 pcchk(&lf.pc,lf.fp.d,npar(&lf.sp),1);
79 haspc(&lf.pc) = !noparcomp(&lf.sp);
80 lf.pc.xtwx.st = JAC_EIGD;
81
82 where = 0;
83 if (mxIsCell(prhs[0])) /* interpret as grid margins */
84 { m = 1;
85 for (i=0; i<lf.fp.d; i++)
86 { fx[i] = mxGetPr(mxGetCell(prhs[0],i));
87 mg[i] = mxGetM(mxGetCell(prhs[0],i));
88 m *= mg[i];
89 }
90 where = 2;
91 }
92 if (mxIsChar(prhs[0]))
93 { mxGetString(prhs[0],wher,16);
94 where = 3; /* data */
95 m = mg[0] = lf.lfd.n;
96 if (wher[0] == 'f') /* or fit points */
97 { where = 4;
98 m = mg[0] = lf.fp.nv;
99 }
100 }
101 if (where==0) /* interpret as numeric vector */
102 { x = mxGetPr(prhs[0]);
103 m = mg[0] = mxGetM(prhs[0]); /* should check mxGetN == d */
104 for (i=0; i<lf.lfd.d; i++)
105 fx[i] = &x[i*m];
106 where=1;
107 }
108
109 plhs[0] = mxCreateDoubleMatrix(m,1,mxREAL); /* fitted values */
110 plhs[1] = mxCreateDoubleMatrix(m,1,mxREAL); /* std errors */
111 fh = mxGetPr(plhs[0]);
112 se = mxGetPr(plhs[1]);
113
114 preplot(&lf,fx,fh,se,band[0],mg,where,wh,dir);
115
116 cv = 0.0;
117 if (band[0] != 'n')
118 { cv = critval(1-lev,kap,nkap,lf.lfd.d,TWO_SIDED,0.0,GAUSS);
119 plhs[2] = mxCreateDoubleMatrix(m,2,mxREAL);
120 lo = mxGetPr(plhs[2]);
121 hi = &lo[m];
122 for (i=0; i<m; i++)
123 { lo[i] = fh[i]-cv*se[i];
124 hi[i] = fh[i]+cv*se[i];
125 }
126 }
127 else plhs[2] = mxCreateDoubleScalar(0.0);
128 }