GammaLib

Jürgen, how should we implement GModelSpectralGauss::mc?
(With Google we didn’t find code that works for a given (min, max) energy range, just simple methods for all energies)

Deil Christoph wrote:

Jürgen, how should we implement GModelSpectralGauss::mc?
(With Google we didn’t find code that works for a given (min, max) energy range, just simple methods for all energies)

I guess you then use a rejection method. Build a do-while loop around the random number generator and exit only if the energy is within the valid range. This can of course be costly in case that the energy range is narrow, but I guess very often min and max will enclose the bulk of the distribution.

Progress update - implemented so far:

- GModelSpectralGauss::flux

Still to implement (in progress):

- GModelSpectralGauss::write (as in GModelSpectralExpPlaw::write)
 - GModelSpectralGauss::flux (as used in Fermi Science Tools)
 - GModelSpectralGauss::mc (no progress yet made)

Also, I noticed the Doxygen docstrings in some cases aren’t formatted correctly - I can correct these tomorrow if neccessary!

One more thing that needs to be done that we forgot today: add GModelSpectralGauss to the user manual:
http://gammalib.sourceforge.net/user_manual/modules/model.html?highlight=curvature#spectral-models
https://github.com/gammalib/gammalib/blob/devel/doc/source/user_manual/modules/model.rst

The GModelSpectralGauss::mc() method needs to be improved.

This shows that currently it sometimes generates negative energies (even though I pass a positive emin) and that it sometimes goes into an infinite loop:

import gammalib
model = gammalib.GModelSpectralGauss()
print(model)
emin = gammalib.GEnergy(1, 'MeV')
emax = gammalib.GEnergy(1e6, 'MeV')
time = gammalib.GTime()
ran = gammalib.GRan()
model.mc(emin, emax, time, ran)
print(model.mc(emin, emax, time, ran)) # call this repeatedly

I was checking on the web: http://www.design.caltech.edu/erik/Misc/Gaussian.html

The algorithm given is

float x1, x2, w, y1, y2;

         do {
                 x1 = 2.0 * ranf() - 1.0;
                 x2 = 2.0 * ranf() - 1.0;
                 w = x1 * x1 + x2 * x2;
         } while ( w >= 1.0 );

         w = sqrt( (-2.0 * log( w ) ) / w );
         y1 = x1 * w;
         y2 = x2 * w;

The implement algorithm is

    // Get uniform value between 0 and 1
    double u = ran.uniform();

    // Do ...
    do {
        double x1;
        double w;
        do{
            x1        = 2.0 * u - 1.0;
            double x2 = 2.0 * u - 1.0;
            w         = x1 * x1 + x2 * x2;
        } while (w >= 1.0);

        // Do ...
        w = std::sqrt( (-2.0 * std::log( w ) ) / w );

        // Compute ...
        double val = x1 * w;

        // Map into [emin,emax] range
        energy = (xmax - xmin) * val + xmin;

    } while ((xmin <= energy) && (energy < xmax));

    // Do ...
    do {
        double x1;
        double w;
        do{
            x1        = 2.0 * ran.uniform() - 1.0;
            double x2 = 2.0 * ran.uniform() - 1.0;
            w         = x1 * x1 + x2 * x2;
        } while (w >= 1.0);

        // Do ...
        w = std::sqrt( (-2.0 * std::log( w ) ) / w );

        // Compute ...
        double val = x1 * w;

        // Map to Gaussian parameters
        energy = m_sigma.value() * val + m_mean.value();

    } while ((xmin <= energy) && (energy < xmax));

Work in progress: https://github.com/ellisowen/gammalib/tree/issue_1014_gradients

Implemented analytical gradients (I still need to check/test this)

Also made changes suggested by Jürgen to the mc algorithm in the GModelSpectralGauss::mc() method (although I may make further changes due to the issues Christoph mentioned)

I also still need to add GModelSpectralGauss to the user manual and work on the pull distribution test(tomorrow)

I’d like now to freeze the code for the gammalib-00-08-00 release. Any progress on the remaining issues? Do you think things will be finished by the end of this week?

Many thanks for your helpful comments! I think from my side there were no further changes or additions I planned to make.

Okay. I guess the only thing still missing is a description of the model in the User Manual. A sub section would be needed in

doc/source/user_manual/modules/model.rst

Could you quickly write something up?

Sure - I will do it this afternoon