Synchronous parallelisation (#67)

williamjameshandley · Lukas Hergt · web-flow · commit 53eb09457589 · 2022-02-03T14:39:27.000Z
* Synchronous parallelisation implemented * Updated version numbers * add missing i for int for the additional boolean sync parameter in weird str describing parameter types (#66) * Updated branch * Added anesthetic to the README * Modification so that now all dead points are saved for anesthetic usage * Bug fix for linear mode setting birth contour * changed huge to logzero in logweights, as this causes inability to resume with intel compilers Co-authored-by: Lukas Hergt <lh561@mrao.cam.ac.uk>
diff --git a/README.rst b/README.rst
@@ -5,7 +5,7 @@
    :target: https://arxiv.org/abs/1506.00171
    :alt: Open-access paper
 
-PolyChord v 1.18.2
+PolyChord v 1.20.0
 
 Will Handley, Mike Hobson & Anthony Lasenby
 
@@ -56,6 +56,14 @@ If any of the above steps fail (this can in general happen for certain Mac OSX v
 
 If you do not have sudo access/virtual environments/anaconda, then appending `--user` to the install command may be necessary.
 
+Post Processing
+===============
+
+We recommend the tool `anesthetic <https://github.com/williamjameshandley/anesthetic>`_ for post-processing your nested sampling runs. A plot gallery can be found `here <http://htmlpreview.github.io/?https://github.com/williamjameshandley/cosmo_example/blob/master/demos/demo.html>`_
+
+
+https://github.com/williamjameshandley/anesthetic
+
 MPI Support
 ===========
 
diff --git a/pypolychord/__init__.py b/pypolychord/__init__.py
@@ -199,6 +199,7 @@ def wrap_prior(cube, theta):
                      settings.write_prior,
                      settings.maximise,
                      settings.compression_factor,
+                     settings.synchronous,
                      settings.base_dir,
                      settings.file_root,
                      settings.grade_frac,
diff --git a/pypolychord/_pypolychord.cpp b/pypolychord/_pypolychord.cpp
@@ -126,7 +126,7 @@ static PyObject *run_pypolychord(PyObject *, PyObject *args)
         
 
     if (!PyArg_ParseTuple(args,
-                "OOOiiiiiiiiddidiiiiiiiiiiidssO!O!O!i:run",
+                "OOOiiiiiiiiddidiiiiiiiiiiidissO!O!O!i:run",
                 &temp_logl,
                 &temp_prior,
                 &temp_dumper,
@@ -154,6 +154,7 @@ static PyObject *run_pypolychord(PyObject *, PyObject *args)
                 &S.write_prior,
                 &S.maximise,
                 &S.compression_factor,
+                &S.synchronous,
                 &base_dir,
                 &file_root,
                 &PyList_Type,
diff --git a/pypolychord/settings.py b/pypolychord/settings.py
@@ -120,6 +120,14 @@ class PolyChordSettings:
         (Default: exp(-1))
         How often to update the files and do clustering
 
+    synchronous : boolean
+        (Default: True)
+        Parallelise with synchronous workers, rather than asynchronous ones.
+        This can be set to False if the likelihood speed is known to be
+        approximately constant across the parameter space. Synchronous
+        parallelisation is less effective than asynchronous by a factor ~O(1)
+        for large parallelisation.
+
     base_dir : string
         (Default: 'chains')
         Where to store output files.
@@ -179,6 +187,7 @@ def __init__(self, nDims, nDerived, **kwargs):
         self.maximise = kwargs.pop('maximise', False)
         self.compression_factor = kwargs.pop('compression_factor',
                                              numpy.exp(-1))
+        self.synchronous = kwargs.pop('synchronous', True)
         self.base_dir = kwargs.pop('base_dir', 'chains')
         self.file_root = kwargs.pop('file_root', 'test')
         self.seed = kwargs.pop('seed', -1)
diff --git a/src/drivers/polychord_CC.cpp b/src/drivers/polychord_CC.cpp
@@ -31,6 +31,7 @@ int main()
 
     settings.feedback      = 1;
     settings.compression_factor  = 0.36787944117144233;
+    settings.synchronous   = true;
 
     settings.boost_posterior= 5.0;
 
diff --git a/src/polychord/c_interface.cpp b/src/polychord/c_interface.cpp
@@ -28,6 +28,7 @@ Settings::Settings(int _nDims,int _nDerived):
     write_prior         {true},
     maximise            {true},
     compression_factor  {0.36787944117144233},
+    synchronous         {true},
     base_dir            {"chains"},
     file_root           {"test"},
     grade_frac          {1.0},
@@ -95,6 +96,7 @@ void run_polychord(
             s.write_prior,
             s.maximise,
             s.compression_factor,
+            s.synchronous,
             s.nDims,
             s.nDerived,
             base_dir,
diff --git a/src/polychord/feedback.f90 b/src/polychord/feedback.f90
@@ -28,8 +28,8 @@ 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.18.2")')
-            write(stdout_unit,'("  release: 7th April 2020")')
+            write(stdout_unit,'("  version: 1.20.0")')
+            write(stdout_unit,'("  release: 1st June 2021")')
             write(stdout_unit,'("    email: wh260@mrao.cam.ac.uk")')
             write(stdout_unit,'("")')
         end if
@@ -40,6 +40,7 @@ subroutine write_opening_statement(settings)
             write(stdout_unit,'("nDims    :",I8)')   settings%nDims
             write(stdout_unit,'("nDerived :",I8)')   settings%nDerived
             if(settings%do_clustering) write(stdout_unit,'("Doing Clustering")')
+            if(settings%synchronous) write(stdout_unit,'("Synchronous parallelisation")')
             if(settings%equals) write(stdout_unit,'("Generating equally weighted posteriors")')
             if(settings%posteriors) write(stdout_unit,'("Generating weighted posteriors")')
             if((settings%equals.or.settings%posteriors).and.settings%cluster_posteriors.and.settings%do_clustering) write(stdout_unit,'("Clustering on posteriors")')
diff --git a/src/polychord/generate.F90 b/src/polychord/generate.F90
@@ -196,8 +196,8 @@ function prior(cube) result(theta)
                     ! Recieve a point from any worker
                     worker_id = catch_point(live_point,mpi_information)
 
-                    ! If its valid, and we need more points, add it to the array
-                    if(live_point(settings%l0)>settings%logzero .and. RTI%nlive(1)<nprior) then
+                    ! If its valid, add it to the array
+                    if(live_point(settings%l0)>settings%logzero) then
 
                         call add_point(live_point,RTI%live,RTI%nlive,1) ! Add this point to the array
 
diff --git a/src/polychord/ini.f90 b/src/polychord/ini.f90
@@ -75,6 +75,7 @@ subroutine read_params(file_name,settings,params,derived_params)
         settings%write_prior        = get_logical(file_name,'write_prior',.false.)
         settings%maximise           = get_logical(file_name,'maximise',.false.)
         settings%compression_factor = get_double(file_name,'compression_factor',exp(-1d0))
+        settings%synchronous        = get_logical(file_name,'synchronous',.true.)
         settings%base_dir           = get_string(file_name,'base_dir','chains')
         settings%file_root          = get_string(file_name,'file_root','test')
         settings%seed               = get_integer(file_name,'seed', -1)
diff --git a/src/polychord/interfaces.F90 b/src/polychord/interfaces.F90
@@ -282,6 +282,7 @@ subroutine polychord_c_interface(&
             write_prior,&
             maximise,&
             compression_factor,&
+            synchronous,&
             nDims,&
             nDerived, &
             base_dir,&
@@ -360,6 +361,7 @@ end subroutine c_dumper
         logical(c_bool), intent(in), value  :: write_prior
         logical(c_bool), intent(in), value  :: maximise
         real(c_double),  intent(in), value  :: compression_factor
+        logical(c_bool), intent(in), value  :: synchronous
         integer(c_int), intent(in), value   :: nDims
         integer(c_int), intent(in), value   :: nDerived
         character(len=1,kind=c_char), intent(in), dimension(STR_LENGTH) :: base_dir
@@ -407,6 +409,7 @@ end subroutine c_dumper
         settings%write_prior         = write_prior
         settings%maximise            = maximise
         settings%compression_factor  = compression_factor         
+        settings%synchronous         = synchronous         
         settings%nDims               = nDims
         settings%nDerived            = nDerived
 
diff --git a/src/polychord/interfaces.h b/src/polychord/interfaces.h
@@ -25,6 +25,7 @@ extern "C" void polychord_c_interface(
         bool,
         bool,
         double,
+        bool,
         int,
         int,
         char*,
diff --git a/src/polychord/interfaces.hpp b/src/polychord/interfaces.hpp
@@ -31,6 +31,7 @@ struct Settings
     bool write_prior;
     bool maximise;
     double compression_factor;
+    bool synchronous;
     std::string base_dir;
     std::string file_root;
     std::vector<double> grade_frac;
diff --git a/src/polychord/nested_sampling.F90 b/src/polychord/nested_sampling.F90
@@ -128,6 +128,12 @@ end subroutine dumper
         ! has reorganised the cluster indices
         integer ::  worker_epoch
         integer ::  administrator_epoch
+
+        ! Nursary for storing babies in synchronous parallel mode
+        real(dp), allocatable, dimension(:,:,:) :: nursary
+        integer :: i_nursary
+        integer, allocatable, dimension(:) :: worker_epochs
+        integer, allocatable, dimension(:,:) :: nlikes
 #endif
 
 
@@ -143,6 +149,7 @@ end subroutine dumper
         ! worker switch
         worker_epoch=0
         administrator_epoch=0
+        i_nursary=0
 #endif
 
         ! Rolling loglikelihood calculation
@@ -217,6 +224,8 @@ end subroutine dumper
         end if
 #ifdef MPI
         call broadcast_integers(num_repeats,mpi_information)
+        allocate(nursary(settings%nTotal,sum(num_repeats), mpi_information%nprocs-1))
+        allocate(worker_epochs(mpi_information%nprocs-1), nlikes(size(settings%grade_dims),mpi_information%nprocs-1))
 #endif
         allocate(baby_points(settings%nTotal,sum(num_repeats)))
 
@@ -252,7 +261,34 @@ end subroutine dumper
 
                     ! Generate a new set of points within the likelihood bound of the late point
                     baby_points = SliceSampling(loglikelihood,prior,settings,logL,seed_point,cholesky,nlike,num_repeats)
+                    baby_points(settings%b0,:) = logL ! Note the moment it is born at
 #ifdef MPI
+                else if(settings%synchronous) then
+                    ! Parallel synchronous mode
+                    ! -------------------------
+
+                    if(i_nursary == 0) then
+                        do worker_id=1,mpi_information%nprocs-1
+                            seed_point = GenerateSeed(settings,RTI,cluster_id)
+                            cholesky = RTI%cholesky(:,:,cluster_id)
+                            logL = RTI%logLp(cluster_id)
+                            call throw_seed(seed_point,cholesky,logL,mpi_information,worker_id,administrator_epoch,.true.)
+                            worker_cluster(worker_id) = cluster_id
+                        end do
+                        do i_worker=1,mpi_information%nprocs-1
+                            i_nursary = catch_babies(baby_points,nlike,worker_epoch,mpi_information)
+                            nursary(:,:,i_nursary) = baby_points
+                            worker_epochs(i_nursary) = worker_epoch
+                            nlikes(:,i_nursary) = nlike
+                        end do
+                        i_nursary = mpi_information%nprocs-1
+                    end if
+                    baby_points = nursary(:,:,i_nursary)
+                    cluster_id = worker_cluster(i_nursary)
+                    worker_epoch = worker_epochs(i_nursary)
+                    nlike = nlikes(:,i_nursary)
+                    i_nursary = i_nursary-1
+
                 else
                     ! Parallel mode
                     ! -------------
@@ -391,18 +427,23 @@ end subroutine dumper
             ! Kill off the final workers.
             ! If we're done, then clean up by receiving the last piece of
             ! data from each node (and throw it away) and then send a kill signal back to it
-            do i_worker=mpi_information%nprocs-1,1,-1
-
-                ! Recieve baby point from worker worker_id
-                worker_id = catch_babies(baby_points,nlike,worker_epoch,mpi_information)
+            if (settings%synchronous) then
+                do worker_id=mpi_information%nprocs-1,1,-1
+                    call throw_seed(seed_point,cholesky,logL,mpi_information,worker_id,administrator_epoch,.false.) 
+                end do
+            else
+                do i_worker=mpi_information%nprocs-1,1,-1
 
-                ! Add the likelihood calls to our counter
-                RTI%nlike = RTI%nlike + nlike
+                    ! Recieve baby point from worker worker_id
+                    worker_id = catch_babies(baby_points,nlike,worker_epoch,mpi_information)
 
-                ! Send kill signal to worker worker_id (note that we no longer care about seed_point, so we'll just use the last one
-                call throw_seed(seed_point,cholesky,logL,mpi_information,worker_id,administrator_epoch,.false.) 
+                    ! Add the likelihood calls to our counter
+                    RTI%nlike = RTI%nlike + nlike
 
-            end do
+                    ! Send kill signal to worker worker_id (note that we no longer care about seed_point, so we'll just use the last one
+                    call throw_seed(seed_point,cholesky,logL,mpi_information,worker_id,administrator_epoch,.false.) 
+                end do
+            end if
 
 
         else !(myrank/=root)
@@ -421,10 +462,13 @@ end subroutine dumper
             ! On the first loop, send a nonsense set of baby_points
             ! to indicate that we're ready to start receiving
 
-            baby_points = 0d0                              ! Avoid sending nonsense
-            baby_points(settings%l0,:) = settings%logzero  ! zero contour to ensure these are all thrown away
-            nlike = 0                                      ! no likelihood calls in this round
-            call throw_babies(baby_points,nlike,worker_epoch,mpi_information)
+            if (.not. settings%synchronous) then
+                baby_points = 0d0                              ! Avoid sending nonsense
+                baby_points(settings%l0,:) = settings%logzero  ! zero contour to ensure these are all thrown away
+                baby_points(settings%b0,:) = settings%logzero  ! zero birth contour
+                nlike = 0                                      ! no likelihood calls in this round
+                call throw_babies(baby_points,nlike,worker_epoch,mpi_information)
+            end if
             wait_time = 0
             slice_time = 0
             time1 = time()
@@ -438,6 +482,7 @@ end subroutine dumper
                 time0 = time()
                 ! 2) Generate a new set of baby points
                 baby_points = SliceSampling(loglikelihood,prior,settings,logL,seed_point,cholesky,nlike,num_repeats)
+                baby_points(settings%b0,:) = logL ! Note the moment it is born at
 
 
                 wait_time = wait_time + time0-time1
diff --git a/src/polychord/run_time_info.f90 b/src/polychord/run_time_info.f90
@@ -717,7 +717,7 @@ function replace_point(settings,RTI,baby_points,cluster_add)
         use utils_module, only: logsumexp,logincexp
         use settings_module, only: program_settings
         use random_module, only: bernoulli_trial
-        use array_module, only: add_point
+        use array_module, only: add_point, reallocate
 
         implicit none
         type(program_settings), intent(in) :: settings !> Program settings
@@ -774,11 +774,14 @@ function replace_point(settings,RTI,baby_points,cluster_add)
                     replace_point = .true.                                                           ! Mark this as a replaced live point
                 end if
                 if (sum(RTI%nlive) < nlive) then
-                    point(settings%b0) = logL                                                        ! Note the moment it is born at
                     call add_point(point,RTI%live,RTI%nlive,cluster_add)                             ! Add the new live point to the live points
                     call find_min_loglikelihoods(settings,RTI)                                       ! Find the new minimum likelihoods
                 end if
             end if
+        else
+            call add_point(point,RTI%dead,RTI%ndead)
+            if (RTI%ndead > size(RTI%logweights)) call reallocate(RTI%logweights,RTI%ndead*2)
+            RTI%logweights(RTI%ndead) = settings%logzero
         end if
 
     end function replace_point
diff --git a/src/polychord/settings.f90 b/src/polychord/settings.f90
@@ -143,6 +143,8 @@ module settings_module
         
         real(dp) :: compression_factor = exp(-1d0)
 
+        logical :: synchronous = .true.
+
     end type program_settings
 
     contains
