Restartable models #111

HDembinski · 2022-12-31T13:57:19Z

An attempt to close #106

This branch is based on #112, which should be merged first.

This is using two approaches. For most models, it is sufficient run some initialization code exactly once, which I ensure with the new _once private method. Concrete models should not implemented __init__ anymore, only _once which is called by MCRun.__init__ exactly once and after the underlying implementation has been initialized.

For Phojet and DPMJet, this does not work (at least I could not get it to work), so I followed the idea expressed in #106. I wrote new base class MCRunRemote, which a brilliant and terrible piece of engineering. Behind the scenes, it runs every call to cross_section and __call__ through a remote process which is created and dies just for the call. The beauty of it: a model which requires this workaround just has to inherit from MCRunRemote instead of MCRun, no other change required!

To maintain the illusion of having a single model and not one that is restarted all the time, I rely strongly on the random-state serialization feature that @jncots added. I further improved this so that the PRNG state of the numpy generator is also saved. Since this so important, I added more tests that check the RNG state serialization and found some issues. @jncots could you please look into this?

Other changes

MCRun.get_stable returns set of particles that were changed to stable
Kinematics is now passed to MCRun.__call__ instead of MCRun.__init__; this makes more sense for models that are restartable. Some models need a maximum energy upfront, this can be passed via init.
MCRun.kinematics property was removed. The kinematics should not be part of the state of the generator, we pass them explicitly to the two functions which need them, MCRun.__call__ and MCRun.cross_section. Less class state is generally better, it simplifies the design and the reasoning about something. It does not matter that the underlying Fortran stores the kinematics, because we abstracted that away in our high-level API.
MCRun.random_state now includes the state of the numpy PRNG
MCRun._composite_plan now generates heavy elements first, as a workaround for DPMJet
MCRun.set_unstable was renamed to MCRun.maydecay

HDembinski · 2023-01-02T23:30:57Z

@afedynitch You said that you are looking into fixing cross-sections for the models. I added a new test test_cross_sections.py which checks the output of all models for various combinations of projectiles and targets. The references for such checks are very rough, of course, because the models seem to disagree a lot especially on nuclear cross-sections, however, the test still checks basic sanity of the model response and I found several bugs this way, too.

HDembinski · 2023-01-02T23:44:55Z

From my side, this PR is ready for review, I need your help to fix the remaining things.

Tests are failing largely because of the new tests for rng_state persistence that I added (that is unrelated to my changes).

HDembinski · 2023-01-02T23:58:00Z

At least one of the tests is stalling, on my own computer and on CI. I haven't figured out which one yet. This means you have to use a keyboard interrupt to complete the tests at the moment. On Windows the tests abort with a memory error, these things could be related, although I did not notice a large increase in RAM on my computer. Parallel computation works differently on Windows and Unix, so it could well be that there is a mistake which Windows reveals more drastically.

Edit: Stalling is now fixed, the issue is well understood.

HDembinski · 2023-01-03T09:14:50Z

The test which runs forever is test_to_hepmc3.py, I am investigating.

src/impy/models/dpmjetIII.py

afedynitch · 2023-01-04T09:40:12Z

src/impy/models/dpmjetIII.py

-            return CrossSectionData(total=stot, elastic=sela, inelastic=stot - sela)
+
+        assert kin.p1.is_hadron
+        assert kin.p2.is_hadron


allowed targets are _projectiles, but in reality it's only Nuclei(), proton and neutron.

asserts are just an executable comment, I put this here to remind myself that we only deal with hadrons at this point in the function. The nuclear targets are handled by the part above.

afedynitch · 2023-01-04T09:42:14Z

src/impy/models/dpmjetIII.py


        # TODO set more cross-sections

    def _set_kinematics(self, kin):
        # Save maximal mass that has been initialized
        # (DPMJET sometimes crashes if higher mass requested than initialized)
+
+        # Relax momentum and energy conservation checks at very high energies
+        if kin.ecm > 5e4:


Can't do that without setting it back. Like this the output will depend on the sequence of requested energies since energy-momentum conservation checks can reject events. Either piece has to stay in init or, one needs to save the default values and reset them if kinematics change.

It cannot be in init.

Ah, I see what you mean, but it does not really matter, I think. We need to use MCRunRemote which always uses a fresh copy of dpmjet when _set_kinematics is called, so it is always called only once per dpmjet instance.

I tried to fix this a bit, but you are the DPMJet expert, so I am relying on your help.

afedynitch · 2023-01-04T09:45:32Z

src/impy/models/dpmjetIII.py

@@ -160,9 +158,9 @@ def _set_kinematics(self, kin):
            self._lib.dt_init(
                -1,
                max(kin.plab, 100.0),


Also the maximal initialization energy cannot be exceeded. So, it has to be saved as well.

dt_init can be only called once, right? At least that was my impression from observing the code. MCRunRemote makes it so that _set_kinematics is called exactly once per instance.

afedynitch · 2023-01-04T09:48:36Z

src/impy/models/epos.py

-
-        super().__init__(seed)
+    def _once(self, ecm_max=1000 * TeV):
+        from impy import debug_level


Does this imply that the debug level cannot change in runtime? It is important if you search for something that goes wrong at event 12345 but don't want to see the output of 12344 events.

Another reason to keep it changeable in runtime is that different debug levels in the models can produce extremely different levels of output.

I don't think that the functionality was there before consistently. So might be no problem, when using this "remote run".

N-th edit: If I change the debug level in impy.debug_level in a notebook or python script, and do a second initialization, it won't be reflected here. Something like this would be great.

impy.debug_level = 0 m1 = Epos() # do something impy.debug_level = 5 m2 = Phojet() # do something different

Another reason to keep it changeable in runtime is that different debug levels in the models can produce extremely different levels of output.

I agree that being able to change debug levels for a model instance is good, but the API we have so far (also before my changes) does not allow this.

If you pass the debug level to fortran in a function that is called only at init, then you cannot change the debug level after init. However, we should be able to access the variable that is used in the fortran code to enable debug output, and then you can change it at any time.

Regarding this:

impy.debug_level = 0 m1 = Epos() # do something impy.debug_level = 5 m2 = Phojet() # do something different

I thought this works, this is why debug_level is imported freshly from top level impy (and to avoid some circular import issue). Did you test it and found it does not work? It does not work, if you call Epos again instead of Phojet, but it should work if it is two different models, and both are created for the first time.

afedynitch · 2023-01-04T09:53:42Z

src/impy/models/phojet.py

+        if self._lib.pho_event(-1, *self._beams)[1]:
+            raise RuntimeError(
+                "initialization failed with the current event kinematics"
+            )


Nope, this is wrong. pho_event(-1, ....) has to be called once for the highest energy requested for a run. Different combinations of particles and energies are set with pho_setpar. Once an event generation is triggered (pho_event(>0, ..)), PHOJET runs the remaining initializations. So pho_event(-1 only once.

Can we call it once with a maximum energy that is never exceeded? like 1000 * TeV?

afedynitch · 2023-01-04T09:59:54Z

src/impy/models/pythia8.py

-
-    def _cross_section(self, kin=None):
+    def _cross_section(self, kin):
+        if (kin.p2.A or 1) > 1:


It's really that bad at nuclei :)

I don't know, perhaps you have to query the cross-section from another variable. At some point I dereference a nullptr in the calls when nuclei are used. I didn't check what exactly happens.

afedynitch · 2023-01-04T10:06:37Z

src/impy/models/sibyll.py

-    def __init__(self, evt_kin, *, seed=None):
-        super().__init__(seed)
+    def _once(self):
+        from impy import debug_level


Same here, no change in run time is an unnecessary design limitation.

I agree, but I did not design the library in this way, it was already like that.

src/impy/models/sibyll.py

afedynitch · 2023-01-04T10:10:10Z

src/impy/models/sibyll.py

    _version = "2.1"
    _library_name = "_sib21"


-class Sibyll23(SIBYLLRun):
+# For some reason, Sibyll23 requires MCRunRemote, but the others don't


Does this mean that only Sibyll23 crashes on multiple inits? In SIBYLL, the initialization should not be called multiple times even if it works. If I remove the check in PHOJET it will also work but will result in undefined or altered behavior.

Yes, you can try it out. Replace MCRunRemote with MCRun further below as base class for Sibyll23 and run the tests. I don't understand why, so it is good if someone else checks this.

…n without remote

HDembinski added 5 commits January 2, 2023 16:25

wip

6cf4162

fix

a85ba34

cleanup

26c769e

bug-fix for epos

9255d40

cleanup

f904713

HDembinski force-pushed the singletons branch from 5d71dc3 to 92773c1 Compare January 2, 2023 16:26

HDembinski marked this pull request as ready for review January 2, 2023 17:20

HDembinski changed the title ~~Turn models into singletons~~ Restartable models Jan 2, 2023

update urqmd reference

8a4d1ff

HDembinski requested review from jncots and afedynitch and removed request for jncots January 2, 2023 23:35

afedynitch reviewed Jan 4, 2023

View reviewed changes

src/impy/models/dpmjetIII.py Outdated Show resolved Hide resolved

afedynitch reviewed Jan 4, 2023

View reviewed changes

src/impy/models/dpmjetIII.py Outdated Show resolved Hide resolved

afedynitch reviewed Jan 4, 2023

View reviewed changes

src/impy/models/sibyll.py Outdated Show resolved Hide resolved

afedynitch reviewed Jan 4, 2023

View reviewed changes

HDembinski added 2 commits January 4, 2023 13:49

fix wrong mass for proton,neutron

af161ce

update urqmd reference

0213fad

HDembinski added 24 commits January 4, 2023 15:38

work in progress

c752368

ported sibyll

1f4982c

ported pythia6

173589e

ported pythia8

baf41b4

move CompositeTarget to util, use process_particle everywhere

b9b6921

port sophia

d41a497

port urqmd

6372c93

ported phojet

eb88290

fixes

ac6fb7d

fix tests, add MCRun.get_stable, add numpy state to rng state

c9adb18

revert name change to test references

40d062e

remove obsolete inits

0f02a8f

IT WORKS

fe45cfe

massive fixes and new tests for cross sections

632da2d

terminate process if it is still alive

e8b6963

kill

34fd277

fix bug in remote control, add reprs

b0555d5

update notebook

698e31c

mark bad combinations for cross_section test, allow MCRunRemote to ru…

2029384

…n without remote

fix remotecall

fd69c9b

fix

74cc715

fix

fda3574

fix

4f662f1

check more stuff

367d90c

HDembinski force-pushed the singletons branch from 5526ca1 to 367d90c Compare January 4, 2023 14:39

HDembinski and others added 4 commits January 4, 2023 17:21

make naming in sibyll more consistent

9b401fe

fix my bug introduced to dpmjet

12ede78

examples

7a2da75

Merge branch 'main' into singletons

b3ab91c

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Restartable models #111

Restartable models #111

HDembinski commented Dec 31, 2022 •

edited

Loading

HDembinski commented Jan 2, 2023

HDembinski commented Jan 2, 2023

HDembinski commented Jan 2, 2023 •

edited

Loading

HDembinski commented Jan 3, 2023

afedynitch Jan 4, 2023

HDembinski Jan 4, 2023

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023

HDembinski Jan 4, 2023

afedynitch Jan 4, 2023

HDembinski Jan 4, 2023 •

edited

Loading

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

afedynitch Jan 4, 2023

HDembinski Jan 4, 2023

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

Restartable models #111

Are you sure you want to change the base?

Restartable models #111

Conversation

HDembinski commented Dec 31, 2022 • edited Loading

Other changes

HDembinski commented Jan 2, 2023

HDembinski commented Jan 2, 2023

HDembinski commented Jan 2, 2023 • edited Loading

HDembinski commented Jan 3, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

afedynitch Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

HDembinski Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

afedynitch Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

HDembinski Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

HDembinski Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

afedynitch Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

HDembinski Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

afedynitch Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

afedynitch Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

HDembinski Jan 4, 2023 • edited Loading

Choose a reason for hiding this comment

HDembinski commented Dec 31, 2022 •

edited

Loading

HDembinski commented Jan 2, 2023 •

edited

Loading

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading

afedynitch Jan 4, 2023 •

edited

Loading

afedynitch Jan 4, 2023 •

edited

Loading

HDembinski Jan 4, 2023 •

edited

Loading