Skip to content

Commit

Permalink
Merge pull request #24 from siadicus/develop
Browse files Browse the repository at this point in the history 
Develop
  • Loading branch information
@jwhite-usgs
jwhite-usgs authored Sep 8, 2019
2 parents 9c289cd + 38deefe commit 50b66ed
Show file tree
Hide file tree
Showing 14 changed files with 2,518 additions and 550 deletions.
1 change: 1 addition & 0 deletions src/programs/pestpp-pso/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ OBJECTS := \
readsubs \
rmsubs \
utilsubs \
parunc \
writesubs
OBJECTS := $(addsuffix $(OBJ_EXT),$(OBJECTS))

Expand Down
192 changes: 162 additions & 30 deletions src/programs/pestpp-pso/bestsubs.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,16 @@ subroutine getgbest(gindex,gbest,gmbest,gpbest)
! specifications:
!----------------------------------------------------------------------------------------
integer,intent(out)::gindex
integer::i,j,ipart,left,right,ineibr,iineibr
integer::i,j,ipart,left,right,ineibr,iineibr,ifeas

double precision,intent(out)::gbest,gmbest,gpbest
double precision::gnbest
!----------------------------------------------------------------------------------------

if (neibr == 1) then
!
if (nptocon > 0) call killrm('Neighborhoods not yet supported with constraints')
!
! find gbest within each neighborhood
gnbest = 0.0d+00
!
Expand Down Expand Up @@ -60,29 +62,89 @@ subroutine getgbest(gindex,gbest,gmbest,gpbest)
!
end if
!
! find gbest of the swarm
!
i = 0
!
do ipart=1,npop
! if using constraints, check for feasibility
if (nptocon > 0) then
!
i = i + 1
ifeas = 0
i = 0
!
do ipart=1,npop
if (feas(ipart) == 1) ifeas = 1
end do
!
if (objopt(ipart) < gbest .or. i == 1) then
gbest = objopt(ipart)
gmbest = objmopt(ipart)
gpbest = objpopt(ipart)
gindex = ipart
! find gbest of the swarm:
! at least one pbest position is feasible
if (ifeas == 1) then
!
do ipart=1,npop
!
if (vioopt(ipart) > 1.0d-16) then
!
cycle
!
else
!
i = i + 1
!
if (objopt(ipart) < gbest .or. i == 1) then
gbest = objopt(ipart)
gindex = ipart
end if
!
end if
!
end do
!
! no pbest position is feasible
else
!
do ipart=1,npop
!
if (vioopt(ipart) < gbest .or. ipart == 1) then
gbest = vioopt(ipart)
gindex = ipart
end if
!
end do
!
end if
!
end do
else
!
do ipart=1,npop
!
if (objopt(ipart) < gbest .or. ipart == 1) then
gbest = objopt(ipart)
gindex = ipart
end if
!
end do
!
end if
!
! ! find gbest of the swarm
! !
! i = 0
! !
! do ipart=1,npop
! !
! i = i + 1
! !
! if (objopt(ipart) < gbest .or. i == 1) then
! gbest = objopt(ipart)
! gmbest = objmopt(ipart)
! gpbest = objpopt(ipart)
! gindex = ipart
! end if
! !
! end do


end subroutine getgbest



subroutine getpbest()
subroutine getpbest(ipart)
!========================================================================================
!==== This subroutine determines the personal best particle position. ====
!==== by Adam Siade ====
Expand All @@ -94,33 +156,103 @@ subroutine getpbest()
implicit none

! specifications:
!----------------------------------------------------------------------------------------
integer::ipart,iparm,igp
!----------------------------------------------------------------------------------------
integer,intent(in)::ipart
integer::iparm,igp
!----------------------------------------------------------------------------------------

! if running standard PSO, set pbest based on composite objective function
do ipart=1,npop
if (obj(ipart) < objopt(ipart)) then
!
if (obj(ipart) < objopt(ipart)) then
objopt(ipart) = obj(ipart)
objmopt(ipart) = objm(ipart)
objpopt(ipart) = objp(ipart)
!
do igp=1,nobsgp
objgpopt(ipart,igp) = objgp(ipart,igp)
end do
!
do iparm=1,npar
pbest(ipart,iparm) = partval(ipart,iparm)
end do
!
end if
!
end subroutine getpbest



subroutine pbestcon(ipart)
!========================================================================================
!==== This subroutine determines the personal best particle positions for either ====
!==== constrained or unconstrained MOPSO. ====
!==== by Adam Siade ====
!========================================================================================
!========================================================================================

use psodat

implicit none

! specifications:
!----------------------------------------------------------------------------------------
integer,intent(in)::ipart
integer::iparm,ipto,igp
!----------------------------------------------------------------------------------------

if (nptocon > 0) then
!
if (feas(ipart) == 0) then
!
objopt(ipart) = obj(ipart)
objmopt(ipart) = objm(ipart)
objpopt(ipart) = objp(ipart)
! particle has never been feasible, set pbest based on violate
if (violate(ipart) < vioopt(ipart)) then
!
vioopt(ipart) = violate(ipart)
!
do iparm=1,npar
pbest(ipart,iparm) = partval(ipart,iparm)
end do
!
objopt(ipart) = obj(ipart)
!
end if
!
do igp=1,nobsgp
objgpopt(ipart,igp) = objgp(ipart,igp)
end do
else
!
do iparm=1,npar
pbest(ipart,iparm) = partval(ipart,iparm)
end do
! particle has been feasible before, but is not right now - do nothing
if (violate(ipart) > 1.0d-16) return
!
! particle is feasible now for the first time, update pbest and zero vioopt
if (vioopt(ipart) > 1.0d-16) then
!
vioopt(ipart) = 0.0d+00
!
do iparm=1,npar
pbest(ipart,iparm) = partval(ipart,iparm)
end do
!
objopt(ipart) = obj(ipart)
!
do igp=1,nobsgp
objgpopt(ipart,igp) = objgp(ipart,igp)
end do
!
! particle is feasible now and has been before, update pbest as usual
else
!
call getpbest(ipart)
!
end if
!
end if
!
end do

else
!
call killrm('Something is wrong with the declaration of constraints')
!
end if

end subroutine getpbest
end subroutine pbestcon



Expand Down
70 changes: 33 additions & 37 deletions src/programs/pestpp-pso/estimation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,23 +7,11 @@ subroutine psoest(basnam)
use psodat

implicit none

! specifications:
!----------------------------------------------------------------------------------------
!----------------------------------------------------------------------------------------
! interfaces
!----------------------------------------------------------------------------------------
interface
subroutine readrst(basnam)
integer::ipart,iparm,irep,iobs,igp,ipto
integer,dimension(:),allocatable::measgp

character(len=100), optional, intent(in)::basnam
character(len=100)::scrc,fnam
end subroutine readrst
end interface
!----------------------------------------------------------------------------------------
! external routines for run management via PANTHER
! external routines for run management via YAMR
!----------------------------------------------------------------------------------------
external rmif_run
integer rmif_run
Expand All @@ -44,8 +32,6 @@ end subroutine readrst
double precision::gbest,gmbest,gpbest,objmin

character(len=100),intent(in)::basnam
integer :: n_seed
integer, dimension(:), allocatable :: a_seed
!----------------------------------------------------------------------------------------
!----------------------------------------------------------------------------------------

Expand Down Expand Up @@ -75,13 +61,8 @@ end subroutine readrst
! write initial conditions to record file
call listini(basnam,gindex)
!
!-- generate random parameter sets and velocities, set pbest, and add corresponding runs to the queue
call random_seed (size = n_seed)
allocate(a_seed(1:n_seed))
a_seed = iseed
call random_seed(put = a_seed)
deallocate(a_seed)

!-- generate random parameter sets and velocities, set pbest, and add corresponding runs to the queue
call random_seed()
!
inpar = 0
!
Expand Down Expand Up @@ -148,8 +129,7 @@ end subroutine readrst
modfail(ipart) = 0
!
! get model run results
irun = ipart-1
call modelrm(0,irun,fail)
call modelrm(0,ipart-1,fail)
!
modfail(ipart) = fail
!
Expand All @@ -175,6 +155,15 @@ end subroutine readrst
objgpopt(ipart,igp) = 1.00d+30
end do
!
if (nptocon > 0) violate(ipart) = 1.0d+30
!
end if
!
if (nptocon > 0) then
!
! since this is the initial condition, vioopt set directly to initial violate
vioopt(ipart) = violate(ipart)
!
end if
!
end do
Expand All @@ -183,7 +172,7 @@ end subroutine readrst
call getgbest(gindex,gbest,gmbest,gpbest)
!
objmin = gbest
!
!
! list output for initial conditions
call listout(0,gbest,gmbest,gpbest,gindex,objmin,iphistp)
!
Expand Down Expand Up @@ -244,8 +233,7 @@ end subroutine readrst
modfail(ipart) = 0
!
!---- get model run results
irun = ipart-1
call modelrm(0,irun,fail)
call modelrm(0,ipart-1,fail)
!
modfail(ipart) = fail
!
Expand All @@ -267,7 +255,19 @@ end subroutine readrst
end do
!
!-- set pbest
call getpbest()
if (nptocon > 0) then
!
do ipart=1,npop
call pbestcon(ipart)
end do
!
else
!
do ipart=1,npop
call getpbest(ipart)
end do
!
end if
!
!-- set gbest
call getgbest(gindex,gbest,gmbest,gpbest)
Expand Down Expand Up @@ -312,6 +312,9 @@ end subroutine readrst
!
objmin = gbest
!
! write pbest and gbest to file
call writebest(gindex,basnam)
!
end do
! end main loop of PSO iterative procedure
!-----------------------------------------
Expand All @@ -330,13 +333,6 @@ end subroutine readrst
!
! execute model run
err = rmif_run()
!
! get model run results
irun = 0
call modelrm(0,irun,fail)
!
! write pbest and gbest to file
call writebest(gindex,basnam)


end subroutine psoest
end subroutine psoest
Loading

0 comments on commit 50b66ed

Please sign in to comment.