Making pymc wrapper for finite element

Questions v5
3

I called the function with “FEmcOperator(h0,h1,h2)” but now they seem like in an infinite loop.Or should I recreate all the pm.model with custom dist? I am sampling the function FEmcOperator as mean (mu) for multivariate normal sampling with covariance, observation “y”, and inputs “a0,a1,a2”.

I also tried to testing code using the guide in this link and it gives me the result that is expected (1x20 numpy array with displacement vector u).

Testing code :

finiteelement_op = fenicsFEmc()
test_out = finiteelement_op(g0,g1,g2)
test_out.eval()

Sampling code :

# Create covariance matrix
covariance = np.eye(len(x_cent)) * perturb

#building pymc model
with pm.Model() as model:
    # Priors
    a0 = pm.Normal('a0', mu=g0, sigma=1.)
    a1 = pm.Normal('a1', mu=g1, sigma=0.5)
    a2 = pm.Normal('a2', mu=g2, sigma=0.2)
    
    h0 = pt.as_tensor_variable(a0)
    h1 = pt.as_tensor_variable(a1)
    h2 = pt.as_tensor_variable(a2)

    # Likelihood
    fe= FEmcOperator(h0,h1,h2)
    likelihood = pm.MvNormal('likelihood', mu=fe, cov = covariance, observed=y)

with model:
    # Use custom number of draws to replace the HMC based defaults
    trace = pm.sample(1000, tune=1000)

# plot the traces
az.plot_trace(trace);