def aggregate_duplicates(X, Y, aggregator="mean", precision=precision):
        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