Subroutine sce

  1 :       subroutine sce ()
  2 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  3 : !
  4 : !  SHUFFLED COMPLEX EVOLUTION
  5 : !    Parameter optimisation algorithm based on a paper by Duan et al. 
  6 : !    (1993, J. Opt. Theory and Appl., 76, 501--521). 
  7 : !
  8 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  9 : !  VERSION 2.1        ---  01 February 1999         
 10 : !  Written by Neil Viney, Centre for Water Research (CWR), The University of WA
 11 : !  Modified by Stan Schymanski, CWR, 05 April 2004 (to run with transpmodel)  
 12 : !  Extended by Stan Schymanski, SESE, 02 June 2006 to follow Muttil & Liong (2004,
 13 : !  Journal of Hydraulic Engineering-Asce 130(12):1202-1205) and Duan et al. (1994, 
 14 : !  Journal of Hydrology 158)
 15 : !        
 16 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 17 : !
 18 : !  This implementation MAXIMISES the objective function, which is calculated
 19 : !  by the model, not by the optimiser.        The optimiser transfers parameter values
 20 : !  to the model subroutine and receives the value of the 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 :       use sce_mod
 42 :       implicit none
 43 : 
 44 :       INTEGER :: i_, stat, first, numcv
 45 :       REAL*8  :: maxcv
 46 :       REAL*8, ALLOCATABLE :: sumvar(:)
 47 :       CHARACTER(24)  :: logdate
 48 :       CHARACTER(300) :: writeformat
 49 : 
 50 : !EXTERNAL compar
 51 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 52 : ! Initialize the random number generator
 53 : ! (standard subroutine, based on the date and time)
 54 :       if (command .ne. 4) then
 55 :         CALL RANDOM_SEED()
 56 :       endif
 57 : 
 58 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 59 : ! WRITE SCREEN HEADER
 60 : 
 61 :       write(*,'(/"SHUFFLED COMPLEX EVOLUTION OPTIMISER")')
 62 :       call fdate(logdate)
 63 :       write(*,'(/"  Run time:   ",a)') logdate
 64 : 
 65 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 66 : ! INITIALIZATION
 67 : 
 68 :       call sce_init()
 69 : 
 70 :       allocate(sumvar(nopt))
 71 : 
 72 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 73 : ! BEGIN MODEL LOOP (EXECUTE s PASSES)
 74 : ! CALL model SUBROUTINE        (OPEN run.log FOR CONSOLE OUTPUT);
 75 : ! CALCULATE OBJECTIVE FUNCTION FOR DEFAULT PARAMETER VALUES
 76 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 77 : ! INSERTED BY STAN TO ALLOW CONTINUATION OF OPTIMSATION FROM PREVIOUSLY SAVED STEP
 78 : 
 79 :       open(kfile_lastloop, file=sfile_lastloop(1:len_trim(sfile_lastloop)), status='old', &
 80 :       & iostat=stat)
 81 :       if (stat .eq. 0) then
 82 :         read(kfile_lastloop,*) ncomp2
 83 :         read(kfile_lastloop,*) nloop
 84 :         read(kfile_lastloop,*) nrun
 85 :         read(kfile_lastloop,*) nsincebest
 86 :         read(kfile_lastloop,loopformat) ofvec(:)
 87 :         do i_ = 1, npar
 88 :           read(kfile_lastloop, loopformat) shufflevar(i_,:)
 89 :         enddo
 90 :         close(kfile_lastloop)
 91 :         bestobj = ofvec(1)
 92 : 
 93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 94 : ! OPEN FILES FOR STORING OBJECTIVE FUNCTION AND PARAMETER VALUES
 95 : 
 96 :         open(kfile_sceout, file=sfile_sceout(1:len_trim(sfile_sceout)), status='old', &
 97 :         & position='append')
 98 :         open(kfile_bestpars, file=sfile_bestpars(1:len_trim(sfile_bestpars)), status='old', &
 99 :         & position='append')
100 :         open(kfile_progress, file=sfile_progress(1:len_trim(sfile_progress)), status='old', &
101 :         & position='append')
102 :         write(kfile_progress,'(/"  NEW Run time:   ",a)') logdate
103 :       else
104 : 
105 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 : ! OPEN AND EMPTY FILE FOR STORING OBJECTIVE FUNCTION AND PARAMETER VALUES
107 : 
108 :         open(kfile_sceout, file=sfile_sceout(1:len_trim(sfile_sceout)), status='replace')
109 :         open(kfile_bestpars, file=sfile_bestpars(1:len_trim(sfile_bestpars)), status='replace')
110 :         open(kfile_progress, file=sfile_progress(1:len_trim(sfile_progress)), status='replace')
111 : 
112 : !       * write file header
113 :         write(kfile_progress,'(/"SHUFFLED COMPLEX EVOLUTION OPTIMISER")')
114 :         write(kfile_progress,'(/"  Run time:   ",a)') logdate
115 : 
116 :         call initialseed
117 : 
118 :         close(kfile_sceout)
119 :         close(kfile_bestpars)
120 :         close(kfile_progress)
121 :         return
122 :       endif
123 : 
124 : ! END OF INSERTION
125 : !
126 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127 : ! BEGIN SCE LOOP
128 : 
129 :       do while (nrun .lt. 20000 .and. nloop .lt. 500)
130 :         nloop = nloop + 1
131 : 
132 : !       * Saving the best OF of the worst complex in worstbest for
133 : !       * assessment of gene pool mixing
134 : 
135 :         first = 1 + (ncomp2 - 1) * mopt
136 :         worstbest = ofvec(first)
137 :         call sortcomp(shufflevar(:,:), shape(shufflevar(:,:)), ofvec(:), size(ofvec(:)))  ! [SORT ENTIRE ARRAYS]
138 : 
139 : !       * WRITE BEST_PARAMETERS FILE FOR PREVIOUS LOOP
140 : 
141 :         writeformat = '("Finished ",i4," main loops'
142 :         writeformat(29:66) = ' --- best objective function =",e12.6)'
143 :         write(*,writeformat) nloop, ofvec(1)
144 :         write(kfile_progress,writeformat) nloop, ofvec(1)
145 :         writeformat = '(/"No improvement in OF for",i5," loops")'
146 :         write(*,writeformat) nsincebest
147 :         write(kfile_progress,writeformat) nsincebest
148 :         nsincebest = nsincebest + 1
149 :     
150 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151 : ! ASSESS CONVERGENCE
152 : ! CHANGED BY STAN TO CALCULATE THE FLUCTUATION RANGE RELATIVELY TO
153 : ! THE FEASIBLE RANGE, INSTEAD OF CV:
154 : 
155 :         numcv = ncomp2 * mopt
156 :         sumvar(:) = sum(shufflevar(optid(:), 1:numcv),2) / numcv  ! mean parameter values
157 :         do i_ = 1, nopt
158 :           cv_(i_) = maxval(abs((shufflevar(optid(i_), 1:numcv) - sumvar(i_))  &
159 :      &          / (parmin(optid(i_)) - parmax(optid(i_)))) * 100.d0)  ! distance from mean in % of feasible range
160 :         enddo
161 :         maxcv = maxval(cv_(:))                ! maximum distance
162 :         writeformat = '("Greatest parameter range: ",f5.2,"%'
163 :         writeformat(38:67) = ' for optimised parameter ",a9)'
164 :         write(*,writeformat) maxcv, parname(optid(maxloc(cv_(:))))
165 :         write(kfile_progress,writeformat) maxcv, parname(optid(maxloc(cv_(:))))
166 :         if (maxcv .ge. resolution) then
167 :           if (nsincebest .le. patience) then
168 :             call writepars
169 :             call run_cce()
170 :             return
171 :           else
172 :             writeformat = '(/"No improvement in OF for",i5," loops"/" '
173 :             writeformat(44:65) = ' About to give up...")'
174 :             write(*,writeformat) nsincebest
175 :             write(kfile_progress,writeformat) nsincebest
176 :             call ck_success()
177 :             if (success .eq. 1) return
178 :           endif
179 :         else
180 :           writeformat = '(/"First Convergence criterion satisfied..."/"'
181 :           writeformat(47:92) = '  parameter ranges are all less than 0.1 %"//)'
182 :           write(*,writeformat)
183 :           write(kfile_progress,writeformat)
184 : 
185 : !         * STAN'S MODIFICATION TO ASSESS SENSITIVITY OF OBJECTIVE
186 : !         * FUNCTION TO EACH PARAMETER:
187 : 
188 :           call ck_success()
189 :           if (success .eq. 1) return
190 :         endif
191 :       enddo
192 : 
193 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
194 : ! END OF BLOCK THAT CAN BE PARALLELISED!
195 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
196 : ! TERMINATE PROGRAM
197 : 
198 :       writeformat = '(/"FAILURE TO CONVERGE..."/"  Number of runs has '
199 :       writeformat(50:96) = 'reached 20000." //"  Program terminated."/,a,a)'
200 :       write(*,writeformat) char(7), char(7)
201 :       write(kfile_progress,writeformat) char(7), char(7)
202 :       call sortcomp(shufflevar(:,:), shape(shufflevar(:,:)), ofvec(:), size(ofvec(:)))  ! [SORT ENTIRE ARRAYS]
203 :       call writepars                              ! PROGRAM STOP
204 :       open(kfile_finalbest, file=sfile_finalbest(1:len_trim(sfile_finalbest)))
205 :       write(kfile_finalbest,outformat) shufflevar(:,1), bestobj
206 :       close(kfile_finalbest)
207 :       close(kfile_sceout)
208 :       close(kfile_bestpars)
209 :       close(kfile_progress)
210 : 
211 :       deallocate(sumvar)
212 : 
213 :       return
214 :       end subroutine sce