X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=sslvrg.f;fp=sslvrg.f;h=221beaf487c2ff8a8f11de3b2da7a8c74fad6280;hb=d0051dc9939d3477bd92b42c86bcd3eda743b955;hp=0000000000000000000000000000000000000000;hpb=7f0cae4f4853cc3f12bc751ee06ea31c7c97496e;p=mothur.git

diff --git a/sslvrg.f b/sslvrg.f
new file mode 100644
index 0000000..221beaf
--- /dev/null
+++ b/sslvrg.f
@@ -0,0 +1,136 @@
+C Output from Public domain Ratfor, version 1.0
+      subroutine sslvrg(penalt,dofoff,x,y,w,ssw, n, knot,nk,coef,
+     *     sz,lev, crit,icrit, lambda, xwy, hs0,hs1,hs2,hs3,
+     *     sg0,sg1,sg2,sg3, abd,p1ip,p2ip,ld4,ldnk,info)
+
+C Purpose :
+C       Compute smoothing spline for smoothing parameter lambda
+C       and compute one of three `criteria' (OCV , GCV , "df match").
+C See comments in ./sbart.f from which this is called
+
+      integer n,nk,icrit,ld4,ldnk,info
+      DOUBLE precision penalt,dofoff,x(n),y(n),w(n),ssw,
+     &     knot(nk+4), coef(nk),sz(n),lev(n), crit, lambda,
+     *     xwy(nk), hs0(nk),hs1(nk),hs2(nk),hs3(nk),
+     *     sg0(nk),sg1(nk),sg2(nk),sg3(nk), abd(ld4,nk),
+     &     p1ip(ld4,nk),p2ip(ldnk,nk)
+
+      EXTERNAL bvalue
+      double precision bvalue
+C local variables
+      double precision vnikx(4,1),work(16)
+      integer i,ileft,j,mflag, lenkno, ilo
+      double precision b0,b1,b2,b3,eps, xv,rss,df, sumw
+c
+c      integer interv
+c      external interv
+
+      lenkno = nk+4
+      ileft = 1
+      eps = 1d-11
+	  ilo = 1
+
+C compute the coefficients coef() of estimated smooth
+
+      do 1 i=1,nk
+         coef(i) = xwy(i)
+         abd(4,i) = hs0(i)+lambda*sg0(i)
+ 1    continue
+
+      do 4 i=1,(nk-1)
+         abd(3,i+1) = hs1(i)+lambda*sg1(i)
+ 4    continue
+
+      do 6 i=1,(nk-2)
+ 6       abd(2,i+2) = hs2(i)+lambda*sg2(i)
+
+      do 8 i=1,(nk-3)
+ 8       abd(1,i+3) = hs3(i)+lambda*sg3(i)
+
+c     factorize banded matrix abd:
+      call dpbfa(abd,ld4,nk,3,info)
+      if(info.ne.0) then
+C        matrix could not be factorized -> ier := info
+         return
+      endif
+c     solve linear system (from factorize abd):
+      call dpbsl(abd,ld4,nk,3,coef)
+
+C     Value of smooth at the data points
+      do 12 i=1,n
+         xv = x(i)
+ 12      sz(i) = bvalue(knot,coef,nk,4,xv,0)
+
+C     Compute the criterion function if requested
+
+      if(icrit .eq. 0)then
+         return
+      else
+C --- Ordinary or Generalized CV or "df match" ---
+
+C     Get Leverages First
+         call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0)
+         do 16 i=1,n
+            xv = x(i)
+            call interv(knot(1), nk+1, xv, ileft, mflag)
+            if(mflag .eq. -1) then
+               ileft = 4
+               xv = knot(4)+eps
+            else if(mflag .eq. 1) then
+               ileft = nk
+               xv = knot(nk+1) - eps
+            endif
+            j=ileft-3
+C           call bspvd(knot,4,1,xv,ileft,4,vnikx,work)
+            call bsplvd(knot,lenkno,4,xv,ileft,work,vnikx,1)
+            b0=vnikx(1,1)
+            b1=vnikx(2,1)
+            b2=vnikx(3,1)
+            b3=vnikx(4,1)
+            lev(i) = (
+     &              p1ip(4,j)*b0**2   + 2.d0*p1ip(3,j)*b0*b1 +
+     *           2.d0*p1ip(2,j)*b0*b2   + 2.d0*p1ip(1,j)*b0*b3 +
+     *              p1ip(4,j+1)*b1**2 + 2.d0*p1ip(3,j+1)*b1*b2 +
+     *           2.d0*p1ip(2,j+1)*b1*b3 +    p1ip(4,j+2)*b2**2 +
+     &           2.d0*p1ip(3,j+2)*b2*b3 +    p1ip(4,j+3)*b3**2
+     &           )*w(i)**2
+ 16      continue
+
+C     Evaluate Criterion
+
+         if(icrit .eq. 1)then
+C     Generalized CV
+            rss = ssw
+            df = 0d0
+            sumw = 0d0
+c	w(i) are sqrt( wt[i] ) weights scaled in ../R/smspline.R such
+c       that sumw =  number of observations with w(i) > 0
+            do 24 i=1,n
+               rss = rss + ((y(i)-sz(i))*w(i))**2
+               df = df + lev(i)
+               sumw = sumw + w(i)**2
+ 24         continue
+
+            crit = (rss/sumw)/((1d0-(dofoff + penalt*df)/sumw)**2)
+c            call dblepr("spar", 4, spar, 1)
+c            call dblepr("crit", 4, crit, 1)
+
+         else if(icrit .eq. 2) then
+C     Ordinary CV
+               crit = 0d0
+               do 30 i = 1,n
+ 30               crit = crit + (((y(i)-sz(i))*w(i))/(1-lev(i)))**2
+               crit = crit/n
+c            call dblepr("spar", 4, spar, 1)
+c            call dblepr("crit", 4, crit, 1)
+            else
+C     df matching
+            crit = 0d0
+            do 32 i=1,n
+ 32            crit = crit+lev(i)
+            crit = 3 + (dofoff-crit)**2
+         endif
+         return
+      endif
+C     Criterion evaluation
+      end