C     =================================================================
C     PROGRAM TO FIND GLOBAL MINIMUM BY REPULSIVE PARTICLE SWARM METHOD
C     WRITTEN BY SK MISHRA, DEPT. OF ECONOMICS, NEHU, SHILLONG (INDIA)
C     -----------------------------------------------------------------
      PARAMETER (N=30,NN=15,MX=50,NSTEP=11,ITRN=10000,NSIGMA=1,ITOP=3)
C     PARAMETER(N=50,NN=25,MX=100,NSTEP=9,ITRN=10000,NSIGMA=1,ITOP=3)
C     PARAMETER (N=100,NN=15,MX=100,NSTEP=9,ITRN=10000,NSIGMA=1,ITOP=3)
C     IN CERTAIN CASES THE ONE OR THE OTHER SPECIFICATION WORKS BETTER
C     DIFFERENT SPECIFICATIONS OF PARAMETERS MAY SUIT DIFFERENT TYPES
C     OF FUNCTIONS OR DIMENSIONS - ONE HAS TO DO SOME TRIAL AND ERROR
C     -----------------------------------------------------------------
C     N = POPULATION SIZE. IN MOST OF THE CASES N=30 IS OK. ITS VALUE
C     MAY BE INCREASED TO 50 OR 100 TOO. THE PARAMETER NN IS THE SIZE OF
C     RANDOMLY CHOSEN NEIGHBOURS. 15 TO 25 (BUT SUFFICIENTLY LESS THAN
C     N) IS A GOOD CHOICE. MX IS THE MAXIMAL SIZE OF DECISION VARIABLES.
C     IN F(X1, X2,...,XM) M SHOULD BE LESS THAN OR EQUAL TO MX. ITRN IS
C     THE NO. OF ITERATIONS. IT MAY DEPEND ON THE PROBLEM. 200(AT LEAST)
C     TO 500 ITERATIONS MAY BE GOOD ENOUGH. BUT FOR FUNCTIONS LIKE
C     ROSENBROCKOR GRIEWANK OF LARGE SIZE (SAY M=30) IT IS NEEDED THAT
C     ITRN IS LARGE, SAY 5000 OR EVEN 10000.
C     SIGMA INTRODUCES PERTURBATION & HELPS THE SEARCH JUMP OUT OF LOCAL
C     OPTIMA. FOR EXAMPLE : RASTRIGIN FUNCTION OF DMENSION 3O OR LARGER
C     NSTEP DOES LOCAL SEARCH BY TUNNELLING AND WORKS WELL BETWEEN 5 AND
C     15, WHICH IS MUCH ON THE HIGHER SIDE.
C     ITOP <=1 (RING); ITOP=2 (RING AND RANDOM); ITOP=>3 (RANDOM)
C     NSIGMA=0 (NO CHAOTIC PERTURBATION);NSIGMA=1 (CHAOTIC PERTURBATION)
C     NOTE THAT NSIGMA=1 NEED NOT ALWAYS WORK BETTER (OR WORSE)
C     SUBROUTINE FUNC( ) DEFINES OR CALLS THE FUNCTION TO BE OPTIMIZED.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      COMMON /KFF/KF,NFCALL,FTIT
      INTEGER IU,IV
      CHARACTER *70 FTIT
      COMMON /XYDATA/X0,Y0,ALIM,NPOINTS
      DIMENSION X0(200),Y0(200),ALIM(2,4) ! IF MORE DATA INCREASE THESE
      DIMENSION X(N,MX),V(N,MX),A(MX),VI(MX)
      DIMENSION XX(N,MX),F(N),V1(MX),V2(MX),V3(MX),V4(MX),BST(MX)
C     A1 A2 AND A3 ARE CONSTANTS AND W IS THE INERTIA WEIGHT.
C     OCCASIONALLY, TINKERING WITH THESE VALUES, ESPECIALLY A3, MAY BE
C     NEEDED.
      DATA A1,A2,A3,W,SIGMA,EPSI /.5D0,.5D0,5.D-04,.5D00,1.D-03,1.D-08/
C     -----------------------------------------------------------------
      M=4 ! NO OF PARAMETERS

C     -----------------------------------------------------------------
      GGBEST=1.D30 ! TO BE USED FOR TERMINATION CRITERION
      LCOUNT=0
      NFCALL=0

      WRITE(*,*)'PROVIDE A FOUR-DIGIT POSITIVE ODD INTEGER, SAY, 1171'
      READ(*,*) IU ! COULD BE ANY OTHER 4 OR 5 DIGIT ODD INTEGER
      FMIN=1.0E30
C     GENERATE N-SIZE POPULATION OF M-TUPLE PARAMETERS X(I,J) RANDOMLY
      DO I=1,N
        DO J=1,M
        CALL RANDOM(RAND)  !GENERATE X WITHIN
        X(I,J)=(RAND-0.5D00)*200   ! GENERATE X WITHIN
C     WE GENERATE RANDOM(-1000,1000). HERE MULTIPLIER IS 200.
C     TINKERING IN SOME CASES MAY BE NEEDED
         ENDDO
        F(I)=1.0D30
      ENDDO
C     INITIALISE VELOCITIES V(I) FOR EACH INDIVIDUAL IN THE POPULATION
      DO I=1,N
      DO J=1,M
      CALL RANDOM(RAND)
       V(I,J)=(RAND-0.5D+00)
C       V(I,J)=RAND
      ENDDO
      ENDDO
      DO 100 ITER=1,ITRN
C     LET EACH INDIVIDUAL SEARCH FOR THE BEST IN ITS NEIGHBOURHOOD
        DO I=1,N
           DO J=1,M
           A(J)=X(I,J)
           VI(J)=V(I,J)
           ENDDO
           CALL LSRCH(A,M,VI,NSTEP,FI)
           IF(FI.LT.F(I)) THEN
            F(I)=FI
            DO IN=1,M
            BST(IN)=A(IN)
            ENDDO
C     F(I) CONTAINS THE LOCAL BEST VALUE OF FUNCTION FOR ITH INDIVIDUAL
C     XX(I,J) IS THE M-TUPLE VALUE OF X ASSOCIATED WITH LOCAL BEST F(I)
             DO J=1,M
             XX(I,J)=A(J)
             ENDDO
             ENDIF
        ENDDO
C      NOW LET EVERY INDIVIDUAL RANDOMLY COSULT NN(<<N) COLLEAGUES AND
C      FIND THE BEST AMONG THEM
      DO I=1,N
C     ------------------------------------------------------------------
      IF(ITOP.GE.3) THEN
C     RANDOM TOPOLOGY ******************************************
C     CHOOSE NN COLLEAGUES RANDOMLY AND FIND THE BEST AMONG THEM
          BEST=1.0D30
           DO II=1,NN
                 CALL RANDOM(RAND)
               NF=INT(RAND*N)+1
                IF(BEST.GT.F(NF)) THEN
                 BEST=F(NF)
                NFBEST=NF
                 ENDIF
            ENDDO
      ENDIF
C----------------------------------------------------------------------
      IF(ITOP.EQ.2) THEN
C     RING + RANDOM TOPOLOGY ******************************************
C     REQUIRES THAT THE SUBROUTINE NEIGHBOR IS TURNED ALIVE
       BEST=1.0D30
          CALL NEIGHBOR(I,N,I1,I3)
          DO II=1,NN
                IF(II.EQ.1) NF=I1
                 IF(II.EQ.2) NF=I
                  IF(II.EQ.3) NF=I3
                      IF(II.GT.3) THEN
                     CALL RANDOM(RAND)
                      NF=INT(RAND*N)+1
                     ENDIF
                  IF(BEST.GT.F(NF)) THEN
                  BEST=F(NF)
                  NFBEST=NF
                 ENDIF
             ENDDO
       ENDIF
C---------------------------------------------------------------------
      IF(ITOP.LE.1) THEN
C     RING TOPOLOGY **************************************************
C     REQUIRES THAT THE SUBROUTINE NEIGHBOR IS TURNED ALIVE
        BEST=1.0D30
           CALL NEIGHBOR(I,N,I1,I3)
              DO II=1,3
              IF (II.NE.I) THEN
             IF(II.EQ.1) NF=I1
              IF(II.EQ.3) NF=I3
                  IF(BEST.GT.F(NF)) THEN
                   BEST=F(NF)
                   NFBEST=NF
                  ENDIF
                  ENDIF
            ENDDO
       ENDIF
C---------------------------------------------------------------------
C     IN THE LIGHT OF HIS OWN AND HIS BEST COLLEAGUES EXPERIENCE, THE
C     INDIVIDUAL I WILL MODIFY HIS MOVE AS PER THE FOLLOWING CRITERION
C     FIRST, ADJUSTMENT BASED ON ONES OWN EXPERIENCE
C     AND OWN BEST EXPERIENCE IN THE PAST (XX(I))
           DO J=1,M
           CALL RANDOM(RAND)
           V1(J)=A1*RAND*(XX(I,J)-X(I,J))
C     THEN BASED ON THE OTHER COLLEAGUES BEST EXPERIENCE WITH WEIGHT W
C     HERE W IS CALLED AN INERTIA WEIGHT 0.01< W < 0.7
C     A2 IS THE CONSTANT NEAR BUT LESS THAN UNITY
           CALL RANDOM(RAND)
           V2(J)=V(I,J)
           IF(F(NFBEST).LT.F(I)) THEN
           V2(J)=A2*W*RAND*(XX(NFBEST,J)-X(I,J))
           ENDIF
C     THEN SOME RANDOMNESS AND A CONSTANT A3 CLOSE TO BUT LESS THAN UNITY
           CALL RANDOM(RAND)
           RND1=RAND
           CALL RANDOM(RAND)
            V3(J)=A3*RAND*W*RND1
C            V3(J)=A3*RAND*W
C     THEN ON PAST VELOCITY WITH INERTIA WEIGHT W
           V4(J)=W*V(I,J)
C     FINALLY A SUM OF THEM
           V(I,J)= V1(J)+V2(J)+V3(J)+V4(J)
           ENDDO
      ENDDO
C     CHANGE X
      DO I=1,N
      DO J=1,M
      RANDS=0.D00
C     ------------------------------------------------------------------
      IF(NSIGMA.EQ.1) THEN
       CALL RANDOM(RAND) ! FOR CHAOTIC PERTURBATION
       IF(DABS(RAND-.5D00).LT.SIGMA) RANDS=RAND-0.5D00
C     SIGMA CONDITIONED RANDS INTRODUCES CHAOTIC ELEMENT IN TO LOCATION
C     IN SOME CASES THIS PERTURBATION HAS WORKED VERY EFFECTIVELY WITH
C     PARAMETER (N=100,NN=15,MX=100,NSTEP=9,ITRN=100000,NSIGMA=1,ITOP=2)
      ENDIF
C     -----------------------------------------------------------------
      X(I,J)=X(I,J)+V(I,J)*(1.D00+RANDS)
      ENDDO
      ENDDO
       DO I=1,N
         IF(F(I).LT.FMIN) THEN
         FMIN=F(I)
         II=I
         DO J=1,M
         BST(J)=XX(II,J)
         ENDDO
         ENDIF
         ENDDO
      IF(LCOUNT.EQ.100) THEN
      LCOUNT=0
      WRITE(*,*)'OPTIMAL SOLUTION UPTO THIS (FUNCTION CALLS=',NFCALL,')'
      WRITE(*,*)'X = ',(BST(J),J=1,M),' MIN F = ',FMIN
C      WRITE(*,*)'NO. OF FUNCTION CALLS = ',NFCALL
      IF(DABS(FMIN-GGBEST).LT.EPSI) THEN
           FMINIM=FMIN
      WRITE(*,*)'COMPUTATION OVER:',FTIT
      WRITE(*,*)'CX=',BST(1),'; CY=',BST(2),'; A=',BST(3),'; B=',BST(4),
     & '; RSQUARE= ',-FMIN
      WRITE(*,*)'PROGAM TERMINATED SUCCESSFULLY. THANK YOU.'
      STOP
           ELSE
           GGBEST=FMIN
           ENDIF
      ENDIF
      LCOUNT=LCOUNT+1
  100 CONTINUE
      FMINIM=FMIN
      WRITE(*,*)'COMPUTATION OVER:',FTIT
      WRITE(*,*)'CX=',BST(1),'; CY=',BST(2),'; A=',BST(3),'; B=',BST(4),
     & '; RSQUARE= ',-FMIN
      WRITE(*,*)'PROGAM TERMINATED SUCCESSFULLY. THANK YOU.'
      END
C     ----------------------------------------------------------------
      SUBROUTINE LSRCH(A,M,VI,NSTEP,FI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
      INTEGER IU,IV
      DIMENSION A(*),B(100),VI(*)
      AMN=1.0D30
      DO J=1,NSTEP
         DO JJ=1,M
         B(JJ)=A(JJ)+(J-(NSTEP/2)-1)*VI(JJ)
         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     -----------------------------------------------------------------
C     THIS SUBROUTINE IS NEEDED IF THE NEIGHBOURHOOD HAS RING TOPOLOGY
C     EITHER PURE OR HYBRIDIZED
       SUBROUTINE NEIGHBOR(I,N,J,K)
       IF(I-1.GE.1 .AND. I.LT.N) THEN
       J=I-1
       K=I+1
       ELSE
       IF(I-1.LT.1) THEN
       J=N-I+1
       K=I+1
       ENDIF
       IF(I.EQ.N) THEN
       J=I-1
       K=1
       ENDIF
       ENDIF
       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 ! INPUT FILE NAME IN WHICH X Y ARE THERE
      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
      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 : REPULSIVE PARTICLE SWARM]'
      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     ------------- 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