Broadcastable pattern error using Metropolis in large graph #1983

junpenglao · 2017-04-02T08:51:46Z

I am reworking on my PyMC3 port of Lee and Wagenmakers' Bayesian Cognitive Modeling book, seems there is a bug running Metropolis sampler if the theano graph is large:

import numpy as np
import theano.tensor as tt
import pymc3 as pm

ns = 20
nc = 9
nq = 30
qcc = np.asarray( [[ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  1., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  1.,  0.,  0.,  1.,  1., -1.],
                   [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  0.,  1.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
                   [ 0.,  1.,  0.,  0.,  0., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
                   [ 0.,  0.,  0.,  1., -1.,  0.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  1.,  1.,  1.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
                   [ 0.,  0.,  1., -1.,  0., -1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  0.,  1., -1., -1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  1., -1., -1.,  0.],
                   [ 0.,  0.,  1.,  0., -1., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  1.,  0., -1., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  1.,  0., -1., -1.,  0.,  0.,  0.]] )
qccmat = np.tile(qcc[np.newaxis, :, :], (ns, 1, 1))
y = np.asarray([[1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 0, 0, 0, 0, 0, 0],
               [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 0, 1, 1, 0, 1, 1, 0],
               [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 0, 0, 0, 1, 0],
               [1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
                1, 1, 0, 1, 0, 1, 1, 1],
               [0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0,
                1, 0, 1, 0, 1, 0, 0, 0],
               [1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
                1, 1, 0, 1, 1, 1, 0, 0],
               [1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
                1, 1, 1, 1, 1, 1, 0, 0],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
                1, 1, 1, 1, 0, 1, 1, 1],
               [1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1,
                1, 1, 0, 0, 1, 1, 1, 1],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 0, 0, 0, 0, 0, 0],
               [0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,
                1, 1, 1, 0, 1, 0, 0, 0],
               [0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1,
                1, 1, 0, 1, 0, 0, 1, 0],
               [1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1,
                1, 1, 0, 1, 1, 0, 0, 0],
               [1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 0, 1, 0, 0, 0],
               [1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1,
                1, 1, 0, 0, 1, 0, 0, 1],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 1, 1, 1, 1, 1],
               [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
                1, 1, 0, 1, 1, 1, 1, 1],
               [1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1,
                0, 1, 1, 0, 1, 1, 0, 0],
               [1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1,
                1, 1, 0, 0, 0, 0, 0, 1],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1,
                0, 1, 1, 1, 0, 1, 1, 1]], dtype=int)

vtmp = np.random.randn(ns, 1, nc)**2
with pm.Model() as model3:
    gamma = pm.Uniform('gamma', lower=.5, upper=1)
    gammat = tt.stack([1-gamma, .5, gamma])
    
    v1 = pm.HalfNormal('v1', sd=1, shape=(ns, 1, nc), testval=vtmp)
    s1 = pm.Deterministic('s1', tt.argsort(v1, axis=2))
    smat2 = tt.tile(s1, (1, nq, 1))  # s[1:nc] <- rank(v[1:nc])
    
    # TTB Model For Each Question
    ttmp = tt.sum(qccmat*tt.power(2, smat2), axis=2)
    tmat = -1*(-ttmp > 0)+(ttmp > 0)+1
    
    yiq = pm.Bernoulli('yiq', p=gammat[tmat], observed=y)
    
    trace3 = pm.sample(1e4, step=pm.Metropolis())

This model compiled without problem, and sample without problem in ADVI and NUTS, however, if I set step=pm.Metropolis() the following error returns:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-87778767da2a> in <module>()
     14     yiq = pm.Bernoulli('yiq', p=gammat[tmat], observed=y)
     15 
---> 16     trace3 = pm.sample(1e4, step=pm.Metropolis())

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/arraystep.py in __new__(cls, *args, **kwargs)
     61                 # If we don't return the instance we have to manually
     62                 # call __init__
---> 63                 step.__init__([var], *args, **kwargs)
     64                 # Hack for creating the class correctly when unpickling.
     65                 step.__newargs = ([var], ) + args, kwargs

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/metropolis.py in __init__(self, vars, S, proposal_dist, scaling, tune, tune_interval, model, mode, **kwargs)
    122 
    123         shared = pm.make_shared_replacements(vars, model)
--> 124         self.delta_logp = delta_logp(model.logpt, vars, shared)
    125         super(Metropolis, self).__init__(vars, shared)
    126 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/metropolis.py in delta_logp(logp, vars, shared)
    469 
    470 def delta_logp(logp, vars, shared):
--> 471     [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
    472 
    473     tensor_type = inarray0.type

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/theanof.py in join_nonshared_inputs(xs, vars, shared, make_shared)
    236     replace.update(shared)
    237 
--> 238     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    239     return xs_special, inarray
    240 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/theanof.py in <listcomp>(.0)
    236     replace.update(shared)
    237 
--> 238     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    239     return xs_special, inarray
    240 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/scan_module/scan_utils.py in clone(output, replace, strict, share_inputs, copy_inputs)
    256                                         [],
    257                                         strict,
--> 258                                         share_inputs)
    259 
    260     return outs

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in rebuild_collect_shared(outputs, inputs, replace, updates, rebuild_strict, copy_inputs_over, no_default_updates)
    230     else:
    231         if isinstance(outputs, Variable):
--> 232             cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over)
    233             cloned_outputs = cloned_v
    234             # computed_list.append(cloned_v)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(
---> 96                     [clone_d[i] for i in owner.inputs], strict=rebuild_strict)
     97                 for old_o, new_o in zip(owner.outputs, clone_d[owner].outputs):
     98                     clone_d.setdefault(old_o, new_o)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/gof/graph.py in clone_with_new_inputs(self, inputs, strict)
    240                     remake_node = True
    241         if remake_node:
--> 242             new_node = self.op.make_node(*new_inputs)
    243             new_node.tag = copy(self.tag).__update__(new_node.tag)
    244         else:

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/tensor/elemwise.py in make_node(self, _input)
    201                         "The broadcastable pattern of the "
    202                         "input is incorrect for this op. Expected %s, got %s."
--> 203                         % (self.input_broadcastable, ib)))
    204                 # else, expected == b or expected is False and b is True
    205                 # Both case are good.

TypeError: The broadcastable pattern of the input is incorrect for this op. Expected (True, False, True, False, True, False), got (True, False, True, False, False, False).

I am on the master branch of pymc3. The background of the model could be found in Chapter 18 of the book.

The text was updated successfully, but these errors were encountered:

junpenglao · 2017-04-02T08:53:51Z

Seems to similar to #1304

junpenglao · 2017-04-16T14:38:43Z

Found a solution to this problem.

instead of doing:

v1 = pm.HalfNormal('v1', sd=1, shape=(ns, 1, nc))
s1 = pm.Deterministic('s1', tt.argsort(v1, axis=2))

create a variable and then reshape it resolved the problem

v1 = pm.HalfNormal('v1', sd=1, shape=ns*nc)
s1 = tt.argsort(v1.reshape((ns, 1, nc)), axis=2)

twiecki · 2017-04-17T07:55:26Z

Glad it's working. I don't see why the above code shouldn't work. It could point to a subtle problem in join_nonshared_inputs for higer dimensional inputs, so I'll reopen this.

lerkoah · 2017-05-02T18:37:24Z

Hi,
I have the same problem for my model.

with rbf_model:
    alpha_model = pm.Normal('alpha', mu=init_alpha, sd=200, shape = numberOfBasis)
    Cx_model = pm.Normal('Cx', mu=init_Cx, sd=20, shape = numberOfBasis)
    Cy_model = pm.Normal('Cy', mu=init_Cy, sd=20, shape = numberOfBasis)
    l_model = pm.Lognormal('l', mu = l, sd=20, shape = 1)
    
    sigma_model = pm.HalfNormal('sigma', tau = np.ones(2), shape=(2,2), testval=init_sigma)

    PHI_Re, PHI_Im = Vobs_function(U, l_model, alpha_model, Cx_model, Cy_model)
    V_model = tt.stack([PHI_Re, PHI_Im], axis = 1)

    

    V_obs = pm.MvNormal('V_obs', mu=V_model, cov=sigma_model, observed=V)
    
    n_samples = 5000
    step = pm.Metropolis()
    trace = pm.sample(n_samples, step)

In my case, my model V_obs is a two dimensional vector, i.e., each measurement has two components. The function Vobs_function() returns this two components and then I concatenate (stack).

I don't know where is the shape-error. Can you help me?

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-63-3118a810fbc9> in <module>()
     47 
     48     n_samples = 5000
---> 49     step = pm.Metropolis()
     50     trace = pm.sample(n_samples, step)
     51 #     estimation = pm.find_MAP()

/root/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py in __new__(cls, *args, **kwargs)
     58                 # If we don't return the instance we have to manually
     59                 # call __init__
---> 60                 step.__init__([var], *args, **kwargs)
     61                 # Hack for creating the class correctly when unpickling.
     62                 step.__newargs = ([var], ) + args, kwargs

/root/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py in __init__(self, vars, S, proposal_dist, scaling, tune, tune_interval, model, **kwargs)
     96 
     97         shared = pm.make_shared_replacements(vars, model)
---> 98         self.delta_logp = delta_logp(model.logpt, vars, shared)
     99         super(Metropolis, self).__init__(vars, shared)
    100 

/root/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py in delta_logp(logp, vars, shared)
    424 
    425 def delta_logp(logp, vars, shared):
--> 426     [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
    427 
    428     tensor_type = inarray0.type

/root/anaconda3/lib/python3.6/site-packages/pymc3/theanof.py in join_nonshared_inputs(xs, vars, shared, make_shared)
    198     replace.update(shared)
    199 
--> 200     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    201     return xs_special, inarray
    202 

/root/anaconda3/lib/python3.6/site-packages/pymc3/theanof.py in <listcomp>(.0)
    198     replace.update(shared)
    199 
--> 200     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    201     return xs_special, inarray
    202 

/root/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan_utils.py in clone(output, replace, strict, share_inputs, copy_inputs)
    256                                         [],
    257                                         strict,
--> 258                                         share_inputs)
    259 
    260     return outs

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in rebuild_collect_shared(outputs, inputs, replace, updates, rebuild_strict, copy_inputs_over, no_default_updates)
    230     else:
    231         if isinstance(outputs, Variable):
--> 232             cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over)
    233             cloned_outputs = cloned_v
    234             # computed_list.append(cloned_v)

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(
---> 96                     [clone_d[i] for i in owner.inputs], strict=rebuild_strict)
     97                 for old_o, new_o in zip(owner.outputs, clone_d[owner].outputs):
     98                     clone_d.setdefault(old_o, new_o)

/root/anaconda3/lib/python3.6/site-packages/theano/gof/graph.py in clone_with_new_inputs(self, inputs, strict)
    240                     remake_node = True
    241         if remake_node:
--> 242             new_node = self.op.make_node(*new_inputs)
    243             new_node.tag = copy(self.tag).__update__(new_node.tag)
    244         else:

/root/anaconda3/lib/python3.6/site-packages/theano/tensor/elemwise.py in make_node(self, _input)
    201                         "The broadcastable pattern of the "
    202                         "input is incorrect for this op. Expected %s, got %s."
--> 203                         % (self.input_broadcastable, ib)))
    204                 # else, expected == b or expected is False and b is True
    205                 # Both case are good.

TypeError: The broadcastable pattern of the input is incorrect for this op. Expected (True,), got (False,).

junpenglao · 2017-05-02T20:35:25Z

hmm, it is a bit difficult to say without the data and the parameters. But I dont think you should use sigma_model = pm.HalfNormal('sigma', tau = np.ones(2), shape=(2,2), testval=init_sigma) as your prior for cov. Try pm.LKJCholeskyCov it should work better

lerkoah · 2017-05-16T04:33:50Z

Using LKJCholeskyCov the problem doesn't disappear.

I think that is a problem from the source specificly for the function join_nonshared_inputs. I'm going to make a fork for trying to solve this issue. For me, It is a persona challenge ahha.

giovannidoni · 2017-08-24T17:23:01Z

I am the same problem here when trying to do some hierarchical stuff.

ids = D.index.values.astype(int)

X = D.values[:,1:-2]

with pm.Model() as maxdiff:
    
    mu = pm.Normal('hyper_alpha', mu=0, sd=5, shape=1)
    sigma = pm.HalfCauchy('hyper_sigma', beta=5)  

    alpha = pm.Normal('alpha', mu=mu, sd=sigma, shape=(dd.ids, dd.attr)) 
    coeff = pm.Normal('beta', mu=mu, sd=sigma, shape=(dd.ids, dd.attr)) 
        
    theta = coeff[ids] * X + alpha[ids]
                     
    p = pm.Deterministic('p', T.nnet.softmax(theta))
    
    y_b = pm.Categorical('yl', p=p, observed=D['best'].values)     
    
    step = pm.Metropolis(vars=[alpha, coeff])
    trace = pm.sample(50000, step)

Mocked input have right dimensions:

X = D.values[:,1:-2]
print X.shape

coeff = np.random.normal(size=(dd.ids, dd.attr))
coeff = coeff[ids]
print coeff.shape

print (coeff * X).shape

Output:

(2400, 12)
(2400, 12)
(2400, 12)

The error traceback is the same as previous comments or at least the error reads analogously. The non hierarchical model works fine. Is it a problem generated by the indexing (which are actually used also in previous cases)?

mfansler · 2018-07-13T18:28:21Z

Same issue (crossposting on SO). It's not clear to me how the workaround for other models in this thread translates to my model implementation. If anyone can help with that, I'd greatly appreciate it!

Model

G = 10    # genes
K = 5     # cell types
C_k = 20  # cells per type

data = sample_data(G, K, C_k)

with pm.Model() as capture_efficiency:

    # Genes expression levels per cell type
    mu_gk = pm.Uniform('mu_gk', 1, 1000, shape=(G, K, 1))

    # Cell capture efficiencies
    p_kc = pm.Beta('p_kc', shape=(1, K, C_k), alpha=5, beta=95)

    # Captured transcripts
    N_gkc = pm.Binomial('N_gkc', shape=(G, K, C_k),
                        n=tt.tensordot(mu_gk, np.ones((C_k, 1)), [[2], [1]]),
                        p=tt.tensordot(np.ones((G, 1)), p_kc, [[1], [0]]),
                        observed=data)

    trace = pm.sample(5000, tune=10000, step=pm.Metropolis())

Sample Data

def sample_data(G=1, K=2, C_k=100):
    mu_gk = np.random.randint(1, 1000, size=(G, K))
    p_kc = np.random.beta(5, 95, (K, C_k))
    N_gkc = np.random.binomial(mu_gk[:, :, np.newaxis], p_kc[np.newaxis, :, :])

    return N_gkc

Error and Stack Trace

Traceback (most recent call last):
File "/Applications/PyCharm.app/Contents/helpers/pydev/pydev_run_in_console.py", line 52, in run_file
pydev_imports.execfile(file, globals, locals) # execute the script
File "/Applications/PyCharm.app/Contents/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
File "/Users/mfansler/Projects/pymc3/intro/capture-efficiency-celltypes.py", line 46, in
trace = pm.sample(5000, tune=10000, step=pm.Metropolis(),
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py", line 65, in new
step.__init__([var], *args, **kwargs)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py", line 136, in init
self.delta_logp = delta_logp(model.logpt, vars, shared)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py", line 624, in delta_logp
[logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/theanof.py", line 245, in join_nonshared_inputs
xs_special = [theano.clone(x, replace, strict=False) for x in xs]
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/theanof.py", line 245, in
xs_special = [theano.clone(x, replace, strict=False) for x in xs]
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/scan_module/scan_utils.py", line 247, in clone
share_inputs)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 232, in rebuild_collect_shared
cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates
clone_v_get_shared_updates(i, copy_inputs_over)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates
clone_v_get_shared_updates(i, copy_inputs_over)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates
clone_v_get_shared_updates(i, copy_inputs_over)
[Previous line repeated 9 more times]
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 96, in clone_v_get_shared_updates
[clone_d[i] for i in owner.inputs], strict=rebuild_strict)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/gof/graph.py", line 246, in clone_with_new_inputs
new_node = self.op.make_node(*new_inputs)
File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/tensor/elemwise.py", line 230, in make_node
% (self.input_broadcastable, ib)))

TypeError: The broadcastable pattern of the input is incorrect for this op. Expected (False, False, True), got (False, False, False).

natalieklein · 2019-07-11T16:04:32Z

I know this thread is older, but I just had the same issue and wanted to share what worked for me. Oddly, what fixed it was not flattening anything (tried that, didn't fix it). It turns out I had one variable that had shape = 1 explicitly specified (because it used to have a different shape, otherwise I normally wouldn't explicitly specify it). Taking that shape = 1 away seems to fix Metropolis sampling.

porterjenkins · 2019-09-13T17:17:21Z

@natalieklein Thanks for posting! That seemed to work for me as well. Changing my shapes from (N x 1) --> (N) seemed to fix the issue for me.

twiecki · 2020-07-27T09:56:41Z

Feel free to reopen if still an issue.

It seems like broadcasting information gets lost when applying `pm.make_shared_replacements`, leading to problems with the metropolis sampler. Potentially related issues below: - #1083 - #1304 - #1983 This fix was previously suggested in the following issue: - #3337 It could be that further adaptations are necessary as indicated in the issue. Strangely, this does not seem to lead to problems when using NUTS.

It seems like broadcasting information gets lost when applying `pm.make_shared_replacements`, leading to problems with the metropolis sampler. Potentially related issues below: - pymc-devs#1083 - pymc-devs#1304 - pymc-devs#1983 This fix was previously suggested in the following issue: - pymc-devs#3337 It could be that further adaptations are necessary as indicated in the issue. Strangely, this does not seem to lead to problems when using NUTS.

junpenglao closed this as completed Apr 2, 2017

junpenglao reopened this Apr 2, 2017

junpenglao changed the title ~~Metropolis sampler failed in large graph~~ Broadcastable pattern error using Metropolis in large graph Apr 2, 2017

junpenglao mentioned this issue Apr 16, 2017

SMC #1923

Closed

junpenglao closed this as completed Apr 16, 2017

twiecki reopened this Apr 17, 2017

junpenglao added the Metropolis label May 8, 2017

twiecki closed this as completed Jul 27, 2020

This was referenced Feb 27, 2021

Broadcasting with metropolis sampler leads to issues where NUTS works perfectly. #4491

Closed

Keep broadcasting information in make_shared_replacements #4492

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Broadcastable pattern error using Metropolis in large graph #1983

Broadcastable pattern error using Metropolis in large graph #1983

junpenglao commented Apr 2, 2017

junpenglao commented Apr 2, 2017 •

edited

Loading

junpenglao commented Apr 16, 2017

twiecki commented Apr 17, 2017

lerkoah commented May 2, 2017 •

edited

Loading

junpenglao commented May 2, 2017

lerkoah commented May 16, 2017 •

edited

Loading

giovannidoni commented Aug 24, 2017

mfansler commented Jul 13, 2018

natalieklein commented Jul 11, 2019

porterjenkins commented Sep 13, 2019

twiecki commented Jul 27, 2020

Broadcastable pattern error using Metropolis in large graph #1983

Broadcastable pattern error using Metropolis in large graph #1983

Comments

junpenglao commented Apr 2, 2017

junpenglao commented Apr 2, 2017 • edited Loading

junpenglao commented Apr 16, 2017

twiecki commented Apr 17, 2017

lerkoah commented May 2, 2017 • edited Loading

junpenglao commented May 2, 2017

lerkoah commented May 16, 2017 • edited Loading

giovannidoni commented Aug 24, 2017

mfansler commented Jul 13, 2018

Model

Sample Data

Error and Stack Trace

natalieklein commented Jul 11, 2019

porterjenkins commented Sep 13, 2019

twiecki commented Jul 27, 2020

junpenglao commented Apr 2, 2017 •

edited

Loading

lerkoah commented May 2, 2017 •

edited

Loading

lerkoah commented May 16, 2017 •

edited

Loading