# to unbundle, sh this file (in an empty directory) echo README 1>&2 sed >README <<'//GO.SYSIN DD README' 's/^-//' -Software for Locally-Weighted Regression 18 August 1992 - -William S. Cleveland -Eric Grosse -Ming-Jen Shyu - -Locally-weighted regression, or loess, is a procedure for estimating a -regression surface by a multivariate smoothing procedure: fitting a -linear or quadratic function of the independent variables in a moving -fashion that is analogous to how a moving average is computed for a -time series. Compared to classical approaches - fitting global -parametric functions - loess substantially increases the domain of -surfaces that can be estimated without distortion. Also, a pleasant -fact about loess is that analogues of the statistical procedures used -in parametric function fitting - for example, ANOVA and t intervals - -involve statistics whose distributions are well approximated by -familiar distributions. - -The follwing files are included in this distribution. - README the instruction file you are reading now - S.h header file - air.c C source for air data example - changes history of changes to loess - depend.ps PostScript figure of how routines are related - ethanol.c C source for ethanol data example - galaxy.c C source for galaxy data example - gas.c C source for gas data example - loess.c C source (high-level loess routines) - loess.h header file for loess_struct and predict_struct - loess.m manual page for user-callable loess routines - loessc.c C source (low-level loess routines) - loessf.f FORTRAN source (low-level loess & predict routines) - loessf.m documentation for FORTRAN source - madeup.c C source for madeup data example - makefile makefile to compile the example codes - misc.c C source (anova, pointwise, and other support routines) - predict.c C source (high-level predict routines) - predict.m manual page for user-callable predict routines - struct.m manual page for loess_struct, pred_struct - supp.f supplemental Fortran loess drivers - -After unpacking these files, just type "make" and if all goes well -you should see output like: - - loess(&gas): - Number of Observations: 22 - Equivalent Number of Parameters: 5.5 - Residual Standard Error: 0.3404 - - loess(&gas_null): - Number of Observations: 22 - Equivalent Number of Parameters: 3.5 - Residual Standard Error: 0.5197 - - predict(gas_fit_E, m, &gas, &gas_pred): - 1.19641 5.06875 0.523682 - - pointwise(&gas_pred, m, coverage, &gas_ci): - 1.98562 4.10981 5.48023 5.56651 3.52761 1.71062 1.47205 - 1.19641 3.6795 5.05571 5.13526 3.14366 1.19693 0.523682 - 0.407208 3.24919 4.63119 4.70401 2.7597 0.683247 -0.424684 - - anova(&gas_null, &gas, &gas_anova): - 2.5531 15.663 10.1397 0.000860102 - -To run other examples, simply type "make galaxy", or "make ethanol", etc. - -If your loader complains about "-llinpack -lblas" in the makefile, change -it to whatever your system prefers for accessing Linpack and the Blas. -If necessary, these Fortran subroutines can be obtained by - mail netlib@netlib.bell-labs.com - send dnrm2 dsvdc dqrdc ddot dqrsl idamax from linpack core. - -A 50 page user guide, in PostScript form, is available by anonymous ftp. - ftp netlib.bell-labs.com - login: anonymous - password: - binary - cd /netlib/a - get cloess.ps.Z - quit - uncompress cloess.ps -This guide describes crucial steps in the proper analysis of data using -loess. Please read it. - -Bug reports are appreciated. Send electronic mail to - ehg@netlib.bell-labs.com -including the words "this is not spam" in the Subject line -or send paper mail to - Eric Grosse - Bell Labs 2T-502 - Murray Hill NJ 07974 -for problems with the Fortran inner core of the algorithm. -The C drivers were written by Ming-Jen Shyu, who left Bell Labs. Eric will -fix problems with them when he can. - -Remember that this is experimental software distributed free of charge -and comes with no warranty! Exercise professional caution. - -Happy Smoothing! - -/* - * The authors of this software are Cleveland, Grosse, and Shyu. - * Copyright (c) 1989, 1992 by AT&T. - * Permission to use, copy, modify, and distribute this software for any - * purpose without fee is hereby granted, provided that this entire notice - * is included in all copies of any software which is or includes a copy - * or modification of this software and in all copies of the supporting - * documentation for such software. - * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED - * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY - * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY - * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. - */ - //GO.SYSIN DD README echo S.h 1>&2 sed >S.h <<'//GO.SYSIN DD S.h' 's/^-//' -#include -#include - -#define Calloc(n,t) (t *)calloc((unsigned)(n),sizeof(t)) -#define Free(p) free((char *)(p)) - -/* the mapping from f77 to C intermediate code -- may be machine dependent - * the first definition satisfies lint's narrowminded preprocessing & should - * stay the same for all implementations. The __STDC__ definition is for - * ANSI standard conforming C compilers. The #else definition should - * generate the version of the fortran subroutine & common block names x - * handed to the local loader; e.g., "x_" in system V, Berkeley & 9th edition - */ - -#ifdef lint -#define F77_SUB(x) x -#define F77_COM(x) x -#else -#ifdef __STDC__ -#define F77_SUB(x) x##_ -#define F77_COM(x) x##_ -#else -#define F77_SUB(x) x/**/_ -#define F77_COM(x) x/**/_ -#endif -#endif - -#define NULL_ENTRY ((int *)NULL) - //GO.SYSIN DD S.h echo air.c 1>&2 sed >air.c <<'//GO.SYSIN DD air.c' 's/^-//' -#include -#include "loess.h" - -struct loess_struct air; -double ozone[] = {3.44821724038273, 3.30192724889463, 2.28942848510666, - 2.6207413942089, 2.84386697985157, 2.66840164872194, 2, - 2.51984209978975, 2.22398009056931, 2.41014226417523, - 2.6207413942089, 2.41014226417523, 3.23961180127748, - 1.81712059283214, 3.10723250595386, 2.22398009056931, 1, - 2.22398009056931, 1.5874010519682, 3.1748021039364, - 2.84386697985157, 3.55689330449006, 4.86294413109428, - 3.33222185164595, 3.07231682568585, 4.14081774942285, - 3.39121144301417, 2.84386697985157, 2.75892417638112, - 3.33222185164595, 2.71441761659491, 2.28942848510666, - 2.35133468772076, 5.12992784003009, 3.65930571002297, - 3.1748021039364, 4, 3.41995189335339, 4.25432086511501, - 4.59470089220704, 4.59470089220704, 4.39682967215818, - 2.15443469003188, 3, 1.91293118277239, 3.63424118566428, - 3.27106631018859, 3.93649718310217, 4.29084042702621, - 3.97905720789639, 2.51984209978975, 4.30886938006377, - 4.7622031559046, 2.71441761659491, 3.73251115681725, - 4.34448148576861, 3.68403149864039, 4, 3.89299641587326, - 3.39121144301417, 2.0800838230519, 2.51984209978975, - 4.9596756638423, 4.46474509558454, 4.79141985706278, - 3.53034833532606, 3.03658897187566, 4.02072575858906, - 2.80203933065539, 3.89299641587326, 2.84386697985157, - 3.14138065239139, 3.53034833532606, 2.75892417638112, - 2.0800838230519, 3.55689330449006, 5.51784835276224, - 4.17933919638123, 4.23582358425489, 4.90486813152402, - 4.37951913988789, 4.39682967215818, 4.57885697021333, - 4.27265868169792, 4.17933919638123, 4.49794144527541, - 3.60882608013869, 3.1748021039364, 2.71441761659491, - 2.84386697985157, 2.75892417638112, 2.88449914061482, - 3.53034833532606, 2.75892417638112, 3.03658897187566, - 2.0800838230519, 2.35133468772076, 3.58304787101595, - 2.6207413942089, 2.35133468772076, 2.88449914061482, - 2.51984209978975, 2.35133468772076, 2.84386697985157, - 3.30192724889463, 1.91293118277239, 2.41014226417523, - 3.10723250595386, 2.41014226417523, 2.6207413942089, - 2.71441761659491}; -double rad_temp_wind[] = {190, 118, 149, 313, 299, 99, 19, 256, 290, 274, 65, - 334, 307, 78, 322, 44, 8, 320, 25, 92, 13, 252, 223, 279, 127, - 291, 323, 148, 191, 284, 37, 120, 137, 269, 248, 236, 175, - 314, 276, 267, 272, 175, 264, 175, 48, 260, 274, 285, 187, - 220, 7, 294, 223, 81, 82, 213, 275, 253, 254, 83, 24, 77, 255, - 229, 207, 192, 273, 157, 71, 51, 115, 244, 190, 259, 36, 212, - 238, 215, 203, 225, 237, 188, 167, 197, 183, 189, 95, 92, 252, - 220, 230, 259, 236, 259, 238, 24, 112, 237, 224, 27, 238, 201, - 238, 14, 139, 49, 20, 193, 191, 131, 223, - 67, 72, 74, 62, 65, 59, 61, 69, 66, 68, 58, 64, 66, 57, 68, - 62, 59, 73, 61, 61, 67, 81, 79, 76, 82, 90, 87, 82, 77, 72, - 65, 73, 76, 84, 85, 81, 83, 83, 88, 92, 92, 89, 73, 81, 80, - 81, 82, 84, 87, 85, 74, 86, 85, 82, 86, 88, 86, 83, 81, 81, - 81, 82, 89, 90, 90, 86, 82, 80, 77, 79, 76, 78, 78, 77, 72, - 79, 81, 86, 97, 94, 96, 94, 91, 92, 93, 93, 87, 84, 80, 78, - 75, 73, 81, 76, 77, 71, 71, 78, 67, 76, 68, 82, 64, 71, 81, - 69, 63, 70, 75, 76, 68, - 7.4, 8, 12.6, 11.5, 8.6, 13.8, 20.1, 9.7, 9.2, 10.9, 13.2, - 11.5, 12, 18.4, 11.5, 9.7, 9.7, 16.6, 9.7, 12, 12, 14.9, 5.7, - 7.4, 9.7, 13.8, 11.5, 8, 14.9, 20.7, 9.2, 11.5, 10.3, 4, 9.2, - 9.2, 4.6, 10.9, 5.1, 6.3, 5.7, 7.4, 14.3, 14.9, 14.3, 6.9, - 10.3, 6.3, 5.1, 11.5, 6.9, 8.6, 8, 8.6, 12, 7.4, 7.4, 7.4, - 9.2, 6.9, 13.8, 7.4, 4, 10.3, 8, 11.5, 11.5, 9.7, 10.3, 6.3, - 7.4, 10.9, 10.3, 15.5, 14.3, 9.7, 3.4, 8, 9.7, 2.3, 6.3, 6.3, - 6.9, 5.1, 2.8, 4.6, 7.4, 15.5, 10.9, 10.3, 10.9, 9.7, 14.9, - 15.5, 6.3, 10.9, 11.5, 6.9, 13.8, 10.3, 10.3, 8, 12.6, 9.2, - 10.3, 10.3, 16.6, 6.9, 14.3, 8, 11.5}; -long n = 111, p = 3; - -main() { - printf("\nloess(&air):\n"); - loess_setup(rad_temp_wind, ozone, n, p, &air); - air.model.span = 0.8; - loess(&air); - loess_summary(&air); - - loess_free_mem(&air); -} //GO.SYSIN DD air.c echo changes 1>&2 sed >changes <<'//GO.SYSIN DD changes' 's/^-//' -CHANGES PLANNED SOMEDAY -1) more vertices in k-d tree for dimension > 2, to get continuity. -2) triangulation based method. ----------------------- - -19 Nov 1987 workspace not big enough for degree=2 - -22 Jan 1988 switched from depth first to breadth first tree build - -14 Mar 1988 lostt.3 extra space needed if (method mod 1000 = 0), - not the documented (method/1000=0) - -28 Apr 1988 l2tr.g vval2 needed to be initialized to 0 - - galaxy smooth needs double precision on vax - -26 May 1988 bbox.g add 10% margin to allow limited extrapolation - -6 June 1988 loess/lostt.f trL wasn't set if method/1000==0 - -10 June 1988 losave, loread - - v(RCOND) 1 / max condition number - -12 June 1988 lofort - -21 June 1988 additional workspace for explicit L - -27 June 1988 workspace checking in lowesf was slightly pessimistic - -30 June 1988 Changed default fdiam to 0. - Added warning messages for memory limits and pseudoinverse. - -4 Aug 1988 bbox.g changed margin from 10% to 0.5%. - -24 Aug 1988 loser documentation should have specified workspace - of size ...+m*n, not ...+m**2. - -Sep 1988 - loess-based approximations of delta1,2. - pseudo-values, so statistics are available with robustness iterations. - reorganize error messages to better fit into S. - sample driver program. - somewhat shorter code generated by ehg170. - -20 Dec 1988 - workspace in loser - -27 Jan 1989 - workspace checking in lostt was a bit pessimistic. - -3 Feb 1989 - l2fit, l2tr: error message should contain sqrt(rho) - -18 Dec 1989 - ehg141, ehg179-ehg181: new delta approximations - -24 Jan 1990 - master copy moved from Sun3/180 to SGI 4D/240S - (no intentional changes) - -25 Jan 1990 - (many routines touched; ehg127 added) cleaned up computational - kernel, added provision for only first dd<=d variables to enter - the distance calculation ("conditionally parametric variables"), - added independent bounds on total and componentwise degree for - local polynomial model, made extrapolation warning message print - a bit more detail. - -14 Mar 1990 - added setLf argument to lowesd; added lowesr, lowesl for resmoothing. - -------------------------------------------------------- -Converting to the new version of loess -5 April 1990 - -Over the past few months, a number of changes have been made to the -loess package, to provide more control over the local model, to allow -conditionally parametric variables, and to return exact statistical -quantities for the blending method. Unlike earlier internal -algorithmic improvements, this round of changes added some extra -arguments in the Fortran calling sequences. The purpose of this note -is to assist in converting programs that called the old version. - -An explicit argument setLf has been added to lowesd(), since it affects -the partitioning of the workspace. To help protect against inadvertent -version mismatches, the version number that lowesd() checks has also -been changed. The componentwise degree and the specification of -conditionally nonparametric variables can be changed from the default -by modifying iv(CDEG) and iv(NDIST). - -The influence matrix L for blending is now explicitly available by -calling a new subroutine lowesl(), but this loses the speed -advantage of blending. A faster, sometimes equivalent method is -to use the influence matrix that carries data values to coefficients -at the vertices of the k-d tree. This information is saved in iv(iv(Lq)) -and v(iv(Lf)), for the afficionado. - -The new subroutine lowesr() takes advantage of Lq and Lf to allow rapid -resmoothing for applications when only y, not x, is subject to change. -------------------------------------------------------- - -7 May 1990 - new delta approximations. - added prior weights to input format for sample driver. - -29 May 1990 - loess,lostt,loser,pseudo moved from Fortran to S. - -11 Jul 1990 - column equilibration, so pseudoinverse is needed less often. - -27 May 1991 -lowesd version 105; increased nvmax,ncmax to max(200,n). -l2fit added ihat=1 (diagL only). -ehg133,lowese removed unused arguments dist,eta. -ehg190,ehg141 changed name to lowesa, slight change to calling sequence. -ehg144 changed name to lowesc -m9rwt changed name to lowesw -pseudo changed name to lowesp - -22 Jul 1991 IMPORTANT BUG FIX! -ehg131 vval2 should be dimensioned 0:d, not 0:8 - -26 Jul 1991 -lowesd change calling sequence to provide tighter memory allocation -diff old/man/internal new/man/internal -< lowesd(105,iv,liv,lv,v,d,n,f,tdeg,setLf) setup workspace -> lowesd(106,iv,liv,lv,v,d,n,f,tdeg,nvmax,setLf) setup workspace -< liv 50+(2^d+6)*max(200,n) -< if setLf, add nf*max(200,n) -< lv 50+(3*d+4)*max(200,n)+(tau+2)*nf -< if setLf, add (d+1)*nf*max(200,n) -> liv 50+(2^d+6)*nvmax -> if setLf, add nf*nvmax -> lv 50+(3*d+4)*nvmax+(tau+2)*nf -> if setLf, add (d+1)*nf*nvmax -> nvmax limit on number of vertices for kd-tree; e.g. max(200,n) - -20 Sep 1991 -sample.f brought in sync with recent loess changes. - -24 Dec 1991 -l2fit.f fixed comment in single precision version - -10 Jan 1992 -ehg197.f new formula for approximating trL, valid for small f - -15 May 1992 -netlib/a/dloess now includes C drivers (written by Ming-Jen Shyu, - adapted from code used inside the S system) - -22 Jun 1992 -ehg191.f Loop 11 ran too far, picking up one more value than necessary. - The value was not used, so the loess computation itself is unaffected, - but on some systems the old code could conceivably cause a reference - to an invalid memory address and abort with a segmentation fault - message. - -23 Jun 1992 -S.h #include , since loessc.c calls floor() and pow(). - -18 Aug 1992 -netlib/a/dloess A new release with bug fixes in all the C drivers, new - example codes, and detail documentations. - -25 Mar 1996 -predict.c fix enormous memory leak. update email address //GO.SYSIN DD changes echo depend.ps 1>&2 sed >depend.ps <<'//GO.SYSIN DD depend.ps' 's/^-//' -%! -/Courier-Bold findfont 10 scalefont setfont -%draw a box -%x y width height box -/box { newpath - /height exch def - /width exch def - /y exch def - /x exch def - x width 2 div sub - y height 2 div sub moveto - width 0 rlineto - 0 height rlineto - width neg 0 rlineto - closepath } def - -%draw a circle -%x y radius circle -/circle { newpath 0 360 arc } def - -%draw an ellipse -%x y width height ellipse -/ellipse { gsave - /height exch def - /width exch def - 1 height width div scale - width height div mul - width 2 div - circle stroke - grestore } def - -%draw a centered label -%x y str -/label { - /str exch def - /y exch def - /x exch def - str stringwidth - pop /width exch def - x width 2 div sub - y 10 3 div sub moveto str show - } def - -%draw a line -%x1 y1 x2 y2 drawline -/drawline { 4 -2 roll moveto lineto stroke } def - -277 684 42 14 box stroke -277 684 (lowesd) label -349 630 42 14 box stroke -349 630 (lowesf) label -205 630 42 14 box stroke -205 630 (lowesb) label -155 565 42 14 box stroke -155 565 (lowesr) label -146 427 42 14 box stroke -146 427 (lowese) label -277 576 42 14 box stroke -277 576 (lowesl) label -203 464 42 14 box stroke -203 464 (lofort) label -81 576 42 14 box stroke -81 576 (losave) label -81 522 42 14 box stroke -81 522 (lohead) label -81 468 42 14 box stroke -81 468 (loread) label -405 540 42 14 box stroke -405 540 (lowesa) label -342 539 42 14 box stroke -342 539 (lowesc) label -92 461 134 434 drawline -124.266363 435.502104 134.000000 434.000000 drawline -134.000000 434.000000 128.592424 442.231532 drawline -81 515 81 475 drawline -77.000000 484.000000 81.000000 475.000000 drawline -81.000000 475.000000 85.000000 484.000000 drawline -81 569 81 529 drawline -77.000000 538.000000 81.000000 529.000000 drawline -81.000000 529.000000 85.000000 538.000000 drawline -289 569 329 546 drawline -319.203959 547.018615 329.000000 546.000000 drawline -329.000000 546.000000 323.191728 553.953865 drawline -154 558 146 434 drawline -142.587739 443.238857 146.000000 434.000000 drawline -146.000000 434.000000 150.571142 442.723799 drawline -188 623 97 583 drawline -103.629564 590.283466 97.000000 583.000000 drawline -97.000000 583.000000 106.848776 582.959760 drawline -204 623 203 471 drawline -199.059296 480.026120 203.000000 471.000000 drawline -203.000000 471.000000 207.059123 479.973490 drawline -214 623 267 583 drawline -257.406670 585.228906 267.000000 583.000000 drawline -267.000000 583.000000 262.225925 591.614419 drawline -199 623 160 572 drawline -162.289620 581.579021 160.000000 572.000000 drawline -160.000000 572.000000 168.644482 576.719420 drawline -220 623 389 547 drawline -379.151237 547.043173 389.000000 547.000000 drawline -389.000000 547.000000 382.432359 554.339352 drawline -202 623 148 434 drawline -146.626394 443.752600 148.000000 434.000000 drawline -148.000000 434.000000 154.318586 441.554831 drawline -348 623 342 546 drawline -338.711268 555.283547 342.000000 546.000000 drawline -342.000000 546.000000 346.687091 554.662054 drawline -353 623 400 547 drawline -391.864262 552.550655 400.000000 547.000000 drawline -400.000000 547.000000 398.668290 556.758409 drawline -267 677 214 637 drawline -218.774075 645.614419 214.000000 637.000000 drawline -214.000000 637.000000 223.593330 639.228906 drawline -286 677 339 637 drawline -329.406670 639.228906 339.000000 637.000000 drawline -339.000000 637.000000 334.225925 645.614419 drawline -showpage //GO.SYSIN DD depend.ps echo ethanol.c 1>&2 sed >ethanol.c <<'//GO.SYSIN DD ethanol.c' 's/^-//' -#include -#include "loess.h" - -struct loess_struct ethanol, ethanol_cp; -struct pred_struct ethanol_pred, ethanol_grid; -struct ci_struct ethanol_ci; -double NOx[] = {3.741, 2.295, 1.498, 2.881, 0.76, 3.12, 0.638, 1.17, 2.358, - 0.606, 3.669, 1, 0.981, 1.192, 0.926, 1.59, 1.806, 1.962, - 4.028, 3.148, 1.836, 2.845, 1.013, 0.414, 0.812, 0.374, 3.623, - 1.869, 2.836, 3.567, 0.866, 1.369, 0.542, 2.739, 1.2, 1.719, - 3.423, 1.634, 1.021, 2.157, 3.361, 1.39, 1.947, 0.962, 0.571, - 2.219, 1.419, 3.519, 1.732, 3.206, 2.471, 1.777, 2.571, 3.952, - 3.931, 1.587, 1.397, 3.536, 2.202, 0.756, 1.62, 3.656, 2.964, - 3.76, 0.672, 3.677, 3.517, 3.29, 1.139, 0.727, 2.581, 0.923, - 1.527, 3.388, 2.085, 0.966, 3.488, 0.754, 0.797, 2.064, 3.732, - 0.586, 0.561, 0.563, 0.678, 0.37, 0.53, 1.9}; -double C_E[] = {12, 12, 12, 12, 12, 9, 9, 9, 12, 12, 12, 12, 15, 18, 7.5, 12, - 12, 15, 15, 9, 9, 7.5, 7.5, 18, 18, 15, 15, 7.5, 7.5, 9, 15, 15, - 15, 15, 15, 9, 9, 7.5, 7.5, 7.5, 18, 18, 18, 18, 9, 9, 9, 9, - 7.5, 7.5, 7.5, 15, 18, 18, 15, 15, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, - 7.5, 18, 18, 18, 12, 12, 9, 9, 9, 15, 15, 15, 15, 15, 7.5, 7.5, - 9, 7.5, 18, 18, 7.5, 9, 12, 15, 18, 18, - 0.907, 0.761, 1.108, 1.016, 1.189, 1.001, 1.231, 1.123, 1.042, - 1.215, 0.93, 1.152, 1.138, 0.601, 0.696, 0.686, 1.072, 1.074, - 0.934, 0.808, 1.071, 1.009, 1.142, 1.229, 1.175, 0.568, 0.977, - 0.767, 1.006, 0.893, 1.152, 0.693, 1.232, 1.036, 1.125, 1.081, - 0.868, 0.762, 1.144, 1.045, 0.797, 1.115, 1.07, 1.219, 0.637, - 0.733, 0.715, 0.872, 0.765, 0.878, 0.811, 0.676, 1.045, 0.968, - 0.846, 0.684, 0.729, 0.911, 0.808, 1.168, 0.749, 0.892, 1.002, - 0.812, 1.23, 0.804, 0.813, 1.002, 0.696, 1.199, 1.03, 0.602, - 0.694, 0.816, 1.037, 1.181, 0.899, 1.227, 1.18, 0.795, 0.99, - 1.201, 0.629, 0.608, 0.584, 0.562, 0.535, 0.655}; -double newdata[] = {7.5, 9.0, 12.0, 15.0, 18.0, 0.6, 0.8, 1.0, 0.8, 0.6}; -double Cmin = 7.5, Cmax = 18.0, Emin = 0.535, Emax = 1.232; -double Cm[7], Em[16], grid[224]; -double tmp, coverage = .99; -long n = 88, p = 2, m = 5, se_fit = FALSE; -int i, j, k; - -main() { - printf("\nloess(ðanol): (span = 0.5)\n"); - loess_setup(C_E, NOx, n, p, ðanol); - ethanol.model.span = 0.5; - loess(ðanol); - loess_summary(ðanol); - - printf("\nloess(ðanol): (span = 0.25)\n"); - ethanol.model.span = 0.25; - loess(ðanol); - loess_summary(ðanol); - - printf("\nloess(ðanol_cp): (span = 0.25)\n"); - loess_setup(C_E, NOx, n, p, ðanol_cp); - ethanol_cp.model.span = 0.25; - ethanol_cp.model.parametric[0] = TRUE; - ethanol_cp.model.drop_square[0] = TRUE; - loess(ðanol_cp); - loess_summary(ðanol_cp); - - printf("\nloess(ðanol_cp): (span = 0.5)\n"); - ethanol_cp.model.span = 0.5; - loess(ðanol_cp); - loess_summary(ðanol_cp); - - printf("\npredict(newdata, m, ðanol, ðanol_pred, %d):\n", se_fit); - predict(newdata, m, ðanol_cp, ðanol_pred, se_fit); - for(i = 0; i < m; i++) - printf("%g ", ethanol_pred.fit[i]); - printf("\n"); - - m = 112; - se_fit = TRUE; - tmp = (Cmax - Cmin) / 6; - for(i = 0; i < 7; i++) - Cm[i] = Cmin + tmp * i; - tmp = (Emax - Emin) / 15; - for(i = 0; i < 16; i++) - Em[i] = Emin + tmp * i; - for(i = 0; i < 16; i++) { - k = i * 7; - for(j = 0; j < 7; j++) { - grid[k + j] = Cm[j]; - grid[m + k + j] = Em[i]; - } - } - predict(grid, m, ðanol_cp, ðanol_grid, se_fit); - pointwise(ðanol_grid, m, coverage, ðanol_ci); - - loess_free_mem(ðanol); - loess_free_mem(ðanol_cp); - pred_free_mem(ðanol_pred); - pred_free_mem(ðanol_grid); -} //GO.SYSIN DD ethanol.c echo galaxy.c 1>&2 sed >galaxy.c <<'//GO.SYSIN DD galaxy.c' 's/^-//' -#include -#include "loess.h" - -struct loess_struct galaxy; -struct pred_struct galaxy_contour, spine_fit, spine_se; -struct ci_struct spine_ci; -double velocity[] = {1769, 1749, 1749, 1758, 1750, 1745, 1750, 1753, 1734, - 1710, 1711, 1709, 1674, 1665, 1680, 1648, 1626, 1581, 1602, - 1558, 1538, 1506, 1521, 1498, 1501, 1491, 1481, 1468, 1455, - 1454, 1456, 1459, 1451, 1465, 1451, 1486, 1433, 1631, 1618, - 1607, 1608, 1601, 1603, 1612, 1607, 1618, 1649, 1595, 1580, - 1574, 1574, 1559, 1578, 1591, 1579, 1588, 1581, 1569, 1572, - 1584, 1565, 1718, 1711, 1710, 1715, 1713, 1717, 1715, 1712, - 1710, 1692, 1669, 1679, 1691, 1647, 1630, 1616, 1576, 1561, - 1558, 1538, 1525, 1509, 1501, 1494, 1489, 1493, 1487, 1495, - 1511, 1505, 1508, 1507, 1513, 1493, 1495, 1736, 1744, 1765, - 1766, 1764, 1715, 1751, 1761, 1763, 1758, 1743, 1738, 1732, - 1734, 1723, 1706, 1665, 1677, 1679, 1601, 1629, 1621, 1574, - 1559, 1540, 1525, 1517, 1506, 1481, 1465, 1468, 1465, 1454, - 1448, 1441, 1441, 1430, 1434, 1445, 1464, 1471, 1442, 1436, - 1434, 1428, 1558, 1563, 1581, 1548, 1572, 1574, 1578, 1576, - 1583, 1584, 1566, 1568, 1577, 1587, 1606, 1593, 1584, 1595, - 1617, 1552, 1597, 1615, 1626, 1626, 1586, 1624, 1600, 1585, - 1738, 1690, 1729, 1719, 1702, 1754, 1741, 1736, 1731, 1725, - 1710, 1673, 1669, 1641, 1675, 1681, 1645, 1594, 1583, 1599, - 1578, 1548, 1543, 1537, 1543, 1519, 1500, 1488, 1486, 1483, - 1481, 1485, 1480, 1479, 1505, 1482, 1481, 1489, 1531, 1533, - 1539, 1526, 1551, 1549, 1532, 1538, 1550, 1536, 1519, 1536, - 1535, 1536, 1533, 1528, 1539, 1546, 1552, 1557, 1573, 1553, - 1576, 1591, 1591, 1624, 1633, 1597, 1605, 1629, 1658, 1664, - 1667, 1671, 1687, 1682, 1668, 1673, 1684, 1668, 1618, 1658, - 1644, 1647, 1642, 1616, 1629, 1610, 1603, 1613, 1603, 1606, - 1603, 1608, 1613, 1616, 1615, 1611, 1580, 1580, 1586, 1591, - 1592, 1562, 1572, 1589, 1588, 1585, 1586, 1573, 1573, 1558, - 1566, 1740, 1704, 1748, 1757, 1775, 1765, 1762, 1752, 1752, - 1753, 1753, 1748, 1730, 1709, 1688, 1687, 1678, 1654, 1634, - 1611, 1590, 1562, 1565, 1541, 1537, 1515, 1498, 1479, 1481, - 1475, 1466, 1461, 1457, 1455, 1452, 1453, 1448, 1469, 1456, - 1448, 1409, 1416, 1429}; -double direction[] = {8.46279, 7.96498, 7.46717, 6.96936, 6.47154, 5.97373, - 5.47592, 4.97811, 4.4803, 3.98249, 3.46303, 2.96522, - 2.46741, 1.9696, 1.47179, 0.973978, 0.476167, -0.021644, - -0.519455, -1.01727, -1.51508, -2.01289, -2.5107, - -3.00851, -3.52797, -4.02578, -4.52359, -5.0214, - -5.51921, -6.01702, -6.51483, -7.01264, -7.51045, - -8.00827, -8.50608, -9.5017, -11.0168, 27.8244, 21.088, - 18.8425, 16.597, 14.3516, 12.1061, 9.86059, 7.61511, - 5.272, 3.02652, 0.781037, -1.46444, -3.70992, -5.95541, - -8.20089, -10.4464, -12.6918, -14.9373, -17.1828, - -19.4283, -21.6738, -23.9193, -26.2624, -28.5078, - 23.8699, 22.3013, 20.7327, 19.1642, 17.5956, 16.027, - 14.3902, 12.8216, 11.253, 9.68438, 8.11578, 6.54718, - 4.97859, 3.40999, 1.8414, 0.272799, -1.2958, -2.86439, - -4.43299, -6.00159, -7.63838, -9.20698, -10.7756, - -12.3442, -13.9128, -15.4814, -17.05, -18.6186, - -20.1872, -21.7557, -23.3243, -24.8929, -26.4615, - -28.0301, -29.6669, 18.4201, 17.5959, 16.7716, 15.9474, - 14.263, 13.4388, 12.6146, 11.7903, 10.9661, 10.1418, - 9.31757, 8.49332, 7.66907, 6.84483, 6.02058, 5.19634, - 4.37209, 3.54784, 2.68776, 1.86351, 1.03927, 0.215021, - -0.609226, -1.43347, -2.25772, -3.08196, -3.90621, - -4.73046, -5.5547, -6.37895, -7.2032, -8.02744, - -8.88752, -9.71177, -10.536, -11.3603, -12.1845, - -13.0088, -13.833, -14.6572, -15.4815, -16.3057, - -17.13, -17.9542, -18.7785, 25.8899, 24.2078, 22.4526, - 20.8436, 19.1615, 17.4794, 15.7972, 14.1151, 12.433, - 10.7509, 9.06879, 7.31354, 5.70456, 3.94931, 2.19406, - 0.511948, -1.09703, -2.77914, -4.46126, -6.07024, - -7.82548, -9.5076, -11.1897, -12.8718, -14.5539, - -16.2361, -23.1108, -24.7198, 1.97596, 1.77531, 1.67498, - 1.57466, 1.47434, 1.37401, 1.27369, 1.17336, 1.07304, - 0.972712, 0.872388, 0.767701, 0.667377, 0.567052, - 0.466727, 0.366403, 0.266078, 0.165754, 0.0654291, - -0.0348955, -0.13522, -0.235545, -0.335869, -0.436194, - -0.536518, -0.636843, -0.74153, -0.841854, -0.942179, - -1.0425, -1.14283, -1.24315, -1.34348, -1.4438, - -1.54413, -1.64445, -1.74478, -1.8451, 24.8532, 23.827, - 22.8007, 21.7298, 20.7036, 19.6773, 18.6511, 16.5539, - 15.5723, 14.546, 13.4752, 12.4489, 11.4227, 10.3964, - 9.37015, 8.3439, 7.31764, 6.29139, 5.26513, 4.23888, - 3.21262, 2.18637, 1.16011, 0.133859, -0.937015, - -1.96327, -2.98953, -4.01578, -5.04204, -6.06829, - -7.04993, -8.07618, -9.14706, -10.1733, -11.1996, - -12.2258, -13.2521, -14.2783, -15.3046, -16.3308, - -17.3571, -18.3834, -19.4096, -20.4359, -21.4621, - -22.4884, 29.4841, 27.0434, 25.0908, 22.6501, 20.4046, - 18.1591, 15.9136, 13.7658, 11.4227, 9.17718, 6.9317, - 4.58859, 2.44074, 0.0976296, -2.05022, -4.19807, - -6.63881, -8.88429, -11.1298, -13.2776, -15.5231, - -17.8662, -20.1117, -22.3572, -24.6027, -26.8481, - -29.0936, 10.8869, 9.39348, 8.91731, 8.39786, 7.92169, - 7.42388, 6.92607, 6.42826, 5.9088, 5.41099, 4.91318, - 4.41537, 3.91756, 3.44139, 2.92193, 2.42412, 1.92631, - 1.4285, 0.93069, 0.432879, -0.0649319, -0.562743, - -1.06055, -1.55837, -2.07782, -2.55399, -3.07344, - -3.57125, -4.06906, -4.56688, -5.06469, -5.5625, - -6.06031, -6.55812, -7.05593, -7.57539, -8.0732, - -8.54937, -9.09046, -9.58827, -10.0428, -10.5406, - -11.0601, - -38.1732, -35.9277, -33.6822, -31.4367, -29.1912, - -26.9458, -24.7003, -22.4548, -20.2093, -17.9638, -15.6207, - -13.3753, -11.1298, -8.88429, -6.63881, -4.39333, -2.14785, - 0.0976296, 2.34311, 4.58859, 6.83407, 9.07955, 11.325, - 13.5705, 15.9136, 18.1591, 20.4046, 22.6501, 24.8955, - 27.141, 29.3865, 31.632, 33.8775, 36.123, 38.3684, 42.8594, - 49.6935, 6.16853, 4.6751, 4.17728, 3.67947, 3.18166, 2.68385, - 2.18604, 1.68823, 1.16877, 0.670963, 0.173152, -0.324659, - -0.822471, -1.32028, -1.81809, -2.3159, -2.81371, -3.31153, - -3.80934, -4.30715, -4.80496, -5.30277, -5.82223, -6.32004, - -25.5974, -23.9153, -22.2332, -20.551, -18.8689, -17.1868, - -15.4316, -13.7494, -12.0673, -10.3852, -8.70311, -7.021, - -5.33888, -3.65677, -1.97466, -0.292541, 1.38957, 3.07169, - 4.7538, 6.43591, 8.19116, 9.87327, 11.5554, 13.2375, 14.9196, - 16.6017, 18.2838, 19.966, 21.6481, 23.3302, 25.0123, 26.6944, - 28.3765, 30.0586, 31.8139, -47.986, -45.8388, -43.6916, - -41.5443, -37.1565, -35.0093, -32.862, -30.7148, -28.5676, - -26.4203, -24.2731, -22.1259, -19.9786, -17.8314, -15.6842, - -13.5369, -11.3897, -9.24245, -7.00185, -4.85462, -2.70738, - -0.560148, 1.58709, 3.73432, 5.88156, 8.02879, 10.176, - 12.3233, 14.4705, 16.6177, 18.765, 20.9122, 23.1528, 25.3, - 27.4473, 29.5945, 31.7417, 33.889, 36.0362, 38.1834, 40.3307, - 42.4779, 44.6251, 46.7724, 48.9196, 24.1427, 22.5741, 20.9373, - 19.437, 17.8684, 16.2998, 14.7312, 13.1626, 11.594, 10.0254, - 8.45678, 6.81998, 5.31959, 3.68279, 2.04599, 0.477399, -1.023, - -2.59159, -4.16019, -5.66059, -7.29738, -8.86598, -10.4346, - -12.0032, -13.5718, -15.1404, -21.5511, -23.0515, -45.2569, - -40.6613, -38.3635, -36.0656, -33.7678, -31.47, -29.1722, - -26.8744, -24.5766, -22.2788, -19.981, -17.5832, -15.2854, - -12.9876, -10.6898, -8.392, -6.09419, -3.79638, -1.49857, - 0.799239, 3.09705, 5.39486, 7.69267, 9.99048, 12.2883, - 14.5861, 16.9838, 19.2816, 21.5794, 23.8773, 26.1751, 28.4729, - 30.7707, 33.0685, 35.3663, 37.6641, 39.9619, 42.2597, 49.8478, - 47.7895, 45.7311, 43.5833, 41.525, 39.4666, 37.4083, 33.2021, - 31.2332, 29.1749, 27.027, 24.9687, 22.9103, 20.852, 18.7936, - 16.7353, 14.6769, 12.6186, 10.5602, 8.50188, 6.44353, 4.38518, - 2.32683, 0.26848, -1.87936, -3.93771, -5.99606, -8.05441, - -10.1128, -12.1711, -14.14, -16.1983, -18.3462, -20.4045, - -22.4629, -24.5212, -26.5796, -28.6379, -30.6962, -32.7546, - -34.8129, -36.8713, -38.9296, -40.988, -43.0463, -45.1047, - 6.53648, 5.99538, 5.5625, 5.0214, 4.52359, 4.02578, 3.52797, - 3.0518, 2.53234, 2.03453, 1.53672, 1.01727, 0.541099, - 0.021644, -0.454523, -0.93069, -1.47179, -1.9696, -2.46741, - -2.94358, -3.44139, -3.96084, -4.45866, -4.95647, -5.45428, - -5.95209, -6.4499, -49.1077, -42.3712, -40.2234, -37.8803, - -35.7324, -33.487, -31.2415, -28.996, -26.6529, -24.4074, - -22.1619, -19.9164, -17.671, -15.5231, -13.18, -10.9345, - -8.68903, -6.44355, -4.19807, -1.95259, 0.292889, 2.53837, - 4.78385, 7.02933, 9.37244, 11.5203, 13.8634, 16.1089, 18.3544, - 20.5998, 22.8453, 25.0908, 27.3363, 29.5818, 31.8272, 34.1704, - 36.4158, 38.5637, 41.0044, 43.2499, 45.3001, 47.5456, - 49.8887}; -double ew[59], ns[99], grid[11682], fit_eval[200], ci_eval[30]; -double tmp, range = 98, coverage = .99; -long n = 323, p = 2, m, se_fit = FALSE; -int i, j, k; - -main() { - printf("\nloess(&galaxy):\n"); - loess_setup(direction, velocity, n, p, &galaxy); - galaxy.model.span = 0.35; - galaxy.model.normalize = FALSE; - galaxy.model.family = "symmetric"; - loess(&galaxy); - loess_summary(&galaxy); - - m = 5841; - tmp = -29.0; - for(i = 0; i < 59; i++) - ew[i] = tmp++; - tmp = -49.0; - for(i = 0; i < 99; i++) - ns[i] = tmp++; - for(i = 0; i < 99; i++) { - k = i * 59; - for(j = 0; j < 59; j++) { - grid[k + j] = ew[j]; - grid[m + k + j] = ns[i]; - } - } - predict(grid, m, &galaxy, &galaxy_contour, se_fit); - - m = 100; - tmp = range / 99; - for(i = 0; i < 100; i++) { - fit_eval[i + 100] = -49 + tmp * i; - fit_eval[i] = fit_eval[i + 100] / (-3.7); - } - predict(fit_eval, m, &galaxy, &spine_fit, se_fit); - - m = 15; - se_fit = TRUE; - tmp = range / 14; - for(i = 0; i < m; i++) { - ci_eval[i + m] = -49 + tmp * i; - ci_eval[i] = fit_eval[i + 100] / (-3.7); - } - predict(ci_eval, m, &galaxy, &spine_se, se_fit); - pointwise(&spine_se, m, coverage, &spine_ci); - - loess_free_mem(&galaxy); - pred_free_mem(&galaxy_contour); - pred_free_mem(&spine_fit); - pred_free_mem(&spine_se); -} //GO.SYSIN DD galaxy.c echo gas.c 1>&2 sed >gas.c <<'//GO.SYSIN DD gas.c' 's/^-//' -/* sample program for the gas data using loess */ - -#include -#include "loess.h" - -struct loess_struct gas, gas_null; -struct pred_struct gas_pred; -struct ci_struct gas_ci; -struct anova_struct gas_anova; -double NOx[] = {4.818, 2.849, 3.275, 4.691, 4.255, 5.064, 2.118, 4.602, - 2.286, 0.97, 3.965, 5.344, 3.834, 1.99, 5.199, 5.283, - 3.752, 0.537, 1.64, 5.055, 4.937, 1.561}; -double E[] = {0.831, 1.045, 1.021, 0.97, 0.825, 0.891, 0.71, 0.801, - 1.074, 1.148, 1, 0.928, 0.767, 0.701, 0.807, 0.902, - 0.997, 1.224, 1.089, 0.973, 0.98, 0.665}; -double gas_fit_E[] = {0.665, 0.949, 1.224}; -double newdata[] = {0.6650000, 0.7581667, 0.8513333, 0.9445000, - 1.0376667, 1.1308333, 1.2240000}; -double coverage = .99; -long i, n = 22, p = 1, m = 3, se_fit = FALSE; - -main() { - printf("\nloess(&gas):\n"); - loess_setup(E, NOx, n, p, &gas); - gas.model.span = 2.0 / 3.0; - loess(&gas); - loess_summary(&gas); - - printf("\nloess(&gas_null):\n"); - loess_setup(E, NOx, n, p, &gas_null); - gas_null.model.span = 1.0; - loess(&gas_null); - loess_summary(&gas_null); - - printf("\npredict(gas_fit_E, m, &gas, &gas_pred, %d):\n", se_fit); - predict(gas_fit_E, m, &gas, &gas_pred, se_fit); - for(i = 0; i < m; i++) - printf("%g ", gas_pred.fit[i]); - printf("\n"); - - m = 7; - se_fit = TRUE; - predict(newdata, m, &gas, &gas_pred, se_fit); - printf("\npointwise(&gas_pred, m, coverage, &gas_ci):\n"); - pointwise(&gas_pred, m, coverage, &gas_ci); - for(i = 0; i < m; i++) - printf("%g ", gas_ci.upper[i]); - printf("\n"); - for(i = 0; i < m; i++) - printf("%g ", gas_ci.fit[i]); - printf("\n"); - for(i = 0; i < m; i++) - printf("%g ", gas_ci.lower[i]); - printf("\n"); - - printf("\nanova(&gas_null, &gas, &gas_anova):\n"); - anova(&gas_null, &gas, &gas_anova); - printf("%g %g %g %g\n", gas_anova.dfn, gas_anova.dfd, - gas_anova.F_value, gas_anova.Pr_F); - - loess_free_mem(&gas); - loess_free_mem(&gas_null); - pred_free_mem(&gas_pred); - pw_free_mem(&gas_ci); -} - - - - //GO.SYSIN DD gas.c echo loess.c 1>&2 sed >loess.c <<'//GO.SYSIN DD loess.c' 's/^-//' -#include "S.h" -#include "loess.h" - -static char *surf_stat; - -void -loess_setup(x, y, n, p, lo) -double *x, *y; -long n, p; -struct loess_struct *lo; -{ - int i, max_kd; - - max_kd = n > 200 ? n : 200; - - lo->in.y = (double *) malloc(n * sizeof(double)); - lo->in.x = (double *) malloc(n * p * sizeof(double)); - lo->in.weights = (double *) malloc(n * sizeof(double)); - for(i = 0; i < (n * p); i++) - lo->in.x[i] = x[i]; - for(i = 0; i < n; i++) { - lo->in.y[i] = y[i]; - lo->in.weights[i] = 1; - } - lo->in.n = n; - lo->in.p = p; - lo->model.span = 0.75; - lo->model.degree = 2; - lo->model.normalize = TRUE; - for(i = 0; i < 8; i++) - lo->model.parametric[i] = lo->model.drop_square[i] = FALSE; - lo->model.family = "gaussian"; - lo->control.surface = "interpolate"; - lo->control.statistics = "approximate"; - lo->control.cell = 0.2; - lo->control.trace_hat = "wait.to.decide"; - lo->control.iterations = 4; - - lo->out.fitted_values = (double *) malloc(n * sizeof(double)); - lo->out.fitted_residuals = (double *) malloc(n * sizeof(double)); - lo->out.pseudovalues = (double *) malloc(n * sizeof(double)); - lo->out.diagonal = (double *) malloc(n * sizeof(double)); - lo->out.robust = (double *) malloc(n * sizeof(double)); - lo->out.divisor = (double *) malloc(p * sizeof(double)); - - lo->kd_tree.parameter = (long *) malloc(7 * sizeof(long)); - lo->kd_tree.a = (long *) malloc(max_kd * sizeof(long)); - lo->kd_tree.xi = (double *) malloc(max_kd * sizeof(double)); - lo->kd_tree.vert = (double *) malloc(p * 2 * sizeof(double)); - lo->kd_tree.vval = (double *) malloc((p + 1) * max_kd * sizeof(double)); -} - -void -loess(lo) -struct loess_struct *lo; -{ - long size_info[2], iterations; - void loess_(); - - size_info[0] = lo->in.p; - size_info[1] = lo->in.n; - - iterations = (!strcmp(lo->model.family, "gaussian")) ? 0 : - lo->control.iterations; - if(!strcmp(lo->control.trace_hat, "wait.to.decide")) { - if(!strcmp(lo->control.surface, "interpolate")) - lo->control.trace_hat = (lo->in.n < 500) ? "exact" : "approximate"; - else - lo->control.trace_hat = "exact"; - } - loess_(lo->in.y, lo->in.x, size_info, lo->in.weights, - &lo->model.span, - &lo->model.degree, - lo->model.parametric, - lo->model.drop_square, - &lo->model.normalize, - &lo->control.statistics, - &lo->control.surface, - &lo->control.cell, - &lo->control.trace_hat, - &iterations, - lo->out.fitted_values, - lo->out.fitted_residuals, - &lo->out.enp, - &lo->out.s, - &lo->out.one_delta, - &lo->out.two_delta, - lo->out.pseudovalues, - &lo->out.trace_hat, - lo->out.diagonal, - lo->out.robust, - lo->out.divisor, - lo->kd_tree.parameter, - lo->kd_tree.a, - lo->kd_tree.xi, - lo->kd_tree.vert, - lo->kd_tree.vval); -} - -void -loess_(y, x_, size_info, weights, span, degree, parametric, drop_square, - normalize, statistics, surface, cell, trace_hat_in, iterations, - fitted_values, fitted_residuals, enp, s, one_delta, two_delta, - pseudovalues, trace_hat_out, diagonal, robust, divisor, - parameter, a, xi, vert, vval) -double *y, *x_, *weights, *span, *cell, *pseudovalues, - *fitted_values, *fitted_residuals, *enp, *s, *one_delta, *two_delta, - *trace_hat_out, *diagonal, *robust, *divisor, *xi, *vert, *vval; -long *size_info, *degree, *parametric, *drop_square, *normalize, - *iterations, *parameter, *a; -char **statistics, **surface, **trace_hat_in; -{ - double *x, *x_tmp, new_cell, trL, delta1, delta2, sum_squares = 0, - *pseudo_resid, *temp, *xi_tmp, *vert_tmp, *vval_tmp, - *diag_tmp, trL_tmp = 0, d1_tmp = 0, d2_tmp = 0, sum, mean; - long i, j, k, p, N, D, sum_drop_sqr = 0, sum_parametric = 0, - setLf, nonparametric = 0, *order_parametric, - *order_drop_sqr, zero = 0, max_kd, *a_tmp, *param_tmp; - int cut, comp(); - char *new_stat; - void condition(); - - D = size_info[0]; - N = size_info[1]; - max_kd = (N > 200 ? N : 200); - *one_delta = *two_delta = *trace_hat_out = 0; - - x = (double *) malloc(D * N * sizeof(double)); - x_tmp = (double *) malloc(D * N * sizeof(double)); - temp = (double *) malloc(N * sizeof(double)); - a_tmp = (long *) malloc(max_kd * sizeof(long)); - xi_tmp = (double *) malloc(max_kd * sizeof(double)); - vert_tmp = (double *) malloc(D * 2 * sizeof(double)); - vval_tmp = (double *) malloc((D + 1) * max_kd * sizeof(double)); - diag_tmp = (double *) malloc(N * sizeof(double)); - param_tmp = (long *) malloc(N * sizeof(long)); - order_parametric = (long *) malloc(D * sizeof(long)); - order_drop_sqr = (long *) malloc(D * sizeof(long)); - if((*iterations) > 0) - pseudo_resid = (double *) malloc(N * sizeof(double)); - - new_cell = (*span) * (*cell); - for(i = 0; i < N; i++) - robust[i] = 1; - for(i = 0; i < (N * D); i++) - x_tmp[i] = x_[i]; - if((*normalize) && (D > 1)) { - cut = ceil(0.100000000000000000001 * N); - for(i = 0; i < D; i++) { - k = i * N; - for(j = 0; j < N; j++) - temp[j] = x_[k + j]; - qsort(temp, N, sizeof(double), comp); - sum = 0; - for(j = cut; j <= (N - cut - 1); j++) - sum = sum + temp[j]; - mean = sum / (N - 2 * cut); - sum = 0; - for(j = cut; j <= (N - cut - 1); j++) { - temp[j] = temp[j] - mean; - sum = sum + temp[j] * temp[j]; - } - divisor[i] = sqrt(sum / (N - 2 * cut - 1)); - for(j = 0; j < N; j++) { - p = k + j; - x_tmp[p] = x_[p] / divisor[i]; - } - } - } - else - for(i = 0; i < D; i++) divisor[i] = 1; - j = D - 1; - for(i = 0; i < D; i++) { - sum_drop_sqr = sum_drop_sqr + drop_square[i]; - sum_parametric = sum_parametric + parametric[i]; - if(parametric[i]) - order_parametric[j--] = i; - else - order_parametric[nonparametric++] = i; - } - for(i = 0; i < D; i++) { - order_drop_sqr[i] = 2 - drop_square[order_parametric[i]]; - k = i * N; - p = order_parametric[i] * N; - for(j = 0; j < N; j++) - x[k + j] = x_tmp[p + j]; - } - if((*degree) == 1 && sum_drop_sqr) { - fprintf(stderr, "Specified the square of a factor predictor to be dropped when degree = 1"); - exit(1); - } - if(D == 1 && sum_drop_sqr) { - fprintf(stderr, "Specified the square of a predictor to be dropped with only one numeric predictor"); - exit(1); - } - if(sum_parametric == D) { - fprintf(stderr, "Specified parametric for all predictors"); - exit(1); - } - for(j = 0; j <= (*iterations); j++) { - new_stat = j ? "none" : *statistics; - for(i = 0; i < N; i++) - robust[i] = weights[i] * robust[i]; - condition(surface, new_stat, trace_hat_in); - setLf = !strcmp(surf_stat, "interpolate/exact"); - loess_raw(y, x, weights, robust, &D, &N, span, degree, - &nonparametric, order_drop_sqr, &sum_drop_sqr, - &new_cell, &surf_stat, fitted_values, parameter, a, - xi, vert, vval, diagonal, &trL, &delta1, &delta2, - &setLf); - if(j == 0) { - *trace_hat_out = trL; - *one_delta = delta1; - *two_delta = delta2; - } - for(i = 0; i < N; i++) - fitted_residuals[i] = y[i] - fitted_values[i]; - if(j < (*iterations)) - F77_SUB(lowesw)(fitted_residuals, &N, robust, temp); - } - if((*iterations) > 0) { - F77_SUB(lowesp)(&N, y, fitted_values, weights, robust, temp, pseudovalues); - - loess_raw(pseudovalues, x, weights, weights, &D, &N, span, - degree, &nonparametric, order_drop_sqr, &sum_drop_sqr, - &new_cell, &surf_stat, temp, param_tmp, a_tmp, xi_tmp, - vert_tmp, vval_tmp, diag_tmp, &trL_tmp, &d1_tmp, &d2_tmp, &zero); - for(i = 0; i < N; i++) - pseudo_resid[i] = pseudovalues[i] - temp[i]; - } - if((*iterations) == 0) - for(i = 0; i < N; i++) - sum_squares = sum_squares + weights[i] * - fitted_residuals[i] * fitted_residuals[i]; - else - for(i = 0; i < N; i++) - sum_squares = sum_squares + weights[i] * - pseudo_resid[i] * pseudo_resid[i]; - *enp = (*one_delta) + 2 * (*trace_hat_out) - N; - *s = sqrt(sum_squares / (*one_delta)); - - free(x); - free(x_tmp); - free(temp); - free(xi_tmp); - free(vert_tmp); - free(vval_tmp); - free(diag_tmp); - free(a_tmp); - free(param_tmp); - free(order_parametric); - free(order_drop_sqr); - if((*iterations) > 0) - free(pseudo_resid); -} - -void -loess_free_mem(lo) -struct loess_struct *lo; -{ - free(lo->in.x); - free(lo->in.y); - free(lo->in.weights); - free(lo->out.fitted_values); - free(lo->out.fitted_residuals); - free(lo->out.pseudovalues); - free(lo->out.diagonal); - free(lo->out.robust); - free(lo->out.divisor); - free(lo->kd_tree.parameter); - free(lo->kd_tree.a); - free(lo->kd_tree.xi); - free(lo->kd_tree.vert); - free(lo->kd_tree.vval); -} - -void -loess_summary(lo) -struct loess_struct *lo; -{ - printf("Number of Observations: %d\n", lo->in.n); - printf("Equivalent Number of Parameters: %.1f\n", lo->out.enp); - if(!strcmp(lo->model.family, "gaussian")) - printf("Residual Standard Error: "); - else - printf("Residual Scale Estimate: "); - printf("%.4f\n", lo->out.s); -} - -void -condition(surface, new_stat, trace_hat_in) -char **surface, *new_stat, **trace_hat_in; -{ - if(!strcmp(*surface, "interpolate")) { - if(!strcmp(new_stat, "none")) - surf_stat = "interpolate/none"; - else if(!strcmp(new_stat, "exact")) - surf_stat = "interpolate/exact"; - else if(!strcmp(new_stat, "approximate")) - { - if(!strcmp(*trace_hat_in, "approximate")) - surf_stat = "interpolate/2.approx"; - else if(!strcmp(*trace_hat_in, "exact")) - surf_stat = "interpolate/1.approx"; - } - } - else if(!strcmp(*surface, "direct")) { - if(!strcmp(new_stat, "none")) - surf_stat = "direct/none"; - else if(!strcmp(new_stat, "exact")) - surf_stat = "direct/exact"; - else if(!strcmp(new_stat, "approximate")) - surf_stat = "direct/approximate"; - } -} - -int -comp(d1, d2) -double *d1, *d2; -{ - if(*d1 < *d2) - return(-1); - else if(*d1 == *d2) - return(0); - else - return(1); -} //GO.SYSIN DD loess.c echo loess.h 1>&2 sed >loess.h <<'//GO.SYSIN DD loess.h' 's/^-//' -/* for the meaning of these fields, see struct.m */ -/* longs are used here so that the codes can be called from S */ - -#define TRUE 1 -#define FALSE 0 - -extern struct loess_struct { - struct { - long n; - long p; - double *y; - double *x; - double *weights; - } in; - struct { - double span; - long degree; - long normalize; - long parametric[8]; - long drop_square[8]; - char *family; - } model; - struct { - char *surface; - char *statistics; - double cell; - char *trace_hat; - long iterations; - } control; - struct { - long *parameter; - long *a; - double *xi; - double *vert; - double *vval; - } kd_tree; - struct { - double *fitted_values; - double *fitted_residuals; - double enp; - double s; - double one_delta; - double two_delta; - double *pseudovalues; - double trace_hat; - double *diagonal; - double *robust; - double *divisor; - } out; -}; - -extern struct pred_struct { - double *fit; - double *se_fit; - double residual_scale; - double df; -}; - -extern struct anova_struct { - double dfn; - double dfd; - double F_value; - double Pr_F; -}; - -extern struct ci_struct { - double *fit; - double *upper; - double *lower; -}; - //GO.SYSIN DD loess.h echo loess.m 1>&2 sed >loess.m <<'//GO.SYSIN DD loess.m' 's/^-//' -NAME - - loess_setup, loess, loess_summary, loess_free_mem, anova - -SYNOPSIS - - #include "loess.h" - double *x, *y; - long n, p; - struct loess_struct *lo, *lo2; - struct anova_struct *aov; - - void loess_setup(x, y, n, p, lo) - - void loess(lo) - - void loess_summary(lo) - - void loess_free_mem(lo) - - void anova(lo, lo2, aov); - -PARAMETERS - - x predictors vector (of length n * p) - The j-th coordinate of the i-th point is in x[i+n*j], - where 0<=j&2 sed >loessc.c <<'//GO.SYSIN DD loessc.c' 's/^-//' -#include "S.h" - -#define min(x,y) ((x) < (y) ? (x) : (y)) -#define max(x,y) ((x) > (y) ? (x) : (y)) -#define GAUSSIAN 1 -#define SYMMETRIC 0 - -static long *iv, liv, lv, tau; -static double *v; - -loess_raw(y, x, weights, robust, d, n, span, degree, nonparametric, - drop_square, sum_drop_sqr, cell, surf_stat, surface, parameter, a, - xi, vert, vval, diagonal, trL, one_delta, two_delta, setLf) -double *y, *x, *weights, *robust, *span, *cell, *surface, *xi, *vert, - *vval, *diagonal, *trL, *one_delta, *two_delta; -long *d, *n, *parameter, *a, *degree, *nonparametric, *drop_square, - *sum_drop_sqr, *setLf; -char **surf_stat; -{ - long zero = 0, one = 1, two = 2, nsing, i, k; - double *hat_matrix, *LL; - - *trL = 0; - loess_workspace(d, n, span, degree, nonparametric, drop_square, - sum_drop_sqr, setLf); - v[1] = *cell; - if(!strcmp(*surf_stat, "interpolate/none")) { - F77_SUB(lowesb)(x, y, robust, &zero, &zero, iv, &liv, &lv, v); - F77_SUB(lowese)(iv, &liv, &lv, v, n, x, surface); - loess_prune(parameter, a, xi, vert, vval); - } - else if (!strcmp(*surf_stat, "direct/none")) { - F77_SUB(lowesf)(x, y, robust, iv, &liv, &lv, v, n, x, - &zero, &zero, surface); - } - else if (!strcmp(*surf_stat, "interpolate/1.approx")) { - F77_SUB(lowesb)(x, y, weights, diagonal, &one, iv, &liv, &lv, v); - F77_SUB(lowese)(iv, &liv, &lv, v, n, x, surface); - nsing = iv[29]; - for(i = 0; i < (*n); i++) *trL = *trL + diagonal[i]; - F77_SUB(lowesa)(trL, n, d, &tau, &nsing, one_delta, two_delta); - loess_prune(parameter, a, xi, vert, vval); - } - else if (!strcmp(*surf_stat, "interpolate/2.approx")) { - F77_SUB(lowesb)(x, y, robust, &zero, &zero, iv, &liv, &lv, v); - F77_SUB(lowese)(iv, &liv, &lv, v, n, x, surface); - nsing = iv[29]; - F77_SUB(ehg196)(&tau, d, span, trL); - F77_SUB(lowesa)(trL, n, d, &tau, &nsing, one_delta, two_delta); - loess_prune(parameter, a, xi, vert, vval); - } - else if (!strcmp(*surf_stat, "direct/approximate")) { - F77_SUB(lowesf)(x, y, weights, iv, &liv, &lv, v, n, x, - diagonal, &one, surface); - nsing = iv[29]; - for(i = 0; i < (*n); i++) *trL = *trL + diagonal[i]; - F77_SUB(lowesa)(trL, n, d, &tau, &nsing, one_delta, two_delta); - } - else if (!strcmp(*surf_stat, "interpolate/exact")) { - hat_matrix = Calloc((*n)*(*n), double); - LL = Calloc((*n)*(*n), double); - F77_SUB(lowesb)(x, y, weights, diagonal, &one, iv, &liv, &lv, v); - F77_SUB(lowesl)(iv, &liv, &lv, v, n, x, hat_matrix); - F77_SUB(lowesc)(n, hat_matrix, LL, trL, one_delta, two_delta); - F77_SUB(lowese)(iv, &liv, &lv, v, n, x, surface); - loess_prune(parameter, a, xi, vert, vval); - Free(hat_matrix); - Free(LL); - } - else if (!strcmp(*surf_stat, "direct/exact")) { - hat_matrix = Calloc((*n)*(*n), double); - LL = Calloc((*n)*(*n), double); - F77_SUB(lowesf)(x, y, weights, iv, liv, lv, v, n, x, - hat_matrix, &two, surface); - F77_SUB(lowesc)(n, hat_matrix, LL, trL, one_delta, two_delta); - k = (*n) + 1; - for(i = 0; i < (*n); i++) - diagonal[i] = hat_matrix[i * k]; - Free(hat_matrix); - Free(LL); - } - loess_free(); -} - -loess_dfit(y, x, x_evaluate, weights, span, degree, nonparametric, - drop_square, sum_drop_sqr, d, n, m, fit) -double *y, *x, *x_evaluate, *weights, *span, *fit; -long *degree, *nonparametric, *drop_square, *sum_drop_sqr, *d, *n, *m; -{ - long zero = 0, one = 1; - - loess_workspace(d, n, span, degree, nonparametric, drop_square, - sum_drop_sqr, &zero); - F77_SUB(lowesf)(x, y, weights, iv, &liv, &lv, v, m, x_evaluate, - &zero, &zero, fit); - loess_free(); -} - -loess_dfitse(y, x, x_evaluate, weights, robust, family, span, degree, - nonparametric, drop_square, sum_drop_sqr, d, n, m, fit, L) -double *y, *x, *x_evaluate, *weights, *robust, *span, *fit, *L; -long *family, *degree, *nonparametric, *drop_square, *sum_drop_sqr, - *d, *n, *m; -{ - long zero = 0, one = 1, two = 2; - - loess_workspace(d, n, span, degree, nonparametric, drop_square, - sum_drop_sqr, &zero); - if(*family == GAUSSIAN) - F77_SUB(lowesf)(x, y, weights, iv, &liv, &lv, v, m, - x_evaluate, L, &two, fit); - else if(*family == SYMMETRIC) - { - F77_SUB(lowesf)(x, y, weights, iv, &liv, &lv, v, m, - x_evaluate, L, &two, fit); - F77_SUB(lowesf)(x, y, robust, iv, &liv, &lv, v, m, - x_evaluate, &zero, &zero, fit); - } - loess_free(); -} -loess_ifit(parameter, a, xi, vert, vval, m, x_evaluate, fit) -double *xi, *vert, *vval, *x_evaluate, *fit; -long *parameter, *a, *m; -{ - loess_grow(parameter, a, xi, vert, vval); - F77_SUB(lowese)(iv, &liv, &lv, v, m, x_evaluate, fit); - loess_free(); -} - -loess_ise(y, x, x_evaluate, weights, span, degree, nonparametric, - drop_square, sum_drop_sqr, cell, d, n, m, fit, L) -double *y, *x, *x_evaluate, *weights, *span, *cell, *fit, *L; -long *degree, *nonparametric, *drop_square, *sum_drop_sqr, *d, *n, *m; -{ - long zero = 0, one = 1; - - loess_workspace(d, n, span, degree, nonparametric, drop_square, - sum_drop_sqr, &one); - v[1] = *cell; - F77_SUB(lowesb)(x, y, weights, &zero, &zero, iv, &liv, &lv, v); - F77_SUB(lowesl)(iv, &liv, &lv, v, m, x_evaluate, L); - loess_free(); -} - -loess_workspace(d, n, span, degree, nonparametric, drop_square, - sum_drop_sqr, setLf) -long *d, *n, *degree, *nonparametric, *drop_square, *sum_drop_sqr, - *setLf; -double *span; -{ - long D, N, tau0, nvmax, nf, version = 106, i; - - D = *d; - N = *n; - nvmax = max(200, N); - nf = min(N, floor(N * (*span))); - tau0 = ((*degree) > 1) ? ((D + 2) * (D + 1) * 0.5) : (D + 1); - tau = tau0 - (*sum_drop_sqr); - lv = 50 + (3 * D + 3) * nvmax + N + (tau0 + 2) * nf; - liv = 50 + ((long)pow((double)2, (double)D) + 4) * nvmax + 2 * N; - if(*setLf) { - lv = lv + (D + 1) * nf * nvmax; - liv = liv + nf * nvmax; - } - iv = Calloc(liv, long); - v = Calloc(lv, double); - - F77_SUB(lowesd)(&version, iv, &liv, &lv, v, d, n, span, degree, - &nvmax, setLf); - iv[32] = *nonparametric; - for(i = 0; i < D; i++) - iv[i + 40] = drop_square[i]; -} - -loess_prune(parameter, a, xi, vert, vval) -double *xi, *vert, *vval; -long *parameter, *a; -{ - long d, vc, a1, v1, xi1, vv1, nc, nv, nvmax, i, j, k; - - d = iv[1]; - vc = iv[3] - 1; - nc = iv[4]; - nv = iv[5]; - a1 = iv[6] - 1; - v1 = iv[10] - 1; - xi1 = iv[11] - 1; - vv1 = iv[12] - 1; - nvmax = iv[13]; - - for(i = 0; i < 5; i++) - parameter[i] = iv[i + 1]; - parameter[5] = iv[21] - 1; - parameter[6] = iv[14] - 1; - - for(i = 0; i < d; i++){ - k = nvmax * i; - vert[i] = v[v1 + k]; - vert[i + d] = v[v1 + vc + k]; - } - for(i = 0; i < nc; i++) { - xi[i] = v[xi1 + i]; - a[i] = iv[a1 + i]; - } - k = (d + 1) * nv; - for(i = 0; i < k; i++) - vval[i] = v[vv1 + i]; -} - -loess_grow(parameter, a, xi, vert, vval) -double *xi, *vert, *vval; -long *parameter, *a; -{ - long d, vc, nc, nv, a1, v1, xi1, vv1, i, j, k; - - d = parameter[0]; - vc = parameter[2]; - nc = parameter[3]; - nv = parameter[4]; - liv = parameter[5]; - lv = parameter[6]; - iv = Calloc(liv, long); - v = Calloc(lv, double); - - iv[1] = d; - iv[2] = parameter[1]; - iv[3] = vc; - iv[5] = iv[13] = nv; - iv[4] = iv[16] = nc; - iv[6] = 50; - iv[7] = iv[6] + nc; - iv[8] = iv[7] + vc * nc; - iv[9] = iv[8] + nc; - iv[10] = 50; - iv[12] = iv[10] + nv * d; - iv[11] = iv[12] + (d + 1) * nv; - iv[27] = 173; - - v1 = iv[10] - 1; - xi1 = iv[11] - 1; - a1 = iv[6] - 1; - vv1 = iv[12] - 1; - - for(i = 0; i < d; i++) { - k = nv * i; - v[v1 + k] = vert[i]; - v[v1 + vc - 1 + k] = vert[i + d]; - } - for(i = 0; i < nc; i++) { - v[xi1 + i] = xi[i]; - iv[a1 + i] = a[i]; - } - k = (d + 1) * nv; - for(i = 0; i < k; i++) - v[vv1 + i] = vval[i]; - - F77_SUB(ehg169)(&d, &vc, &nc, &nc, &nv, &nv, v+v1, iv+a1, - v+xi1, iv+iv[7]-1, iv+iv[8]-1, iv+iv[9]-1); -} - -loess_free() -{ - Free(v); - Free(iv); -} - -/* begin ehg's FORTRAN-callable C-codes */ - -void -F77_SUB(ehg182)(i) - int *i; -{ - char *mess, mess2[50]; - switch(*i){ -case 100: mess="wrong version number in lowesd. Probably typo in caller."; break; -case 101: mess="d>dMAX in ehg131. Need to recompile with increased dimensions."; break; -case 102: mess="liv too small. (Discovered by lowesd)"; break; -case 103: mess="lv too small. (Discovered by lowesd)"; break; -case 104: mess="span too small. fewer data values than degrees of freedom."; break; -case 105: mess="k>d2MAX in ehg136. Need to recompile with increased dimensions."; break; -case 106: mess="lwork too small"; break; -case 107: mess="invalid value for kernel"; break; -case 108: mess="invalid value for ideg"; break; -case 109: mess="lowstt only applies when kernel=1."; break; -case 110: mess="not enough extra workspace for robustness calculation"; break; -case 120: mess="zero-width neighborhood. make span bigger"; break; -case 121: mess="all data on boundary of neighborhood. make span bigger"; break; -case 122: mess="extrapolation not allowed with blending"; break; -case 123: mess="ihat=1 (diag L) in l2fit only makes sense if z=x (eval=data)."; break; -case 171: mess="lowesd must be called first."; break; -case 172: mess="lowesf must not come between lowesb and lowese, lowesr, or lowesl."; break; -case 173: mess="lowesb must come before lowese, lowesr, or lowesl."; break; -case 174: mess="lowesb need not be called twice."; break; -case 175: mess="need setLf=.true. for lowesl."; break; -case 180: mess="nv>nvmax in cpvert."; break; -case 181: mess="nt>20 in eval."; break; -case 182: mess="svddc failed in l2fit."; break; -case 183: mess="didnt find edge in vleaf."; break; -case 184: mess="zero-width cell found in vleaf."; break; -case 185: mess="trouble descending to leaf in vleaf."; break; -case 186: mess="insufficient workspace for lowesf."; break; -case 187: mess="insufficient stack space"; break; -case 188: mess="lv too small for computing explicit L"; break; -case 191: mess="computed trace L was negative; something is wrong!"; break; -case 192: mess="computed delta was negative; something is wrong!"; break; -case 193: mess="workspace in loread appears to be corrupted"; break; -case 194: mess="trouble in l2fit/l2tr"; break; -case 195: mess="only constant, linear, or quadratic local models allowed"; break; -case 196: mess="degree must be at least 1 for vertex influence matrix"; break; -case 999: mess="not yet implemented"; break; -default: sprintf(mess=mess2,"Assert failed; error code %d\n",*i); break; - } - Recover(mess,NULL_ENTRY); /* in /usr/s/current/src/qpe/debug.c */ -} - -void -F77_SUB(ehg183)(s,i,n,inc) - char *s; - int *i, *n, *inc; -{ - char mess[4000], num[20]; - int j; - strcpy(mess,s); - for (j=0; j<*n; j++) { - sprintf(num," %d",i[j * *inc]); - strcat(mess,num); - } - strcat(mess,"\n"); - Warning(mess,NULL_ENTRY); -} - -void -F77_SUB(ehg184)(s,x,n,inc) - char *s; - double *x; - int *n, *inc; -{ - char mess[4000], num[30]; - int j; - strcpy(mess,s); - for (j=0; j<*n; j++) { - sprintf(num," %.5g",x[j * *inc]); - strcat(mess,num); - } - strcat(mess,"\n"); - Warning(mess,NULL_ENTRY); -} //GO.SYSIN DD loessc.c echo loessf.f 1>&2 sed >loessf.f <<'//GO.SYSIN DD loessf.f' 's/^-//' - subroutine ehg126(d,n,vc,x,v,nvmax) - integer d,execnt,i,j,k,n,nv,nvmax,vc - DOUBLE PRECISION machin,alpha,beta,mu,t - DOUBLE PRECISION v(nvmax,d),x(n,d) - DOUBLE PRECISION D1MACH - external D1MACH - save machin,execnt - data execnt /0/ -c MachInf -> machin - execnt=execnt+1 - if(execnt.eq.1)then - machin=D1MACH(2) - end if -c fill in vertices for bounding box of $x$ -c lower left, upper right - do 3 k=1,d - alpha=machin - beta=-machin - do 4 i=1,n - t=x(i,k) - alpha=min(alpha,t) - beta=max(beta,t) - 4 continue -c expand the box a little - mu=0.005D0*max(beta-alpha,1.d-10*max(DABS(alpha),DABS(beta))+1. - +d-30) - alpha=alpha-mu - beta=beta+mu - v(1,k)=alpha - v(vc,k)=beta - 3 continue -c remaining vertices - do 5 i=2,vc-1 - j=i-1 - do 6 k=1,d - v(i,k)=v(1+mod(j,2)*(vc-1),k) - j=DFLOAT(j)/2.D0 - 6 continue - 5 continue - return - end - subroutine ehg125(p,nv,v,vhit,nvmax,d,k,t,r,s,f,l,u) - logical i1,i2,match - integer d,execnt,h,i,i3,j,k,m,mm,nv,nvmax,p,r,s - integer f(r,0:1,s),l(r,0:1,s),u(r,0:1,s),vhit(nvmax) - DOUBLE PRECISION t - DOUBLE PRECISION v(nvmax,d) - external ehg182 - save execnt - data execnt /0/ - execnt=execnt+1 - h=nv - do 3 i=1,r - do 4 j=1,s - h=h+1 - do 5 i3=1,d - v(h,i3)=v(f(i,0,j),i3) - 5 continue - v(h,k)=t -c check for redundant vertex - match=.false. - m=1 -c top of while loop - 6 if(.not.match)then - i1=(m.le.nv) - else - i1=.false. - end if - if(.not.(i1))goto 7 - match=(v(m,1).eq.v(h,1)) - mm=2 -c top of while loop - 8 if(match)then - i2=(mm.le.d) - else - i2=.false. - end if - if(.not.(i2))goto 9 - match=(v(m,mm).eq.v(h,mm)) - mm=mm+1 - goto 8 -c bottom of while loop - 9 m=m+1 - goto 6 -c bottom of while loop - 7 m=m-1 - if(match)then - h=h-1 - else - m=h - if(vhit(1).ge.0)then - vhit(m)=p - end if - end if - l(i,0,j)=f(i,0,j) - l(i,1,j)=m - u(i,0,j)=m - u(i,1,j)=f(i,1,j) - 4 continue - 3 continue - nv=h - if(.not.(nv.le.nvmax))then - call ehg182(180) - end if - return - end - integer function ehg138(i,z,a,xi,lo,hi,ncmax) - logical i1 - integer d,execnt,i,j,nc,ncmax - integer a(ncmax),hi(ncmax),lo(ncmax) - DOUBLE PRECISION xi(ncmax),z(8) - save execnt - data execnt /0/ - execnt=execnt+1 -c descend tree until leaf or ambiguous - j=i -c top of while loop - 3 if(a(j).ne.0)then - i1=(z(a(j)).ne.xi(j)) - else - i1=.false. - end if - if(.not.(i1))goto 4 - if(z(a(j)).le.xi(j))then - j=lo(j) - else - j=hi(j) - end if - goto 3 -c bottom of while loop - 4 ehg138=j - return - end - subroutine ehg106(il,ir,k,nk,p,pi,n) - integer execnt,i,ii,il,ir,j,k,l,n,nk,r - integer pi(n) - DOUBLE PRECISION t - DOUBLE PRECISION p(nk,n) - save execnt - data execnt /0/ - execnt=execnt+1 -c find the $k$-th smallest of $n$ elements -c Floyd+Rivest, CACM Mar '75, Algorithm 489 - l=il - r=ir -c top of while loop - 3 if(.not.(l.lt.r))goto 4 -c to avoid recursion, sophisticated partition deleted -c partition $x sub {l..r}$ about $t$ - t=p(1,pi(k)) - i=l - j=r - ii=pi(l) - pi(l)=pi(k) - pi(k)=ii - if(t.lt.p(1,pi(r)))then - ii=pi(l) - pi(l)=pi(r) - pi(r)=ii - end if -c top of while loop - 5 if(.not.(i.lt.j))goto 6 - ii=pi(i) - pi(i)=pi(j) - pi(j)=ii - i=i+1 - j=j-1 -c top of while loop - 7 if(.not.(p(1,pi(i)).lt.t))goto 8 - i=i+1 - goto 7 -c bottom of while loop - 8 continue -c top of while loop - 9 if(.not.(t.lt.p(1,pi(j))))goto 10 - j=j-1 - goto 9 -c bottom of while loop - 10 goto 5 -c bottom of while loop - 6 if(p(1,pi(l)).eq.t)then - ii=pi(l) - pi(l)=pi(j) - pi(j)=ii - else - j=j+1 - ii=pi(r) - pi(r)=pi(j) - pi(j)=ii - end if - if(j.le.k)then - l=j+1 - end if - if(k.le.j)then - r=j-1 - end if - goto 3 -c bottom of while loop - 4 return - end - subroutine ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,r - +cond,sing,sigma,u,e,dgamma,qraux,work,tol,dd,tdeg,cdeg,s) - integer column,d,dd,execnt,i,i3,i9,info,inorm2,j,jj,jpvt,k,kernel, - +n,nf,od,sing,tdeg - integer cdeg(8),psi(n) - double precision machep,f,i1,i10,i2,i4,i5,i6,i7,i8,rcond,rho,scal, - +tol - double precision g(15),sigma(15),u(15,15),e(15,15),b(nf,k),colnor( - +15),dist(n),eta(nf),dgamma(15),q(d),qraux(15),rw(n),s(0:od),w(nf), - +work(15),x(n,d),y(n) - external ehg106,ehg182,ehg184,dqrdc,dqrsl,dsvdc - integer idamax - external idamax - double precision d1mach - external d1mach - double precision ddot - external ddot - save machep,execnt - data execnt /0/ -c colnorm -> colnor -c E -> g -c MachEps -> machep -c V -> e -c X -> b - execnt=execnt+1 - if(execnt.eq.1)then - machep=d1mach(4) - end if -c sort by distance - do 3 i3=1,n - dist(i3)=0 - 3 continue - do 4 j=1,dd - i4=q(j) - do 5 i3=1,n - dist(i3)=dist(i3)+(x(i3,j)-i4)**2 - 5 continue - 4 continue - call ehg106(1,n,nf,1,dist,psi,n) - rho=dist(psi(nf))*max(1.d0,f) - if(.not.(0.lt.rho))then - call ehg182(120) - end if -c compute neighborhood weights - if(kernel.eq.2)then - do 6 i=1,nf - if(dist(psi(i)).lt.rho)then - i1=dsqrt(rw(psi(i))) - else - i1=0 - end if - w(i)=i1 - 6 continue - else - do 7 i3=1,nf - w(i3)=dsqrt(dist(psi(i3))/rho) - 7 continue - do 8 i3=1,nf - w(i3)=dsqrt(rw(psi(i3))*(1-w(i3)**3)**3) - 8 continue - end if - if(dabs(w(idamax(nf,w,1))).eq.0)then - call ehg184('at ',q,dd,1) - call ehg184('radius ',rho,1,1) - if(.not..false.)then - call ehg182(121) - end if - end if -c fill design matrix - column=1 - do 9 i3=1,nf - b(i3,column)=w(i3) - 9 continue - if(tdeg.ge.1)then - do 10 j=1,d - if(cdeg(j).ge.1)then - column=column+1 - i5=q(j) - do 11 i3=1,nf - b(i3,column)=w(i3)*(x(psi(i3),j)-i5) - 11 continue - end if - 10 continue - end if - if(tdeg.ge.2)then - do 12 j=1,d - if(cdeg(j).ge.1)then - if(cdeg(j).ge.2)then - column=column+1 - i6=q(j) - do 13 i3=1,nf - b(i3,column)=w(i3)*(x(psi(i3),j)-i6)**2 - 13 continue - end if - do 14 jj=j+1,d - if(cdeg(jj).ge.1)then - column=column+1 - i7=q(j) - i8=q(jj) - do 15 i3=1,nf - b(i3,column)=w(i3)*(x(psi(i3),j)-i7)*(x(psi(i3), - +jj)-i8) - 15 continue - end if - 14 continue - end if - 12 continue - k=column - end if - do 16 i3=1,nf - eta(i3)=w(i3)*y(psi(i3)) - 16 continue -c equilibrate columns - do 17 j=1,k - scal=0 - do 18 inorm2=1,nf - scal=scal+b(inorm2,j)**2 - 18 continue - scal=dsqrt(scal) - if(0.lt.scal)then - do 19 i3=1,nf - b(i3,j)=b(i3,j)/scal - 19 continue - colnor(j)=scal - else - colnor(j)=1 - end if - 17 continue -c singular value decomposition - call dqrdc(b,nf,nf,k,qraux,jpvt,work,0) - call dqrsl(b,nf,nf,k,qraux,eta,work,eta,eta,work,work,1000,info) - do 20 i9=1,k - do 21 i3=1,k - u(i3,i9)=0 - 21 continue - 20 continue - do 22 i=1,k - do 23 j=i,k - u(i,j)=b(i,j) - 23 continue - 22 continue - call dsvdc(u,15,k,k,sigma,g,u,15,e,15,work,21,info) - if(.not.(info.eq.0))then - call ehg182(182) - end if - tol=sigma(1)*(100*machep) - rcond=min(rcond,sigma(k)/sigma(1)) - if(sigma(k).le.tol)then - sing=sing+1 - if(sing.eq.1)then - call ehg184('Warning. pseudoinverse used at',q,d,1) - call ehg184('neighborhood radius',dsqrt(rho),1,1) - call ehg184('reciprocal condition number ',rcond,1,1) - else - if(sing.eq.2)then - call ehg184('There are other near singularities as well.' - +,rho,1,1) - end if - end if - end if -c compensate for equilibration - do 24 j=1,k - i10=colnor(j) - do 25 i3=1,k - e(j,i3)=e(j,i3)/i10 - 25 continue - 24 continue -c solve least squares problem - do 26 j=1,k - if(tol.lt.sigma(j))then - i2=ddot(k,u(1,j),1,eta,1)/sigma(j) - else - i2=0.d0 - end if - dgamma(j)=i2 - 26 continue - do 27 j=0,od -c bug fix 2006-07-04 for k=1, od>1. (thanks btyner@gmail.com) - if(j.lt.k)then - s(j)=ddot(k,e(j+1,1),15,dgamma,1) - else - s(j)=0.d0 - end if - 27 continue - return - end - subroutine ehg131(x,y,rw,trl,diagl,kernel,k,n,d,nc,ncmax,vc,nv,nvm - +ax,nf,f,a,c,hi,lo,pi,psi,v,vhit,vval,xi,dist,eta,b,ntol,fd,w,vval2 - +,rcond,sing,dd,tdeg,cdeg,lq,lf,setlf) - logical setlf - integer identi,d,dd,execnt,i1,i2,j,k,kernel,n,nc,ncmax,nf,ntol,nv, - +nvmax,sing,tdeg,vc - integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax),lo(ncm - +ax),pi(n),psi(n),vhit(nvmax) - double precision f,fd,rcond,trl - double precision lf(0:d,nvmax,nf),b(*),delta(8),diagl(n),dist(n),e - +ta(nf),rw(n),v(nvmax,d),vval(0:d,nvmax),vval2(0:d,nvmax),w(nf),x(n - +,d),xi(ncmax),y(n) - external ehg126,ehg182,ehg139,ehg124 - double precision dnrm2 - external dnrm2 - save execnt - data execnt /0/ -c Identity -> identi -c X -> b - execnt=execnt+1 - if(.not.(d.le.8))then - call ehg182(101) - end if -c build $k$-d tree - call ehg126(d,n,vc,x,v,nvmax) - nv=vc - nc=1 - do 3 j=1,vc - c(j,nc)=j - vhit(j)=0 - 3 continue - do 4 i1=1,d - delta(i1)=v(vc,i1)-v(1,i1) - 4 continue - fd=fd*dnrm2(d,delta,1) - do 5 identi=1,n - pi(identi)=identi - 5 continue - call ehg124(1,n,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v,vhit,nvmax, - +ntol,fd,dd) -c smooth - if(trl.ne.0)then - do 6 i2=1,nv - do 7 i1=0,d - vval2(i1,i2)=0 - 7 continue - 6 continue - end if - call ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k,dist,di - +st,eta,b,d,w,diagl,vval2,nc,vc,a,xi,lo,hi,c,vhit,rcond,sing,dd,tde - +g,cdeg,lq,lf,setlf,vval) - return - end - subroutine ehg133(n,d,vc,nvmax,nc,ncmax,a,c,hi,lo,v,vval,xi,m,z,s) - integer d,execnt,i,i1,m,nc,ncmax,nv,nvmax,vc - integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax) - double precision delta(8),s(m),v(nvmax,d),vval(0:d,nvmax),xi(ncmax - +),z(m,d) - double precision ehg128 - external ehg128 - save execnt - data execnt /0/ - execnt=execnt+1 - do 3 i=1,m - do 4 i1=1,d - delta(i1)=z(i,i1) - 4 continue - s(i)=ehg128(delta,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval) - 3 continue - return - end - subroutine ehg140(iw,i,j) - integer execnt,i,j - integer iw(i) - save execnt - data execnt /0/ - execnt=execnt+1 - iw(i)=j - return - end - subroutine ehg141(trl,n,deg,k,d,nsing,dk,delta1,delta2) - integer d,deg,dk,k,n,nsing - external ehg176 - double precision ehg176 - double precision corx,delta1,delta2,trl,z - double precision c(48), c1, c2, c3, c4 -c coef, d, deg, del - data c / .2971620d0,.3802660d0,.5886043d0,.4263766d0,.3346498d0,.6 - +271053d0,.5241198d0,.3484836d0,.6687687d0,.6338795d0,.4076457d0,.7 - +207693d0,.1611761d0,.3091323d0,.4401023d0,.2939609d0,.3580278d0,.5 - +555741d0,.3972390d0,.4171278d0,.6293196d0,.4675173d0,.4699070d0,.6 - +674802d0,.2848308d0,.2254512d0,.2914126d0,.5393624d0,.2517230d0,.3 - +898970d0,.7603231d0,.2969113d0,.4740130d0,.9664956d0,.3629838d0,.5 - +348889d0,.2075670d0,.2822574d0,.2369957d0,.3911566d0,.2981154d0,.3 - +623232d0,.5508869d0,.3501989d0,.4371032d0,.7002667d0,.4291632d0,.4 - +930370d0 / - if(deg.eq.0) dk=1 - if(deg.eq.1) dk=d+1 - if(deg.eq.2) dk=dfloat((d+2)*(d+1))/2.d0 - corx=dsqrt(k/dfloat(n)) - z=(dsqrt(k/trl)-corx)/(1-corx) - if(nsing .eq. 0 .and. 1 .lt. z) call ehg184('Chernobyl! trLn',trl,1,1) - z=min(1.0d0,max(0.0d0,z)) - c4=dexp(ehg176(z)) - i=1+3*(min(d,4)-1+4*(deg-1)) - if(d.le.4)then - c1=c(i) - c2=c(i+1) - c3=c(i+2) - else - c1=c(i)+(d-4)*(c(i)-c(i-3)) - c2=c(i+1)+(d-4)*(c(i+1)-c(i-2)) - c3=c(i+2)+(d-4)*(c(i+2)-c(i-1)) - endif - delta1=n-trl*dexp(c1*z**c2*(1-z)**c3*c4) - i=i+24 - if(d.le.4)then - c1=c(i) - c2=c(i+1) - c3=c(i+2) - else - c1=c(i)+(d-4)*(c(i)-c(i-3)) - c2=c(i+1)+(d-4)*(c(i+1)-c(i-2)) - c3=c(i+2)+(d-4)*(c(i+2)-c(i-1)) - endif - delta2=n-trl*dexp(c1*z**c2*(1-z)**c3*c4) - return - end - subroutine lowesc(n,l,ll,trl,delta1,delta2) - integer execnt,i,j,n - double precision delta1,delta2,trl - double precision l(n,n),ll(n,n) - double precision ddot - external ddot - save execnt - data execnt /0/ - execnt=execnt+1 -c compute $LL~=~(I-L)(I-L)'$ - do 3 i=1,n - l(i,i)=l(i,i)-1 - 3 continue - do 4 i=1,n - do 5 j=1,i - ll(i,j)=ddot(n,l(i,1),n,l(j,1),n) - 5 continue - 4 continue - do 6 i=1,n - do 7 j=i+1,n - ll(i,j)=ll(j,i) - 7 continue - 6 continue - do 8 i=1,n - l(i,i)=l(i,i)+1 - 8 continue -c accumulate first two traces - trl=0 - delta1=0 - do 9 i=1,n - trl=trl+l(i,i) - delta1=delta1+ll(i,i) - 9 continue -c $delta sub 2 = "tr" LL sup 2$ - delta2=0 - do 10 i=1,n - delta2=delta2+ddot(n,ll(i,1),n,ll(1,i),1) - 10 continue - return - end - subroutine ehg169(d,vc,nc,ncmax,nv,nvmax,v,a,xi,c,hi,lo) - integer d,execnt,i,j,k,mc,mv,nc,ncmax,nv,nvmax,p,vc - integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),novhit(1) - DOUBLE PRECISION v(nvmax,d),xi(ncmax) - external ehg125,ehg182 - integer ifloor - external ifloor - save execnt - data execnt /0/ - execnt=execnt+1 -c as in bbox -c remaining vertices - do 3 i=2,vc-1 - j=i-1 - do 4 k=1,d - v(i,k)=v(1+mod(j,2)*(vc-1),k) - j=ifloor(DFLOAT(j)/2.D0) - 4 continue - 3 continue -c as in ehg131 - mc=1 - mv=vc - novhit(1)=-1 - do 5 j=1,vc - c(j,mc)=j - 5 continue -c as in rbuild - p=1 -c top of while loop - 6 if(.not.(p.le.nc))goto 7 - if(a(p).ne.0)then - k=a(p) -c left son - mc=mc+1 - lo(p)=mc -c right son - mc=mc+1 - hi(p)=mc - call ehg125(p,mv,v,novhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k), - +c(1,p),c(1,lo(p)),c(1,hi(p))) - end if - p=p+1 - goto 6 -c bottom of while loop - 7 if(.not.(mc.eq.nc))then - call ehg182(193) - end if - if(.not.(mv.eq.nv))then - call ehg182(193) - end if - return - end - DOUBLE PRECISION function ehg176(z) - DOUBLE PRECISION z(*) - integer d,vc,nv,nc - integer a(17), c(2,17) - integer hi(17), lo(17) - DOUBLE PRECISION v(10,1) - DOUBLE PRECISION vval(0:1,10) - DOUBLE PRECISION xi(17) - DOUBLE PRECISION ehg128 - data d,vc,nv,nc /1,2,10,17/ - data a(1) /1/ - data hi(1),lo(1),xi(1) /3,2,0.3705D0/ - data c(1,1) /1/ - data c(2,1) /2/ - data a(2) /1/ - data hi(2),lo(2),xi(2) /5,4,0.2017D0/ - data c(1,2) /1/ - data c(2,2) /3/ - data a(3) /1/ - data hi(3),lo(3),xi(3) /7,6,0.5591D0/ - data c(1,3) /3/ - data c(2,3) /2/ - data a(4) /1/ - data hi(4),lo(4),xi(4) /9,8,0.1204D0/ - data c(1,4) /1/ - data c(2,4) /4/ - data a(5) /1/ - data hi(5),lo(5),xi(5) /11,10,0.2815D0/ - data c(1,5) /4/ - data c(2,5) /3/ - data a(6) /1/ - data hi(6),lo(6),xi(6) /13,12,0.4536D0/ - data c(1,6) /3/ - data c(2,6) /5/ - data a(7) /1/ - data hi(7),lo(7),xi(7) /15,14,0.7132D0/ - data c(1,7) /5/ - data c(2,7) /2/ - data a(8) /0/ - data c(1,8) /1/ - data c(2,8) /6/ - data a(9) /0/ - data c(1,9) /6/ - data c(2,9) /4/ - data a(10) /0/ - data c(1,10) /4/ - data c(2,10) /7/ - data a(11) /0/ - data c(1,11) /7/ - data c(2,11) /3/ - data a(12) /0/ - data c(1,12) /3/ - data c(2,12) /8/ - data a(13) /0/ - data c(1,13) /8/ - data c(2,13) /5/ - data a(14) /0/ - data c(1,14) /5/ - data c(2,14) /9/ - data a(15) /1/ - data hi(15),lo(15),xi(15) /17,16,0.8751D0/ - data c(1,15) /9/ - data c(2,15) /2/ - data a(16) /0/ - data c(1,16) /9/ - data c(2,16) /10/ - data a(17) /0/ - data c(1,17) /10/ - data c(2,17) /2/ - data vval(0,1) /-9.0572D-2/ - data v(1,1) /-5.D-3/ - data vval(1,1) /4.4844D0/ - data vval(0,2) /-1.0856D-2/ - data v(2,1) /1.005D0/ - data vval(1,2) /-0.7736D0/ - data vval(0,3) /-5.3718D-2/ - data v(3,1) /0.3705D0/ - data vval(1,3) /-0.3495D0/ - data vval(0,4) /2.6152D-2/ - data v(4,1) /0.2017D0/ - data vval(1,4) /-0.7286D0/ - data vval(0,5) /-5.8387D-2/ - data v(5,1) /0.5591D0/ - data vval(1,5) /0.1611D0/ - data vval(0,6) /9.5807D-2/ - data v(6,1) /0.1204D0/ - data vval(1,6) /-0.7978D0/ - data vval(0,7) /-3.1926D-2/ - data v(7,1) /0.2815D0/ - data vval(1,7) /-0.4457D0/ - data vval(0,8) /-6.4170D-2/ - data v(8,1) /0.4536D0/ - data vval(1,8) /3.2813D-2/ - data vval(0,9) /-2.0636D-2/ - data v(9,1) /0.7132D0/ - data vval(1,9) /0.3350D0/ - data vval(0,10) /4.0172D-2/ - data v(10,1) /0.8751D0/ - data vval(1,10) /-4.1032D-2/ - ehg176=ehg128(z,d,nc,vc,a,xi,lo,hi,c,v,nv,vval) - end - subroutine lowesa(trl,n,d,tau,nsing,delta1,delta2) - integer d,dka,dkb,execnt,n,nsing,tau - double precision alpha,d1a,d1b,d2a,d2b,delta1,delta2,trl - external ehg141 - save execnt - data execnt /0/ - execnt=execnt+1 - call ehg141(trl,n,1,tau,d,nsing,dka,d1a,d2a) - call ehg141(trl,n,2,tau,d,nsing,dkb,d1b,d2b) - alpha=dfloat(tau-dka)/dfloat(dkb-dka) - delta1=(1-alpha)*d1a+alpha*d1b - delta2=(1-alpha)*d2a+alpha*d2b - return - end - subroutine ehg191(m,z,l,d,n,nf,nv,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vv - +al2,lf,lq) - integer lq1,d,execnt,i,i1,i2,j,m,n,nc,ncmax,nf,nv,nvmax,p,vc - integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax) - double precision l(m,n),lf(0:d,nvmax,nf),v(nvmax,d),vval2(0:d,nvma - +x),xi(ncmax),z(m,d),zi(8) - double precision ehg128 - external ehg128 - save execnt - data execnt /0/ - execnt=execnt+1 - do 3 j=1,n - do 4 i2=1,nv - do 5 i1=0,d - vval2(i1,i2)=0 - 5 continue - 4 continue - do 6 i=1,nv -c linear search for i in Lq - lq1=lq(i,1) - lq(i,1)=j - p=nf -c top of while loop - 7 if(.not.(lq(i,p).ne.j))goto 8 - p=p-1 - goto 7 -c bottom of while loop - 8 lq(i,1)=lq1 - if(lq(i,p).eq.j)then - do 9 i1=0,d - vval2(i1,i)=lf(i1,i,p) - 9 continue - end if - 6 continue - do 10 i=1,m - do 11 i1=1,d - zi(i1)=z(i,i1) - 11 continue - l(i,j)=ehg128(zi,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval2) - 10 continue - 3 continue - return - end - subroutine ehg196(tau,d,f,trl) - integer d,dka,dkb,execnt,tau - double precision alpha,f,trl,trla,trlb - external ehg197 - save execnt - data execnt /0/ - execnt=execnt+1 - call ehg197(1,tau,d,f,dka,trla) - call ehg197(2,tau,d,f,dkb,trlb) - alpha=dfloat(tau-dka)/dfloat(dkb-dka) - trl=(1-alpha)*trla+alpha*trlb - return - end - subroutine ehg197(deg,tau,d,f,dk,trl) - integer d,deg,dk,tau - double precision trl, f - dk = 0 - if(deg.eq.1) dk=d+1 - if(deg.eq.2) dk=dfloat((d+2)*(d+1))/2.d0 - g1 = (-0.08125d0*d+0.13d0)*d+1.05d0 - trl = dk*(1+max(0.d0,(g1-f)/f)) - return - end - subroutine ehg192(y,d,n,nf,nv,nvmax,vval,lf,lq) - integer d,execnt,i,i1,i2,j,n,nf,nv,nvmax - integer lq(nvmax,nf) - DOUBLE PRECISION i3 - DOUBLE PRECISION lf(0:d,nvmax,nf),vval(0:d,nvmax),y(n) - save execnt - data execnt /0/ - execnt=execnt+1 - do 3 i2=1,nv - do 4 i1=0,d - vval(i1,i2)=0 - 4 continue - 3 continue - do 5 i=1,nv - do 6 j=1,nf - i3=y(lq(i,j)) - do 7 i1=0,d - vval(i1,i)=vval(i1,i)+i3*lf(i1,i,j) - 7 continue - 6 continue - 5 continue - return - end - DOUBLE PRECISION function ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax - +,vval) - logical i10,i2,i3,i4,i5,i6,i7,i8,i9 - integer d,execnt,i,i1,i11,i12,ig,ii,j,lg,ll,m,nc,ncmax,nt,nv,nvmax - +,ur,vc - integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),t(20) - DOUBLE PRECISION ge,gn,gs,gw,gpe,gpn,gps,gpw,h,phi0,phi1,psi0,psi1 - +,s,sew,sns,v0,v1,xibar - DOUBLE PRECISION g(0:8,256),g0(0:8),g1(0:8),v(nvmax,d),vval(0:d,nv - +max),xi(ncmax),z(d) - external ehg182,ehg184 - save execnt - data execnt /0/ - execnt=execnt+1 -c locate enclosing cell - nt=1 - t(nt)=1 - j=1 -c top of while loop - 3 if(.not.(a(j).ne.0))goto 4 - nt=nt+1 -c bug fix 2006-07-18 (thanks, btyner@gmail.com) - if(z(a(j)).le.xi(j))then - i1=lo(j) - else - i1=hi(j) - end if - t(nt)=i1 - if(.not.(nt.lt.20))then - call ehg182(181) - end if - j=t(nt) - goto 3 -c bottom of while loop - 4 continue -c tensor - do 5 i12=1,vc - do 6 i11=0,d - g(i11,i12)=vval(i11,c(i12,j)) - 6 continue - 5 continue - lg=vc - ll=c(1,j) - ur=c(vc,j) - do 7 i=d,1,-1 - h=(z(i)-v(ll,i))/(v(ur,i)-v(ll,i)) - if(h.lt.-.001D0)then - call ehg184('eval ',z,d,1) - call ehg184('lowerlimit ',v(ll,1),d,nvmax) - else - if(1.001D0.lt.h)then - call ehg184('eval ',z,d,1) - call ehg184('upperlimit ',v(ur,1),d,nvmax) - end if - end if - if(-.001D0.le.h)then - i2=(h.le.1.001D0) - else - i2=.false. - end if - if(.not.i2)then - call ehg182(122) - end if - lg=DFLOAT(lg)/2.D0 - do 8 ig=1,lg -c Hermite basis - phi0=(1-h)**2*(1+2*h) - phi1=h**2*(3-2*h) - psi0=h*(1-h)**2 - psi1=h**2*(h-1) - g(0,ig)=phi0*g(0,ig)+phi1*g(0,ig+lg)+(psi0*g(i,ig)+psi1*g(i, - +ig+lg))*(v(ur,i)-v(ll,i)) - do 9 ii=1,i-1 - g(ii,ig)=phi0*g(ii,ig)+phi1*g(ii,ig+lg) - 9 continue - 8 continue - 7 continue - s=g(0,1) -c blending - if(d.eq.2)then -c ----- North ----- - v0=v(ll,1) - v1=v(ur,1) - do 10 i11=0,d - g0(i11)=vval(i11,c(3,j)) - 10 continue - do 11 i11=0,d - g1(i11)=vval(i11,c(4,j)) - 11 continue - xibar=v(ur,2) - m=nt-1 -c top of while loop - 12 if(m.eq.0)then - i4=.true. - else - if(a(t(m)).eq.2)then - i3=(xi(t(m)).eq.xibar) - else - i3=.false. - end if - i4=i3 - end if - if(.not.(.not.i4))goto 13 - m=m-1 -c voidp junk - goto 12 -c bottom of while loop - 13 if(m.ge.1)then - m=hi(t(m)) -c top of while loop - 14 if(.not.(a(m).ne.0))goto 15 - if(z(a(m)).le.xi(m))then - m=lo(m) - else - m=hi(m) - end if - goto 14 -c bottom of while loop - 15 if(v0.lt.v(c(1,m),1))then - v0=v(c(1,m),1) - do 16 i11=0,d - g0(i11)=vval(i11,c(1,m)) - 16 continue - end if - if(v(c(2,m),1).lt.v1)then - v1=v(c(2,m),1) - do 17 i11=0,d - g1(i11)=vval(i11,c(2,m)) - 17 continue - end if - end if - h=(z(1)-v0)/(v1-v0) -c Hermite basis - phi0=(1-h)**2*(1+2*h) - phi1=h**2*(3-2*h) - psi0=h*(1-h)**2 - psi1=h**2*(h-1) - gn=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0) - gpn=phi0*g0(2)+phi1*g1(2) -c ----- South ----- - v0=v(ll,1) - v1=v(ur,1) - do 18 i11=0,d - g0(i11)=vval(i11,c(1,j)) - 18 continue - do 19 i11=0,d - g1(i11)=vval(i11,c(2,j)) - 19 continue - xibar=v(ll,2) - m=nt-1 -c top of while loop - 20 if(m.eq.0)then - i6=.true. - else - if(a(t(m)).eq.2)then - i5=(xi(t(m)).eq.xibar) - else - i5=.false. - end if - i6=i5 - end if - if(.not.(.not.i6))goto 21 - m=m-1 -c voidp junk - goto 20 -c bottom of while loop - 21 if(m.ge.1)then - m=lo(t(m)) -c top of while loop - 22 if(.not.(a(m).ne.0))goto 23 - if(z(a(m)).le.xi(m))then - m=lo(m) - else - m=hi(m) - end if - goto 22 -c bottom of while loop - 23 if(v0.lt.v(c(3,m),1))then - v0=v(c(3,m),1) - do 24 i11=0,d - g0(i11)=vval(i11,c(3,m)) - 24 continue - end if - if(v(c(4,m),1).lt.v1)then - v1=v(c(4,m),1) - do 25 i11=0,d - g1(i11)=vval(i11,c(4,m)) - 25 continue - end if - end if - h=(z(1)-v0)/(v1-v0) -c Hermite basis - phi0=(1-h)**2*(1+2*h) - phi1=h**2*(3-2*h) - psi0=h*(1-h)**2 - psi1=h**2*(h-1) - gs=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0) - gps=phi0*g0(2)+phi1*g1(2) -c ----- East ----- - v0=v(ll,2) - v1=v(ur,2) - do 26 i11=0,d - g0(i11)=vval(i11,c(2,j)) - 26 continue - do 27 i11=0,d - g1(i11)=vval(i11,c(4,j)) - 27 continue - xibar=v(ur,1) - m=nt-1 -c top of while loop - 28 if(m.eq.0)then - i8=.true. - else - if(a(t(m)).eq.1)then - i7=(xi(t(m)).eq.xibar) - else - i7=.false. - end if - i8=i7 - end if - if(.not.(.not.i8))goto 29 - m=m-1 -c voidp junk - goto 28 -c bottom of while loop - 29 if(m.ge.1)then - m=hi(t(m)) -c top of while loop - 30 if(.not.(a(m).ne.0))goto 31 - if(z(a(m)).le.xi(m))then - m=lo(m) - else - m=hi(m) - end if - goto 30 -c bottom of while loop - 31 if(v0.lt.v(c(1,m),2))then - v0=v(c(1,m),2) - do 32 i11=0,d - g0(i11)=vval(i11,c(1,m)) - 32 continue - end if - if(v(c(3,m),2).lt.v1)then - v1=v(c(3,m),2) - do 33 i11=0,d - g1(i11)=vval(i11,c(3,m)) - 33 continue - end if - end if - h=(z(2)-v0)/(v1-v0) -c Hermite basis - phi0=(1-h)**2*(1+2*h) - phi1=h**2*(3-2*h) - psi0=h*(1-h)**2 - psi1=h**2*(h-1) - ge=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0) - gpe=phi0*g0(1)+phi1*g1(1) -c ----- West ----- - v0=v(ll,2) - v1=v(ur,2) - do 34 i11=0,d - g0(i11)=vval(i11,c(1,j)) - 34 continue - do 35 i11=0,d - g1(i11)=vval(i11,c(3,j)) - 35 continue - xibar=v(ll,1) - m=nt-1 -c top of while loop - 36 if(m.eq.0)then - i10=.true. - else - if(a(t(m)).eq.1)then - i9=(xi(t(m)).eq.xibar) - else - i9=.false. - end if - i10=i9 - end if - if(.not.(.not.i10))goto 37 - m=m-1 -c voidp junk - goto 36 -c bottom of while loop - 37 if(m.ge.1)then - m=lo(t(m)) -c top of while loop - 38 if(.not.(a(m).ne.0))goto 39 - if(z(a(m)).le.xi(m))then - m=lo(m) - else - m=hi(m) - end if - goto 38 -c bottom of while loop - 39 if(v0.lt.v(c(2,m),2))then - v0=v(c(2,m),2) - do 40 i11=0,d - g0(i11)=vval(i11,c(2,m)) - 40 continue - end if - if(v(c(4,m),2).lt.v1)then - v1=v(c(4,m),2) - do 41 i11=0,d - g1(i11)=vval(i11,c(4,m)) - 41 continue - end if - end if - h=(z(2)-v0)/(v1-v0) -c Hermite basis - phi0=(1-h)**2*(1+2*h) - phi1=h**2*(3-2*h) - psi0=h*(1-h)**2 - psi1=h**2*(h-1) - gw=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0) - gpw=phi0*g0(1)+phi1*g1(1) -c NS - h=(z(2)-v(ll,2))/(v(ur,2)-v(ll,2)) -c Hermite basis - phi0=(1-h)**2*(1+2*h) - phi1=h**2*(3-2*h) - psi0=h*(1-h)**2 - psi1=h**2*(h-1) - sns=phi0*gs+phi1*gn+(psi0*gps+psi1*gpn)*(v(ur,2)-v(ll,2)) -c EW - h=(z(1)-v(ll,1))/(v(ur,1)-v(ll,1)) -c Hermite basis - phi0=(1-h)**2*(1+2*h) - phi1=h**2*(3-2*h) - psi0=h*(1-h)**2 - psi1=h**2*(h-1) - sew=phi0*gw+phi1*ge+(psi0*gpw+psi1*gpe)*(v(ur,1)-v(ll,1)) - s=(sns+sew)-s - end if - ehg128=s - return - end - integer function ifloor(x) - DOUBLE PRECISION x - ifloor=x - if(ifloor.gt.x) ifloor=ifloor-1 - end - DOUBLE PRECISION functionDSIGN(a1,a2) - DOUBLE PRECISION a1, a2 - DSIGN=DABS(a1) - if(a2.ge.0)DSIGN=-DSIGN - end - subroutine ehg136(u,lm,m,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,o - +d,o,ihat,w,rcond,sing,dd,tdeg,cdeg,s) - integer identi,d,dd,execnt,i,i1,ihat,info,j,k,kernel,l,lm,m,n,nf,o - +d,sing,tdeg - integer cdeg(8),psi(n) - double precision f,i2,rcond,scale,tol - double precision o(m,n),sigma(15),e(15,15),g(15,15),b(nf,k),dist(n - +),eta(nf),dgamma(15),q(8),qraux(15),rw(n),s(0:od,m),u(lm,d),w(nf), - +work(15),x(n,d),y(n) - external ehg127,ehg182,dqrsl - double precision ddot - external ddot - save execnt - data execnt /0/ -c V -> g -c U -> e -c Identity -> identi -c L -> o -c X -> b - execnt=execnt+1 - if(.not.(k.le.nf-1))then - call ehg182(104) - end if - if(.not.(k.le.15))then - call ehg182(105) - end if - do 3 identi=1,n - psi(identi)=identi - 3 continue - do 4 l=1,m - do 5 i1=1,d - q(i1)=u(l,i1) - 5 continue - call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,rcon - +d,sing,sigma,e,g,dgamma,qraux,work,tol,dd,tdeg,cdeg,s(0,l)) - if(ihat.eq.1)then -c $L sub {l,l} = -c V sub {1,:} SIGMA sup {+} U sup T -c (Q sup T W e sub i )$ - if(.not.(m.eq.n))then - call ehg182(123) - end if -c find $i$ such that $l = psi sub i$ - i=1 -c top of while loop - 6 if(.not.(l.ne.psi(i)))goto 7 - i=i+1 - if(.not.(i.lt.nf))then - call ehg182(123) - end if - goto 6 -c bottom of while loop - 7 do 8 i1=1,nf - eta(i1)=0 - 8 continue - eta(i)=w(i) -c $eta = Q sup T W e sub i$ - call dqrsl(b,nf,nf,k,qraux,eta,eta,eta,eta,eta,eta,1000,info - +) -c $gamma = U sup T eta sub {1:k}$ - do 9 i1=1,k - dgamma(i1)=0 - 9 continue - do 10 j=1,k - i2=eta(j) - do 11 i1=1,k - dgamma(i1)=dgamma(i1)+i2*e(j,i1) - 11 continue - 10 continue -c $gamma = SIGMA sup {+} gamma$ - do 12 j=1,k - if(tol.lt.sigma(j))then - dgamma(j)=dgamma(j)/sigma(j) - else - dgamma(j)=0.d0 - end if - 12 continue -c voidp junk -c voidp junk - o(l,1)=ddot(k,g(1,1),15,dgamma,1) - else - if(ihat.eq.2)then -c $L sub {l,:} = -c V sub {1,:} SIGMA sup {+} -c ( U sup T Q sup T ) W $ - do 13 i1=1,n - o(l,i1)=0 - 13 continue - do 14 j=1,k - do 15 i1=1,nf - eta(i1)=0 - 15 continue - do 16 i1=1,k - eta(i1)=e(i1,j) - 16 continue - call dqrsl(b,nf,nf,k,qraux,eta,eta,work,work,work,work - +,10000,info) - if(tol.lt.sigma(j))then - scale=1.d0/sigma(j) - else - scale=0.d0 - end if - do 17 i1=1,nf - eta(i1)=eta(i1)*(scale*w(i1)) - 17 continue - do 18 i=1,nf - o(l,psi(i))=o(l,psi(i))+g(1,j)*eta(i) - 18 continue - 14 continue - end if - end if - 4 continue - return - end - subroutine ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k,d - +ist,phi,eta,b,od,w,diagl,vval2,ncmax,vc,a,xi,lo,hi,c,vhit,rcond,si - +ng,dd,tdeg,cdeg,lq,lf,setlf,s) - logical setlf - integer identi,d,dd,execnt,i,i2,i3,i5,i6,ii,ileaf,info,j,k,kernel, - +l,n,nc,ncmax,nf,nleaf,nv,nvmax,od,sing,tdeg,vc - integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax),leaf(2 - +56),lo(ncmax),phi(n),pi(n),psi(n),vhit(nvmax) - DOUBLE PRECISION f,i1,i4,i7,rcond,scale,term,tol,trl - DOUBLE PRECISION lf(0:d,nvmax,nf),sigma(15),u(15,15),e(15,15),b(nf - +,k),diagl(n),dist(n),eta(nf),DGAMMA(15),q(8),qraux(15),rw(n),s(0:o - +d,nv),v(nvmax,d),vval2(0:d,nv),w(nf),work(15),x(n,d),xi(ncmax),y(n - +),z(8) - external ehg127,ehg182,DQRSL,ehg137 - DOUBLE PRECISION ehg128 - external ehg128 - DOUBLE PRECISION DDOT - external DDOT - save execnt - data execnt /0/ -c V -> e -c Identity -> identi -c X -> b - execnt=execnt+1 -c l2fit with trace(L) - if(.not.(k.le.nf-1))then - call ehg182(104) - end if - if(.not.(k.le.15))then - call ehg182(105) - end if - if(trl.ne.0)then - do 3 i5=1,n - diagl(i5)=0 - 3 continue - do 4 i6=1,nv - do 5 i5=0,d - vval2(i5,i6)=0 - 5 continue - 4 continue - end if - do 6 identi=1,n - psi(identi)=identi - 6 continue - do 7 l=1,nv - do 8 i5=1,d - q(i5)=v(l,i5) - 8 continue - call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,rcon - +d,sing,sigma,u,e,DGAMMA,qraux,work,tol,dd,tdeg,cdeg,s(0,l)) - if(trl.ne.0)then -c invert $psi$ - do 9 i5=1,n - phi(i5)=0 - 9 continue - do 10 i=1,nf - phi(psi(i))=i - 10 continue - do 11 i5=1,d - z(i5)=v(l,i5) - 11 continue - call ehg137(z,vhit(l),leaf,nleaf,d,nv,nvmax,ncmax,vc,a,xi,lo - +,hi,c,v) - do 12 ileaf=1,nleaf - do 13 ii=lo(leaf(ileaf)),hi(leaf(ileaf)) - i=phi(pi(ii)) - if(i.ne.0)then - if(.not.(psi(i).eq.pi(ii)))then - call ehg182(194) - end if - do 14 i5=1,nf - eta(i5)=0 - 14 continue - eta(i)=w(i) -c $eta = Q sup T W e sub i$ - call DQRSL(b,nf,nf,k,qraux,eta,work,eta,eta,work,wo - +rk,1000,info) - do 15 j=1,k - if(tol.lt.sigma(j))then - i4=DDOT(k,u(1,j),1,eta,1)/sigma(j) - else - i4=0.D0 - end if - DGAMMA(j)=i4 - 15 continue - do 16 j=1,d+1 -c bug fix 2006-07-15 for k=1, od>1. (thanks btyner@gmail.com) - if(j.le.k)then - vval2(j-1,l)=DDOT(k,e(j,1),15,DGAMMA,1) - else - vval2(j-1,l)=0 - end if - 16 continue - do 17 i5=1,d - z(i5)=x(pi(ii),i5) - 17 continue - term=ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval2 - +) - diagl(pi(ii))=diagl(pi(ii))+term - do 18 i5=0,d - vval2(i5,l)=0 - 18 continue - end if - 13 continue - 12 continue - end if - if(setlf)then -c $Lf sub {:,l,:} = V SIGMA sup {+} U sup T Q sup T W$ - if(.not.(k.ge.d+1))then - call ehg182(196) - end if - do 19 i5=1,nf - lq(l,i5)=psi(i5) - 19 continue - do 20 i6=1,nf - do 21 i5=0,d - lf(i5,l,i6)=0 - 21 continue - 20 continue - do 22 j=1,k - do 23 i5=1,nf - eta(i5)=0 - 23 continue - do 24 i5=1,k - eta(i5)=u(i5,j) - 24 continue - call DQRSL(b,nf,nf,k,qraux,eta,eta,work,work,work,work,10 - +000,info) - if(tol.lt.sigma(j))then - scale=1.D0/sigma(j) - else - scale=0.D0 - end if - do 25 i5=1,nf - eta(i5)=eta(i5)*(scale*w(i5)) - 25 continue - do 26 i=1,nf - i7=eta(i) - do 27 i5=0,d - if(i5.lt.k)then - lf(i5,l,i)=lf(i5,l,i)+e(1+i5,j)*i7 - else - lf(i5,l,i)=0 - end if - 27 continue - 26 continue - 22 continue - end if - 7 continue - if(trl.ne.0)then - if(n.le.0)then - trl=0.D0 - else - i3=n - i1=diagl(i3) - do 28 i2=i3-1,1,-1 - i1=diagl(i2)+i1 - 28 continue - trl=i1 - end if - end if - return - end - subroutine dqrdc(x,ldx,n,p,qraux,jpvt,work,job) - integer ldx,n,p,job - integer jpvt(1) - double precision x(ldx,1),qraux(1),work(1) -c -c dqrdc uses householder transformations to compute the qr -c factorization of an n by p matrix x. column pivoting -c based on the 2-norms of the reduced columns may be -c performed at the users option. -c -c on entry -c -c x double precision(ldx,p), where ldx .ge. n. -c x contains the matrix whose decomposition is to be -c computed. -c -c ldx integer. -c ldx is the leading dimension of the array x. -c -c n integer. -c n is the number of rows of the matrix x. -c -c p integer. -c p is the number of columns of the matrix x. -c -c jpvt integer(p). -c jpvt contains integers that control the selection -c of the pivot columns. the k-th column x(k) of x -c is placed in one of three classes according to the -c value of jpvt(k). -c -c if jpvt(k) .gt. 0, then x(k) is an initial -c column. -c -c if jpvt(k) .eq. 0, then x(k) is a free column. -c -c if jpvt(k) .lt. 0, then x(k) is a final column. -c -c before the decomposition is computed, initial columns -c are moved to the beginning of the array x and final -c columns to the end. both initial and final columns -c are frozen in place during the computation and only -c free columns are moved. at the k-th stage of the -c reduction, if x(k) is occupied by a free column -c it is interchanged with the free column of largest -c reduced norm. jpvt is not referenced if -c job .eq. 0. -c -c work double precision(p). -c work is a work array. work is not referenced if -c job .eq. 0. -c -c job integer. -c job is an integer that initiates column pivoting. -c if job .eq. 0, no pivoting is done. -c if job .ne. 0, pivoting is done. -c -c on return -c -c x x contains in its upper triangle the upper -c triangular matrix r of the qr factorization. -c below its diagonal x contains information from -c which the orthogonal part of the decomposition -c can be recovered. note that if pivoting has -c been requested, the decomposition is not that -c of the original matrix x but that of x -c with its columns permuted as described by jpvt. -c -c qraux double precision(p). -c qraux contains further information required to recover -c the orthogonal part of the decomposition. -c -c jpvt jpvt(k) contains the index of the column of the -c original matrix that has been interchanged into -c the k-th column, if pivoting was requested. -c -c linpack. this version dated 08/14/78 . -c g.w. stewart, university of maryland, argonne national lab. -c -c dqrdc uses the following functions and subprograms. -c -c blas daxpy,ddot,dscal,dswap,dnrm2 -c fortran dabs,dmax1,min0,dsqrt -c -c internal variables -c - integer j,jp,l,lp1,lup,maxj,pl,pu - double precision maxnrm,dnrm2,tt - double precision ddot,nrmxl,t - logical negj,swapj -c -c - pl = 1 - pu = 0 - if (job .eq. 0) go to 60 -c -c pivoting has been requested. rearrange the columns -c according to jpvt. -c - do 20 j = 1, p - swapj = jpvt(j) .gt. 0 - negj = jpvt(j) .lt. 0 - jpvt(j) = j - if (negj) jpvt(j) = -j - if (.not.swapj) go to 10 - if (j .ne. pl) call dswap(n,x(1,pl),1,x(1,j),1) - jpvt(j) = jpvt(pl) - jpvt(pl) = j - pl = pl + 1 - 10 continue - 20 continue - pu = p - do 50 jj = 1, p - j = p - jj + 1 - if (jpvt(j) .ge. 0) go to 40 - jpvt(j) = -jpvt(j) - if (j .eq. pu) go to 30 - call dswap(n,x(1,pu),1,x(1,j),1) - jp = jpvt(pu) - jpvt(pu) = jpvt(j) - jpvt(j) = jp - 30 continue - pu = pu - 1 - 40 continue - 50 continue - 60 continue -c -c compute the norms of the free columns. -c - if (pu .lt. pl) go to 80 - do 70 j = pl, pu - qraux(j) = dnrm2(n,x(1,j),1) - work(j) = qraux(j) - 70 continue - 80 continue -c -c perform the householder reduction of x. -c - lup = min0(n,p) - do 200 l = 1, lup - if (l .lt. pl .or. l .ge. pu) go to 120 -c -c locate the column of largest norm and bring it -c into the pivot position. -c - maxnrm = 0.0d0 - maxj = l - do 100 j = l, pu - if (qraux(j) .le. maxnrm) go to 90 - maxnrm = qraux(j) - maxj = j - 90 continue - 100 continue - if (maxj .eq. l) go to 110 - call dswap(n,x(1,l),1,x(1,maxj),1) - qraux(maxj) = qraux(l) - work(maxj) = work(l) - jp = jpvt(maxj) - jpvt(maxj) = jpvt(l) - jpvt(l) = jp - 110 continue - 120 continue - qraux(l) = 0.0d0 - if (l .eq. n) go to 190 -c -c compute the householder transformation for column l. -c - nrmxl = dnrm2(n-l+1,x(l,l),1) - if (nrmxl .eq. 0.0d0) go to 180 - if (x(l,l) .ne. 0.0d0) nrmxl = dsign(nrmxl,x(l,l)) - call dscal(n-l+1,1.0d0/nrmxl,x(l,l),1) - x(l,l) = 1.0d0 + x(l,l) -c -c apply the transformation to the remaining columns, -c updating the norms. -c - lp1 = l + 1 - if (p .lt. lp1) go to 170 - do 160 j = lp1, p - t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) - call daxpy(n-l+1,t,x(l,l),1,x(l,j),1) - if (j .lt. pl .or. j .gt. pu) go to 150 - if (qraux(j) .eq. 0.0d0) go to 150 - tt = 1.0d0 - (dabs(x(l,j))/qraux(j))**2 - tt = dmax1(tt,0.0d0) - t = tt - tt = 1.0d0 + 0.05d0*tt*(qraux(j)/work(j))**2 - if (tt .eq. 1.0d0) go to 130 - qraux(j) = qraux(j)*dsqrt(t) - go to 140 - 130 continue - qraux(j) = dnrm2(n-l,x(l+1,j),1) - work(j) = qraux(j) - 140 continue - 150 continue - 160 continue - 170 continue -c -c save the transformation. -c - qraux(l) = x(l,l) - x(l,l) = -nrmxl - 180 continue - 190 continue - 200 continue - return - end - integer function idamax(n,dx,incx) -c -c finds the index of element having max. absolute value. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dmax - integer i,incx,ix,n -c - idamax = 0 - if( n .lt. 1 ) return - idamax = 1 - if(n.eq.1)return - if(incx.eq.1)go to 20 -c -c code for increment not equal to 1 -c - ix = 1 - dmax = dabs(dx(1)) - ix = ix + incx - do 10 i = 2,n - if(dabs(dx(ix)).le.dmax) go to 5 - idamax = i - dmax = dabs(dx(ix)) - 5 ix = ix + incx - 10 continue - return -c -c code for increment equal to 1 -c - 20 dmax = dabs(dx(1)) - do 30 i = 2,n - if(dabs(dx(i)).le.dmax) go to 30 - idamax = i - dmax = dabs(dx(i)) - 30 continue - return - end - subroutine lowesb(xx,yy,ww,diagl,infl,iv,liv,lv,wv) - logical infl,setlf - integer execnt - integer iv(*) - DOUBLE PRECISION trl - DOUBLE PRECISION diagl(*),wv(*),ww(*),xx(*),yy(*) - external ehg131,ehg182,ehg183 - integer ifloor - external ifloor - save execnt - data execnt /0/ - execnt=execnt+1 - if(.not.(iv(28).ne.173))then - call ehg182(174) - end if - if(iv(28).ne.172)then - if(.not.(iv(28).eq.171))then - call ehg182(171) - end if - end if - iv(28)=173 - if(infl)then - trl=1.D0 - else - trl=0.D0 - end if - setlf=(iv(27).ne.iv(25)) - call ehg131(xx,yy,ww,trl,diagl,iv(20),iv(29),iv(3),iv(2),iv(5),iv( - +17),iv(4),iv(6),iv(14),iv(19),wv(1),iv(iv(7)),iv(iv(8)),iv(iv(9)), - +iv(iv(10)),iv(iv(22)),iv(iv(27)),wv(iv(11)),iv(iv(23)),wv(iv(13)), - +wv(iv(12)),wv(iv(15)),wv(iv(16)),wv(iv(18)),ifloor(iv(3)*wv(2)),wv - +(3),wv(iv(26)),wv(iv(24)),wv(4),iv(30),iv(33),iv(32),iv(41),iv(iv( - +25)),wv(iv(34)),setlf) - if(iv(14).lt.iv(6)+DFLOAT(iv(4))/2.D0)then - call ehg183('Warning. k-d tree limited by memory; nvmax=',iv(14 - +),1,1) - else - if(iv(17).lt.iv(5)+2)then - call ehg183('Warning. k-d tree limited by memory. ncmax=',iv - +(17),1,1) - end if - end if - return - end - subroutine lowesd(versio,iv,liv,lv,v,d,n,f,ideg,nvmax,setlf) - logical setlf - integer bound,d,execnt,i,i1,i2,ideg,j,liv,lv,n,ncmax,nf,nvmax,vc,v - +ersio - integer iv(liv) - double precision f - double precision v(lv) - external ehg182 - integer ifloor - external ifloor - save execnt - data execnt /0/ -c version -> versio - execnt=execnt+1 - if(.not.(versio.eq.106))then - call ehg182(100) - end if - iv(28)=171 - iv(2)=d - iv(3)=n - vc=2**d - iv(4)=vc - if(.not.(0.lt.f))then - call ehg182(120) - end if - nf=min(n,ifloor(n*f)) - iv(19)=nf - iv(20)=1 - if(ideg.eq.0)then - i1=1 - else - if(ideg.eq.1)then - i1=d+1 - else - if(ideg.eq.2)then - i1=dfloat((d+2)*(d+1))/2.d0 - end if - end if - end if - iv(29)=i1 - iv(21)=1 - iv(14)=nvmax - ncmax=nvmax - iv(17)=ncmax - iv(30)=0 - iv(32)=ideg - if(.not.(ideg.ge.0))then - call ehg182(195) - end if - if(.not.(ideg.le.2))then - call ehg182(195) - end if - iv(33)=d - do 3 i2=41,49 - iv(i2)=ideg - 3 continue - iv(7)=50 - iv(8)=iv(7)+ncmax - iv(9)=iv(8)+vc*ncmax - iv(10)=iv(9)+ncmax - iv(22)=iv(10)+ncmax -c initialize permutation - j=iv(22)-1 - do 4 i=1,n - iv(j+i)=i - 4 continue - iv(23)=iv(22)+n - iv(25)=iv(23)+nvmax - if(setlf)then - iv(27)=iv(25)+nvmax*nf - else - iv(27)=iv(25) - end if - bound=iv(27)+n - if(.not.(bound-1.le.liv))then - call ehg182(102) - end if - iv(11)=50 - iv(13)=iv(11)+nvmax*d - iv(12)=iv(13)+(d+1)*nvmax - iv(15)=iv(12)+ncmax - iv(16)=iv(15)+n - iv(18)=iv(16)+nf - iv(24)=iv(18)+iv(29)*nf - iv(34)=iv(24)+(d+1)*nvmax - if(setlf)then - iv(26)=iv(34)+(d+1)*nvmax*nf - else - iv(26)=iv(34) - end if - bound=iv(26)+nf - if(.not.(bound-1.le.lv))then - call ehg182(103) - end if - v(1)=f - v(2)=0.05d0 - v(3)=0.d0 - v(4)=1.d0 - return - end - subroutine lowese(iv,liv,lv,wv,m,z,s) - integer execnt,m - integer iv(*) - double precision s(m),wv(*),z(m,1) - external ehg133,ehg182 - save execnt - data execnt /0/ - execnt=execnt+1 - if(.not.(iv(28).ne.172))then - call ehg182(172) - end if - if(.not.(iv(28).eq.173))then - call ehg182(173) - end if - call ehg133(iv(3),iv(2),iv(4),iv(14),iv(5),iv(17),iv(iv(7)),iv(iv( - +8)),iv(iv(9)),iv(iv(10)),wv(iv(11)),wv(iv(13)),wv(iv(12)),m,z,s) - return - end - subroutine lowesf(xx,yy,ww,iv,liv,lv,wv,m,z,l,ihat,s) - logical i1 - integer execnt,ihat,m,n - integer iv(*) - double precision l(m,*),s(m),wv(*),ww(*),xx(*),yy(*),z(m,1) - external ehg182,ehg136 - save execnt - data execnt /0/ - execnt=execnt+1 - if(171.le.iv(28))then - i1=(iv(28).le.174) - else - i1=.false. - end if - if(.not.i1)then - call ehg182(171) - end if - iv(28)=172 - if(.not.(iv(14).ge.iv(19)))then - call ehg182(186) - end if - call ehg136(z,m,m,iv(3),iv(2),iv(19),wv(1),xx,iv(iv(22)),yy,ww,iv( - +20),iv(29),wv(iv(15)),wv(iv(16)),wv(iv(18)),0,l,ihat,wv(iv(26)),wv - +(4),iv(30),iv(33),iv(32),iv(41),s) - return - end - subroutine lowesl(iv,liv,lv,wv,m,z,l) - integer execnt,m,n - integer iv(*) - double precision l(m,*),wv(*),z(m,1) - external ehg182,ehg191 - save execnt - data execnt /0/ - execnt=execnt+1 - if(.not.(iv(28).ne.172))then - call ehg182(172) - end if - if(.not.(iv(28).eq.173))then - call ehg182(173) - end if - if(.not.(iv(26).ne.iv(34)))then - call ehg182(175) - end if - call ehg191(m,z,l,iv(2),iv(3),iv(19),iv(6),iv(17),iv(4),iv(iv(7)), - +wv(iv(12)),iv(iv(10)),iv(iv(9)),iv(iv(8)),wv(iv(11)),iv(14),wv(iv( - +24)),wv(iv(34)),iv(iv(25))) - return - end - subroutine lowesr(yy,iv,liv,lv,wv) - integer execnt - integer iv(*) - DOUBLE PRECISION wv(*),yy(*) - external ehg182,ehg192 - save execnt - data execnt /0/ - execnt=execnt+1 - if(.not.(iv(28).ne.172))then - call ehg182(172) - end if - if(.not.(iv(28).eq.173))then - call ehg182(173) - end if - call ehg192(yy,iv(2),iv(3),iv(19),iv(6),iv(14),wv(iv(13)),wv(iv(34 - +)),iv(iv(25))) - return - end - subroutine lowesw(res,n,rw,pi) - integer identi,execnt,i,i1,n,nh - integer pi(n) - double precision cmad,rsmall - double precision res(n),rw(n) - external ehg106 - integer ifloor - external ifloor - double precision d1mach - external d1mach - save execnt - data execnt /0/ -c Identity -> identi - execnt=execnt+1 -c tranliterated from Devlin's ratfor -c find median of absolute residuals - do 3 i1=1,n - rw(i1)=dabs(res(i1)) - 3 continue - do 4 identi=1,n - pi(identi)=identi - 4 continue - nh=ifloor(dfloat(n)/2.d0)+1 -c partial sort to find 6*mad - call ehg106(1,n,nh,1,rw,pi,n) - if((n-nh)+1.lt.nh)then - call ehg106(1,nh-1,nh-1,1,rw,pi,n) - cmad=3*(rw(pi(nh))+rw(pi(nh-1))) - else - cmad=6*rw(pi(nh)) - end if - rsmall=d1mach(1) - if(cmad.lt.rsmall)then - do 5 i1=1,n - rw(i1)=1 - 5 continue - else - do 6 i=1,n - if(cmad*0.999d0.lt.rw(i))then - rw(i)=0 - else - if(cmad*0.001d0.lt.rw(i))then - rw(i)=(1-(rw(i)/cmad)**2)**2 - else - rw(i)=1 - end if - end if - 6 continue - end if - return - end - subroutine lowesp(n,y,yhat,pwgts,rwgts,pi,ytilde) - integer identi,execnt,i2,i3,i5,m,n - integer pi(n) - double precision c,i1,i4,mad - double precision pwgts(n),rwgts(n),y(n),yhat(n),ytilde(n) - external ehg106 - integer ifloor - external ifloor - save execnt - data execnt /0/ -c Identity -> identi - execnt=execnt+1 -c median absolute deviation - do 3 i5=1,n - ytilde(i5)=dabs(y(i5)-yhat(i5))*dsqrt(pwgts(i5)) - 3 continue - do 4 identi=1,n - pi(identi)=identi - 4 continue - m=ifloor(dfloat(n)/2.d0)+1 - call ehg106(1,n,m,1,ytilde,pi,n) - if((n-m)+1.lt.m)then - call ehg106(1,m-1,m-1,1,ytilde,pi,n) - mad=(ytilde(pi(m-1))+ytilde(pi(m)))/2 - else - mad=ytilde(pi(m)) - end if -c magic constant - c=(6*mad)**2/5 - do 5 i5=1,n - ytilde(i5)=1-((y(i5)-yhat(i5))**2*pwgts(i5))/c - 5 continue - do 6 i5=1,n - ytilde(i5)=ytilde(i5)*dsqrt(rwgts(i5)) - 6 continue - if(n.le.0)then - i4=0.d0 - else - i3=n - i1=ytilde(i3) - do 7 i2=i3-1,1,-1 - i1=ytilde(i2)+i1 - 7 continue - i4=i1 - end if - c=n/i4 -c pseudovalues - do 8 i5=1,n - ytilde(i5)=yhat(i5)+(c*rwgts(i5))*(y(i5)-yhat(i5)) - 8 continue - return - end - subroutine ehg124(ll,uu,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v,vhi - +t,nvmax,fc,fd,dd) - logical i1,i2,i3,leaf - integer d,dd,execnt,fc,i4,inorm2,k,l,ll,m,n,nc,ncmax,nv,nvmax,p,u, - +uu,vc,lower,upper,check,offset - integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),pi(n),vhit(nvmax) - DOUBLE PRECISION diam,fd - DOUBLE PRECISION diag(8),sigma(8),v(nvmax,d),x(n,d),xi(ncmax) - external ehg125,ehg106,ehg129 - integer IDAMAX - external IDAMAX - save execnt - data execnt /0/ - execnt=execnt+1 - p=1 - l=ll - u=uu - lo(p)=l - hi(p)=u -c top of while loop - 3 if(.not.(p.le.nc))goto 4 - do 5 i4=1,dd - diag(i4)=v(c(vc,p),i4)-v(c(1,p),i4) - 5 continue - diam=0 - do 6 inorm2=1,dd - diam=diam+diag(inorm2)**2 - 6 continue - diam=DSQRT(diam) - if((u-l)+1.le.fc)then - i1=.true. - else - i1=(diam.le.fd) - end if - if(i1)then - leaf=.true. - else - if(ncmax.lt.nc+2)then - i2=.true. - else - i2=(nvmax.lt.nv+DFLOAT(vc)/2.D0) - end if - leaf=i2 - end if - if(.not.leaf)then - call ehg129(l,u,dd,x,pi,n,sigma) - k=IDAMAX(dd,sigma,1) - m=DFLOAT(l+u)/2.D0 - call ehg106(l,u,m,1,x(1,k),pi,n) - -c bug fix from btyner@gmail.com 2006-07-20 - offset = 0 - 7 if(((m+offset).ge.u).or.((m+offset).lt.l))then - goto 8 - else - if(offset.lt.0)then - lower = l - check = m + offset - upper = check - else - lower = m + offset + 1 - check = lower - upper = u - end if - call ehg106(lower,upper,check,1,x(1,k),pi,n) - if(x(pi(m + offset),k).eq.x(pi(m+offset+1),k))then - offset = (-offset) - if(offset.ge.0)then - offset = offset + 1 - end if - goto 7 - else - m = m + offset - goto 8 - end if - end if - - 8 if(v(c(1,p),k).eq.x(pi(m),k))then - leaf=.true. - else - leaf=(v(c(vc,p),k).eq.x(pi(m),k)) - end if - end if - if(leaf)then - a(p)=0 - else - a(p)=k - xi(p)=x(pi(m),k) -c left son - nc=nc+1 - lo(p)=nc - lo(nc)=l - hi(nc)=m -c right son - nc=nc+1 - hi(p)=nc - lo(nc)=m+1 - hi(nc)=u - call ehg125(p,nv,v,vhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k),c( - +1,p),c(1,lo(p)),c(1,hi(p))) - end if - p=p+1 - l=lo(p) - u=hi(p) - goto 3 -c bottom of while loop - 4 return - end - subroutine ehg129(l,u,d,x,pi,n,sigma) - integer d,execnt,i,k,l,n,u - integer pi(n) - DOUBLE PRECISION machin,alpha,beta,t - DOUBLE PRECISION sigma(d),x(n,d) - DOUBLE PRECISION D1MACH - external D1MACH - save machin,execnt - data execnt /0/ -c MachInf -> machin - execnt=execnt+1 - if(execnt.eq.1)then - machin=D1MACH(2) - end if - do 3 k=1,d - alpha=machin - beta=-machin - do 4 i=l,u - t=x(pi(i),k) - alpha=min(alpha,x(pi(i),k)) - beta=max(beta,t) - 4 continue - sigma(k)=beta-alpha - 3 continue - return - end - subroutine ehg137(z,kappa,leaf,nleaf,d,nv,nvmax,ncmax,vc,a,xi,lo,h - +i,c,v) - integer d,execnt,nc,ncmax,nleaf,p,stackt - integer a(ncmax),hi(ncmax),leaf(256),lo(ncmax),pstack(20) - DOUBLE PRECISION xi(ncmax),z(d) - external ehg182 - save execnt - data execnt /0/ -c stacktop -> stackt - execnt=execnt+1 -c find leaf cells affected by $z$ - stackt=0 - p=1 - nleaf=0 -c top of while loop - 3 if(.not.(0.lt.p))goto 4 - if(a(p).eq.0)then -c leaf - nleaf=nleaf+1 - leaf(nleaf)=p -c Pop - if(stackt.ge.1)then - p=pstack(stackt) - else - p=0 - end if - stackt=max(0,stackt-1) - else - if(z(a(p)).eq.xi(p))then -c Push - stackt=stackt+1 - if(.not.(stackt.le.20))then - call ehg182(187) - end if - pstack(stackt)=hi(p) - p=lo(p) - else - if(z(a(p)).le.xi(p))then - p=lo(p) - else - p=hi(p) - end if - end if - end if - goto 3 -c bottom of while loop - 4 if(.not.(nleaf.le.256))then - call ehg182(185) - end if - return - end //GO.SYSIN DD loessf.f echo loessf.m 1>&2 sed >loessf.m <<'//GO.SYSIN DD loessf.m' 's/^-//' -*************************************************************** -* LOESS smoothing scattered data in one or more variables * -* documentation of Fortran routines * -* Cleveland, Devlin, Grosse, Shyu * -*************************************************************** - -1. The typical program would call lowesd, set tolerances in iv,v if - desired, then call lowesb and lowese. -2. To save the k-d tree, call lowesd, lowesb and then losave; subsequent - programs would call lohead, set liv and lv, then call loread and lowese. -3. For statistics, get diagL and then call lowesa or get the full hat - matrix and call lowesc. Robustness iterations can take advantage of - lowesw and lowesp. - -lowesd(106,iv,liv,lv,v,d,n,f,tdeg,nvmax,setLf) setup workspace -lowesf(x,y,w,iv,liv,lv,v,m,z,L,hat,s) slow smooth at z -lowesb(x,y,w,diagL,infl,iv,liv,lv,v) build k-d tree -lowesr(y,iv,liv,lv,v) rebuild with new data values - (does not change y) -lowese(iv,liv,lv,v,m,z, s) evaluate smooth at z -lowesl(iv,liv,lv,v,m,z, L) explicit hat matrix, - which maps from y to z -lofort(iunit,iv,liv,lv,v) save k-d tree as Fortran -losave(iunit,iv,liv,lv,v) save k-d tree in file -lohead(iunit,d,vc,nc,nv) read d,vc,nc,nv from file - liv = 50+(vc+3)*nc determine space - lv = 50+(2*d+1)*nv+nc requirements -loread(iunit,d,vc,nc,nv,iv,liv,lv,v) finish reading k-d tree, - ready for lowese -lowesa(trL,n,d,tau,nsing, del1,del2) approximate delta -lowesc(n,L,LL, trL,del1,del2) exact delta -lowesp(n,y,yhat,w,rw, pi,ytilde) pseudo-values -lowesw(res,n, rw,pi) robustness weights - -=== arguments === -d number of independent variables [integer] (called "p" elsewhere) -del1,del2 delta1, delta2 -diagL diagonal of hat matrix, only set if infl=.true. (n) -f fraction of points to use in local smooth (called "alpha" elsewhere) -fc don't refine cells with less than fc*n points; ordinarily=.05 -hat is hat matrix desired? [integer] - 0 = none - 1 = diagonal only - 2 = full matrix -infl is diagonal of hat matrix desired? [logical] -iunit Fortran unit number for i/o -iv workspace (liv) -L hat matrix (m,n) [real] - in lowesf, only computed if hat nonzero; if hat=1 only size (n) -LL workspace (n,n) -liv 50+(2^d+4)*nvmax+2*n - if setLf, add nf*nvmax -lv 50+(3*d+3)*nvmax+n+(tau0+2)*nf - if setLf, add (d+1)*nf*nvmax -m number of points to smooth at; ordinarily=n -n number of observations -nf min(n,floor(n*f)) -nsing if 0, print warning in lowesa when trL0. -4 RCOND reciprocal condition number -... 49 reserved for future use -iv(V1) v vertices (nv,d) -iv(VV1) vval vertex values (0:d,nv) -iv(XI1) xi cut values (nc) -------------------------eval only needs workspace up to here -iv(WORK1) workspace (n) l2fit:dist -iv(WORK2) workspace (nf) l2fit:eta -iv(WORK3) workspace (dim,nf) l2fit:X -iv(VV2) vval2 workspace ((d+1)*nv) pseudo-vval for trL computation -iv(LF) Lf hat matrix (data to vertex) ((d+1)*nvmax*nf) -iv(WORK4) workspace (nf) l2fit:w - -Internal routine names have been hidden as follows: -ehg106 select q-th smallest by partial sorting -ehg124 rbuild -ehg125 cpvert -ehg126 bbox -ehg127 l2fit,l2tr computational kernel -ehg128 eval -ehg129 spread -ehg131 lowesb after workspace expansion -ehg133 lowese after workspace expansion -ehg134 abort by calling S Recover function -ehg136 l2fit with hat matrix L -ehg137 vleaf -ehg138 descend -ehg139 l2tr -ehg140(w,i,j) w(i)=j used when w is declared real, but should store an int -ehg141 delta1,2 from trL -ehg142 robust iteration -ehg144 now called lowesc -ehg152 like ehg142, but for lowesf -ehg167 kernel for losave -ehg168 kernel for loread -ehg169 compute derived k-d tree information -ehg170 generate Fortran -ehg176,ehg177,ehg178,ehg179,ehg180,ehg181 loeval for delta -ehg182 ehgdie(i) -ehg183 warning(message,i,n,inc) -ehg184 warning(message,x,n,inc) -ehg190 now called lowesa, with slight change in calling sequence -ehg191 lowesl after workspace expansion -ehg192 lowesr after workspace expansion -ehg196(tau,d,f,trl) trL approximation -ehg197 for deg=1,2 -m9rwt now called lowesw -pseudo now called lowesp //GO.SYSIN DD loessf.m echo madeup.c 1>&2 sed >madeup.c <<'//GO.SYSIN DD madeup.c' 's/^-//' -#include -#include "loess.h" - -struct loess_struct madeup, madeup_new, madeup2; -struct pred_struct madeup_pred; -struct ci_struct madeup_ci; -struct anova_struct madeup_anova; -double one_two[] = {-0.957581198938384, -2.80954937859791, -0.696510845605909, - 3.45100038854536, 0.509259838818042, 0.557854035598286, - 0.0525817201002309, -2.05064423770094, -1.11567547099143, - -1.18366549451658, 0.511958575232253, 0.334364244817592, - -2.05706205756846, -0.121896645718402, 0.54423804521097, - 0.600501641888935, 0.531074442421607, 0.495400347786053, - -1.60860176347294, 0.277370954937718, 0.290464363258084, - 0.579894254111128, -0.290441177117614, 1.30622601704777, - -0.482897816720494, -0.716423394441349, 0.742412540254878, - -0.91161346344296, 1.27943556865527, -0.189153217811851, - 0.592292730243945, 0.952415888511291, 0.491436176457309, - -0.30568088056691, -0.363871357644093, -0.285424162901343, - -0.0372094292657342, -0.923529247741133, 1.13805430719146, - -1.33122338081553, 0.55123448290722, -0.852726057929887, - 1.19687530878469, 0.498781686408254, 0.320179856418398, - 0.21244678210441, 1.00935803951191, -0.900989007058962, - 1.13216444413294, 0.0188670824356362, 0.424169515300288, - -0.19862121711326, 0.955170163539181, 0.948320512371124, - 0.473848149342783, -0.699121305560135, -0.612853026250685, - 0.580431200426044, 1.27799640925722, 0.806797458367235, - -1.03855925707748, 1.00866312622584, -0.578256568822387, - -0.323244575961333, -0.756301997657156, 1.38635211208482, - 0.722419488760045, -1.2160777034384, -0.498279906600592, - 0.726247405185, -0.260119271608589, -0.741134528045221, - -0.184110574491516, 0.307761674659839, 0.464568227698959, - -0.25253136951752, -0.486503680414154, 0.426634057655542, - -1.30396915580526, 0.0671486396913438, 1.77117635735777, - 0.907249468179712, 0.432349548721498, 1.41989705188111, - -0.413389471016361, 2.44202481656431, 0.0411377482323225, - 0.509505377681864, -0.282743502058313, 0.179881630718384, - -1.18808328118875, 0.98265314676344, -1.04288590077335, - 1.18136543233696, -0.398339818481707, -1.33556478800344, - -0.502789555455575, 0.484761653956289, -0.806445812279308, - 1.41207651978306, -0.878873945799123, -0.935197083131863, - -0.33925477332393439, 0.16449721487453731, 1.3700178698345999, - -1.4946841727166, 1.3805047732704381, - 0.88508389905048512, 0.83560940141892148, 0.89623509727336315, - -1.289541425794579, 0.2332028995229195, 1.183197522953588, - -0.85793361589157902, -1.33423445483833, -0.9233512315474407, - 0.76914556896670361, -0.37794794349382183, 0.059114341211622581, - -1.8706153553475069, -0.67786838062170507, 0.038184754648735768, - 0.37530087746353391, 0.96471695952212921, 0.69505105492152874, - -0.34214020737803602, -1.1454631827640021, -0.99324551114161375, - -0.13057284978088679, 1.213711380869505, 0.29124075688915307, - 1.106890512068581, 0.94957063346615733, 0.46367541051066768, - 0.45572327290248621, 0.39878553409592049, -0.015849431703916221, - -1.3973725035064171, 0.7700624622976332, 0.083291190129894818, - 0.53179773252409901, 0.049727349788233177, -0.73414037626738005, - -0.96348659055127073, 0.57356064323574374, -0.28194211032947131, - -0.59450289683584279, 0.77026173196827941, 1.0739830028467161, - -0.61570603602075391, -0.084794357704615464, -0.49163022652120109, - -1.526968705617602, -0.19688130817103111, 0.1656534453607213, - 0.19835657518696179, 0.97492977599052544, -0.95484796495550817, - 0.58847390467129868, -0.42688317000127768, 0.1771186872105201, - -0.91644209647809238, -1.8851386926119349, 0.086893856222760746, - 0.45630642515021741, 0.17428542070878469, -0.0013077214871275221, - -0.00058541929918550742, 0.28402285608099398, -0.36567881757010029, - -0.54886653165173238, 0.8578476816688223, 0.69909448655308448, - -0.14002628501260239, 1.332454137144605, 1.6017946938719501, - 0.01241549637061686, 0.24342918633361621, 1.0773689561938919, - 1.8592463357601141, 0.18590984985424869, 0.033342258305766252, - 0.6130082357970067, 1.068594886375418, -0.68330464261374424, - -0.12882583544682871, -1.6555248021907429, 0.013086014377651681, - 0.062454455755349927, 0.77304176654886514, 0.12704646649909671, - 0.40865153244567209, 1.195437623807228, -0.18555786800092541, - -1.299714084101439, 0.89967540292281434, -0.033647925669371137, - -1.5446015243088369, 0.65520298400478949, -0.71393501757996425}; -double response[] = {14.4535533874191, 6.62282520910778, 13.6714139876233, - 14.1975175236874, 12.8605301149348, 12.5228556826206, 14.2146384581959, - 7.9242642010286, 12.5069380013745, 13.7342047122325, 14.7108554131065, - 13.5962229304995, 5.89001909002711, 13.5586535685782, 14.0431671811957, - 13.9313910018427, 13.2189198447833, 17.0905598230825, 15.1993220372035, - 13.2616669404325, 15.7606359467964, 12.0838552528602, 14.344906985408, - 12.6094936116173, 11.9329594317628, 13.4086741328164, 13.7007653532941, - 13.0133656112894, 15.794998892751, 14.600198458049, 16.2757508936254, - 11.5643493993645, 14.8090225170414, 12.9823612913134, 15.003502495484, - 14.7373366435951, 15.7476765061616, 11.6745084114309, 14.047278212178, - 14.6669170934119, 13.8062403198314, 13.6111487435938, 13.3471486192318, - 14.2251519152709, 14.7188461068404, 14.2172164843947, 14.4180584862351, - 14.7196335400403, 12.799715984732, 13.9330377247579, 15.2646032349699, - 14.6603872891079, 9.73869078623634, 14.4434243169553, 14.4172837909381, - 15.1845379738711, 13.3449384473427, 15.3729427547467, 13.8115544407009, - 15.103777322749, 15.3838341258708, 14.368611819712, 12.525202176137, - 14.3250330647389, 15.2596577477861, 13.0045474727206, 14.515987797507, - 15.176981889542, 14.9241874861469, 13.872430121229, 15.3953655496863, - 13.4280761187509, 15.2034304536162, 14.1866308929129, 13.3058326261246, - 14.0746238485616, 14.1030921763152, 13.49966901054, 11.5846746059002, - 14.2648911116312, 14.88561614061, 13.9672969505607, 16.604679813678, - 10.3676055239145, 14.7434725924834, 16.3088265042892, 14.1086733681544, - 13.5909878288487, 14.6745463058857, 15.2940472804827, 14.6867226502357, - 13.6114224063955, 11.9702698734486, 13.8841573398, 15.0717757159234, - 12.5898155750775, 13.8187450898422, 14.2453171289186, 14.4065299197652, - 14.3479407847109}; -double newdata1[] = {-2.5, 0., 2.5, 0., 0., 0.}; -double newdata2[] = {-0.5, 0.5, 0., 0.}; -double coverage = .99; -long n = 100, p = 2, m = 3, se_fit = FALSE; -int i; - -main() { - printf("\nloess(&madeup):\n"); - loess_setup(one_two, response, n, p, &madeup); - madeup.model.span = 0.5; - loess(&madeup); - loess_summary(&madeup); - - printf("\nloess(&madeup_new):\n"); - loess_setup(one_two, response, n, p, &madeup_new); - madeup_new.model.span = 0.8; - madeup_new.model.drop_square[0] = TRUE; - madeup_new.model.parametric[0] = TRUE; - loess(&madeup_new); - loess_summary(&madeup_new); - - printf("\nloess(&madeup_new) (family = symmetric):\n"); - madeup_new.model.family = "symmetric"; - loess(&madeup_new); - loess_summary(&madeup_new); - - printf("\nloess(&madeup_new) (normalize = FALSE):\n"); - madeup_new.model.normalize = FALSE; - loess(&madeup_new); - loess_summary(&madeup_new); - - printf("\npredict(newdata1, m, &madeup, &madeup_pred, %d):\n", se_fit); - predict(newdata1, m, &madeup, &madeup_pred, se_fit); - printf("%g %g %g\n", madeup_pred.fit[0], - madeup_pred.fit[1], madeup_pred.fit[2]); - - m = 2; - se_fit = TRUE; - printf("\npredict(newdata2, m, &madeup, &madeup_pred, %d):\n", se_fit); - predict(newdata2, m, &madeup, &madeup_pred, se_fit); - printf("%g %g\n", madeup_pred.fit[0], madeup_pred.fit[1]); - printf("%g %g\n", madeup_pred.se_fit[0], madeup_pred.se_fit[1]); - printf("%g\n", madeup_pred.residual_scale); - printf("%g\n", madeup_pred.df); - - printf("\npointwise(&madeup_pred, m, coverage, &madeup_ci):\n"); - pointwise(&madeup_pred, m, coverage, &madeup_ci); - for(i = 0; i < m; i++) - printf("%g ", madeup_ci.upper[i]); - printf("\n"); - for(i = 0; i < m; i++) - printf("%g ", madeup_ci.fit[i]); - printf("\n"); - for(i = 0; i < m; i++) - printf("%g ", madeup_ci.lower[i]); - printf("\n"); - - loess_setup(one_two, response, n, p, &madeup2); - madeup2.model.span = 0.75; - loess(&madeup2); - - printf("\nanova(&madeup2, &madeup, &madeup_anova):\n"); - anova(&madeup2, &madeup, &madeup_anova); - printf("%g %g %g %g\n", madeup_anova.dfn, madeup_anova.dfd, - madeup_anova.F_value, madeup_anova.Pr_F); - - loess_free_mem(&madeup); - loess_free_mem(&madeup2); - loess_free_mem(&madeup_new); - pred_free_mem(&madeup_pred); - pw_free_mem(&madeup_ci); -} //GO.SYSIN DD madeup.c echo makefile 1>&2 sed >makefile <<'//GO.SYSIN DD makefile' 's/^-//' -CFLAGS = -cckr -O -FFLAGS = -O -OBJ = loessc.o loess.o predict.o misc.o loessf.o - -gas: gas.x - gas.x -gas.x: gas.o $(OBJ) - cc -o gas.x gas.o $(OBJ) -llinpack -lblas -lm -lF77 - -madeup: madeup.x - madeup.x -madeup.x: madeup.o $(OBJ) - cc -o madeup.x madeup.o $(OBJ) -llinpack -lblas -lm -lF77 - -ethanol: ethanol.x - ethanol.x -ethanol.x: ethanol.o $(OBJ) - cc -o ethanol.x ethanol.o $(OBJ) -llinpack -lblas -lm -lF77 - -air: air.x - air.x -air.x: air.o $(OBJ) - cc -o air.x air.o $(OBJ) -llinpack -lblas -lm -lF77 - -galaxy: galaxy.x - galaxy.x -galaxy.x: galaxy.o $(OBJ) - cc -o galaxy.x galaxy.o $(OBJ) -llinpack -lblas -lm -lF77 - -clean: - rm -f *.o *.x core //GO.SYSIN DD makefile echo misc.c 1>&2 sed >misc.c <<'//GO.SYSIN DD misc.c' 's/^-//' -#include "S.h" -#include "loess.h" - -/* If your compiler is so ancient it doesn't recognize void, say -#define void -*/ - -void -anova(one, two, out) -struct loess_struct *one, *two; -struct anova_struct *out; -{ - double one_d1, one_d2, one_s, two_d1, two_d2, two_s, - rssdiff, d1diff, tmp, pf(); - int max_enp; - - one_d1 = one->out.one_delta; - one_d2 = one->out.two_delta; - one_s = one->out.s; - two_d1 = two->out.one_delta; - two_d2 = two->out.two_delta; - two_s = two->out.s; - - rssdiff = fabs(one_s * one_s * one_d1 - two_s * two_s * two_d1); - d1diff = fabs(one_d1 - two_d1); - out->dfn = d1diff * d1diff / fabs(one_d2 - two_d2); - max_enp = (one->out.enp > two->out.enp); - tmp = max_enp ? one_d1 : two_d1; - out->dfd = tmp * tmp / (max_enp ? one_d2 : two_d2); - tmp = max_enp ? one_s : two_s; - out->F_value = (rssdiff / d1diff) / (tmp * tmp); - out->Pr_F = 1 - pf(out->F_value, out->dfn, out->dfd); -} - -void -pointwise(pre, m, coverage, ci) -struct pred_struct *pre; -long m; -double coverage; -struct ci_struct *ci; -{ - double t_dist, limit, fit, qt(); - int i; - - ci->fit = (double *) malloc(m * sizeof(double)); - ci->upper = (double *) malloc(m * sizeof(double)); - ci->lower = (double *) malloc(m * sizeof(double)); - - t_dist = qt(1 - (1 - coverage)/2, pre->df); - for(i = 0; i < m; i++) { - limit = pre->se_fit[i] * t_dist; - ci->fit[i] = fit = pre->fit[i]; - ci->upper[i] = fit + limit; - ci->lower[i] = fit - limit; - } -} - -void -pw_free_mem(ci) -struct ci_struct *ci; -{ - free(ci->fit); - free(ci->upper); - free(ci->lower); -} - -double -pf(q, df1, df2) -double q, df1, df2; -{ - double ibeta(); - - return(ibeta(q*df1/(df2+q*df1), df1/2, df2/2)); -} - -double -qt(p, df) -double p, df; -{ - double t, invibeta(); - - t = invibeta(fabs(2*p-1), 0.5, df/2); - return((p>0.5?1:-1) * sqrt(t*df/(1-t))); -} - -/**********************************************************************/ - /* - * Incomplete beta function. - * Reference: Abramowitz and Stegun, 26.5.8. - * Assumptions: 0 <= x <= 1; a,b > 0. - */ -#define DOUBLE_EPS 2.2204460492503131E-16 -#define IBETA_LARGE 1.0e30 -#define IBETA_SMALL 1.0e-30 - -double -ibeta(x, a, b) -double x, a, b; -{ - int flipped = 0, i, k, count; - double I, temp, pn[6], ak, bk, next, prev, factor, val; - - if (x <= 0) - return(0); - if (x >= 1) - return(1); - - /* use ibeta(x,a,b) = 1-ibeta(1-x,b,a) */ - if ((a+b+1)*x > (a+1)) { - flipped = 1; - temp = a; - a = b; - b = temp; - x = 1 - x; - } - - pn[0] = 0.0; - pn[2] = pn[3] = pn[1] = 1.0; - count = 1; - val = x/(1.0-x); - bk = 1.0; - next = 1.0; - do { - count++; - k = count/2; - prev = next; - if (count%2 == 0) - ak = -((a+k-1.0)*(b-k)*val)/ - ((a+2.0*k-2.0)*(a+2.0*k-1.0)); - else - ak = ((a+b+k-1.0)*k*val)/ - ((a+2.0*k)*(a+2.0*k-1.0)); - pn[4] = bk*pn[2] + ak*pn[0]; - pn[5] = bk*pn[3] + ak*pn[1]; - next = pn[4] / pn[5]; - for (i=0; i<=3; i++) - pn[i] = pn[i+2]; - if (fabs(pn[4]) >= IBETA_LARGE) - for (i=0; i<=3; i++) - pn[i] /= IBETA_LARGE; - if (fabs(pn[4]) <= IBETA_SMALL) - for (i=0; i<=3; i++) - pn[i] /= IBETA_SMALL; - } while (fabs(next-prev) > DOUBLE_EPS*prev); - factor = a*log(x) + (b-1)*log(1-x); - factor -= gamma(a+1) + gamma(b) - gamma(a+b); - I = exp(factor) * next; - return(flipped ? 1-I : I); -} - -/* - * Rational approximation to inverse Gaussian distribution. - * Absolute error is bounded by 4.5e-4. - * Reference: Abramowitz and Stegun, page 933. - * Assumption: 0 < p < 1. - */ - -static double num[] = { - 2.515517, - 0.802853, - 0.010328 -}; - -static double den[] = { - 1.000000, - 1.432788, - 0.189269, - 0.001308 -}; - -double -invigauss_quick(p) -double p; -{ - int lower; - double t, n, d, q; - - if(p == 0.5) - return(0); - lower = p < 0.5; - p = lower ? p : 1 - p; - t = sqrt(-2 * log(p)); - n = (num[2]*t + num[1])*t + num[0]; - d = ((den[3]*t + den[2])*t + den[1])*t + den[0]; - q = lower ? n/d - t : t - n/d; - return(q); -} - -/* - * Inverse incomplete beta function. - * Assumption: 0 <= p <= 1, a,b > 0. - */ - -double -invibeta(p, a, b) -double p, a, b; -{ - int i; - double ql, qr, qm, qdiff; - double pl, pr, pm, pdiff; - double invibeta_quick(), ibeta(); - -/* MEANINGFUL(qm);*/ - qm = 0; - if(p == 0 || p == 1) - return(p); - - /* initialize [ql,qr] containing the root */ - ql = qr = invibeta_quick(p, a, b); - pl = pr = ibeta(ql, a, b); - if(pl == p) - return(ql); - if(pl < p) - while(1) { - qr += 0.05; - if(qr >= 1) { - pr = qr = 1; - break; - } - pr = ibeta(qr, a, b); - if(pr == p) - return(pr); - if(pr > p) - break; - } - else - while(1) { - ql -= 0.05; - if(ql <= 0) { - pl = ql = 0; - break; - } - pl = ibeta(ql, a, b); - if(pl == p) - return(pl); - if(pl < p) - break; - } - - /* a few steps of bisection */ - for(i = 0; i < 5; i++) { - qm = (ql + qr) / 2; - pm = ibeta(qm, a, b); - qdiff = qr - ql; - pdiff = pm - p; - if(fabs(qdiff) < DOUBLE_EPS*qm || fabs(pdiff) < DOUBLE_EPS) - return(qm); - if(pdiff < 0) { - ql = qm; - pl = pm; - } else { - qr = qm; - pr = pm; - } - } - - /* a few steps of secant */ - for(i = 0; i < 40; i++) { - qm = ql + (p-pl)*(qr-ql)/(pr-pl); - pm = ibeta(qm, a, b); - qdiff = qr - ql; - pdiff = pm - p; - if(fabs(qdiff) < 2*DOUBLE_EPS*qm || fabs(pdiff) < 2*DOUBLE_EPS) - return(qm); - if(pdiff < 0) { - ql = qm; - pl = pm; - } else { - qr = qm; - pr = pm; - } - } - - /* no convergence */ - return(qm); -} - -/* - * Quick approximation to inverse incomplete beta function, - * by matching first two moments with the Gaussian distribution. - * Assumption: 0 < p < 1, a,b > 0. - */ - -double -invibeta_quick(p, a, b) -double p, a, b; -{ - double x, m, s, fmax(), fmin(), invigauss_quick(); - - x = a + b; - m = a / x; - s = sqrt((a*b) / (x*x*(x+1))); - return(fmax(0.0, fmin(1.0, invigauss_quick(p)*s + m))); -} - -static double -fmin(a, b) -double a, b; -{ - return(a < b ? a : b); -} - -static double -fmax(a, b) -double a, b; -{ - return(a > b ? a : b); -} - -long -max(a, b) -long a, b; -{ - return(a > b ? a : b); -} - -typedef double doublereal; -typedef int integer; - -void -Recover(a, b) -char *a; -int *b; -{ - printf(a); - exit(1); -} - -void -Warning(a, b) -char *a; -int *b; -{ - printf(a); -} - -/* d1mach may be replaced by Fortran code: - mail netlib@netlib.bell-labs.com - send d1mach from core. -*/ - -#include - -doublereal F77_SUB(d1mach) (i) -integer *i; -{ - switch(*i){ - case 1: return DBL_MIN; - case 2: return DBL_MAX; - case 3: return DBL_EPSILON/FLT_RADIX; - case 4: return DBL_EPSILON; - case 5: return log10(FLT_RADIX); - default: Recover("Invalid argument to d1mach()", 0L); - } -} - //GO.SYSIN DD misc.c echo predict.c 1>&2 sed >predict.c <<'//GO.SYSIN DD predict.c' 's/^-//' -#include "S.h" -#include "loess.h" - -void -predict(eval, m, lo, pre, se) -double *eval; -long m, se; -struct loess_struct *lo; -struct pred_struct *pre; -{ - long size_info[3]; - void pred_(); - - pre->fit = (double *) malloc(m * sizeof(double)); - pre->se_fit = (double *) malloc(m * sizeof(double)); - pre->residual_scale = lo->out.s; - pre->df = (lo->out.one_delta * lo->out.one_delta) / lo->out.two_delta; - - size_info[0] = lo->in.p; - size_info[1] = lo->in.n; - size_info[2] = m; - - pred_(lo->in.y, lo->in.x, eval, size_info, &lo->out.s, - lo->in.weights, - lo->out.robust, - &lo->model.span, - &lo->model.degree, - &lo->model.normalize, - lo->model.parametric, - lo->model.drop_square, - &lo->control.surface, - &lo->control.cell, - &lo->model.family, - lo->kd_tree.parameter, - lo->kd_tree.a, - lo->kd_tree.xi, - lo->kd_tree.vert, - lo->kd_tree.vval, - lo->out.divisor, - &se, - pre->fit, - pre->se_fit); -} - -void -pred_(y, x_, new_x, size_info, s, weights, robust, span, degree, - normalize, parametric, drop_square, surface, cell, family, - parameter, a, xi, vert, vval, divisor, se, fit, se_fit) -double *y, *x_, *new_x, *weights, *robust, *span, *cell, *fit, *s, - *xi, *vert, *vval, *divisor, *se_fit; -long *size_info, *degree, *normalize, *parametric, *drop_square, - *parameter, *a, *se; -char **surface, **family; -{ - double *x, *x_tmp, *x_evaluate, *L, new_cell, z, tmp, *fit_tmp, - *temp, sum, mean; - long N, D, M, sum_drop_sqr = 0, sum_parametric = 0, - nonparametric = 0, *order_parametric, *order_drop_sqr; - int i, j, k, p, cut, comp(); - - D = size_info[0]; - N = size_info[1]; - M = size_info[2]; - - x = (double *) malloc(N * D * sizeof(double)); - x_tmp = (double *) malloc(N * D * sizeof(double)); - x_evaluate = (double *) malloc(M * D * sizeof(double)); - L = (double *) malloc(N * M * sizeof(double)); - order_parametric = (long *) malloc(D * sizeof(long)); - order_drop_sqr = (long *) malloc(D * sizeof(long)); - temp = (double *) malloc(N * D * sizeof(double)); - - for(i = 0; i < (N * D); i++) - x_tmp[i] = x_[i]; - for(i = 0; i < D; i++) { - k = i * M; - for(j = 0; j < M; j++) { - p = k + j; - new_x[p] = new_x[p] / divisor[i]; - } - } - if(!strcmp(*surface, "direct") || se) { - for(i = 0; i < D; i++) { - k = i * N; - for(j = 0; j < N; j++) { - p = k + j; - x_tmp[p] = x_[p] / divisor[i]; - } - } - } - j = D - 1; - for(i = 0; i < D; i++) { - sum_drop_sqr = sum_drop_sqr + drop_square[i]; - sum_parametric = sum_parametric + parametric[i]; - if(parametric[i]) - order_parametric[j--] = i; - else - order_parametric[nonparametric++] = i; - } - for(i = 0; i < D; i++) { - order_drop_sqr[i] = 2 - drop_square[order_parametric[i]]; - k = i * M; - p = order_parametric[i] * M; - for(j = 0; j < M; j++) - x_evaluate[k + j] = new_x[p + j]; - k = i * N; - p = order_parametric[i] * N; - for(j = 0; j < N; j++) - x[k + j] = x_tmp[p + j]; - } - for(i = 0; i < N; i++) - robust[i] = weights[i] * robust[i]; - - if(!strcmp(*surface, "direct")) { - if(*se) { - loess_dfitse(y, x, x_evaluate, weights, robust, - !strcmp(*family, "gaussian"), span, degree, - &nonparametric, order_drop_sqr, &sum_drop_sqr, - &D, &N, &M, fit, L); - } - else { - loess_dfit(y, x, x_evaluate, robust, span, degree, - &nonparametric, order_drop_sqr, &sum_drop_sqr, - &D, &N, &M, fit); - } - } - else { - loess_ifit(parameter, a, xi, vert, vval, &M, x_evaluate, fit); - if(*se) { - new_cell = (*span) * (*cell); - fit_tmp = (double *) malloc(M * sizeof(double)); - loess_ise(y, x, x_evaluate, weights, span, degree, - &nonparametric, order_drop_sqr, &sum_drop_sqr, - &new_cell, &D, &N, &M, fit_tmp, L); - free(fit_tmp); - } - } - if(*se) { - for(i = 0; i < N; i++) { - k = i * M; - for(j = 0; j < M; j++) { - p = k + j; - L[p] = L[p] / weights[i]; - L[p] = L[p] * L[p]; - } - } - for(i = 0; i < M; i++) { - tmp = 0; - for(j = 0; j < N; j++) - tmp = tmp + L[i + j * M]; - se_fit[i] = (*s) * sqrt(tmp); - } - } - free(x); - free(x_tmp); - free(x_evaluate); - free(L); - free(order_parametric); - free(order_drop_sqr); - free(temp); -} - -void -pred_free_mem(pre) -struct pred_struct *pre; -{ - free(pre->fit); - free(pre->se_fit); -} - - - - - - //GO.SYSIN DD predict.c echo predict.m 1>&2 sed >predict.m <<'//GO.SYSIN DD predict.m' 's/^-//' -NAME - predict, pointwise, pred_free_mem, pw_free_mem - -SYNOPSIS - #include "loess.h" - - double *eval, coverage; - long m, se; - struct loess_struct *lo; - struct predict_struct *pre; - struct ci_struct *ci; - - void predict(eval, m, lo, pre, se) - - void pointwise(pre, m, coverage, ci) - - void pred_free_mem(pre) - - void pw_free_mem(ci) - -PARAMETERS - - eval a vector of length m * p specifying the values of the - predictors at which the evaluation is to be carried out. - The j-th coordinate of the i-th point is in eval[i+m*j], - where 0<=j&2 sed >struct.m <<'//GO.SYSIN DD struct.m' 's/^-//' -(All default values mentioned here are set by loess_setup().) - -struct loess_struct *lo; - -in - n: number of observations. - - p: number of numeric predictors. - - y: vector of response (length n). - - x: vector of predictors, of length (n * p). - The j-th coordinate of the i-th point is in x[i+n*j], - where 0<=j&2 sed >supp.f <<'//GO.SYSIN DD supp.f' 's/^-//' - subroutine ehg182(i) - integer i - if(i.eq.100) print *,'wrong version number in lowesd. Probably ty - +po in caller.' - if(i.eq.101) print *,'d>dMAX in ehg131. Need to recompile with in - +creased dimensions.' - if(i.eq.102) print *,'liv too small. (Discovered by lowesd)' - if(i.eq.103) print *,'lv too small. (Discovered by lowesd)' - if(i.eq.104) print *,'alpha too small. fewer data values than deg - +rees of freedom.' - if(i.eq.105) print *,'k>d2MAX in ehg136. Need to recompile with i - +ncreased dimensions.' - if(i.eq.106) print *,'lwork too small' - if(i.eq.107) print *,'invalid value for kernel' - if(i.eq.108) print *,'invalid value for ideg' - if(i.eq.109) print *,'lowstt only applies when kernel=1.' - if(i.eq.110) print *,'not enough extra workspace for robustness ca - +lculation' - if(i.eq.120) print *,'zero-width neighborhood. make alpha bigger' - if(i.eq.121) print *,'all data on boundary of neighborhood. make a - +lpha bigger' - if(i.eq.122) print *,'extrapolation not allowed with blending' - if(i.eq.123) print *,'ihat=1 (diag L) in l2fit only makes sense if - + z=x (eval=data).' - if(i.eq.171) print *,'lowesd must be called first.' - if(i.eq.172) print *,'lowesf must not come between lowesb and lowe - +se, lowesr, or lowesl.' - if(i.eq.173) print *,'lowesb must come before lowese, lowesr, or l - +owesl.' - if(i.eq.174) print *,'lowesb need not be called twice.' - if(i.eq.180) print *,'nv>nvmax in cpvert.' - if(i.eq.181) print *,'nt>20 in eval.' - if(i.eq.182) print *,'svddc failed in l2fit.' - if(i.eq.183) print *,'didnt find edge in vleaf.' - if(i.eq.184) print *,'zero-width cell found in vleaf.' - if(i.eq.185) print *,'trouble descending to leaf in vleaf.' - if(i.eq.186) print *,'insufficient workspace for lowesf.' - if(i.eq.187) print *,'insufficient stack space' - if(i.eq.188) print *,'lv too small for computing explicit L' - if(i.eq.191) print *,'computed trace L was negative; something is - +wrong!' - if(i.eq.192) print *,'computed delta was negative; something is wr - +ong!' - if(i.eq.193) print *,'workspace in loread appears to be corrupted' - if(i.eq.194) print *,'trouble in l2fit/l2tr' - if(i.eq.195) print *,'only constant, linear, or quadratic local mo - +dels allowed' - if(i.eq.196) print *,'degree must be at least 1 for vertex influen - +ce matrix' - if(i.eq.999) print *,'not yet implemented' - print *,'Assert failed, error code ',i - stop - end - subroutine ehg183(s,i,n,inc) - character*(*) s - integer n, inc, i(inc,n), j - print *,s,(i(1,j),j=1,n) - end - subroutine ehg184(s,x,n,inc) - character*(*) s - integer n, inc, j - double precision x(inc,n) - print *,s,(x(1,j),j=1,n) - end - subroutine losave(iunit,iv,liv,lv,v) - integer execnt,iunit,liv,lv - integer iv(liv) - DOUBLE PRECISION v(lv) - external ehg167 - save execnt - data execnt /0/ - execnt=execnt+1 - call ehg167(iunit,iv(2),iv(4),iv(5),iv(6),iv(14),v(iv(11)),iv(iv(7 - +)),v(iv(12)),v(iv(13))) - return - end - subroutine ehg167(iunit,d,vc,nc,nv,nvmax,v,a,xi,vval) - integer iunit,d,vc,nc,nv,a(nc),magic,i,j - DOUBLE PRECISION v(nvmax,d),xi(nc),vval(0:d,nv) - write(iunit,*)d,nc,nv - do 10 i=1,d -10 write(iunit,*)v(1,i),v(vc,i) - j = 0 - do 20 i=1,nc - if(a(i).ne.0)then - write(iunit,*)a(i),xi(i) - else - write(iunit,*)a(i),j - end if -20 continue - do 30 i=1,nv -30 write(iunit,*)(vval(j,i),j=0,d) - end - subroutine lohead(iunit,d,vc,nc,nv) - integer iunit,d,vc,nc,nv - read(iunit,*)d,nc,nv - vc = 2**d - end - subroutine loread(iunit,d,vc,nc,nv,iv,liv,lv,v) - integer bound,d,execnt,iunit,liv,lv,nc,nv,vc - integer iv(liv) - DOUBLE PRECISION v(lv) - external ehg168,ehg169,ehg182 - save execnt - data execnt /0/ - execnt=execnt+1 - iv(28)=173 - iv(2)=d - iv(4)=vc - iv(14)=nv - iv(17)=nc - iv(7)=50 - iv(8)=iv(7)+nc - iv(9)=iv(8)+vc*nc - iv(10)=iv(9)+nc - bound=iv(10)+nc - if(.not.(bound-1.le.liv))then - call ehg182(102) - end if - iv(11)=50 - iv(13)=iv(11)+nv*d - iv(12)=iv(13)+(d+1)*nv - bound=iv(12)+nc - if(.not.(bound-1.le.lv))then - call ehg182(103) - end if - call ehg168(iunit,d,vc,nc,nv,nv,v(iv(11)),iv(iv(7)),v(iv(12)),v(iv - +(13))) - call ehg169(d,vc,nc,nc,nv,nv,v(iv(11)),iv(iv(7)),v(iv(12)),iv(iv(8 - +)),iv(iv(9)),iv(iv(10))) - return - end - subroutine ehg168(iunit,d,vc,nc,nv,nvmax,v,a,xi,vval) - integer iunit,d,vc,nc,nv,a(nc),magic,i,j - DOUBLE PRECISION v(nvmax,d),xi(nc),vval(0:d,nv) - do 10 i=1,d -10 read(iunit,*)v(1,i),v(vc,i) - do 20 i=1,nc -20 read(iunit,*)a(i),xi(i) - do 30 i=1,nv -30 read(iunit,*)(vval(j,i),j=0,d) - end - subroutine ehg170(k,d,vc,nv,nvmax,nc,ncmax,a,c,hi,lo,v,vval,xi) - integer d,execnt,i,j,nc,ncmax,nv,nvmax,vc - integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax) - double precision v(nvmax,d),vval(0:d,nvmax),xi(ncmax) - save execnt - data execnt /0/ - execnt=execnt+1 - write(k,*)' double precision function loeval(z)' - write(k,50)d - write(k,*)' integer d,vc,nv,nc' - write(k,51)nc,vc,nc - write(k,52)nc,nc - write(k,53)nv,d - write(k,54)d,nv - write(k,55)nc - write(k,56) - write(k,57)d,vc,nv,nc -50 format(' double precision z(',i2,')') -51 format(' integer a(',i5,'), c(',i3,',',i5,')') -52 format(' integer hi(',i5,'), lo(',i5,')') -53 format(' double precision v(',i5,',',i2,')') -54 format(' double precision vval(0:',i2,',',i5,')') -55 format(' double precision xi(',i5,')') -56 format(' double precision ehg128') -57 format(' data d,vc,nv,nc /',i2,',',i3,',',i5,',',i5,'/') - do 3 i=1,nc - write(k,58)i,a(i) -58 format(' data a(',i5,') /',i5,'/') - if(a(i).ne.0)then - write(k,59)i,i,i,hi(i),lo(i),xi(i) -59 format(' data hi(',i5,'),lo(',i5,'),xi(',i5,') /', - $ i5,',',i5,',',1pe15.6,'/') - end if - do 4 j=1,vc - write(k,60)j,i,c(j,i) -60 format(' data c(',i3,',',i5,') /',i5,'/') - 4 continue - 3 continue - do 5 i=1,nv - write(k,61)i,vval(0,i) -61 format(' data vval(0,',i5,') /',1pe15.6,'/') - do 6 j=1,d - write(k,62)i,j,v(i,j) -62 format(' data v(',i5,',',i2,') /',1pe15.6,'/') - write(k,63)j,i,vval(j,i) -63 format(' data vval(',i2,',',i5,') /',1pe15.6,'/') - 6 continue - 5 continue - write(k,*)' loeval=ehg128(z,d,nc,vc,a,xi,lo,hi,c,v,nv,vval)' - write(k,*)' end' - return - end - subroutine lofort(iunit,iv,liv,lv,wv) - integer execnt,iunit - integer iv(*) - DOUBLE PRECISION wv(*) - external ehg170 - save execnt - data execnt /0/ - execnt=execnt+1 - call ehg170(iunit,iv(2),iv(4),iv(6),iv(14),iv(5),iv(17),iv(iv(7)), - +iv(iv(8)),iv(iv(9)),iv(iv(10)),wv(iv(11)),wv(iv(13)),wv(iv(12))) - return - end //GO.SYSIN DD supp.f