docs: def aggregate_duplicates(X, Y, aggregator="mean", precision=precision):
          """ A function that will attempt to collapse duplicates in
              domain space, X, by aggregating values over the range space,
              Y.
              @ In, X, an m-by-n array of values specifying m
              n-dimensional samples
              @ In, Y, a m vector of values specifying the output
              responses corresponding to the m samples specified by X
              @ In, aggregator, an optional string or callable object that
              specifies what type of aggregation to do when duplicates are
              found in the domain space. Default value is mean meaning the
              code will calculate the mean range value over each of the
              unique, duplicated samples.
              @ In, precision, an optional positive integer specifying how
              many digits numbers should be rounded to in order to
              determine if they are unique or not.
              @ Out, (unique_X, aggregated_Y), a tuple where the first
              value is an m'-by-n array specifying the unique domain
              samples and the second value is an m' vector specifying the
              associated range values. m' <= m.
          """
          if callable(aggregator):
              pass
          elif "min" in aggregator.lower():
              aggregator = np.min
          elif "max" in aggregator.lower():
              aggregator = np.max
          elif "median" in aggregator.lower():
              aggregator = np.median
          elif aggregator.lower() in ["average", "mean"]:
              aggregator = np.mean
          elif "first" in aggregator.lower():
  
              def aggregator(x):
                  return x[0]
  
          elif "last" in aggregator.lower():
  
              def aggregator(x):
                  return x[-1]
  
          else:
              warnings.warn(
                  'Aggregator "{}" not understood. Skipping sample '
                  "aggregation.".format(aggregator)
              )
              return X, Y
  
          is_y_multivariate = Y.ndim > 1
  
          X_rounded = X.round(decimals=precision)
          unique_xs = np.unique(X_rounded, axis=0)
  
          old_size = len(X_rounded)
          new_size = len(unique_xs)
          if old_size == new_size:
              return X, Y
  
          if not is_y_multivariate:
              Y = np.atleast_2d(Y).T
  
          reduced_y = np.empty((new_size, Y.shape[1]))
  
          warnings.warn(
              "Domain space duplicates caused a data reduction. "
              + "Original size: {} vs. New size: {}".format(old_size, new_size)
          )
          for col in range(Y.shape[1]):
              for i, distinct_row in enumerate(unique_xs):
                  filtered_rows = np.all(X_rounded == distinct_row, axis=1)
                  reduced_y[i, col] = aggregator(Y[filtered_rows, col])
  
          if not is_y_multivariate:
              reduced_y = reduced_y.flatten()
  
          return unique_xs, reduced_y