Duplicated information in netiter.combine_results

[Sorah-Darkhat]

Good morning Johannes,

Lately I have been a bit annoying, since I'm having a lot of issues with my model as it is quite expensive, so I want to apologize, and thank you for always replying fast!

Currently, I have decided to downgrade my model's precision to at least get a posterior and understand how to better constrain my model to get results in reasonable times. For this reason I decreased the number of live points to 400, yet this produces a small number of effective samples, so I increased the min_ess to 10.000.

I assume this is the correct way to go about this, since this will then increase the number of live points near the posterior region, rather than sampling from the prior.

My issue is that I'm running my model in a cluster using 128 nodes each with 128 cores for a total of 16384 cores.
Each node has 256 Gb or ram, and this is plenty for the sampling, but during the save of the results the code is always killed by what I assumed was the OOM, what surprised me was that it worked just fine if I decreased the number of cores, so I went into the code and realized that the issue was in a mpi_comm.gather, I then removed the gather and used the simple send and recv, this worked but then it crashed during the bcast, so I also changed the bcast into a send and recv loop of 512MB chunks to obey the MPI memory limit, once again this worked, but then it still crashed.

So, to prevent every core of having the full saved_logwt_bs and logZ_bs, I gathered them in rank 0, and tried to bcast the results dictionary from rank 0 as I assumed this would be lighter.

However, during what I assume is the rank order test, it iteratively generates the results dictionary and stores it into memory, which is mostly fine if the results dictionary is small, the issue is that when I use 1 core, the results dictionary is relatively light, as well as saved_logwt_bs and logZ_bs, but when I run with 16384 cores, these files are in the order of the GB, which then result in massive dictionaries, and in the code crashing.

I assume the information must then be duplicated across the cores, because I'm resuming an already concluded run, and with 1 to 128 cores, everything is fine, but when I start to go into the thousand cores things get funky on the ram usage and crash.

Once again, sorry for troubling you, but I was wondering if this is something you could help me look into since I'm not as familiar in what the code is actually doing.

With the best regards,

João

[JohannesBuchner]

Hi João,

I'm glad you are looking at the code and bringing it to these very large compute jobs. It seems to me that you are getting a good understanding of what the code is doing.

I can see that the memory could be an issue. The object which does the nested sampling integration on a graph (allowing for dynamic number of live points), MultiCounter, is built here: https://github.com/JohannesBuchner/UltraNest/blob/master/ultranest/integrator.py#L2604 nbootstraps is not allowed to be smaller than 1. You don't actually need that many bootstraps if you are running with so many cores. Maybe we should allow some cores to not do bootstraps at all.

In combine_results, the cores share the results: https://github.com/JohannesBuchner/UltraNest/blob/master/ultranest/netiter.py#L880
I can see that the data then might get large (number of iterations x number of cores x ...).

I did not expect that that might be an issue in computing clusters, I thought they'd have loads of memory.

For sure, some or most calculations could also only be done at rank 0, or the first num_bootstraps ranks, ideally.
In fact it is a problem for MLFriends here: https://github.com/JohannesBuchner/UltraNest/blob/master/ultranest/integrator.py#L388 and https://github.com/JohannesBuchner/UltraNest/blob/master/ultranest/integrator.py#L2141
A bootstrap is recomputed at least once per core, and then the worst case across all calculations is taken. This means the more cores, the less efficient the MLFriends region proposal. In these three places, the max(1, nbootstraps // self.mpi_size) should probably be replaced by something like max(1, nbootstraps // self.mpi_size) if self.mpi_rank < nbootstraps else 0, and the called function needs to be checked if it can handle zero number of bootstraps and return something sensible for that case (zero radius).

[Sorah-Darkhat]

Hm...!

Thank you for the reply! I have also just noticed that the try/except block that leads to raising np.linalg.LinAlgError in the region.compute_enlargement, depends on the bootstraps, I'm now wondering if the reason my model was previously failing to previously build the mlfriends region is also a side effect of the number of cores I'm currently using, I will make a short run using fewer cores to see if this is actually the reason.

I will see if I can make all of this work in my cluster, which does have a rather small ram for the number of cores in each node.

[JohannesBuchner]

Yes, that sounds likely, good catch.

[Sorah-Darkhat]

The number of cores didn’t seem to fix the region exception.

I also managed to work around the RAM issue by getting a bit creative. Instead of launching 16 384 MPI tasks, I launch 128, one per node, and then vectorize my likelihood call by distributing the calculations across the cores within each node. This seems to work fine.

For the regular ReactiveNestedSampler, this is probably a simple fix for HPCs with low ram.

The only catch is with the SliceSampler, while it parallelizes fine when evaluating 16 384 calls at once, the vectorization only assigns one call per MPI task, so I am limited to evaluating 128 calls at a time.

[JohannesBuchner]

That's a good solution. You can adjust ndraw_min and ndraw_max as well.

You could use the ultranest.popstepsampler.PopulationSimpleSliceSampler or PopulationSliceSampler to vectorize the slice sampler on that level, but SliceSampler is probably better currently. The PopulationSimpleSliceSampler is not expected to converge, but give quick look results. However, recently there is an idea to make it unbiased, it just needs a small change described at #170 and - of course - testing.

There is another issue with ingesting PopulationSliceSampler results in the right order that PR #167 tries to fix when running with a high number of MPI cores, but users have reported the fix has a problem. I could not figure out why though, it does not make sense to me, the code should work. Maybe you can have a look at the patch as well, it might help you too.

[Sorah-Darkhat]

Hello Johannes,
Sorry for reopening this issue, but is SliceSampler vectorization supposed to return vectors of 1? I assumed it was, but from your reply, I'm now wondering if it is supposed to return more.

I have ndraw_min set to 512, yet it only draws 1, is it possible to make it draw the 512?

Thank you once again!

[EDIT]
When sampling the live points from the prior, it does samples more than 1, but during the actual run it returns vectors of 1.

[JohannesBuchner]

__next__ (the iterator) is supposed to return tuples of u, p, L, nc . The first three contain either None or information for one point.
https://github.com/JohannesBuchner/UltraNest/blob/master/ultranest/popstepsampler.py#L358

Internally it advances multiple, but in each call the nested sampling integrator expects to ingest (up to) one point per node. Only 1 node needs to succeed for the nested sampling iteration to advance, so the slice sampling population / efficiency should be chosen accordingly.

In slice samplers, __next__ only returns a point after nsteps have been completed, which can take a long time. Probably about popsize/(nsteps*ncallperstep) points are done at the bottom of a __next__ call.

[Sorah-Darkhat]

My current issue, is that I changed the code even further, and now I'm testing some very heavy likelihoods, but I'm stopping the calculation mid way for most samples:

If your loglikelihood can only decrease with data (all terms are negative, for example by normalizing to the highest possible density, or because it is -0.5*chi²), then you can stop computing early, as you suggest. Your likelihood could access sampler.Lmin, which gives the current likelihood threshold. Making the likelihood a class, for which you can set the sampler after you create the sampler, might be helpful. Likelihood points below sampler.Lmin are rejected by the sampler anyhow -- their value is inconsequential.

The issue is that, while some stop very early, taking only 2 or 3 seconds to compute the likelihood, my new likelihood takes around 300 to 900 seconds, depending on the specific parameters. This would be fine, because I'm okay getting a very rough posterior so I can use a small number of live points as well as a small number of nsteps, ( currently I'm not too worried about the quality ). But the issue is that even if 90% of the samples are only taking 2/3 seconds to complete, the cores end up waiting for the samples that take 900 seconds, and this is slowlying the inference down tremendously.

An easy fix for me, would be to just make the cores do vectors of samples, so that they can keep computing instead of waiting for the samples that are computing the whole likelihood. This should be in principle possible, because ultranest already allows me to use a very large number of cores, my idea was to use the same cores, but to vectorize them:

I also managed to work around the RAM issue by getting a bit creative. Instead of launching 16 384 MPI tasks, I launch 128, one per node, and then vectorize my likelihood call by distributing the calculations across the cores within each node. This seems to work fine.

So that they keep computing until they find a point that is above Lmin, and then they send that point or those points for the chains to progress. I might be misunderstanding how the slice sampler works, because it seems possible given the fact that ultranest already accepts a lot of cores.

My alternative would be to oversubscribe the number of cores, and just use an ungodly number of mpi tasks, but thats bound to lead to ram issues, or other problems of the sort.

This would not be an issue on the reactive nested sampler, because there the vectorization is simple, and if I submit 1000 samples, some are bound to take 900 seconds but most will take 1 or 2, and it evens out, but sadly it is much slower than the slice sampler.

I was wondering if it would be possible to unblock the cores, and let them keep computing instead of waiting for the ones taking 900 seconds? Vectorization would fix this, but the popstepsampler seems also much slower than the slice sampler, so I was trying to get a work around for the slice sampler.

[JohannesBuchner]

What you suggest sounds very complicated, would have a lot of race conditions, and probably will not give you a speed-up factor comparable to 300s/3s.
You could modify the SliceSampler to perform your quick check, and if rejected already goes ahead with the next proposal. There is something along those lines for the "region_filter" functionality.