`import gctree` causes `FloatingPointError` in independent matplotlib/numpy calls

[wsdewitt]

Certain matplotlib calls—seemingly specific to computing log-scales axes—fail after importing gctree, even without any calls to the gctree package.

The following MWE has been demonstrated on ARM Mac, Intel Mac, and x86_64 Linux.

Plot a log-scaled color bar before gctree import:

import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import LogNorm

plt.colorbar(ScalarMappable(norm=LogNorm(vmin=1, vmax=10)))

Plot a log-scaled color bar after gctree import:

import gctree

Same color bar command as above:

plt.colorbar(ScalarMappable(norm=LogNorm(vmin=1, vmax=10)))

---------------------------------------------------------------------------
FloatingPointError                        Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 plt.colorbar(ScalarMappable(norm=LogNorm(vmin=1, vmax=10)))

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/matplotlib/pyplot.py:2109, in colorbar(mappable, cax, ax, **kw)
   2104     if mappable is None:
   2105         raise RuntimeError('No mappable was found to use for colorbar '
   2106                            'creation. First define a mappable such as '
   2107                            'an image (with imshow) or a contour set ('
   2108                            'with contourf).')
-> 2109 ret = gcf().colorbar(mappable, cax=cax, ax=ax, **kw)
   2110 return ret

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/matplotlib/figure.py:1210, in FigureBase.colorbar(self, mappable, cax, ax, use_gridspec, **kw)
   1206 NON_COLORBAR_KEYS = ['fraction', 'pad', 'shrink', 'aspect', 'anchor',
   1207                      'panchor']
   1208 cb_kw = {k: v for k, v in kw.items() if k not in NON_COLORBAR_KEYS}
-> 1210 cb = cbar.Colorbar(cax, mappable, **cb_kw)
   1212 if not userax:
   1213     self.sca(current_ax)

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/matplotlib/colorbar.py:494, in Colorbar.__init__(self, ax, mappable, cmap, norm, alpha, values, boundaries, orientation, ticklocation, extend, spacing, ticks, format, drawedges, filled, extendfrac, extendrect, label)
    492 else:
    493     self.formatter = format  # Assume it is a Formatter or None
--> 494 self.draw_all()
    496 if isinstance(mappable, contour.ContourSet) and not mappable.filled:
    497     self.add_lines(mappable)

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/matplotlib/colorbar.py:583, in Colorbar.draw_all(self)
    581 self.vmin, self.vmax = self._boundaries[self._inside][[0, -1]]
    582 # Compute the X/Y mesh.
--> 583 X, Y = self._mesh()
    584 # draw the extend triangles, and shrink the inner axes to accommodate.
    585 # also adds the outline path to self.outline spine:
    586 self._do_extends()

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/matplotlib/colorbar.py:1165, in Colorbar._mesh(self)
   1163 norm.vmin = self.vmin
   1164 norm.vmax = self.vmax
-> 1165 y, _ = self._proportional_y()
   1166 # Use the vmin and vmax of the colorbar, which may not be the same
   1167 # as the norm. There are situations where the colormap has a
   1168 # narrower range than the colorbar and we want to accommodate the
   1169 # extra contours.
   1170 if (isinstance(norm, (colors.BoundaryNorm, colors.NoNorm))
   1171         or self.boundaries is not None):
   1172     # not using a norm.

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/matplotlib/colorbar.py:1280, in Colorbar._proportional_y(self)
   1278         yscaled = y
   1279 else:
-> 1280     y = self.norm(self._boundaries.copy())
   1281     y = np.ma.filled(y, np.nan)
   1282     # the norm and the scale should be the same...

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/matplotlib/colors.py:1543, in make_norm_from_scale.<locals>.Norm.__call__(self, value, clip)
   1541     raise ValueError("Invalid vmin or vmax")
   1542 t_value -= t_vmin
-> 1543 t_value /= (t_vmax - t_vmin)
   1544 t_value = np.ma.masked_invalid(t_value, copy=False)
   1545 return t_value[0] if is_scalar else t_value

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/numpy/ma/core.py:4331, in MaskedArray.__itruediv__(self, other)
   4326 """
   4327 True divide self by other in-place.
   4328 
   4329 """
   4330 other_data = getdata(other)
-> 4331 dom_mask = _DomainSafeDivide().__call__(self._data, other_data)
   4332 other_mask = getmask(other)
   4333 new_mask = mask_or(other_mask, dom_mask)

File ~/miniconda3/envs/replay/lib/python3.9/site-packages/numpy/ma/core.py:851, in _DomainSafeDivide.__call__(self, a, b)
    849 a, b = np.asarray(a), np.asarray(b)
    850 with np.errstate(invalid='ignore'):
--> 851     return umath.absolute(a) * self.tolerance >= umath.absolute(b)

FloatingPointError: underflow encountered in multiply

Conda environment

name: test
channels:
  - conda-forge
  - defaults
dependencies:
  - python=3.9
  - gctree

[wsdewitt]

h/t to @willdumm, who guessed this is probably due to the global change in numpy numerical error handling that occurs at module level in gctree.branching_processes: https://github.com/matsengrp/gctree/blob/28390fb87512b0535e4ea04399e78e8a10f22899/gctree/branching_processes.py#L38 Indeed, if I replace import gctree with np.seterr(all="raise") I get the same FloatingPointError.

We should remove this global np.seterr, and instead wrap numerical code in the np.errstate context manager to avoid these pesky side effects.

[willdumm]

That’s good to know about! It looks like a function wrapper is also provided that does the same thing, if that seems like it makes sense. I’ll think carefully about which code needs errors turned on.