C C Copyright (C) 1998--2020 The R Core Team C The authors of this software are Cleveland, Grosse, and Shyu. C Copyright (c) 1989, 1992 by AT&T. C Permission to use, copy, modify, and distribute this software for any C purpose without fee is hereby granted, provided that this entire notice C is included in all copies of any software which is or includes a copy C or modification of this software and in all copies of the supporting C documentation for such software. C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED C WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C altered by B.D. Ripley to C C remove unused variables C make phi in ehg139 double precision to match calling sequence C pass integer not logical from C C C by M. Maechler, renaming ehg182() to loesswarn() C Note that loesswarn(errormsg_code) is in ./loessc.c subroutine ehg126(d,n,vc,x,v,nvmax) integer d,execnt,i,j,k,n,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 c initialize d1mach(2) === DBL_MAX: 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) c Integer division would do here j=INT(DBLE(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,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 loesswarn 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 loesswarn(180) end if return end integer function ehg138(i,z,a,xi,lo,hi,ncmax) logical i1 integer i,j,ncmax integer a(ncmax),hi(ncmax),lo(ncmax) DOUBLE PRECISION xi(ncmax),z(8) 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) c Partial sorting of p(1, il:ir) returning the sort indices pi() only c such that p(1, pi(k)) is correct c implicit none c Arguments c Input: integer il,ir,k,nk,n DOUBLE PRECISION p(nk,n) c using only p(1, pi(*)) c Output: integer pi(n) c Variables DOUBLE PRECISION t integer i,ii,j,l,r c find the $k$-th smallest of $n$ elements c Floyd+Rivest, CACM Mar '75, Algorithm 489 l=il r=ir c while (l < r ) 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, + rcond,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) integer idamax double precision d1mach, ddot external ehg106,loesswarn,ehg184,dqrdc,dqrsl,dsvdc external idamax, d1mach, 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 c initialize d1mach(4) === 1 / DBL_EPSILON === 2^52 : 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(rho .le. 0)then call loesswarn(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(1),dd,1) call ehg184('radius ',rho,1,1) if(.not..false.)then call loesswarn(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 c FIXME: this has i = 3 vs bound 2 in a ggplot2 test 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 loesswarn(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('pseudoinverse used at',q(1),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, + nvmax,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,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(ncmax),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), + eta(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,loesswarn,ehg139,ehg124 double precision dnrm2 external dnrm2 c Identity -> identi c X -> b if(.not.(d.le.8))then call loesswarn(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, + dist,eta,b,d,w,diagl,vval2,nc,vc,a,xi,lo,hi,c,rcond, + sing,dd,tdeg,cdeg,lq,lf,setlf,vval) return end c called from lowese() only : subroutine ehg133(d,vc,nvmax,ncmax, a,c,hi,lo, v,vval,xi,m,z,s) integer d,vc,nvmax,ncmax, m integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax) double precision v(nvmax,d),vval(0:d,nvmax),xi(ncmax),z(m,d),s(m) c Var double precision delta(8) integer i,i1 double precision ehg128 external ehg128 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 i,j integer iw(i) iw(i)=j return end subroutine ehg141(trl,n,deg,k,d,nsing,dk,delta1,delta2) double precision trl,delta1,delta2 integer n,deg,k,d,nsing,dk double precision c(48), c1, c2, c3, c4, corx,z,zz(1) integer i external ehg176 double precision ehg176 c coef, d, deg, del data c / .2971620d0,.3802660d0,.5886043d0,.4263766d0,.3346498d0, +.6271053d0,.5241198d0,.3484836d0,.6687687d0,.6338795d0,.4076457d0, +.7207693d0,.1611761d0,.3091323d0,.4401023d0,.2939609d0,.3580278d0, +.5555741d0,.3972390d0,.4171278d0,.6293196d0,.4675173d0,.4699070d0, +.6674802d0,.2848308d0,.2254512d0,.2914126d0,.5393624d0,.2517230d0, +.3898970d0,.7603231d0,.2969113d0,.4740130d0,.9664956d0,.3629838d0, +.5348889d0,.2075670d0,.2822574d0,.2369957d0,.3911566d0,.2981154d0, +.3623232d0,.5508869d0,.3501989d0,.4371032d0,.7002667d0,.4291632d0, +.4930370d0 / if(deg.eq.0) dk=1 if(deg.eq.1) dk=d+1 if(deg.eq.2) dk=int(dble((d+2)*(d+1))/2.d0) corx=dsqrt(k/dble(n)) z=(dsqrt(k/trl)-corx)/(1-corx) if(nsing .eq. 0 .and. 1 .lt. z) then call ehg184('Chernobyl! trLn',trl,1,1) endif z=min(1.0d0,max(0.0d0,z)) c R fix zz(1)=z c4=dexp(ehg176(zz)) 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 i,j,n double precision delta1,delta2,trl double precision l(n,n),ll(n,n) double precision ddot external ddot 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,vc,nc,ncmax,nv,nvmax integer a(ncmax), c(vc,ncmax), hi(ncmax), lo(ncmax) DOUBLE PRECISION v(nvmax,d),xi(ncmax) integer novhit(1),i,j,k,mc,mv,p external ehg125,loesswarn integer ifloor external ifloor 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(DBLE(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 loesswarn(193) end if if(.not.(mv.eq.nv))then call loesswarn(193) end if return end DOUBLE PRECISION function ehg176(z) c DOUBLE PRECISION z(*) c 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 external 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 n,d,tau,nsing double precision trl, delta1,delta2 integer dka,dkb double precision alpha,d1a,d1b,d2a,d2b external ehg141 call ehg141(trl,n,1,tau,d,nsing,dka,d1a,d2a) call ehg141(trl,n,2,tau,d,nsing,dkb,d1b,d2b) alpha=dble(tau-dka)/dble(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, + vval2,lf,lq) c Args integer m,d,n,nf,nv,ncmax,nvmax,vc double precision z(m,d), l(m,n), xi(ncmax), v(nvmax,d), + vval2(0:d,nvmax), lf(0:d,nvmax,nf) integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),lo(ncmax),hi(ncmax) c Var integer lq1,i,i1,i2,j,p double precision zi(8) double precision ehg128 external ehg128 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,tau double precision alpha,f,trl,trla,trlb external ehg197 call ehg197(1,d,f,dka,trla) call ehg197(2,d,f,dkb,trlb) alpha=dble(tau-dka)/dble(dkb-dka) trl=(1-alpha)*trla+alpha*trlb return end subroutine ehg197(deg,d,f,dk,trl) integer deg,d,dk double precision f, trl double precision g1 dk = 0 if(deg.eq.1) dk=d+1 if(deg.eq.2) dk=int(dble((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,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) 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) c implicit none c Args integer d,ncmax,nvmax,vc integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax) DOUBLE PRECISION z(d),xi(ncmax),v(nvmax,d), vval(0:d,nvmax) c Vars logical i2,i3,i4,i5,i6,i7,i8,i9,i10 integer i,i1,i11,i12,ig,ii,j,lg,ll,m,nt,ur integer 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) external loesswarn,ehg184 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 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 loesswarn(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(1),d,1) call ehg184('lowerlimit ',v(ll,1),d,nvmax) else if(1.001D0.lt.h)then call ehg184('eval ',z(1),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 loesswarn(122) end if lg=int(DBLE(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=int(x) if(ifloor.gt.x) ifloor=ifloor-1 end c DSIGN is unused, causes conflicts on some platforms c DOUBLE PRECISION function DSIGN(a1,a2) c DOUBLE PRECISION a1, a2 c DSIGN=DABS(a1) c if(a2.ge.0)DSIGN=-DSIGN c end c ehg136() is the workhorse of lowesf(.) c n = number of observations c m = number of x values at which to evaluate c f = span c nf = min(n, floor(f * n)) subroutine ehg136(u,lm,m,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b, + od,o,ihat,w,rcond,sing,dd,tdeg,cdeg,s) integer identi,d,dd,i,i1,ihat,info,j,k,kernel,l,lm,m,n,nf, + od,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,loesswarn,dqrsl double precision ddot external ddot c V -> g c U -> e c Identity -> identi c L -> o c X -> b if(k .gt. nf-1) call loesswarn(104) if(k .gt. 15) call loesswarn(105) 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, + rcond,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 loesswarn(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 loesswarn(123) c next line is not in current dloess goto 7 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 c called from lowesb() ... compute fit ..?..?... c somewhat similar to ehg136 subroutine ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k, + dist,phi,eta,b,od,w,diagl,vval2,ncmax,vc,a,xi,lo,hi,c, + rcond,sing,dd,tdeg,cdeg,lq,lf,setlf,s) logical setlf integer identi,d,dd,i,i2,i3,i5,i6,ii,ileaf,info,j,k,kernel, + l,n,ncmax,nf,nleaf,nv,nvmax,od,sing,tdeg,vc integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax), + leaf(256),lo(ncmax),pi(n),psi(n) 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:od,nv),v(nvmax,d),vval2(0:d,nv),w(nf),work(15), + x(n,d),xi(ncmax),y(n),z(8) DOUBLE PRECISION phi(n) external ehg127,loesswarn,DQRSL,ehg137 DOUBLE PRECISION ehg128 external ehg128 DOUBLE PRECISION DDOT external DDOT c V -> e c Identity -> identi c X -> b c l2fit with trace(L) if(k .gt. nf-1) call loesswarn(104) if(k .gt. 15) call loesswarn(105) 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, + rcond,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,leaf,nleaf,d,ncmax,a,xi, + lo,hi) do 12 ileaf=1,nleaf do 13 ii=lo(leaf(ileaf)),hi(leaf(ileaf)) i=int(phi(pi(ii))) if(i.ne.0)then if(.not.(psi(i).eq.pi(ii)))then call loesswarn(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, + work,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.0d0 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 loesswarn(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, + 10000,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 lowesb(xx,yy,ww,diagl,infl,iv,wv) integer infl integer iv(*) DOUBLE PRECISION xx(*),yy(*),ww(*),diagl(*),wv(*) c Var DOUBLE PRECISION trl logical setlf integer ifloor external ifloor external ehg131,loesswarn,ehg183 if(.not.(iv(28).ne.173))then call loesswarn(174) end if if(iv(28).ne.172)then if(.not.(iv(28).eq.171))then call loesswarn(171) end if end if iv(28)=173 if(infl.ne.0)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)+DBLE(iv(4))/2.D0)then call ehg183('k-d tree limited by memory; nvmax=', + iv(14),1,1) else if(iv(17).lt.iv(5)+2)then call ehg183('k-d tree limited by memory. ncmax=', + iv(17),1,1) end if end if return end c lowesd() : Initialize iv(*) and v(1:4) c ------ called only by loess_workspace() in ./loessc.c subroutine lowesd(iv, liv,lv, v, d,n,f,ideg,nf,nvmax, setlf) integer liv,lv, d,n, ideg,nf,nvmax, setlf c setlf {Rboolean}: if(true) need L [nf x nvmax] matrices integer iv(liv) double precision f, v(lv) c had nf = min(n,ifloor(n*f)) integer bound,i,i1,i2,j,ncmax,vc external loesswarn integer ifloor external ifloor c c unnecessary initialization of i1 to keep g77 -Wall happy c i1 = 0 cc version -> versio c if(.not.(versio.eq.106))then c call loesswarn(100) c end if iv(28)=171 iv(2)=d iv(3)=n vc=2**d iv(4)=vc if(.not.(0.lt.f))then call loesswarn(120) end if 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=int(dble((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 loesswarn(195) end if if(.not.(ideg.le.2))then call loesswarn(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.ne.0)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 loesswarn(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.ne.0)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 loesswarn(103) end if v(1)=f v(2)=0.05d0 v(3)=0.d0 v(4)=1.d0 return end subroutine lowese(iv,wv,m,z,s) integer m integer iv(*) double precision s(m),wv(*),z(m,1) external ehg133,loesswarn if(.not.(iv(28).ne.172))then call loesswarn(172) end if if(.not.(iv(28).eq.173))then call loesswarn(173) end if call ehg133(iv(2),iv(4),iv(14),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 c "direct" (non-"interpolate") fit aka predict() : subroutine lowesf(xx,yy,ww,iv,wv,m,z,l,ihat,s) integer m,ihat c m = number of x values at which to evaluate integer iv(*) double precision xx(*),yy(*),ww(*),wv(*),z(m,1),l(m,*),s(m) logical i1 external loesswarn,ehg136 if(171.le.iv(28))then i1=(iv(28).le.174) else i1=.false. end if if(.not.i1)then call loesswarn(171) end if iv(28)=172 if(.not.(iv(14).ge.iv(19)))then call loesswarn(186) end if c do the work; in ehg136() give the argument names as they are there: c ehg136(u,lm,m, n, d, nf, f, x, psi, y ,rw, call ehg136(z,m,m,iv(3),iv(2),iv(19),wv(1),xx,iv(iv(22)),yy,ww, c kernel, k, dist, eta, b, od,o,ihat, + iv(20),iv(29),wv(iv(15)),wv(iv(16)),wv(iv(18)),0,l,ihat, c w, rcond,sing, dd, tdeg,cdeg, s) + wv(iv(26)),wv(4),iv(30),iv(33),iv(32),iv(41),s) return end c Called either from loess_raw() only for case [surf_stat = "interpolate/exact"], or c from loess_ise() {used only when 'se = TRUE' and surface = "interpolate"} subroutine lowesl(iv,wv,m,z,l) integer m integer iv(*) double precision l(m,*), wv(*), z(m,1) external loesswarn,ehg191 if(.not.(iv(28).ne.172))then call loesswarn(172) end if if(.not.(iv(28).eq.173))then call loesswarn(173) end if if(.not.(iv(26).ne.iv(34)))then call loesswarn(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 c Not used subroutine lowesr(yy,iv,wv) integer iv(*) DOUBLE PRECISION yy(*),wv(*) external loesswarn,ehg192 if(.not.(iv(28).ne.172))then call loesswarn(172) end if if(.not.(iv(28).eq.173))then call loesswarn(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 c Update Robustness Weights -- called via .Fortran() from R's simpleLoess() in ../R/loess.R c subroutine lowesw(res,n,rw,pi) c Tranliterated from Devlin's ratfor c implicit none c Args integer n double precision res(n), rw(n) integer pi(n) c Var integer identi,i,i1,nh double precision cmad,rsmall integer ifloor double precision d1mach external ehg106 external ifloor external d1mach c Identity -> identi 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(dble(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 c Compute Pseudovalues -- called via .Fortran() from R's simpleLoess() in ../R/loess.R c in the case of robustness iterations subroutine lowesp(n,y,yhat,pwgts,rwgts,pi,ytilde) integer n integer pi(n) double precision y(n),yhat(n),pwgts(n),rwgts(n),ytilde(n) c Var double precision c,i1,i4,mad integer i2,i3,i,m external ehg106 integer ifloor external ifloor c median absolute deviation (using partial sort): do 3 i=1,n ytilde(i)=dabs(y(i)-yhat(i))*dsqrt(pwgts(i)) pi(i) = i 3 continue m=ifloor(dble(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 i=1,n ytilde(i)= 1 - ((y(i)-yhat(i))**2 * pwgts(i))/c 5 continue do 6 i=1,n ytilde(i)=ytilde(i)*dsqrt(rwgts(i)) 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 i=1,n ytilde(i)=yhat(i) + (c*rwgts(i))*(y(i)-yhat(i)) 8 continue return end subroutine ehg124(ll,uu,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v, + vhit,nvmax,fc,fd,dd) integer ll,uu,d,n,nv,nc,ncmax,vc,nvmax,fc,dd integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),pi(n),vhit(nvmax) DOUBLE PRECISION fd, v(nvmax,d),x(n,d),xi(ncmax) logical i1,i2,leaf integer i4,inorm2,k,l,m,p,u, upper, lower, check, offset DOUBLE PRECISION diam,diag(8),sigma(8) external ehg125,ehg106,ehg129 integer IDAMAX external IDAMAX 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+DBLE(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=int(DBLE(l+u)/2.D0) call ehg106(l,u,m,1,x(1,k),pi,n) c all ties go with hi son c top of while loop c bug fix from btyner@gmail.com 2006-07-20 offset = 0 7 if(((m+offset).ge.u).or.((m+offset).lt.l))goto 8 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 c bottom of while loop 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 c initialize d1mach(2) === DBL_MAX: 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 c {called only from ehg127} purpose...?... subroutine ehg137(z,leaf,nleaf,d,ncmax,a,xi,lo,hi) integer d,nleaf,ncmax integer leaf(256),a(ncmax),hi(ncmax),lo(ncmax),pstack(20) DOUBLE PRECISION z(d),xi(ncmax) integer p,stackt external loesswarn c stacktop -> stackt 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 loesswarn(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 loesswarn(185) end if return end C-- For Error messaging, call the "a" routines at the bottom of ./loessc.c : subroutine ehg183(s, i, n, inc) character s*(*) integer i, n, inc call ehg183a(s, len(s), i, n, inc) end subroutine ehg184(s, x, n, inc) character s*(*) double precision x integer n, inc call ehg184a(s, len(s), x, n, inc) end