Using Latin hypercube sampler for Prior.draw (#521)

DevarshiChoudhury · thjsal · web-flow · commit 2af6ad12a166 · 2025-05-12T13:02:45.000+03:00
* Using Latin hypercube sampler instead for `Prior.draw`

* Testing sampling that sampling does not crash in the continous integration tests

---------

Co-authored-by: Tuomo Salmi &lt;t.h.j.salmi@uva.nl&gt;
Co-authored-by: Tuomo Salmi &lt;tuomo.salmi@helsinki.fi&gt;
diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml
@@ -35,7 +35,7 @@ jobs:
           python=${{ matrix.python-version }}
     - name: mamba installs
       run: |
-        micromamba install flake8 pytest
+        micromamba install flake8 pytest h5py "ultranest>=4.0"
     - name: Set up GCC
       run: |
          sudo apt-get install -y gcc
@@ -49,6 +49,23 @@ jobs:
         flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    #- name: multinest and pymultinest install (seems working but the pytest below still cannot find the MultiNest library files)
+    #  run: |
+    #   sudo apt-get install cmake libblas-dev liblapack-dev libatlas-base-dev
+    #   sudo apt install mpich
+    #   git clone https://github.com/farhanferoz/MultiNest.git multinest
+    #   cd multinest/MultiNest_v3.12_CMake/multinest/
+    #   mkdir build
+    #   cd build
+    #   CC=gcc FC=mpif90 CXX=g++ cmake -DCMAKE_{C,CXX}_FLAGS="-O3 -march=native -funroll-loops" -DCMAKE_Fortran_FLAGS="-O3 -march=native -funroll-loops" ..
+    #   make
+    #   cd ../../../../
+    #   git clone https://github.com/JohannesBuchner/PyMultiNest.git pymultinest
+    #   cd pymultinest
+    #   python setup.py install
+    #   cd ..
+    #   export LD_LIBRARY_PATH=$(pwd)/multinest/MultiNest_v3.12_CMake/multinest/lib/:$LD_LIBRARY_PATH
+    #   echo $LD_LIBRARY_PATH
     - name: install X-PSI package
       run: |
         python setup.py install
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -33,6 +33,43 @@ and this project adheres to
 .. Attribution
 .. ^^^^^^^^^^^
 
+
+[v3.0.6] - 2025-05-12
+~~~~~~~~~~~~~~~~~~~~~
+
+Summary
+^^^^^^^
+
+* Drawing samples from prior to now use Latin hypercube sampler (LHS) to more uniformly draw samples over high-dimensions.
+
+Added
+^^^^^
+
+* Continuous integration test for sampling with UltraNest added.
+
+Changed
+^^^^^^^
+
+* `Prior.draw` modified to use `scipy.stats.qmc.LatinHypercube` instead of `numpy.random.rand` for drawing samples from prior distribution.
+* `Prior.estimate_hypercube_frac` modified to accept the kwarg `LHS_seed` to enable reproducible LHS.
+* `Prior.unit_hypercube_frac` modified to accept the kwarg `LHS_seed` to enable reproducible LHS.
+* `NestedSampler.__call__` accepts `LHS_seed` as a kwarg to passs to `Prior.unit_hypercube_frac`
+* `examples/examples_fast/Sampling.ipynb`: `runtime_params` includes a fixed `LHS_seed`.
+* `examples.examples_fast.Modules.main_IS_likelihood.runtime_params` includes a fixed `LHS_seed`.
+* `examples.examples_fast.Modules.main_IS_prior.runtime_params` includes a fixed `LHS_seed`.
+* `examples.examples_fast.Modules.main.runtime_params` includes a fixed `LHS_seed`.
+* `examples.examples_modeling_tutorial.TestRun_BB.runtime_params` includes a fixed `LHS_seed`.
+* `examples.examples_modeling_tutorial.TestRun_Num.runtime_params` includes a fixed `LHS_seed`.
+* `examples.examples_modeling_tutorial.TestRun_NumBeam.runtime_params` includes a fixed `LHS_seed`.
+* `examples.examples_modeling_tutorial.TestRun_PolNum_split_inference.runtime_params` includes a fixed `LHS_seed`.
+
+Attribution
+^^^^^^^^^^^
+
+Devarshi Choudhury
+Tuomo Salmi
+
+
 [v3.0.5] - 2025-03-19
 ~~~~~~~~~~~~~~~~~~~~~
 
@@ -42,7 +79,7 @@ Summary
 * Issues with phase bin interpolations were fixed when using the time-invariant version of the Everywhere class.
 
 Fixed
-^^^^^^^
+^^^^^
 
 * Phase bin interpolations are now skipped if having only one phase bin during the likelihood calculation (``xpsi/likelihoods/default_background_marginalisation.pyx``) or during the data synthesisation (``xpsi/tools/synthesise.pyx``). Previously, these resulted in errors.
 
@@ -106,6 +143,7 @@ Changed
 
 Attribution
 ^^^^^^^^^^^
+
 Devarshi Choudhury,
 Tuomo Salmi,
 Bas Dorsman
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -80,7 +80,7 @@
 # The short X.Y version.
 version = '3.0-b'
 # The full version, including alpha/beta/rc tags.
-release = '3.0.5'
+release = '3.0.6'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/examples/examples_fast/Modules/main.py b/examples/examples_fast/Modules/main.py
@@ -155,6 +155,7 @@
                       'evidence_tolerance': 0.1,
                       'max_iter': 100, #-1,
                       'seed' : 0, # Fixing the seed
+                      'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation
                       'verbose': True}
 
 
diff --git a/examples/examples_fast/Modules/main_IS_likelihood.py b/examples/examples_fast/Modules/main_IS_likelihood.py
@@ -155,6 +155,7 @@
                       'evidence_tolerance': 0.1,
                       'max_iter': 100, #-1,
                       'seed' : 0, # Fixing the seed
+                      'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation
                       'verbose': True}
 
 
diff --git a/examples/examples_fast/Modules/main_IS_prior.py b/examples/examples_fast/Modules/main_IS_prior.py
@@ -152,6 +152,7 @@
                       'evidence_tolerance': 0.1,
                       'max_iter': 100, #-1,
                       'seed' : 0, # Fixing the seed
+                      'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation
                       'verbose': True}
 
 
diff --git a/examples/examples_fast/Sampling.ipynb b/examples/examples_fast/Sampling.ipynb
@@ -781,7 +781,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": null,
    "metadata": {
     "scrolled": true
    },
@@ -804,7 +804,8 @@
     "#                       'wrapped_params': wrapped_params,\n",
     "#                       'evidence_tolerance': 0.1,\n",
     "#                       'max_iter': -1,\n",
-    "#                       'seed' : 0, # Fixing the seed \n",
+    "#                       'seed' : 0, # Fixing the seed\n",
+    "#                       'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation', \n",
     "#                       'verbose': True}\n",
     "\n",
     "# xpsi.Sample.nested(likelihood, prior, **runtime_params)\n",
diff --git a/examples/examples_modeling_tutorial/TestRun_BB.py b/examples/examples_modeling_tutorial/TestRun_BB.py
@@ -435,6 +435,7 @@ def transform(self, p, **kwargs):
                   'evidence_tolerance': 0.5,
                   'seed': 7,
                   'max_iter': 100,# manual termination condition for short test
+                  'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation
                   'verbose': True}
 
 xpsi.set_phase_interpolant('Akima')
diff --git a/examples/examples_modeling_tutorial/TestRun_Num.py b/examples/examples_modeling_tutorial/TestRun_Num.py
@@ -616,6 +616,7 @@ def transform(self, p, **kwargs):
                   'wrapped_params': wrapped_params,
                   'evidence_tolerance': 0.5,
                   'seed': 7,
+                  'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation
                   'max_iter': 100, # manual termination condition for short test
                   'verbose': True}
 
diff --git a/examples/examples_modeling_tutorial/TestRun_NumBeam.py b/examples/examples_modeling_tutorial/TestRun_NumBeam.py
@@ -586,6 +586,7 @@ def transform(self, p, **kwargs):
                   'wrapped_params': wrapped_params,
                   'evidence_tolerance': 0.5,
                   'seed': 7,
+                  'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation
                   'max_iter': 100, # manual termination condition for short test
                   'verbose': True}
 
diff --git a/examples/examples_modeling_tutorial/TestRun_PolNum_split_inference.py b/examples/examples_modeling_tutorial/TestRun_PolNum_split_inference.py
@@ -724,6 +724,7 @@ def transform(self, p, **kwargs):
                   'wrapped_params': wrapped_params,
                   'evidence_tolerance': 0.5,
                   'seed': 7,
+                  'LHS_seed': 42, # Fixing the LHS seed for hypercube fraction estimation
                   'max_iter': 100, # manual termination condition for short test
                   'verbose': True}
 
diff --git a/setup.py b/setup.py
@@ -271,7 +271,7 @@ def EXTENSION(modname):
 
     setup(
         name = 'xpsi',
-        version = '3.0.5',
+        version = '3.0.6',
         author = 'The X-PSI Core Team',
         author_email = 'A.L.Watts@uva.nl',
         url = 'https://github.com/xpsi-group/xpsi',
diff --git a/xpsi/NestedSampler.py b/xpsi/NestedSampler.py
@@ -48,7 +48,7 @@ def __init__(self,
             self._prior = prior
 
     @make_verbose('Commencing integration', 'Integration completed')
-    def __call__(self, **kwargs):
+    def __call__(self, LHS_seed=None, **kwargs):
         """ Integrate.
 
         :param kwargs:
@@ -57,7 +57,7 @@ def __call__(self, **kwargs):
         """
 
         kwargs.setdefault('sampling_efficiency', 0.8)
-        kwargs['sampling_efficiency'] /= self._prior.unit_hypercube_frac
+        kwargs['sampling_efficiency'] /= self._prior.unit_hypercube_frac(LHS_seed=LHS_seed)
 
         yield 'Sampling efficiency set to: %.4f.'%kwargs['sampling_efficiency']
 
diff --git a/xpsi/Prior.py b/xpsi/Prior.py
@@ -9,6 +9,8 @@
 from abc import ABCMeta, abstractmethod
 from xpsi.ParameterSubspace import ParameterSubspace
 
+from scipy.stats import qmc
+
 class Prior(ParameterSubspace, metaclass=ABCMeta):
     """ The joint prior distribution of parameters (including hyperparameters).
 
@@ -220,7 +222,7 @@ def index(self, name):
         return len(self) + self.derived_names.index(name)
 
     @make_verbose('Drawing samples from the joint prior','Samples drawn')
-    def draw(self, ndraws, transform=False):
+    def draw(self, ndraws, transform=False, LHS_seed=None):
         """ Draw samples uniformly from the prior via inverse sampling.
 
         :param int ndraws: Number of draws.
@@ -229,7 +231,12 @@ def draw(self, ndraws, transform=False):
 
         """
 
-        h = _np.random.rand(int(ndraws), len(self))
+        # Create a Latin Hypercube sampler
+        sampler = qmc.LatinHypercube(d=len(self), 
+                                     seed=LHS_seed)
+
+        # Generate Latin Hypercube samples in the range [0, 1]
+        h = sampler.random(n=ndraws)
 
         if transform:
             try:
@@ -263,7 +270,7 @@ def draw(self, ndraws, transform=False):
                 else:
                     redraw = ndraws
                 redraw -= finite_counter
-                h = _np.random.rand(int(redraw)+1, len(self))
+                h = sampler.random(n=int(redraw) + 1)
                 index = 0
                 if transform:
                     try:
@@ -284,7 +291,7 @@ def draw(self, ndraws, transform=False):
     @make_verbose('Estimating fractional hypervolume of the unit hypercube '
                   'with finite prior density:',
                   'Fractional hypervolume estimated')
-    def estimate_hypercube_frac(self, ndraws=5):
+    def estimate_hypercube_frac(self, ndraws=5, LHS_seed=None):
         """
         Estimate using Monte Carlo integration the fractional hypervolume
         within a unit hypercube at which prior density is finite.
@@ -299,29 +306,28 @@ def estimate_hypercube_frac(self, ndraws=5):
             yield ('Requiring %.E draws from the prior support '
                    'for Monte Carlo estimation' % ndraws)
 
-            self._unit_hypercube_frac = self.draw(ndraws)[1]
+            self._unit_hypercube_frac = self.draw(ndraws, 
+                                                  LHS_seed=LHS_seed)[1]
         else:
             self._unit_hypercube_frac = None
 
         if _size > 1:
             self._unit_hypercube_frac = _comm.bcast(self._unit_hypercube_frac,
                                                     root=0)
-
         yield ('The support occupies an estimated '
                '%.1f%% of the hypervolume within the unit hypercube'
                % (self._unit_hypercube_frac*100.0))
 
         yield self._unit_hypercube_frac
 
-    @property
-    def unit_hypercube_frac(self):
+    def unit_hypercube_frac(self, LHS_seed=None):
         """ Get the fractional hypervolume with finite prior density. """
-
         try:
             return self._unit_hypercube_frac
         except AttributeError:
             try:
-                self.estimate_hypercube_frac(self.__draws_from_support__)
+                self.estimate_hypercube_frac(self.__draws_from_support__, 
+                                             LHS_seed=LHS_seed)
             except AttributeError:
                 print('Cannot locate method for estimating fraction.')
             else:
diff --git a/xpsi/__init__.py b/xpsi/__init__.py
@@ -4,7 +4,7 @@
 
 """
 
-__version__ = "3.0.5"
+__version__ = "3.0.6"
 __author__ = "The X-PSI Core Team"
 
 try:
diff --git a/xpsi/tests/test_sampling.py b/xpsi/tests/test_sampling.py
