C     MAIN PROGRAM : PROVIDES TO USE DIFFERENTIAL EVOLUTION METHOD TO
C                    COMPUTE COMPOSITE INDEX INDICES
C     BY MAXIMIZING SUM OF (SQUARES, OR ABSOLUTES, OR MINIMUM) OF
C     CORRELATION OF THE INDEX WITH THE CONSTITUENT VARIABLES. THE MAX
C     SUM OF SQUARES IS THE PRINCIPAL COMPONENT INDEX. IT ALSO PRIVIDES
C     TO OBTAIN MAXIMUM ENTROPY ABSOLUTE CORRELATION INDICES.
C     PRODUCT MOMENT AS WELL AS ABSOLUTE CORRELATION (BRADLEY, 1985) MAY
C     BE USED. PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NORTH-EASTERN
C     HILL UNIVERSITY, SHILLONG (INDIA)
C     -----------------------------------------------------------------
C         ADJUST THE PARAMETERS SUITABLY IN SUBROUTINES DE
C     WHEN THE PROGRAM ASKS FOR PARAMETERS, FEED THEM SUITABLY
C     -----------------------------------------------------------------
      PROGRAM DERPSINDEX
      PARAMETER(NOB=50,MVAR=5)!CHANGE THE PARAMETERS HERE AS NEEDED.
C     -----------------------------------------------------------------
C     NOB=NO. OF CASES AND MVAR=NO. OF VARIABLES
C     TO BE ADJUSTED IN SUBROUTINE CORD(M,X,F) ALSO: STATEMENT 931
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS & TITLE
      CHARACTER *30  METHOD(1)
      CHARACTER *70 FTIT
      CHARACTER *40 INFILE,OUTFILE
      COMMON /CORDAT/CDAT(NOB,MVAR),QIND(NOB),R(MVAR),ENTROPY,NORM,NCOR
      COMMON /XBASE/XBAS
      COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED)
      COMMON /GETRANK/MRNK
      INTEGER IU,IV
      DIMENSION XX(3,50),KKF(3),MM(3),FMINN(3),XBAS(1000,50)
      DIMENSION ZDAT(NOB,MVAR+1),FRANK(NOB),RMAT(MVAR+1,MVAR+1)
      DIMENSION X(50)! X IS THE DECISION VARIABLE X IN F(X) TO MINIMIZE
C     M IS THE DIMENSION OF THE PROBLEM, KF IS TEST FUNCTION CODE AND
C     FMIN IS THE MIN VALUE OF F(X) OBTAINED FROM DE OR RPS
      WRITE(*,*)'====================     WARNING    =============== '
      WRITE(*,*)'ADJUST PARAMETERS IN SUBROUTINES DE IF NEEDED '
      NOPT=1 ! OPTIMIZATION BY DE METHOD
      WRITE(*,*)'=================================================== '
      METHOD(1)=' : DIFFERENTIAL EVOLUTION'
C     INITIALIZATION. THIS XBAS WILL BE USED TO
C     INITIALIZE THE POPULATION.
      WRITE(*,*)' '
      WRITE(*,*)'FEED RANDOM NUMBER SEED,NORM,ENTROPY,NCOR'
      WRITE(*,*)'SEED[ANY 4-DIGIT NUMBER]; NORM[1,2,3]; ENTROPY[0,1]; &
     &NCOR[0,1]'
      WRITE(*,*)' '
      WRITE(*,*)'NORM(1)=ABSOLUTE;NORM(2)=PCA-EUCLIDEAN;NORM(3)=MAXIMIN'
      WRITE(*,*)'ENTROPY(0)=MAXIMIZES NORM;ENTROPY(1)=MAXIMIZES ENTROPY'
      WRITE(*,*)'NCOR(0)=PRODUCT MOMENT; NCOR(1)=ABSOLUTE CORRELATION'
      READ(*,*) IU,NORM,ENTROPY,NCOR
      WRITE(*,*)'WANT RANK SCORE OPTIMIZATION? YES(1); NO(OTHER THAN 1)'
      READ(*,*) MRNK
      WRITE(*,*)'INPUT FILE TO READ DATA:YOUR DATA MUST BE IN THIS FILE'
      WRITE(*,*)'CASES (NOB) IN ROWS ; VARIABLES (MVAR) IN COLUMNS'
      READ(*,*) INFILE
      WRITE(*,*)'SPECIFY THE OUTPUT FILE TO STORE THE RESULTS'
      READ(*,*) OUTFILE
      OPEN(9, FILE=OUTFILE)
      OPEN(7,FILE=INFILE)
      DO I=1,NOB
      READ(7,*),CDA,(CDAT(I,J),J=1,MVAR)
      ENDDO
      CLOSE(7)
      DO I=1,NOB
      DO J=1,MVAR
      ZDAT(I,J+1)=CDAT(I,J)
      ENDDO
      ENDDO
      WRITE(*,*)'DATA HAS BEEN READ. WOULD YOU UNITIZE VARIABLES? [YES=1
     & ELSE NO UNITIZATION]'
      WRITE(*,*)'UNITIZE MEANS TRANSFORMATION FROM X(I,J) TO UNITIZED X'
      WRITE(*,*)'[X(I,J)-MIN(X(.,J))]/[MAX(X(.,J))-MIN(X(.,J))]'
      READ(*,*) NUN
      IF(NUN.EQ.1) THEN
      DO J=1,MVAR
      CMIN=CDAT(1,J)
      CMAX=CDAT(1,J)
      DO I=2,NOB
      IF(CMIN.GT.CDAT(I,J)) CMIN=CDAT(I,J)
      IF(CMAX.LT.CDAT(I,J)) CMAX=CDAT(I,J)
      ENDDO
      DO I=1,NOB
      CDAT(I,J)=(CDAT(I,J)-CMIN)/(CMAX-CMIN)
      ENDDO
      ENDDO
      ENDIF
C     -----------------------------------------------------
      WRITE(*,*)' '
      WRITE(*,*)'FEED RANDOM NUMBER SEED [4-DIGIT ODD INTEGER] TO BEGIN'
      READ(*,*) IU
C     THIS XBAS WILL BE USED AS INITIAL X
      DO I=1,1000
      DO J=1,50
      CALL RANDOM(RAND)
      XBAS(I,J)=RAND ! RANDOM NUMBER BETWEEN (0, 1)
      ENDDO
      ENDDO
      WRITE(*,*)' *****************************************************'
C     ------------------------------------------------------------------
      DO I=1,NOPT
      IF(I.EQ.1) THEN
      WRITE(*,*)'==== WELCOME TO DE/RPS PROGRAM FOR INDEX CONSTRUCTION'
      CALL DE(M,X,FMINDE,Q0,Q1) !CALLS DE AND RETURNS OPTIMAL X AND FMIN
      IF(KF.EQ.1) THEN
      WRITE(9,*)'DIFFERENTIAL EVALUATION OPTIMIZATION RESULTS'
      RSUM1=0.D0
      RSUM2=0.D0
      DO J=1,MVAR
      RSUM1=RSUM1+DABS(R(J))
      RSUM2=RSUM2+DABS(R(J))**2
      ENDDO
      WRITE(9,*)'CORRELATION OF INDEX WITH CONSTITUENT VARIABLES'
      WRITE(9,*)(R(J),J=1,MVAR)
      WRITE(9,*)'SUM OF ABS (R)=',RSUM1,';   SUM OF SQUARE(R)=',RSUM2
      WRITE(9,*)'THE INDEX OR SCORE OF DIFFERENT CASES'
      DO II=1,NOB
      WRITE(9,*)QIND(II)
      FRANK(II)=QIND(II)
      ENDDO
      ENDIF
      FMIN=FMINDE
      ENDIF
C     ------------------------------------------------------------------
      DO J=1,M
      XX(I,J)=X(J)
      ENDDO
      KKF(I)=KF
      MM(I)=M
      FMINN(I)=FMIN
      ENDDO
      WRITE(*,*)' '
      WRITE(*,*)' '
      WRITE(*,*)'---------------------- FINAL RESULTS=================='
      DO I=1,NOPT
      WRITE(*,*)'FUNCT CODE=',KKF(I),'  FMIN=',FMINN(I),' : DIM=',MM(I)
      WRITE(*,*)'OPTIMAL DECISION VARIABLES : ',METHOD(I)
      WRITE(*,*) 'WEIGHTS ARE AS FOLLOWS --------------'
      WRITE(9,*) 'WEIGHTS ARE AS FOLLOWS --------------'
      WRITE(9,*)(XX(I,J),J=1,M)
      WRITE(*,*)(XX(I,J),J=1,M)
      WRITE(*,*)'/////////////////////////////////////////////////////'
      ENDDO
      WRITE(*,*)'OPTIMIZATION PROGRAM ENDED'
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'MEASURE OF EQUALITY/INEQUALITY'
      WRITE(*,*)'DE: BEFORE AND AFTER OPTIMIZATION = ',Q0,Q1
      WRITE(*,*)' '
      WRITE(*,*)'RESULTS STORED IN FILE= ',OUTFILE
      WRITE(*,*)'OPEN BY MSWORD OR EDIT OR ANY OTHER EDITOR'
      WRITE(*,*)' '
      WRITE(*,*)'NOTE:VECTORS OF CORRELATIONS & INDEX(BOTH TOGETHER) ARE
     & IDETERMINATE FOR SIGN &   MAY BE MULTIPLED BY (-1) IF NEEDED'
      WRITE(*,*)'THAT IS IF R(J) IS TRANSFORMED TO -R(J) FOR ALL J THEN
     &THE INDEX(I) TOO IS         TRANSFORMED TO -INDEX(I) FOR ALL I'
      WRITE(9,*)' '
      WRITE(9,*)'NOTE: VECTORS OF CORRELATIONS AND INDEX (BOTH TOGETHER)
     & ARE IDETERMINATE FOR SIGN AND MAY BE MULTIPLED BY (-1) IF NEEDED'
      WRITE(9,*)'THAT IS IF R(J) IS TRANSFORMED TO -R(J) FOR ALL J THEN
     &THE INDEX(I) TOO IS TRANSFORMED TO -INDEX(I) FOR ALL I'
      CALL DORANK(FRANK,NOB)
      DO I=1,NOB
      ZDAT(I,1)=FRANK(I)
      ENDDO
      CALL CORREL(ZDAT,NOB,MVAR+1,RMAT)
      WRITE(9,*)'-------------------- CORRELATION MATRIX --------------'
      WRITE(*,*)'-------------------- CORRELATION MATRIX --------------'
      DO I=1,MVAR+1
      WRITE(9,1)(RMAT(I,J),J=1,MVAR+1)
      WRITE(*,1)(RMAT(I,J),J=1,MVAR+1)
      ENDDO
    1 FORMAT(10F8.4)
      WRITE(9,*)'=================================================== '
      WRITE(*,*)'=================================================== '
      WRITE(9,*)'VARIABLES: 1ST IS THE INDEX AND 2ND THE RANK OF INDEX'
      WRITE(*,*)'VARIABLES: 1ST IS THE INDEX AND 2ND THE RANK OF INDEX'
      WRITE(9,*)'=================================================== '
      WRITE(*,*)'=================================================== '
      DO I=1,NOB
      WRITE(9,2)I,QIND(I),(ZDAT(I,J),J=1,MVAR+1)
      WRITE(*,2)I,QIND(I),(ZDAT(I,J),J=1,MVAR+1)
      ENDDO
    2 FORMAT(I3,F12.4,13F5.0)
      SR2=0.D0
      DO J=2,MVAR+1
      SR2=SR2+RMAT(1,J)**2
      ENDDO
      WRITE(9,*)'SUM OF SQUARE OF CORRELATION R(RANK(INDEX),VAR)=',SR2
      WRITE(*,*)'SUM OF SQUARE OF CORRELATION R(RANK(INDEX),VAR)=',SR2
      CLOSE(9)
      WRITE(*,*) 'THE JOB IS OVER'
      END
C     -----------------------------------------------------------------
      SUBROUTINE DE(M,A,FBEST,G0,G1)
C     PROGRAM: "DIFFERENTIAL EVOLUTION ALGORITHM" OF GLOBAL OPTIMIZATION
C     THIS METHOD WAS PROPOSED BY R. STORN AND K. PRICE IN 1995. REF --
C     "DIFFERENTIAL EVOLUTION - A SIMPLE AND EFFICIENT ADAPTIVE SCHEME
C     FOR GLOBAL OPTIMIZATION OVER CONTINUOUS SPACES" : TECHNICAL REPORT
C     INTERNATIONAL COMPUTER SCIENCE INSTITUTE, BERKLEY, 1995.
C     PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA)
C     -----------------------------------------------------------------
C     PROGRAM DE
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! TYPE DECLARATION
      PARAMETER(NMAX=500,MMAX=50) ! MAXIMUM DIMENSION PARAMETERS
      PARAMETER (RX1=0.D0, RX2=0.D0) ! TO BE ADJUSTED SUITABLY, IF NEEDED
C     RX1 AND RX2 CONTROL THE SCHEME OF CROSSOVER. (0 <= RX1 <= RX2) <=1
C     RX1 DETERMINES THE UPPER LIMIT OF SCHEME 1 (AND LOWER LIMIT OF
C     SCHEME 2; RX2 IS THE UPPER LIMIT OF SCHEME 2 AND LOWER LIMIT OF
C     SCHEME 3. THUS RX1 = .2 AND RX2 = .8 MEANS 0-20% SCHEME1, 20 TO 80
C     PERCENT SCHEME 2 AND THE REST (80 TO 100 %) SCHEME 3.
C     PARAMETER(NCROSS=2) ! CROSS-OVER SCHEME (NCROSS <=0 OR =1 OR =>2)
      PARAMETER(IPRINT=100,EPS=1.D-08)!FOR WATCHING INTERMEDIATE RESULTS
C     IT PRINTS THE INTERMEDIATE RESULTS AFTER EACH IPRINT ITERATION AND
C     EPS DETERMINES ACCURACY FOR TERMINATION. IF EPS= 0, ALL ITERATIONS
C     WOULD BE UNDERGONE EVEN IF NO IMPROVEMENT IN RESULTS IS THERE.
C     ULTIMATELY "DID NOT CONVERGE" IS REOPORTED.
      COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED)
      INTEGER IU,IV      ! FOR RANDOM NUMBER GENERATION
      COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS * TITLE
      COMMON /XBASE/XBAS
      CHARACTER *70 FTIT ! TITLE OF THE FUNCTION
C     -----------------------------------------------------------------
C     THE PROGRAM REQUIRES INPUTS FROM THE USER ON THE FOLLOWING ------
C     (1) FUNCTION CODE (KF), (2) NO. OF VARIABLES IN THE FUNCTION (M);
C     (3) N=POPULATION SIZE (SUGGESTED 10 TIMES OF NO. OF VARIABLES, M,
C         FOR SMALLER PROBLEMS N=100 WORKS VERY WELL);
C     (4) PCROS = PROB. OF CROSS-OVER (SUGGESTED : ABOUT 0.85 TO .99);
C     (5) FACT = SCALE (SUGGESTED 0.5 TO .95 OR  1, ETC);
C     (6) ITER = MAXIMUM NUMBER OF ITERATIONS PERMITTED (5000 OR MORE)
C     (7) RANDOM NUMBER SEED (4 DIGITS INTEGER)
C     ----------------------------------------------------------------
      DIMENSION X(NMAX,MMAX),Y(NMAX,MMAX),A(MMAX),FV(NMAX)
      DIMENSION IR(3),XBAS(1000,50)
C     ----------------------------------------------------------------
C     ------- SELECT THE FUNCTION TO MINIMIZE AND ITS DIMENSION -------
      CALL FSELECT(KF,M,FTIT)
C     SPECIFY OTHER PARAMETERS ---------------------------------------
      WRITE(*,*)'POPULATION SIZE [N] AND NO. OF ITERATIONS [ITER] ?'
      WRITE(*,*)'SUGGESTED : N => 100 OR =>10.M; ITER 10000 OR SO'
      READ(*,*) N,ITER
      WRITE(*,*)'CROSSOVER PROBABILITY [PCROS] AND SCALE [FACT] ?'
      WRITE(*,*)'SUGGESTED : PCROS ABOUT 0.9; FACT=.5 OR LARGER BUT <=1'
      READ(*,*) PCROS,FACT
      WRITE(*,*)'RANDOM NUMBER SEED ?'
      WRITE(*,*)'A FOUR-DIGIT POSITIVE ODD INTEGER, SAY, 1171'
      READ(*,*) IU
      NFCALL=0 ! INITIALIZE COUNTER FOR FUNCTION CALLS
      GBEST=1.D30 ! TO BE USED FOR TERMINATION CRITERION
C     INITIALIZATION : GENERATE X(N,M) RANDOMLY
      DO I=1,N
      DO J=1,M
C      CALL RANDOM(RAND) ! GENERATES INITION X WITHIN
C      X(I,J)=(RAND-.5D00)*2000 ! GENERATES INITION X WITHIN
C     RANDOM NUMBERS BETWEEN -RRANGE AND +RRANGE (BOTH EXCLUSIVE)
      X(I,J)=XBAS(I,J)! TAKES THESE NUMBERS FROM THE MAIN PROGRAM
      ENDDO
      ENDDO
      WRITE(*,*)'COMPUTING --- PLEASE WAIT '
      IPCOUNT=0
      DO 100 ITR=1,ITER  ! ITERATION BEGINS
C     -----------------------------------------------------------------
C     EVALUATE ALL X FOR THE GIVEN FUNCTION
      DO I=1,N
      DO J=1,M
      A(J)=X(I,J)
      ENDDO
      CALL FUNC(A,M,F)
C     STORE FUNCTION VALUES IN FV VECTOR
      FV(I)=F
      ENDDO
      IF(ITR.EQ.1) CALL GINI(FV,N,G0)
C     ----------------------------------------------------------------
C     FIND THE FITTEST (BEST) INDIVIDUAL AT THIS ITERATION
                FBEST=FV(1)
                KB=1
                DO IB=2,N
                     IF(FV(IB).LT.FBEST) THEN
                     FBEST=FV(IB)
                     KB=IB
                     ENDIF
                ENDDO
C     BEST FITNESS VALUE = FBEST : INDIVIDUAL X(KB)
C     -----------------------------------------------------------------
C     GENERATE OFFSPRINGS
      DO I=1,N    ! I LOOP BEGINS
C     INITIALIZE CHILDREN IDENTICAL TO PARENTS; THEY WILL CHANGE LATER
           DO J=1,M
           Y(I,J)=X(I,J)
           ENDDO
C     SELECT RANDOMLY THREE OTHER INDIVIDUALS
   20     DO IRI=1,3  ! IRI LOOP BEGINS
          IR(IRI)=0

          CALL RANDOM(RAND)
           IRJ=INT(RAND*N)+1
C     CHECK THAT THESE THREE INDIVIDUALS ARE DISTICT AND OTHER THAN I
           IF(IRI.EQ.1.AND.IRJ.NE.I) THEN
           IR(IRI)=IRJ
           ENDIF
           IF(IRI.EQ.2.AND.IRJ.NE.I.AND.IRJ.NE.IR(1)) THEN
           IR(IRI)=IRJ
           ENDIF
       IF(IRI.EQ.3.AND.IRJ.NE.I.AND.IRJ.NE.IR(1).AND.IRJ.NE.IR(2)) THEN
            IR(IRI)=IRJ
            ENDIF
            ENDDO   ! IRI LOOP ENDS
C     CHECK IF ALL THE THREE IR ARE POSITIVE (INTEGERS)
          DO IX=1,3
          IF(IR(IX).LE.0) THEN
          GOTO 20  ! IF NOT THEN REGENERATE
          ENDIF
          ENDDO
C     THREE RANDOMLY CHOSEN INDIVIDUALS DIFFERENT FROM I AND DIFFERENT
C     FROM EACH OTHER ARE IR(1),IR(2) AND IR(3)
C     ===================== RANDOMIZATION OF NCROSS ===================
C     RANDOMIZES NCROSS
      NCROSS=0
      CALL RANDOM(RAND)
      IF(RAND.GT.RX1) NCROSS=1 ! IF RX1=>1, SCHEME 2 NEVER IMPLEMENTED
      IF(RAND.GT.RX2) NCROSS=2 ! IF RX2=>1, SCHEME 3 NEVER IMPLEMENTED

C     ---------------------- SCHEME 1 ----------------------------------
C      NO CROSS OVER, ONLY REPLACEMENT THAT IS PROBABILISTIC
          IF(NCROSS.LE.0) THEN
          DO J=1,M      ! J LOOP BEGINS
          CALL RANDOM(RAND)
          IF(RAND.LE.PCROS) THEN ! REPLACE IF RAND < PCROS
          A(J)=X(IR(1),J)+(X(IR(2),J)-X(IR(3),J))*FACT ! CANDIDATE CHILD
          ENDIF
          ENDDO  ! J LOOP ENDS
          ENDIF
C     ----------------------- SCHEME 2 ---------------------------------
C     THE STANDARD CROSSOVER SCHEME
C     CROSSOVER SCHEME (EXPONENTIAL) SUGGESTED BY KENNETH PRICE IN HIS
C     PERSONAL LETTER TO THE AUTHOR (DATED SEPTEMBER 29, 2006)
      IF(NCROSS.EQ.1) THEN
        CALL RANDOM(RAND)
        JR=INT(RAND*M)+1
        J=JR
    2   A(J)=X(IR(1),J)+FACT*(X(IR(2),J)-X(IR(3),J))
        J=J+1
        IF(J.GT.M) J=1
        IF(J.EQ.JR) GOTO 10
        CALL RANDOM(RAND)
        IF(PCROS.LE.RAND) GOTO 2
    6   A(J)=X(I,J)
        J=J+1
        IF(J.GT.M) J=1
        IF (J.EQ.JR) GOTO 10
        GOTO 6
   10   CONTINUE
      ENDIF
C     ------------------------ SCHEME 3 --------------------------------
C     ESPECIALLY SUITABLE TO NON-DECOMPOSABLE (NON-SEPERABLE) FUNCTIONS
C     CROSSOVER SCHEME (NEW) SUGGESTED BY KENNETH PRICE IN HIS
C     PERSONAL LETTER TO THE AUTHOR (DATED OCTOBER 18, 2006)
      IF(NCROSS.GE.2) THEN
          CALL RANDOM(RAND)
          IF(RAND.LE.PCROS) THEN
             CALL NORMAL(RN)
             DO J=1,M
             A(J)=X(I,J)+(X(IR(1),J)+ X(IR(2),J)-2*X(I,J))*RN
             ENDDO
           ELSE
            DO J=1,M
            A(J)=X(I,J)+(X(IR(1),J)- X(IR(2),J))! FACT ASSUMED TO BE 1
            ENDDO
           ENDIF
      ENDIF
C     -----------------------------------------------------------------
          CALL FUNC(A,M,F) ! EVALUATE THE OFFSPRING
          IF(F.LT.FV(I)) THEN ! IF BETTER, REPLACE PARENTS BY THE CHILD
          FV(I)=F
          DO J=1,M
          Y(I,J)=A(J)
          ENDDO
          ENDIF
      ENDDO   ! I LOOP ENDS
      DO I=1,N
      DO J=1,M
      X(I,J)=Y(I,J) ! NEW GENERATION IS A MIX OF BETTER PARENTS AND
C                     BETTER CHILDREN
      ENDDO
      ENDDO
      IPCOUNT=IPCOUNT+1
      IF(IPCOUNT.EQ.IPRINT) THEN
      DO J=1,M
      A(J)=X(KB,J)
      ENDDO
      WRITE(*,*)(X(KB,J),J=1,M),'  FBEST UPTO NOW = ',FBEST
      WRITE(*,*)'TOTAL NUMBER OF FUNCTION CALLS =',NFCALL
           IF(DABS(FBEST-GBEST).LT.EPS) THEN
           WRITE(*,*) FTIT
           WRITE(*,*)'COMPUTATION OVER'
           GOTO 999
           ELSE
           GBEST=FBEST
           ENDIF
      IPCOUNT=0
      ENDIF
C     ----------------------------------------------------------------
  100 ENDDO   ! ITERATION ENDS : GO FOR NEXT ITERATION, IF APPLICABLE
C     ----------------------------------------------------------------
      WRITE(*,*)'DID NOT CONVERGE. REDUCE EPS OR RAISE ITER OR DO BOTH'
      WRITE(*,*)'INCREASE N, PCROS, OR SCALE FACTOR (FACT)'
  999 CALL FUNC(A,M,FBEST)
      CALL GINI(FV,N,G1)
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE NORMAL(R)
C     PROGRAM TO GENERATE N(0,1) FROM RECTANGULAR RANDOM NUMBERS
C     IT USES BOX-MULLER VARIATE TRANSFORMATION FOR THIS PURPOSE.
C     -----------------------------------------------------------------
C     ----- BOX-MULLER METHOD BY GEP BOX AND ME MULLER (1958) ---------
C     BOX, G. E. P. AND MULLER, M. E. "A NOTE ON THE GENERATION OF
C     RANDOM NORMAL DEVIATES." ANN. MATH. STAT. 29, 610-611, 1958.
C     IF U1 AND U2 ARE UNIFORMLY DISTRIBUTED RANDOM NUMBERS (0,1),
C     THEN X=[(-2*LN(U1))**.5]*(COS(2*PI*U2) IS N(0,1)
C     ALSO,  X=[(-2*LN(U1))**.5]*(SIN(2*PI*U2) IS N(0,1)
C     PI = 4*ARCTAN(1.0)= 3.1415926535897932384626433832795
C     2*PI = 6.283185307179586476925286766559
C     -----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      INTEGER IU,IV
C     -----------------------------------------------------------------
      CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1]
      U1=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1]
      CALL RANDOM(RAND) ! INVOKES RANDOM TO GENERATE UNIFORM RAND [0, 1]
      U2=RAND ! U1 IS UNIFORMLY DISTRIBUTED [0, 1]
      R=DSQRT(-2.D0*DLOG(U1))
      R=R*DCOS(U2*6.283185307179586476925286766559D00)
C     R=R*DCOS(U2*6.28318530718D00)
      RETURN
      END
C     -----------------------------------------------------------------
C     RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE)
      SUBROUTINE RANDOM(RAND1)
       DOUBLE PRECISION  RAND1
       COMMON /RNDM/IU,IV
      INTEGER IU,IV
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAND=IV
       IU=IV
       RAND=RAND*0.4656613E-09
       RAND1= DBLE(RAND)
       RETURN
       END
C     -----------------------------------------------------------------
      SUBROUTINE GINI(F,N,G)
      PARAMETER (K=1) !K=1 GINI COEFFICENT; K=2 COEFFICIENT OF VARIATION
C     THIS PROGRAM COMPUTES MEASURE OF INEQUALITY
C     IF K =1 GET THE GINI COEFFICIENT. IF K=2 GET COEFF OF VARIATIONE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(*)
      S=0.D0
      DO I=1,N
      S=S+F(I)
      ENDDO
      S=S/N
      H=0.D00
      DO I=1,N-1
      DO J=I+1,N
      H=H+(DABS(F(I)-F(J)))**K
      ENDDO
      ENDDO
      H=(H/(N**2))**(1.D0/K)! FOR K=1 H IS MEAN DEVIATION;
C                             FOR K=2 H IS STANDARD DEVIATION
      WRITE(*,*)'MEASURES OF DISPERSION AND CENTRAL TENDENCY = ',G,S
      G=DEXP(-H)! G IS THE MEASURE OF EQUALITY (NOT GINI OR CV)
C     G=H/DABS(S) !IF S NOT ZERO, K=1 THEN G=GINI, K=2 G=COEFF VARIATION
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE FSELECT(KF,M,FTIT)
C     THE PROGRAM REQUIRES INPUTS FROM THE USER ON THE FOLLOWING ------
C     (1) FUNCTION CODE (KF), (2) NO. OF VARIABLES IN THE FUNCTION (M);
      CHARACTER *70 TIT(100),FTIT
      NFN=1
      KF=1
      WRITE(*,*)'----------------------------------------------------'
      DATA TIT(1)/'CONSTRUCTION OF INDEX FROM M VARIABLES '/
C     -----------------------------------------------------------------
      DO I=1,NFN
      WRITE(*,*)TIT(I)
      ENDDO
      WRITE(*,*)'----------------------------------------------------'
      WRITE(*,*)'SPECIFY NO. OF VARIABLES [MVAR] HERE ALSO ?'
      READ(*,*) M
      FTIT=TIT(KF) ! STORE THE NAME OF THE CHOSEN FUNCTION IN FTIT
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE FUNC(X,M,F)
C     TEST FUNCTIONS FOR GLOBAL OPTIMIZATION PROGRAM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /KFF/KF,NFCALL,FTIT
      INTEGER IU,IV
      DIMENSION X(*)
      CHARACTER *70 FTIT
      NFCALL=NFCALL+1 ! INCREMENT TO NUMBER OF FUNCTION CALLS
C     KF IS THE CODE OF THE TEST FUNCTION
      IF(KF.EQ.1) THEN
      CALL CORD(M,X,F)
      RETURN
      ENDIF
C     =================================================================
      WRITE(*,*)'FUNCTION NOT DEFINED. PROGRAM ABORTED'
      STOP
      END
C     ------------------------------------------------------------------
      SUBROUTINE CORD(M,X,F)
      PARAMETER (NOB=50,MVAR=5)! CHANGE THE PARAMETERS HERE AS NEEDED.
C     ------------------------------------------------------------------
C     NOB=NO. OF OBSERVATIONS (CASES) & MVAR= NO. OF VARIABLES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /CORDAT/CDAT(NOB,MVAR),QIND(NOB),R(MVAR),ENTROPY,NORM,NCOR
      COMMON /GETRANK/MRNK
      INTEGER IU,IV
      DIMENSION X(*),Z(NOB,2)
      DO I=1,M
      IF(X(I).LT.-1.0D0.OR.X(I).GT.1.0D0) THEN
      CALL RANDOM(RAND)
      X(I)=(RAND-0.5D0)*2
      ENDIF
      ENDDO
      XNORM=0.D0
      DO J=1,M
      XNORM=XNORM+X(J)**2
      ENDDO
      XNORM=DSQRT(XNORM)
      DO J=1,NOB
      X(J)=X(J)/XNORM
      ENDDO
C     CONSTRUCT INDEX
      DO I=1,NOB
      QIND(I)=0.D0
      DO J=1,M
      QIND(I)=QIND(I)+CDAT(I,J)*X(J)
      ENDDO
      ENDDO
C     ------------------------------------------------------------------
      !FIND THE RANK OF QIND
      IF(MRNK.EQ.1) CALL DORANK(QIND,NOB)
C     ------------------------------------------------------------------
C     COMPUTE CORRELATIONS
      DO I=1,NOB
      Z(I,1)=QIND(I)
      ENDDO
      DO J=1,M
      DO I=1,NOB
      Z(I,2)=CDAT(I,J)
      ENDDO
      IF(NCOR.EQ.0) THEN
      CALL CORLN(Z,NOB,RHO)
      ELSE
      CALL CORA(Z,NOB,RHO)
      ENDIF
      R(J)=RHO
      ENDDO
      IF(ENTROPY.EQ.0.D0) THEN
C     ------------------ MAXIMIN SOLUTION ----------------------------
      IF(NORM.GT.2) THEN
       F=DABS(R(1))
       DO J=2,M
       IF(F.GT.DABS(R(J))) F= DABS(R(J))
       ENDDO
       ENDIF
C     ------------------ FOR  NORM =1 OR 2 ---------------------------
      IF(NORM.LE.2) THEN
      F=0.D0
      DO J=1,M
      F=F+DABS(R(J))**NORM
      ENDDO
      ENDIF
C     --------------------------------------------------------------
      ELSE
      IF(ENTROPY.NE.0.D0) THEN
C     ENTROPY MAXIMIZATION
      ENT=0.0D0
      DO J=1,M
      ENT=ENT+DABS(R(J))
      ENDDO
      F=ENT*DLOG(ENT)
      DO J=1,M
      FX=DABS(R(J))
      F=F+ (FX/ENT)*DLOG(FX/ENT)
      ENDDO
      ENDIF
      ENDIF
C     ------------------------------------------------------------------
      F=-F
      RETURN
      END
      SUBROUTINE CORLN(Z,NOB,RHO)
C     NOB = NO. OF CASES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(NOB,2),AV(2),SD(2)
      DO J=1,2
      AV(J)=0.D0
      SD(J)=0.D0
      DO I=1,NOB
      AV(J)=AV(J)+Z(I,J)
      SD(J)=SD(J)+Z(I,J)**2
      ENDDO
      ENDDO
      DO J=1,2
      AV(J)=AV(J)/NOB
      SD(J)=DSQRT(SD(J)/NOB-AV(J)**2)
      ENDDO
C      WRITE(*,*)'AV AND SD ', AV(1),AV(2),SD(1),SD(2)
      RHO=0.D0
      DO I=1,NOB
      RHO=RHO+(Z(I,1)-AV(1))*(Z(I,2)-AV(2))
      ENDDO
      RHO=(RHO/NOB)/(SD(1)*SD(2))
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE CORA(Z,N,R)
C     COMPUTING ABSOLUTE CORRELATION MATRIX
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(N,2),X(N),Y(N)
C     ------------------------------------------------------------------
C     PUT Z INTO X AND Y
      DO I=1,N
      X(I)=Z(I,1)
      Y(I)=Z(I,2)
      ENDDO
C     ARRANGE X ANY IN AN ASCENDING ORDER
      DO I=1,N-1
      DO II=I+1,N
      IF(X(I).GT.X(II)) THEN
      TEMP=X(I)
      X(I)=X(II)
      X(II)=TEMP
      ENDIF
      IF(Y(I).GT.Y(II)) THEN
      TEMP=Y(I)
      Y(I)=Y(II)
      Y(II)=TEMP
      ENDIF
      ENDDO
      ENDDO
C     FIND MEDIAN
      IF(INT(N/2).EQ.N/2.D0) THEN
      XMED=(X(N/2)+X(N/2+1))/2.D0
      YMED=(Y(N/2)+Y(N/2+1))/2.D0
      ENDIF
      IF(INT(N/2).NE.N/2.D0) THEN
      XMED=X(N/2+1)
      YMED=Y(N/2+1)
      ENDIF
C     SUBTRACT RESPECTIVE MEDIANS FROM X AND Y AND FIND ABS DEVIATIONS
      VX=0.D0
      VY=0.D0
      DO I=1,N
      X(I)=X(I)-XMED
      Y(I)=Y(I)-YMED
      VX=VX+DABS(X(I))
      VY=VY+DABS(Y(I))
      ENDDO
C     SCALE THE VARIABLES X AND Y SUCH THAT VX=VY
      IF(VX.EQ.0.D0.OR.VY.EQ.0.D0) THEN
      R=0.D0
      RETURN
      ENDIF
      DO I=1,N
      X(I)=X(I)*VY/VX
      ENDDO
C     COMPUTE CORRELATION COEFFICIENT
      VZ=0.D0
      R=0.D0
      DO I=1,N
      VZ=VZ+DABS(X(I))+DABS(Y(I))
      R=R+DABS(X(I)+Y(I))-DABS(X(I)-Y(I))
      ENDDO
      R=R/VZ
      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE DORANK(X,N)! N IS THE NUMBER OF OBSERVATIONS
      PARAMETER (MXD=1000)! MXD IS MAX DIMENSION FOR TEMPORARY VARIABLES
      ! THAT ARE LOCAL AND DO NOT GO TO THE INVOKING PROGRAM
      ! X IS THE VARIABLE TO BE SUBSTITUTED BY ITS RANK VALUES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N),ID(MXD)
      ! GENERATE ID
      DO I=1,N
      ID(I)=I
      ENDDO

      DO I=1,N-1
      DO II=I+1,N
      IF(X(I).GT.X(II)) THEN ! ARRANGE ACCORDING TO ASCENDING ORDER OF X
      T=X(I)
      X(I)=X(II)
      X(II)=T
      IT=ID(I)
      ID(I)=ID(II)
      ID(II)=IT
      ENDIF
      ENDDO
      ENDDO

      DO I=1,N
      X(I)=I+0.D0
      ENDDO

      DO I=1,N-1
      DO II=I+1,N
      IF(ID(I).GT.ID(II)) THEN!ARRANGE ACCORDING TO ASCENDNG ORDER OF ID
      T=X(I)
      X(I)=X(II)
      X(II)=T
      IT=ID(I)
      ID(I)=ID(II)
      ID(II)=IT
      ENDIF
      ENDDO
      ENDDO
      RETURN
      END
C     ----------------------------------------------------------------
      SUBROUTINE CORREL(X,N,M,RMAT)
      PARAMETER (NMX=30)!DO NOT CHANGE UNLESS NO. OF VARIABLES EXCEED 30
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N,M),RMAT(M,M),AV(NMX),SD(NMX)
      DO J=1,M
      AV(J)=0.D0
      SD(J)=0.D0
      DO I=1,N
      AV(J)=AV(J)+X(I,J)
      SD(J)=SD(J)+X(I,J)**2
      ENDDO
      AV(J)=AV(J)/N
      SD(J)=DSQRT(SD(J)/N-AV(J)**2)
      ENDDO
      DO J=1,M
      DO JJ=1,M
      RMAT(J,JJ)=0.D0
      DO I=1,N
      RMAT(J,JJ)=RMAT(J,JJ)+X(I,J)*X(I,JJ)
      ENDDO
      ENDDO
      ENDDO
      DO J=1,M
      DO JJ=1,M
      RMAT(J,JJ)=RMAT(J,JJ)/N-AV(J)*AV(JJ)
      RMAT(J,JJ)=RMAT(J,JJ)/(SD(J)*SD(JJ))
      ENDDO
      ENDDO
      RETURN
      END