I think I should extend the answer with more specifics.

```
with pm.Model() as model:
a = pm.Uniform('a', 0, 60)
R = pm.Uniform('R', 0, 60)
N0 = pm.Uniform('N0', -2, 2)
f = pm.Uniform('f', -1, 1)
d = pm.Uniform('d', 0, 50)
B0 = pm.Uniform('B0', 0, 1000)
alpha = pm.Uniform('alpha', -np.pi/2, np.pi/2)
xc = pm.Uniform('xc', 0, 113)
yc = pm.Uniform('yc', 0, 53)
```

#This is the problem I must solve, strategy #1 = using pm.Deterministic. Please note that xv and yv are constant matrices define outside the function. So the output mag is a matrix of the same dimension.

```
def torusmodel(input):
a, R, N0, f, d, B0, alpha, xc, yc = input
x=(xv-xc+0.5)*np.cos(alpha) + (yv-yc+0.5)*np.sin(alpha)
y=-(xv-xc+0.5)*np.sin(alpha) + (yv-yc+0.5)*np.cos(alpha)
xr=np.sqrt(x**2+d**2) - R
rho=np.sqrt(xr**2 + y**2)
u=np.sqrt(x**2 + d**2)
costh=xr/rho
Nt=N0*(1+f*(rho/a)**2)
mag = (x*(1+g*(costh)) - 2*Nt*d*y/u)*(-1)*B0*np.exp((-1)*(rho/a)**2)/u
return mag
with model:
torusmodel_pm = pm.Deterministic('torusmodel_pm', torusmodel(a,R,N0,f,g,d,B0,alpha,xc,yc))
with model:
observations = pm.Normal( "obs", mu=torusmodel_pm, tau=1/0.1, observed=data.values)
```

#########################

#This is strategy #2 = using decorators. Please note that xv and yv are constant matrices define outside the function. So the output mag is a matrix of the same dimension.

```
@T.compile.ops.as_op(itypes=[T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector, T.dvector],otypes=[T.dvector])
def torusmodel(input):
a, R, N0, f, d, B0, alpha, xc, yc = input
x=(xv-xc+0.5)*np.cos(alpha) + (yv-yc+0.5)*np.sin(alpha)
y=-(xv-xc+0.5)*np.sin(alpha) + (yv-yc+0.5)*np.cos(alpha)
xr=np.sqrt(x**2+d**2) - R
rho=np.sqrt(xr**2 + y**2)
u=np.sqrt(x**2 + d**2)
costh=xr/rho
Nt=N0*(1+f*(rho/a)**2)
mag = (x*(1+g*(costh)) - 2*Nt*d*y/u)*(-1)*B0*np.exp((-1)*(rho/a)**2)/u
return mag
with model:
```

`observations = pm.Normal( "obs", mu=torusmodel, tau=1/0.1, observed=data.values)`

Both approaches fail, with obscure errors. At this point, I am not sure if this is related to: (1) the definition of the deterministic function or (2) the fact that I am building matrices of PyMC3 random variables. In any case, this simple approach do work.

```
# this dummy example that works fine
def torusmodel_dummy(a,R,N0,f,g,d,B0,alpha,xc,yc):
return a + R + N0 + f + g + d + B0 + alpha + xc + yc
with model:
torusmodel_pm = pm.Deterministic('torusmodel_pm', torusmodel_dummy(a,R,N0,f,g,d,B0,alpha,xc,yc))
```

Any feedback is welcomed!