////////////////////////////////////////////////////////////////////////////////
//																			  //
// AUTHOR:      Tevfik Metin Sezgin	                        				  //
//				Massachusetts Institute of Technology						  //
//				Department of Electrical Engineering and Computer Science	  //
//				Artificial Intelligence Laboratory							  //
//																			  //
// E-MAIL:		mtsezgin@ai.mit.edu, mtsezgin@mit.edu						  //
//																			  //
// COPYRIGHT:   Tevfik Metin Sezgin                    						  //
//              All rights reserved. This code can not be copied, modified,   //
//              or distributed in whole or partially without the written 	  //
//				permission of the author. Also see the COPYRIGHT file.		  //
//																			  //
////////////////////////////////////////////////////////////////////////////////

      PROGRAM SAMPLE

C  ODRPACK ARGUMENT DEFINITIONS
C      ==> FCN      NAME OF THE USER SUPPLIED FUNCTION SUBROUTINE
C      ==> N        NUMBER OF OBSERVATIONS 
C      ==> M        COLUMNS OF DATA IN THE EXPLANATORY VARIABLE
C      ==> NP       NUMBER OF PARAMETERS
C      ==> NQ       NUMBER OF RESPONSES PER OBSERVATION
C     <==> BETA     FUNCTION PARAMETERS
C      ==> Y        RESPONSE VARIABLE
C      ==> LDY      LEADING DIMENSION OF ARRAY Y
C      ==> X        EXPLANATORY VARIABLE
C      ==> LDX      LEADING DIMENSION OF ARRAY X
C      ==> WE       "EPSILON" WEIGHTS
C      ==> LDWE     LEADING DIMENSION OF ARRAY WE
C      ==> LD2WE    SECOND DIMENSION OF ARRAY WE
C      ==> WD       "DELTA" WEIGHTS
C      ==> LDWD     LEADING DIMENSION OF ARRAY WD
C      ==> LD2WD    SECOND DIMENSION OF ARRAY WD
C      ==> IFIXB    INDICATORS FOR "FIXING" PARAMETERS (BETA)
C      ==> IFIXX    INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X)
C      ==> LDIFX    LEADING DIMENSION OF ARRAY IFIXX
C      ==> JOB      TASK TO BE PERFORMED
C      ==> NDIGIT   GOOD DIGITS IN SUBROUTINE FCN RESULTS
C      ==> TAUFAC   TRUST REGION INITIALIZATION FACTOR
C      ==> SSTOL    SUM OF SQUARES CONVERGENCE CRITERION
C      ==> PARTOL   PARAMETER CONVERGENCE CRITERION
C      ==> MAXIT    MAXIMUM NUMBER OF ITERATIONS
C      ==> IPRINT   PRINT CONTROL
C      ==> LUNERR   LOGICAL UNIT FOR ERROR REPORTS
C      ==> LUNRPT   LOGICAL UNIT FOR COMPUTATION REPORTS
C      ==> STPB     STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT BETA
C      ==> STPD     STEP SIZES FOR FINITE DIFFERENCE DERIVATIVES WRT DELTA
C      ==> LDSTPD   LEADING DIMENSION OF ARRAY STPD
C      ==> SCLB     SCALE VALUES FOR PARAMETERS BETA
C      ==> SCLD     SCALE VALUES FOR ERRORS DELTA IN EXPLANATORY VARIABLE
C      ==> LDSCLD   LEADING DIMENSION OF ARRAY SCLD
C     <==> WORK     DOUBLE PRECISION WORK VECTOR
C      ==> LWORK    DIMENSION OF VECTOR WORK
C     <==  IWORK    INTEGER WORK VECTOR
C      ==> LIWORK   DIMENSION OF VECTOR IWORK
C     <==  INFO     STOPPING CONDITION 
 
C  PARAMETERS SPECIFYING MAXIMUM PROBLEM SIZES HANDLED BY THIS DRIVER
C     MAXN          MAXIMUM NUMBER OF OBSERVATIONS 
C     MAXM          MAXIMUM NUMBER OF COLUMNS IN EXPLANATORY VARIABLE
C     MAXNP         MAXIMUM NUMBER OF FUNCTION PARAMETERS
C     MAXNQ         MAXIMUM NUMBER OF RESPONSES PER OBSERVATION

C  PARAMETER DECLARATIONS AND SPECIFICATIONS
      INTEGER    LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE,
     +           LIWORK,LWORK,MAXM,MAXN,MAXNP,MAXNQ
      PARAMETER (MAXM=2,MAXN=100,MAXNP=2,MAXNQ=2,MAXALL=1000,
     +           LDY=MAXN,LDX=MAXN,
     +           LDWE=MAXN,LD2WE=MAXNQ,LDWD=MAXN,LD2WD=1,
     +           LDIFX=MAXN,LDSCLD=1,LDSTPD=1,
     +           LWORK=18 + 11*MAXNP + MAXNP**2 + MAXM + MAXM**2 + 
     +                 4*MAXN*MAXNQ + 6*MAXN*MAXM + 2*MAXN*MAXNQ*MAXNP +  
     +                 2*MAXN*MAXNQ*MAXM + MAXNQ**2 + 
     +                 5*MAXNQ + MAXNQ*(MAXNP+MAXM) + LDWE*LD2WE*MAXNQ,
     +           LIWORK=20+MAXNP+MAXNQ*(MAXNP+MAXM))

C  VARIABLE DECLARATIONS 
      INTEGER          I,INFO,IPRINT,J,JOB,L,LUNERR,LUNRPT,M,MAXIT,N,
     +                 NDIGIT,NP,NQ,NWS,NPTS,NTMP
      INTEGER          IFIXB(MAXNP),IFIXX(LDIFX,MAXM),IWORK(LIWORK)   
      DOUBLE PRECISION PARTOL,SSTOL,TAUFAC
      DOUBLE PRECISION BETA(MAXNP),SCLB(MAXNP),SCLD(LDSCLD,MAXM),
     +                 STPB(MAXNP),STPD(LDSTPD,MAXM),
     +                 WD(LDWD,LD2WD,MAXM),WE(LDWE,LD2WE,MAXNQ),
     +                 WORK(LWORK),X(LDX,MAXM),Y(LDY,MAXNQ),
     +                 XALL(MAXALL,MAXM),YALL(MAXALL,MAXNQ)
      EXTERNAL         FCN


C  SPECIFY DEFAULT VALUES FOR DODRC ARGUMENTS
      WE(1,1,1)  = -1.0D0
      WD(1,1,1)  = -1.0D0
      IFIXB(1)   = -1
      IFIXX(1,1) = -1
      JOB        = -1
      NDIGIT     = -1
      TAUFAC     = -1.0D0
      SSTOL      = -1.0D0
      PARTOL     = -1.0D0
      MAXIT      = -1
      IPRINT     = -1
      LUNERR     = -1
      LUNRPT     = -1
      STPB(1)    = -1.0D0
      STPD(1,1)  = -1.0D0
      SCLB(1)    = -1.0D0
      SCLD(1,1)  = -1.0D0
 
C  SET UP ODRPACK REPORT FILES
      LUNERR  =   0
      LUNRPT  =   0

      M          = 1
      NP         = 2
      NQ         = 1
      BETA(1)    = 1.0
      BETA(2)    = 1.0
C  READ PROBLEM DATA
C
C  N    number of input points to DODRC
C  NPTS number of input points, first number in the input file
C  M    set to 1 (only one explanatory data column)
C  NP   Number of parameters set to 2 (slope and y intercept)
C  NQ   Number of responses per observation set to 1 (one y for each x)
C  NWS  window size, second number in the input file
C  
C
C
      OPEN (UNIT=5,FILE='parameters')
      READ (5,FMT=*) NPTS,NWS
      DO 10 I=1,NPTS
         READ (5,FMT=*) (XALL(I,J),J=1,M),(YALL(I,L),L=1,NQ)
 10   CONTINUE

C  SPECIFY TASK AS EXPLICIT ORTHOGONAL DISTANCE REGRESSION
C                  WITH CENTRAL DIFFERENCE DERIVATIVES
C                  COVARIANCE MATRIX CONSTRUCTED WITH RECOMPUTED DERIVATIVES
C                  DELTA INITIALIZED BY USER
C                  NOT A RESTART
C  AND INDICATE LONG INITIAL REPORT
C               NO ITERATION REPORTS
C               LONG FINAL REPORT
      JOB     = 01010
      IPRINT  = 0000

      OPEN  (11, FILE='result' )
C  SET WEIGHTS
C  COMPUTE SOLUTION

      DO 13 I=2,NWS
         N    = I
C        WRITE (*,*) "I    =", I
C        WRITE (*,*) "N    =", N
C        WRITE (*,*) "NPTS =", NPTS
         DO 15 J=1,I
            X(J,1)=XALL(J,1) 
            Y(J,1)=YALL(J,1) 
 15      CONTINUE

         BETA(1)    = 1.0
         BETA(2)    = 1.0
         
         CALL DODRC(FCN,
     +        N,M,NP,NQ,
     +        BETA,
     +        Y,LDY,X,LDX,
     +        WE,LDWE,LD2WE,WD,LDWD,LD2WD,
     +        IFIXB,IFIXX,LDIFX,
     +        JOB,NDIGIT,TAUFAC,
     +        SSTOL,PARTOL,MAXIT,
     +        IPRINT,LUNERR,LUNRPT,
     +        STPB,STPD,LDSTPD,
     +        SCLB,SCLD,LDSCLD,
     +        WORK,LWORK,IWORK,LIWORK,
     +        INFO)
         
         WRITE (*,*) BETA(1), ' , ', BETA(2)
         WRITE (11,*) BETA(1), "   ", BETA(2)
 13   CONTINUE

      N       = NWS*2 + 1

      DO 18 I=0,NPTS-N
C        WRITE (*,*) "I    =", I
C        WRITE (*,*) "N    =", N
C        WRITE (*,*) "NPTS =", NPTS
         DO 20 J=1,N
            X(J,1)=XALL(I+J,1) 
            Y(J,1)=YALL(I+J,1) 
 20      CONTINUE

         BETA(1)    = 1.0
         BETA(2)    = 1.0
         
         CALL DODRC(FCN,
     +        N,M,NP,NQ,
     +        BETA,
     +        Y,LDY,X,LDX,
     +        WE,LDWE,LD2WE,WD,LDWD,LD2WD,
     +        IFIXB,IFIXX,LDIFX,
     +        JOB,NDIGIT,TAUFAC,
     +        SSTOL,PARTOL,MAXIT,
     +        IPRINT,LUNERR,LUNRPT,
     +        STPB,STPD,LDSTPD,
     +        SCLB,SCLD,LDSCLD,
     +        WORK,LWORK,IWORK,LIWORK,
     +        INFO)
         
         WRITE (*,*) BETA(1), ' , ', BETA(2)
         WRITE (11,*) BETA(1), "   ", BETA(2)
 18   CONTINUE

      NTMP = NWS
      DO 23 NTMP=0,NWS-2
         I    = NWS-NTMP
         N    = I
C        WRITE (*,*) "I    =", I
C        WRITE (*,*) "N    =", N
C        WRITE (*,*) "NPTS =", NPTS
         DO 25 J=1,I
            X(J,1)=XALL(NPTS-N+J,1) 
            Y(J,1)=YALL(NPTS-N+J,1) 
 25      CONTINUE

         BETA(1)    = 1.0
         BETA(2)    = 1.0
         
         CALL DODRC(FCN,
     +        N,M,NP,NQ,
     +        BETA,
     +        Y,LDY,X,LDX,
     +        WE,LDWE,LD2WE,WD,LDWD,LD2WD,
     +        IFIXB,IFIXX,LDIFX,
     +        JOB,NDIGIT,TAUFAC,
     +        SSTOL,PARTOL,MAXIT,
     +        IPRINT,LUNERR,LUNRPT,
     +        STPB,STPD,LDSTPD,
     +        SCLB,SCLD,LDSCLD,
     +        WORK,LWORK,IWORK,LIWORK,
     +        INFO)
         
         WRITE (*,*) BETA(1), ' , ', BETA(2)
         WRITE (11,*) BETA(1), "   ", BETA(2)
 23   CONTINUE

      END


      SUBROUTINE FCN(N,M,NP,NQ,
     +               LDN,LDM,LDNP,
     +               BETA,XPLUSD,
     +               IFIXB,IFIXX,LDIFX,
     +               IDEVAL,F,FJACB,FJACD,
     +               ISTOP)

C  SUBROUTINE ARGUMENTS
C      ==> N        NUMBER OF OBSERVATIONS
C      ==> M        NUMBER OF COLUMNS IN EXPLANATORY VARIABLE
C      ==> NP       NUMBER OF PARAMETERS
C      ==> NQ       NUMBER OF RESPONSES PER OBSERVATION
C      ==> LDN      LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING N
C      ==> LDM      LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING M
C      ==> LDNP     LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING NP
C      ==> BETA     CURRENT VALUES OF PARAMETERS
C      ==> XPLUSD   CURRENT VALUE OF EXPLANATORY VARIABLE, I.E., X + DELTA
C      ==> IFIXB    INDICATORS FOR "FIXING" PARAMETERS (BETA)
C      ==> IFIXX    INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X)
C      ==> LDIFX    LEADING DIMENSION OF ARRAY IFIXX
C      ==> IDEVAL   INDICATOR FOR SELECTING COMPUTATION TO BE PERFORMED
C     <==  F        PREDICTED FUNCTION VALUES
C     <==  FJACB    JACOBIAN WITH RESPECT TO BETA
C     <==  FJACD    JACOBIAN WITH RESPECT TO ERRORS DELTA
C     <==  ISTOP    STOPPING CONDITION, WHERE
C                     0 MEANS CURRENT BETA AND X+DELTA WERE
C                       ACCEPTABLE AND VALUES WERE COMPUTED SUCCESSFULLY
C                     1 MEANS CURRENT BETA AND X+DELTA ARE
C                       NOT ACCEPTABLE;  ODRPACK SHOULD SELECT VALUES  
C                       CLOSER TO MOST RECENTLY USED VALUES IF POSSIBLE 
C                    -1 MEANS CURRENT BETA AND X+DELTA ARE
C                       NOT ACCEPTABLE;  ODRPACK SHOULD STOP

C  INPUT ARGUMENTS, NOT TO BE CHANGED BY THIS ROUTINE:
      INTEGER          I,IDEVAL,ISTOP,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ
      DOUBLE PRECISION BETA(NP),XPLUSD(LDN,M)
      INTEGER          IFIXB(NP),IFIXX(LDIFX,M)
C  OUTPUT ARGUMENTS:
      DOUBLE PRECISION F(LDN,NQ),FJACB(LDN,LDNP,NQ),FJACD(LDN,LDM,NQ)
C  LOCAL VARIABLES
      INTRINSIC        ATAN2,EXP,SQRT


C  CHECK FOR UNACCEPTABLE VALUES FOR THIS PROBLEM
C  COMPUTE PREDICTED VALUES
      IF (MOD(IDEVAL,10).GE.1) THEN
         DO 100 I = 1,N
            F(I,1) = BETA(1)*XPLUSD(I,1) + BETA(2)
  100    CONTINUE
      END IF

      RETURN
      END
