Interacting with TreeMesh object before finalization causes __repr__ and _repr_html_ to later error

[williamjsdavis]

This is some odd behavior I encountered when interacting with a TreeMesh object before finalizing it.

Working example

The code here is modified from the documentation: link

from discretize import TreeMesh
from discretize.utils import mkvc
import numpy as np

dx = 5  # minimum cell width (base mesh cell width) in x
dy = 5  # minimum cell width (base mesh cell width) in y

x_length = 300.0  # domain width in x
y_length = 300.0  # domain width in y

# Compute number of base mesh cells required in x and y
nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))
nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))

# Define the base mesh
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
mesh = TreeMesh([hx, hy], x0="CC")

# Note this print statement
# print(mesh.h_gridded)

# Refine surface topography
xx = mesh.nodes_x
yy = -3 * np.exp((xx**2) / 100**2) + 50.0
pts = np.c_[mkvc(xx), mkvc(yy)]
padding = [[0, 2], [0, 2]]
mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False)

# Refine mesh near points
xx = np.array([0.0, 10.0, 0.0, -10.0])
yy = np.array([-20.0, -10.0, 0.0, -10])
pts = np.c_[mkvc(xx), mkvc(yy)]
mesh.refine_points(pts, padding_cells_by_level=[2, 2], finalize=False)

mesh.finalize()

mesh.__repr__

The final display of mesh gives the following output:

<bound method TreeMesh.__repr__ of 
QuadTreeMesh: 14.01% filled

Level : Number of cells               Mesh Extent               Cell Widths    
-----------------------           min     ,     max            min   ,   max   
  2   :        2             ---------------------------   --------------------
  3   :       28          x:    -160.0    ,    160.0          5.0    ,    80.0   
  4   :       50          y:    -160.0    ,    160.0          5.0    ,    80.0   
  5   :       166      
  6   :       328      
-----------------------
Total :       574      >

Erroring example

Running the exact same code, with the print statement not commented out, i.e.,

from discretize import TreeMesh
from discretize.utils import mkvc
import numpy as np

dx = 5  # minimum cell width (base mesh cell width) in x
dy = 5  # minimum cell width (base mesh cell width) in y

x_length = 300.0  # domain width in x
y_length = 300.0  # domain width in y

# Compute number of base mesh cells required in x and y
nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))
nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))

# Define the base mesh
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
mesh = TreeMesh([hx, hy], x0="CC")

# Note this print statement
print(mesh.h_gridded)

# Refine surface topography
xx = mesh.nodes_x
yy = -3 * np.exp((xx**2) / 100**2) + 50.0
pts = np.c_[mkvc(xx), mkvc(yy)]
padding = [[0, 2], [0, 2]]
mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False)

# Refine mesh near points
xx = np.array([0.0, 10.0, 0.0, -10.0])
yy = np.array([-20.0, -10.0, 0.0, -10])
pts = np.c_[mkvc(xx), mkvc(yy)]
mesh.refine_points(pts, padding_cells_by_level=[2, 2], finalize=False)

mesh.finalize()

mesh.__repr__

Causes the following error:

ValueError: zero-size array to reduction operation minimum which has no identity

Full stacktrace:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/Documents/Projects/vermelhos-pomdp/.venv/lib/python3.12/site-packages/IPython/core/formatters.py:770, in PlainTextFormatter.__call__(self, obj)
    763 stream = StringIO()
    764 printer = pretty.RepresentationPrinter(stream, self.verbose,
    765     self.max_width, self.newline,
    766     max_seq_length=self.max_seq_length,
    767     singleton_pprinters=self.singleton_printers,
    768     type_pprinters=self.type_printers,
    769     deferred_pprinters=self.deferred_printers)
--> 770 printer.pretty(obj)
    771 printer.flush()
    772 return stream.getvalue()

File ~/Documents/Projects/vermelhos-pomdp/.venv/lib/python3.12/site-packages/IPython/lib/pretty.py:394, in RepresentationPrinter.pretty(self, obj)
    391 for cls in _get_mro(obj_class):
    392     if cls in self.type_pprinters:
    393         # printer registered in self.type_pprinters
--> 394         return self.type_pprinters[cls](obj, self, cycle)
    395     else:
    396         # deferred printer
    397         printer = self._in_deferred_types(cls)

File ~/Documents/Projects/vermelhos-pomdp/.venv/lib/python3.12/site-packages/IPython/lib/pretty.py:794, in _repr_pprint(obj, p, cycle)
    792 """A pprint that just redirects to the normal repr function."""
    793 # Find newlines and replace them with p.break_()
--> 794 output = repr(obj)
    795 lines = output.splitlines()
    796 with p.group():

File ~/Documents/Projects/vermelhos-pomdp/.venv/lib/python3.12/site-packages/discretize/tree_mesh.py:302, in TreeMesh.__repr__(self)
    300 h_display.append("-" * (len(h_display[0])))
    301 h_gridded = self.h_gridded
--> 302 mins = np.min(h_gridded, axis=0)
    303 maxs = np.max(h_gridded, axis=0)
    304 for dim in range(self.dim):

File ~/Documents/Projects/vermelhos-pomdp/.venv/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:3302, in min(a, axis, out, keepdims, initial, where)
   3190 @array_function_dispatch(_min_dispatcher)
   3191 def min(a, axis=None, out=None, keepdims=np._NoValue, initial=np._NoValue,
   3192         where=np._NoValue):
   3193     """
   3194     Return the minimum of an array or minimum along an axis.
   3195 
   (...)
   3300     6
   3301     """
-> 3302     return _wrapreduction(a, np.minimum, 'min', axis, None, out,
   3303                           keepdims=keepdims, initial=initial, where=where)

File ~/Documents/Projects/vermelhos-pomdp/.venv/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:86, in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs)
     83         else:
     84             return reduction(axis=axis, out=out, **passkwargs)
---> 86 return ufunc.reduce(obj, axis, dtype, out, **passkwargs)

ValueError: zero-size array to reduction operation minimum which has no identity

Additional comments

The same behavior happens with the _repr_html_ method of TreeMesh.

[santisoler]

Hi @williamjsdavis, thanks for opening this Issue.

I couldn't reproduce the error with the scripts you provided, but I suspect what issue you are trying to report. Trying to print a TreeMesh before it gets finalized raises an error. For example:

from discretize import TreeMesh
from discretize.utils import mkvc
import numpy as np

dx = 5  # minimum cell width (base mesh cell width) in x
dy = 5  # minimum cell width (base mesh cell width) in y

x_length = 300.0  # domain width in x
y_length = 300.0  # domain width in y

# Compute number of base mesh cells required in x and y
nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))
nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))

# Define the base mesh
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
mesh = TreeMesh([hx, hy], x0="CC")

# Print the TreeMesh before we finilize it
print(mesh)

Traceback (most recent call last):
  File "/home/santi/git/simpeg/should-error.py", line 21, in <module>
    print(mesh)
  File "/home/santi/.miniforge3/envs/simpeg-dev/lib/python3.11/site-packages/discretize/tree_mesh.py", line 303, in __repr__
    mins = np.min(h_gridded, axis=0)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/santi/.miniforge3/envs/simpeg-dev/lib/python3.11/site-packages/numpy/_core/fromnumeric.py", line 3042, in min
    return _wrapreduction(a, np.minimum, 'min', axis, None, out,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/santi/.miniforge3/envs/simpeg-dev/lib/python3.11/site-packages/numpy/_core/fromnumeric.py", line 86, in _wrapreduction
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: zero-size array to reduction operation minimum which has no identity

Is this the behavior you intended to report?

If so, I think the __repr__() shouldn't raise an error just because the mesh is not finalized.

[williamjsdavis]

I just double checked running my code and I reproduced my original error and working case. I was running these code blocks in jupyter notebook cells; I don't know if that makes a difference? Is there any additional information I can provide to help you reproduce my error?

I was able to reproduce your example error. It is related to the behavior I was intended to report. The difference is that, in my example, even just printing mesh.h_gridded will cause the __repr__() to error after the mesh is finalized.

[santisoler]

Thanks for the swift reply.

I opened a #392 to fix the error I described: in that PR the __repr__ method will not error out when the mesh is not finalized.

I can reproduce the error when I run your script in a Jupyter Notebook. I think the problem is not coming from the __repr__, but by accessing the h_gridded property before the mesh is finalized. For example, the following script creates the same mesh, access the h_gridded property and then finalize the mesh:

from discretize import TreeMesh
from discretize.utils import mkvc
import numpy as np

dx = 5  # minimum cell width (base mesh cell width) in x
dy = 5  # minimum cell width (base mesh cell width) in y

x_length = 300.0  # domain width in x
y_length = 300.0  # domain width in y

# Compute number of base mesh cells required in x and y
nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))
nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))

# Define the base mesh
hx = [(dx, nbcx)]
hy = [(dy, nbcy)]
mesh = TreeMesh([hx, hy], x0="CC")

# Access h_gridded before finalizing
mesh.h_gridded

# Refine surface topography
xx = mesh.nodes_x
yy = -3 * np.exp((xx**2) / 100**2) + 50.0
pts = np.c_[mkvc(xx), mkvc(yy)]
padding = [[0, 2], [0, 2]]
mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False)

# Refine mesh near points
xx = np.array([0.0, 10.0, 0.0, -10.0])
yy = np.array([-20.0, -10.0, 0.0, -10])
pts = np.c_[mkvc(xx), mkvc(yy)]
mesh.refine_points(pts, padding_cells_by_level=[2, 2], finalize=False)
mesh.finalize()

# h_gridded is empty even after finalizing the mesh
print(f"{mesh.h_gridded=}")

Even after finalizing the mesh, the h_gridded is still empty:

mesh.h_gridded=array([], shape=(0, 2), dtype=float64)

I think the h_gridded should be updated after the mesh gets finalized, or it should not be cached if it hasn't been finalized yet. Once computed, it's cached in a private self._h_gridded attribute:

discretize/discretize/_extensions/tree_ext.pyx

Lines 2058 to 2086 in 60f3af1

	@property
	def h_gridded(self):
	"""Gridded cell dimensions.
	This property returns a numpy array of shape (n_cells, dim)
	containing the dimensions of the cells along each axis
	direction in order. E.g. the columns of *h_gridded* for a
	3D tree mesh would be ordered [hx,hy,hz].
	Returns
	-------
	(n_cells, dim) numpy.ndarray of float
	Gridded cell dimensions
	"""
	if self._h_gridded is not None:
	return self._h_gridded
	cdef np.float64_t[:, :] gridCH
	cdef np.int64_t ii, ind, dim
	cdef np.float64_t len
	cdef int epc = 4 if self._dim==3 else 2
	dim = self._dim
	self._h_gridded = np.empty((self.n_cells, dim), dtype=np.float64)
	gridCH = self._h_gridded
	for cell in self.tree.cells:
	ind = cell.index
	for ii in range(dim):
	gridCH[ind, ii] = cell.edges[ii*epc].length
	return self._h_gridded

)

Maybe we could return an empty array and not cache anything if the mesh is not finalized.

cc @jcapriot

[jcapriot]

Related to #293

We need to add guards on most TreeMesh operations that will immediately error if the TreeMesh is not finalized

[williamjsdavis]

@santisoler Yes this is the erroring behavior that I was experiencing, thanks for opening that PR!