## [R] translation of lowess_module.f90 
## Dan Gavin; Mar 1, 2008
## Fortran code is commented out

## module lowess
## contains

## subroutine lowess_curve(nctrl,xin,yin,ntarg,xt,wintype="w",order=0, nrobust=0,width=300,yfit)
## manages computation of a lowess "curve"

lowess_curve<-function(xin,yin,xt,width=300,wintype="w",order=0,nrobust=0){

    ## xin=array of x values
    ## yin=array of y values
    ## xt=array of target x values to fit loess.  Use seq(low,high,step).
    ## wintype= "s" to use obs in width proportion of data
    ## wintype= "w" to use obs within plus or minus window half-width (width)
    ## order=0,1,2 (median, linear, quadratic)
    ## nrobust=number of robustifying iterations
    
    ## implicit none
    ## integer             :: nctrl,ntarg
    ## real(8)             :: xin(nctrl),yin(nctrl),xt(ntarg),width,yfit(ntarg),y0
    ## real(8)             :: xs(nctrl),ys(nctrl),targ,fitval,lwmean1,lwmean2,newwidth
    ## integer             :: order,nrobust
    ## integer             :: np 
    ## character(1)        :: wintype
    ## logical             :: monotonic,yfitnan
    ## integer             :: i,ii,iii,j,k
    
    y0 <- 0.0
    xs<-xin
    ys<-yin
    ntarg<-length(xt)
    nctrl<-length(xin)
    yfit<-array(NA,ntarg)
    stderr<-array(NA,ntarg)
    
    #!write (*,*) nctrl,ntarg
    #!write (*,'(8f10.2)') xs
    #!write (*,'(8f10.4)') ys
    #!write (*,'(8f10.2)') xt
    #!pause
    
    #! sort the input data by x-values if necessary
    # monotonic=.true.
    # do i=2,nctrl
    #    if (xs(i).lt.xs(i-1)) then
    #        monotonic=.false.
    #        ! write (*,'("xs(i),xs(i-1) non-monotonic: ",i6,2f14.6)') i,xs(i),xs(i-1)
    #    end if
    #end do
    #if (.not. monotonic) then
    #    ! write (*,'(a)') "Note:  input x's were not monotonic"
    #    ! write (10,'(a)') "Note:  input x's were not monotonic"
    #    call sortxy(nctrl,xin,yin,xs,ys)  !sortxy replaced by R command 'order'
    #end if
    yin<-yin[order(xin)]
    xs<-xs[order(xin)]
    ys<-ys[order(xin)]
    xin<-xin[order(xin)]

    # ! loop over target points, calling localfit for each target point
    for(ii in 1:ntarg){    # do ii=1,ntarg  
      targ<-xt[ii]     #    targ=xt(ii)  
      #  call localfit(nctrl,xs,ys,targ,wintype,width,order,nrobust,fitval,lwmean1)
      #    yfit(ii)=fitval

      fitted<-localfit(nctrl,xs,ys,targ,wintype,width,order,nrobust)
      yfit[ii]<-fitted[1]
      stderr[ii]<-fitted[2]
      }  #end do  
    #    !write (*,*) ii,targ,yfit(ii),lwmean1
    return(data.frame(xt,yfit,stderr))
    }  #end subroutine lowess_curve

# subroutine localfit(nctrl,xs,ys,targval,wintype,width,order,nrobust,fitval)
# ! gets a lowess-type fitted value at a particular target point
localfit <- function(nctrl,xs,ys,targval,wintype,width,order,nrobust){
    
    #implicit none
    #integer             :: nctrl
    #real(8)             :: targval,width,fitval,lwmean
    #real(8)             :: xs(nctrl),ys(nctrl),x(nctrl),y(nctrl)
    #real(8)             :: w(nctrl),wsum,ywm,ym,ystd
    #real(8)             :: b0,b1,b2,yhat(nctrl),resid(nctrl),mse,stderr,rsq
    #real(8)             :: dmax,diff,u,ul,ll
    #real(8)             :: r(nctrl),median,rmedx6,ru,rw(nctrl)
    #integer             :: order,nrobust
    #integer             :: np,np2,npts,low,high,nonzero
    #character(1)        :: wintype
    #integer             :: i,ii,iii,j
    #!write (*,*) nctrl
    #!write (*,'(8f10.2)') xs
    #!write (*,'(8f10.4)') ys
    #!write (*,'(8f10.2)') targval
    
    x<-array(NA,nctrl)
    y<-array(NA,nctrl)
    r<-array(NA,nctrl)
    rw<-array(NA,nctrl)
    
    #! select data to use for a specific fitted value
    if(wintype=="s"){     #if (wintype.eq.'s') then
      #! find observations in width proportion of data
      npts<-floor(width*nctrl)  #npts=int(width*nctrl)
      np2<-npts/2   #np2=npts/2
      low<-ii-np2   #low=ii-np2
      high<-ii+np2   #high=ii+np2
      if(low<1){    #if (low.lt.1) then
        low<-1    #low=1
        high<-low+npts    #high=low+npts
        }  # end if
      if(high>nctrl){   #if (high.gt.nctrl) then
        high<-nctrl    #high=nctrl
        low<-nctrl-npts    #low=nctrl-npts
        }   #end if
      np<-high-low+1    #np=high-low+1
      i<-0     #i=0
      for(iii in low:high){      #do iii=low,high
        i<-i+1    #i=i+1
        x[i]<-xs[iii]   #x(i)=xs(iii)
        y[i]<-ys[iii]    #y(i)=ys(iii)
        }   #end do
        
      # ! get distance weights
      dmax<-0    #dmax=0.0d0
      for(i in 1:np){    #do i=1,np
        diff<-x[i]-targval   #diff=x(i)-targval
        if(abs(diff)>dmax) dmax<-abs(diff)    #if (dabs(diff).gt.dmax) dmax=dabs(diff)
        }  #end do
      
 w<-array(NA,np)
     for(i in 1:np){     #do i=1,np
        u<-abs(x[i]-targval)/dmax    #u=dabs(x(i)-targval)/dmax
        if(u>1) u<-1.0     #if (u.gt.1.0d0) u=1.0d0
        w[i]<-(1-u^3)^3   #w(i)=(1.0d0-u**3.0d0)**3.0d0
        }  #end do
          
      } else {    #else if (wintype.eq.'w') then  #assume 'w'
      #! find observations within plus or minus window half-width (width)
      ul<-targval+width    #ul=targval+width
      ll<-targval-width  #ll=targval-width
      i<-0   #i=0
      for(iii in 1:nctrl){   #do iii=1,nctrl
      	if(!is.na(xs[iii])){
        if(xs[iii]>ll && xs[iii]<ul){    #if (xs(iii).ge.ll .and. xs(iii).le.ul) then
          i<-i+1   #i=i+1
          x[i]<-xs[iii]    #x(i)=xs(iii)
          y[i]<-ys[iii]
          #! write (*,*) i,x(i),y(i)
          }  #end if
          }
        }   #end do
      np<-i   #np=i   
    	if(np==0) return(c(NA,NA))  #if no points in window, return nulls
    	
      #!write (*,*) np
      #! get distance weights
      dmax<-width    #dmax=width
      w<-array(NA,np)
      for(i in 1:np){     #do i=1,np
        u<-abs(x[i]-targval)/dmax   #u=dabs(x(i)-targval)/dmax
        if(u>1) u<-1.0    #if (u.gt.1.0d0) u=1.0d0
        w[i]<-(1-u^3)^3    #w(i)=(1.0d0-u**3.0d0)**3.0d0
        }  #end do         
      }    #end if

    # ! get a locally weighted mean
    # call wm(np,y,w,lwmean,resid)
    resid<-wm_resid(np,y,w)
    
    fitval<-NA     # fitval=0.0d0/0.0d0 ! set to NaN initially
    seval<-NA

    # ! set initial robustness weights equal to distance weights
    nonzero<-0  #nonzero=0
    for(i in 1:np){   #do i=1,np
      rw[i]<-w[i]   #rw(i)=w(i)
      if(rw[i] != 0.0) nonzero <- nonzero+1 #if (rw(i).ne.0.0d0) nonzero=nonzero+1
      #!write (*,*) i,w(i),rw(i)
      }  #end do

    #! robustness iterations
    for(j in 0:nrobust){    #do j=0,nrobust

      #! lowess regressions
      #select case (order)
            
      if(order==0){   #case (0)  #! (distance) weighted mean
        fitval <- wm(np,y,rw)  ## call wm(np,y,rw,ywm,resid)
        seval <- wm_resid(np,y,w)
           # fitval=ywm
        } #end if
      if(order==1){    #case (1)   #! distance-weighted linear regression
        if(nonzero>3){   #if (nonzero.gt.3) then
          #call fitline(np,x,y,rw,b0,b1,yhat,resid,mse,stderr,rsq)
          fitout <- fitline(np,x,y,rw)
          fitval <- fitout$b0+fitout$b1*targval
          seval <- fitout$stderr
          }  #end if
        }  #end if
      if(order==2){   #case (2)
        #! distance-weighted quadratic regression
        if(nonzero>8){     #if (nonzero.gt.8) then
          #call fitquad(np,x,y,rw,b0,b1,b2,yhat,resid,mse,stderr,rsq)
          fitout <- fitquad(np,x,y,rw)
          #fitval=b0+b1*targval+b2*targval*targval
          fitval<-fitout$b0+fitout$b1*targval+fitout$b2*targval*targval
          seval<-fitout$stderr
          }   #end if
        }  #end if
      #case default
      #    write (*,*) "Polynomial order not supported...."
      #    stop 
      if(order>2){
        print("Polynomial order not supported....")
        break
        }
              
      #end select
        
      #! robustness weights
      if(nrobust!=0 && j<nrobust){   #if (nrobust.ne.0 .and. j.lt.nrobust) then
    
        #! median of absolute values of current residuals
        #call medianx(np,dabs(resid),median)
        rmedx6<-median(abs(resid))*6    #rmedx6=median*6.0d0
        
        #! bisquare weights (r) and new robustness weights (rw)
        nonzero<-0    #nonzero=0
        for(i in 1:np){      #do i=1,np
          ru<-abs(resid[i]/rmedx6)  #ru=dabs(resid(i)/rmedx6)    
          if(ru>1) ru<-1.0   #if (ru.gt.1.0d0) ru=1.0d0
          r[i]<-(1-ru^2)^2   #r(i)=(1.0d0-ru**2.0d0)**2.0d0
          rw[i]<-w[i]*r[i]    #rw(i)=w(i)*r(i)
          if(rw[i] != 0) nonzero <- nonzero+1  #if (rw(i).ne.0.0d0) nonzero=nonzero+1
          }  #end do
          
        } #end if           
      }  #end do (next j , robust step)
        
    #!write (*,*) lwmean,fitval
    #!pause
    return(c(fitval,seval))
    }    #end subroutine localfit

# subroutine wm(npts,y,w,ywm,dev)
wm <- function(npts,y,w){  
    #! weighted mean  --this function returns the weighted mean

    #implicit none

    #integer     :: npts
    #real(8)     :: y(npts),w(npts),dev(npts)
    #real(8)     :: wsum,ywm,ym
    #integer     :: i
    dev<-array(NA,npts)

    ym<-0.0    #ym=0.0d0
    ywm<-0.0   #ywm=0.0d0
    wsum<-0.0  #wsum=0.0d0
    for(i in 1:npts){   #do i=1,npts
      ym<-ym+y[i]  #ym=ym+y(i)
      ywm<-ywm+y[i]*w[i]   #ywm=ywm+y(i)*w(i)
      wsum<-wsum+w[i]    #wsum=wsum+w(i)
      }  #end do
    if(wsum>0.0){      #if (wsum.gt.0.0) then
      ywm<-ywm/wsum
      } else {
      ywm<-ym/npts   #ywm=ym/dble(npts)
      }  #end if

    return(ywm)   
    }   #end subroutine wm
    
wm_resid <- function(npts,y,w){  
    #! weighted mean  --this function returns the deviations

    #implicit none

    #integer     :: npts
    #real(8)     :: y(npts),w(npts),dev(npts)
    #real(8)     :: wsum,ywm,ym
    #integer     :: i
    dev<-array(NA,npts)

    ym<-0.0    #ym=0.0d0
    ywm<-0.0   #ywm=0.0d0
    wsum<-0.0  #wsum=0.0d0
    for(i in 1:npts){   #do i=1,npts
      ym<-ym+y[i]  #ym=ym+y(i)
      ywm<-ywm+y[i]*w[i]   #ywm=ywm+y(i)*w(i)
      wsum<-wsum+w[i]    #wsum=wsum+w(i)
      }  #end do
    if(wsum>0.0){      #if (wsum.gt.0.0) then
      ywm<-ywm/wsum
      } else {
      ywm<-ym/npts   #ywm=ym/dble(npts)
      }  #end if

    for(i in 1:npts){    #do i=1,npts
      dev[i]<-y[i]-ywm    #dev(i)=y(i)-ywm
      }  #end do
    
    return(dev)   
    }   #end subroutine wm


#subroutine fitline(npts,x,y,w,b0,b1,yhat,resid,mse,stderr,rsq)
fitline <- function(npts,x,y,w){
    #! simple weighted linear regression

    #implicit none
    #integer         :: npts
    #real(8)         :: x(npts),y(npts),w(npts),b0,b1,yhat(npts),resid(npts)
    #real(8)         :: xm,ym,wsum,sxx,sxy,sse,sst,mse,stderr,sratio,rsq
    #integer         :: i
    yhat<-array(NA,npts)
    resid<-array(NA,npts)

    #xm=0.0d0; ym=0.0d0; wsum=0.0d0; sxx=0.0d0; sxy=0.0d0
    xm<-0
    ym<-0
    wsum<-0
    sxx<-0
    sxy<-0
    
    for(i in 1:npts){    #do i=1,npts
      wsum<-wsum+w[i]  #wsum=wsum+w(i)
      xm<-xm+x[i]*w[i]  #xm=xm+x(i)*w(i)
      ym<-ym+y[i]*w[i]   #ym=ym+y(i)*w(i)
      #!write (*,*) i,x(i),y(i),w(i)
      }  #end do
    xm<-xm/wsum  #xm=xm/wsum
    ym<-ym/wsum   #ym=ym/wsum
    #!write (*,*) xm,ym,wsum

    for(i in 1:npts){   #do i=1,npts
      sxx<-sxx+((x[i]-xm)^2.0)*w[i]  #sxx=sxx+((x(i)-xm)**2.0d0)*w(i)
      sxy<-sxy+(x[i]-xm)*(y[i]-ym)*w[i]   #sxy=sxy+(x(i)-xm)*(y(i)-ym)*w(i)
      }  #end do

    b1<-sxy/sxx    #b1=sxy/sxx
    b0<-ym-(b1*xm)  #b0=ym-(b1*xm)

    sse<-0.0
    sst<-0.0    #sse=0.0d0; sst=0.0d0
    for(i in 1:npts){   #do i=1,npts
      yhat[i]<-b0+b1*x[i]   #yhat(i)=b0+b1*x(i)
      resid[i]<-y[i]-yhat[i]  #resid(i)=y(i)-yhat(i)
      sse<-sse+(resid[i]*sqrt(w[i]))^2.0  #sse=sse+(resid(i)*dsqrt(w(i)))**2.0d0
      sst<-sst+((y[i]-ym)*sqrt(w[i]))^2.0  #sst=sst+((y(i)-ym)*dsqrt(w(i)))**2.0d0
      } #end do

    mse<-sse/(npts-2)  #mse=sse/dble(npts-2)
    stderr<-sqrt(mse)  #stderr=dsqrt(mse)
    rsq <- -1.0        #rsq=-1.0d0
    sratio<-sse/sst   #sratio=sse/sst
    #if (sratio.le.1.0d0 .and. sratio.ge.0.0d0) rsq=1.0d0-sratio
    if(sratio<1.0 && sratio>0.0) rsq<-1.0-sratio
    
    return(data.frame(b0,b1,mse,stderr,rsq))
    }  #end subroutine fitline


#subroutine fitquad(npts,x1,y,w,b0,b1,b2,yhat,resid,mse,stderr,rsq)
#! simple weighted quadratic regression
fitquad <- function(npts,x1,y,w){
    #implicit none

    #integer         :: npts
    #real(8)         :: x1(npts),x2(npts),y(npts),w(npts),b0,b1,b2,yhat(npts),resid(npts)
    #real(8)         :: x1m,x2m,ym,wsum,sx1x1,sx2x2,sx1x2,sx1y,sx2y,sse,sst,mse,stderr,sratio,rsq
    #real(8)         :: bx1y,bx2y,bx2x1,bx1x2
    #integer         :: i
    x2<-array(NA,npts)
    yhat<-array(NA,npts)
    resid<-array(NA,npts)

    #x1m=0.0d0; x2m=0.0d0; ym=0.0d0; wsum=0.0d0
    x1m<-0.0
    x2m<-0.0
    ym<-0.0
    wsum<-0.0
    sx1x1<-0.0
    sx2x2<-0.0
    sx1x2<-0.0
    sx1y<-0.0
    sx2y<-0.0
    
    for(i in 1:npts){    #do i=1,npts
      x2[i]<-x1[i]^2  #x2(i)=x1(i)**2.0d0
      wsum<-wsum+w[i]   #wsum=wsum+w(i)
      x1m<-x1m+x1[i]*w[i]  #x1m=x1m+x1(i)*w(i)
      x2m<-x2m+x2[i]*w[i]   #x2m=x2m+x2(i)*w(i)
      ym<-ym+y[i]*w[i]    #ym=ym+y(i)*w(i)
      }    #end do
    x1m<-x1m/wsum     #x1m=x1m/wsum
    x2m<-x2m/wsum     #x2m=x2m/wsum
    ym<-ym/wsum     #ym=ym/wsum

    for(i in 1:npts){     #do i=1,npts
      sx1x1<-sx1x1+((x1[i]-x1m)^2)*w[i]  #sx1x1=sx1x1+((x1(i)-x1m)**2.0d0)*w(i)
      sx2x2<-sx2x2+((x2[i]-x2m)^2)*w[i]  #sx2x2=sx2x2+((x2(i)-x2m)**2.0d0)*w(i)
      sx1x2<-sx1x2+(x1[i]-x1m)*(x2[i]-x2m)*w[i]   #sx1x2=sx1x2+(x1(i)-x1m)*(x2(i)-x2m)*w(i)
      sx1y<-sx1y+(x1[i]-x1m)*(y[i]-ym)*w[i]     #sx1y=sx1y+(x1(i)-x1m)*(y(i)-ym)*w(i)
      sx2y<-sx2y+(x2[i]-x2m)*(y[i]-ym)*w[i]    #sx2y=sx2y+(x2(i)-x2m)*(y(i)-ym)*w(i)
      }  #end do

    bx1y<-sx1y/sx1x1    #bx1y=sx1y/sx1x1    ! b12
    bx2y<-sx2y/sx2x2    #bx2y=sx2y/sx2x2    ! b13
    bx2x1<-sx1x2/sx2x2  #bx2x1=sx1x2/sx2x2  ! b32
    bx1x2<-sx1x2/sx1x1  #bx1x2=sx1x2/sx1x1  ! b23

    b1<-(bx1y-bx2y*bx1x2)/(1.0-bx1x2*bx2x1)   #b1=(bx1y-bx2y*bx1x2)/(1.0d0-bx1x2*bx2x1)
    b2<-(bx2y-bx1y*bx2x1)/(1.0-bx1x2*bx2x1)   #b2=(bx2y-bx1y*bx2x1)/(1.0d0-bx1x2*bx2x1)

    b0<-ym-(b1*x1m)-(b2*x2m)   #b0=ym-(b1*x1m)-(b2*x2m)

    sse<-0.0   #sse=0.0d0; sst=0.0d0
    sst<-0.0
    for(i in 1:npts){    #do i=1,npts
      yhat[i]<-b0+b1*x1[i]+b2*x2[i]   #yhat(i)=b0+b1*x1(i)+b2*x2(i)
      resid[i]<-y[i]-yhat[i]   #resid(i)=y(i)-yhat(i)
      sse<-sse+(resid[i]*sqrt(w[i]))^2   #sse=sse+(resid(i)*dsqrt(w(i)))**2.0d0
      sst<-sst+((y[i]-ym)*sqrt(w[i]))^2    #sst=sst+((y(i)-ym)*dsqrt(w(i)))**2.0d0
      }  #end do

    mse<-sse/(npts-2)   #mse=sse/dble(npts-2)
    stderr<-sqrt(mse)   #stderr=dsqrt(mse)
    rsq<- -1.0    #rsq=-1.0d0
    sratio<-sse/sst   #sratio=sse/sst
    if(sratio<1 && sratio>0) rsq <- 1-sratio  #if (sratio.le.1.0d0 .and. sratio.ge.0.0d0) rsq=1.0d0-sratio

    return(data.frame(b0,b1,b2,yhat,mse,stderr,rsq))
    }  #end subroutine fitquad
    

# end module lowess