PROGRAM BESTSUBH2
      IMPLICIT INTEGER(A-Z)
      REAL*8 X(1000,100),XS(1000,100),Y(1000),R(100,100),TIMEA,
     1       MN(100),TIMEB,SD(100),SDY,MEANY,SUM,R2,RR(100,100),
     1       R2MIN, R2STORE(100),XT(1000,20),R2BEST(5),RDUM
      INTEGER S(100),SBEST(5,100),V(100,2),ORDER(100),SEL(100),
     1        SS(5,100),PERM(100),SDUM(100),SSEL(5,100)
C
C ########################################################################
C AUG 21, 2008: BEST SUBSETS REGRESSION - HIERARCHICALLY WELL-FORMULATED
C               SUBSETS FOR 2ND-ORDER POLYNOMIAL REEGRESSION USING
C               A BRANCH AND BOUND ALGORITHM
C     ***** STORE THE 5 BEST SUBSETS FOR EACH VALUE OF Q **********
C ########################################################################
C
      OPEN(5, FILE='PARAM2') ! N = OBJECTS, H = CANDIDATES, Q = SIZE, R2
      OPEN(1, FILE='DATAX.PRN')       ! DATA FILE - X VARIABLES
      OPEN(4, FILE='DATAY.PRN')       ! DATA FILE - Y VARIABLES
      OPEN(3, FILE='RESULT.OUT')      ! BRIEF OUTPUT FILE
      READ(5,*) N,H,QMIN,QMAX,(R2BEST(J),J=1,5)
C
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEA = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
      READ(4,*) (Y(I),I=1,N)
      DO I = 1,N
        READ(1,*) (X(I,J),J=1,H)
      END DO
      DO J = 1,H
        MN(J) = 0.0
        DO I = 1,N
          MN(J) = MN(J) + X(I,J)
        END DO
        MN(J) = MN(J)/DFLOAT(N)
        DO I = 1,N
          X(I,J) = X(I,J)/MN(J)
        END DO
      END DO
      P = H
      DO J = 1,H
        P = P+1
        DO I = 1,N
          X(I,P) = X(I,J)**2
        END DO
      END DO
      DO I1 = 1,H-1
        DO J1 = I1+1,H
          P = P + 1
          DO I = 1,N
            X(I,P) = X(I,I1)*X(I,J1)
          END DO
        END DO
      END DO
      DO J = 1,P
        MN(J) = 0.
        SD(J) = 0.
        DO I = 1,N
          MN(J) = MN(J) + X(I,J)
        END DO
        MN(J) = MN(J)/DFLOAT(N)
        DO I = 1,N
          SD(J) = SD(J) + (X(I,J)-MN(J))**2
        END DO
        SD(J) = SD(J)/DFLOAT(N-1)
        SD(J) = SD(J)**0.5
        DO I = 1,N
          XS(I,J) = (X(I,J)-MN(J))/SD(J)
        END DO
      END DO

C
C #####################################################################
C  COMPUTE THE CORRELATION MATRIX
C #####################################################################
C
      MEANY=0.0D0                       ! COMPUTE TSS FOR DEPENDENT
      DO I = 1,N                         
        MEANY = MEANY+Y(I)
      END DO
      MEANY=MEANY/DFLOAT(N)
      SDY = 0.
      DO I = 1,N
        SDY = SDY + (Y(I)-MEANY)**2
      END DO
      SDY = SDY/DFLOAT(N-1)
      SDY = SDY**0.5
      DO I = 1,N
        XS(I,P+1) = (Y(I)-MEANY)/SDY
      END DO
C
      DO K = 1,P+1
        DO L = 1,P+1
          SUM = 0.
          DO I = 1,N
            SUM = SUM + XS(I,K)*XS(I,L)
          END DO
          R(K,L) = SUM/DFLOAT(N-1)
        END DO
      END DO
C
      DO I = 1,P
        ICOL = I
        CALL SWEEP(P+1,ICOL,R)
      END DO
C
C      WRITE(*,150) 1. - R(P+1,P+1)
C
  150 FORMAT(' FULL MODEL R SQUARED WITH Q = P IS',F8.5)
C
      DO I = 1,H
        SEL(I) = 0
        ORDER(I) = 0
      END DO
      DO III = 1,H
        R2MIN = 9999.
        DO 620 JJJ = 1,H
          IF(SEL(JJJ).EQ.1) GO TO 620
          DO I = 1,P+1
            DO J = 1,P+1
              RR(I,J) = R(I,J)
            END DO
          END DO
          CALL SWEEP(P+1,JJJ,RR)
          CALL SWEEP(P+1,JJJ+H,RR)
          HH = 2*H
          DO 626 I = 1,H-1
            DO 627 J = I+1,H
              HH = HH + 1
              IF(SEL(I).EQ.1.OR.SEL(J).EQ.1) GO TO 627
              IF(JJJ.EQ.I.OR.JJJ.EQ.J) THEN
                CALL SWEEP(P+1,HH,RR)
              END IF
 627        CONTINUE
 626      CONTINUE
          IF(1-RR(P+1,P+1).LT.R2MIN) THEN
            R2MIN = 1-RR(P+1,P+1)
            JSEL = JJJ
          END IF
 620    CONTINUE
        ORDER(III) = JSEL
        SEL(JSEL) = 1
        PERM(III) = JSEL
        JSET = JSEL
        PERM(III+H) = JSEL+H
        CALL SWEEP(P+1,JSEL,R)
        CALL SWEEP(P+1,JSEL+H,R)
        HH = 2*H
        DO 624 I = 1,H-1
          DO 625 J = I+1,H
            HH = HH + 1
            IF(I.NE.JSET.AND.J.NE.JSET) GO TO 625
            IF(I.NE.JSET.AND.SEL(I).EQ.1) GO TO 625
            IF(J.NE.JSET.AND.SEL(J).EQ.1) GO TO 625
            CALL SWEEP(P+1,HH,R)
 625      CONTINUE
 624    CONTINUE
      END DO
C
C      WRITE(*,621) (ORDER(I),I=1,H)
C 621  FORMAT(20I3)
C
      DO J = 1,H
        DO I = 1,N
          XT(I,J) = X(I,ORDER(J))
        END DO
      END DO
      DO J = 1,H
        DO I = 1,N
          X(I,J) = XT(I,J)
        END DO
      END DO
      P = H
      DO J = 1,H
        P = P+1
        DO I = 1,N
          X(I,P) = X(I,J)**2
        END DO
      END DO
      DO I1 = 1,H-1
        DO J1 = I1+1,H
          P = P + 1
          DO I = 1,N
            X(I,P) = X(I,I1)*X(I,J1)
          END DO
          NNN = 2*H
          DO II = 1,H-1
            DO JJ = II+1,H
              NNN = NNN + 1
           IF((II.EQ.ORDER(I1).AND.JJ.EQ.ORDER(J1)).OR.(II.EQ.ORDER(J1)
     1          .AND.JJ.EQ.ORDER(I1))) THEN
                PERM(P) = NNN
              END IF
            END DO
          END DO
        END DO
      END DO
      DO J = 1,P
        MN(J) = 0.
        SD(J) = 0.
        DO I = 1,N
          MN(J) = MN(J) + X(I,J)
        END DO
        MN(J) = MN(J)/DFLOAT(N)
        DO I = 1,N
          SD(J) = SD(J) + (X(I,J)-MN(J))**2
        END DO
        SD(J) = SD(J)/DFLOAT(N-1)
        SD(J) = SD(J)**0.5
        DO I = 1,N
          XS(I,J) = (X(I,J)-MN(J))/SD(J)
        END DO
      END DO
      DO I = 1,N
        XS(I,P+1) = (Y(I)-MEANY)/SDY
      END DO
C
      DO I = 1,P+1
        DO J = 1,P+1
          R(I,J) = 0
        END DO
      END DO
C
      DO K = 1,P+1
        DO L = 1,P+1
          SUM = 0.
          DO I = 1,N
            SUM = SUM + XS(I,K)*XS(I,L)
          END DO
          RR(K,L) = SUM/DFLOAT(N-1)
        END DO
      END DO
C
      DO I = 1,P
        ICOL = I
        CALL SWEEP(P+1,ICOL,RR)
      END DO
C
      WRITE(*,150) 1. - RR(P+1,P+1)
      WRITE(3,150) 1. - RR(P+1,P+1)
C
      DO I = H+1,2*H
        V(I,1) = I - H
      END DO
      HH = 2*H
      DO I = 1,H-1
        DO J = 1+I,H
          HH = HH + 1
          V(HH,1) = I
          V(HH,2) = J
        END DO
      END DO
C
      DO 900 IJKL = QMIN,QMAX
C
      Q = IJKL
C
      DO I = 1,P+1
        DO J = 1,P+1
          R(I,J) = RR(I,J)
        END DO
      END DO
C      
C ###################################################################
C BEGIN THE BRANCH AND BOUND ALGORITHM
C ###################################################################
C
C   ######################
C   STEP 0: INITIALIZE
C   ######################
C
      IP=0
      DO K = 1,Q
        S(K)=0
      END DO
C
C   ###############################
C   STEP 1: INCREMENT SEARCH DEPTH
C   ###############################
C
 100  IP=IP+1
      IF(IP.EQ.1) THEN
        S(IP) = 1
        CALL SWEEP(P+1,S(IP),R)
      ELSE
        S(IP) = S(IP-1)+1
        CALL SWEEP(P+1,S(IP),R)
      END IF
C
C   ###############################
C   STEP 2: FEASIBILITY
C   ###############################
C
 200  IF(P-S(IP).LT.P-Q-IP) GO TO 500
C
C   ###############################
C   STEP 3: SUBOPTIMALITY
C   ###############################
C
 300  R2 = 1 - R(P+1,P+1)
      IF(R2.LT.R2BEST(5)) GO TO 500
C
C   ###############################
C   STEP 4: UPDATE INCUMBENT
C   ###############################
C
 400  IF(IP.NE.P-Q) GO TO 100
      DO K = P,S(IP)+1,-1
        DO L = 1,IP
          LL = S(L)
          DO M = 1,2
            IF(V(K,M).EQ.LL) GO TO 500
          END DO
        END DO
      END DO
      R2BEST(5) = R2
      DO I = 1,IP
        SBEST(5,I)=S(I)
      END DO
      DO I = 1,4
        IF(R2BEST(5).GT.R2BEST(I)) THEN
          RDUM = R2BEST(5)
          DO K = 1,IP
            SDUM(K) = SBEST(5,K)
          END DO
          DO J = 5,I+1,-1
            R2BEST(J) = R2BEST(J-1)
            DO K = 1,IP
              SBEST(J,K) = SBEST(J-1,K)
            END DO
          END DO
          R2BEST(I) = RDUM
          DO K = 1,IP
            SBEST(I,K) = SDUM(K)
          END DO
          GO TO 500
        END IF
      END DO
C
C   ###############################
C   STEP 5: WELL FORMULATED TEST
C   ###############################
C
 500  IF(S(IP).GT.H) THEN
        DO L = 1,IP-1
          LL = S(L)
          DO M = 1,2
            IF(V(S(IP),M).EQ.LL) THEN
              CALL SWEEP(P+1,S(IP),R)
              S(IP) = 0
              IP = IP - 1
              IF(IP.EQ.0) GO TO 800
              GO TO 500
            END IF
          END DO
        END DO
      END IF
C
C   ###############################
C   STEP 6: BRANCH RIGHT
C   ###############################
C
 600  IF(S(IP).GE.P) GO TO 700
      CALL SWEEP(P+1,S(IP),R)
      S(IP) = S(IP) + 1
      CALL SWEEP(P+1,S(IP),R)
      GO TO 200
C
C   ###############################
C   STEP 7: FATHOM: DEPTH RETRACTION
C   ###############################
C
 700  CALL SWEEP(P+1,S(IP),R)
      S(IP)=0
      IP = IP-1
      IF(IP.GT.0) GO TO 500
C
 800  DO K = 1,5
        DO I = 1,P
          SSEL(K,I) = 0
        END DO
        DO I = 1,P-Q
          SSEL(K,SBEST(K,I)) = 1
        END DO
      END DO
      DO K = 1,5
        ISEL = 0
        DO I = 1,P
          IF(SSEL(K,I).EQ.0) THEN
            ISEL = ISEL + 1
            SS(K,ISEL) = PERM(I)
          END IF
        END DO
      END DO
      DO J = 1,5
        WRITE(*,152) Q,J,R2BEST(J),(SS(J,I),I=1,Q)
        WRITE(3,152) Q,J,R2BEST(J),(SS(J,I),I=1,Q)
      END DO
 900  CONTINUE
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEB = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
      WRITE(*,154) TIMEB-TIMEA
      WRITE(3,154) TIMEB-TIMEA
 152  FORMAT(' AT Q =',I3,' SOL = ',I2,' RSQUARED IS ',F7.4,55I3)
 154  FORMAT(' TOTAL COMPUTATION TIME IS',F10.2,' SECONDS')
      END
C
      SUBROUTINE SWEEP(N,ICOL,R)
      IMPLICIT INTEGER(A-Z)
      REAL*8 R(100,100),PIVOT,G(100,100)
      PIVOT = 1./R(ICOL,ICOL)
      DO I = 1,N
        DO J = 1,N
          G(I,J) = R(I,J) - PIVOT*R(I,ICOL)*R(ICOL,J)
        END DO
      END DO
      DO J = 1,N
        G(ICOL,J) = R(ICOL,J)*PIVOT
      END DO
      DO I = 1,N
        G(I,ICOL) = -R(I,ICOL)*PIVOT
      END DO
      G(ICOL,ICOL) = PIVOT
      DO I = 1,N
        DO J = 1,N
          R(I,J) = G(I,J)
        END DO
      END DO
      RETURN
      END