C ----------------------------------------------------------------- C PROGRAM: LOGSPIRAL - FOR FITTING LOG SPIRAL TO DATA WITH DISPLACED C ORIGIN BY DIFFERENTIAL EVOLUTION METHOD OF GLOBAL OPTIMIZATION. C ------------------------------------------------------------------ C "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 SUBROUTINE DE IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! TYPE DECLARATION PARAMETER(NMAX=500,MMAX=50) ! MAXIMUM DIMENSION PARAMETERS PARAMETER (RX1=0.0, RX2=0.8) ! 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 /XYDATA/X0,Y0,ALIM,NPOINTS DIMENSION X0(200),Y0(200),ALIM(2,4) ! IF MORE DATA INCREASE THESE 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 ---------------------------------------------------------------- M=4 ! NO. OF PARAMETERS C ----------------------------------------------------------------- C SPECIFY OTHER PARAMETERS --------------------------------------- WRITE(*,*)'POPULATION SIZE [N] AND NO. OF ITERATIONS [ITER] ?' WRITE(*,*)'SUGGESTED : N => 40 ; 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)*200 ! 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 IF(DABS(FBEST-GBEST).LT.EPS) THEN WRITE(*,*) FTIT WRITE(*,*)'NO. OF VARIABLES =', M WRITE(*,*)'======================================================' WRITE(*,*)'CX=',A(1),'; CY=',A(2),'; A=',A(3),'; B=',A(4), & '; RSQUARE= ',-FBEST WRITE(*,*)'COMPUTATION OVER. THANK YOU' STOP ELSE GBEST=FBEST ENDIF IPCOUNT=0 ENDIF C ---------------------------------------------------------------- 100 ENDDO ! ITERATION ENDS : GO FOR NEXT ITERATION, IF APPLICABLE C ---------------------------------------------------------------- WRITE(*,*)'======================================================' WRITE(*,*)'CX=',A(1),'; CY=',A(2),'; A=',A(3),'; B=',A(4), & '; RSQUARE= ',-FBEST 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(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 : DIFFERENTIAL EVOLUTION]' 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