Scheduled service maintenance on November 22


On Friday, November 22, 2024, between 06:00 CET and 18:00 CET, GIN services will undergo planned maintenance. Extended service interruptions should be expected. We will try to keep downtimes to a minimum, but recommend that users avoid critical tasks, large data uploads, or DOI requests during this time.

We apologize for any inconvenience.

mexpp.c 3.8 KB

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