Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement RANLUX random number generator in egs++ #292

Closed
wants to merge 1 commit into from
Closed

Implement RANLUX random number generator in egs++ #292

wants to merge 1 commit into from

Conversation

seirios
Copy link

@seirios seirios commented May 20, 2017

I noticed a lack of the RANLUX random number generator in egs++, so I implemented a class EGS_Ranlux for it based on GSL's integer implementation (gsl_rng_ranlux), which is pretty much the same as the original Fortran implementation in EGSnrc. I haven't done extensive testing, but it seems to be working fine.

I think the next step could be to implement it using Lüscher's floating-point (SSE optimized) code from his own implementation or based on GSL's floating-point implementation (gsl_rng_ranlxs* and gsl_rng_ranlxd*).

@seirios
Copy link
Author

seirios commented May 20, 2017

Sorry for creating a new pull request, but I changed it to be based off my develop branch as well. Couldn't find a way to do it otherwise.

@ftessier
Copy link
Member

ftessier commented Jun 1, 2017

Thank you @seirios! Can you comment on the source code from a licence standpoint? Did you write it, say, based on Lüscher's 1994 paper or egsnrc.mortran, or did you work off Lüscher's code or the GSL implementation?

@seirios
Copy link
Author

seirios commented Jun 3, 2017

I based the code on the GSL implementation. I also looked into egsnrc.mortran, but I'm not sure if any bits made it to the code.

@ftessier ftessier self-assigned this Jun 9, 2017
@mainegra
Copy link
Contributor

mainegra commented May 7, 2020

@ftessier what ever happened to this PR? Were there issues with Open Source licenses? @seirios did you test this implementation and compared it with the current high-resolution modification of ranmar? I am asking now because I think we should make ranlux or the high-resolution version the default random number generator in EGSnrc. There are situations were ranmar simply gives wrong results, albeit not many. But for robustness it is a better approach to have someone switch to the original ranmar intentionally than having users unaware of this issue obtaining wrong results. This is espcially important as more simulations move to low energies and small scales.

@crcrewso
Copy link
Contributor

crcrewso commented May 7, 2020

This is a naive question, but Linux RAND and GCC Rand have seen a lot of improvements over the last few years, As such is RANLUX still necessary? (For backwards consistency? for matching reality?) If so, could we include the whole GSL so we're not maintaining our own code in parallel?

@mainegra
Copy link
Contributor

mainegra commented May 7, 2020

@crcrewso I agree with that. The only issue would be if there is a lot of repetition in terms of the tools that EGSnrc uses. So, one would have to comb through the system removing duplicates such as the random number generator. And, it is also a matter of keeping EGSnrc as lean as possible without dependency to other libraries which might be easily available on some OSs but not in other. And might also have version incompatibilities. See for instance the issues with the QT libraries for our GUIs. Although QT has gotten much better since QT4 and QT5. Anyway, I think this topic should be open for debate.

@rfmthoms
Copy link

rfmthoms commented May 7, 2020 via email

@ftessier
Copy link
Member

ftessier commented May 11, 2020

Agreed that the default generator should be "high-resolution", with decent performance. But ranmar and ranlux are such old generators, no? They don't even make in on any modern list of the optimal rngs these days. If we are going to move to a new default, I suggest what seems to be the current state of the art, the xoshiro generators, specifically the xoshiro256++ generator, clocking in at ~0.86 nanosecond per 64 bits (on an Intel i7-8700B CPU @3.20GHz using gcc 8.3.0), with a periodicity of 2^256 - 1 and no failures on the BigCrush test suite. There is no issue of copyright since the code is placed in the public domain. The code is only about 60 lines in C.

Read the xoshiro paper by Blackman and Vigna: https://arxiv.org/pdf/1805.01407.pdf

@mainegra
Copy link
Contributor

@ftessier the issue of avoiding a potential bug should be separate from the inclusion of a new, better rng in EGSnrc. Setting the default rng to our best available option is a simple modification that avoids serious errors in some situations. Modernizing our rng's should be a new feature PR.

@ftessier
Copy link
Member

Agreed that we should make ranlux the default on the mortran backend. But inasmuch as there is currently no implementation of ranlux on the egs++ side, then implementing ranlux or xoshiro equivalent in terms of simplicity etc. I propose we:

  • test and merge this ranlux implementation by @seirios
  • make ranlux the default on the mortran and egs++ side
  • implement the xoshiro rng in a separate pull request, indeed; and make it the default?

My only concern with switching to ranlux is that simulations will be significantly slower (what is the relative efficiency of ranmar vs ranlux?). That is why I am tempted to jump to xoshiro as the default as soon as possible, to provide both a better and much faster rng.

@ftessier
Copy link
Member

About the "high resolution", this does not have to do with the quality of the rng, but the number of bits. Why not then stick with ranmar, but set the default to "high resolution"? Is that an acceptable default?

@mainegra
Copy link
Contributor

mainegra commented May 11, 2020

The question remains: Who is going to bell the cat? ;-) Meanwhile we can avoid potential errors by simple changing the current defaults. In terms of time penalty, I observed only a 4% cpu time increase when using the high resolution rng with egs_app for a simple photon-only calculation. This is not too onerous in my opinion. But there should be a more detailed assessment.

@ftessier
Copy link
Member

ftessier commented May 11, 2020

Yes, indeed. We change the default to high resolution ranmar! I thought we were talking about changing the default to ranlux... But of course the high resolution ranmar ought to be the default, thanks for the timing test. Upon reviewing my correspondance, I realize that @MartinMartinov had pointed it out in 2015: he solved his discrete sampling problem at small scale by turning on the high resolution option of ranmar. Go for that! So defer this discussion to #594 for changing the default to ranmar high-resolution everywhere, and continue this thread here solely for the egs++ implementation of ranlux by @seirios .

@mainegra
Copy link
Contributor

mainegra commented Jun 8, 2020

About the "high resolution", this does not have to do with the quality of the rng, but the number of bits. Why not then stick with ranmar, but set the default to "high resolution"? Is that an acceptable default?

@ftessier I haven't looked inside the Mortran ranmar implementation to check whether there is the option to set it to high resolution! When I have a chance I will and then perhaps create a PR setting the high-resolution ranmar as default all across EGSnrc/BEAMnrc.

Update: The Mortran ranmar implementation does not have a hires option. It should be easily implementedin a similar fashion as in the C++ side, but in that case it might be better to go with ranlux for both Mortran and C++ apps taking advantage of this PR.

@ftessier
Copy link
Member

ftessier commented Jun 9, 2020

I will merge the @seirios implementation of ranlux on the egs++ side, and make ranlux the default generator on both the mortran and the egs++ frontends, pending ranlux efficiency tests. I will also implement the "high-resolution" (48 bits) ranmar on the mortran side (although I don't like how it is built from two sequential 24-bit numbers from the generator...). Finally, I will propose the xoshiro256++ generator via another pull-request.

@mainegra
Copy link
Contributor

mainegra commented Jun 9, 2020

I like that @ftessier !
Why not xoshiro256+? It's 15% faster with similar properties!

@ftessier
Copy link
Member

ftessier commented Jun 9, 2020

Yes, you are right, we could implement xoshiro256+ for floating point numbers (53 mantissa bits), and perhaps xoshiro256++ if we need 64-bit integers (I don't think we do!). Although the low linearity of the low bits are not considered a major concern in general, according to the authors, so yes let's go with xoshiro256+!

@seirios
Copy link
Author

seirios commented Jun 9, 2020

Hi there, sorry for the late reply.
@mainegra I tested this RANLUX implementation with a small EGS++ program that computed dose deposition in a voxelized grid, and it gave very similar results to DOSXYZnrc for that use case. I haven't compared it to high-res ranmar, as it's been a while since I worked with EGSnrc.
FYI: The reason for this PR is that I like RANLUX for its mathematically proven properties, although it's slower than other generators, which may be undesirable in practice.
Finally, to add my grain of salt, have you looked at the MIXMAX generator (https://mixmax.hepforge.org)? According to ROOT's documentation it's fast (6 ns/call), and it's the default in CLHEP/Geant4, so it may be worth looking at.
Cheers!

@ftessier
Copy link
Member

ftessier commented Jun 9, 2020

Thank you @seirios for this contribution, and the reference to the mixmax generator, will look at it! In another discussion (pymedphys/pymedphys#399), @SimonBiggs mentioned an article reporting (transient) initialization issues with ranlux (and other generators); I was surprised about this, and must admit I have never given much thought to rng initialization. I don't think it matters for EGSnrc (it is only a transient problem), but thought I'd mention it here for the record: https://dl.acm.org/doi/10.1145/1276927.1276928 (interestingly, mixmax addressed the sequential initial seeds explicitly).

@mainegra
Copy link
Contributor

mainegra commented Jun 9, 2020

@seirios thanks for the comments. I also am very fond of ranlux! Used it quite a bit back in the 90's when it was considered state-of-the-art!
But the xoshiro256 generation of rngs seem to have much better statistical qualities and are much faster! Xoshiro256+ is clocked at 0.76 ns per rn!
@ftessier I support your suggestion of going ranlux on both ends (Mortran and C++) and adding the note about the faster, riskier options. And we should do a mini timing study of ranlux with luxury level 1 (default).

@ftessier
Copy link
Member

ftessier commented Jun 10, 2020

Oops, as I am reviewing this code, and as @dworogers pointed out (private communication), this ranlux implementation is a single-precision version. I propose that we implement a 64-bit version instead (or both if we want). Our current related issue is that people working with nano- to micro-scale simulations are having issued with the floating-point grid being too coarse. Here is a recent review on High-Quality Random Number Generators that may prove useful.

@ftessier
Copy link
Member

ftessier commented Jun 10, 2020

I propose we implement the ranlux generator, single and double precision, directly with Lüscher's code (GPL), and use appropriate bindings to interact with it from the mortran and egs++ frontends. And perhaps implement the mixmax generator, after reading the linked above. Or ranlux++ (from Sibidanov in Vicotria!), but given the required assembly codes, it does not seem generally portable (yet): https://github.com/sibidanov/ranluxpp

At ~30 ns per double-precision, I don't think we make ranlux the default generator. I still side with xoshiro256+ as a default, and ranlux as a reference to check simulations where the rng quality is called into question.

Comments welcome!

@mainegra
Copy link
Contributor

@ftessier I see a larger project on RNGs in EGSnrc shaping up .... any computer science student interested ?

@ftessier
Copy link
Member

Another great reference (blog post) on ranlux: https://christoph-conrads.name/faster-ranlux-pseudo-random-number-generators/. I was surprised to learn there that the standard 24-bit ranlux implementation is about twice as slow on X86-64 as a proper 64-bit implementation (see the first table). Really good background information in there too!

@ftessier ftessier added this to the Release 2022 milestone Mar 25, 2021
@ftessier
Copy link
Member

@seirios your develop branch for this pull request has been deleted, apparently. Hence I will close the PR. Please reopen it if you still have the original branch for your EGS_Ranlux class somewhere! As mentioned in the discussion, the implementation was single-precision, but it would remain useful as a starting point for the double-precision version!

@seirios
Copy link
Author

seirios commented Apr 21, 2021

Hi @ftessier , sorry for that. I moved to GitLab some time ago and removed my EGSnrc fork from here. I've now resubmitted the pull request here #702. Cheers!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants