C     CONSTRUCTION OF ROBUST COMPOSTE INDICES -------------------------
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.0, RX2=0.9) ! 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     ------------------------------  NOTE -----------------------------
C     [RX1=0,RX2=1] (PURE EXPONENTIAL CROSSOVER) IS BEST IN MOST CASES
C     ------------------------------------------------------------------
C     NOTE:(NCROSS=2) ! CROSS-OVER SCHEME (NCROSS <=0 OR =1 OR =>2)
      PARAMETER(IPRINT=500,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
      CHARACTER *70 FTIT ! TITLE OF THE FUNCTION
      COMMON /RINDEX/VARS(30,6),FINDEX(30),RVAL(6),NSTART
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),IR(3)
C     ----------------------------------------------------------------
      NSTART=0
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 !SEED OF RANDOM NUMBER (4-DIGIT ODD NATURAL NUMBER)

      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
      CALL RANDOM(RAND) ! GENERATES INITION X WITHIN
      X(I,J)=(RAND-.5D00)*2000 ! GENERATES INITION X WITHIN
C     RANDOM NUMBERS BETWEEN -1000 AND 1000 (BOTH EXCLUSIVE)
      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
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)
    1   JR=INT(RAND*M)+1
        J=JR
    2   A(J)=X(IR(1),J)+FACT*(X(IR(2),J)-X(IR(3),J))
    3   J=J+1
        IF(J.GT.M) J=1
    4   IF(J.EQ.JR) GOTO 10
    5   CALL RANDOM(RAND)
        IF(PCROS.LE.RAND) GOTO 2
    6   A(J)=X(I,J)
    7   J=J+1
        IF(J.GT.M) J=1
    8   IF (J.EQ.JR) GOTO 10
    9   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
C     ----------------------------------------------------------------
           WRITE(*,*)'CORRELATION COEFFICIENTS ARE AS FOLLOWS'
           WRITE(*,*)(RVAL(J),J=1,M)

C     -----------------------------------------------------------------
           IF(DABS(FBEST-GBEST).LT.EPS) THEN
           WRITE(*,*) FTIT
           WRITE(*,*)'NO. OF VARIABLES =', M
           WRITE(*,*)'COMPUTATION OVER. THANK YOU'
           IF(KF.EQ.1) OPEN(9,FILE='CINDEXB.TXT')
           IF(KF.GE.2.AND.KF.LE.4) OPEN(9,FILE='CINDEXP.TXT')
           IF(KF.EQ.5) OPEN(9,FILE='CINDEXS.TXT')
           DO I=1,30
           WRITE(9,117)(VARS(I,J),J=1,6),FINDEX(I)
           ENDDO
  117 FORMAT(7F11.5)
           WRITE(*,*)'FINAL CORRELATION COEFFICIENTS ARE AS FOLLOWS'
           WRITE(*,*)(RVAL(J),J=1,M)
           CLOSE(9)


           STOP
           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)'
      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(RAND)
       IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
       COMMON /RNDM/IU,IV
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAND=IV
       IU=IV
       RAND=RAND*0.4656613E-09
       RETURN
       END
C     -----------------------------------------------------------------
      SUBROUTINE FSELECT(KF,M,FTIT)
      PARAMETER(NFUNCT=5)
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
      WRITE(*,*)'----------------------------------------------------'
      DATA TIT(1)/'KF=1 CINDEX (BRADLEY) FUNCTION: 6-VARIABLES 6=?'/
      DATA TIT(2)/'KF=2 CINDEX (PCA) FUNCTION: 6-VARIABLES 6=?'/
      DATA TIT(3)/'KF=3 RANK CINDEX (PCA) FUNCTION: 6-VARIABLES 6=?'/
      DATA TIT(4)/'KF=4 SIGNUM CINDEX (PCA) FUNCTION: 6-VARIABLES 6=?'/
      DATA TIT(5)/'KF=5 SHEVLYAKOV CINDEX FUNCTION: 6-VARIABLES 6=?'/
C     -----------------------------------------------------------------
      DO I=1,NFUNCT
      WRITE(*,*)TIT(I)
      ENDDO
      WRITE(*,*)'----------------------------------------------------'
      WRITE(*,*)'FUNCTION CODE [KF] AND NO. OF VARIABLES [M] ?'
      READ(*,*) KF,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
      PI=4.D+00*DATAN(1.D+00)! DEFINING THE VALUE OF PI
      NFCALL=NFCALL+1 ! INCREMENT TO NUMBER OF FUNCTION CALLS
C     KF IS THE CODE OF THE TEST FUNCTION
C     ------------------------------------------------------------------
      IF(KF.EQ.1) THEN
C     CINDEX (BRADLEY) FUNCTION
      CALL  CINDEX(M,X,F)
      RETURN
      ENDIF
      IF(KF.GE.2.AND.KF.LE.5) THEN
C     CINDEX (PCA) FUNCTION
      CALL  CINDEX(M,X,F)
      RETURN
      ENDIF
C     -----------------------------------------------------------------
C      IF(KF.EQ.99) THEN
C     *****  FUNCTION
C      CALL SUBROUTINE(M,X,F)
C      RETURN
C      ENDIF
C     =================================================================
      WRITE(*,*)'FUNCTION NOT DEFINED. PROGRAM ABORTED'
      STOP
      END
C     >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      SUBROUTINE CINDEX(M,W,F)
      PARAMETER (N=30,MMAX=6,ITYPE=0)!ITYPE=0 MAX SUM, ELSE MAXIMIN
C     IF ITYPE=0 OBTAIN MAX SUM CORREL, ELSE MAXIMIN CORRELATION
      PARAMETER(LP=2) ! ABSOLUTE LP=1; SQUARED LP=2
C     KF=1, LP=1 MAXIMIZES SUM OF ABS BRADLEY CORRELATIONS
C     KF=2, LP=1 MAXIMIZES SUM OF ABS PEARSON CORRELATIONS
C     KF=2, LP=2 MAXIMIZES SUM OF SQUARED PEARSONION CORRELATIONS
C     -----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N,MMAX),Y(N),R(MMAX),W(*)
      COMMON /RNDM/IU,IV
      COMMON /RINDEX/VARS(30,6),FINDEX(30),RVAL(6),NSTART
      COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS * TITLE
      CHARACTER *70 FTIT ! TITLE OF THE FUNCTION
C     CONSTRUCT INDEX Y(N) AS A WEIGHTED SUM OF X(N,M)
      IF(NSTART.EQ.0) THEN
      OPEN(9,FILE='RDAT.TXT')
      DO I=1,N
      READ(9,*)(VARS(I,J),J=1,MMAX)
      ENDDO
      CLOSE(9)
      NSTART=1
      ENDIF
      DO I=1,N
      DO J=1,MMAX
      X(I,J)=VARS(I,J)
      ENDDO
      ENDDO
      SW=0.D0
      DO J=1,M

      IF(DABS(W(J)).GT.1.D0) THEN
      CALL RANDOM(RAND)
      W(J)=(RAND-0.5D0)*2
      ENDIF
      SW=SW+W(J)
      ENDDO
      DO J=1,M
      W(J)=W(J)/SW
      ENDDO
C     ------------------------------------------------------------------
      DO I=1,N
      Y(I)=0.D0
      DO J=1,M
      Y(I)=Y(I)+VARS(I,J)*W(J)
      ENDDO
      ENDDO
      CALL BRADLEY(X,Y,N,M,R)
      FX=0.D0
      RMIN=1.D0
      DO J=1,M
      RVAL(J)=R(J)
      RJ=DABS(R(J))
      FX=FX+RJ**LP ! IF LP=1 ABSOLUTE, IF LP=2 THEN SQUARED (PCA TYPE)
      IF(RJ.LT.RMIN) RMIN=RJ
      ENDDO
      IF(ITYPE.EQ.0) F=-FX ! MAXIMIZES SUM OF ABSOLUTE CORRELATION
      IF(ITYPE.NE.0) F=-RMIN ! MAXIMIZES MINIMAL ABSOLUTE CORRELATION
C     OBTAIN STANDARDIZED INDEX
      YMAX=Y(1)
      YMIN=Y(1)
      DO I=2,N
      IF(Y(I).GT.YMAX) YMAX=Y(I)
      IF(Y(I).LT.YMIN) YMIN=Y(I)
      ENDDO
      DO I=1,N
C      FINDEX(I)=(Y(I)-YMIN)/(YMAX-YMIN)
      FINDEX(I)=Y(I)
      ENDDO
      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE BRADLEY(X,Y,N,M,R) ! -----------------------------------
C     SUBROUTINE BRADLEY:FINDS BRADLEY COEFF. OF ABSOLUTE CORRELATION
C     REFERENCE: BRADLEY, C. (1985) "THE ABSOLUTE CORRELATION",
C     THE MATHEMATICAL GAZETTE, 69 (447), PP. 12-17.
C     X(N,M) ARE THE GIVEN VARIATES IN N OBSERVATIONS; Y(N) IS FIXED
C     R IS VECTOR OF COEF OF ABSOLUTE CORRELATION BETWEEN Y(N) & X(N,M)
      PARAMETER (NMAX=30,MMAX=6)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N,M),Y(N),Z1(NMAX),Z2(NMAX),R(MMAX)
      COMMON /KFF/KF,NFCALL,FTIT ! FUNCTION CODE, NO. OF CALLS * TITLE
      COMMON /RNDM/IU,IV
      CHARACTER *70 FTIT ! TITLE OF THE FUNCTION
C     STORE Y(N) IN Z1()
      DO I=1,N
      Z1(I)=Y(I)
      ENDDO
      IF(KF.EQ.3) THEN
      CALL RANK(Z1,N)
      ENDIF
      IF(KF.EQ.4) THEN
      CALL MEDIAN(Z1,N,A,V)
      DO I=1,N
      IF(DABS(Z1(I)-A).GT.1.0D-05) THEN
      Z1(I)=(Z1(I)-A)/DABS(Z1(I)-A)
      ELSE
      Z1(I)=1.D0
      ENDIF
      ENDDO
      ENDIF
C     IF KF =1 THEN FIND MEDIAN AND MEAN DEVIATION OF Z1
C     IF KF =2 THEN FIND ARITHMETIC AND STANDARD DEVIATION OF Z1
      IF(KF.EQ.1) CALL MEDIAN(Z1,N,A,V)
      IF(KF.EQ.2.OR.KF.EQ.3.OR.KF.EQ.4) CALL MEAN(Z1,N,A,V)
C     STANDARDIZE Z1 TO HAVE ZERO MEAN AND UNIT MEAN DEVIATION
      IF(KF.NE.5) THEN
       DO I=1,N
       Z1(I)=(Z1(I)-A)/V
       ENDDO
      ENDIF

C     CORRELATION BETWEEN Y(N) AND X(N,1), X(N,2), X(N,3) ETC
C     -----------------------------------------------------------------
      DO J=1,M ! VARIABLE 1 TO M
      DO I=1,N
      Z2(I)=X(I,J)
      ENDDO
      IF(KF.EQ.3) THEN
      CALL RANK(Z2,N)
      ENDIF
      IF(KF.EQ.4) THEN
      CALL MEDIAN(Z2,N,A,V)
      DO I=1,N
      IF(DABS(Z2(I)-A).GT.1.0D-05) THEN
      Z2(I)=(Z2(I)-A)/DABS(Z2(I)-A)
      ELSE
      Z2(I)=1.D0
      ENDIF
      ENDDO
      ENDIF
C     IF KF =1 THEN FIND MEDIAN AND MEAN DEVIATION OF Z2
C     IF KF =2 THEN FIND ARITHMETIC MEAN AND STD DEVIATION OF Z2
      IF(KF.EQ.1) CALL MEDIAN(Z2,N,A,V)
      IF(KF.EQ.2.OR.KF.EQ.3.OR.KF.EQ.4) CALL MEAN(Z2,N,A,V)
C     STANDARDIZE Z1 AND Z2 TO HAVE ZERO MEAN AND EQUAL MEAN DEVIATION
       IF(KF.NE.5) THEN
       DO I=1,N
       Z2(I)=(Z2(I)-A)/V
       ENDDO
       ENDIF
      IF(KF.EQ.1) THEN
C     OBTAIN BRADLEY COEFFICIENT OF ABSOLUTE CORRELATION
      S1=0.D0
      S2=0.D0
      DO I=1,N
      S1=S1+DABS(Z1(I)+Z2(I))-DABS(Z1(I)-Z2(I))
      S2=S2+DABS(Z1(I))+DABS(Z2(I))
      ENDDO
      R(J)=S1/S2 ! STORE CORRELATION IN VECTOR R(M), M=1,2,3,...,M
      ENDIF

      IF(KF.EQ.2.OR.KF.EQ.3.OR.KF.EQ.4) THEN
C     OBTAIN PEARSONION PRODUCT MOMENT COEFFICIENT OF CORRELATION
      S1=0.D0
      DO I=1,N
      S1=S1+Z1(I)*Z2(I)
      ENDDO
      R(J)=S1/N ! STORE CORRELATION IN VECTOR R(M), M=1,2,3,...,M
      ENDIF
      IF(KF.EQ.5) THEN
      CALL SHEVLYAKOV(Z1,Z2,N,COR)
      R(J)=COR
      ENDIF
      ENDDO ! THE J LOOP
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE MEDIAN(X,N,A,V) ! ------------------------------------
C     SUBROUTINE MEDIAN : FINDS MEDIAN (A) AND MEAN DEVIATION (V) OF A
C     GIVEN VARIATE, VARIATE X(N)
      PARAMETER (NMAX=30)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NMAX),Z(NMAX)
      COMMON /RNDM/IU,IV
C     STORE X IN Z
      DO I=1,N
      Z(I)=X(I)
      ENDDO
C     ARRANGE Z IN AN ASCENDING ORDER
      DO I=1,N-1
      DO J=I+1,N
      IF(Z(I).GT.Z(J)) THEN ! EXCHANGE
      TEMP=Z(I)
      Z(I)=Z(J)
      Z(J)=TEMP
      ENDIF
      ENDDO
      ENDDO
      K=(N+1)/2 ! K IS OBTAINED AS INT((N+1)/2.0D0)
      A=(Z(K)+Z(N+1-K))/2.D0 ! GIVES MEDIAN FOR ODD AS WELL AS EVEN N
C     FIND MEAN DEVIATION
      V=0.D0
      DO I=1,N
      V=V+DABS(Z(I)-A) ! A IS MEDIAN
      ENDDO
      V=V/N ! V IS MEAN DEVIATION FROM MEDIAN
C     WRITE(*,*)'MEDIAN =',A,'  MEAN DEVIATION =',V
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE MEAN(X,N,A,V) ! ------------------------------------
C     SUBROUTINE MEDIAN : FINDS MEDIAN (A) AND MEAN DEVIATION (V) OF A
C     GIVEN VARIATE, VARIATE X(N)
      PARAMETER (NMAX=30)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NMAX),Z(NMAX)
      COMMON /RNDM/IU,IV
C     STORE X IN Z
      DO I=1,N
      Z(I)=X(I)
      ENDDO
      A=0.D0
      V=0.D0
      DO I=1,N
      A=A+Z(I)
      V=V+Z(I)**2
      ENDDO
      A=A/N  ! A IS THE ARITHMETIC MEAN
      V=DSQRT(V/N-A**2) ! V IS STANDARD DEVIATION
C     WRITE(*,*)'ARITHMETIC MEAN =',A,'  STANDARD DEVIATION =',V
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE RANK(X,N)
      PARAMETER (NMAX=1000)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N), SL(NMAX)
      DO I=1,N
      SL(I)=DFLOAT(I)
      ENDDO
      DO I=1,N-1
      DO II=I+1,N
      IF(X(I).LT.X(II)) THEN
      T=X(I)
      X(I)=X(II)
      X(II)=T
      T=SL(I)
      SL(I)=SL(II)
      SL(II)=T
      ENDIF
      ENDDO
      ENDDO
      DO I=1,N
      X(I)=I
      ENDDO
      DO I=1,N-1
      DO II=I+1,N
      IF(SL(I).GT.SL(II)) THEN
      T=X(I)
      X(I)=X(II)
      X(II)=T
      T=SL(I)
      SL(I)=SL(II)
      SL(II)=T
      ENDIF
      ENDDO
      ENDDO
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE SHEVLYAKOV(Z1,Z2,N,COR)
C     ---------------   SHEVLYAKOV CORRELATION   --------------------
      PARAMETER(NMX=1000, M=2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z1(N),Z2(N),Z(NMX,2),X(NMX,2),XX(NMX),R(2,2)
      DIMENSION U(NMX),V(NMX)
      DO I=1,N
      Z(I,1)=Z1(I)
      Z(I,2)=Z2(I)
      ENDDO

      DO J=1,M
      DO I=1,N
      XX(I)=Z(I,J)
      ENDDO
      CALL MEDIAN(XX,N,AMED,AMDEV)
      DO I=1,N
      XX(I)=DABS(Z(I,J)-AMED)
      ENDDO
      CALL MEDIAN(XX,N,HMED,HMDEV)
      DO I=1,N
      X(I,J)=(Z(I,J)-AMED)/HMED
      ENDDO
      ENDDO
      DO J=1,M
      DO JJ=1,M
      R(J,JJ)=0.D0
      DO I=1,N
      U(I)=DABS(X(I,J)+X(I,JJ))
      V(I)=DABS(X(I,J)-X(I,JJ))
      ENDDO
      CALL MEDIAN(U,N,UMED,UMDEV)
      CALL MEDIAN(V,N,VMED,VMDEV)
      R(J,JJ)=(UMED**2-VMED**2)/(UMED**2+VMED**2)
      ENDDO
      ENDDO
C      DO J=1,M
C      WRITE(*,1)(R(J,JJ),JJ=1,M)
C      ENDDO
C      WRITE(*,*)'----------------------------------------------------'
    1  FORMAT(8F9.5)
      COR=R(1,2)
      RETURN
      END