Parallelization issues with pocoMC

[lisaotten]

When employing the pocoMC package for Bayesian Inference runs using tellurium for modeling, we have encountered issues with parallelization. Using multiprocess(ing), we noticed a very big discrepancy between the results obtained by parallelized and non-parallelized runs (see also attached corner plots). Both runs run through smoothly without any error messages or other large differences. I have been able to reproduce the results of both runs separately multiple times both on a HPC and my personal notebook. Besides the change of parallelization, all other parameters are kept exactly the same. Changes in the number of parallel kernels do not seem to change these results. The results of the non-parallelized run are what we would expect as the correct results from experience.

I have attached a self-contained code including the environment I am running it on. The config-file has an option under bayesian_inference to turn parallelization on/off as well as specify the number of kernels.

parallel.pdf
not_parallel.pdf
Bayesian_Transporter.zip

[matthiaskoenig]

Hi @lisaotten,

I had a quick look at the plot. Not sure the plots are really different.
You have to set the axes of the two plots identical to do a better comparison visually.
Often you get single rare samples which are far outside resulting in very different axes ranges (if the axes are adapted automatically). I assume you have just 1-2 outlier (very unlikely samples) in the one run resulting in very small distributions visually due to changes of the axes limits.

You should do some actual test if the distributions are different, e.g. by comparing the modes of your multi-dimensional distributions or using something like:
EFECT – A Method and Metric to Assess the Reproducibility of Stochastic Simulation Studies
T.J. Sego, Matthias König, Luis L. Fonseca, Baylor Fain, Adam C. Knapp, Krishna Tiwari, Henning Hermjakob, Herbert M. Sauro, James A. Glazier, Reinhard C. Laubenbacher, Rahuman S. Malik-Sheriff
arXiv:2406.16820 (preprint). doi:10.48550/arXiv.2406.16820

Hope this helps.
TDLR: most likely a plotting issue, not a sampling issue

[lisaotten]

Hi @matthiaskoenig,

Thanks for your reply!
Both plots actually have the exact same axis ranges. I agree with you that the parallel results are much broader than the results from the non-parallel run, which is where my problem lies. I have reproduced these results multiple times with very similar results both on a HPC and my personal notebook.

[luciansmith]

The first thing I can think of is that if the seeds are being set from the system clock (which they are by default), you might be getting the same seed on multiple runs? I can imagine reasons for both the parallel and non-parallel runs to end up this way, so you might want to try explicitly setting the seed for each run manually, to ensure that they're all unique.

You could also examine the individual results to see if this is actually happening or not.

[lisaotten]

I played around with the seed quite a bit in finding a possible source for the errors, but even setting the seed manually resulted in the same distributions.

[lisaotten]

I have attached a combined corner plot that maybe illustrates the issues a little better: The black line corresponds to a non-parallel run, while the other three lines correspond to parallel runs using the multiprocess, multiprocessing and pathos packages. They were all created using the same fixed seed at the start of the runs.

The last parameter in the corner plot corresponds to the deviation of the parameter fits from the data points that we are trying to analyze and is much larger for the parallel run. This clearly indicates that the parallel runs fit the data much worse.

[hsauro]

I wonder if it is possible to create a small problem to look at? I can see there appears to be a difference. What happens if you run it again, do you get a similar corner plot?

…

On Wed, Sep 18, 2024 at 9:39 AM lisaotten ***@***.***> wrote: I have attached a combined corner plot that maybe illustrates the issues a little better: The black line corresponds to a non-parallel run, while the other three lines correspond to parallel runs using the multiprocess, multiprocessing and pathos packages. They were all created using the same fixed seed at the start of the runs. The last parameter in the corner plot corresponds to the deviation of the parameter fits from the data points that we are trying to analyze and is much larger for the parallel run. This clearly indicates that the parallel runs fit the data much worse. corner_16Dpococheck_3.png (view on web) <https://urldefense.com/v3/__https://github.com/user-attachments/assets/3656420b-7719-4161-b5ff-cbecb31a2313__;!!K-Hz7m0Vt54!gttSNTitQRpvNxt7NVNEJMih4al3YU4A6XOICHbJTWgYKRSHGcXs3IuufcYhNZu4FC7hF4Kiu2goWAJu6zfQoC1Vq0uYCA$> — Reply to this email directly, view it on GitHub <https://urldefense.com/v3/__https://github.com/sys-bio/tellurium/issues/591*issuecomment-2358945406__;Iw!!K-Hz7m0Vt54!gttSNTitQRpvNxt7NVNEJMih4al3YU4A6XOICHbJTWgYKRSHGcXs3IuufcYhNZu4FC7hF4Kiu2goWAJu6zfQoC10I_7x7Q$>, or unsubscribe <https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AAIBSDRLGAZ5HYL34HTUPTDZXGUDNAVCNFSM6AAAAABMA2544KVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNJYHE2DKNBQGY__;!!K-Hz7m0Vt54!gttSNTitQRpvNxt7NVNEJMih4al3YU4A6XOICHbJTWgYKRSHGcXs3IuufcYhNZu4FC7hF4Kiu2goWAJu6zfQoC3RPX1GBw$> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- Herbert Sauro, Professor Director: NIH Center for model reproducibility University of Washington, Bioengineering 206-685-2119, www.sys-bio.org, http://reproduciblebiomodels.org/ Mobile: 206-880-8093 ***@***.*** Books: http://books.analogmachine.org/

[hsauro]

I’ve been thinking more about this. Would it be possible to create a simple one variable model, eg -> x ->. With simple mass action on each step? That would create a one element corner plot. This would be easier to debug.

[luciansmith]

@lisaotten : We have some time again to look at this--can you explain what it is that you are varying and what you are collecting across runs? I'd like to try to pare down your example to something small that can be more easily tested.

[lisaotten]

Hi @luciansmith,

Thank you for following up on this. I am sorry for not diving into this more, I haven't been working on the project in the last few month as other projects have taken priority.

In the example that I send, we are inferring 16 parameters: 11 rate constants for the transporter and 5 additional nuisance parameters (concentrations and noise). In the config file you can find the names, ranges and the ground truth of these parameters:

- name: "log10_k1_f"
  bounds: [6,12]
  nominal: 10
- name: "log10_k1_r"
  bounds: [-1, 5]
  nominal: 3
- name: "log10_k2_f"
  bounds: [-2, 4]
  nominal: 2
- name: "log10_k2_r"
  bounds: [-2, 4]
  nominal: 2
- name: "log10_k3_f"
  bounds: [3,9]
  nominal: 7
- name: "log10_k3_r"
  bounds: [-1, 5]
  nominal: 3
- name: "log10_k4_f"
  bounds: [-1, 5]
  nominal: 3
- name: "log10_k4_r"
  bounds: [6,12]
  nominal: 10
- name: "log10_k5_f"
  bounds: [-2,4]
  nominal: 2
- name: "log10_k5_r"
  bounds: [-2,4]
  nominal: 2
- name: "log10_k6_f"
  bounds: [-1,5]
  nominal: 3
- name: "H_out_scale"
  bounds: [0.8,1.2]
  nominal: 1
- name: "S_out_scale"
  bounds: [0.8,1.2]
  nominal: 1
- name: "protein_conc_scale"
  bounds: [0.8,1.2]
  nominal: 1
- name: "bias_scale"
  bounds: [0.8,1.2]
  nominal: 1
- name: "log10_sigma"
  bounds: [-11, -9]
  nominal: -10.5

It might be feasible to set some of these parameters to their true value to reduce the dimensionality of the problem.
I am also attaching the publication associated with the code, in case this is helpful.

transporter-manuscript-final-30jan.pdf

Again, thank you for looking into this. Please let me know, if I can be of any assistance.

[luciansmith]

OK, looking through your code, here's what I think you can do to ensure that roadrunner is acting independently in each of your runs. Basically, you'll create a roadrunner object, then save its state, then pass the state to the simulate_assay function instead of the object.

Every time you load the SBML model, you then convert it into a 'state':

    rr_model = te.loadSBMLModel(model_file)
    rr_state = rr_model.saveStateS()

Then replace all your rr_model arguments with rr_state, and replace your ssme_function with the attached file:

ssme_function.py.txt

Basically, the new function takes a rr_state as an input instead of a model, then runs:

    rr = roadrunner.RoadRunner()
    rr.loadStateS(rr_state)

to set up a new roadrunner instance.

[lisaotten]

Thank you! I will test this out.