Make live point generation deterministic

README.rst

@@ -5,7 +5,7 @@
  :target: https://arxiv.org/abs/1506.00171
  :alt: Open-access paper
 
-PolyChord v 1.20.2
+PolyChord v 1.21.0
 
 Will Handley, Mike Hobson & Anthony Lasenby
 

pypolychord/__init__.py

@@ -1,3 +1,3 @@
-__version__ = "1.20.2"
+__version__ = "1.21.0"
 from pypolychord.settings import PolyChordSettings
 from pypolychord.polychord import run_polychord

src/polychord/feedback.f90

@@ -28,7 +28,7 @@ subroutine write_opening_statement(settings)
  write(stdout_unit,'("")')
  write(stdout_unit,'("PolyChord: Next Generation Nested Sampling")')
  write(stdout_unit,'("copyright: Will Handley, Mike Hobson & Anthony Lasenby")')
- write(stdout_unit,'(" version: 1.20.2")')
+ write(stdout_unit,'(" version: 1.21.0")')
  write(stdout_unit,'(" release: 1st June 2021")')
  write(stdout_unit,'(" email: wh260@mrao.cam.ac.uk")')
  write(stdout_unit,'("")')

src/polychord/generate.F90

@@ -61,15 +61,15 @@ end function GenerateSeed
  subroutine GenerateLivePoints(loglikelihood,prior,settings,RTI,mpi_information)
  use settings_module, only: program_settings
  use random_module, only: random_reals
- use utils_module, only: write_phys_unit,DB_FMT,fmt_len,minpos,time
+ use utils_module, only: write_phys_unit,DB_FMT,fmt_len,minpos,time,sort_doubles
  use calculate_module, only: calculate_point
  use read_write_module, only: phys_live_file, prior_info_file
  use feedback_module, only: write_started_generating,write_finished_generating,write_generating_live_points
  use run_time_module, only: run_time_info,initialise_run_time_info, find_min_loglikelihoods
  use array_module, only: add_point
  use abort_module
 #ifdef MPI
- use mpi_module, only: mpi_bundle,is_root,linear_mode,throw_point,catch_point,more_points_needed,sum_integers,sum_doubles,request_point,no_more_points
+ use mpi_module, only: mpi_bundle,is_root,linear_mode,throw_point,catch_point,more_points_needed,sum_integers,sum_doubles,request_point,no_more_points,request_live_point,live_point_needed
 #else
  use mpi_module, only: mpi_bundle,is_root,linear_mode
 #endif
@@ -112,13 +112,15 @@ function prior(cube) result(theta)
 
  integer :: nlike ! number of likelihood calls
  integer :: nprior, ndiscarded
+ integer :: ngenerated ! use to track order points are generated in
 
  real(dp) :: time0,time1,total_time
  real(dp),dimension(size(settings%grade_dims)) :: speed
 
 
  ! Initialise number of likelihood calls to zero here
  nlike = 0
+ ngenerated = 1
 
 
  if(is_root(mpi_information)) then
@@ -190,6 +192,14 @@ function prior(cube) result(theta)
 
 
  active_workers=mpi_information%nprocs-1 ! Set the number of active processors to the number of workers
+ do worker_id=1,active_workers
+ ! Request a point from any worker
+ live_point(settings%h0:settings%h1) = random_reals(settings%nDims) ! Generate a random hypercube coordinate
+ ! use the time as an ordering identifier, cheat by using the birth contour
+ live_point(settings%b0) = ngenerated
+ ngenerated = ngenerated+1
+ call request_live_point(live_point,mpi_information,worker_id)
+ end do
 
  do while(active_workers>0) 
 
@@ -215,31 +225,36 @@ function prior(cube) result(theta)
 
 
  if(RTI%nlive(1)<nprior) then
- call request_point(mpi_information,worker_id) ! If we still need more points, send a signal to have another go
+ live_point(settings%h0:settings%h1) = random_reals(settings%nDims) ! Generate a random hypercube coordinate
+ ! use the time as a unique identifier, cheat by using the birth contour
+ live_point(settings%b0) = ngenerated
+ ngenerated = ngenerated+1
+ call request_live_point(live_point,mpi_information,worker_id)
  else
  call no_more_points(mpi_information,worker_id) ! Otherwise, send a signal to stop
  active_workers=active_workers-1 ! decrease the active worker counter
  end if
 
  end do
 
+ ! sort live points by order the prior samples were generated in
+ RTI%live(:,:RTI%nlive(1),:) = RTI%live(:,sort_doubles(RTI%live(settings%b0,:RTI%nlive(1),1)),:)
+ ! restore birth contour to logzero
+ RTI%live(settings%b0,:RTI%nlive(1),:) = settings%logzero
 
 
 
  else
 
  ! The workers simply generate and send points until they're told to stop by the administrator
- do while(.true.)
-
- live_point(settings%h0:settings%h1) = random_reals(settings%nDims) ! Generate a random hypercube coordinate
+
+ do while(live_point_needed(live_point,mpi_information))
  time0 = time()
  call calculate_point( loglikelihood, prior, live_point, settings,nlike) ! Compute physical coordinates, likelihoods and derived parameters
  ndiscarded=ndiscarded+1
  time1 = time()
- live_point(settings%b0) = settings%logzero
  if(live_point(settings%l0)>settings%logzero) total_time = total_time + time1-time0
  call throw_point(live_point,mpi_information) ! Send it to the root node
- if(.not. more_points_needed(mpi_information)) exit ! If we've recieved a kill signal, then exit this loop
 
  end do
  end if

src/polychord/mpi_utils.F90

@@ -693,6 +693,70 @@ function more_points_needed(mpi_information)
 
  end function more_points_needed
 
+ !> Request specific live points
+ ! administrator ----> worker
+ ! request_live_point(live_point) point_needed -> true
+ ! no_more_points (defined above) point_needed -> false
+ ! 
+
+ !> Request live point
+ !!
+ !! This subroutine is used by the root node to request a specific live point
+ subroutine request_live_point(live_point,mpi_information,worker_id)
+ implicit none
+ type(mpi_bundle), intent(in) :: mpi_information
+ integer, intent(in) :: worker_id !> Worker to request a new point from
+ real(dp), intent(in), dimension(:) :: live_point !> The live point to be sent
+
+
+ call MPI_SEND( &!
+ live_point, &! live point being sent
+ size(live_point), &!
+ MPI_DOUBLE_PRECISION, &! sending doubles
+ worker_id, &! process id to send to
+ tag_gen_request, &! continuation tag
+ mpi_information%communicator,&! mpi handle
+ mpierror &! error flag
+ )
+
+ end subroutine request_live_point
+
+ !> See if another point is needed
+ !!
+ !! This subroutine is used by the root node to request a specific live point
+ function live_point_needed(live_point,mpi_information)
+ use abort_module
+ implicit none
+ type(mpi_bundle), intent(in) :: mpi_information
+ real(dp),intent(out),dimension(:) :: live_point !> live point to throw
+
+ integer, dimension(MPI_STATUS_SIZE) :: mpistatus ! status identifier
+
+ logical :: live_point_needed !> Whether we need more points or not
+
+ call MPI_RECV( &!
+ live_point, &! live point recieved
+ size(live_point), &!
+ MPI_DOUBLE_PRECISION, &! receiving doubles
+ mpi_information%root, &! root node
+ MPI_ANY_TAG, &! mpi tag
+ mpi_information%communicator,&! mpi handle
+ mpistatus, &! status identifier
+ mpierror &! error flag
+ )
+
+ ! If we've recieved a kill signal, then exit this loop
+ if(mpistatus(MPI_TAG) == tag_gen_stop ) then
+ live_point_needed = .false.
+ else if(mpistatus(MPI_TAG) == tag_gen_request) then
+ live_point_needed = .true.
+ else
+ call halt_program('generate error: unrecognised tag')
+ end if
+
+ end function live_point_needed
+
+
 #endif