Bizarre dtype error

[dfm]

I came across some strange bugginess in the likelihood calculation when the dtypes of the data are inconsistent. Here are two examples:

First, if the dtypes are not the same for the different bands, numba throws an error:

import isochrones
import numpy as np

mod = isochrones.SingleStarModel(isochrones.get_ichrone("mist", bands=["G", "BP"]),
                                 G=(np.float32(12.5), np.float32(0.1)),
                                 BP=(np.float64(12.5), np.float64(0.1)))

np.random.seed(42)
pars = np.random.rand(mod.n_params)
mod.mnest_prior(pars, None, None)
mod.lnpost(pars)

Executing this, I get the following error:

---------------------------------------------------------------------------
TypingError                               Traceback (most recent call last)
<ipython-input-36-83c12f9ac8af> in <module>
      9 pars = np.random.rand(mod.n_params)
     10 mod.mnest_prior(pars, None, None)
---> 11 mod.lnpost(pars)

~/anaconda3/lib/python3.6/site-packages/isochrones/starmodel.py in lnpost(self, p, **kwargs)
    507         if not np.isfinite(lnpr):
    508             return -np.inf
--> 509         return lnpr + self.lnlike(p, **kwargs)
    510 
    511     def lnlike(self, p, **kwargs):

~/anaconda3/lib/python3.6/site-packages/isochrones/starmodel.py in lnlike(self, pars)
   1511                              *self.ic.model_grid.interp.index_columns,
   1512                              self.ic.bc_grid.interp.grid,
-> 1513                              *self.ic.bc_grid.interp.index_columns)
   1514 
   1515         if 'parallax' in self.kwargs:

~/anaconda3/lib/python3.6/site-packages/numba/dispatcher.py in _compile_for_args(self, *args, **kws)
    346                 e.patch_message(msg)
    347 
--> 348             error_rewrite(e, 'typing')
    349         except errors.UnsupportedError as e:
    350             # Something unsupported is present in the user code, add help info

~/anaconda3/lib/python3.6/site-packages/numba/dispatcher.py in error_rewrite(e, issue_type)
    313                 raise e
    314             else:
--> 315                 reraise(type(e), e, None)
    316 
    317         argtypes = []

~/anaconda3/lib/python3.6/site-packages/numba/six.py in reraise(tp, value, tb)
    656             value = tp()
    657         if value.__traceback__ is not tb:
--> 658             raise value.with_traceback(tb)
    659         raise value
    660 

TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Invalid use of Function(<built-in function getitem>) with argument(s) of type(s): ((float32, float64), int64)
 * parameterized
In definition 0:
    All templates rejected with literals.
In definition 1:
    All templates rejected without literals.
In definition 2:
    All templates rejected with literals.
In definition 3:
    All templates rejected without literals.
In definition 4:
    All templates rejected with literals.
In definition 5:
    All templates rejected without literals.
In definition 6:
    All templates rejected with literals.
In definition 7:
    All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at /Users/dforeman/anaconda3/lib/python3.6/site-packages/isochrones/likelihood.py (88)

File "../../../../../anaconda3/lib/python3.6/site-packages/isochrones/likelihood.py", line 88:
def star_lnlike(pars, index_order,
    <source elided>
    for i in range(len(mag_vals)):
        val = mag_vals[i]
        ^

This is not usually a problem with Numba itself but instead often caused by
the use of unsupported features or an issue in resolving types.

To see Python/NumPy features supported by the latest release of Numba visit:
http://numba.pydata.org/numba-doc/dev/reference/pysupported.html
and
http://numba.pydata.org/numba-doc/dev/reference/numpysupported.html

For more information about typing errors and how to debug them visit:
http://numba.pydata.org/numba-doc/latest/user/troubleshoot.html#my-code-doesn-t-compile

If you think your code should work with Numba, please report the error message
and traceback, along with a minimal reproducer at:
https://github.com/numba/numba/issues/new

Second, I would expect the following to give the same answer to at least floating point precision, but I get and error of 0.04 which (while small) is much larger than it should be and I find that it makes real differences in the results I get!

import isochrones
import numpy as np

mod1 = isochrones.SingleStarModel(isochrones.get_ichrone("mist", bands=["G", "BP"]),
                                 G=(np.float32(12.5), np.float32(0.1)),
                                 BP=(np.float32(12.5), np.float32(0.1)))

mod2 = isochrones.SingleStarModel(isochrones.get_ichrone("mist", bands=["G", "BP"]),
                                 G=(np.float64(12.5), np.float64(0.1)),
                                 BP=(np.float64(12.5), np.float64(0.1)))

np.random.seed(42)
pars = np.random.rand(mod1.n_params)
mod1.mnest_prior(pars, None, None)
mod1.lnpost(pars) - mod2.lnpost(pars)

[timothydmorton]

To handle the first, we can just cast to float64 on initialization to avoid that. The second, however, doesn't seem to have anything to do with dtypes(!):

import isochrones
import numpy as np

mod1 = isochrones.SingleStarModel(isochrones.get_ichrone("mist", bands=["G", "BP"]),
                                 G=(12.5, 0.1),
                                 BP=(12.5, 0.1))

mod2 = isochrones.SingleStarModel(isochrones.get_ichrone("mist", bands=["G", "BP"]),
                                 G=(12.5, 0.1),
                                 BP=(12.5, 0.1))

mod3 = isochrones.SingleStarModel(isochrones.get_ichrone("mist", bands=["G", "BP"]),
                                 G=(12.5, 0.1),
                                 BP=(12.5, 0.1))

mod4 = isochrones.SingleStarModel(isochrones.get_ichrone("mist", bands=["G", "BP"]),
                                 G=(12.5, 0.1),
                                 BP=(12.5, 0.1))

np.random.seed(42)
pars = np.random.rand(mod1.n_params)
mod1.mnest_prior(pars, None, None)
print(mod2.lnpost(pars) - mod1.lnpost(pars))
print(mod3.lnpost(pars) - mod2.lnpost(pars))
print(mod4.lnpost(pars) - mod3.lnpost(pars))

gives

-0.04686564768883272
0.0
0.0

In fact, it seems as though something in the call to mnest_prior is changing something internally, as the outlier appears to be whichever object's mnest_prior function is called. eep.

[timothydmorton]

This difference also comes from the prior, not the likelihood. With the different dtypes, the difference in lnlike is on the order of 1e-5.

[timothydmorton]

Fixed on master.