program WeightRandomSamplingReservoir

implicit none

integer :: i,j

integer*4 timeArray(3) ! Holds the hour, minute, and second

integer, dimension(100,2) :: V ! Population

real, dimension(10,2) :: R ! Reservoir

integer, dimension(1) :: MinIndex
integer, dimension(2) :: R_dim

real :: ui,ki
real :: ki_min

!----------------------------------------------------------------------------------------!

call itime(timeArray) ! Generate sequence of random number

i = rand ( timeArray(1)+timeArray(2)+timeArray(3) )

R_dim = shape(R)

!----------------------------------------------------------------------------------------!
! Population V: (vi,wi)

write(6,*) '!----- Population -----!'

V = 0
do i = 1, 100
V(i,1) = i
if( i < 90 )then
V(i,2) = 1
else
V(i,2) = 2
end if
write(6,*) 'Data, Weight: ',V(i,1),V(i,2)
end do

!----------------------------------------------------------------------------------------!
! Weighted Random Sampling

R = 9999.0
ki_min = 9999.0

ui = rand(0)
R(1,1) = V(1,1)
R(1,2) = ui**(1.0/V(1,2))
ki_min = R(1,2)
do i = 2, R_dim(1)
ui = rand(0)
R(i,1) = V(i,1)
R(i,2) = ui**(1.0/V(i,2))
if( R(i,2) < ki_min ) ki_min = R(i,2)
end do

do i = R_dim(1)+1, 100
ui = rand(0)
ki = ui**(1.0/V(i,2))

!write(6,*) V(i,1),V(i,2),ki
!write(6,*) (R(j,1), j=1,10)
!write(6,*) (R(j,2), j=1,10)
!write(6,*) minloc(R(:,2)),'dimension',shape(minloc(R)),MinIndex(1)
!write(6,*) "------------"

MinIndex = minloc(R(:,2))

if( ki > ki_min )then
j = MinIndex(1)
R(j,1) = V(i,1)
R(j,2) = ki

MinIndex = minloc(R(:,2))
j = MinIndex(1)
ki_min = R(j,2)
end if

!if( rand(0) < 0.2 )then
!write(6,*) D(i,1),D(i,2)
!end if
end do

!----------------------------------------------------------------------------------------!
! Print results

write(6,*) '!----- Sampling -----!'

do i=1,R_dim(1)

write(6,*) R(i,1), R(i,2)

end do

!----------------------------------------------------------------------------------------!

end program WeightRandomSamplingReservoir


#### Infos

Main author / Auteur principal:
Bazinga

Created / Créé: 2017-05-12

Edited / Modifié: 2017-05-12