Solving TypeError / ValueError for integrate.quad

[marenaud]

Hello hello,

~So I took a quick look. I have a few questions though. integrate.quad is going to except that the function you pass to it (ie xxint) returns a scalar, but right now it seems to return an array. I had to fill in a few functions in there so make it work, so maybe I added the wrong things. Can you explain more what xxint is supposed to do?~

See below, maybe (hopefully!) that's all that's needed?

[marenaud]

For the optimize.newton part, if I write it as

## The code in the other file works well if kf is a function of theta. The functions to find kf are:
from scipy import optimize
import numpy as np

# bunch of functions omitted

def interpolate_fermi_surface(a, e0, t0, t1, t2):
    """ Interpolate the tight binding energy to find the Fermi
    wavevector at theta. Find zeros of the energy in each direction. """
    kf = lambda theta: optimize.newton(
        polar_energy, x0=np.pi / (2 * a), args=(theta, a, e0, t0, t1, t2), tol=1e-4 / a
    )
    return kf

a = 1.0
e0 = 1.0
t0 = 1.0
t1 = 1.0
t2 = 1.0

kf_theta = interpolate_fermi_surface(a, e0, t0, t1, t2)
print(kf_theta(1.0))

Then it works fine. Is the issue that you want it to work if theta is an array and return one kf value per theta? You could do something like

def interpolate_fermi_surface(a, e0, t0, t1, t2):
    """ Interpolate the tight binding energy to find the Fermi
    wavevector at theta. Find zeros of the energy in each direction. """
    kf = lambda theta_array: np.array(
        [
            optimize.newton(
                polar_energy,
                x0=np.pi / (2 * a),
                args=(theta, a, e0, t0, t1, t2),
                tol=1e-4 / a,
            )
            for theta in theta_array
        ]
    )
    return kf

kf_theta = interpolate_fermi_surface(a, e0, t0, t1, t2)
print(kf_theta([1.0, 2.0]))
[1.51003341 1.56126305]

I suppose that's a little ugly and it might be slow, so you could also do something like

def interpolate_fermi_surface(a, e0, t0, t1, t2):
    """ Interpolate the tight binding energy to find the Fermi
    wavevector at theta. Find zeros of the energy in each direction. """
    kf = lambda theta: optimize.newton(
        polar_energy, x0=np.pi / (2 * a), args=(theta, a, e0, t0, t1, t2), tol=1e-4 / a,
    )
    return np.vectorize(kf)

a = 1.0
e0 = 1.0
t0 = 1.0
t1 = 1.0
t2 = 1.0

kf_theta = interpolate_fermi_surface(a, e0, t0, t1, t2)
print(kf_theta([1.0, 2.0]))

The other file with integrate.quad seems a bit tougher and I think I'd need more context to figure out how to make it work properly.

[If-iNewtonsMore]

Thank you! I've tried this and it still gets stuck on the integrate.quad part because kf_theta ends up being class np.vectorizerather than a function.

For integrate.quad, I'd like to integrate xxint (which is a function of a bunch of other functions) with respect to phi from 0 to infinity (it's currently integrating to 100 until the code works properly). If I were to define sxx(phi, theta, T, B) rather than sxx(theta, T, B) I'll have to pass an argument for phi later on which I suspect will be problematic since it would be an infinite array. That's why I initially used the lambda function, which also doesn't work.

The code after xxint is:

def sigma_xx(T, B): outer_integrand = lambda theta: sxx(theta, T, B) xx, _ = integrate.quad(outer_integrand, 0, np.pi, args=(T, B) prefactors = 2*E*C2/(B*4*np.power(np.pi, 3)) return xx * prefactors

which I suspect will give the same problem

[marenaud]

Hmm in principle all that should matter is that kf_theta is callable. If we look at xxint, it'll return an array for each value of theta rather than a scalar, which throws off integrate.quad, so I assume you'd also want to integrate xxint element-wise for every value of theta?

To first order, is there anything that prevents us from doing the same thing as with kf_theta? i.e. either vectorize or explicitly call integrate.quad for each theta in the array.

At first glance, f_zeta is the problematic one since it assumes theta is an array whereas you'd want to run the whole thing with one value of theta at a time. It doesn't depend on any of your integration variables though, so you could rewrite f_zeta to work for a single value of theta, i. e.

kf = np.arange(5)
theta_array = np.arange(5) * np.pi / 5

def f_zeta(theta):
    N = len(kf)
    zeta_array = np.empty(N)
    for i in range(0, N):
        dtheta = theta[i] - theta[i - 1]
        dkf = np.log(kf[i]) - np.log(kf[i - 1])
        zeta_array[i] = np.arctan(dkf / dtheta)
    return zeta_array

f_zeta_interp = interpolate.interp1d(
    theta_array, f_zeta(theta_array), bounds_error=False, fill_value="extrapolate"
)

There's probably a lot of context I'm missing as to what has to be an array and a scalar at what point of the code. You won't really be able to get around the hard constraint that integrate.quad expects the integrand to be scalar. So for example, something like this will work and run, but it's a mess and I have no idea if I'm understanding the data types at each part of the way :P

import numpy as np
from scipy import integrate, interpolate

kf = np.arange(5)
theta_array = np.arange(5) * np.pi / 5

def f_zeta(theta):
    N = len(kf)
    zeta_array = np.empty(N)
    for i in range(0, N):
        dtheta = theta[i] - theta[i - 1]
        dkf = np.log(kf[i]) - np.log(kf[i - 1])
        zeta_array[i] = np.arctan(dkf / dtheta)
    return zeta_array

f_zeta_interp = interpolate.interp1d(
    theta_array, f_zeta(theta_array), bounds_error=False, fill_value="extrapolate"
)

kf_interp = interpolate.interp1d(
    theta_array, kf, bounds_error=False, fill_value="extrapolate"
)

def func_gamma(T):
    return T

def func_wct(T, B):
    return T + B

def vx(theta):
    return np.cos(theta - f_zeta_interp(theta))

def vy(theta):
    return np.sin(theta - f_zeta_interp(theta))

def G2(theta, T, B):
    return (func_gamma(T) * np.sin(4 * theta) + 4 * theta) / (4 * func_wct(T, B))

def xxint(theta, phi, T, B):
    return vx(theta - phi) * np.exp(G2(theta - phi, T, B)) * np.exp(-G2(theta, T, B))

def xyint(phi, theta, T, B):
    return vy(theta - phi) * np.exp(G2(theta - phi, T, B)) * np.exp(-G2(theta, T, B))

def sxx(theta, T, B):
    integral = integrate.quad(lambda phi: xxint(theta, phi, T, B), 0, 100)[0]
    ans = kf_interp(theta) * vx(theta) * integral
    print(ans)
    return ans

def sigma_xx(T, B):
    outer_integrand = lambda theta: sxx(theta, T, B)
    xx, _ = integrate.quad(outer_integrand, 0, np.pi)
    prefactors = 2 * E * C2 / (B * 4 * np.power(np.pi, 3))
    return xx * prefactors

sigma_xx(1.0, 1.0)

[If-iNewtonsMore]

The integration still seems to be throwing errors or just giving up. Just for more contect, the full code (including the calculation of the kf array is:

`import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate import tight_binding as tb from scipy import optimize from scipy.optimize import curve_fit from numba import njit

def sctif():

Aiso=4
Aani=0.004
Biso=2
Bani=0.02

E=1.602176487e-19;
C= 13.26*10**(-10)/2;
C2=2*np.pi/C;

TB model for kf

a  = 3.77e-10 # FROM TOM CROFT ARXIV

# Some pyplot to create a six panel figure
fig, ax=plt.subplots(1)

# doping
dopings = 0.20

##the lee-hone tight binding parameters
e0, t0, t1, t2 = tb.lee_hone_parameters(dopings)

    # Define a load of points over the BZ
kx = np.linspace(-np.pi, np.pi, 500)/a
ky = np.linspace(-np.pi, np.pi, 500)/a
# and meshgrid them so we can go through each point in k-space one at a time
KX, KY = np.meshgrid(kx, ky)

# Compute the enery surface 
ENERGY = tb.energy(KX, KY, a, e0, t0, t1, t2)

# Define the range of thetas
thetas = np.linspace(0, np.pi*2, 100)

# This returns all of the relevant quantities as a function of theta
th, kf, vf, kx, ky, vx, vy, lx, ly, af, al = tb.compute_all(thetas, a, e0, t0, t1, t2)

# Plot 10 energy contours
ax.contour(KX, KY, ENERGY, levels=10, alpha=0.2)
# Plot the Fermi surface (contour at zero where energy=e0)
ax.contour(KX, KY, ENERGY, 0)

# Plot computed kx, ky to check it correctly matches the 0 contour as expected
ax.plot(kx, ky, 'k--', linewidth=0.5)

    # The l-curve (only contains anisotropy due to v_F)
# Multiplied by 1e4 for visibility
ax.plot(lx*1e4, ly*1e4, 'r-', linewidth=0.5)

plt.show()

###########

@njit()
def energy(kx, ky, a, e0, t0, t1, t2):
     """Tight binding energy as a function of kx and ky"""
     return     e0 - 2*t0*(np.cos(kx*a) + np.cos(ky*a)) - 4*t1*(np.cos(kx*a) * np.cos(ky*a)) - 2*t2*(np.cos(2*kx*a) + np.cos(2*ky*a))

@njit()
def polar_energy(r, theta, a, e0, t0, t1, t2):
    """ A wrapper for the tight binding energy in
        ar coordinates."""

    kx, ky = polar_to_cartesian(r, theta)
    return energy(kx, ky, a, e0, t0, t1, t2)

@njit()
def polar_to_cartesian(r, theta):
    """ Utility function for converting from polar to
    cartesian coordinates """
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    return x, y

def interpolate_fermi_surface(a, e0, t0, t1, t2):
    """ Interpolate the tight binding energy to find the Fermi
    evector at theta. Find zeros of the energy in each direction. """
    kf = lambda theta: optimize.newton(
        polar_energy, x0=np.pi / (2 * a), args=(theta, a, e0, t0, t1, t2), tol=1e-4 / a,)
    return np.vectorize(kf)

func_kf = interpolate_fermi_surface(a, e0, t0, t1, t2)
# print(type(func_kf))

alpha=lambda T: Aiso+Biso*T**2;
beta=lambda T: Aani+Bani*T;

func_gamma=lambda T: beta(T)/(2*alpha(T)+beta(T))
func_wct=lambda T, B:((1-func_gamma(T))/alpha(T))*(B/45)###I'm guessing the 45 is the max field they used

# func_kf=lambda theta: k00*(1+k40*np.cos(4*theta))   This is kf for Tl2201
func_tau=lambda theta, T, B: func_wct(T,B)/(1+func_gamma(T)*np.cos(4*theta))

def make_zeta(kf, dphi=1e-8):
    f = lambda phi: np.log(kf(phi))
    g = lambda phi: np.arctan( (f(phi+dphi) - f(phi-dphi))/(2*dphi) )
    return g

f_zeta = make_zeta(func_kf)

def vx(theta):
    return func_kf(theta)*np.cos(theta - f_zeta(theta))

def vy(theta):
    return func_kf(theta)*np.sin(theta - f_zeta(theta))

def G2(theta,T,B):
    return (func_gamma(T)*np.sin(4*theta)+4*theta)/(4*func_wct(T, B))

def xxint(phi,theta,T,B):
    return vx(theta-phi)*np.exp(G2(theta-phi, T, B))*np.exp(-G2(theta, T, B))

def xyint(phi,theta,T, B):
    return vy(theta-phi)*np.exp(G2(theta-phi, T, B))*np.exp(-G2(theta, T, B))

# def sxx(theta, T, B):
#     integrand = lambda phi: xxint
#     print(integrand)
#     integral, _ = integrate.quad(integrand, 0, 100)
#     ans = vx(theta)*integral
#     return ans

def sxx(theta, T, B):
    integral = integrate.quad(lambda phi: xxint(phi, theta, T, B), 0, 100)[0]
    ans = vx(theta)*integral
    return ans

# def sigma_xx(phi,theta,T, B):
#     outer_integrand = lambda theta: sxx(theta, T, B)
#     xx, _ = integrate.quad(outer_integrand, 0, np.pi)[0]
#     prefactors = 2*E*C2/(B*4*np.power(np.pi, 3))
#     return xx * prefactors

# def sxy(theta, T, B):
#     integrand = lambda phi: xyint(theta, phi, T, B)
#     intr,_ = integrate.quad(integrand, 0, 100)[0]#np.inf
#     ans = vx(theta)*intr
#     return ans

# print(type(sxx))

# def sigma_xy(T, B):
#   outer_integrand = lambda theta: sxy(theta, T, B)
#   xy, _ = integrate.quad(outer_integrand, 0, np.pi)[0]
#   prefactors = 2*E*C2/(B*4*np.power(np.pi, 3))
#   return xy * prefactors

B, T = 1, 1
theta = np.linspace(0, np.pi*2, 100)
zeta = f_zeta(theta)

    # Plot intermediate quantitites
#   #
fig, axes = plt.subplots(2, 3)

axes[0, 0].plot(theta, func_kf(theta))
axes[1, 0].plot(theta, vx(theta))
axes[1, 0].plot(theta, vy(theta))
axes[0, 1].plot(theta, zeta)
axes[1, 1].plot(theta, func_tau(theta, 10, 10))
axes[0, 2].plot(theta, sxx(theta,10, 10))

plt.show()

#   # Plot some results
# #     #
# fig, axes = plt.subplots(2, 2)
# # # start = time.time()

# for T in [0, 25, 50, 75, 100, 125, 150]:
    # for B in np.arange(3, 45, 1):
# #          #print(time.time()-start, T, B)
            # xx = sigma_xx(T, B)
            # xy = sigma_xy(T, B)

            # Rxx = xx/(xx*xx + xy*xy)
            # Rxy = xy/(xx*xx + xy*xy)
            # Rh = xy/(xx**2)

            # axes[0,0].plot(B, (xx),'k.')
            # axes[0,1].plot(B, (xy*1e8), 'k.')
            # axes[0,1].set_ylim([-1,10])
            # axes[1,0].plot(B, Rxx*1e8, 'k.')
            # axes[1,1].plot(B, Rxy*1e8,'k.')
            # print(xy)

# plt.show()

if name == "main": sctif() `

I'd rather have kf defined as a function of theta rather than an array of values, which is also why I've changed the f_zeta() function a little. Regardless of whether I'm using arrays or functions, the sxx() integral is being the most problematic. Eventually I'd like this to run from 0-> infinity (maybe when the weather is more agreeable and my poor laptop can run for fifteen thousand years without overheating...)

I think the important thing here is kf. kf is essentially a shape. In the first iteration of the code it was a very simple shape given by kf = k0*(1+k40*cos(4*theta) which was nice because it worked, however the shape is now a much more complicated sum of stupidly high order sines and cosines and still isn't accurate; the actual shape is given by interpolate_fermi_surface

When I run the above code I get the following error:

`run approx_SCTIF.py Traceback (most recent call last):

File "approx_SCTIF.py", line 256, in sctif()

File "approx_SCTIF.py", line 221, in sctif axes[0, 2].plot(theta, sxx(theta,10, 10))

File "/approx_SCTIF.py", line 174, in sxx integral = integrate.quad(lambda phi: xxint(phi, theta, T, B), 0, 100)[0]

File "/home/caitlin/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 352, in quad points)

File "/home/caitlin/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 463, in _quad return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)

TypeError: only size-1 arrays can be converted to Python scalars`

If I run your suggested code using numerical values for kf and theta (from the interpolate bit), it gives up:

run test_SCTIF.py test_SCTIF.py:165: IntegrationWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. integral = integrate.quad(lambda phi: xxint(theta, phi, T, B), 0, 100)[0] /test_SCTIF.py:165: IntegrationWarning: The integral is probably divergent, or slowly convergent. integral = integrate.quad(lambda phi: xxint(theta, phi, T, B), 0, 100)[0] /test_SCTIF.py:173: IntegrationWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. xx, _ = integrate.quad(outer_integrand, 0, np.pi)

So I think the best way forward is to try and solve the integral for when kf and zeta are functions of theta (i.e solve the size-1 array nightmare!!) Perhaps it's just a definition problem or something; I can't quite see why it wouldn't work for the function-defined kf, other than there being an issue with having an np.vectorize class...