Scipy.root problem using PyMC3

Hi, i’m currently working with PyMC3 and got an error when using scipy.root(). This rises this error: ‘error: Result from function call is not a proper array of floats.’. When i call the function dl(x,L,C) alone (out of PyMC3), giving L and C fixed, this works fine. Then, i think the problem appears when PyMC3 interacts with scipy.root, since in that line appears the error.

# Model function

def function(q,x,L,C):
    return np.divide(1,1+x)-np.divide(yt(q,L,C),yt(1.,L,C))
def YL(x,L,C):
    return root(function,1.,args=(x,L,C))['x'][0]

def yt(Y,L,C):
    return Y*np.sqrt(np.divide(1-np.divide(1.,3.)*L*Y*np.sqrt(C+Y),1-L*Y*np.sqrt(C+Y)))

def integrand(u,L,C):

    return np.divide(u,np.sqrt(C+u)*yt(u,L,C))

def dl(x,L,C):
    return c*(1+x)*np.sqrt(C)*np.divide(yt(1.,L,C),100.*np.sqrt(h2OmegaR))*(quad(integrand,YL(x,L,C),1.,args=(L,C))[0])

# define the model/function to be fitted.

with pymc3.Model(): 
    L = pymc3.Normal('L', 0., 1., testval= 0.505)
    C = pymc3.Normal('C', 0.00010, 0.0003, testval= 0.00020)
    #M = pymc3.Normal('M', -21, -19, testval= -19.5)
    Delta_m = pymc3.Deterministic('Delta_m',5.*np.log10(dl(x,L,C)) + 25. - 19.26)
    y = pymc3.Normal('y', mu=Delta_m, tau=1.0/np.power(f_error,2), observed = f)
    start = pymc3.find_MAP()
    step = pymc3.NUTS()
    trace = pymc3.sample(100,start=start)

PyMC3 variables are theano tensors, which means you can not apply numpy/scipy functions on them directly.

The most direct (and many cases slow) way to go around it is to wrap the numpy/scipy function as a theano function:
a pymc3 example:

The ideal way is to rewrite the numpy/scipy function as a theano function. It could be a bit involve at times, but you can have a look at the doc for an example (which also finding the root of some function):