custom photon fields / SOPHIA update #220
Conversation
…ed from reference to pointer due to problems with swig.
…onal version but it works jet.
… is an analytical solution for small distances from the symmetry axis (MagneticBottle) . Also added a magnetic field in which charged particles perform a gyration motion in a small and a big scale (GyroField). The right comment to these two classes is pending.
…er version of a candidate
|
Now that step 0 is taken, there is the agreement to introduce a second version of SOPHIA that includes my modifications in parallel to its original version. This version is agreed to be toggled by some boolean switch in the PhotoPionProduction, defaulting to original SOPHIA. Yet I am having trouble to actually link this second version of SOPHIA to CRPropa and after going through all sorts of attempts I am asking for a few thoughts of yours. Objective: State of progress How SOPHIA is implemented / linked (main CMakeLists) (sophia.h) (sophia_interface.f) (SOPHIA CMakeLists) (import / use in PPP) in all these files I have implemented additional entries of the same syntax yet mentioning an object called SOPHIAMOD instead, including a new directory called sophiaMod in the /libs folder. I added an import in the PPP Problem Do you have any educated guesses what the solution could be? Otherwise I would be happy to provide a more detailed description of this problem. There is a branch in my repository that currently has these settings. Also, the PPP has got a version of sophiaeventMod (toggled out) that if enabled would produce the error described. |
|
Give me one day - I'll write a CMake option for it from which we would be able to switch between sophia and sophiamod. |
|
Sorry for the delay. I somehow repeated the procedure you described and it worked for me. Or to be more specific:
Here you can see the diff in my branch. I can import and use both modules in python in parallel after doing this. |
|
Thanks a real lot to @adundovi ! A branch of mine now contains a modified version of SOPHIA following the suggestions of @TobiasWinchen where only unused parts of the code have been removed. No indentation / other style specifications / goto removal or new functions and the like have been realized. This is to tackle step 1 of the agreed roadmap. It should be shown that nothing changed to SOPHIA by removing unused parts of the code. This modified version can be used via a new temporary module called At this point, how shall we proceed? I could create a pull request with this modification or conduct further tests on the scientific identity of the two versions. So far I have conducted Kolmogorov-Smirnov tests on each of SOPHIAs output particles between two samples 1) drawn from the original SOPHIA in the master branch (for comparison of the KS parameters) and 2) drawn from master and SOPHIAmod. Results suggest identity between all of the produced secondary output particles, i.e. photons, electrons and neutrinos (hadrons not yet tested). I am happy to share any setup or tests I have made. Any opinions? |
|
Tests are certainly welcome if not necessary. I suggest that you write them in the form of unit tests and place them in It would be great if we could control the SOPHIA random number generator, i.e. to fix its seed so that both modules (the old and the new) produce deterministic results for given inputs. In that way, we would avoid time and resource consuming unit tests which are based on statistical comparisons. Moreover, if you have some theoretical estimates which should be easy to reproduce with SOPHIA, they could be included, too. Both modules should pass all these tests. After ending the basic testing step, we should proceed by including your PPPmod in the master repo following manually testing in different "large-scale scenarios". If everything goes well we should make PPPmod as the default PPP module, while the current one label as Legacy and after some time to put a deprecated notice and remove it from the source. |
|
There is progress with two aspects:
Note that you should get exactly the same output as SOPHIAs RNG always starts with the same seed. I observe that if you renamed sophiaevent_ to sophiaeventMod_ in libs/sophiaMod/sophiaMod.h, libs/sophiaMod/sophiaMod_interface.f and include/crpropa/module/PhotoPionProductionMod.cpp, there will be an error upon compilation stating "libcrpropa.so: undefined reference to sophiaeventMod_". Do you have any thoughts on this? |
That might be problematic, as one usually assumes that the events are independent of each other. With this bug, to my understanding combining the results of two independent CRPropa runs does not reduce the MonteCarlo error for the combined results. Although this helps us now to test the code we should definetly change that behavior in the future, i.e. using the random seed from the CRPropa Random class also in the SOPHIA RNG or replace the SOPHIA RNG by the Random-class. |
|
So far, following your @adundovi approach, I can tell that the compiler randomly links to either sophia or sophiaMod. Hence, PhotoPionProduction and PhotoPionProductionMod both use the same version of either sophia or sophiaMod and there is no separate access. |
|
Re: the linker problems I've only taken a brief look at this, so feel free to disregard what I'm saying if you've already thought these points through. As far as I can tell, the functions that This is because the linker only resolves functions by their name; it has no information on where a function is meant to come from. (Consider that (Even so, it is a little bit confusing that the compiler didn't report this as an error, but it seems problematic nonetheless.) |
|
Thanks a lot @antoniusf for looking into this. I indeed have considered this and e.g. renaming |
This is the error I would have expected to occur with unchanged function names... Did you change all names or just a few ones? (To get rid of this error, you would definitely have to change all of them, or at least all that are exported... which are quite a few, it appears??) In any case, I think that renaming all of the exported functions (and constants) is probably the most technically simple option. Unfortunately, this appears to be quite a bit of work, so here's some other possibilities I can think of right now: – I believe that CRPropa currently uses static linking for sophia. In principle, it might be possible to split up the static linking process such that PhotoPionProduction/Mod are linked first against their respective versions of sophia (such that all function names are resolved), and only then linked together with the rest of CRPropa. However, I do not know how to make gcc or cmake do this – maybe you could build PhotoPionProduction in a separate folder, with its own CMakeLists, and then integrate that somehow? – The other option might be to manually load the correct version of sophia at runtime. While this would mean that you don't have to muck around with gcc, it is not platform-independent and requires you to write extra code to load each function. (see here for what this would look like) Edit: platform independence may not even be that big of a deal, since most people probably run on linux?? And you could set it up such that the platform-specific code is only required if you want to test the new sophia. |
|
Update on renaming the functions: Since only the Currently, However, before proceeding with any of the approaches above, it might be prudent to re-evaluate if the capability to switch between sophia and sophiaMod at run-time (as opposed to compile-time) is worth the extra effort and added complexity. Almost every approach I can think of either complicates the build process (step-by-step linking), causes issues with platform independence (dynamic linking, symbol hiding for static libraries), or introduces the possibility for (additional) bugs in the sophia code (renaming all exported symbols). This is in addition to the time necessary to make each of the approaches work, which shouldn't be neglected either. |
|
Sorry for the long absence on this topic. I'm looking into it. |
|
Hi all!
from crpropa import *
ppp = PhotoPionProduction()
c = Candidate()
c.current.setId(nucleusId(1, 1))
c.current.setEnergy(100 * EeV)
c.setCurrentStep(1000 * Mpc)
ppp.process(c)
pppmod = PhotoPionProductionMod()
c = Candidate()
c.current.setId(nucleusId(1, 1))
c.current.setEnergy(100 * EeV)
c.setCurrentStep(1000 * Mpc)
pppmod.process(c)
extern "C" {
void sophiaevent_(int& channel, double& inputenergy, double momentum[][2000],
int id[], int& n, double& redshift, int& photonbackground, double& maxz,
int&, double[], double[]);
}
extern "C" {
void sophiaeventmod_(int& channel, double& inputenergy, double momentum[][2000],
int id[], int& n, double& redshift, int& photonbackground, double& maxz,
int&, double[], double[]);
}
I tested this idea by putting different print output in sophiaevent_ and sophiaeventmod_, and the above python script gives: |
|
Thanks a lot @adundovi for the proposal! Could you push a working example of what you propose? That would help me a lot! |
|
Thanks a lot @adundovi for pushing your example. Based on your changes, I am able / intend to continue with step 1 of the agreed road map above which is to clear up SOPHIA to introduce custom photon fields into CRPropa. This sub-step is to demonstrate that, despite having an internal RNG:
Once we agree on this step, further changes to SOPHIA are much easier to validate as we would not have to rely on statistical tests as @TobiasWinchen has pointed out. The modifications can be found in this branch of my repository. This modification adds the following features:
How non-randomness can be shown:
The output should be: eps = 1.4145363205544225E-003
Can you verify this behaviour? If so, I suggest cleaning out unused parts of the SOPHIA code to make it way more readable and show that we will get the very same output after this. Then, we could move on to manipulating physics inside SOPHIA. |
|
@Froehliche-Kernschmelze OK, I can confirm the behaviour with your git branch: |
|
Thanks for your confirmation. Now I suggest to either (or both!)
Which is the most viable way for you? |
|
If I understand you correctly, the first point is about cleaning/refactoring the FORTRAN code, and the second one is about to introduce new features / change the interface, right? As far as I'm concerned, you can do whatever you want and what you find necessary in the modified PPP as long as we have both versions of PPP available in parallel (the new and the original). However, to ease testing of the modified/new version (when you go beyond cleaning the code, and start to introduce new features/change the interface), you should provide an additional interface/wrapper which, by construction, would use the same usage logic, and produce identical output as the old procedure. If this approach is possible, of course. If not, then we have to do the 1st and the 2nd point separately, in which you will have to wait with the inclusion of the 2nd until the 1st is well tested. |
|
Exactly, the first would introduce changes to the original SOPHIA code by cutting away code that never gets used by CRPropa. This procedure may be lengthy time-wise but worth the review effort as this would reduce the code by approximately a factor of 2, i.e. ~8k lines of code. However, this step can still be done after having introduced the parallel version of SOPHIA. Given that time is the scarcer resource, I would directly go to modifying the newly added interface for sophiamod (leaving the original completely unchanged and in parallel). The idea behind this step is to show that the photon sampling for SOPHIA can be done externally - and if done correctly, exactly reproduces the results of original SOPHIA. Thus follows a detailed yet straight forward explanation of how I intend to do this:
Here comes the critical part:
the output should be eps = 1.0000000000000000E-003 In another terminal (to re-set the RNG counter as both versions call sophia_legacy.f in the end) which should yield sophiaeventmod version This works for each combination of parameters that lie over the threshold CM energy for photo-pion production Summary of this step:
Can this be verified by anyone? |
|
I've checked your changes. And it works as you presented, except one detail, you should probably change the data prefix in the mod version of PPP to use the same data files as the original PPP: if (haveRedshiftDependence)
initRate(getDataPath("PhotoPionProduction/rate_" + fname.replace(0, 3, "IRBz") + ".txt"));
else
initRate(getDataPath("PhotoPionProduction/rate_" + fname + ".txt"));
BTW how do you plan to integrate |
|
Thanks for looking into this! My suggestions on how to proceed / integrate are the following:
Do your concerns evolve around how to provide eps? Because since we are at the same page up to here, this would be the next step(s). These could look like the following
These points are pushed in another branch but I am happy to provide this step-by-step. I can also elaborate on how custom photon sampling would work in CRPropa before continuing. Do you have a preferred way to move on? |
|
Ciao @Froehliche-Kernschmelze , |
|
Does anyone of @adundovi, @MHoerbe or @TobiasWinchen know if there is still something in this PR that was not merged with different pull requests? I tried to quickly go through the discussion again but could come to a fast conclusion. So any hint would be appreciated. |
I remember that we kept this PR open but eventually agreed to never merge it. Our outline was:
NOTE: over the years more features came to life and I had documented their implementation in my PhD thesis. Afaik the following changes are in fact implemented now:
Super good work so far; I am very happy to see this come to life in the main branch. What remains, partly with respect to this PR, is
TL;DR
|
|
Thanks @MHoerbe for the very comprehensive summary. I like the idea to discuss open points of this PR and then move the interesting ones into issues. I will bring it up in the next CRPropa meeting and make the necessary changes afterwards. |
Dear all,
with this pull request I want to provide the first of three parts that I wish to contribute to CRPropa. This approach follows the comments made in #218
Content:
This extension is meant to introduce variable photon fields to CRPropa. Previously, there has only been the CMB and several IRB models. This photon field extension basically follows three concepts:
The photon field may
Changelog (most severe changes):
~/src/PhotonBackground.cpp
~/include/PhotonBackground.h
There are now eight new slots PF1..PF8 that may be included in the simulation. A class has been added to draw a random photon from the photon field that takes the primary, the nucleon/gamma crossection and the photon field's shape into account. This method and all those that it depends on have been inspired by SOPHIA.
Photon fields are now contained in photon field files in ~/share/crpropa/Scaling/ . For detailed information on the format and use of these fields, there is a Wiki page in my forked repository.
~/libs/sophia/sophia_interface.f
The main point of these changes to SOPHIA is to run the code on other photon fields than the CMB and IRB Kneiske (as only these two are implemented so far in SOPHIA). SOPHIA samples photons from these fields with its own photon sampling routine which in this version has been removed and replaced by the sampling method mentioned in the section above. For this, the input and output parameters of SOPHIA changed from
to
Furthermore, SOPHIA has recieved:
~/src/module/PhotonPionProduction.cpp
~/include/module/PhotonPionProduction.h
The flag "haveRedshiftDependence" and all its associated functions have been removed as now all photon fields are represented on grids that are redshift-dependent. The signature of SOPHIA has been updated. The treatment of secondary particles being produced in SOPHIA has been rewritten in a more intuitively readable way. Also, this particle treatment did not allow for more than one hadron to be added among all secondaries produced by SOPHIA. This is valid for the CMB and IRB models as not eanough energy in the center of mass system would be present to produce nucleon/anti-nucleon pairs. However, with more energetic photon fields, these particles may occur. I would like to stress that no kind of anti-nucleon can be produced in the current implementation, even if it were produced.
Tests
As photon sampling is removed from SOPHIA and now provided within CRPropa, SOPHIA's sampling routines and its CRPropa counterparts are expected to produce similar results. Therefore, the following tests have been run according to their hierarchy in the procedure of sampling photons:
Hierarchy of sampling functions (highest to lowest):
Reference simulation setup:
crossection is the deepest layer to be tested. Note that the SOPHIA function that returns the nucleon-gamma xsection does not depend on the nucleon's energy itself but rather on the fact whether or not it deals with a proton or a neutron.
The reference simulation called the crossection function roughly 500.000 times. The sampled photon's energies have been injected into the newly written C++ version of crossection producing another set of ca 500.000 crossection values. These two samples have been inserted into the scipy.stats.ks_2samp Kolmogorov-Smirnov test. A result of D=0.0 and p=1.0 corresponds to total identity in this implementation of the Kolmogorov-Smirnov test. The test achieved a result of D=0.0011, p=0.9183. Note that the functions Ef, Pl and breitwigner are only used by crossection.
In an analogous way, the function gaussInt has been tested, which integrates the function functs over all center-of-mass energies. the function functs is only called by gaussint, whereas functs is the only function calling crossection. The KS test yielded a result of D=9.61e-5, p=0.9999999999 between the SOPHIA output and the C++ version.
prob_eps calculates the non-normalized probability to encounter a photon of the given photon field. Normalization is achieved by dividing the integral of prob_eps over its entire photon energy range. prob_eps is the only function calling gaussInt and get_photonDensity. For better comparison, get_photonDensity has been replaced by an analytical expression returning the photon density of the CMB, i.e. a blackbody with a temprature of 2.73 K. The results of SOPHIA vs C++ yield D=9.61e-5, p=0.9999999999. Results on grid data (and runtimes) are discussed in a later section.
Finally, sample_eps calls prob_eps to sample a photon. This happens in two steps: a) calculate the normalization and the photon energy with which one most likely gets an interaction. b) monte-carlo sample a photon from the field, use norm factor and max probability for this. For comparison, the SOPHIA output has been tested vs the C++ analytical expression for the photon density and vs the grid readout. The KS test yields:
Insights: the grid version and the analytic version are able to produce nearly identical outputs, however less similarity is achieved compared to SOPHIA. In return, SOPHIA operates a dedicated sampling routine that only works for the CMB this making several assumptions about the shape of the field in beforehand. The C++ version do not require any information that is not encoded in the photon field file. One might be interested in the fact that SOPHIA uses its own random number generator whereas the C++ routines have been run with std::srand().
Insights: The C++ small grid / analytic versions perform roughly 20% slower than SOPHIA. A larger grid may reach the analytic version at the expense of runtime. All C++ versions perform similarly but do not reach identity with the SOPHIA version. The two points coming to my mind are:
a) the C++ version of sampling has been done on CRPropa's random number generator, SOPHIA uses its own one
b) SOPHIA operates a dedicated sampling method that only works on the CMB (yet another one with IRB_Kneiske) making assumptions about the CMB shape in beforehand. The C++ version does not rely on such information.
Important notes:
Thanks in advance and best wishes,
Mario