-
Notifications
You must be signed in to change notification settings - Fork 4
/
mo_anneal.f90
1563 lines (1356 loc) · 70.6 KB
/
mo_anneal.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
MODULE mo_anneal
! This module is minimizing a cost function via Simulated Annealing
! and is part of the JAMS Fortran library.
! Written Juliane Mai, Mar 2012 : module implementation
! Modified Juliane Mai, May 2012 : anneal: sp version
! Juliane Mai, May 2012 : anneal: documentation
! Juliane Mai, May 2012 : GetTemperature: sp and dp version
! Juliane Mai, Jun 2012 : weighted parameter selection
! Juliane Mai, Aug 2012 : - function anneal instead of subroutine
! - using new module get_timeseed as default for seeding
! - new optional for minimization or maximization
! - fixed parameter ranges possible instead of interface range
! Juliane Mai, Nov 2012 : history of achieved objective function values as optional out
! only in anneal but not anneal_valid
! Juliane Mai, Jan 2013 : - including DDS features in anneal, i.e. reflection at parameter boundaries,
! different parameter selection modes (one, all, neighborhood), and
! different parameter pertubation modes (flexible r=dR (anneal version) or
! constant r=0.2 (dds version))
! - remove sp versions
! - fixed and flexible parameter ranges are now in one function
! using optional arguments
! - undef_funcval instead of anneal_valid function
! Juliane Mai, Feb 2013 : - xor4096 optionals combined in save_state
! Matthias Cuntz, May 2018 : Change in cost function must be better 1e-15 instead of 0.0 to accept
! new state because of numerics of different compilers
! License
! -------
! This file is part of the JAMS Fortran package, distributed under the MIT License.
!
! Copyright (c) 2012-2013 Juliane Mai
!
! 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.
USE mo_kind, ONLY: i4, i8, dp
USE mo_utils, ONLY: le, ge
USE mo_xor4096, ONLY: get_timeseed, xor4096, xor4096g, n_save_state
IMPLICIT NONE
PUBLIC :: anneal ! Minimization of a cost function via Simaulated Annealing
PUBLIC :: GetTemperature ! Function which returns the optimal initial temperature for
! ! a given acceptance ratio and initial parameter set
! ------------------------------------------------------------------
! NAME
!> \brief anneal
! PURPOSE
!> \details Optimizes a user provided cost function using the Simulated Annealing strategy.
! CALLING SEQUENCE
! para = (/ 1.0_dp , 2.0_dp /)
! prange(1,:) = (/ 0.0_dp, 10.0_dp /)
! prange(2,:) = (/ 0.0_dp, 50.0_dp /)
!
! user defined function cost_dp which calculates the cost function value for a
! parameter set (the interface given below has to be used for this function!)
! parabest = anneal(cost_dp, para, prange)
! see also test folder for a detailed example
! INTENT(IN)
!> \param[in] "INTERFACE :: cost_dp" interface calculating the
!> cost function at a given point
!> \param[in] "REAL(DP), DIMENSION(:) :: para" initial parameter set
! INTENT(INOUT)
! None
!
! INTENT(OUT)
! None
!
! INTENT(IN), OPTIONAL
!> \param[in] "REAL(DP), DIMENSION(size(para),2), optional :: prange"
!> lower and upper bound per parameter
! OR
!> \param[in] "INTERFACE, optional :: prange_func" interface calculating the
!> feasible range for a parameter
!> at a certain point, if ranges are variable
!> \param[in] "REAL(DP), optional :: temp" initial temperature \n
!> DEFAULT: Get_Temperature
!> \param[in] "REAL(DP), optional :: DT" geometrical decreement of temperature \n
!> 0.7<DT<0.999 \n
!> DEFAULT: 0.9_dp
!> \param[in] "INTEGER(I4), optional :: nITERmax" maximal number of iterations \n
!> will be increased by 10% if stopping criteria of
!> acceptance ratio or epsilon decreement of cost
!> function is not fullfilled \n
!> DEFAULT: 1000_i4
!> \param[in] "INTEGER(I4), optional :: LEN" Length of Markov Chain \n
!> DEFAULT: MAX(250_i4, size(para,1))
!> \param[in] "INTEGER(I4), optional :: nST" Number of consecutive LEN steps \n
!> DEFAULT: 5_i4
!> \param[in] "REAL(DP), optional :: eps" Stopping criteria of epsilon decreement of
!> cost function \n
!> DEFAULT: 0.0001_dp
!> \param[in] "REAL(DP), optional :: acc" Stopping criteria for Acceptance Ratio \n
!> acc <= 0.1_dp \n
!> DEFAULT: 0.1_dp
!> \param[in] "INTEGER(I4/I8), DIMENSION(3), optional :: seeds"
!> Seeds of random numbers used for random parameter
!> set generation \n
!> DEFAULT: dependent on current time
!> \param[in] "LOGICAL, optional :: printflag"
!> If .true. detailed command line output is written \n
!> DEFAULT: .false.
!> \param[in] "LOGICAL, DIMENSION(size(para)), optional :: maskpara"
!> maskpara(i) = .true. --> parameter is optimized \n
!> maskpara(i) = .false. --> parameter is discarded
!> from optimiztaion \n
!> DEFAULT: .true.
!> \param[in] "REAL(DP), DIMENSION(size(para)), optional :: weight"
!> vector of weights per parameter: \n
!> gives the frequency of parameter to be chosen for
!> optimization (will be scaled to a CDF internally) \n
!> eg. [1,2,1] --> parameter 2 is chosen twice as
!> often as parameter 1 and 2 \n
!> DEFAULT: weight = 1.0_dp
!> \param[in] "INTEGER(I4), optional :: changeParaMode" which and how many param.are changed in one step \n
!> 1 = one parameter \n
!> 2 = all parameter \n
!> 3 = neighborhood parameter \n
!> DEFAULT: 1_i4
!> \param[in] "LOGICAL, optional :: reflectionFlag" if new parameter values are Gaussian
!> distributed and reflected (.true.) or
!> uniform in range (.false.) \n
!> DEFAULT: .false.
!> \param[in] "LOGICAL, optional :: pertubFlexFlag" if pertubation of Gaussian distributed
!> parameter values is constant
!> at 0.2 (.false.) or
!> depends on dR (.true.) \n
!> DEFAULT: .true.
!> \param[in] "LOGICAL, optional :: maxit" maximizing (.true.) or minimizing (.false.)
!> a function \n
!> DEFAULT: .false. (minimization)
!> \param[in] "REAL(DP), optional :: undef_funcval" objective function value defining invalid
!> model output, e.g. -999.0_dp \n
!> \param[in] "character(len=*) , optional :: tmp_file" file with temporal output
! INTENT(INOUT), OPTIONAL
! None
! INTENT(OUT), OPTIONAL
!> \param[out] "REAL(DP), optional :: funcbest" minimized value of cost function
!> \param[out] "REAL(DP), DIMENSION(:,:), ALLOCATABLE, optional
!> :: history" returns a vector of achieved objective
! after ith model evaluation
! RETURN
!> \return real(dp) :: parabest(size(para)) — Parameter set minimizing the cost function.
! RESTRICTIONS
!> \note Either fixed parameter range (prange) OR flexible parameter range (function interface prange_func)
!> has to be given in calling sequence. \n
!> Only double precision version available. \n
!> If single precision is needed not only DP has to be replaced by SP
!> but also I8 of save_state (random number variables) has to be replaced by I4. \n
!> ParaChangeMode > 1 is not applied in GetTemperature.
!> For Temperature estimation always only one single parameter is changed (ParaChangeMode=1)
!> which should give theoretically always the best estimate. \n
!> Cost and prange_func are user defined functions. See interface definition.
! EXAMPLE
! see test_mo_anneal
! LITERATURE
! (1) S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi.
! Optimization by simulated annealing.
! Science, 220:671-680, 1983.
! (2) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.
! Equation of state calculations by fast computing machines.
! J. Chem. Phys., 21:1087-1092, June 1953.
!
! (3) B. A. Tolson, and C. A. Shoemaker.
! Dynamically dimensioned search algorithm for computationally efficient watershed model calibration.
! WRR, 43(1), W01413, 2007.
!
! HISTORY
! Written Samaniego, Jan 2000 : Created
! Samaniego, Mar 2003 : Re-heating
! Juliane Mai, March 2012 : modular version
! Modified Juliane Mai, May 2012 : sp version
! Juliane Mai, May 2012 : documentation
! ------------------------------------------------------------------
INTERFACE anneal
MODULE PROCEDURE anneal_dp
END INTERFACE anneal
! ------------------------------------------------------------------
! NAME
!> \brief GetTemperature
! PURPOSE
!> \details Determines an initial temperature for Simulated Annealing achieving
! certain acceptance ratio.
! CALLING SEQUENCE
! para = (/ 1.0_dp , 2.0_dp /)
! acc_goal = 0.95_dp
! prange(1,:) = (/ 0.0_dp, 10.0_dp /)
! prange(2,:) = (/ 0.0_dp, 50.0_dp /)
!
! user defined function cost_dp which calculates the cost function value for a
! parameter set (the interface given below has to be used for this function!)
! temp = GetTemperature(para, cost_dp, acc_goal, prange)
! see also test folder for a detailed example
! INTENT(IN)
!> \param[in] "REAL(DP), DIMENSION(:) :: paraset" initial (valid) parameter set
!> \param[in] "INTERFACE :: cost_dp" interface calculating the
!> cost function at a given point
!> \param[in] "REAL(DP) :: acc_goal" Acceptance Ratio which has to be achieved
! INTENT(INOUT)
! None
! INTENT(OUT)
! None
! INTENT(IN), OPTIONAL
!> \param[in] "REAL(DP), DIMENSION(size(para),2), optional :: prange"
!> lower and upper bound per parameter
! OR
!> \param[in] "INTERFACE, optional :: prange_func" interface calculating the
!> feasible range for a parameter
!> at a certain point, if ranges are variable
!> \param[in] "INTEGER(I4), optional :: samplesize" number of iterations the estimation of temperature
!> is based on \n
!> DEFAULT: Max(20_i4*n,250_i4)
!> \param[in] "LOGICAL, DIMENSION(size(para)), optional :: maskpara"
!> maskpara(i) = .true. --> parameter is optimized \n
!> maskpara(i) = .false. --> parameter is discarded
!> from optimiztaion \n
!> DEFAULT: .true.
!> \param[in] "INTEGER(I4/I8), DIMENSION(2), optional :: seeds"
!> Seeds of random numbers used for random parameter
!> set generation \n
!> DEFAULT: dependent on current time
!> \param[in] "LOGICAL, optional :: printflag"
!> If .true. detailed command line output is written
!> DEFAULT: .false.
!> \param[in] "REAL(DP), DIMENSION(size(para,1)), optional :: weight"
!> vector of weights per parameter \n
!> gives the frequency of parameter to be chosen for
!> optimization (will be scaled to a CDF internally) \n
!> eg. [1,2,1] --> parameter 2 is chosen twice as
!> often as parameter 1 and 2 \n
! DEFAULT: weight = 1.0_dp
!> \param[in] "LOGICAL, optional :: maxit" minimizing (.false.) or maximizing (.true.)
!> a function \n
!> DEFAULT: .false. (minimization)
!> \param[in] "REAL(DP), optional :: undef_funcval" objective function value defining invalid
!> model output, e.g. -999.0_dp
! INTENT(INOUT), OPTIONAL
! None
! INTENT(OUT), OPTIONAL
! None
! RETURN
!> \return real(dp) :: temperature — Temperature achieving a certain
!> acceptance ratio in Simulated Annealing
! RESTRICTIONS
!> \note Either fixed parameter range (prange) OR flexible parameter range (function interface prange_func)
!> has to be given in calling sequence. \n
!> Only double precision version available.
!> If single precision is needed not only DP has to be replaced by SP
!> but also I8 of save_state (random number variables) has to be replaced by I4. \n
!> ParaChangeMode > 1 is not applied in GetTemperature.
!> For Temperature estimation always only one single parameter is changed (ParaChangeMode=1)
!> which should give theoretically always the best estimate. \n
!> Cost and prange_func are user defined functions. See interface definition.
! EXAMPLE
! see test_mo_anneal
! LITERATURE
! (1) Walid Ben-Ameur
! Compututing the Initial Temperature of Simulated Annealing
! Comput. Opt. and App. (2004).
! HISTORY
! Written Juliane Mai, May 2012
! ------------------------------------------------------------------
INTERFACE GetTemperature
MODULE PROCEDURE GetTemperature_dp
END INTERFACE GetTemperature
PRIVATE
INTERFACE Generate_Neighborhood_weight
MODULE PROCEDURE Generate_Neighborhood_weight_dp
END INTERFACE Generate_Neighborhood_weight
! ------------------------------------------------------------------
CONTAINS
FUNCTION anneal_dp( cost, para, & ! obligatory
prange, prange_func, & ! optional IN: exactly one of both
temp, Dt, nITERmax, Len, nST, eps, acc, & ! optional IN
seeds, printflag, maskpara, weight, & ! optional IN
changeParaMode, reflectionFlag, & ! optional IN
pertubFlexFlag, maxit, undef_funcval, & ! optional IN
tmp_file, & ! optional IN
funcbest, history & ! optional OUT
) &
result(parabest)
IMPLICIT NONE
INTERFACE
FUNCTION cost(paraset)
! calculates the cost function at a certain parameter set paraset
use mo_kind
REAL(DP), DIMENSION(:), INTENT(IN) :: paraset
REAL(DP) :: cost
END FUNCTION cost
END INTERFACE
INTERFACE
SUBROUTINE prange_func(paraset,iPar,rangePar)
! gives the range (min,max) of the parameter iPar at a certain parameter set paraset
use mo_kind
REAL(DP), DIMENSION(:), INTENT(IN) :: paraset
INTEGER(I4), INTENT(IN) :: iPar
REAL(DP), DIMENSION(2), INTENT(OUT) :: rangePar
END SUBROUTINE prange_func
END INTERFACE
optional :: prange_func
REAL(DP), DIMENSION(:), INTENT(IN) :: para
! ! initial parameter
REAL(DP), DIMENSION(size(para,1)) :: parabest
! ! parameter set minimizing objective
! optionals
REAL(DP), OPTIONAL, DIMENSION(size(para,1),2), INTENT(IN) :: prange
! ! lower and upper limit per parameter
REAL(DP), OPTIONAL, INTENT(IN) :: temp
! ! starting temperature
! ! (DEFAULT: Get_Temperature)
REAL(DP), OPTIONAL, INTENT(IN) :: Dt
! ! geometrical decreement, 0.7<DT<0.999
! ! (DEFAULT: 0.9)
INTEGER(I4), OPTIONAL, INTENT(IN) :: nITERmax
! ! maximal number of iterations
! ! (DEFAULT: 1000)
INTEGER(I4), OPTIONAL, INTENT(IN) :: Len
! ! Length of Markov Chain,
! ! DEFAULT: max(250, size(para,1))
INTEGER(I4), OPTIONAL, INTENT(IN) :: nST
! ! Number of consecutive LEN steps
! ! (DEFAULT: 5)
REAL(DP), OPTIONAL, INTENT(IN) :: eps
! ! epsilon decreement of cost function
! ! (DEFAULT: 0.01)
REAL(DP), OPTIONAL, INTENT(IN) :: acc
! ! Acceptance Ratio, <0.1 stopping criteria
! ! (DEFAULT: 0.1)
INTEGER(I8), OPTIONAL, INTENT(IN) :: seeds(3)
! ! Seeds of random numbers
! ! (DEFAULT: Get_timeseed)
LOGICAL, OPTIONAL, INTENT(IN) :: printflag
! ! If command line output is written (.true.)
! ! (DEFAULT: .false.)
LOGICAL, OPTIONAL, DIMENSION(size(para,1)), INTENT(IN) :: maskpara
! ! true if parameter will be optimized
! ! false if parameter is discarded in optimization
! ! (DEFAULT: .true.)
REAL(DP), OPTIONAL, DIMENSION(size(para,1)), INTENT(IN) :: weight
! ! vector of weights per parameter
! ! gives the frequency of parameter to be
! ! chosen for optimization
! ! (DEFAULT: uniform)
INTEGER(I4), OPTIONAL, INTENT(IN) :: changeParaMode
! ! which and how many param. are changed in a step
! ! 1 = one parameter
! ! 2 = all parameter
! ! 3 = neighborhood parameter
! ! (DEFAULT: 1_i4)
LOGICAL, OPTIONAL, INTENT(IN) :: reflectionFlag
! ! if new parameter values are selected normal
! ! distributed and reflected (.true.) or
! ! uniform in range (.false.)
! ! (DEFAULT: .false.)
LOGICAL, OPTIONAL, INTENT(IN) :: pertubFlexFlag
! ! if pertubation of normal distributed parameter
! ! values is constant 0.2 (.false.) or
! ! depends on dR (.true.)
! ! (DEFAULT: .true.)
LOGICAL, OPTIONAL, INTENT(IN) :: maxit
! ! Maximization or minimization of function
! ! maximization = .true., minimization = .false.
! ! (DEFAULT: .false.)
REAL(DP), OPTIONAL, INTENT(IN) :: undef_funcval
! ! objective function value occuring if
! ! parameter set leads to invalid model results,
! ! e.g. -9999.0_dp
! ! (DEFAULT: not present)
CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tmp_file
! ! file for temporal output
REAL(DP), OPTIONAL, INTENT(OUT) :: funcbest
! ! minimized value of cost function
! ! (DEFAULT: not present)
REAL(DP), OPTIONAL, DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) :: history
! ! returns a vector of achieved objective
! ! after ith model evaluation
! ! (DEFAULT: not present)
! local variables
INTEGER(I4) :: n ! Number of parameters
integer(I4) :: iPar, par ! counter for parameter
real(DP), DIMENSION(2) :: iParRange ! parameter's range
REAL(DP) :: T_in ! Temperature
REAL(DP) :: DT_IN ! Temperature decreement
INTEGER(I4) :: nITERmax_in ! maximal number of iterations
INTEGER(I4) :: LEN_IN ! Length of Markov Chain
INTEGER(I4) :: nST_in ! Number of consecutive LEN steps
REAL(DP) :: eps_in ! epsilon decreement of cost function
REAL(DP) :: acc_in ! Acceptance Ratio, <0.1 stopping criteria
INTEGER(I8) :: seeds_in(3) ! Seeds of random numbers
LOGICAL :: printflag_in ! If command line output is written
LOGICAL :: coststatus ! Checks status of cost function value,
! ! i.e. is parameter set is feasible
LOGICAL, DIMENSION(size(para,1)) :: maskpara_in ! true if parameter will be optimized
INTEGER(I4), DIMENSION(:), ALLOCATABLE :: truepara ! indexes of parameters to be optimized
REAL(DP), DIMENSION(size(para,1)) :: weight_in ! CDF of parameter chosen for optimization
REAL(DP), DIMENSION(size(para,1)) :: weightUni ! uniform CDF of parameter chosen for optimization
REAL(DP), DIMENSION(size(para,1)) :: weightGrad ! linear combination of weight and weightUni
INTEGER(I4) :: changeParaMode_inin ! 1=one, 2=all, 3=neighborhood
LOGICAL :: reflectionFlag_inin ! true=gaussian distributed and reflected,
! ! false=uniform distributed in range
LOGICAL :: pertubFlexFlag_inin ! variance of normal distributed parameter values
! ! .true. = depends on dR, .false. = constant 0.2
REAL(DP) :: maxit_in ! maximization = -1._dp, minimization = 1._dp
REAL(DP) :: costbest ! minimized value of cost function
REAL(DP), DIMENSION(:,:), ALLOCATABLE :: history_out ! best cost function value after k iterations
REAL(DP), DIMENSION(:,:), ALLOCATABLE :: history_out_tmp ! helper vector
type paramLim
real(DP) :: min ! minimum value
real(DP) :: max ! maximum value
real(DP) :: new ! new state value
real(DP) :: old ! old state value
real(DP) :: best ! best value found
real(DP) :: dMult ! sensitivity multiplier for parameter search
end type paramLim
type (paramLim), dimension (size(para,1)), target :: gamm ! Parameter
! for random numbers
real(DP) :: RN1, RN2, RN3 ! Random numbers
integer(I8), dimension(n_save_state) :: save_state_1 ! optional arguments for restarting RN stream 1
integer(I8), dimension(n_save_state) :: save_state_2 ! optional arguments for restarting RN stream 2
integer(I8), dimension(n_save_state) :: save_state_3 ! optional arguments for restarting RN stream 3
! for dds parameter selection
logical, dimension(size(para,1)) :: neighborhood ! selected parameter in neighborhood
real(dp) :: pertubationR ! neighborhood pertubation size parameter
! for SA
integer(I4) :: idummy,i
real(DP) :: NormPhi
real(DP) :: ac_ratio, pa
integer(i4), dimension(:,:), allocatable :: iPos_iNeg_history ! stores iPos and iNeg of nST last Markov Chains
real(DP) :: fo, fn, df, fBest
real(DP) :: rho, fInc, fbb, dr
real(DP) :: T0, DT0
real(DP), parameter :: small = -700._DP
integer(I4) :: j, iter, kk
integer(I4) :: Ipos, Ineg
integer(I4) :: iConL, iConR, iConF
integer(I4) :: iTotalCounter ! includes reheating for final conditions
integer(I4) :: iTotalCounterR ! counter of interations in one reheating
logical :: iStop
logical :: ldummy
! CHECKING OPTIONALS
n = size(para,1)
! either range or rangfunc has to be given
if ( present(prange) .eqv. present(prange_func) ) then
stop 'anneal: Either range or prange_func has to be given'
end if
if (present(Dt)) then
if ( (Dt .lt. 0.7_dp) .or. (Dt .gt. 0.999_dp) ) then
stop 'Input argument DT must lie between 0.7 and 0.999'
else
DT_IN = Dt
end if
else
DT_IN = 0.9_dp
endif
if (present(nITERmax)) then
if (nITERmax .lt. 1_I4) then
stop 'Input argument nITERmax must be greater 0'
else
nITERmax_in = nITERmax
end if
else
nITERmax_in = 1000_i4
endif
if (present(Len)) then
if (Len .lt. Max(20_i4*n,250_i4)) then
write(*,*) 'WARNING: Input argument LEN should be greater than Max(250,20*N), N=number of parameters'
LEN_IN = Len
else
LEN_IN = Len
end if
else
LEN_IN = Max(20_i4*n,250_i4)
endif
idummy = nITERmax_in/LEN_IN+1_i4
allocate(history_out(idummy,2))
if (present(nST)) then
if (nST .lt. 0_i4) then
stop 'Input argument nST must be greater than 0'
else
nST_in = nST
end if
else
nST_in = 5_i4
endif
allocate(iPos_iNeg_history(nST_in,2))
iPos_iNeg_history = 0_i4
if (present(eps)) then
if ( le(eps, 0.0_dp) ) then
stop 'Input argument eps must be greater than 0'
else
eps_in = eps
end if
else
eps_in = 0.0001_dp
endif
if (present(acc)) then
if ( le(acc, 0.0_dp) .or. ge(acc, 1.0_dp) ) then
stop 'Input argument acc must lie between 0.0 and 1.0'
else
acc_in = acc
end if
else
acc_in = 0.1_dp
endif
if (present(seeds)) then
seeds_in = seeds
else
! Seeds depend on actual time
call get_timeseed(seeds_in)
endif
if (present(printflag)) then
printflag_in = printflag
else
printflag_in = .false.
endif
if (present(maskpara)) then
if (count(maskpara) .eq. 0_i4) then
stop 'Input argument maskpara: At least one element has to be true'
else
maskpara_in = maskpara
end if
else
maskpara_in = .true.
endif
if (present(maxit)) then
ldummy = maxit
if (maxit) then
maxit_in = -1._dp
else
maxit_in = 1._dp
end if
else
ldummy = .false.
maxit_in = 1._dp
endif
if (present(temp)) then
if ( (temp .lt. 0.0_dp) ) then
stop 'Input argument temp must be greater then zero'
else
T_in = temp
end if
else
if (present(prange_func)) then
if (present(undef_funcval)) then
T_in = GetTemperature( para, cost, 0.95_dp, prange_func=prange_func, &
maskpara=maskpara_in, samplesize=2_i4*LEN_IN, &
seeds=seeds_in(1:2), printflag=printflag_in, &
maxit=ldummy,undef_funcval=undef_funcval)
else
T_in = GetTemperature( para, cost, 0.95_dp, prange_func=prange_func, &
maskpara=maskpara_in, samplesize=2_i4*LEN_IN, &
seeds=seeds_in(1:2), printflag=printflag_in, &
maxit=ldummy)
end if
else
if (present(undef_funcval)) then
T_in = GetTemperature( para, cost, 0.95_dp, prange=prange, &
maskpara=maskpara_in, samplesize=2_i4*LEN_IN, &
seeds=seeds_in(1:2), printflag=printflag_in, &
maxit=ldummy,undef_funcval=undef_funcval)
else
T_in = GetTemperature( para, cost, 0.95_dp, prange=prange, &
maskpara=maskpara_in, samplesize=2_i4*LEN_IN, &
seeds=seeds_in(1:2), printflag=printflag_in, &
maxit=ldummy)
end if
end if
end if
if (present(changeParaMode)) then
changeParaMode_inin = changeParaMode
else
changeParaMode_inin = 1_i4
end if
if (present(reflectionFlag)) then
reflectionFlag_inin = reflectionFlag
else
reflectionFlag_inin = .false.
end if
if (present(pertubFlexFlag)) then
pertubFlexFlag_inin = pertubFlexFlag
else
pertubFlexFlag_inin = .true.
end if
! Temporal file writing
if(present(tmp_file)) then
open(unit=999,file=trim(adjustl(tmp_file)), action='write', status = 'unknown')
write(999,*) '# settings :: general'
write(999,*) '# nIterations iseed'
write(999,*) nITERmax_in, seeds_in
write(999,*) '# settings :: anneal specific'
write(999,*) '# sa_tmp'
write(999,*) T_in
write(999,*) '# iter bestf (bestx(j),j=1,nopt)'
close(999)
end if
! INITIALIZATION
allocate ( truepara(count(maskpara_in)) )
idummy = 0_i4
do i=1,N
if ( maskpara_in(i) ) then
idummy = idummy+1_i4
truepara(idummy) = i
end if
end do
if (printflag_in) then
write(*,*) 'Following parameters will be optimized: ',truepara
end if
weight_in = 0.0_dp
weightUni = 0.0_dp
if (present(weight)) then
where ( maskpara_in(:) )
weight_in(:) = weight(:)
weightUni(:) = 1.0_dp
end where
else
where ( maskpara_in(:) )
weight_in(:) = 1.0_dp
weightUni(:) = 1.0_dp
end where
endif
! scaling the weights
weight_in = weight_in/sum(weight_in)
weightUni = weightUni/sum(weightUni)
! cummulating the weights
do i=2,n
weight_in(i) = weight_in(i) + weight_in(i-1)
weightUni(i) = weightUni(i) + weightUni(i-1)
end do
call xor4096 (seeds_in(1), RN1, save_state=save_state_1)
call xor4096 (seeds_in(2), RN2, save_state=save_state_2)
call xor4096g(seeds_in(3), RN3, save_state=save_state_3)
seeds_in = 0_i8
! Start Simulated Annealing routine
gamm(:)%dmult = 1.0_dp
gamm(:)%new = para(:)
gamm(:)%old = para(:)
gamm(:)%best = para(:)
NormPhi = -9999.9_dp
T0= T_in
DT0= DT_IN
! Generate and evaluate the initial solution state
fo = cost(gamm(:)%old) * maxit_in
if ( abs(fo) .lt. tiny(0.0_dp) ) fo = 0.0000001_dp * maxit_in
file_write: if (present(tmp_file)) then
open(unit=999,file=trim(adjustl(tmp_file)), action='write', position='append', recl=(n+2)*30)
if (.not. ldummy) then
write(999,*) '0', fo, gamm(:)%old
else
write(999,*) '0', -fo, gamm(:)%old
end if
close(999)
end if file_write
! initialize counters /var for new SA
iConL= 0_i4
iConR= 1_i4
iConF= 0_i4
iTotalCounter= 0_i4
iTotalCounterR = 0_i4
fn = 0.0_DP
ac_ratio = 1.0_DP
iStop= .TRUE.
iter = 0_i4
! restoring initial T param.
DT_IN= DT0
! Storing the best solution so far
NormPhi = fo * maxit_in
fo = fo/NormPhi
fBest = 1.0_dp * maxit_in
if (printflag_in) then
print '(A15,E15.7,A4,E15.7,A4)', ' start NSe = ', fBest*maxit_in , ' ( ',fBest*normPhi*maxit_in,' ) '
print '(I8, 2I5, 4E15.7)', 1_i4, 0_i4, 0_i4, 1._dp, T_in, fo*normPhi*maxit_in, fBest*normPhi*maxit_in
end if
! ****************** Stop Criterium *******************
! Repeat until the % reduction of nST consecutive Markov
! chains (LEN) of the objective function (f) <= epsilon
! ******************************************************
loopTest: do while (iStop)
iter = iter + 1_i4
Ipos= 0_i4
Ineg= 0_i4
fbb= fBest
!LEN_IN = int( real(iter,dp)*sqrt(real(iter,dp)) / nITER + 1.5*real(N,dp),i4)
! Repeat LEN times with feasible solution
j=1
loopLEN: do while (j .le. LEN_IN)
iTotalCounterR = iTotalCounterR + 1_i4
iTotalCounter = iTotalCounter + 1_i4
! Generate a random subsequent state and evaluate its objective function
! (1) Generate new parameter set
dR=(1.0_DP - real(iTotalCounterR,dp) / real(nITERmax_in,dp))**2.0_DP
if ( .not. ( ge(dR , 0.05_DP) ) .and. &
(iTotalCounterR <= int(real(nIterMax_in,dp)/3._dp*4_dp,i4))) then
dR = 0.05_DP
end if
if (pertubFlexFlag_inin) then
pertubationR = max(dR/5.0_dp,0.05_dp)
else
pertubationR = 0.2_dp
end if
! gradual weights:
! acc_ratio close to 1 --> weight
! acc_ratio close to 0 --> weight uniform
! gradual weighting is linear combination of these two weights
weightGrad(:) = weightUni(:)+(ac_ratio)*(weight_in(:)-weightUni(:))
select case(changeParaMode_inin)
case(1_i4) ! only one parameter is changed
! (1a) Change one parameter traditionally
call xor4096(seeds_in(2), RN2, save_state=save_state_2)
iPar=1_i4
!
do while (weightGrad(iPar) .lt. RN2)
iPar = iPar + 1_i4
end do
if (present(prange_func)) then
call prange_func(gamm(:)%old, iPar, iParRange)
else
iParRange(1) = prange(iPar,1)
iParRange(2) = prange(iPar,2)
endif
gamm(iPar)%min = iParRange(1)
gamm(iPar)%max = iParRange(2)
if (reflectionFlag_inin) then
call xor4096g(seeds_in(3), RN3, save_state=save_state_3)
gamm(iPar)%new = parGen_dds_dp( gamm(iPar)%old, pertubationR, &
gamm(iPar)%min, gamm(iPar)%max,RN3)
else
call xor4096(seeds_in(2), RN2, save_state=save_state_2)
gamm(iPar)%new = parGen_anneal_dp( gamm(iPar)%old, dR, &
gamm(iPar)%min, gamm(iPar)%max,RN2)
end if
case(2_i4) ! all parameter are changed
do par=1, size(truepara)
iPar = truepara(par)
if (present(prange_func)) then
call prange_func(gamm(:)%old, iPar, iParRange)
else
iParRange(1) = prange(iPar,1)
iParRange(2) = prange(iPar,2)
endif
gamm(iPar)%min = iParRange(1)
gamm(iPar)%max = iParRange(2)
if (reflectionFlag_inin) then
call xor4096g(seeds_in(3),RN3, save_state=save_state_3)
gamm(iPar)%new = parGen_dds_dp( gamm(iPar)%old, pertubationR, &
gamm(iPar)%min, gamm(iPar)%max,RN3)
else
call xor4096(seeds_in(2),RN2, save_state=save_state_2)
gamm(iPar)%new = parGen_anneal_dp( gamm(iPar)%old, dR, &
gamm(iPar)%min, gamm(iPar)%max,RN2)
end if
end do
case(3_i4) ! parameter in neighborhood are changed
! Generate new neighborhood
call generate_neighborhood_weight_dp( truepara, weightGrad, save_state_1, &
iTotalCounter, nIterMax_in, neighborhood)
!
! change parameter in neighborhood
do iPar=1, n
if (neighborhood(iPar)) then
! find range of parameter
if (present(prange_func)) then
call prange_func(gamm(:)%old, iPar, iParRange)
else
iParRange(1) = prange(iPar,1)
iParRange(2) = prange(iPar,2)
endif
gamm(iPar)%min = iParRange(1)
gamm(iPar)%max = iParRange(2)
!
if (reflectionFlag_inin) then
! generate gaussian distributed new parameter value which is reflected if out of bound
call xor4096g(seeds_in(3),RN3, save_state=save_state_3)
gamm(iPar)%new = parGen_dds_dp( gamm(iPar)%old, pertubationR, &
gamm(iPar)%min, gamm(iPar)%max,RN3)
else
! generate new parameter value uniform distributed in range (no reflection)
call xor4096(seeds_in(2),RN2, save_state=save_state_2)
gamm(iPar)%new = parGen_anneal_dp( gamm(iPar)%old, dR, &
gamm(iPar)%min, gamm(iPar)%max,RN2)
end if
end if
end do
end select
! (2) Calculate new objective function value
fn = cost(gamm(:)%new) * maxit_in
coststatus = .true.
if (present(undef_funcval)) then
if ( abs(fn*maxit_in-undef_funcval) .lt. tiny(1.0_dp) ) then
coststatus = .false.
end if
end if
feasible: if (coststatus) then ! feasible parameter set
fn = fn/normPhi
! Change in cost function value
df = fn-fo
! analyze change in the objective function: df
!if (df < 0.0_DP) then
if (df < (-1.e-15)) then
! accept the new state
Ipos = Ipos+1_i4
fo = fn
gamm(:)%old = gamm(:)%new
! keep best solution
if (fo < fBest) then
fBest = fo
gamm(:)%best = gamm(:)%new
endif
else
if ( df > eps_in ) then
rho = -df/T_in
if (rho < small) then
pa = 0.0_DP
else
pa = EXP(rho)
end if
!
call xor4096(seeds_in(1), RN1, save_state=save_state_1)
!
if (pa > RN1) then
! accept new state with certain probability
Ineg = Ineg+1_i4
fo = fn
! save old state
gamm(:)%old = gamm(:)%new
end if
end if
end if
j=j+1
else
! function value was not valid
iTotalCounterR = iTotalCounterR - 1_i4
iTotalCounter = iTotalCounter - 1_i4
end if feasible !valid parameter set
end do loopLEN
! estimate acceptance ratio
ac_ratio=(real(Ipos,dp) + real(Ineg,dp))/real(LEN_IN,dp)
if (modulo(iTotalCounter,LEN_IN*nST_in)/LEN_IN .gt. 0_i4) then
idummy = modulo(iTotalCounter,LEN_IN*nST_in)/LEN_IN
else
idummy = nST_in
end if
iPos_iNeg_history(idummy,:) = (/ iPos, iNeg /)
! store current best in history vector
history_out(iTotalCounter/LEN_IN,1) = real(iTotalCounter,dp)
history_out(iTotalCounter/LEN_IN,2) = maxit_in * fBest*normPhi
file_write2: if (present(tmp_file)) then
open(unit=999,file=trim(adjustl(tmp_file)), action='write', position='append',recl=(n+2)*30)
write(999,*) iTotalCounter, maxit_in * fBest*normPhi, gamm(:)%best
close(999)
end if file_write2
if (printflag_in) then
print '(I8, 2I5, E15.7, 3E15.7)', ITotalCounter, Ipos, Ineg, ac_ratio, T_in, &
fo*NormPhi*maxit_in, fBest*NormPhi*maxit_in
end if