improving multidimentional integrals?

[jsitarek]

I think it would be good to improve how the integrals are being done in agnpy, for the moment integration is done with simple traperoid integration in each dimention one by one over a predefined grid. There is no check if the integration has converged, and the points are not selected more densely in the regions where the integrant is varying fast. It would be good to use a dedicated multidimentional integration function with a speficied accuracy.

Connected with it, for the case of BLR sometimes the integrant function goes singular. Since the default binning is 50 bins over 5 decades in distance when you select e.g. r = 0.1 * R_line one of the points falls on top of the BLR shell. Proper multidimentional integral should hopefully solve this automatically

[cosimoNigro]

Hi @jsitarek,

I already thought about this, this is why I left a integrator parameter in all the functions performing SED or absorption calculation.

At the moment you can basically choose between two methods: the np.trapz we actually use in most of the functions and this trapz_loglog I implemented in utils/math.py.
It is based on the same function implementation in naima
https://github.com/zblz/naima/blob/master/src/naima/utils.py#L291
and in gammapy
https://github.com/gammapy/gammapy/blob/master/gammapy/utils/integrate.py#L8
I tried to make it a more efficient but did not really manage.
The gammapy developers tell me a log-log integration allows you to have the same precision with a lower number of grid points.
It is not adaptive though.

There is actually a test function checking the two integration methods (np.trapz and trapz_loglog) against each other, for synch and inverse Compton processes:

I just wanted to say that I already started to think about this 😉, but let's discuss more!

[jsitarek]

Hi @cosimoNigro
I realized one more problem with those integrals. They take really increadible amount of system memory.
If you want to compute EC SED over 500 gamma values, 100 mu values and 50 phi values for 100 energies, even with 4 bytes per number it would already end up in 1 GB of memory, and I think in reality python is using much more, because during tests I was running out of memory of a 16 GB machine.
I think we should consider some integration method that takes as a parameter not a grid of points, but a function and boundaries. This will take much less memory and most probably will be also much faster.