Using Latin hypercube sampler for Prior.draw

[DevarshiChoudhury]

Addresses feature request #520.

I have checked its expected application on the modelling notebook where we print the outputs of Prior.draw. The prior distributions looked ~identical, perhaps even slightly more uniformly distributed than before as seen by eye.

[thjsal]

Looks good to me. However, I find one thing worse with change. Looks like you cannot re-produce exactly the example results anymore every time you run them. I guess that is because you should fix somehow the seed of this new Latin hypercube sampler.

For example, if running TestRun_BB.py twice in a row, I get first time:

The support occupies an estimated 25.4% of the hypervolume within the unit hypercube...
Fractional hypervolume estimated.
Sampling efficiency set to: 3.1504.

And second time:

The support occupies an estimated 24.8% of the hypervolume within the unit hypercube...
Fractional hypervolume estimated.
Sampling efficiency set to: 3.2264.

Previously, fixing the numpy seed in the main script was enough to guarantee the numbers won't change here when executing several times.

[thjsal]

Perhaps you could fix the seed for this hypercube sampler in all the example files, so that they remain very re-producible. Adding something similar to this line:

xpsi/examples/examples_modeling_tutorial/TestRun_BB.py

Line 16 in d38db0b

np.random.seed(xpsi._rank+10)

[DevarshiChoudhury]

@thjsal it wasn't as straightforward a fix. I've added a LHS_seed arg in Prior.draw, Prior.estimate_hypercube_frac and Prior.unit_hypercube_frac. A LHS_seed is then added in runtime_params dict in the example files for Sample.nested to parse as a kwarg, passing it to NestedSampler as a kwarg which then uses it to call Prior.unit_hypercube_frac. I think this should make the examples fully reproducible although I haven't tested it myself. Also, if no LHS_seed is passed, it'll just default to a random seed.

[thjsal]

Seems like it is not working for me. I get first (when running TestRun_BB.py):

  File "/home/saltuomo/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-3.0.5-py3.12-linux-x86_64.egg/xpsi/NestedSampler.py", line 59, in __call__
    LHS_seed = LHS_seed if LHS_seed is not None else None
                           ^^^^^^^^
UnboundLocalError: cannot access local variable 'LHS_seed' where it is not associated with a value

That I can fix by doing LHS_seed = kwargs.setdefault('LHS_seed', None) instead of LHS_seed = LHS_seed if LHS_seed is not None else None

Then I get

  File "/home/saltuomo/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-3.0.5-py3.12-linux-x86_64.egg/xpsi/Prior.py", line 235, in draw
    sampler = qmc.LatinHypercube(d=len(self), 
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: LatinHypercube.__init__() got an unexpected keyword argument 'rng'

Seems like the "rng" input parameter was used to be named "seed" in Scipy version before 1.15. Do you think we should force X-PSI to use Scipy versions newer than this (it also requires Python newer than 3.10)? Anyway, I changed that to "seed", but now I get:

  File "/home/saltuomo/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-3.0.5-py3.12-linux-x86_64.egg/xpsi/NestedSampler.py", line 62, in __call__
    kwargs['sampling_efficiency'] /= self._prior.unit_hypercube_frac(LHS_seed=LHS_seed)
                                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: 'float' object is not callable

I don't know yet how to solve this one.

[DevarshiChoudhury]

Thanks for fixing the dict related issue. I had assumed that since the env files only mention scipy and not a version number, it would automatically install the latest version (v1.15.2) but I suppose not since I have v1.12.0 in my system. Do you envision any problem associated any other problem with updating to v1.15 for Scipy?

However, looking at your error, it also seems that the error is taking place at an earlier/later point, i.e. at self._prior.unit_hypercube_frac(LHS_seed=LHS_seed) and not qmc.LatinHypercube(d=len(self), different from where the rng was getting stuck, so I guess using the seed argument worked and now the issue is somewhere else?

It's not immediately obvious to me why the second error arises though

[thjsal]

I think it may use the latest Scipy version if creating the X-PSI environment now. But probably the version v1.15 was not yet available when we last time created that environment. I don't know if upgrading to that version would cause issues. Looks like at least the Github CI-tests are already using that version without issues. But would be good to add another CI-test for drawing samples, and calling the sampler (this is not tested currently). If we do that, we should also stop testing Python 3.8 and Python 3.9 as those would not work anymore. Alternatively, we could call the parameter "seed", since it should still work, probably just giving a deprecation warning for now.

Regarding the latest error, I found that is caused because giving input arguments for a property function is not allowed. The error goes away if not declaring the unit_hypercube_frac() as a property function in Prior.py.

Next error is from pymultinest which is not happy that it got an additional LHS_seed parameter.

  File "/home/saltuomo/.conda/envs/xpsi_py3/lib/python3.12/site-packages/xpsi-3.0.5-py3.12-linux-x86_64.egg/xpsi/NestedSampler.py", line 72, in __call__
    _ = pymultinest.solve(self._likelihood,
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/saltuomo/.conda/envs/xpsi_py3/lib/python3.12/site-packages/pymultinest/solve.py", line 71, in solve
    run(**kwargs)
TypeError: run() got an unexpected keyword argument 'LHS_seed'

[thjsal]

Okay, pushed now fixes to all the errors mentioned above. The last one (pymultinest error) was solved by just not having LHS_seed defined in the kwargs that are sent to pymultinest.

Now the test runs seem to work as intended. But perhaps should still check/think if the removal of the "@Property" definition from unit_hypercube_frac() function was completely fine, and whether we should force to upgrade to a newer Scipy.

I am also planning to set up (maybe in a couple of days) one or two new Github CI-tests to make sure the sampling works also with newest packages.

[thjsal]

I also changed LHS_seed=LHS_seed if LHS_seed is not None else None to simply LHS_seed=LHS_seed, since I could not understand why the longer version would be needed (it just sets the seed to None if it already is None).

[DevarshiChoudhury]

The @property allows the class method to be treated as a class attribute, where it can be directly called as has been done in NestedSampler as self._prior.unit_hypercube_frac. If I understand correctly, removing that would imply the need to call it as self._prior.unit_hypercube_frac() instead. Had you made that change to make it work?
I also checked & noticed that the unit_hypercube_frac has only been used in Prior & NestedSampler so if things are running without the @property, I imagine it should be fine.
We can also just force to upgrade to a newer Scipy. Perhaps something to bring up this Friday to discuss with others.

For the pymultinest error, if you didn't pass it via the runtime_params kwarg, how did you pass it?

[thjsal]

Had you made that change to make it work?
I did not change that, but you had changed it when trying to pass the LHS_seed argument to the self._prior.unit_hypercube_frac function in NestedSampler.py (you were calling it like a function). But I did not check other places where this function is used.

For the pymultinest error, if you didn't pass it via the runtime_params kwarg, how did you pass it?
I still passed it via the runtime_params but I defined the call function of NestedSampler to be def __call__(self, LHS_seed=None, **kwargs): instead of def __call__(self, **kwargs):. Then LHS_seed is not part of the kwargs and is not passed to pymultinest.

[DevarshiChoudhury]

Ah right, I was looking at the main branch where it was still called as a method.

In any case, let's hear the opinion of others on Python & SciPy updates, then make the appropriate changes before merging the pull request

[thjsal]

I made another pull request for the CI tests to be merged to this branch (#525 ). There I confirmed that indeed if using the "rng" argument instead of "seed", the tests only pass if using python3.10 or newer: https://github.com/xpsi-group/xpsi/actions/runs/13408191787

But they all pass if using "seed".

The exact value of the estimated prior volume using by either "seed" or "rng" is though different and needs to be changed depending on which one is used. Now I assume "seed" is used.