Subroutine simplex

942 : 
943 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
944 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
945 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
946 : 
947 :       subroutine simplex (invar, objfun)
948 :       use vom_sce_mod
949 : 
950 :       implicit none
951 : 
952 : !     * Declarations
953 :       REAL*8, DIMENSION(vom_npar,qopt), INTENT(inout) :: invar
954 :       REAL*8, DIMENSION(qopt),          INTENT(inout) :: objfun
955 : 
956 : !     * Definitions
957 :       INTEGER       :: l_, kk
958 :       INTEGER       :: ii, jj
959 :       INTEGER       :: max_j_
960 :       REAL*8        :: newobjfun
961 :       REAL*8        :: minj, rangej
962 :       REAL*8, DIMENSION(vom_npar) :: centroid  ! Centroid of parameter sets for simplex procedure
963 :       REAL*8, DIMENSION(vom_npar) :: newpoint  ! New parameter set resulting from simplex procedure
964 :       REAL*8, ALLOCATABLE :: ranarr_simplex(:)  ! Array of random numbers
965 : 
966 :       allocate(ranarr_simplex(nopt))
967 : 
968 :       do l_ = 1, i_nsimp
969 : 
970 : !         * REFLECTION STEP
971 : 
972 :           evolution = 'reflection'
973 :           newpoint(:) = invar(:,qopt)
974 :           centroid(optid(:)) = 1.d0 / (qopt - 1) * SUM(invar(optid(:), 1:qopt - 1), 2)
975 :           newpoint(optid(:)) = 2.d0 * centroid(optid(:)) - invar(optid(:),qopt)
976 :           if (MINVAL(newpoint(:) - parmin(:)) .lt. 0.d0 .or.           &
977 :      &        MAXVAL(newpoint(:) - parmax(:)) .gt. 0.d0) then
978 : 
979 : !           * MUTATION STEP
980 : !           * NB: MUTATION IS BASED ON THE SMALLEST HYPERCUBE OF THE SUBCOMPLEX,
981 : !           * NOT THE SMALLEST HYPERCUBE OF THE COMPLEX AS SUGGESTED BY DUAN ET AL.
982 : 
983 :             do ii = 1, nopt
984 :               call random_number(ranscal)  ! SCALAR RANDOM NUMBER
985 :               jj = optid(ii)
986 :               minj = parmin(jj)
987 :               rangej = parmax(jj) - minj
988 :               newpoint(jj) = minj + rangej * ranscal  ! UNIFORMLY DISTRIBUTED SAMPLE
989 :             enddo
990 :             evolution = 'mutation'
991 :           endif
992 : 
993 :           kk = 1
994 :         do while (kk .le. 3)
995 : 
996 : !         * CALCULATE OBJECTIVE FUNCTION
997 :           call runmodel(newpoint(:), newobjfun)
998 : 
999 :             jj = 1
1000 :             if (kk .eq. 1) then
1001 :               if (newobjfun .gt. objfun(qopt)) then
1002 :                 max_j_ = qopt
1003 :                 kk = 4
1004 :               else
1005 : !               * CONTRACTION STEP
1006 :                 newpoint(optid(:)) = (centroid(optid(:)) + invar(optid(:),qopt)) * 0.5d0
1007 :                 evolution = 'contraction'
1008 :               endif
1009 :             endif
1010 : 
1011 :             if (kk .eq. 2) then
1012 :               if (newobjfun .gt. objfun(qopt)) then
1013 :                 max_j_ = qopt
1014 :                 kk = 4
1015 :               else
1016 : !               * ANOTHER MUTATION STEP
1017 :                 do ii = 1, nopt
1018 :                   call random_number(ranarr_simplex(:))  ! RANDOM ARRAY OF SIZE 1:nopt
1019 :                   jj = optid(ii)
1020 :                   minj = MINVAL(invar(jj,:))
1021 :                   rangej = MAXVAL(invar(jj,:)) - minj
1022 :                   newpoint(jj) = minj + rangej * ranarr_simplex(ii)  ! UNIFORMLY DISTRIBUTED SAMPLE
1023 :                 enddo
1024 :                 evolution = 'mutation'
1025 :               endif
1026 :             endif
1027 : 
1028 :             if (kk .eq. 3) then
1029 :               if (newobjfun .gt. objfun(qopt-1)) then
1030 :                 max_j_ = qopt - 1
1031 :                 kk = 4
1032 :               elseif (newobjfun .gt. 0.d0) then  ! REPLACE ANYWAY
1033 :                 objfun(qopt) = newobjfun
1034 :                 invar(:,qopt) = newpoint(:)
1035 :               endif
1036 :             endif
1037 : 
1038 :             if (kk .eq. 4) then
1039 : !            * SORT objfun HERE IN CASE alpha > 1
1040 :               do while (newobjfun .le. objfun(jj) .and. jj .le. max_j_)
1041 :                 jj = jj + 1
1042 :               enddo
1043 :               if (jj .lt. qopt) then
1044 :                 objfun(jj+1:qopt) = objfun(jj:qopt-1)
1045 :                 invar(:,jj+1:qopt) = invar(:,jj:qopt-1)
1046 :               endif
1047 :               if (jj .le. qopt) then
1048 :                 objfun(jj) = newobjfun
1049 :                 invar(:,jj) = newpoint(:)
1050 :               endif
1051 :             endif
1052 : 
1053 :             kk = kk + 1
1054 : 
1055 :           enddo
1056 : 
1057 : !         * UPDATE 'currentbest' IF APPROPRIATE
1058 : 
1059 :           if (newobjfun .gt. bestobj) then
1060 :             bestobj = newobjfun
1061 :             bestincomp = newobjfun
1062 :             call write_lastbest(invar(:,1), vom_npar, bestobj, 0)
1063 :             nsincebest = 0
1064 :           elseif (newobjfun .gt. bestincomp) then
1065 :             bestincomp = newobjfun
1066 :           endif
1067 : 
1068 :       enddo
1069 : 
1070 :       deallocate(ranarr_simplex)
1071 : 
1072 :       return
1073 :       end subroutine simplex