C     -----------------------------------------------------------------
C     BARTER ALGORITHM
C     -----------------------------------------------------------------
C     *************** THIS METHOD IS PROPOSED BY SK MISHRA ************
C     PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA)
C     -----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! TYPE DECLARATION
      PARAMETER(IPRINT=500, EPS=1.D-08)
      COMMON /RNDM/IU,IV ! RANDOM NUMBER GENERATION (IU = 4-DIGIT SEED)
      COMMON /XYDATA/X0,Y0,ALIM,NPOINTS
      DIMENSION X0(200),Y0(200),ALIM(2,4) ! IF MORE DATA INCREASE THESE
      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
C     -----------------------------------------------------------------
C     ----------------------- VERY IMPORTANT --------------------------
C     FOR THE BARTER METHOD THE DIMENSION M MUST EXCEED UNITY : OR M =>2
C     ----------------------------------------------------------------
      DIMENSION X(500,10),FV(500),A(50),B(50),BEST(50)
C     ----------------------------------------------------------------
      M=4  ! no. of parameters

C     ----------------------------------------------------------------
      NBART=0 ! NO. OF TIMES  THE BARTER WAS SUCCESSFUL
      NTRY=0  ! NO. OF TIMES  ATTEMPT WAS MADE TO BARTER
      N=M*10
      IF(N.LT.50) N=50
      ITER=10000
      WRITE(*,*)'RANDOM NUMBER SEED ?'
      READ(*,*) IU
C     IU=7113 ! COULD BE ANY OTHER 4 OR 5 DIGIT ODD INTEGER
C     -----------------------------------------------------------------
      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)   ! GENERATE INITIAL X WITHIN
      X(I,J)=(RAND-.5D00)*200   ! GENERATE INITIAL X WITHIN
C     RANDOM NUMBERS BETWEEN -RRANGE AND +RRANGE (BOTH EXCLUSIVE)
      ENDDO
      ENDDO
      IPCOUNT=0
      DO 100 ITR=1,ITER  ! ITERATION BEGINS

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,FA)
      FV(I)=FA
      CALL SEARCH(A,M,FI)
C     STORE FUNCTION VALUES IN FV VECTOR
      IF(FI.LT.FV(I)) THEN
      DO J=1,M
      X(I,J)=A(J)
      ENDDO
      FV(I)=FI
      ENDIF
      ENDDO
C     ----------------------------------------------------------------
      DO I=1,N
      NTRY=NTRY+1
C     CHOOSE IB TH INDIVIDUAL RANDOMLY
2      CALL RANDOM(RAND)
       IB=INT(RAND*N)+1 !THE RANDOM INDIVIDUAL
      IF(IB.EQ.I) GOTO 2
C     IB=2*M+1
C     STORE ITH IN A AND RANDOMLY SELECTED INDIVIDUAL IN B


         DO J=1,M
         A(J)=X(I,J)
         B(J)=X(IB,J) ! OF THE INDIVIDUAL RANDOMLY SELECTED
         ENDDO
C     CHOSE AN INDEX BETWEEN 1 AND M RANDOMLY

          CALL RANDOM(RAND)
          JA=INT(RAND*M)+1
C     CHOOSE ANOTHER INDEX RANDOMLY : MUST BE DIFFERENT FROM JA
   1       CALL RANDOM(RAND)
          JB=INT(RAND*M)+1
          IF(JA.EQ.JB) GOTO 1
C     EXCHANGE A(JA) WITH B(JB)
           TEMP1=A(JA)
           TEMP2=B(JB)
            CALL NORMAL(RN)! OBTAIN STANDARD NORMAL  RANDOM NUMBER
           A(JB)=A(JB)+RN*TEMP2
           B(JB)=B(JB)-RN*TEMP2
           A(JA)=A(JA)-RN*TEMP1
           B(JA)=B(JA)+RN*TEMP1
C     EVALUATE A AND B VECTORS
           CALL FUNC(A,M,FA)
           CALL FUNC(B,M,FB)
C     CHECK IF FA < FV(I) AND FB < FV(IB)
           IF(FA.LT.FV(I) .AND. FB.LT.FV(IB)) THEN
           NBART=NBART+1
           FV(I)=FA
           FV(IB)=FB

           DO J=1,M
           X(I,J)=A(J)
           X(IB,J)=B(J)
           ENDDO
           ENDIF

      ENDDO
C     ----------------------------------------------------------------
C     FIND THE  BEST
      FBEST=1.D30
      DO I=1,N
      IF(FV(I).LT.FBEST) THEN
      FBEST=FV(I)
      KB=I
      ENDIF
      ENDDO
      DO J=1,M
      BEST(J)=X(KB,J)
      ENDDO
C     ----------------------------------------------------------------
      IPCOUNT=IPCOUNT+1
      IF(IPCOUNT.EQ.IPRINT) THEN
      WRITE(*,*)(BEST(J),J=1,M),'  FBEST UPTO NOW = ',FBEST
      WRITE(*,*)'TOTAL NUMBER OF FUNCTION CALLS =',NFCALL
           IF(DABS(FBEST-GBEST).LT.EPS) THEN
           WRITE(*,*) FTIT
           PSUCCESS=NBART/DFLOAT(NTRY)*100.0
      WRITE(*,*) 'NO. OF TIMES AN ATTEPT TO BARTER WAS MADE = ', NTRY
      WRITE(*,*) 'NO. OF TIMES BARTER WAS SUCCESSFUL =' , NBART
      WRITE(*,*) 'PERCENTAGE SUCCESS OF BARTER ATTEMPTS = ', PSUCCESS
      WRITE(*,*)'CX=',BEST(1),' CY=',BEST(2),' A=',BEST(3),' B=',BEST(4)
     * ,' RSQUARE=',-FBEST
      WRITE(*,*)'PROGAM TERMINATED SUCCESSFULLY. THANK YOU.'
      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  -- MAY MAKE IT ZERO  --'
      WRITE(*,*)'RAISE ITER OR DO BOTH OR INCREASE N, POPULATION SIZE'
      PSUCCESS=NBART/DFLOAT(NTRY)*100.0
      WRITE(*,*) 'NO. OF TIMES AN ATTEPT TO BARTER WAS MADE = ', NTRY
      WRITE(*,*) 'NO. OF TIMES BARTER WAS SUCCESSFUL =' , NBART
      WRITE(*,*) 'PERCENTAGE SUCCESS OF BARTER ATTEMPTS = ', PSUCCESS
      WRITE(*,*)'CX=',BEST(1),' CY=',BEST(2),' A=',BEST(3),' B=',BEST(4)
     * ,' RSQUARE=',-FBEST
      WRITE(*,*)'PROGAM TERMINATED SUCCESSFULLY. THANK YOU.'
      END
C     ----------------------------------------------------------------
      SUBROUTINE SEARCH(A,M,FI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /XYDATA/X0,Y0,ALIM,NPOINTS
      DIMENSION X0(200),Y0(200),ALIM(2,4) ! IF MORE DATA INCREASE THESE
      INTEGER IU,IV
      DIMENSION A(*),B(10)
      NSTEP=11
      AMN=1.0D30
      DO J=1,NSTEP
         DO JJ=1,M
         CALL RANDOM(RAND)
         B(JJ)=A(JJ)+(J-(NSTEP/2)-1)*RAND*0.0001D00
         ENDDO
      CALL FUNC(B,M,FI)
        IF(FI.LT.AMN) THEN
        AMN=FI
        DO JJ=1,M
        A(JJ)=B(JJ)
        ENDDO
        ENDIF
      ENDDO
      FI=AMN
      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     -----------------------------------------------------------------

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
       RAND=REAL(RAND1)
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAND=IV
       IU=IV
       RAND=RAND*0.4656613E-09
       RAND1= RAND
       RETURN
       END
C     -----------------------------------------------------------------
      SUBROUTINE FUNC(A,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 A(*)
      CHARACTER *70 FTIT
      CHARACTER *30 INFILE
      COMMON /XYDATA/X0,Y0,ALIM,NPOINTS
      DIMENSION X0(200),Y0(200),ALIM(2,4) ! IF MORE DATA INCREASE THESE
C     -----------------------------------------------------------------
      NFCALL=NFCALL+1 ! INCREMENT TO NUMBER OF FUNCTION CALLS
C     KF IS THE CODE OF THE TEST FUNCTION
      IF(NFCALL.EQ.1) THEN
C     ----------------------------------------------------------------
      FTIT='LOGARITHMIC SPIRAL'  ! TITLE OF THE FUNCTION TO FIT
      KF=1 ! FUNCTION CODE = 1
C      M=4 ! FOUR PARAMETERS CX, CY, A AND B
C     CX AND CY ARE DISPLACEMENT IN X AND Y,
C     A AND B ARE AS IN  SPIRAL R=A*EXP(B*THETA)
C     ---------- GUESS THE LIMITS ON PARAMETERS CX, CY, A AND B --------
      WRITE(*,*)'PROGRAM TO FIT A LOGARITHMIC SPIRAL TO GIVEN DATA'
      WRITE(*,*)'PROGRAM BY PROF. SK MISHRA'
      WRITE(*,*)'DEPARTMENT OF ECONOMICS, NEHU, SHILLONG (INDIA)'
      WRITE(*,*)'[METHOD OF OPTIMIZATION : BARTER ALGORITHM]'
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'PROVIDE LOWER LIMITS ON CX, CY, A AND B'
      READ(*,*)(ALIM(1,J),J=1,M) ! LOWER LIMITS ON PARAMETERS
      WRITE(*,*)'PROVIDE UPPER LIMITS ON CX, CY, A AND B'
      READ(*,*)(ALIM(2,J),J=1,M) ! UPPER LIMITS ON PARAMETERS
C     -----------------------------------------------------------------
C     ------------- READ X Y DATA FROM A SPECIFIED FILE ----------------
      WRITE(*,*)'INPUT FILE NAME ?'
      READ(*,*) INFILE
      OPEN(7,FILE=INFILE)
      READ(7,*) NPOINTS ! NO OF OBSERVATIONS OR POINTS (X,Y) IN DATA
      DO I=1,NPOINTS
      READ(7,*) X0(I),Y0(I)
C     WRITE(*,*) X0(I),Y0(I)
      ENDDO
      CLOSE(7)
      WRITE(*,*)'COMPUTING --- PLEASE WAIT '
      ENDIF
C     -----------------------------------------------------------------
C     CHECK IF PARAMETERS ARE WITHIN LIMITS. IF NOT, BRING THEM WITHIN
      DO J=1,M
      IF(A(J).LT.ALIM(1,J).OR. A(J).GT.ALIM(2,J)) THEN
      CALL RANDOM(RAND)
      A(J)=RAND*(ALIM(2,J)-ALIM(1,J))+ALIM(1,J)
      ENDIF
      ENDDO
      CALL FUNCT(A,M,F)
      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE FUNCT(A,M,F)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      COMMON /XYDATA/X0,Y0,ALIM,NPOINTS
      DIMENSION X0(200),Y0(200),ALIM(2,4) ! IF MORE DATA INCREASE THESE
      DIMENSION A(*),X(200),Y(200),R(200),Q(200),AK(200),ANG(200)
      DIMENSION D(2),C(2,2)
C       WRITE(*,*)'ENTERED IN FUNCT**'
C     EVALUATE THE FUNCTION
      N=NPOINTS
      PI=4.0*DATAN(1.D0)
C     SET X AND Y TO ORIGINAL X0 AND Y0
      DO 12 I=1,N
      X(I)=X0(I)
      Y(I)=Y0(I)
   12 CONTINUE
      DO 1 I=1,N
      X(I)=X(I)-A(1) !SUBSTRACT DISPLACEMENT  CX = A(1)
      Y(I)=Y(I)-A(2) !SUBSTRACT DISPLACEMENT  CY = A(2)
      R(I)=DSQRT(X(I)**2+Y(I)**2)
C      WRITE(*,*) X(I),Y(I),R(I)
C      READ(*,*) AAAAA
    1 CONTINUE
      DO 6 I=1,N
      Q(I)=1
      IF((X(I).LT.0.0).AND.(Y(I).GE.0.0)) Q(I)=2
      IF((X(I).LT.0.0).AND.(Y(I).LT.0.0)) Q(I)=3
      IF((X(I).GE.0.0).AND.(Y(I).LT.0.0)) Q(I)=4
    6 CONTINUE
      DO 7 I=1,N-1
      DO 8 II=I+1,N
      IF(R(I).GT.R(II)) THEN
      TEMP=Q(I)
      Q(I)=Q(II)
      Q(II)=TEMP
      TEMP=X(I)
      X(I)=X(II)
      X(II)=TEMP
      TEMP=Y(I)
      Y(I)=Y(II)
      Y(II)=TEMP
      TEMP=R(I)
      R(I)=R(II)
      R(II)=TEMP
      ENDIF
    8 CONTINUE
    7 CONTINUE
      KK=0
      AK(1)=KK
      DO 5 I=2,N
      AK(I)=KK
      IF(Q(I).LT.Q(I-1)) THEN
      KK=KK+1
      AK(I)=KK
      ENDIF
    5 CONTINUE
      DO 9 I=1,N
      IF(DABS(X(I)).LE.1.0E-30) THEN
      IF(Y(I).GT.0) ANGO=PI/2.0
      IF(Y(I).LT.0) ANGO=3.0*PI/2.0
      ENDIF
      IF(DABS(X(I)).GT.1.0E-30) ANGO=DATAN(Y(I)/X(I))
      ANG(I)=INT(Q(I)/2)*PI+ANGO
C      ANG(I)=ATAN(Y(I)/X(I))
C      IF(Q(I).EQ.2) ANG(I)=PI+ANG(I)
C      IF(Q(I).EQ.3) ANG(I)=PI+ANG(I)
C      IF(Q(I).EQ.4) ANG(I)=2*PI+ANG(I)
      ANG(I)=ANG(I)+2.0*PI*AK(I)
    9 CONTINUE
C     LEAST SQUARES ESTIMATION
      SSR=0
      SR=0
      SANG=0
      SSANG=0
      SRANG=0
      DO 10 I=1,N
      RL=DLOG(R(I))
      SR=SR+RL
      SANG=SANG+ANG(I)
      SSANG=SSANG+ANG(I)**2
      SRANG=SRANG+RL*ANG(I)
      SSR=SSR+RL**2
   10 CONTINUE
      D(1)=SR
      D(2)=SRANG
      C(1,1)=N
      C(1,2)=SANG
      C(2,1)=SANG
      C(2,2)=SSANG
C     MATRIX INVERSION OF C
      TEMP=C(1,1)
      C(1,1)=C(2,2)
      C(2,2)=TEMP
      C(1,2)=-C(1,2)
      C(2,1)=-C(2,1)
      DET=C(1,1)*C(2,2)-C(1,2)*C(2,1)
      C(1,1)=C(1,1)/DET
      C(1,2)=C(1,2)/DET
      C(2,1)=C(2,1)/DET
      C(2,2)=C(2,2)/DET
C     INVERTED MATRIX C-INVERSE IS POST-MULTIPLIED BY D
      AA=C(1,1)*D(1)+C(1,2)*D(2)
      BB=C(2,1)*D(1)+C(2,2)*D(2)
C     ROOT MEAN SQUARES OF ERRORS
      SMS=0.0
      DO 11 I=1,N
      SMS=SMS+(LOG(R(I))-AA-BB*ANG(I))**2
  11  CONTINUE
C      RMS=SQRT(SMS/(N))
      VARR=SSR/N-(SR/N)**2
      RSQUARE=1.0-(SMS/N)/VARR
C      TA=AA/(RMS*SQRT(C(1,1)))
C      TB=BB/(RMS*SQRT(C(2,2)))
C     TAKING ANTILOG(AA)
      AA= EXP(AA)
      F=RSQUARE
      ALPHA=AA
C     ----------------- MINIMIZE -R_SQUARE ----------------------------
      F=-F
C     -----------------------------------------------------------------
      BETA=BB
      A(3)=ALPHA  ! AS IN R=ALPHA*EXP(BETA*THETA))
      A(4)=BETA   ! AS IN R=ALPHA*EXP(BETA*THETA))
      RETURN
      END