GammaLib

The Monte Carlo sampling in GModelSpectralSuperExpPlaw::mc is done using a rejection method that draws random energies for a power law and compares then the super-exponential cut-off model to the power law to decide whether to keep or to reject the energy. Since for the first energy bins the sampling is done deeply in the cut-off of the model, the energy is very often rejected before being accepted. Below some debugging of the Monte Carlo method, where samples indicates the number of random energies that are drawn before one is accepted. This number can go up to millions!

GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=35.7830766905839 GeV samples=786164
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=40.3034500835591 GeV samples=345242
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=46.6816418047913 GeV samples=116586
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=36.3483683005087 GeV samples=284844
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=37.2933148694958 GeV samples=170091
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=31.566807167513 GeV samples=71186
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=31.0769872610411 GeV samples=159230
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=31.3044072005203 GeV samples=45231
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=37.9303234539904 GeV samples=324109
GModelSpectralSuperExpPlaw::mc(30 GeV,72.3259983091833 GeV,0 s (TT)): energy=45.909979097982 GeV samples=6397

I implemented a speed-up by simply re-normalising the random number generator used for rejection:

double acceptance_fraction;
double inv_ecut  = 1.0 / m_ecut.value();
double norm_emin = std::exp(- std::pow(emin.MeV() * inv_ecut, m_index2.value()));
double norm_emax = std::exp(- std::pow(emax.MeV() * inv_ecut, m_index2.value()));
double norm      = (norm_emin > norm_emax) ? norm_emin : norm_emax;
do {
    ...
    acceptance_fraction = std::exp(- std::pow(eng * inv_ecut, m_index2.value()));
} while (ran.uniform() * norm > acceptance_fraction);

That this works is illustrated in the figures below, that show on the left the old results and on the rights the results after implementing the sampling with improved efficiency. The red line indicates the analytical model, the blue histogram the sampled energies.

That the code actually works also for rising power laws is illustrated in the figure below.

I did the same implementation of the exponentially cut-off power law GModelSpectralExpPlaw. Below the same comparison: left is before the code change, centre is after the code change, right is for a power-law that increases with energy. Speed-up is impressive and everything looks fine.

And the final verification for the GModelSpectralExpInvPlaw model, panels below are the same as before. Everything looks okay.

The code was merged into bugfix release 1.7.4 and devel. Close issue now.