Compute numerical integration setting the axis over I want to integrate Python

Question

I would like to numerically integrate an array d dimensional for instance named c.

This hyper-surface has to be integrated using some axis with a specific increment.

Let's say that these particular axis are:

• x[1]
• x[2]
• x[d]

I wrote a function that compute it in 2 d:

from numpy import*
import scipy.integrate as scint


def int2d(c,x,y):
    g=[]
    a=arange(0,size(y))
    for i in a:
        g.append(scint.simps(c[i],x))   
    return scint.simps(g,y)

This works,

How do I extend it to a multidimensional input array?

I need it because I would like to compute the hyper volume of some function of a histogram.

Answer 1

You can use recursion.

Also from scipy.integrate.simps it says the assumed axis is the last axis.

import numpy
import scipy.integrate

def intNd(c,axes):
    ''' c is the array
    axes is a list of the corresponding coordinates
    ''' 
    assert len(c.shape) == len(axes)
    assert all([c.shape[i] == axes[i].shape[0] 
                for i in range(len(axes))])
    if len(axes) == 1:
        return scipy.integrate.simps(c,axes[0])
    else:
        return intNd(scipy.integrate.simps(c,axes[-1]),axes[:-1])

You may also consider improving the efficiency for integrating large complex arrays by integrating along the longest dimension first. However. not integrating along the last axis, which is the fastest, may impose some penalty. I did not look into these details myself.