-
Notifications
You must be signed in to change notification settings - Fork 6
/
mo_sce.f90
1951 lines (1820 loc) · 87.4 KB
/
mo_sce.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! ToDo: be able to switch of writing of restart files (e.g. restartfile1='')
!> \file mo_sce.f90
!> \brief Shuffled Complex Evolution optimization algorithm.
!> \details Optimization algorithm using Shuffled Complex Evolution strategy.
!> Original version 2.1 of Qingyun Duan (1992) rewritten in Fortran 90.
!> \authors Juliane Mai
!> \date Feb 2013
MODULE mo_sce
! This module is the Shuffled Complex Evolution optimization algorithm.
! Written Juliane Mai, Feb 2013
! License
! -------
! This file is part of the JAMS Fortran package, distributed under the MIT License.
!
! Copyright (c) 2011-2014 Juliane Mai, Matthias Cuntz - mc (at) macu (dot) de
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
IMPLICIT NONE
PUBLIC :: sce ! sce optimization
! ------------------------------------------------------------------
! NAME
! sce
! PURPOSE
!> \brief Shuffled Complex Evolution (SCE) algorithm for global optimization.
!
!> \details Shuffled Complex Evolution method for global optimization\n
!> -- version 2.1\n
!> \n
!> by Qingyun Duan\n
!> Department of Hydrology & Water Resources\n
!> University of Arizona, Tucson, AZ 85721\n
!> (602) 621-9360, email: duan@hwr.arizona.edu\n
!> \n
!> Written by Qingyun Duan, Oct 1990.\n
!> Revised by Qingyun Duan, Aug 1991.\n
!> Revised by Qingyun Duan, Apr 1992.\n
!> \n
!> Re-written by Juliane Mai, Feb 2013.\n
!> \n
!> Statement by Qingyun Duan:\n
!> ----------------------------------\n
!> \n
!> This general purpose global optimization program is developed at
!> the Department of Hydrology & Water Resources of the University
!> of Arizona. Further information regarding the SCE-UA method can
!> be obtained from Dr. Q. Duan, Dr. S. Sorooshian or Dr. V.K. Gupta
!> at the address and phone number listed above. We request all
!> users of this program make proper reference to the paper entitled
!> 'Effective and Efficient Global Optimization for Conceptual
!> Rainfall-runoff Models' by Duan, Q., S. Sorooshian, and V.K. Gupta,
!> Water Resources Research, Vol 28(4), pp.1015-1031, 1992.\n
!> \n
!> The function to be minimized is the first argument of DDS and must be defined as \n
!> \code
!> function functn(p)
!> use mo_kind, only: dp
!> implicit none
!> real(dp), dimension(:), intent(in) :: p
!> real(dp) :: functn
!> end function functn
!> \endcode
!
! INTENT(IN)
!> \param[in] "real(dp) :: functn(p)" Function on which to search the optimum
!> \param[in] "real(dp) :: pini(:)" inital value of decision variables
!> \param[in] "real(dp) :: prange(size(pini),2)" Min/max range of decision variables
!
! INTENT(INOUT)
! None
! INTENT(OUT)
! None
!
! INTENT(IN), OPTIONAL
!> \param[in] "integer(i8), optional :: mymaxn" max no. of trials allowed before optimization is terminated\n
!> DEFAULT: 1000_i8
!> \param[in] "logical, optional :: mymaxit" maximization (.true.) or minimization (.false.) of function\n
!> DEFAULT: false
!> \param[in] "integer(i4), optional :: mykstop" number of shuffling loops in which the criterion value must
!> change by given percentage before optimiz. is terminated\n
!> DEFAULT: 10_i4
!> \param[in] "real(dp), optional :: mypcento" percentage by which the criterion value must change in
!> given number of shuffling loops\n
!> DEFAULT: 0.0001_dp
!> \param[in] "real(dp), optional :: mypeps" optimization is terminated if volume of complex has
!> converged to given percentage of feasible space\n
!> DEFAULT: 0.001_dp
!> \param[in] "integer(i8), optional :: myseed" initial random seed\n
!> DEFAULT: get_timeseed
!> \param[in] "integer(i4), optional :: myngs" number of complexes in the initial population\n
!> DEFAULT: 2_i4
!> \param[in] "integer(i4), optional :: mynpg" number of points in each complex\n
!> DEFAULT: 2*n+1
!> \param[in] "integer(i4), optional :: mynps" number of points in a sub-complex\n
!> DEFAULT: n+1
!> \param[in] "integer(i4), optional :: mynspl" number of evolution steps allowed for each complex before
!> complex shuffling\n
!> DEFAULT: 2*n+1
!> \param[in] "integer(i4), optional :: mymings" minimum number of complexes required, if the number of
!> complexes is allowed to reduce as the
!> optimization proceeds\n
!> DEFAULT: ngs = number of complexes in initial population
!> \param[in] "integer(i4), optional :: myiniflg" flag on whether to include the initial point in population\n
!> 0, not included\n
!> 1, included (DEFAULT)
!> \param[in] "integer(i4), optional :: myprint" flag for controlling print-out after each shuffling loop\n
!> 0, print information on the best point of the population\n
!> 1, print information on every point of the population\n
!> 2, no printing (DEFAULT)
!> 3, same as 0 but print progress '.' on every function call
!> 4, same as 1 but print progress '.' on every function call
!> \param[in] "logical, optional :: mymask(size(pini))"
!> parameter included in optimization (true) or discarded (false)
!> DEFAULT: .true.
!> \param[in] "real(dp), optional :: myalpha" parameter for reflection of points in complex\n
!> DEFAULT: 0.8_dp
!> \param[in] "real(dp), optional :: mybeta" parameter for contraction of points in complex\n
!> DEFAULT: 0.45_dp
!> \param[in] "character(len=*), optional :: tmp_file" if given: write results after each evolution loop
!> to temporal output file of that name\n
!> # of headlines: 7\n
!> format: '# nloop icall ngs1 bestf worstf ... \n
!> ... gnrng (bestx(j),j=1,nn)'
!> \param[in] "character(len=*), optional :: popul_file" if given: write whole population to file of that name\n
!> # of headlines: 1 \n
!> format: #_evolution_loop, xf(i), (x(i,j),j=1,nn)\n
!> total number of lines written <= neval <= mymaxn\n
!> \param[in] "logical, optional :: popul_file_append" if true, append to existing population file (default: false)\n
!> \param[in] "logical, optional :: parallel" sce runs in parallel (true) or not (false)
!> parallel sce should only be used if model/ objective
!> is not parallel
!> DEFAULT: .false.
!> \param[in] "logical, optional :: restart" if .true.: restart former sce run from restart_file
!> \param[in] "character(len=*), optional :: restart_file" file name for read/write of restart file
!> (default: mo_sce.restart)
!
! INTENT(INOUT), OPTIONAL
! None
!
! INTENT(OUT), OPTIONAL
!> \param[out] "real(dp), optional :: bestf" the best value of the function.
!> \param[out] "integer(i8), optional :: neval" number of function evaluations needed.
!> \param[out] "real(dp), optional, allocatable :: history(:)" the history of best function values,
!> history(neval)=bestf
!
! RETURN
!> \return real(dp) :: bestx(size(pini)) — The parameters of the point which is estimated
!> to minimize/maximize the function.
!
! RESTRICTIONS
!> \note Maximal number of parameters is 1000.\n
!> SCE is OpenMP enabled on the loop over the complexes.\n
!> OMP_NUM_THREADS > 1 does not give reproducible results even when seeded!
!
! EXAMPLE
! use mo_opt_functions, only: griewank
!
! prange(:,1) = (/ -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0, -600.0 /)
! prange(:,2) = (/ 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0, 600.0 /)
! pini = (/ -.226265E+01, -.130187E+01, -.151219E+01, 0.133983E+00, 0.988159E+00, &
! -.495074E+01, -.126574E+02, 0.572684E+00, 0.303864E+01, 0.343031E+01 /)
!
! opt = sce(griewank, pini, prange)
!
! -> see also example in test directory
! LITERATURE
! Duan, Q., S. Sorooshian, and V.K. Gupta -
! Effective and Efficient Global Optimization for Conceptual Rainfall-runoff Models
! Water Resources Research, Vol 28(4), pp.1015-1031, 1992.
! HISTORY
!> \author Juliane Mai
!> \date Feb 2013
! Modified Juliane Mai, Matthias Cuntz, Jul 2013 - OpenMP
! - NaN and Inf in objective function
! Juliane Mai, Oct 2013 - added peps as optional argument
! - allow for masked parameters
! - write population to file --> popul_file
! - write intermediate results to file --> tmp_file
! - flag parallel introduced
! Matthias Cuntz, Nov 2013 - progress dots
! - use iso_fortran_env
! - treat functn=NaN as worse function value in cce
! Matthias Cuntz, May 2014 - sort -> orderpack
! Matthias Cuntz, May 2014 - popul_file_append
! Matthias Cuntz, May 2014 - sort with NaNs
! Matthias Cuntz, Juliane Mai Feb 2015 - restart
! Matthias Cuntz Mar 2015 - use is_finite and special_value from mo_utils
! Juliane Mai Apr 2015 - handling of array-like variables in restart-namelists
! Matthias Cuntz Feb 2018 - made counter local in contained subroutines/functions
! Matthias Cuntz Mar 2018 - mask parameters with degenerated ranges, e.g. upper<lower
! Matthias Cuntz Mar 2018 - use only masked in normalized geometric range of parameters
! ------------------------------------------------------------------
PRIVATE
! chkcst ! check if a trial point satisfies all constraints.
! comp ! reduces the number of complexes and points in a population
! parstt ! checks for parameter convergence
! getpnt ! generated a new point within its range
! cce ! generates new point(s) from sub-complex
! ------------------------------------------------------------------
CONTAINS
function sce(functn,pini,prange, & ! IN
mymaxn,mymaxit,mykstop,mypcento,mypeps,myseed, & ! Optional IN
myngs,mynpg,mynps,mynspl,mymings,myiniflg,myprint, & ! Optional IN
mymask,myalpha, mybeta, & ! Optional IN
tmp_file, popul_file, & ! Optional IN
popul_file_append, & ! Optional IN
parallel, restart, restart_file, & ! OPTIONAL IN
bestf,neval,history & ! Optional OUT
) result(bestx)
use mo_kind, only: i4, i8, dp
use mo_orderpack, only: sort
use mo_string_utils, only: num2str, compress
use mo_xor4096, only: get_timeseed, n_save_state, xor4096, xor4096g
!$ use omp_lib, only: OMP_GET_THREAD_NUM, OMP_GET_NUM_THREADS
use iso_fortran_env, only: output_unit, error_unit
use mo_utils, only: is_finite, special_value
implicit none
!
interface
function functn(paraset)
! calculates the cost function at a certain parameter set paraset
use mo_kind, only: dp
real(dp), dimension(:), intent(in) :: paraset
real(dp) :: functn
end function functn
end interface
real(dp), dimension(:), intent(in) :: pini ! initial parameter set
real(dp), dimension(:,:), intent(in) :: prange ! lower and upper bound on parameters
!
! algorithmic control & convergence parameter
integer(i8), optional, intent(in) :: mymaxn ! max no. of trials allowed before optimization is terminated
! ! DEFAULT: 1000_i8
logical, optional, intent(in) :: mymaxit ! maximization (.true.) or minimization (.false.) of function
! ! DEFAULT: false
integer(i4), optional, intent(in) :: mykstop ! number of shuffling loops in which the criterion value must
! ! change by given percentage before optimiz. is terminated
! ! DEFAULT: 10_i4
real(dp), optional, intent(in) :: mypcento ! percentage by which the criterion value must change in
! ! given number of shuffling loops
! ! DEFAULT: 0.0001_dp
real(dp), optional, intent(in) :: mypeps ! optimization is terminated if volume of complex has
! ! converged to given percentage of feasible space
! ! DEFAULT: 0.001_dp
integer(i8), optional, intent(in) :: myseed ! initial random seed
! ! DEFAULT: get_timeseed
integer(i4), optional, intent(in) :: myngs ! number of complexes in the initial population
! ! DEFAULT: 2_i4
integer(i4), optional, intent(in) :: mynpg ! number of points in each complex
! ! DEFAULT: 2*n+1
integer(i4), optional, intent(in) :: mynps ! number of points in a sub-complex
! ! DEFAULT: n+1
integer(i4), optional, intent(in) :: mynspl ! number of evolution steps allowed for each complex before
! ! complex shuffling
! ! DEFAULT: 2*n+1
integer(i4), optional, intent(in) :: mymings ! minimum number of complexes required, if the number of
! ! complexes is allowed to reduce as the
! ! optimization proceeds
! ! DEFAULT: ngs = number of complexes in initial population
integer(i4), optional, intent(in) :: myiniflg ! flag on whether to include the initial point in population
! ! 0, not included
! ! 1, included (DEFAULT)
integer(i4), optional, intent(in) :: myprint ! flag for controlling print-out after each shuffling loop
! ! 0, print information on the best point of the population
! ! 1, print information on every point of the population
! ! 2, no printing (DEFAULT)
logical, optional, intent(in), &
dimension(size(pini,1)) :: mymask ! parameter included in optimization (true) or discarded (false)
! ! DEFAULT: .true.
real(dp), optional, intent(in) :: myalpha ! parameter for reflection of points in complex
! ! DEFAULT: 0.8_dp
real(dp), optional, intent(in) :: mybeta ! parameter for contraction of points in complex
! ! DEFAULT: 0.45_dp
character(len=*), optional, intent(in) :: tmp_file ! file for temporal output: write results after evolution loop
! ! # of headlines: 7
! ! format: '# nloop icall ngs1 bestf worstf ...
! ! ... gnrng (bestx(j),j=1,nn)'
character(len=*), optional, intent(in) :: popul_file ! file for temporal output: writes whole population
! ! # of headlines: 1
! ! format: #_evolution_loop, xf(i), (x(i,j),j=1,nn)
! ! total number of lines written <= neval <= mymaxn
logical, optional, intent(in) :: popul_file_append ! if true, append to popul_file
logical, optional, intent(in) :: parallel ! sce runs in parallel (true) or not (false)
! ! parallel sce should only be used if model/ objective
! ! is not parallel
! ! DEAFULT: .false.
logical, optional, intent(in) :: restart ! if .true., start from restart file
character(len=*), optional, intent(in) :: restart_file ! restart file name (default: mo_sce.restart)
real(dp), optional, intent(out) :: bestf ! function value of bestx(.)
integer(i8), optional, intent(out) :: neval ! number of function evaluations
real(dp), optional, &
dimension(:), allocatable, intent(out) :: history ! history of best function values after each iteration
real(dp), dimension(size(pini,1)) :: bestx ! best point at current shuffling loop (is RETURN)
!
! optionals transfer
real(dp), dimension(size(pini,1)) :: bl ! lower bound on parameters
real(dp), dimension(1000) :: dummy_bl ! dummy lower bound (only for namelist)
real(dp), dimension(size(pini,1)) :: bu ! upper bound on parameters
real(dp), dimension(1000) :: dummy_bu ! dummy upper bound (only for namelist)
integer(i8) :: maxn ! max no. of trials allowed before optimization is terminated
logical :: maxit ! minimization (false) or maximization (true)
integer(i4) :: kstop ! number of shuffling loops in which the criterion value
! ! must change
real(dp) :: pcento ! percentage by which the criterion value must change
real(dp) :: peps ! optimization is terminated if volume of complex has
! ! converged to given percentage of feasible space
integer(i4) :: ngs ! number of complexes in the initial population
integer(i4) :: npg ! number of points in each complex
integer(i4) :: nps ! number of points in a sub-complex
integer(i4) :: nspl ! number of evolution steps allowed for each complex before
! ! complex shuffling
integer(i4) :: mings ! minimum number of complexes required
integer(i4) :: iniflg ! flag on whether to include the initial point in population
integer(i4) :: iprint ! flag for controlling print-out after each shuffling loop
logical :: idot ! controls progress report .
logical, dimension(size(pini,1)) :: maskpara ! mask(i) = .true. --> parameter i will be optimized
! ! mask(i) = .false. --> parameter i will not be optimized
logical, dimension(1000) :: dummy_maskpara ! dummy maskpara (only for namelist)
real(dp) :: alpha ! parameter for reflection of points in complex
real(dp) :: beta ! parameter for contraction of points in complex
logical :: itmp_file ! if temporal results wanted
#ifdef __pgiFortran__
character(len=299) :: istmp_file ! local copy of file for temporal results output
#else
character(len=1024) :: istmp_file ! local copy of file for temporal results output
#endif
logical :: ipopul_file ! if temporal population wanted
#ifdef __pgiFortran__
character(len=299) :: ispopul_file ! local copy of file for temporal population output
#else
character(len=1024) :: ispopul_file ! local copy of file for temporal population output
#endif
real(dp), dimension(:), allocatable :: history_tmp ! history of best function values after each iteration
real(dp), dimension(:), allocatable :: ihistory_tmp ! local history for OpenMP
real(dp) :: bestf_tmp ! function value of bestx(.)
logical :: parall ! if sce is used parallel or not
logical :: irestart ! if restart wanted
character(len=1024) :: isrestart_file ! local copy of restart file name
!
! local variables
integer(i4) :: nopt ! number of parameters to be optimized
integer(i4) :: nn ! total number of parameters
integer(i4) :: npt ! total number of points in initial population (npt=ngs*npg)
real(dp) :: fpini ! function value at initial point
real(dp), dimension(:,:), allocatable :: x ! coordinates of points in the population
real(dp), dimension(:), allocatable :: xf ! function values of x
real(dp), dimension(size(pini,1)) :: xx ! coordinates of a single point in x
real(dp), dimension(1000) :: dummy_xx ! dummy xx (only for namelist)
real(dp), dimension(:,:), allocatable :: cx ! coordinates of points in a complex
real(dp), dimension(:), allocatable :: cf ! function values of cx
real(dp), dimension(:,:), allocatable :: s ! coordinates of points in the current sub-complex !simplex
real(dp), dimension(:), allocatable :: sf ! function values of s
real(dp), dimension(:), allocatable :: worstx ! worst point at current shuffling loop
real(dp) :: worstf ! function value of worstx(.)
real(dp), dimension(:), allocatable :: xnstd ! standard deviation of parameters in the population
real(dp) :: gnrng ! normalized geometric mean of parameter ranges
integer(i4), dimension(:), allocatable :: lcs ! indices locating position of s(.,.) in x(.,.)
real(dp), dimension(:), allocatable :: bound ! length of range of ith variable being optimized
real(dp), dimension(:), allocatable :: unit ! standard deviation of initial parameter distributions
integer(i4) :: ngs1 ! number of complexes in current population
integer(i4) :: ngs2 ! number of complexes in last population
real(dp), dimension(:), allocatable :: criter ! vector containing the best criterion values of the last
! ! (kstop+1) shuffling loops
integer(i4) :: ipcnvg ! flag indicating whether parameter convergence is reached
! ! (i.e., check if gnrng is less than 0.001)
! ! 0, parameter convergence not satisfied
! ! 1, parameter convergence satisfied
integer(i4) :: nloop
integer(i4) :: loop
integer(i4) :: ii, jj, kk
integer(i4) :: lpos
logical :: lpos_ok ! for selction of points based on triangular
! ! probability distribution
integer(i4) :: npt1
integer(i8) :: icall ! counter for function evaluations
integer(i8) :: icall_merk, iicall ! local icall for OpenMP
integer(i4) :: igs
integer(i4) :: k1, k2
real(dp) :: denomi ! for checking improvement of last steps
real(dp) :: timeou ! for checking improvement of last steps
character(4), dimension(:), allocatable :: xname ! parameter names: "p1", "p2", "p3", ...
character(512) :: format_str1 ! format string
character(512) :: format_str2 ! format string
! for random numbers
real(dp) :: rand ! random number
integer(i8), dimension(:), allocatable :: iseed ! initial random seed
integer(i8), dimension(:,:), allocatable :: save_state_unif ! save state of uniform stream
integer(i8), dimension(:,:), allocatable :: save_state_gauss ! save state of gaussian stream
real(dp), dimension(:), allocatable :: rand_tmp ! random number
integer(i4) :: ithread, n_threads ! OMP or not
integer(i8), dimension(n_save_state) :: itmp
real(dp), dimension(:,:), allocatable :: xtmp ! tmp array for complex reduction
real(dp), dimension(:), allocatable :: ftmp ! %
real(dp) :: large ! for treating NaNs
logical :: ipopul_file_append
integer(i4) :: nonan ! # of non-NaN in history_tmp
real(dp), dimension(:), allocatable :: htmp ! tmp storage for history_tmp
namelist /restartnml1/ &
bestx, dummy_bl, dummy_bu, maxn, maxit, kstop, pcento, peps, ngs, npg, nps, nspl, &
mings, iniflg, iprint, idot, dummy_maskpara, alpha, beta, bestf_tmp, &
parall, nopt, nn, npt, fpini, dummy_xx, worstf, gnrng, &
ngs1, ngs2, nloop, npt1, icall, &
n_threads, ipopul_file_append, itmp_file, istmp_file, &
ipopul_file, ispopul_file
namelist /restartnml2/ &
history_tmp, x, xf, worstx, xnstd, &
bound, unit, criter, xname, iseed, save_state_unif, &
save_state_gauss
if (present(restart)) then
irestart = restart
else
irestart = .false.
endif
if (present(restart_file)) then
isrestart_file = restart_file
else
isrestart_file = 'mo_sce.restart'
endif
if (.not. irestart) then
if (present(parallel)) then
parall = parallel
else
parall = .false.
end if
if (parall) then
! OpenMP or not
n_threads = 1
!$ write(output_unit,*) '--------------------------------------------------'
!$OMP parallel
!$ n_threads = OMP_GET_NUM_THREADS()
!$OMP end parallel
!$ write(output_unit,*) ' SCE is parellel with ',n_threads,' threads'
!$ write(output_unit,*) '--------------------------------------------------'
else
n_threads = 1
end if
! One random number chain per OpenMP thread
allocate(rand_tmp(n_threads))
allocate(iseed(n_threads))
allocate(save_state_unif(n_threads,n_save_state))
allocate(save_state_gauss(n_threads,n_save_state))
if (present(mymask)) then
if (.not. any(mymask)) stop 'mo_sce: all parameters are masked --> none will be optimized'
maskpara = mymask
else
maskpara = .true.
end if
! number of parameters to optimize
nopt = count(maskpara,dim=1)
! total number of parameters
nn = size(pini,1)
! input checking
if (size(prange,dim=1) .ne. size(pini,1)) then
stop 'mo_sce: prange has not matching rows'
end if
if (size(prange,dim=2) .ne. 2) then
stop 'mo_sce: two colums expected for prange'
end if
bl(:) = prange(:,1)
bu(:) = prange(:,2)
do ii=1, nn
if( ((bu(ii)-bl(ii)) .lt. tiny(1.0_dp) ) .and. maskpara(ii) ) then
maskpara(ii) = .false.
write(error_unit,*) 'para #',ii,' :: range = ( ',bl(ii),' , ',bu(ii),' )'
! stop 'mo_sce: inconsistent or too small parameter range'
end if
end do
!
! optionals checking
if (present(mymaxn)) then
if (mymaxn .lt. 2_i4) stop 'mo_sce: maxn has to be at least 2'
maxn = mymaxn
else
maxn = 1000_i8
end if
if (present(mymaxit)) then
maxit = mymaxit
else
maxit = .false.
end if
if (present(mykstop)) then
if (mykstop .lt. 1_i4) stop 'mo_sce: kstop has to be at least 1'
kstop = mykstop
else
kstop = 10_i4
end if
if (present(mypcento)) then
if (mypcento .lt. 0_dp) stop 'mo_sce: pcento should be positive'
pcento = mypcento
else
pcento = 0.0001_dp
end if
if (present(mypeps)) then
if (mypeps .lt. 0_dp) stop 'mo_sce: peps should be positive'
peps = mypeps
else
peps = 0.001_dp
end if
if (present(myseed)) then
if (myseed .lt. 1_i8) stop 'mo_sce: seed should be non-negative'
forall(ii=1:n_threads) iseed(ii) = myseed + (ii-1)*1000_i8
else
call get_timeseed(iseed)
end if
if (present(myngs)) then
if (myngs .lt. 1_i4) stop 'mo_sce: ngs has to be at least 1'
ngs = myngs
else
ngs = 2_i4
end if
if (present(mynpg)) then
if (mynpg .lt. 3_i4) stop 'mo_sce: npg has to be at least 3'
npg = mynpg
else
npg = 2*nopt + 1
end if
if (present(mynps)) then
if (mynps .lt. 2_i4) stop 'mo_sce: nps has to be at least 2'
nps = mynps
else
nps = nopt + 1_i4
end if
if (present(mynspl)) then
if (mynspl .lt. 3_i4) stop 'mo_sce: nspl has to be at least 3'
nspl = mynspl
else
nspl = 2*nopt + 1
end if
if (present(mymings)) then
if (mymings .lt. 1_i4) stop 'mo_sce: mings has to be at least 1'
mings = mymings
else
mings = ngs ! no reduction of complexes
end if
if (present(myiniflg)) then
if ( (myiniflg .ne. 1_i4) .and. (myiniflg .ne. 0_i4)) stop 'mo_sce: iniflg has to be 0 or 1'
iniflg = myiniflg
else
iniflg = 1_i4
end if
idot = .false.
if (present(myprint)) then
if ( (myprint .lt. 0_i4) .or. (myprint .gt. 4_i4)) stop 'mo_sce: iprint has to be between 0 and 4'
if (myprint > 2_i4) then
iprint = myprint-3
idot = .true.
else
iprint = myprint
endif
else
iprint = 2_i4 ! no printing
iprint = 0_i4
end if
if (present(myalpha)) then
alpha = myalpha
else
alpha = 0.8_dp
end if
if (present(mybeta)) then
beta = mybeta
else
beta = 0.45_dp
end if
if (present(popul_file_append)) then
ipopul_file_append = popul_file_append
else
ipopul_file_append = .false.
endif
if (present(tmp_file)) then
itmp_file = .true.
istmp_file = tmp_file
open(999, file=trim(adjustl(istmp_file)), action='write', status = 'unknown')
write(999,*) '# settings :: general'
write(999,*) '# nIterations seed'
write(999,*) maxn, iseed
write(999,*) '# settings :: sce specific'
write(999,*) '# sce_ngs sce_npg sce_nps'
write(999,*) ngs, npg, nps
write(999,*) '# nloop icall ngs1 bestf worstf gnrng (bestx(j),j=1,nn)'
close(999)
else
itmp_file = .false.
istmp_file = ''
end if
if (present(popul_file)) then
ipopul_file = .true.
ispopul_file = popul_file
else
ipopul_file = .false.
ispopul_file = ''
endif
if (ipopul_file .and. (.not. ipopul_file_append)) then
open(999, file=trim(adjustl(ispopul_file)), action='write', status = 'unknown')
write(999,*) '# xf(i) (x(i,j),j=1,nn)'
close(999)
end if
! allocation of arrays
allocate(x(ngs*npg,nn))
allocate(xf(ngs*npg))
allocate(worstx(nn))
allocate(xnstd(nn))
allocate(bound(nn))
allocate(unit(nn))
allocate(criter(kstop+1))
allocate(xname(nn))
allocate(history_tmp(maxn+3*ngs*nspl))
allocate(xtmp(npg,nn))
allocate(ftmp(npg))
if (maxit) then
large = -0.5_dp*huge(1.0_dp)
else
large = 0.5_dp*huge(1.0_dp)
endif
history_tmp(:) = large
criter(:) = large
!
! initialize variables
do ii=1,nn
xname(ii) = compress('p'//num2str(ii,'(I3.3)'))
end do
nloop = 0_i4
loop = 0_i4
igs = 0_i4
!
! compute the total number of points in initial population
npt = ngs * npg
ngs1 = ngs
npt1 = npt
!
if (iprint .lt. 2) then
write(output_unit,*) '=================================================='
write(output_unit,*) 'ENTER THE SHUFFLED COMPLEX EVOLUTION GLOBAL SEARCH'
write(output_unit,*) '=================================================='
end if
!
! Print seed
if (iprint .lt. 2) then
write (*,*) ' Seeds used : ',iseed
end if
! initialize random number generator: Uniform
call xor4096(iseed, rand_tmp, save_state=save_state_unif)
! initialize random number generator: Gaussian
iseed = iseed + 1000_i8
call xor4096g(iseed, rand_tmp, save_state=save_state_gauss)
iseed = 0_i8
!
! compute the bound for parameters being optimized
do ii = 1, nn
bound(ii) = bu(ii) - bl(ii)
unit(ii) = 1.0_dp
end do
!--------------------------------------------------
! compute the function value of the initial point
!--------------------------------------------------
! function evaluation will be counted later...
if (idot) write(output_unit,'(A1)') '.'
if (.not. maxit) then
fpini = functn(pini)
history_tmp(1) = fpini
else
fpini = -functn(pini)
history_tmp(1) = -fpini
end if
! print the initial point and its criterion value
bestx = pini
bestf_tmp = fpini
if (iprint .lt. 2) then
write(output_unit,*) ''
write(output_unit,*) ''
write(output_unit,*) '*** PRINT THE INITIAL POINT AND ITS CRITERION VALUE ***'
call write_best_final()
end if
if (itmp_file) then ! initial tmp file
open(999, file=trim(adjustl(istmp_file)), action='write', position='append', recl=(nn+6)*30)
if (.not. maxit) then
write(999,*) 0,1,ngs1,fpini,fpini,1.0_dp, pini
else
write(999,*) 0,1,ngs1,-fpini,-fpini,1.0_dp, pini
end if
close(999)
end if
!
! generate an initial set of npt1 points in the parameter space
! if iniflg is equal to 1, set x(1,.) to initial point pini(.)
if (iniflg .eq. 1) then
do ii = 1, nn
x(1,ii) = pini(ii)
end do
xf(1) = fpini
!
! else, generate a point randomly and set it equal to x(1,.)
else
itmp = save_state_unif(1,:)
call getpnt(1,bl(1:nn),bu(1:nn),unit(1:nn),pini(1:nn),maskpara,itmp,xx)
save_state_unif(1,:) = itmp
do ii=1, nn
x(1,ii) = xx(ii)
end do
if (idot) write(output_unit,'(A1)') '.'
if (.not. maxit) then
xf(1) = functn(xx)
else
xf(1) = -functn(xx)
end if
end if
!
! count function evaluation of the first point
icall = 1_i8
! if (icall .ge. maxn) return
!
if (parall) then
! -----------------------------------------------------------------------
! Parallel version of complex-loop
! -----------------------------------------------------------------------
! generate npt1-1 random points distributed uniformly in the parameter
! space, and compute the corresponding function values
ithread = 1
!$OMP parallel default(shared) private(ii,jj,ithread,xx)
!$OMP do
do ii = 2, npt1
!$ ithread = OMP_GET_THREAD_NUM() + 1
itmp = save_state_unif(ithread,:)
call getpnt(1,bl(1:nn),bu(1:nn),unit(1:nn),pini(1:nn),maskpara,itmp,xx)
save_state_unif(ithread,:) = itmp
do jj = 1, nn
x(ii,jj) = xx(jj)
end do
if (idot) write(output_unit,'(A1)') '.'
if (.not. maxit) then
xf(ii) = functn(xx)
history_tmp(ii) = xf(ii) ! min(history_tmp(ii-1),xf(ii)) --> will be sorted later
else
xf(ii) = -functn(xx)
history_tmp(ii) = -xf(ii) ! max(history_tmp(ii-1),-xf(ii)) --> will be sorted later
end if
end do
!$OMP end do
!$OMP end parallel
else
! -----------------------------------------------------------------------
! Non-Parallel version of complex-loop
! -----------------------------------------------------------------------
! generate npt1-1 random points distributed uniformly in the parameter
! space, and compute the corresponding function values
ithread = 1
do ii = 2, npt1
call getpnt(1,bl(1:nn),bu(1:nn),unit(1:nn),pini(1:nn),maskpara,save_state_unif(ithread,:),xx)
do jj = 1, nn
x(ii,jj) = xx(jj)
end do
if (idot) write(output_unit,'(A1)') '.'
if (.not. maxit) then
xf(ii) = functn(xx)
icall = icall + 1_i8
history_tmp(icall) = xf(ii) ! min(history_tmp(icall-1),xf(ii)) --> will be sorted later
else
xf(ii) = -functn(xx)
icall = icall + 1_i8
history_tmp(icall) = -xf(ii) ! max(history_tmp(icall-1),-xf(ii)) --> will be sorted later
end if
if (icall .ge. maxn) then
npt1 = ii
if (iprint .lt. 2) then
! maximum trials reached
call write_termination_case(1)
call write_best_final()
end if
exit
end if
end do
end if
icall = int(npt1,i8)
!
! arrange the points in order of increasing function value
if (maxit) then
large = minval(xf(1:npt1))
large = merge(0.9_dp*large, 1.1_dp*large, large>0._dp)
else
large = maxval(xf(1:npt1))
large = merge(1.1_dp*large, 0.9_dp*large, large>0._dp)
endif
xf(1:npt1) = merge(xf(1:npt1), large, is_finite(xf(1:npt1))) ! NaN and Infinite
! sort does not work with NaNs
! -> get history_tmp w/o NaN, sort it, and set the rest to NaN
nonan = size(pack(history_tmp(1:npt1), mask=is_finite(history_tmp(1:npt1))))
if (nonan /= npt1) then
allocate(htmp(nonan))
htmp(1:nonan) = pack(history_tmp(1:npt1), mask=is_finite(history_tmp(1:npt1)))
call sort(htmp(1:nonan))
history_tmp(1:nonan) = htmp(1:nonan)
history_tmp(nonan+1:npt1) = special_value(1.0_dp, 'IEEE_QUIET_NAN')
deallocate(htmp)
else
call sort(history_tmp(1:npt1))
endif
call sort_matrix(x(1:npt1,1:nn),xf(1:npt1))
!
! record the best and worst points
do ii = 1, nn
bestx(ii) = x(1,ii)
worstx(ii) = x(npt1,ii)
end do
bestf_tmp = xf(1)
worstf = xf(npt1)
!
! compute the parameter range for the initial population
call parstt(x(1:npt1,1:nn),bound,peps,maskpara,xnstd,gnrng,ipcnvg)
!
! write currently best x and xf to temporal file
if (itmp_file) then
call write_best_intermediate(.true.)
end if
! write population to file
if (ipopul_file) then
call write_population(.true.)
end if
!
! print the results for the initial population
print: if (iprint .lt. 2) then
! write currently best x and xf to screen
call write_best_intermediate(.false.)
! write population on screen
if (iprint .eq. 1) then
call write_population(.false.)
end if
end if print
!
! Maximum number of function evaluations reached?
if (icall .ge. maxn) then
if (iprint .lt. 2) then
! maximum trials reached
call write_termination_case(1)
call write_best_final()
end if
call set_optional()
! -----------------------
! Abort subroutine
! -----------------------
return
end if
!
! Feasible parameter space converged?
if (ipcnvg .eq. 1) then
if (iprint .lt. 2) then
! converged because feasible parameter space small
call write_termination_case(3)
call write_best_final()
end if
call set_optional()
! -----------------------
! Abort subroutine
! -----------------------
return
end if
!
! transfer all array-like variables in namelist to fixed-size dummy-arrays
dummy_maskpara = .false.
dummy_maskpara(1:size(pini,1)) = maskpara
dummy_bl = -9999.0_dp
dummy_bl(1:size(pini,1)) = bl
dummy_bu = -9999.0_dp
dummy_bu(1:size(pini,1)) = bu
dummy_xx = -9999.0_dp
dummy_xx(1:size(pini,1)) = xx
!
! write restart
open(999, file=isrestart_file, status='unknown', action='write', delim='quote')
write(999, restartnml1)
write(999, restartnml2)
close(999)
endif ! restart or not
if (irestart) then
! read 1st namelist with allocated/scalar variables
open(999, file=isrestart_file, status='old', action='read', delim='quote')
read(999, nml=restartnml1)
close(999)
! transfer all array-like variables in namelist to fixed-size dummy-arrays
maskpara = dummy_maskpara(1:size(pini,1))
bl = dummy_bl(1:size(pini,1))
bu = dummy_bu(1:size(pini,1))
xx = dummy_xx(1:size(pini,1))
! allocate global arrays
allocate(rand_tmp(n_threads))
allocate(iseed(n_threads))
allocate(save_state_unif(n_threads,n_save_state))
allocate(save_state_gauss(n_threads,n_save_state))
allocate(x(ngs*npg,nn))
allocate(xf(ngs*npg))
allocate(worstx(nn))
allocate(xnstd(nn))
allocate(bound(nn))
allocate(unit(nn))
allocate(criter(kstop+1))
allocate(xname(nn))
allocate(history_tmp(maxn+3*ngs*nspl))
allocate(xtmp(npg,nn))
allocate(ftmp(npg))
endif
!
! begin the main loop
mainloop:do while (icall .lt. maxn)
! read variables from restart files
open(999, file=isrestart_file, status='old', action='read', delim='quote')
read(999, nml=restartnml1)
read(999, nml=restartnml2)
close(999)
!
! transfer all array-like variables in namelist to fixed-size dummy-arrays
maskpara = dummy_maskpara(1:size(pini,1))
bl = dummy_bl(1:size(pini,1))
bu = dummy_bu(1:size(pini,1))
xx = dummy_xx(1:size(pini,1))
!
nloop = nloop + 1
!
if (iprint .lt. 2) then
write(output_unit, *) ''
write(output_unit,'(A28,I4)') ' *** Evolution Loop Number ',nloop
end if
!
! begin loop on complexes
! <beta> loop from duan(1993)
if (parall) then
! -----------------------------------------------------------------------
! Parallel version of complex-loop
! -----------------------------------------------------------------------
ithread = 1