-
Notifications
You must be signed in to change notification settings - Fork 0
/
simplex.f90
159 lines (127 loc) · 4.16 KB
/
simplex.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
!
! Simplex redesign
! Alberto Garcia, March 2007
!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------
program simplex
use precision, only: dp
use vars_module, only: nvars, var, read_vars, constrained
use vars_module, only: generate_subs_file,get_subs_file
use toms_minimizer, only: energy=>objective_function,delayed_submit
use m_amoeba, only: amoeba
implicit none
real(dp), dimension(:,:), allocatable, save :: p
real(dp), dimension(:), allocatable, save :: y, x_final, pmin
real(dp) :: ftol ! = 0.0001_dp
real(dp) :: lambda ! = 0.4_dp
real(dp) :: lambda_min, lambda_factor
integer :: mem_stat, iter, itmax, iu, imin, imin_dummy(1)
integer :: i
logical :: file_exists
call read_vars("VARS") ! Sets nvars
allocate(p(nvars+1,nvars),stat=mem_stat)
if (mem_stat /=0 ) stop "alloc error"
allocate(y(nvars+1),stat=mem_stat)
if (mem_stat /=0 ) stop "alloc error"
allocate(x_final(nvars),stat=mem_stat)
if (mem_stat /=0 ) stop "alloc error"
allocate(pmin(nvars),stat=mem_stat)
if (mem_stat /=0 ) stop "alloc error"
!--------------------------------------------------------------------
call io_assign(iu)
open(iu,file="PARAMS",form="formatted",status="old")
read(iu,*) lambda
read(iu,*) lambda_factor
read(iu,*) lambda_min
read(iu,*) itmax
read(iu,*) ftol
close(iu)
!--------------------------------------------------------------------
!
! Initial guess at a random point in the corresponding
! interval, or whatever was specified by the user in the
! VARS file (this should probably be made selectable in
! the PARAMS file).
!
inquire(file="initial.sed", exist=file_exists)
if(file_exists) then
print *, "inital file exists"
call get_subs_file(p(1,:),"initial")
else
print *, "initial file doesn't exist"
do i = 1, nvars
p(1,i) = var(i)%x
enddo
call generate_subs_file(p(1,:),"initial")
endif
!
! do i = 1, nvars
! p(1,i) = var(i)%x
! enddo
! call generate_subs_file(p(1,:),"initial")
do
call generate_subs_file(p(1,:),"best_so_far")
if (lambda < lambda_min) then
print *, "Lambda < ", lambda_min,". Convergence presumably reached"
exit
endif
print *, "Start of an amoeba cycle: [lambda= ", lambda, "]"
!
! Simple-minded handling of constraints
!
do i=1, nvars
p(1+i,:) = p(1,:)
p(1+i,i) = constrained(p(1+i,i) + lambda*var(i)%range,i)
enddo
!$omp parallel do shared(p,y,nvars) private(i)
!This section will be repeated with the periodic restarts
!We can submit and wait for all these to be done
! do i=1, nvars+1
! y(i) = energy(p(i,:))
! enddo
!my edit
print *, "call delayed submit"
call delayed_submit(p,y)
print *, "end delayed submit"
!$omp end parallel do
print *, "Points in simplex and values:"
do i=1, nvars+1
print *, p(i,:), " --- ", y(i)
enddo
print *, "call amoeba"
call amoeba(p,y,nvars+1,nvars,nvars,ftol,energy,iter,req_itmax=itmax)
print *, "end amoeba"
print *, "Amoeba cycle finished in ", iter, " iterations."
print *, "Points in simplex and values:"
do i=1, nvars+1
print *, p(i,:), " --- ", y(i)
enddo
!
! Make sure that the minimum is at the first point, so
! that restarts build the new simplex around the minimum.
!
imin_dummy = minloc(y)
imin = imin_dummy(1)
pmin(:) = p(imin,:) ! Swapping...
p(imin,:) = p(1,:)
p(1,:) = pmin(:)
lambda = lambda * lambda_factor ! Next simplex will be smaller
enddo
!
! Take averages of positions in final simplex
!
print *, "Final (average) point:"
do i=1, nvars
x_final(i) = sum(p(:,i)) / (nvars + 1)
var(i)%x = x_final(i)
print *, "Final (average)", trim(var(i)%name), " :", var(i)%x
enddo
call generate_subs_file(x_final,"final_average")
print *, "Final (actual minimum):"
do i=1, nvars
var(i)%x = pmin(i)
print *, "Final", trim(var(i)%name), " :", var(i)%x
enddo
call generate_subs_file(pmin(:),"final")
deallocate(p,y,x_final,pmin)
End Program simplex