code: def make_cut(lb, ub, RVStart, RVend, vol):
      params, subdiams, diam, probmass = test_cuboids(
          lb, ub, RVstart, RVend, vol
      )
      SOLVED_PARAMETERS, SUB_DIAMETERS = params, subdiams
      TOTAL_DIAMETERS, PROBABILITY_MASS = diam, probmass
  
      #####XXX: probably needs to be some sort of loop ######
      # get region with largest probability mass
      region = probmass.index(max(probmass))
      NEW_SLICES = [region, region + 1]
  
      # get direction with largest subdiameter
      direction = subdiams[region].index(max(subdiams[region]))
  
      # XXX: should check 'params' for 'no more improvement possilble'
      # XXX: ... if so, check all the directions in that region
      # XXX: ... if still so, then select region with next largest probmass
      #####XXX: probably needs to be some sort of loop ######
  
      # get the midpoint
      cutvalue = 0.5 * (ub[region][direction] + lb[region][direction])
  
      # modify bounds to include cut plane
      l = lb[: region + 1]
      l += [
          lb[region][:direction]
          + [cutvalue]
          + lb[region][direction + 1 :]
      ]
      l += lb[region + 1 :]
      u = ub[:region]
      u += [
          ub[region][:direction]
          + [cutvalue]
          + ub[region][direction + 1 :]
      ]
      u += ub[region:]
  
      return l, u
  
file_path: mystic/examples4/TEST_surrogate_cut.py
imports: from TEST_surrogate_diam import *
  from mystic.math.stats import volume, prob_mass, mean, mcdiarmid_bound
  from mystic.math.integrate import integrate
  from mystic.math.samples import sample
label_desc: Use image node_id 4 for calling a global function with example usage: make_cut(lb, ub, RVStart, RVend, vol) and returns: l,  u