Error when trying to pass 'bounds' to scipy curve_fit

Question

I've written a Python script that performs the fitting of the function to the experimental data. The code returns: TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe'' if I pass bounds to curve_fit function from SciPy library. However, without bounds, the code executes and returns fitting parameters.

What is the source of the error and how to correctly pass the bounds to curve_fit function?

I have tried different, including non-zero, values for bounds, and defining it as bounds = (np.array([0, 0,0,0]), np.array([2500, 2500, 2500, 1000])) - did not help.

Here is the code I'm using:

import pandas as pd
import numpy as np
from mpmath import zeta
from scipy.optimize import curve_fit
import os

# defining the function which should describe experimental data:

def sigmaWAL(B, l, lfi, lso, N):
    q = 1.602*10**-19 
    h = 6.626*10**-34
    hbar = h/(2*np.pi) # set of physical constants
    
    lb = np.sqrt(hbar/(4*q*np.abs(B))) # magnetic length
    
    l = l*10**-9
    lfi = lfi*10**-9
    lso = lso*10**-9 # characteristic coherence lengths of the system, recalculated to nm
    
# transforming mpmath.zeta to be compatible with numpy:
    
    npzeta = np.frompyfunc(zeta,2,1)
    
# defining "sub-functions" for the final result

    zeta1 = npzeta(0.5, 0.5 + (lb**2)/(l**2))
    zeta2 = npzeta(0.5, 0.5 + (lb**2)/(lfi**2))
    zeta3 = npzeta(0.5, 0.5 + 4*(lb**2)/(lso**2)+(lb**2)/(lfi**2))

    return N*((q**2)/(4*np.pi*h*lb))*(2*zeta1 + zeta2 - 3*zeta3)

# subtracting the value in the near-zero magnetic field:

def sigmaWAL_normalized(B, l, lfi, lso, N):
    return sigmaWAL(B, l, lfi, lso, N) - sigmaWAL(0.0001, l, lfi, lso, N)

# defining the fitting function, which takes pandas DataFrame with columns 'B' and 'dSIGMA' as a parameter:

def WAL_fit(data):
    Bexp = data['B'] #x values for fitting
    dS = data['dSIGMA'] #y values for fitting
    p0 = [10, 50, 10, 1] # initial guess for parameters

    parameters, covariance = curve_fit(sigmaWAL_normalized,Bexp, dS, p0, bounds = ([0, 0, 0, 0], [2500, 2500, 2500, 1000]), maxfev = 100000)    
    
    l = parameters[0]
    lfi = parameters[1]
    lso = parameters[2]
    N = parameters[3]

    return parameters, covariance, sigmaWAL_normalized(Bexp, l, lfi, lso, N)
    
# iterating over files with experimental data:

for file in os.listdir(os.fsencode('data')):
    
    name = 'data/' + os.fsdecode(file)
    panda = pd.read_csv(name) #reading data to variable 'panda'
    fitting = WAL_fit(panda)

The experimental data does not include B = 0, it spans from B = 0.002 to B = 1 with step of 0.004. dSIGMA takes negative values.

The fitting function is defined as eq. 1 in H. Nakamura et al., Nature Communications 11, 1161 (2020).

Here is the link to the datafiles