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.

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