How to find out which values a function f is called when using scipy.integrate.quad?

Question

I have a function f(x) = 1/x^2 and I evaluate the integral from 0 to 1 using scipy.integrate.quad. The scipy.integrate.quad is an adaptive integration routine, and I would like to know in which regions of [0,1] is the function f evaluated. So, for which values of x is the function f called when making an estimate of the integral?

I'm thinking of using a global variable and appending the x values called to it to keep track of the which x values get used. However, I'm not too familiar on how to do this and any help would be greatly appreciated.

The plan is to then plot a histogram to see what regions in the interval [0,1] are evaluated the most.

Answer 1

You could use a decorator class that saves each value x before evaluating the function:

class MemoizePoints:
    def __init__(self, fun):
        self.fun = fun
        self.points = []

    def __call__(self, x, *args):
        self.points.append(x)
        return self.fun(x, *args)

@MemoizePoints
def fun(x):
    return 1 / x**2

quad(fun, a = 1e-6, b = 1.0)

Then, fun.points contains all the x values for which the function f was evaluated.

Here, the decorator @MemoizePoints is just syntactic sugar for

def fun(x):
    return 1 / x**2

fun = MemoizePoints(fun)