Unexpected index error when calling `find_MAP`

[malmaud]

This code generates an error:

import pymc3 as pm

prior=np.array([.5,.5])
emissions=theano.shared(np.array([[.1,.9],[.9,.1]]))

with pm.Model() as m:
    state = pm.Categorical('state', prior)
    observation = pm.Categorical('observation', emissions[state,:], observed=np.array(0))
    map_state = pm.find_MAP(vars=[state])

But the code runs fine when the observed value is changed:

import pymc3 as pm

prior=np.array([.5,.5])
emissions=theano.shared(np.array([[.1,.9],[.9,.1]]))

with pm.Model() as m:
    state = pm.Categorical('state', prior)
    observation = pm.Categorical('observation', emissions[state,:], observed=np.array(1))
    map_state = pm.find_MAP(vars=[state])

The error is

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-226-3b3e234e64a3> in <module>()
      7     state = pm.Categorical('state', prior)
      8     observation = pm.Categorical('observation', emissions[state,:], observed=np.array(0))
----> 9     map_state = pm.find_MAP(vars=[state])

/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/tuning/starting.pyc in find_MAP(start, vars, fmin, return_raw, model, *args, **kwargs)
     84             start), fprime=grad_logp_o, *args, **kwargs)
     85     else:
---> 86         r = fmin(logp_o, bij.map(start), *args, **kwargs)
     87 
     88     if isinstance(r, tuple):

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in fmin_powell(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback, direc)
   2369             'return_all': retall}
   2370 
-> 2371     res = _minimize_powell(func, x0, args, callback=callback, **opts)
   2372 
   2373     if full_output:

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_powell(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, direc, return_all, **unknown_options)
   2454             fx2 = fval
   2455             fval, x, direc1 = _linesearch_powell(func, x, direc1,
-> 2456                                                  tol=xtol * 100)
   2457             if (fx2 - fval) > delta:
   2458                 delta = fx2 - fval

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _linesearch_powell(func, p, xi, tol)
   2262     def myfunc(alpha):
   2263         return func(p + alpha*xi)
-> 2264     alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol)
   2265     xi = alpha_min*xi
   2266     return squeeze(fret), p + xi, xi

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in brent(func, args, brack, tol, full_output, maxiter)
   2001     options = {'xtol': tol,
   2002                'maxiter': maxiter}
-> 2003     res = _minimize_scalar_brent(func, brack, args, **options)
   2004     if full_output:
   2005         return res['x'], res['fun'], res['nit'], res['nfev']

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_scalar_brent(func, brack, args, xtol, maxiter, **unknown_options)
   2033                   full_output=True, maxiter=maxiter)
   2034     brent.set_bracket(brack)
-> 2035     brent.optimize()
   2036     x, fval, nit, nfev = brent.get_result(full_output=True)
   2037     return OptimizeResult(fun=fval, x=x, nit=nit, nfev=nfev,

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in optimize(self)
   1839         # set up for optimization
   1840         func = self.func
-> 1841         xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info()
   1842         _mintol = self._mintol
   1843         _cg = self._cg

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in get_bracket_info(self)
   1813         ### carefully DOCUMENT any CHANGES in core ##
   1814         if brack is None:
-> 1815             xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
   1816         elif len(brack) == 2:
   1817             xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in bracket(func, xa, xb, args, grow_limit, maxiter)
   2194         fa, fb = fb, fa
   2195     xc = xb + _gold * (xb - xa)
-> 2196     fc = func(*((xc,) + args))
   2197     funcalls = 3
   2198     iter = 0

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in myfunc(alpha)
   2261     """
   2262     def myfunc(alpha):
-> 2263         return func(p + alpha*xi)
   2264     alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol)
   2265     xi = alpha_min*xi

/storage/malmaud/anaconda2/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in function_wrapper(*wrapper_args)
    290     def function_wrapper(*wrapper_args):
    291         ncalls[0] += 1
--> 292         return function(*(wrapper_args + args))
    293 
    294     return ncalls, function_wrapper

/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/tuning/starting.pyc in logp_o(point)
     74 
     75     def logp_o(point):
---> 76         return nan_to_high(-logp(point))
     77 
     78     def grad_logp_o(point):

/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/blocking.pyc in __call__(self, x)
    121 
    122     def __call__(self, x):
--> 123         return self.fa(self.fb(x))

/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/model.pyc in __call__(self, state)
    471 
    472     def __call__(self, state):
--> 473         return self.f(**state)
    474 
    475 

/storage/malmaud/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.pyc in __call__(self, *args, **kwargs)
    869                     node=self.fn.nodes[self.fn.position_of_error],
    870                     thunk=thunk,
--> 871                     storage_map=getattr(self.fn, 'storage_map', None))
    872             else:
    873                 # old-style linkers raise their own exceptions

/storage/malmaud/anaconda2/lib/python2.7/site-packages/theano/gof/link.pyc in raise_with_op(node, thunk, exc_info, storage_map)
    312         # extra long error message in that case.
    313         pass
--> 314     reraise(exc_type, exc_value, exc_trace)
    315 
    316 

/storage/malmaud/anaconda2/lib/python2.7/site-packages/theano/compile/function_module.pyc in __call__(self, *args, **kwargs)
    857         t0_fn = time.time()
    858         try:
--> 859             outputs = self.fn()
    860         except Exception:
    861             if hasattr(self.fn, 'position_of_error'):

IndexError: index out of bounds
Apply node that caused the error: Subtensor{int64}(TensorConstant{(2,) of 0.5}, ScalarFromTensor.0)
Toposort index: 3
Inputs types: [TensorType(float64, vector), Scalar(int64)]
Inputs shapes: [(2,), ()]
Inputs strides: [(8,), ()]
Inputs values: [array([ 0.5,  0.5]), 2]
Outputs clients: [[Elemwise{Composite{(Switch((GE(i0, i1) * LE(i0, i2)), log(i3), i4) + log((i5 / i6)))}}[(0, 6)](state, TensorConstant{0}, TensorConstant{1}, Subtensor{int64}.0, TensorConstant{-inf}, Subtensor{int64, int64}.0, Sum{acc_dtype=float64}.0)]]

Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
  File "/storage/malmaud/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py", line 2723, in run_cell
    interactivity=interactivity, compiler=compiler, result=result)
  File "/storage/malmaud/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py", line 2825, in run_ast_nodes
    if self.run_code(code, result):
  File "/storage/malmaud/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py", line 2885, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-226-3b3e234e64a3>", line 7, in <module>
    state = pm.Categorical('state', prior)
  File "/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/distributions/distribution.py", line 27, in __new__
    return model.Var(name, dist, data)
  File "/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/model.py", line 285, in Var
    var = FreeRV(name=name, distribution=dist, model=self)
  File "/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/model.py", line 514, in __init__
    self.logp_elemwiset = distribution.logp(self)
  File "/storage/malmaud/anaconda2/lib/python2.7/site-packages/pymc3/distributions/discrete.py", line 410, in logp
    a = tt.log(p[value])

HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

Thanks for your help!

[ColCarroll]

This is super weird. I'm looking into it -- there are various combinations of observations that do or don't work...

observed=np.array([0, 1]) #error
observed=np.array([0, 1, 1]) #fine!
observed=np.array([0, 0, 1, 1]) #fine!
observed=np.array([0, 0, 0, 1, 1]) #error

at the very least, the output should be consistent.

[ColCarroll]

So this is fun: fmin_powellsets a bracket to start optimizing in a fairly inscrutable manner, searching through floats, but theano is casting everything to int32, and using them as indexes in emissions. For whatever reason, if you initialize at 1, fmin_powell first evaluates logp at -1.618034, which rounds to -1, which pymc3 handles as 1e+100. Initializing at 0 has fmin_powell start at 2.6180339999999998, which rounds to 2, which throws an IndexError.

Trying to figure out the right place to handle this, since we probably should error immediately, or do specific handling for categorical variables -- @fonnesbeck?

[fonnesbeck]

I would not have expected Powell's method to work with a binary variable. For a problem like this, I think brute force is the way forward (scipy.optimize.brute). It will required a tweak to find_MAP, which I can add.

[fonnesbeck]

Using the PR above, running:

map_state = pm.find_MAP(vars=[state], fmin=pm.starting.optimize.brute, ranges=[slice(0,2)])

will work. Since there are only two states, brute force is just fine. Would have been nice to automatically find the range of the variable, but I could not think of a good way to do this.

[ColCarroll]

1474

[springcoil]

Wow, that's a bizarre bug. I'll keep an eye out for any weird behaviour like that.

I think specific handling for categorical variables would be a good idea.

On Sun, Oct 23, 2016 at 4:15 PM, Colin notifications@github.com wrote:

So this is fun: fmin_powellsets a bracket to start optimizing in a fairly inscrutable manner, searching through floats, but theano is casting everything to int32, and using them as indexes in emissions. For whatever reason, if you initialize at 1, fmin_powell first evaluates logp at -1.618034, which rounds to -1, which pymc3 handles as 1e+100. Initializing at 0 has fmin_powell start at 2.6180339999999998, which rounds to 2, which throws an IndexError.

Trying to figure out the right place to handle this, since we probably should error immediately, or do specific handling for categorical variables -- @fonnesbeck https://github.com/fonnesbeck?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/pymc-devs/pymc3/issues/1469#issuecomment-255594170, or mute the thread https://github.com/notifications/unsubscribe-auth/AA8DiAhkWun10wfXwrUeIAr3kwjTe_R1ks5q23oQgaJpZM4KcVPl .

Peadar Coyle Skype: springcoilarch www.twitter.com/springcoil peadarcoyle.wordpress.com