c----------------- SUN/PC version, April 13,1993------------ program lenns implicit double precision(a-h,o-z) parameter(kmax=8,jdmax=10,nxmax=2025,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 integer tim1,tim2,time common /size/ k,jd,nx,m,jforce common /xdata/ xt,xlag external objfun tim1=time() c----------- Initialize random number generator iseed=time() 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 16 continue 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) c------------- if there was no convergence, try again if(inform.gt.0) goto 16 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 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 tim2=time() open(12,file=outname,status='old',access='append') write(12,*) '@ending time= ',tim2, ' elapsed time',tim2-tim1,'@' 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=2025,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,jd 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=2025,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