-
Notifications
You must be signed in to change notification settings - Fork 0
/
test_gauss.f90~
142 lines (137 loc) · 4.39 KB
/
test_gauss.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
program test_gauss
implicit none
integer,parameter :: dp=8,n=2,maxiter=10000
character(60) :: task
real(kind=dp) :: f,c_t,x(n),m(n),eps
real(kind=dp) :: Q(n,n),r,xmin(n),fmin
integer :: i,counter(2)
logical :: iprint
task='start'
eps=1e-6
x=3._dp
counter=0
iprint = .true.
if (iprint) print *,'*********************start************************'
do i = 1,maxiter
if (task(1:5).eq.'new_x') call fun(x,n,f)
if (task.eq.'new_x_accepted') then
if (iprint) print *,'accepted'
counter(1)=counter(1)+1
else
if (iprint) print *,'rejected'
counter(2)=counter(2)+1
endif
if (iprint) then
print *,'---------------------------------'
print *,'loop',i
print *,'x',x,'f',f,'ct',c_t
print *,'m',m
print *,'xmin',xmin,'fmin',fmin
print *,'r',r
print *,'C',r**2*matmul(Q,transpose(Q))
print *,'---------------------------------'
endif
call GaussAdapt(n,x,f,m,Q,r,c_t,xmin,fmin,task)
if (task(1:4).eq.'stop') then
print *,task
exit
endif
if (abs(fmin).lt.eps) then
if (iprint) print *,'min found'
exit
endif
enddo
print *,'=============================================================='
print *,'final x',xmin,'final f',fmin
print *,'accepted',counter(1),'rejected',counter(2),'iterations',counter(2)+counter(1)
end program
!----------------------------objective function---------------------------------
subroutine fun(x,n,f)
use random
implicit none
integer,intent(in) :: n
real(kind=dp) :: x(n),a,f
f = (x(1)-1._dp)**2 + (x(2)-3._dp)**2
!a=sqrt(2._dp)
!f = 100*(x(2)-x(1)**2)**2+(a-x(1))
end subroutine fun
!-----------------------Gaussian Adaption Algorithm-----------------------------
subroutine GaussAdapt(n,x,f,m,Q,r,c_t,xmin,fmin,task)
use random
implicit none
character(60),intent(inout) :: task
integer,intent(in) :: n
real(kind=dp),intent(inout) :: f,c_t,x(n),m(n)
real(kind=dp),intent(inout) :: Q(n,n),r
real(kind=dp) :: Q_old(n,n),dC(n,n),eta(n),Id(n,n)
real(kind=dp) :: wm,wc,wt,beta,p,fe,fc,dx(n)
real(kind=dp) :: fmin,xmin(n),eig_val(n),work(3*n-1)
real(kind=dp) :: detQ,dQ(n,n)
integer :: i,j,info
!----------------------------set GauAdapt prameter------------------------------
Id = 0._dp; do i = 1,n; Id(i,i) = 1._dp; end do
wm=exp(1._dp)*n
wc=(n+1._dp)**2/(log(n+1._dp))
wt=exp(1._dp)*n
beta=1._dp/wc
p=1._dp/exp(1._dp)
fe=1._dp + beta*(1._dp-p)
! print *,'f_e',fe
fc=1._dp-beta*p
! print *,'f_c',fc
!-------------------------------------------------------------------------------
if (task(1:5).eq.'new_x') then
if (f .lt. c_t) then
task='new_x_accepted'
!----------------------------save current best----------------------------------
if (f .lt. fmin) then
fmin = f
xmin = x
endif
!------------------------------adapt parameters---------------------------------
r=fe*r
c_T = (1-1/wt)*c_T+1/wt*f
m = (1-1/wm) * m + 1/wm*x
Q_old=Q
!-----------------------------update convariance--------------------------------
call DGER(n,n,1._dp,dx,int(1),dx,int(1),dQ,n)
dC=(1._dp-1._dp/wc)*Id+1._dp/wc*dQ
dQ=dC
call DSYEV( 'V', 'L', n, dQ, n, eig_val, WORK, 3*n-1, INFO )
if (info.ne.0) task='stop_eigendecomp'
! print *,'dQ',dQ
! print *,'eig_val',eig_val
do i=1,n; Id(i,i)=Id(i,i)*sqrt(eig_val(i)); enddo
Q=matmul(dQ,matmul(Id,transpose(dQ)))
!--------------------------------normalize dQ-----------------------------------
detQ=product(eig_val)
detQ=detQ**(1D+0/n)
dQ=dQ/detQ
! print *,detQ,'detQ'
!--------------------------------update Q---------------------------------------
Q=matmul(Q_old,dQ)
else
!---------------rejected lower step size dont adopt cov. + mean-----------------
task='new_x_rejected'
r=fc*r
endif
!------------------------------------sample-------------------------------------
do i=1,n; eta(i)=random_normal(); enddo
! print *,'eta',eta
x=m
call DGEMV('N',n,n,r,Q,n,eta,int(1),1._dp,x,int(1))
if (any(isnan(x))) task='stop_sampling'
! if (r.lt.1e-15) task='stop_r'
endif
if (task(1:5).eq.'start') then
Q = 0._dp; do i = 1,n; Q(i,i) = 1._dp; end do
! print *,'ini C',C
r=1._dp
m = x
task='new_x'
call fun(x,n,f)
c_t=f
fmin=f
xmin=x
endif
end subroutine