Subroutine sce_main

  1 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2 : !
  3 : ! SHUFFLED COMPLEX EVOLUTION
  4 : !   Parameter optimisation algorithm based on a paper by Duan et al.
  5 : !   (1993, J. Opt. Theory and Appl., 76, 501--521).
  6 : !
  7 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  8 : ! VERSION 2.1 ---  01 February 1999
  9 : ! Written by Neil Viney, Centre for Water Research (CWR), The University of WA
 10 : ! Modified by Stan Schymanski, CWR, 05 April 2004 (to run with transpmodel)
 11 : ! Extended by Stan Schymanski, SESE, 02 June 2006 to follow Muttil & Liong
 12 : ! (2004, Journal of Hydraulic Engineering-Asce 130(12):1202-1205) and
 13 : ! Duan et al. (1994, Journal of Hydrology 158)
 14 : !
 15 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 16 : !
 17 : ! This implementation MAXIMISES the objective function, which is
 18 : ! calculated by the model, not by the optimiser. The optimiser transfers
 19 : ! parameter values to the model subroutine and receives the value of the
 20 : ! objective function.
 21 : !
 22 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 23 : !
 24 : ! Copyright (C) 2008 Stan Schymanski
 25 : !
 26 : !   This program is free software: you can redistribute it and/or modify
 27 : !   it under the terms of the GNU General Public License as published by
 28 : !   the Free Software Foundation, either version 3 of the License, or
 29 : !   (at your option) any later version.
 30 : !
 31 : !   This program is distributed in the hope that it will be useful,
 32 : !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 33 : !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 34 : !   GNU General Public License for more details.
 35 : !
 36 : !   You should have received a copy of the GNU General Public License
 37 : !   along with this program. If not, see .
 38 : !
 39 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 40 : 
 41 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 42 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 43 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 44 : 
 45 :       subroutine sce_main ()
 46 :       use vom_sce_mod
 47 :       implicit none
 48 : 
 49 :       INTEGER             :: run_initialseed
 50 :       INTEGER             :: ii, first, numcv
 51 :       REAL*8              :: maxcv
 52 :       CHARACTER(300)      :: writeformat
 53 :       CHARACTER(len=135)  :: msg
 54 :       INTEGER             :: tmp2(2)
 55 :       CHARACTER(len=9)    :: tmp3(1)
 56 : 
 57 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 58 : ! INITIALIZATION
 59 : 
 60 :       call sce_init(run_initialseed)
 61 : 
 62 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 63 : ! BEGIN MODEL LOOP (EXECUTE s PASSES)
 64 : ! CALL model SUBROUTINE        (OPEN run.log FOR CONSOLE OUTPUT);
 65 : ! CALCULATE OBJECTIVE FUNCTION FOR DEFAULT PARAMETER VALUES
 66 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 67 : 
 68 :       if (run_initialseed == 1) then
 69 :         call initialseed
 70 :         return
 71 :       endif
 72 : 
 73 : !
 74 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 75 : ! BEGIN SCE LOOP
 76 : 
 77 :       do while (nrun .lt. 20000 .and. nloop .lt. 500)
 78 : 
 79 :           nloop = nloop + 1
 80 : 
 81 : !         * Saving the best OF of the worst complex in worstbest for
 82 : !         * assessment of gene pool mixing
 83 : 
 84 :           first = 1 + (ncomp2 - 1) * mopt
 85 :           worstbest = ofvec(first)
 86 : !         * [SORT ENTIRE ARRAYS]
 87 : !         * use temporary variable to prevent warning in ifort
 88 :           tmp2(:) = SHAPE(shufflevar(:,:))
 89 :           call sortcomp(shufflevar(:,:), tmp2(:), ofvec(:), SIZE(ofvec(:)))
 90 : 
 91 : !         * WRITE BEST_PARAMETERS FILE FOR PREVIOUS LOOP
 92 : 
 93 :           writeformat = '("Finished ",i4," main loops'
 94 :           writeformat(29:66) = ' --- best objective function =",e12.6)'
 95 :           write(msg,writeformat) nloop, ofvec(1)
 96 :           write(kfile_progress,*) TRIM(msg)
 97 :           write(kfile_progress,*) " "
 98 :           writeformat = '("No improvement in OF for",i5," loops")'
 99 :           write(msg,writeformat) nsincebest
100 :           write(kfile_progress,*) TRIM(msg)
101 :           nsincebest = nsincebest + 1
102 : 
103 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104 : ! ASSESS CONVERGENCE
105 : ! CHANGED BY STAN TO CALCULATE THE FLUCTUATION RANGE RELATIVELY TO
106 : ! THE FEASIBLE RANGE, INSTEAD OF CV:
107 : 
108 :           numcv = ncomp2 * mopt
109 :           sumvar(:) = SUM(shufflevar(optid(:), 1:numcv),2) / numcv  ! mean parameter values
110 :           do ii = 1, nopt
111 : !           * distance from mean in % of feasible range
112 :             cv_(ii) = MAXVAL(ABS((shufflevar(optid(ii), 1:numcv)       &
113 :      &              - sumvar(ii)) / (parmin(optid(ii))                 &
114 :      &              - parmax(optid(ii)))) * 100.d0)
115 :           enddo
116 :           maxcv = MAXVAL(cv_(:))        ! maximum distance
117 :           writeformat = '("Greatest parameter range: ",f5.2,"%'
118 :           writeformat(38:67) = ' for optimised parameter ",a9)'
119 : !         * use temporary variable to prevent warning in ifort
120 :           tmp3(:) = parname(optid(MAXLOC(cv_(:))))
121 :           write(msg,writeformat) maxcv, tmp3
122 :           write(kfile_progress,*) TRIM(msg)
123 :           if (maxcv .ge. i_resolution) then
124 :             if (nsincebest .le. i_patience) then
125 :               call writepars()
126 : 
127 :               call run_cce()
128 : 
129 :               return
130 :             else
131 :               write(kfile_progress,*) " "
132 :               writeformat = '("No improvement in OF for",i5," loops")'
133 :               write(msg,writeformat) nsincebest
134 :               write(kfile_progress,*) TRIM(msg)
135 :               write(kfile_progress,*) "  About to give up..."
136 :             endif
137 :           else
138 :             write(kfile_progress,*) " "
139 :             write(kfile_progress,*) "First Convergence criterion satisfied..."
140 :             write(kfile_progress,*) "  parameter ranges are all less than 0.1 %"
141 :             write(kfile_progress,*) " "
142 :             write(kfile_progress,*) " "
143 :           endif
144 : 
145 : !       * STAN'S MODIFICATION TO ASSESS SENSITIVITY OF OBJECTIVE
146 : !       * FUNCTION TO EACH PARAMETER:
147 : 
148 :           call ck_success()
149 : 
150 :           if (success .eq. 1) then
151 :             !close(kfile_sceout)
152 :             close(kfile_bestpars)
153 :             if (kfile_progress .ne. 6) close(kfile_progress)
154 :             exit
155 :           endif
156 :       enddo
157 : 
158 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159 : ! TERMINATE PROGRAM
160 : 
161 :         if (success .ne. 1) then
162 :           write(kfile_progress,*) " "
163 :           write(kfile_progress,*) "FAILURE TO CONVERGE..."
164 :           write(kfile_progress,*) "  Number of runs has reached 20000."
165 :           write(kfile_progress,*) " "
166 :           write(kfile_progress,*) "  Program terminated."
167 :           write(msg,'(A,A)') CHAR(7), CHAR(7)
168 :           write(kfile_progress,*) TRIM(msg)
169 : !         * [SORT ENTIRE ARRAYS]
170 : !         * use temporary variable to prevent warning in ifort
171 :           tmp2(:) = SHAPE(shufflevar(:,:))
172 :           call sortcomp(shufflevar(:,:), tmp2(:), ofvec(:), SIZE(ofvec(:)))
173 :           call writepars()              ! PROGRAM STOP
174 :           call write_lastbest(shufflevar(:,1), vom_npar, bestobj, 1)
175 :             !close(kfile_sceout)
176 :             close(kfile_bestpars)
177 :           if (kfile_progress .ne. 6) close(kfile_progress)
178 :         endif
179 : 
180 :       return
181 :       end subroutine sce_main