Separate photon sampling from SOPHIA

[MHoerbe]

Dear all,

this PR extracts the internal SOPHIA photon sampling methods and re-implements them in CRPropa.

This is a necessary preparatory step towards the implementation of variable photon fields as discussed in #220 (sometimes referred to as "step zero" in the agreed road map). The changes are:

• SOPHIA's sampling methods have been implemented in C++
• PhotoPionProduction now calls the C++ versions of SOPHIA's photon sampling methods
• This allowed for SOPHIA's signature to be changed to more meaningful parameters
• A direct SOPHIA interface has been added which could be used outside of a simulation

This PR comes with scripts which

• run minimal simulations of the PhotoPionProduction on a wider range of parameters
• analyse the similarity of the introduced changes vs. the master branch

To make a review easier:

• 5c650ec provides the same functionalities as all later commits
• 77d0f60 deletes other parts of SOPHIA which are not used under any circumstances in CRPropa (Note that no non-linear behaviour of the code was touched with this)
• The tests above compare simulation data (i.e. at the highest hierarchy level). However, there are also tests available for each single new method which is introduced with this PR if this would become important at some point.

Important remarks:

• Please consider that the CRPropa naming style has not been met at each point to make comparisons to the SOPHIA methods feasible. Also the unit convention is broken internally in these methods. I would most eagerly recommend to amend the naming issues once the changes themselves have been reviewed.
• The scripts for testing are found in testScripts.zip. Note that the simulation script needs to be changed for each branch in which you are running the script (one for master, one for these changes, description provided).

Please do get in touch at any given point!

[TobiasWinchen]

Nice work - there are still some readability issues that I noted. Also, please try to make the review not harder than necessary by changing lines without reason. This is in particular true for the old fortran code, where you changed several thousands of lines, and a lot are e.g. changes of capitalizations of print to PRINT. This does not help in understanding what happened.

[MHoerbe]

Thanks lot @TobiasWinchen for the fast reply! I have implemented your feedback on all your points up to my best understanding.

Regarding the changes in SOPHIA: Sorry, I should have elaborated more on what actually changed.

• @adundovi and me concluded that, since we are modifying SOPHIA anyways, we should aim for the maximal amount of code reduction we are able to achieve given the change of the CRPropa/SOPHIA interface.
• the steps of modifying SOPHIA therefore evolve as follows:

1. Minimal solution: out-comment calls to internal sampling routines
2. Middle solution: delete all photon sampling routines from SOPHIA (as they are not used)
3. Maximal solution: delete all unused code from SOPHIA (checked for GOTO violations, comprehensive tests provided) and improve readability by proper indentations and styles

What you were seeing here was the maximal solution and I totally see how this is confusing at first glance. For the sake of a more efficient review, I will change sophia_interface.f to the minimal solution to ease this step of iteration while hoping to get to the maximal solution later on.

[TobiasWinchen]

Thanks for addressing my comments. There is IMO one major issue left with the gauss integration I missed before.

[adundovi]

* @adundovi and me concluded that, since we are modifying SOPHIA anyways, we should aim for the maximal amount of code reduction we are able to achieve given the change of the CRPropa/SOPHIA interface.

@TobiasWinchen I confirm this. The idea is to clear the SOPHIA code as much as possible from unneeded parts, so that we can focus on what really matters and proceed with fixing and improving PPP. The current implementation of PPP with SOPHIA is not only slow (cannot be parallelized beyond 8+ threads) and unmaintainable (totally unreadable), but also lacking the sampling on other photon fields which some users consider as a bug (because it was not explicitly stated). Thus, Mario's work here is more like a bug fix, not a new feature implementation. Therefore, keeping the old state of the code in parallel in this case is not reasonable. After this commit, or how Mario calls it, the zeroth step, we can proceed with new contributions which can be checked systematically against the old one.

[adundovi]

Tobias covered all needed changes that I can think of already. My only suggestion is that we start to use more strict types in interfaces (such as the enum class) across the code to make the code more readable, and to try catching as many bugs as we can during the compile time, not run time.

[TobiasWinchen]

On Tuesday 7 January 2020 10:53:08 CET Andrej Dundovic wrote: > * @adundovi and me concluded that, since we are modifying SOPHIA > anyways, we should aim for the maximal amount of code reduction we > are able to achieve given the change of the CRPropa/SOPHIA interface. @TobiasWinchen I confirm this. The idea is to clear the SOPHIA code as much as possible from unneeded parts, so that we can focus on what really matters and proceed with fixing and improving PPP. The current implementation of PPP with SOPHIA is not only slow (cannot be parallelized beyond 8+ threads) and unmaintainable (totally unreadable), but also lacking the sampling on other photon fields which some users consider as a bug (because it was not explicitly stated). Thus, Mario's work here is more like a bug fix, not a new feature implementation. Therefore, keeping the old state of the code in parallel in this case is not reasonable. After this commit, or how Mario calls it, the zeroth step, we can proceed with new contributions which can be checked systematically against the old one.

I agree that we should remove unused code, but changing e.g. capitalizations from print to PRINT at the same time just blows up the commit and makes the review harder.

[adundovi]

I agree that we should remove unused code, but changing e.g. capitalizations from print to PRINT at the same time just blows up the commit and makes the review harder.

Oh, sorry, I, of course, agree with that.

[MHoerbe]

Thanks all again for your useful comments. I am happy to resolve the integration part and consider enums for sophia particles. May I also emphasise that, based on your comments, the current version of the suggested SOPHIA code refrains from any cosmetics and only out-comments code which is not used as it is related to internal photon sampling (i.e. to be deleted later or now). The interface itself in sophiaevent(...) , of course, contains a very small amount of newly written code.

[adundovi]

@Froehliche-Kernschmelze I have just tried to write you more code snippets for the special SophiaParticle type and realized that with the current state of the implementation where sophiaEvent() and performInteraction() methods are separated while doing the same PID conversion, it makes little sense to do what I suggested. So please, skip the suggestion for the moment. We'll do it later when we merge sophiaEvent() into performInteraction().

[adundovi]

I'm fine with this PR now. Some parts will obviously be reworked once sampling on other field models are introduced, but as an intermediate step, these changes are fine.

[MHoerbe]

Thanks again to you two. Before we merge,

1. shall I summarise here what the KS tests yield on the respective observables or do you prefer running the scripts from above on your systems as well first?
2. Should I remove the out-commented SOPHIA parts or leave them as "comments"?
3. If at all, would you like to have the rest of "unused SOPHIA" removed in a later step or another PR? This would reduce SOPHIA to roughly 9k lines. The same test scripts could show here that nothing changed, too

[adundovi]

1. shall I summarise here what the KS tests yield on the respective observables or do you prefer running the [scripts from above](https://github.com/CRPropa/CRPropa3/files/3983762/testScripts.zip) on your systems as well first?

Are those the same tests you already presented in the last CRPropa call? If yes, this is OK, we don't need to go through these checks again.

2. Should I remove the out-commented SOPHIA parts or leave them as "comments"?

Yes, please. Just don't change "print" to "PRINT" or similar things that are complicating the review process.

3. If at all, would you like to have the rest of "unused SOPHIA" removed in a later step or another PR? This would reduce SOPHIA to roughly 9k lines. The same test scripts could show here that nothing changed, too

Is this the old sampling part? If yes, remove them too as we agreed. If this is something else, keep it for the next PR.

[MHoerbe]

1. shall I summarise here what the KS tests yield on the respective observables or do you prefer running the [scripts from above](https://github.com/CRPropa/CRPropa3/files/3983762/testScripts.zip) on your systems as well first?

Are those the same tests you already presented in the last CRPropa call? If yes, this is OK, we don't need to go through these checks again.

The test scripts above are identical to the tests I had presented, however, on a wider selection of parameters. The proposed tests conduct tests on all combinations of:

• CMB, IRB_Kneiske04 (i.e. all photon fields)
• primary protons and neutrons (i.e. all input particle types)
• primary particle energies of log10(E / GeV) = {7, 8, 9, 10, 11, 12, 13, 14} (i.e. a selection of some energies from above the threshold up to quite high energies)

and compare the results of this PR to the master branch via KS tests.

2. Should I remove the out-commented SOPHIA parts or leave them as "comments"?

Yes, please. Just don't change "print" to "PRINT" or similar things that are complicating the review process.

3. If at all, would you like to have the rest of "unused SOPHIA" removed in a later step or another PR? This would reduce SOPHIA to roughly 9k lines. The same test scripts could show here that nothing changed, too

Is this the old sampling part? If yes, remove them too as we agreed. If this is something else, keep it for the next PR.

All right, I will delete the sampling-related parts in this case. Everything else will be subject to another PR.

[adundovi]

@TobiasWinchen shell we merge it now?