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