Theano bracket nesting level error

Questions v3
1

I’m running into an issue with Theano when building a custom hierarchical PyMC3 model using a black box likelihood. When running the model with too many subjects, I receive the following error:

'Compilation failed (return status=1): /Users/matthewmurrow/.theano/compiledir_macOS-10.16-x86_64-i386-64bit-i386-3.8.8-64/tmplbcx0v1h/mod.cpp:35298:32: fatal error: bracket nesting level exceeded maximum of 256. if (!PyErr_Occurred())

The relevant part of the code is shown below. N_subjects is the number of individuals with data I am including in the model. The logl_fun callouts in the pm.Potentials reference the custom likelihood function itself. From looking at other posts, the issue seems to be related to the python loop I’m doing in the PyMC function. Is this correct, and if so, will a Theano scan loop solve the problem?

with pm.Model() as hm:

# prior distributions
mean_mu = pm.Normal('mean_mu0', mu = 1.0, sigma = 2.0) 
sd_mu = pm.Gamma('sd_mu0', alpha = 1.0, beta = 1.0)
mean_mu1 = pm.Normal('mean_mu1', mu = 1.0, sigma = 1.0) 
sd_mu1 = pm.Gamma('sd_mu1', alpha = 1.0, beta = 1.0)
tnd = pm.Uniform("tnd", lower = 0.0, upper = 0.5)
w = pm.Uniform("w", lower = 0.3, upper = 0.7)
b = pm.Uniform("b", lower = 0.1, upper = 1.0)
    
# group level distributions
mu = pm.Normal('mu_c0', mu = mean_mu, sigma = sd_mu, shape = N_subjects)
mu1 = pm.Normal('mu_c1', mu = mean_mu1, sigma = sd_mu1, shape = N_subjects)
    
psi_temp = [0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
psi_temp[0] = w
psi_temp[2] = tnd
psi_pymc = tt.as_tensor_variable(psi_temp)

for ii in range(N_subjects):
    phi_pymc[ii] = tt.as_tensor_variable([mu[ii], 0.0, 0.0, 1.0, b, 0.0, 0.0, np.log10(1.0), np.log10(1.0), np.log10(1.0), 1.0, 0.0, 0.0, 1.0])
    phi_pymc1[ii] = tt.as_tensor_variable([mu1[ii], 0.0, 0.0, 1.0, b, 0.0, 0.0, np.log10(1.0), np.log10(1.0), np.log10(1.0), 1.0, 0.0, 0.0, 1.0])
    
    pm.Potential('c0' + str(ii), logl_fun(phi_pymc[ii], psi_pymc, ii))
    pm.Potential('c1' + str(ii), logl_fun1(phi_pymc1[ii], psi_pymc, ii))
    
step = pm.DEMetropolisZ(tune = 'lambda')
samples = 50000
tune = 0
chains = 2
cores = 2

idata = pm.sample(draws = samples, tune = tune, step = step, compute_convergence_checks = False, chains = chains, cores = cores, return_inferencedata = True, discard_tuned_samples = False)
idata.to_netcdf('test' + '.nc')
2

Yes, you can only nest 256 levels with indexing in Theano, so you will need a scan if you can.

3

With the nesting issue, will I still be able to reference lists outside of the loop? i.e. if I turn the loop into a scan, will I still be able to call mu[ii], for example, without incurring the nesting issue?

4

Yes, you will still be able to index the output.

5

So, I’ve tinkered with this a bit over the last couple days, and I ran into something that confused me a bit. I was trying to avoid using scan by eliminating the loop entirely within the model context. The code is as follows.

with pm.Model() as hm0:

    # prior distributions
    mean_mu = pm.Normal('mean_mu0', mu = 1.0, sigma = 2.0) 
    sd_mu = pm.Gamma('sd_mu0', alpha = 1.0, beta = 1.0)
    mean_mu1 = pm.Normal('mean_mu1', mu = 1.0, sigma = 1.0) 
    sd_mu1 = pm.Gamma('sd_mu1', alpha = 1.0, beta = 1.0)
    tnd = pm.Uniform("tnd", lower = 0.0, upper = 0.5)
    w = pm.Uniform("w", lower = 0.3, upper = 0.7)
    b = pm.Uniform("b", lower = 0.1, upper = 1.0)
        
    # group level distributions
    mu = pm.Normal('mu_c0', mu = mean_mu, sigma = sd_mu, shape = N_subjects)
    mu1 = pm.Normal('mu_c1', mu = mean_mu1, sigma = sd_mu1, shape = N_subjects)
        
    psi_temp = [w, 1.0, tnd, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    psi_pymc = tt.as_tensor_variable(psi_temp)
    
for ii in range(N_subjects):
    phi_pymc[ii*N_phi] = mu
    phi_pymc[1 + ii*N_phi] = 0.0
    phi_pymc[2 + ii*N_phi] = 0.0
    phi_pymc[3 + ii*N_phi] = 1.0
    phi_pymc[4 + ii*N_phi] = b
    phi_pymc[5 + ii*N_phi] = 0.0
    phi_pymc[6 + ii*N_phi] = 0.0
    phi_pymc[7 + ii*N_phi] = np.log10(1.0)
    phi_pymc[8 + ii*N_phi] = np.log10(1.0)
    phi_pymc[9 + ii*N_phi] = np.log10(1.0)
    phi_pymc[10 + ii*N_phi] = 1.0
    phi_pymc[11 + ii*N_phi] = 0.0
    phi_pymc[12 + ii*N_phi] = 0.0
    phi_pymc[13 + ii*N_phi] = 1.0
    
    phi_pymc1[ii*N_phi] = mu1[ii]
    phi_pymc1[1 + ii*N_phi] = 0.0
    phi_pymc1[2 + ii*N_phi] = 0.0
    phi_pymc1[3 + ii*N_phi] = 1.0
    phi_pymc1[4 + ii*N_phi] = b
    phi_pymc1[5 + ii*N_phi] = 0.0
    phi_pymc1[6 + ii*N_phi] = 0.0
    phi_pymc1[7 + ii*N_phi] = np.log10(1.0)
    phi_pymc1[8 + ii*N_phi] = np.log10(1.0)
    phi_pymc1[9 + ii*N_phi] = np.log10(1.0)
    phi_pymc1[10 + ii*N_phi] = 1.0
    phi_pymc1[11 + ii*N_phi] = 0.0
    phi_pymc1[12 + ii*N_phi] = 0.0
    phi_pymc1[13 + ii*N_phi] = 1.0
    
phi_pymc = tt.as_tensor_variable(phi_pymc)
phi_pymc1 = tt.as_tensor_variable(phi_pymc1)
N_sub = tt.as_tensor_variable(np.int32(N_subjects))
    
with hm0:
    
    pm.Potential('c0', logl_fun(phi_pymc, psi_pymc, N_sub ))
    pm.Potential('c1', logl_fun1(phi_pymc1, psi_pymc, N_sub ))
        
    step = pm.DEMetropolisZ(tune = 'lambda')
    samples = 50000
    tune = 0
    chains = 2
    cores = 2
    
    idata = pm.sample(draws = samples, tune = tune, step = step, compute_convergence_checks = False, chains = chains, cores = cores, return_inferencedata = True, discard_tuned_samples = False)
    idata.to_netcdf('test' + '.nc')

The major change is that I create a one dimensional list containing all the parameters (phi_pymc and phi_pymc1) which I then convert into a tensor object outside the model context. I then shove this large list (hundreds of elements long) into the potential. Since this eliminated the loop, I thought it would solve the problem, but I still get the same nesting error.

To clarify, does the nesting error pop up whenever you need to convert a large python object into theano, or is it just a looping thing?

6

It’s the loop thing.