///
/// Element-wise functor that performs the operation:
/// ```
/// scale[i] * (x[i] + shift[i])
/// ```
/// For each component `i` in a vector dataset.
///
template <typename T> class ScaleShift {
  public:
    using vector_type = std::vector<double>;
    using ptr_type = std::shared_ptr<vector_type>;

    ScaleShift(const ptr_type& scale, const ptr_type& shift)
        : scale_{scale}
        , shift_{shift}
        , modified_buffer_(scale->size()) {
        if ((*scale_).size() != (*shift_).size()) {
            throw ANNEXCEPTION("Scale and shift mismatch!");
        }
    }

    ScaleShift(const vector_type& scale, const vector_type& shift)
        : ScaleShift{
              std::make_shared<vector_type>(scale), std::make_shared<vector_type>(shift)} {}

    ///
    /// Return the vector dimensionality that this operator is constructed to operate on.
    ///
    size_t size() const { return scale_->size(); }

    ///
    /// Subtract the stored bias from the data span and return a data-span to the modified
    /// data. Does not modify its argument.
    ///
    /// The returned span is valid so long as:
    /// * The parent `RemoveBias` is live.
    /// * `operator()` is not called on another data.
    ///
    /// Pre-conditions:
    /// * `data.size() == size()`
    ///
    std::span<const T> operator()(const std::span<const T>& data) {
        assert(data.size() == size());
        assert(shift_->size() == scale_->size());
        assert(modified_buffer_.size() == size());

        const auto& shift = *shift_;
        const auto& scale = *scale_;
        for (size_t i = 0, imax = size(); i < imax; ++i) {
            modified_buffer_[i] =
                static_cast<T>(scale[i] * (static_cast<double>(data[i]) + shift[i]));
        }
        return lib::as_const_span(modified_buffer_);
    }

  private:
    // Generally, don't modify `bias_` because it's meant to be shared among thread-local
    // copies.
    std::shared_ptr<std::vector<double>> scale_;
    std::shared_ptr<std::vector<double>> shift_;
    std::vector<T> modified_buffer_;
};

// Determine the average value for each component.
// Remove this bias from each component and return a distance functor that is able to
// lazily re-apply the bias.
struct VectorBias : public DatasetPreOpBase {
    static std::string name() { return "preop-vector-bias"; }

    template <typename T> using element_type_t = typename T::element_type;
    using misc_type = std::vector<double>;

    ///
    /// Compute the mean of each dimension in the dataset.
    /// Return a tuple consisting of:
    /// (1) A distance function that preserves the semantics of the original distance
    ///     function but can be applied to a dataset that has had the per-dimension mean
    ///     removed from each element.
    /// (2) A copyable map operator that operates subtracts the componentwise mean from
    ///     an entry in the dataset.
    ///
    template <data::ImmutableMemoryDataset Data, svs::threads::ThreadPool Pool>
    std::tuple<ScaleShift<element_type_t<Data>>, misc_type>
    operator()(const Data& data, Pool& pool) const {
        using T = element_type_t<Data>;

        // Compute the component-wise mean of the dataset.
        // Negate the medioid to get the bias we've applied to the dataset.
        std::vector<double> means_f64 = utils::compute_medioid(data, pool);
        auto negative_means = std::vector<double>(means_f64.begin(), means_f64.end());
        range::negate(negative_means);

        auto ones = std::vector<double>(negative_means.size(), 1.0);
        return std::make_tuple(
            ScaleShift<T>(std::move(ones), std::move(negative_means)), std::move(means_f64)
        );
    }
};

auto ones = std::vector<double>(negative_means.size(), 1.0);

template <size_t Bits, size_t Extent> class MinRange : public CVStorage {
  public:
    using encoding_type = Encoding<Unsigned, Bits>;
    using return_type = ScaledBiasedVector<Bits, Extent>;
    static constexpr float min_s = encoding_type::min();
    static constexpr float max_s = encoding_type::max();

    ///
    /// Construct a MinRange encoder that uses per-vector constants for its scaling
    /// parameters.
    ///
    MinRange(lib::MaybeStatic<Extent> size = {})
        : CVStorage{}
        , size_{size} {}

    ///
    /// Construct a MinRange encoder that uses pre-determined minimum and maximum values
    /// for its scaling parameters.
    ///
    MinRange(float min, float max, lib::MaybeStatic<Extent> size = {})
        : CVStorage{}
        , discover_extrema_{false}
        , min_{min}
        , max_{max}
        , size_{size} {}

    size_t size() const { return size_; }

    std::pair<float, float> extrema(lib::AnySpanLike auto data) {
        if (!discover_extrema_) {
            return std::make_pair(min_, max_);
        }

        float min{std::numeric_limits<float>::max()};
        float max{std::numeric_limits<float>::min()};
        for (size_t i = 0; i < data.size(); ++i) {
            float val = static_cast<float>(data[i]);
            min = std::min(min, val);
            max = std::max(max, val);
        }
        return std::make_pair(min, max);
    }

    // Compression Operator.
    return_type operator()(lib::AnySpanLike auto data, selector_t selector = 0) {
        // Find the maximum absolute value of the data to encode.
        auto [min, max] = extrema(data);

        float bias = min;
        float decompressor = 1;
        float compressor = 1;
        if (max != min) {
            decompressor = (max - min) / max_s;
            compressor = max_s / (max - min);
        }

        auto cv = view<Unsigned, Bits, Extent>(size_);
        for (size_t i = 0; i < data.size(); ++i) {
            auto v = lib::narrow<float>(data[i]);
            auto compressed = crunch(compressor, v - bias, min_s, max_s);
            cv.set(compressed, i);
        }

        return return_type{decompressor, bias, selector, cv};
    }

  private:
    bool discover_extrema_ = true;
    float min_ = 0;
    float max_ = 0;
    [[no_unique_address]] lib::MaybeStatic<Extent> size_ = {};
};

auto [min, max] = extrema(data);

size_t insert_inner(value_type neighbor) {
    const auto start = begin();
    // Binary search to the first location where `distance` is less than the stored
    // neighbor.
    auto pos = std::lower_bound(
        start,
        end(),
        neighbor.distance(),
        [&](const value_type& other, const float& d) {
            return !compare_(d, other.distance());
        }
    );

    // Because repeat ids can exist, we have to search until we're sure no repeat will
    // be found.
    //
    // To do that, we start one before the insertion position.
    // Because each instance of repeat ids should have the same distance, we only need
    // to look at `ids` until the buffer elements have a distance less than the
    // distance of the current node.
    if (pos != start) {
        auto back = pos;
        do {
            --back;
            const auto& candidate = *back;
            if (compare_(candidate.distance(), neighbor.distance())) {
                break;
            } else if (candidate.id_ == neighbor.id()) {
                return size() + 1;
            }
        } while (back != start);
    }

    // We're explicitly avoiding moving the underlying range in this implementaion,
    // so we don't need to worry about iterator invalidation.
    //
    // Nevertheless, probably safe to compute the actual index because messing around
    // with the internal state of the buffer.
    size_t i = pos - start;
    unsafe_insert(neighbor, pos);
    size_ = std::min(size_ + 1, capacity_);
    best_unvisited_ = std::min(best_unvisited_, i);
    return i;
}

// We're explicitly avoiding moving the underlying range in this implementaion,
// so we don't need to worry about iterator invalidation.
//
// Nevertheless, probably safe to compute the actual index because messing around
// with the internal state of the buffer.
size_t i = pos - start;
unsafe_insert(neighbor, pos);
size_ = std::min(size_ + 1, capacity_);
best_unvisited_ = std::min(best_unvisited_, i);
return i;