c----------------- SUN/PC version, April 13,1993------------
        program lenns
        implicit double precision(a-h,o-z)
        parameter(kmax=8,jdmax=10,nxmax=525,lmax=65)
        parameter(npmax=1+kmax*(jdmax+2),twopi=6.2831853)
        dimension xt(nxmax),xdat(nxmax),xlag(nxmax,jdmax),
     *      theta(npmax),svexp(lmax),qrexp(lmax),
     *      der(lmax,nxmax),tsav(20,npmax),rsav(20),xraw(nxmax)
        dimension h0(npmax,npmax)
        character*35 fname,outname
        common /size/ k,jd,nx,m,jforce
        common /xdata/ xt,xlag
        external objfun

c----------- Initialize random number generator
	write(*,*) 'enter seed for random numbers'
	write(*,*) 'integer between 1 and 150000'
	read(*,*) iseed
        
        rn=uran(-abs(mod(iseed,15000)))

c----------- Set initial Jacobian to identity matrix
	call setidmat(h0,npmax)
 
c----------- Read in parameters and minimization routine options
        open(11,file='lenns.par',status='old')
        read(11,10) fname
        read(11,10) outname
10      format(a35)
        read(11,*) nx,jt1,ltr,ifreq,jd1,jd2,k1,k2,kcrit,l1,l2,
     *             jtpsav,alf,jforce,jper,cgcv,scale,ftol1,ftol2,
     *             itmax1,itmax2,maxstep1,maxstep2,ibrent,fltol
        close(11)
	iprint=itmax1+itmax2

c----------- Check if jt1 is OK, and fix if necessary
	if (jtpsav.lt.0) then
	    jt1min=jd2*l2
	else
            jt1min=(jd2-1)*l2+jtpsav
	endif
	if (jt1.lt.jt1min) jt1=jt1min

c------------ read data series into xdat; set nx=total series length
       call getx(fname,xdat,nxmax,nx,ifreq,ltr)
	
c---------- Set up the output file
	call makeout(fname,outname,nx,jt1,ltr,ifreq,k1,k2,cgcv,
     *          kcrit,jtpsav,alf,scale,jforce,jper,ftol1,ftol2,
     *          itmax1,itmax2,maxstep1,maxstep2,ibrent,fltol)

c------------ Standardize values to predict: SD=1, mean=0 if alf=0
       call sdev(xdat,jt1+1,nx,sd,xbar)
       if (alf.eq.0.) then
         do 150 j=1,nx
           xdat(j)=xdat(j)-xbar
150      continue
       endif
       call dscal(nx,1./sd,xdat,1)

c------------------ START MAIN LOOP on time-delay and #lags
       do 1010 ldelay=l1,l2
       do 1005 jd=jd1,jd2

	jtp=jtpsav
	if (jtpsav.le.0) jtp=ldelay

        call recontp(xraw,xdat,xlag,nxmax,jdmax,nx,ldelay,jd,jt1,m,jtp)
c --------------- m is the # of state vectors (# rows in xlag)
	do 200 j=1,m
	    xt(j)=xraw(j)-alf*xlag(j,1)	    	
200	continue

c----------------- Add sin & cos terms if fitting forced model 
	if (jforce.gt.0) then
	    xper=float(jper)
	    do 210 j=1,m
		ajper=float(mod(j,jper))/xper
		xlag(j,jd+1)=cos(twopi*ajper)
		xlag(j,jd+2)=sin(twopi*ajper)
210	    continue
	endif
		
c---------------- Initialize markers for best # units 'so far'
        kbic=k1
        bicopt=10000.
	kgcv=k1
	gcvopt=10000.
	mbic=0
	mgcv=0
c------------------ START LOOP on # units
        do 1004 k=k1,k2

c ---------------- decide if it's time to quit increasing # units
          if (k.gt.(2+kbic)) mbic=mbic+1
	  if (k.gt.(2+kgcv)) mgcv=mgcv+1
	  if (kcrit.eq.0) kquit=mbic
	  if (kcrit.eq.1) kquit=mgcv
	  if (kcrit.eq.2) kquit=mbic*mgcv
	  if (kquit.gt.0) goto 1004
          npar=1+k*(jd+2+2*jforce)
          if (2*npar.gt.m) goto 1004
          if (npar.gt.npmax) goto 1004

          do 15 i=1,20
              rsav(i)=10000.
15        continue

c------------------ START LOOP on repeated fits
           do 99 irep=1,251

c----------------- initialize parameter choices
            eps=scale*10**( .01*(dfloat(irep)-126.) )
            call tinit(theta,npar,eps,100)

c------------------ do a high-tolerance fit
	    ibrent1=0	
            call bfgsfm(objfun,npar,npmax,theta,h0,ftol1,fltol,
     *              itmax1,maxstep1,fnew,iter,iprint,inform,ibrent1)

c            call npsol(npar,nclin,ncnln,nrowa,nrowj,nrowr,a,bl,bu,
c     *              confun,objfun,inform,iter,istate,c,cjac,clamda,
c     *              fnew,g,r,theta,iw,leniw,w,lenw)


            xm=float(m)
            bic=1.419+dlog(fnew)+0.5*float(npar)*dlog(xm)/xm
            imod=mod(irep-1,5)
            if(imod.eq.0) then 
                open (13,file='lenns.tmp',status='unknown',
     *                access='append')
                write(13,97) ldelay,jd,k,eps,fnew,bic,iter,inform
                write(*,97) ldelay,jd,k,eps,fnew,bic,iter,inform
                close(13)   
            endif
            imax=idamax(20,rsav,1)
            rmax=rsav(imax)
            if (fnew.lt.rmax) then
                rsav(imax)=fnew
                do 85 jpar=1,npar
                    tsav(imax,jpar)=theta(jpar)
85              continue
            endif
99       continue
c---------------------- END LOOP on repeated fits
                
         do 999 i=1,20
            do 997 jpar=1,npar
                    theta(jpar)=tsav(i,jpar)
997         continue

            call bfgsfm(objfun,npar,npmax,theta,h0,ftol2,fltol,
     *               itmax2,maxstep2,fnew,iter,iprint,inform,ibrent)

c            call npsol(npar,nclin,ncnln,nrowa,nrowj,nrowr,a,bl,bu,
c     *              confun,objfun,inform,iter,istate,c,cjac,clamda,
c     *              fnew,g,r,theta,iw,leniw,w,lenw)
                xm=float(m)
                bic=1.419+dlog(fnew)+0.5*float(npar)*dlog(xm)/xm
		gcv=(fnew/(1.-(cgcv*float(npar)/xm)))**2
                if (bic.lt.bicopt) then
                    kbic=k
                    bicopt=bic
                endif
		if (gcv.lt.gcvopt) then
		    kgcv=k
		    gcvopt=gcv
		endif

c------------------ Compute estimated LE if #lags > 0
		if (jd.gt.0) then
                    call dnetdx(npar,theta,der,ldelay,jtp,alf)
                    nlgs=jtp+(jd-1)*ldelay
                    call liapqr(nlgs,lmax,der,m,svexp,qrexp)
                    ebar=0.
                    ebarqr=0.
                    leq=igcdiv(ldelay,jtp)
                    do 686 jl=1,leq
                        ebar=ebar+svexp(jl)
                        ebarqr=ebarqr+qrexp(jl)
686                 continue
                    ebar=ebar/dfloat(leq)
                    ebarqr=ebarqr/dfloat(leq)
		else
		    ebar=-99.999
		    ebarqr=-77.777
		endif
                open(12,file=outname,status='old',access='append')
                write(12,98) ldelay,jd,k,fnew,bic,gcv,npar,iter,
     *                       inform,ebar,ebarqr
		close(12)
97          format(1x,3i3,2x,f7.3,2(2x,f9.6),2x,i5,i3)
98          format(1x,3i3,2x,f9.6,2(2x,f7.4),i3,i7,i3,2(2x,f7.3))

999         continue
           
1004    continue
c---------------------- END LOOP on #units
1005   continue
1010   continue
c---------------------- END MAIN LOOP on time delay and  #lags
       open(12,file=outname,status='old',access='append')
       write(12,*) '@ LENNS finished @'
       close(12)
       stop
       end

        subroutine dnetdx(npar,theta,der,ldelay,jtp,alf)
c
c       Input:
c           npar  = # of variables   (integer)
c           theta = Parameter values (vector of length nvar)
c           ldelay= time delay in attractor reconstruction
c           jtp   = prediction time
c           alf   = linear coefficient in most recent x(t) value
c        On output:
c        Der(i,j) contains the derivative of f(x_1,x_2,... ;theta)
c        with respect to x_i, evaluated at X=x( x(t-1), x(t-2),... ;theta)
c
        implicit double precision(a-h,o-z)
        parameter(kmax=8,jdmax=10,nxmax=525,lmax=65)
        parameter(npmax=1+kmax*(jdmax+2))
	parameter(zero=0.d0,one=1.d0,hf=0.5d0,two=2.d0)
        common /size/ k,jd,nx,m,jforce
        common /xdata/ xt,xlag
        dimension gama(kmax,jdmax),der(lmax,nxmax),h(jdmax),
     *   xt(nxmax),xlag(nxmax,jdmax),theta(npar),w(kmax,jdmax)
        do 15 i=1,lmax
            do 10 j=1,nxmax
                der(i,j)=zero
10          continue
15      continue
        m0=k+1
        j0=2*k+1
	jdf=jd+2*jforce
        do 30 iunit=1,k
            loc0=j0+(iunit-1)*jdf
            do 20 lag=1,jdf
                gama(iunit,lag)=theta(loc0+lag)
c               ---------- theta(1+iunit)=beta_i ---------
                w(iunit,lag)=theta(1+iunit)*gama(iunit,lag)
20          continue
30      continue
        do 100 it=1,m
            do 45 iunit=1,k
c               ----------  theta(m0+iunit)=mu_i --------
                uin=theta(m0+iunit)
                do 40 lag=1,jdf
                    uin=uin+gama(iunit,lag)*xlag(it,lag)
40              continue
                auin=abs(uin)
                dij=(two+auin +hf*(auin**2))
                h(iunit)=two*(one+auin)/(dij**2)
45          continue
            do 60 j=1,jdf
                djt=0.
                do 55 iunit=1,k
                    djt=djt+w(iunit,j)*h(iunit)
55              continue
                der(jtp+(j-1)*ldelay,it)=djt
60          continue
            der(jtp,it)=alf+der(jtp,it)
100     continue
        return
        end
  

       subroutine tinit(theta,npar,eps,nreps)
       implicit double precision(a-h,o-z)
       parameter(kmax=8,jdmax=10,nxmax=525,lmax=65)
       parameter(npmax=1+kmax*(jdmax+2))
       dimension theta(npmax),ti(npmax),g(npmax)
       twoeps=2.*eps
       fmin=10000.
       do 30 i=1,nreps
           do 10 j=1,npar
               ti(j)=twoeps*(uran(1)-0.5)
10         continue
           call objfun(0,npar,ti,fi,g,1)
           if(fi.lt.fmin) then
               call dcopy(npar,ti,1,theta,1)
               fmin=fi
           endif
30     continue
       return
       end

        subroutine setidmat(h,n)
	implicit real*8(a-h,o-z)
	dimension h(n,n)
	do 5 i=1,n
            do 3 j=1,n
               h(i,j)=0.
3           continue
            h(i,i)=1.
5       continue
	return
	end

	subroutine makeout(fname,outname,nx,jt1,ltr,ifreq,k1,k2,cgcv,
     *          kcrit,jtpsav,alf,scale,jforce,jper,ftol1,ftol2,
     *          itmax1,itmax2,maxstep1,maxstep2,ibrent,fltol)
	implicit double precision(a-h,o-z)
	character*35 fname,outname
        open(12,file=outname,status='unknown')
        write(12,100) fname
100     format(1x,'@ ',a35,'Output from LENNS.FOR')
        write(12,*) 'Dependent variable data scaled to Std Dev ==1'
	write(12,*) 'Parameter values:'
        write(12,110) nx,jt1,ltr,ifreq
110     format(i4,2x,3i3,9x,'<-- nx,jt1,ltr,ifreq')
        write(12,112) k1,k2,kcrit,cgcv
112	format(3(i3,2x),f6.3,3x,'<-- kmin,kmax,kcrit,cgcv')
	write(12,114) jtpsav,alf,scale
114	format(i3,2(2x,f5.2),7x,'<-- T_p,alfa,scale')
	write(12,115) jforce, jper
115	format(i3,2x,i5,14x,'<-- jforce, jper')
	write(12,116) ftol1,ftol2
116     format(e10.3,2x,e10.3,2x,'<-- ftol1,ftol2')
	write(12,118) itmax1,itmax2,maxstep1,maxstep2
118     format(1x,4(i5,2x),'<-- itmax1,itmax2,maxstep1,maxstep2')
	if (ibrent.gt.0) write(12,*) 'Polished with Brent,fltol= ',fltol
        write(12,*) 'output l,d,k,RMS err,BIC,GCV,#params,iters,
     *info,LE1(SV),LE1(QR) @'
        close(12)
	return
	end
