$NOFLOATCALLS
C     This program interpolates a function surface from data
C     it reads off an input file, and it allows interactive
C     exploration of the form of the interpolated surface.
C     The input file is ASCII, and free format, but must present
C     data in a certain order, and in three groups, each beginning on
C     a separate line (see below).  
C     NELTS:  the number of arguments to the function being interpolated --
C             the number of parameters in a likelihood application.
C     NOBS:   the number of observed points to be read in.
C     SCALE:  a NELTSx1 vector describing the size of the change in
C             argument which would be large enough to substantially 
C             reduce the likelihood that a second order Taylor series 
C             expansion will still apply.
C     B:      the array of observed argument values
C     FVAL:   the vector of observed function values
C     FBASE:  a "normal" level for the function, to which the interpolation
C             reverts at argument values far from any observed points.
C     FINCIN: a good increment to expect to achieve on a single iteration
C             from a point fairly far from the optimum.  
C     CRIT:   change in function value at which routine stops, assuming
C             convergence.  
C     The program allows you to change FBASE and SCALE, to display the input
C     function values in sorted order, to display the whole B and FVAL
C     matrices, and to compute the interpolated value and gradient at any
C     value of the argument.  After each function evaluation, it prompts you
C     as to whether you desire a line search.  This will involve the program
C     computing a Newton-Raphson step toward a maximum of the interpolated
C     function, if possible and if you tell it to attempt to use the Hessian.
C     otherwise it prompts you for a guess at the attainable function value
C     and computes a corresponding step along the gradient direction.  These
C     computed steps need not always result in an increase in the function
C     value.  In any case, you have the option of altering the step size along
C     the initial search direction, and a small enough step size will always
C     increase the interpolated function value.  The program gives you the
C     option of chaining into a new Newton-Raphson or gradient search starting
C     from the result of its current line search.
C     --------------------------------------------------------------
C      Additions 4/24/85
C      It is important that the scale factors be roughly correct in order
C      to get good results.  They should be chosen to be the same order of
C      magnitude as the differences among paramter values along the
C      corresponding dimension.  Ordinarily it's best to begin by
C      multiplying the scale factors by shrinking or expansion factors at
C      the start, experimenting until the likelihood value is
C      approximately maximized.  On very simple functions this likelihood
C      may turn out to be unbounded; this is likely if the function is
C      nearly quadratic or has no local maxima on the grid given.  In that
C      case be wary of inaccuracy if you make the scale factors very
C      large.  Smaller scale factors with slightly lower likelihood
C      will still give good interpolation in this kind of situation.  You
C      can detect inaccuracy by comparing interpolated values at input
C      data points to actual values.  The previous version of this
C      program did not calculate likelihood correctly.
C      
C      The DIST1 distance function uses a normal density as the filtering
C      kernel for the underlying stochastic process.  The new DIST2 uses a
C      "product of independent splines" form for the autocovariance
C      function.  It tends to do less well on very smooth functions, but
C      because it cuts dependence to zero faster with distance may work
C      better on functions which behave differently in different regions.
C     ----------------------------------------------------------------
      PROGRAM SMTHTS
      IMPLICIT REAL*8(A-H,O-Z)
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV, DET, FBASE, FIT
      INTEGER NELTS,NOBS,MAXELT,MAXOBS, MOBSYM
      LOGICAL LNFLAG
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C---------------------------------------------
      REAL*8 SCALE
      COMMON /SCALE/ SCALE(MAXELT)
      REAL*8 A(MAXELT),GVEC(MAXELT),NFBASE, CRIT, FINCIN
     & ,FINC, X, SMTH, F, FACTOR
      LOGICAL ANSYES,ABORT
      INTEGER ORDER(MAXOBS), IERR, IS5, IELT, I, IOBS, II, NIT
     & ,IT, J, N1, N2
      CHARACTER FNAME*80,C*1,UPCASE
C---------------------------------------------------
C     Set IERR to 1, since output file for parameters is not set yet.
      IERR = 1      
3     WRITE(*,1000) 'Input file name: '
1000  FORMAT(1X,A)
      READ(*,1001) FNAME
1001  FORMAT(A)
      OPEN(5,FILE=FNAME,STATUS='OLD',IOSTAT=IS5)
      IF (IS5.GT.0) THEN
         WRITE(*,*) 'Error opening ',FNAME
         GO TO 3
      ELSEIF (IS5.LT.0) THEN
         WRITE(*,*) 'Empty file ', FNAME
         GOTO 3
      ENDIF
      READ(5,*) CRIT 
      READ(5,*) NELTS
      READ(5,*) (SCALE(IELT),IELT=1,NELTS), FBASE, FINCIN
      WRITE(*,*) 'Transform to log parameters?'
      IF(ANSYES()) THEN 
         LNFLAG=.TRUE.
      ELSE
         LNFLAG = .FALSE.
      ENDIF
      DO 70 I=1,NELTS
70       SCALE(I)=1./SCALE(I)
      do 700 NOBS=1,maxobs 
         READ(5,*,END=701) (B(IELT,NOBS),IELT=1,NELTS),FVAL(NOBS)
700   continue
      NOBS=MAXOBS+1
701   CONTINUE
      NOBS=NOBS-1
      IF(LNFLAG) THEN
         DO 600 IOBS=1,NOBS
               CALL LOGB(B(1,IOBS),NELTS)
600      CONTINUE
      ENDIF
      CALL MMNSCL(FVAL,FBASE,NOBS)
      CALL BUBBLE(FVAL,ORDER,NOBS)
      CALL STBCVN
      WRITE(*,*) '2 x Likelihood = ', FIT
      II=0
C      DO 30 I=1,NOBS
C         WRITE(*,*) (TBCVNV(II+J),J=1,I)
C         II=II+I
C30    CONTINUE
11    WRITE(*,*) 'Function eval., Scale change, Print, Crit. change,'
      Write(*,1000) 'Display, Order(data points), Auto, or Quit: '
      READ(*,1001) C
      C=UPCASE(C)
C------------------------- Convergence Criterion Change --------
C
      IF(C .EQ. 'C') THEN
         WRITE(*,*) 'Current criterion is: ', CRIT,' New one: '
         READ(*,*) CRIT
C-------------------------- Print  -----------------------------
C
      ELSEIF (C .EQ. 'P') THEN
         WRITE(*,*) 'Open a new parameter file? '
         IF (ANSYES()) THEN
            WRITE(*,*) 'File name: '
            READ(*,'(A12)') FNAME
            OPEN(2,FILE=FNAME,STATUS='NEW',IOSTAT=IERR)
            IF(IERR.GT.0) THEN
               WRITE(*,*) 'Can''t open file  ',FNAME
            ENDIF
         ENDIF
         IF(IERR .EQ. 0 ) THEN
            IF(LNFLAG) THEN
               WRITE(2,'(1X,5G15.7)') (EXP(A(I)), I=1,NELTS)
            ELSE
               WRITE(2,'(1X,5G15.7)') (A(I),I=1,NELTS)
            ENDIF
         ELSE
            WRITE(*,*) 'Nothing written because there''s no file open.'
         ENDIF
C------------------------------ auto search ------------------	
      ELSEIF (C.EQ.'A') THEN
         WRITE(*,*) 'Start at data point?'
         IOBS=0
          IF(ANSYES()) THEN
             WRITE(*,*) 'Which one?'
             READ(*,*) IOBS
             IF (IOBS.GE.1) THEN
                CALL MATMOV(B(1,IOBS),A,NELTS)
             ENDIF
          ENDIF
          IF(IOBS.LT.1) THEN
             WRITE(*,'(1X,I4,A)') NELTS, '-parameter vector: '
             READ(*,*) (A(IELT),IELT=1,NELTS)
             IF (LNFLAG)  CALL LOGB(A,NELTS)
          ENDIF
          WRITE (*,*) 'Iteration bound: '
          READ( *,*) NIT
          IT=0
          FINC=FINCIN
300       CONTINUE
             X=SMTH(A)
             CALL GRAD(A,GVEC)
             WRITE(*,*) 'Function value: ',X
             CALL AUTNWT(A,F,FINC,ABORT)
             WRITE(*,*) IT,' Coeff. vector, function value: '
             IF(LNFLAG) THEN
                 WRITE(*,*) (EXP(A(J)),J=1,NELTS), F
             ELSE
                 WRITE(*,*) (A(J),J=1,NELTS),F
             ENDIF
C           Check convergence criteria.
             IF(ABORT) WRITE(*,*) 'Aborted search.'
             IF(IT.GE.NIT) WRITE(*,*) 'Iteration bound hit.'
             IF (IT.GE.NIT .OR. ABS(F-X).LT.CRIT.OR.ABORT) THEN
                 CALL GRAD(A,GVEC)               
             ELSE
               IT=IT+1
               GOTO 300
             ENDIF
C------------------------------ Function eval ------------------	
      ELSEIF (C.EQ.'F') THEN
10      CONTINUE
          WRITE(*,*) 'Start at a data point?'
          IF(ANSYES()) THEN
             WRITE(*,*) 'Which one? (0 if you don''t know): '
             READ(*,*) IOBS
             IF (IOBS.LT.1) THEN
                WRITE(*,*) 'Use display or order option.'
                GO TO 200
C               200: end of the 'F' block
             ELSE
                CALL MATMOV(B(1,IOBS),A,NELTS)
             ENDIF
          ELSE
             WRITE(*,'(1X,I4,A)') NELTS, '-parameter vector: '
             READ(*,*) (A(IELT),IELT=1,NELTS)
             IF(LNFLAG) CALL LOGB(A,NELTS)
          ENDIF
9         CONTINUE
             X=SMTH(A)
             CALL GRAD(A,GVEC)
             WRITE(*,*) 'Function value: ',X
             WRITE(*,*) 'Gradient: ',(GVEC(I),I=1,NELTS)
             WRITE(*,1000) 'Line search?'
             IF (ANSYES()) THEN
                CALL NEWT(A,ABORT)
                IF (.NOT.ABORT) GO TO 9
C               NEWT will have reset A, usually, so a new GRAD
C               evaluation will be helpful.
             ENDIF
200       WRITE(*,*) 'Input another point? (y/n): '
          IF(ANSYES()) GO TO 10
C------------------------------ Scale change ----------------------
      ELSEIF (C.EQ.'S') THEN
          WRITE(*,*) 'Existing Scales and base fcn value: '
          WRITE(*,*) (1./SCALE(J),J=1,NELTS), FBASE
		WRITE(*,*) '2 x Likelihood = ', FIT
          WRITE(*,*) 'Any changes?(y/n): '
          IF(ANSYES()) THEN
              WRITE(*,*) 'Multiply all scales by common factor?'
              IF(ANSYES()) THEN
                  WRITE(*,*) 'Factor: '
                  READ(*,*) FACTOR
                  IF (FACTOR.NE.1.) CALL SCALAR(SCALE,1./FACTOR,NELTS)
              ELSE
                  WRITE(*,*) 'New Scales: '
                  READ(*,*) (SCALE(J),J=1,NELTS)
                  DO 80 I=1,NELTS
80                   SCALE(I)=1./SCALE(I)
              ENDIF
              WRITE(*,*) 'Base change from ',FBASE, ' ?: '
              IF(ANSYES()) THEN
                  WRITE(*,*) 'Replacement: '
                  READ(*,*) NFBASE
                  CALL MMNSCL(FVAL,NFBASE-FBASE,NOBS)
                  FBASE=NFBASE
              ENDIF
              CALL STBCVN
	        WRITE(*,*) '2 x Likelihood = ', FIT
          ENDIF
C-------------------------- Display ----------------------
      ELSEIF (C.EQ.'D') THEN
          WRITE(*,*) 'All? (y/n): '
          IF(ANSYES()) THEN
              N1=1
              N2=NOBS
          ELSE
              WRITE(*,*) 'Range of observation numbers: '
              READ(*,*) N1,N2
          ENDIF
          IF(N1.GE.1.AND.N2.LE.NOBS) THEN
              DO 40 IOBS=N1,N2
                IF(LNFLAG) THEN
                  WRITE(*,*) IOBS,
     1               (EXP(B(J,IOBS)),J=1,NELTS),FVAL(IOBS)+FBASE
                ELSE
                  WRITE(*,*) IOBS,
     1               (B(J,IOBS),J=1,NELTS),FVAL(IOBS)+FBASE
                ENDIF
40            CONTINUE
          ENDIF
C------------------------ Order data ----------------------
      ELSEIF (C.EQ.'O') THEN
          WRITE(*,'(1X,3(I3,'': '',I3,1X,G16.8))')
     1     (J,ORDER(J),FVAL(ORDER(J))+FBASE,J=1,NOBS)
      ELSEIF (C.EQ.'Q') THEN
          STOP
      ENDIF
      GOTO 11
      END

      SUBROUTINE LOGB(A,NELTS)
      REAL*8 A(*)
      INTEGER I
      DO 10 I=1,NELTS
         A(I) = LOG(A(I))
10    CONTINUE
      RETURN
      END
