Fitting of mpmath.zeta using scipy.optimize.curve_fit with bounds encounters NaNs

Question

This is a follow-up of my previous question. More information, including actual working example, can be found there. After doing some research I decided to post a new question rather than edit the previous one.
I have written a python script that performs a fitting of the function based on mpmath.zeta to the experimental data, using scipy.optimize.curve_fit. When I try to pass the bounds to curve_fit, I encounter the following error:
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''
Without the bounds, code excecutes and returns some (sometimes even reasonable) results.

I have tried to trace the error through the source code of curve_fit and least_squares. The error occurs in the line 836 of least_squares, which writes:

if not np.all(np.isfinite(f0)):
    raise ValueError("Residuals are not finite in the initial point.")

It means that np.isfinite(f0) gets the parameter type it cannot handle. The parameter f0 is a np.array type, and type(f0[1]) returns <class 'mpmath.ctx_mp_python.mpf'>. After re-defining f0 = np.array(f0).astype(float), the error dissapears, but I don't want to modify the source of least_squares. Here is what I've figured out so far:

1. least_squares is called by curve_fit in the following way:

least_squares(func, p0, jac=jac, bounds=bounds, method=method, **kwargs)

with
p0 = np.array([10, 50, 10, 1]) # my initial guess for parameters
bounds = (np.array([0.0, 0.0,0.0,0.0]), np.array([2500.0, 2500.0, 2500.0,1000.0])) # my bounds for fitting parameters
2. func is defined by curve_fit as func = _wrap_func(f, xdata, ydata, transform), where

def _wrap_func(func, xdata, ydata, transform):
    if transform is None:
        def func_wrapped(params):
            return func(xdata, *params) - ydata
    elif transform.ndim == 1:
        def func_wrapped(params):
            return transform * (func(xdata, *params) - ydata)
    else:
        # Chisq = (y - yd)^T C^{-1} (y-yd)
        # transform = L such that C = L L^T
        # C^{-1} = L^{-T} L^{-1}
        # Chisq = (y - yd)^T L^{-T} L^{-1} (y-yd)
        # Define (y-yd)' = L^{-1} (y-yd)
        # by solving
        # L (y-yd)' = (y-yd)
        # and minimize (y-yd)'^T (y-yd)'
        def func_wrapped(params):
            return solve_triangular(transform, func(xdata, *params) - ydata, lower=True)
    return func_wrapped

with xdata, ydata being my experimental values, and f - function that I pass to curve_fit
3. In my case,

def f(B, l, lfi, lso, N):
    q = 1.602*10**-19 #elementary charge in Coulomb
    h = 6.626*10**-34 #Planck constant in J*s (NOT hbar!)
    hbar = h/(2*np.pi) #reduced Planck constant
    
    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 ith 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)

4. The parameter f0, which is passed to np.isfinite(), is defined as:

def fun_wrapped(x):
        return np.atleast_1d(fun(x, *args, **kwargs))
f0 = fun_wrapped(p0)

Since p0 is just an array of floats or integers, it seems to me that the problem is that fun_wrapped(p0) gives an array of mpmath-type objects. Here are my questions:

1. What does _wrap_func(func, xdata, ydata, transform) actually do?
2. Why does f0 contain mpmath-type objects, not floats?
3. Can I overcome this issue without making modifications to least_squares source code?
4. Why is the error present only when I pass bounds to curve_fit?

I tried to make this question as clear and concise as possible, but the issue seems to be quite complex - I am happy to provide more information or edit it if something is not clear.