Using scipy to perform discrete integration of the sample

Question

I am trying to port from labview to python.

In labview there is a function "Integral x(t) VI" that takes a set of samples as input, performs a discrete integration of the samples and returns a list of values (the areas under the curve) according to Simpsons rule.

I tried to find an equivalent function in scipy, e.g. scipy.integrate.simps, but those functions return the summed integral across the set of samples, as a float.

How do I get the list of integrated values as opposed to the sum of the integrated values?

Am I just looking at the problem the wrong way around?

Typecast

Answer 1

I think you may be using scipy.integrate.simps slightly incorrectly. The area returned by scipy.integrate.simps is the total area under y (the first parameter passed). The second parameter is optional, and are sample values for the x-axis (the actual x values for each of the y values). ie:

>>> import numpy as np
>>> import scipy
>>> a=np.array([1,1,1,1,1])
>>> scipy.integrate.simps(a)
4.0
>>> scipy.integrate.simps(a,np.array([0,10,20,30,40]))
40.0

I think you want to return the areas under the same curve between different limits? To do that you pass the part of the curve you want, like this:

>>> a=np.array([0,1,1,1,1,10,10,10,10,0])
>>> scipy.integrate.simps(a)
44.916666666666671
>>> scipy.integrate.simps(a[:5])
3.6666666666666665
>>> scipy.integrate.simps(a[5:])
36.666666666666664

Dermatology

yeah that is what I wanted to achieve. I guess I was expecting there would be a 'quicker' way using the magic of maths. As an aside, is there a reason you coerce the Python list to a numpy.array? – Typecast

@Typecast - If you have a bunch of start and end points you could do it like this:[scipy.integrate.simps(y[s:e]) for s,e in [(s1,e1),(s2,e2),...etc.]]. It is likely faster to use numpy arrays, and for some functions they are neccessary, particularly when using special axis features etc. But in this example it works on regular lists :) – Dermatology

Thanks! I am not so at sea that I cannot use a loop, but thanks for the tip. And the numpy arrays thing is as I thought. – Typecast

Answer 2

There is only one method in SciPy that does cumulative integration which is scipy.integrate.cumtrapz() which does what you want as long as you don't specifically need to use the Simpson rule or another method. For that, you can as suggested always write the loop on your own.

Hartwig

Answer 3

If you want to integrate discrete data, you can just use the trapezoidal rule. As was stated in the previous answer, you can use scipy.integrate.cumtrapz(), but actually, you can write your own simple variant:

def integ(x,y):
    return x[:-1], (x[:-1] - x[1:]) * (y[1:] + y[:-1])/2

This should basically do the same, but somehow I get a too-smooth line with scipy.integrate.cumtrapz() (which is strange because I checked its source code, and there is no smoothing applied at all), so I found the direct solution above to be more useful.

Delve