Better infinite bound handling in Python API plot methods

[jeinstei]

Description

Currently, the plot methods in the Python API rely on setting defaults or raising exceptions in cases where geometry is unbounded. in model.plot() the bbox currently is set by self.geometry.bounding_box, while in Plot.from_geometry(), an exception is raised. If we add an argument to these calls for a [width, height], we can immediately generate plots with user-generated extents for the basis or field instead of relying on a possibly unbounded dimension.

Alternatives

This seems possibly similar to the discussions raised in #2632 where there are inconsistencies in BoundingBox handling between APIs. I'm not sure what the desired alternatives would be besides creating custom overrides. Possibly changing the bounding_box property or adding a bounding_box_restricted that handles infinite extents in a reproducible manner.

Compatibility

This should add new arguments to some plotting APIs, but should also not change existing usage as they would have well-handled default values.

[pshriwise]

@jeinstei I believe model.plot will accept a width argument that is an Iterable of two ints via the **kwargs. These keyword arguments are passed on to the Universe.plot method.

As for Plot.from_geometry() I think that would be a welcome addition!

[jeinstei]

Thanks @pshriwise -- I see the flow down on model.plot -> geometry.plot. I'll give that a try.

[paulromano]

For the Model.plot() method, one can already set the origin and width as arguments, which are taken as keyword arguments that are passed on to Universe.plot. For Plot.from_geometry, the main purpose of that classmethod is to automatically determine the origin/width, so if it can't do that based on the bounding box, raising an exception is indeed the intended behavior.

[jeinstei]

Once I figure out a solid workflow for this, I will close out my (hopefully) unnecessary PR, and maybe instead convert this to a plotting example in the documentation then. Thanks for the assistance.

[jeinstei]

I'm currently kit-bashing with a workaround function to create at least a plot in a programmatic manner. I know I could speed up this code with some fancy array work, but this is a MWE. I'll have a PR with some documentation updates for plots.rst and plots.py for review and debate this week.

def determine_extents_from_model(model):
    """Takes a model with infinte bounds and finds extents automatically

    Parameters
    ----------
    model : openmc.Model
        Model with infinite bounds to determine extent of

    Returns
    -------
    work : openmc.BoundingBox
        Bounding box (lower_left, upper_right) extents of model

    width : (float, float, float)
        Calculated width along each axis
        
    origin : (float, float, float)
        Calculated center of found extents

    Raises
    ------
    ValueError
        If calculated extents is still unbounded

    """
    base = model.geometry.bounding_box
    work = openmc.BoundingBox(base[0], base[1])
    
    # Determine indices where extents are infinite
    base_infinite = np.isinf(base)
    if not base_infinite.any():
        # Calculate width and origin of original extents
        width = base[1] - base[0]
        origin = (width/2) + base[0]
        return work, width, origin
        
    active_tests = np.argwhere(base_infinite)
    
    bboxes = []
    for k,v in model.geometry.get_all_cells().items():
        bboxes.append(np.asarray(v.bounding_box))
    
    # Loop over bounding boxes for indices to test and determine max/min
    # NOTE this could be done using masked array math, etc, but is not needed here
    for ind in active_tests:
        idx = tuple(ind)
        for box in bboxes:
            bval = work[idx]
            tval = box[idx]
    
            # Skip testing if test value is infinite itself
            if np.isinf(tval):
                continue
    
            # If base value is infinite, save initial value
            if np.isinf(bval):
                work[idx] = tval
                continue
            
            if ind[0] == 0: # lower_left
                if tval < bval:
                    work[idx] = tval
            else: # upper_right
                if tval > bval:
                    work[idx] = tval

    if np.isinf(work).any():
        raise ValueError("Unable to determine bounded extents from model cells")
    
    # Calculate width and origin of new extents
    width = work[1] - work[0]
    origin = (width/2) + work[0]
    
    return work, width, origin