"argument must be a callable function" when using scipy.integrate.quad

Question

I've written some code that enables me to take data from DICOM files and separate the data into the individual FID signals.

from pylab import *
import dicom
import numpy as np
import scipy as sp

plan=dicom.read_file("1.3.46.670589.11.38085.5.22.3.1.4792.2013050818105496124")
all_points = array(plan.SpectroscopyData)
cmplx_data = all_points[0::2] -1j*all_points[1::2]
frames = int(plan.NumberOfFrames)
fid_pts = len(cmplx_data)/frames
del_t = plan.AcquisitionDuration / (frames * fid_pts)

fid_list = []
for fidN in arange(frames):
    offset = fidN * fid_pts 
    current_fid = cmplx_data[offset:offset+fid_pts]
    fid_list.append(current_fid)

I'd now like to quantify some of this data so, after applying a Fourier transform and shift I've tried to use Scipy's quad function and encountered the following error:

spec = fftshift(fft(fid_list[0])))
sp.integrate.quad(spec, 660.0, 700.0)



error                                     Traceback (most recent call last)
/home/dominicc/Experiments/In Vitro/Glu, Cr/Phantom 1: 10mM Cr, 5mM Glu/WIP_SV_PRESS_ME_128TEs_5ms_spacing_1828/<ipython-input-107-17cb50e45927> in <module>()
----> 1 sp.integrate.quad(fid, 660.0, 700.0)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
243     if type(args) != type(()): args = (args,)
244     if (weight is None):
--> 245         retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
246     else:
247         retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
307     if points is None:
308         if infbounds == 0:
--> 309             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
310         else:
311             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

error: First argument must be a callable function.

Could somebody please suggest a way to make this object callable? I'm still not particularly clear after reading What is a "callable" in Python?. After reading that I tried

def fid():
    spec = fftshift(fft(fid_list[0])) 
    return spec

Which returned the same error.

Any help would be greatly appreciated, thanks.

Answer 1

As a beginner being interest in FFT of QM wave nature for years, I'd like to contribute a simple coding to the coming Python beginners who might be asking the same on Stack Exchange.

The error syntax seeing is related to how I am attempting to integrate the FFT_psik (i.e. an array of FFT results) over the frequency domain k. The issue here is that np.abs(FFT_psik)**2 is an array (not a callable scalar function) , and when I am attempting to multiply an array by array inside the integrand, things get messy for the quad function mathematically.

I had ever encountered the same syntax problem, can be outlined here:

def integrand(k_value):
    return k_value**2 * np.abs(FFT_psik)**2  # FFT_psik is an array of FFT of psi

FFT_psik_2nd_moment_quad, error_quad = quad(integrand, k.min(), k.max(), full_output=True)
print(f"FFT_psik_2nd_moment_quad: {FFT_psik_2nd_moment_quad}")
print(f"Estimated error from quad: {error_quad}")

To resolve it, one needs to interpolate np.abs(FFT_psik)**2 so that given any k value, quad can compute the scalar function value for any k values. The interp1d function from scipy.interpolate is a convenient way to fix the problem.

First, make np.abs(FFT_psik)**2 callable scalar function using interp1d. Next, rewrite my quad integrand as a scalar function and use the interpolation to integrate with quad. Here's an example showing how it works fine with resolving this:

#!/usr/bin/python
import numpy as np
from numpy.fft import fft, fftshift
from scipy.integrate import simps, quad
from scipy.interpolate import interp1d

def psi_x(x, α):
    A = (2*α/np.pi)**(1/4)
    return A * np.exp(-α * x**2)

# Define x, k arrays and constant α
α = 3
x1 = -20  # Adjust based on the desired Δk size
x2 = 20   # Adjust based on the desired Δk size
N = 1000  # Increase for better accuracy in x-domain
x = np.linspace(x1, x2, N)  # Adjust the range and points as needed
Δx = x[1] - x[0]
k = fftshift(np.fft.fftfreq(x.size, Δx))  # frequency grid array in K-space
Δk = k[1] - k[0]  # the desired Δk size is smaller enough to catch the varying of FFT spectrum ROI

# Calculate the wave function in x-space and its FFT spectrum in k-space
psi = psi_x(x, α)
FFT_psik = fftshift(fft(psi))

### Calculate the 2nd moments using 'np.sum' method and 'simps' integration
# FFT_psik_2nd_moment = np.sum(k*k*np.abs(FFT_psik)**2) / np.sum(np.abs(FFT_psik)**2)  # no good for my use case, error is unpredictable large
FFT_psik_2nd_moment = np.sum(k**2 * np.abs(FFT_psik)**2) * Δk  # reliable one
FFT_psik_2nd_momentS = simps(k**2 * np.abs(FFT_psik)**2, k)

print("Second Moment Estimation Alternatives: \n")
print(f"FFT_psik_2nd_moment of np.sum(^2) : {FFT_psik_2nd_moment}")
print(f"Using 'simps' integration method  : {FFT_psik_2nd_momentS} \n")

### Using Gaussian quadrature to compute the 2nd moment of FFT_psik. One needs to define inegrand as a scalar function (not array) before calling 'quad'
# 1. Interpolate |FFT_psik|^2 making it callable 'scalar function'
FFT_psik_inter = interp1d(k, np.abs(FFT_psik)**2, kind='cubic', fill_value="extrapolate")

# 2. Define the integrand for 'quad'
def integrand(k_val):
    return k_val**2 * FFT_psik_inter(k_val)

# 3. Integrate using the saclar integrand() with quad
FFT_psik_2nd_moment_quad, error_quad = quad(integrand, k.min(), k.max())

print("Using Gaussian quad numerical integration: \n")
print(f"FFT_psik_2nd_moment by quad method: {FFT_psik_2nd_moment_quad}")
print(f"Estimated error from Gaussian quad: {error_quad}")  # the side benefit of using quad is providing the reliable 'estimated error'

The actual outputs in my Python environment are as follows:

Second Moment Estimation Alternatives: 

FFT_psik_2nd_moment of np.sum(^2) : 47.399363716988454
Using 'simps' integration method  : 47.39936371698599 

Using Gaussian quad numerical integration: 

FFT_psik_2nd_moment by quad method: 47.39936294876803
Estimated error from Gaussian quad: 6.517835867177056e-07

By interpolating, one may effectively create a callable scalar function which returns the squared magnitude of FFT_psik for arbitrary k values. This resolves the issues when using quad.

Realizing this might be arriving quite late for docar who asked, being late for even after a decade, hope that this could be beneficial for newcomers to Python mathematics, as it has been for me.

Answer 2

Since you are integrating a function defined only on a grid (i.e., an array of numbers), you need to use one of the routines from "Integration, given fixed samples".

In fact, the first argument of quad() must be a function, something you can call (a "callable"). For instance, you can do:

>>> from scipy.integrate import quad
>>> def f(x):
...     return x**2
... 
>>> quad(f, 0, 1)
(0.33333333333333337, 3.700743415417189e-15)

This can be done instead by sampling f():

>>> from scipy.integrate import simps
>>> x_grid = numpy.linspace(0, 1, 101)  # 101 numbers between 0 and 1
>>> simps(x_grid**2, dx=x_grid[1]-x_grid[0])  # x_grid**2 is an *array*.  dx=0.01 between x_grid values
0.33333333333333337

Your situation is like in this second example: you integrate a sampled function, so it is natural to use one or cumtrapz(), simps() or romb().