sparse_vector.h

// Copyright 2010-2022 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// Classes to represent sparse vectors.
//
// The following are very good references for terminology, data structures,
// and algorithms:
//
// I.S. Duff, A.M. Erisman and J.K. Reid, "Direct Methods for Sparse Matrices",
// Clarendon, Oxford, UK, 1987, ISBN 0-19-853421-3,
// http://www.amazon.com/dp/0198534213.
//
//
// T.A. Davis, "Direct methods for Sparse Linear Systems", SIAM, Philadelphia,
// 2006, ISBN-13: 978-0-898716-13, http://www.amazon.com/dp/0898716136.
//
//
// Both books also contain a wealth of references.

#ifndef OR_TOOLS_LP_DATA_SPARSE_VECTOR_H_
#define OR_TOOLS_LP_DATA_SPARSE_VECTOR_H_

#include <algorithm>
#include <cstring>
#include <memory>
#include <string>

#include "absl/strings/str_format.h"
#include "ortools/base/integral_types.h"
#include "ortools/base/logging.h"  // for CHECK*
#include "ortools/graph/iterators.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/lp_data/permutation.h"
#include "ortools/util/return_macros.h"

namespace operations_research {
namespace glop {

template <typename IndexType>
class SparseVectorEntry;

// --------------------------------------------------------
// SparseVector
// --------------------------------------------------------
// This class allows to store a vector taking advantage of its sparsity.
// Space complexity is in O(num_entries).
// In the current implementation, entries are stored in a first-in order (order
// of SetCoefficient() calls) when they are added; then the "cleaning" process
// sorts them by index (and duplicates are removed: the last entry takes
// precedence).
// Many methods assume that the entries are sorted by index and without
// duplicates, and DCHECK() that.
//
// Default copy construction is fully supported.
//
// This class uses strong integer types (i.e. no implicit cast to/from other
// integer types) for both:
// - the index of entries (eg. SparseVector<RowIndex> is a SparseColumn,
//   see ./sparse_column.h).
// - the *internal* indices of entries in the internal storage, which is an
//   entirely different type: EntryType.
// This class can be extended with a custom iterator/entry type for the
// iterator-based API. This can be used to extend the interface with additional
// methods for the entries returned by the iterators; for an example of such
// extension, see SparseColumnEntry in sparse_column.h. The custom entries and
// iterators should be derived from SparseVectorEntry and SparseVectorIterator,
// or at least provide the same public and protected interface.
//
// TODO(user): un-expose this type to client; by getting rid of the
// index-based APIs and leveraging iterator-based APIs; if possible.
template <typename IndexType,
          typename IteratorType = VectorIterator<SparseVectorEntry<IndexType>>>
class SparseVector {
 public:
  typedef IndexType Index;

  typedef StrictITIVector<Index, Fractional> DenseVector;
  typedef Permutation<Index> IndexPermutation;

  using Iterator = IteratorType;
  using Entry = typename Iterator::Entry;

  SparseVector();

  // NOTE(user): STL uses the expensive copy constructor when relocating
  // elements of a vector, unless the move constructor exists *and* it is marked
  // as noexcept. However, the noexcept annotation is banned by the style guide,
  // and the only way to get it is by using the default move constructor and
  // assignment operator generated by the compiler.
  SparseVector(const SparseVector& other);
#if !defined(_MSC_VER)
  SparseVector(SparseVector&& other) = default;
#endif

  SparseVector& operator=(const SparseVector& other);
#if !defined(_MSC_VER)
  SparseVector& operator=(SparseVector&& other) = default;
#endif

  // Read-only API for a given SparseVector entry. The typical way for a
  // client to use this is to use the natural range iteration defined by the
  // Iterator class below:
  //   SparseVector<int> v;
  //   ...
  //   for (const SparseVector<int>::Entry e : v) {
  //     LOG(INFO) << "Index: " << e.index() << ", Coeff: " << e.coefficient();
  //   }
  //
  // Note that this can only be used when the vector has no duplicates.
  //
  // Note(user): using either "const SparseVector<int>::Entry&" or
  // "const SparseVector<int>::Entry" yields the exact same performance on the
  // netlib, thus we recommend to use the latter version, for consistency.
  Iterator begin() const;
  Iterator end() const;

  // Clears the vector, i.e. removes all entries.
  void Clear();

  // Clears the vector and releases the memory it uses.
  void ClearAndRelease();

  // Reserve the underlying storage for the given number of entries.
  void Reserve(EntryIndex new_capacity);

  // Returns true if the vector is empty.
  bool IsEmpty() const;

  // Cleans the vector, i.e. removes zero-values entries, removes duplicates
  // entries and sorts remaining entries in increasing index order.
  // Runs in O(num_entries * log(num_entries)).
  void CleanUp();

  // Returns true if the entries of this SparseVector are in strictly increasing
  // index order and if the vector contains no duplicates nor zero coefficients.
  // Runs in O(num_entries). It is not const because it modifies
  // possibly_contains_duplicates_.
  bool IsCleanedUp() const;

  // Swaps the content of this sparse vector with the one passed as argument.
  // Works in O(1).
  void Swap(SparseVector* other);

  // Populates the current vector from sparse_vector.
  // Runs in O(num_entries).
  void PopulateFromSparseVector(const SparseVector& sparse_vector);

  // Populates the current vector from dense_vector.
  // Runs in O(num_indices_in_dense_vector).
  void PopulateFromDenseVector(const DenseVector& dense_vector);

  // Appends all entries from sparse_vector to the current vector; the indices
  // of the appended entries are increased by offset. If the current vector
  // already has a value at an index changed by this method, this value is
  // overwritten with the value from sparse_vector.
  // Note that while offset may be negative itself, the indices of all entries
  // after applying the offset must be non-negative.
  void AppendEntriesWithOffset(const SparseVector& sparse_vector, Index offset);

  // Returns true when the vector contains no duplicates. Runs in
  // O(max_index + num_entries), max_index being the largest index in entry.
  // This method allocates (and deletes) a Boolean array of size max_index.
  // Note that we use a mutable Boolean to make subsequent call runs in O(1).
  bool CheckNoDuplicates() const;

  // Same as CheckNoDuplicates() except it uses a reusable boolean vector
  // to make the code more efficient. Runs in O(num_entries).
  // Note that boolean_vector should be initialized to false before calling this
  // method; It will remain equal to false after calls to CheckNoDuplicates().
  // Note that we use a mutable Boolean to make subsequent call runs in O(1).
  bool CheckNoDuplicates(StrictITIVector<Index, bool>* boolean_vector) const;

  // Defines the coefficient at index, i.e. vector[index] = value;
  void SetCoefficient(Index index, Fractional value);

  // Removes an entry from the vector if present. The order of the other entries
  // is preserved. Runs in O(num_entries).
  void DeleteEntry(Index index);

  // Sets to 0.0 (i.e. remove) all entries whose fabs() is lower or equal to
  // the given threshold.
  void RemoveNearZeroEntries(Fractional threshold);

  // Same as RemoveNearZeroEntries, but the entry magnitude of each row is
  // multiplied by weights[row] before being compared with threshold.
  void RemoveNearZeroEntriesWithWeights(Fractional threshold,
                                        const DenseVector& weights);

  // Moves the entry with given Index to the first position in the vector. If
  // the entry is not present, nothing happens.
  void MoveEntryToFirstPosition(Index index);

  // Moves the entry with given Index to the last position in the vector. If
  // the entry is not present, nothing happens.
  void MoveEntryToLastPosition(Index index);

  // Multiplies all entries by factor.
  // i.e. entry.coefficient *= factor.
  void MultiplyByConstant(Fractional factor);

  // Multiplies all entries by its corresponding factor,
  // i.e. entry.coefficient *= factors[entry.index].
  void ComponentWiseMultiply(const DenseVector& factors);

  // Divides all entries by factor.
  // i.e. entry.coefficient /= factor.
  void DivideByConstant(Fractional factor);

  // Divides all entries by its corresponding factor,
  // i.e. entry.coefficient /= factors[entry.index].
  void ComponentWiseDivide(const DenseVector& factors);

  // Populates a dense vector from the sparse vector.
  // Runs in O(num_indices) as the dense vector values have to be reset to 0.0.
  void CopyToDenseVector(Index num_indices, DenseVector* dense_vector) const;

  // Populates a dense vector from the permuted sparse vector.
  // Runs in O(num_indices) as the dense vector values have to be reset to 0.0.
  void PermutedCopyToDenseVector(const IndexPermutation& index_perm,
                                 Index num_indices,
                                 DenseVector* dense_vector) const;

  // Performs the operation dense_vector += multiplier * this.
  // This is known as multiply-accumulate or (fused) multiply-add.
  void AddMultipleToDenseVector(Fractional multiplier,
                                DenseVector* dense_vector) const;

  // WARNING: BOTH vectors (the current and the destination) MUST be "clean",
  // i.e. sorted and without duplicates.
  // Performs the operation accumulator_vector += multiplier * this, removing
  // a given index which must be in both vectors, and pruning new entries whose
  // absolute value are under the given drop_tolerance.
  void AddMultipleToSparseVectorAndDeleteCommonIndex(
      Fractional multiplier, Index removed_common_index,
      Fractional drop_tolerance, SparseVector* accumulator_vector) const;

  // Same as AddMultipleToSparseVectorAndDeleteCommonIndex() but instead of
  // deleting the common index, leave it unchanged.
  void AddMultipleToSparseVectorAndIgnoreCommonIndex(
      Fractional multiplier, Index removed_common_index,
      Fractional drop_tolerance, SparseVector* accumulator_vector) const;

  // Applies the index permutation to all entries: index = index_perm[index];
  void ApplyIndexPermutation(const IndexPermutation& index_perm);

  // Same as ApplyIndexPermutation but deletes the index if index_perm[index]
  // is negative.
  void ApplyPartialIndexPermutation(const IndexPermutation& index_perm);

  // Removes the entries for which index_perm[index] is non-negative and appends
  // them to output. Note that the index of the entries are NOT permuted.
  void MoveTaggedEntriesTo(const IndexPermutation& index_perm,
                           SparseVector* output);

  // Returns the coefficient at position index.
  // Call with care: runs in O(number-of-entries) as entries may not be sorted.
  Fractional LookUpCoefficient(Index index) const;

  // Note this method can only be used when the vector has no duplicates.
  EntryIndex num_entries() const {
    DCHECK(CheckNoDuplicates());
    return EntryIndex(num_entries_);
  }

  // Returns the first entry's index and coefficient; note that 'first' doesn't
  // mean 'entry with the smallest index'.
  // Runs in O(1).
  // Note this method can only be used when the vector has no duplicates.
  Index GetFirstIndex() const {
    DCHECK(CheckNoDuplicates());
    return GetIndex(EntryIndex(0));
  }
  Fractional GetFirstCoefficient() const {
    DCHECK(CheckNoDuplicates());
    return GetCoefficient(EntryIndex(0));
  }

  // Like GetFirst*, but for the last entry.
  Index GetLastIndex() const {
    DCHECK(CheckNoDuplicates());
    return GetIndex(num_entries() - 1);
  }
  Fractional GetLastCoefficient() const {
    DCHECK(CheckNoDuplicates());
    return GetCoefficient(num_entries() - 1);
  }

  // Allows to loop over the entry indices like this:
  //   for (const EntryIndex i : sparse_vector.AllEntryIndices()) { ... }
  // TODO(user): consider removing this, in favor of the natural range
  // iteration.
  ::util::IntegerRange<EntryIndex> AllEntryIndices() const {
    return ::util::IntegerRange<EntryIndex>(EntryIndex(0), num_entries_);
  }

  // Returns true if this vector is exactly equal to the given one, i.e. all its
  // index indices and coefficients appear in the same order and are equal.
  bool IsEqualTo(const SparseVector& other) const;

  // An exhaustive, pretty-printed listing of the entries, in their
  // internal order. a.DebugString() == b.DebugString() iff a.IsEqualTo(b).
  std::string DebugString() const;

 protected:
  // Adds a new entry to the sparse vector, growing the internal buffer if
  // needed. It does not set may_contain_duplicates_ to true.
  void AddEntry(Index index, Fractional value) {
    DCHECK_GE(index, 0);
    // Grow the internal storage if there is no space left for the new entry. We
    // increase the size to max(4, 1.5*current capacity).
    if (num_entries_ == capacity_) {
      // Reserve(capacity_ == 0 ? EntryIndex(4)
      //                        : EntryIndex(2 * capacity_.value()));
      Reserve(capacity_ == 0 ? EntryIndex(4)
                             : EntryIndex(2 * capacity_.value()));
      DCHECK_LT(num_entries_, capacity_);
    }
    const EntryIndex new_entry_index = num_entries_;
    ++num_entries_;
    MutableIndex(new_entry_index) = index;
    MutableCoefficient(new_entry_index) = value;
  }

  // Resizes the sparse vector to a smaller size, without re-allocating the
  // internal storage.
  void ResizeDown(EntryIndex new_size) {
    DCHECK_GE(new_size, 0);
    DCHECK_LE(new_size, num_entries_);
    num_entries_ = new_size;
  }

  // Read-only access to the indices and coefficients of the entries of the
  // sparse vector.
  Index GetIndex(EntryIndex i) const {
    DCHECK_GE(i, 0);
    DCHECK_LT(i, num_entries_);
    return index_[i.value()];
  }
  Fractional GetCoefficient(EntryIndex i) const {
    DCHECK_GE(i, 0);
    DCHECK_LT(i, num_entries_);
    return coefficient_[i.value()];
  }

  // Mutable access to the indices and coefficients of the entries of the sparse
  // vector.
  Index& MutableIndex(EntryIndex i) {
    DCHECK_GE(i, 0);
    DCHECK_LT(i, num_entries_);
    return index_[i.value()];
  }
  Fractional& MutableCoefficient(EntryIndex i) {
    DCHECK_GE(i, 0);
    DCHECK_LT(i, num_entries_);
    return coefficient_[i.value()];
  }

  // The internal storage of the sparse vector. Both the indices and the
  // coefficients are stored in the same buffer; the first
  // sizeof(Index)*capacity_ bytes are used for storing the indices, the
  // following sizeof(Fractional)*capacity_ bytes contain the values. This
  // representation ensures that for small vectors, both the indices and the
  // coefficients are in the same page/cache line.
  // We use a single buffer for both arrays. The amount of data copied during
  // relocations is the same in both cases, and it is much smaller than the cost
  // of an additional allocation - especially when the vectors are small.
  // Moreover, using two separate vectors/buffers would mean that even small
  // vectors would be spread across at least two different cache lines.
  std::unique_ptr<char[]> buffer_;
  EntryIndex num_entries_;
  EntryIndex capacity_;

  // Pointers to the first elements of the index and coefficient arrays.
  Index* index_;
  Fractional* coefficient_;

  // This is here to speed up the CheckNoDuplicates() methods and is mutable
  // so we can perform checks on const argument.
  mutable bool may_contain_duplicates_;

 private:
  // Actual implementation of AddMultipleToSparseVectorAndDeleteCommonIndex()
  // and AddMultipleToSparseVectorAndIgnoreCommonIndex() which is shared.
  void AddMultipleToSparseVectorInternal(
      bool delete_common_index, Fractional multiplier, Index common_index,
      Fractional drop_tolerance, SparseVector* accumulator_vector) const;
};

// --------------------------------------------------------
// SparseVectorEntry
// --------------------------------------------------------

// A reference-like class that points to a certain element of a sparse data
// structure that stores its elements in two parallel arrays. The main purpose
// of the entry class is to support implementation of iterator objects over the
// sparse data structure.
// Note that the entry object does not own the data, and it is valid only as
// long as the underlying sparse data structure; it may also be invalidated if
// the underlying sparse data structure is modified.
template <typename IndexType>
class SparseVectorEntry {
 public:
  using Index = IndexType;

  Index index() const { return index_[i_.value()]; }
  Fractional coefficient() const { return coefficient_[i_.value()]; }

 protected:
  // Creates the sparse vector entry from the given base pointers and the index.
  // We accept the low-level data structures rather than a SparseVector
  // reference to make it possible to use the SparseVectorEntry and
  // SparseVectorIterator classes also for other data structures using the same
  // internal data representation.
  // Note that the constructor is intentionally made protected, so that the
  // entry can be created only as a part of the construction of an iterator over
  // a sparse data structure.
  SparseVectorEntry(const Index* indices, const Fractional* coefficients,
                    EntryIndex i)
      : i_(i), index_(indices), coefficient_(coefficients) {}

  // The index of the sparse vector entry represented by this object.
  EntryIndex i_;
  // The index and coefficient arrays of the sparse vector.
  // NOTE(user): Keeping directly the index and the base pointers gives the
  // best performance with a tiny margin of the options:
  // 1. keep the base pointers and an index of the current entry,
  // 2. keep pointers to the current index and the current coefficient and
  //    increment both when moving the iterator.
  // 3. keep a pointer to the sparse vector object and the index of the current
  //    entry.
  const Index* index_;
  const Fractional* coefficient_;
};

template <typename IndexType, typename IteratorType>
IteratorType SparseVector<IndexType, IteratorType>::begin() const {
  return Iterator(this->index_, this->coefficient_, EntryIndex(0));
}

template <typename IndexType, typename IteratorType>
IteratorType SparseVector<IndexType, IteratorType>::end() const {
  return Iterator(this->index_, this->coefficient_, num_entries_);
}

// --------------------------------------------------------
// SparseVector implementation
// --------------------------------------------------------
template <typename IndexType, typename IteratorType>
SparseVector<IndexType, IteratorType>::SparseVector()
    : num_entries_(0),
      capacity_(0),
      index_(nullptr),
      coefficient_(nullptr),
      may_contain_duplicates_(false) {}

template <typename IndexType, typename IteratorType>
SparseVector<IndexType, IteratorType>::SparseVector(const SparseVector& other) {
  PopulateFromSparseVector(other);
}

template <typename IndexType, typename IteratorType>
SparseVector<IndexType, IteratorType>&
SparseVector<IndexType, IteratorType>::operator=(const SparseVector& other) {
  PopulateFromSparseVector(other);
  return *this;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::Clear() {
  num_entries_ = EntryIndex(0);
  may_contain_duplicates_ = false;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::ClearAndRelease() {
  capacity_ = EntryIndex(0);
  num_entries_ = EntryIndex(0);
  index_ = nullptr;
  coefficient_ = nullptr;
  buffer_.reset();
  may_contain_duplicates_ = false;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::Reserve(EntryIndex new_capacity) {
  if (new_capacity <= capacity_) return;
  // Round up the capacity to a multiple of four. This way, the start of the
  // coefficient array will be aligned to 16-bytes, provided that the buffer
  // used for storing the data is aligned in that way.
  if (new_capacity.value() & 3) {
    new_capacity += EntryIndex(4 - (new_capacity.value() & 3));
  }

  const size_t index_buffer_size = new_capacity.value() * sizeof(Index);
  const size_t value_buffer_size = new_capacity.value() * sizeof(Fractional);
  const size_t new_buffer_size = index_buffer_size + value_buffer_size;
  std::unique_ptr<char[]> new_buffer(new char[new_buffer_size]);
  IndexType* const new_index = reinterpret_cast<Index*>(new_buffer.get());
  Fractional* const new_coefficient =
      reinterpret_cast<Fractional*>(new_index + new_capacity.value());

  // Avoid copying the data if the vector is empty.
  if (num_entries_ > 0) {
    // NOTE(user): We use memmove instead of std::copy, because the latter
    // leads to naive copying code when used with strong ints (a loop that
    // copies a single 32-bit value in each iteration), and as of 06/2016,
    // memmove is 3-4x faster on Haswell.
    std::memmove(new_index, index_, sizeof(IndexType) * num_entries_.value());
    std::memmove(new_coefficient, coefficient_,
                 sizeof(Fractional) * num_entries_.value());
  }
  std::swap(buffer_, new_buffer);
  index_ = new_index;
  coefficient_ = new_coefficient;
  capacity_ = new_capacity;
}

template <typename IndexType, typename IteratorType>
bool SparseVector<IndexType, IteratorType>::IsEmpty() const {
  return num_entries_ == EntryIndex(0);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::Swap(SparseVector* other) {
  std::swap(buffer_, other->buffer_);
  std::swap(num_entries_, other->num_entries_);
  std::swap(capacity_, other->capacity_);
  std::swap(may_contain_duplicates_, other->may_contain_duplicates_);
  std::swap(index_, other->index_);
  std::swap(coefficient_, other->coefficient_);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::CleanUp() {
  // TODO(user): Implement in-place sorting of the entries and cleanup. The
  // current version converts the data to an array-of-pairs representation that
  // can be sorted easily with std::stable_sort, and the converts the sorted
  // data back to the struct-of-arrays implementation.
  // The current version is ~20% slower than the in-place sort on the
  // array-of-struct representation. It is not visible on GLOP benchmarks, but
  // it increases peak memory usage by ~8%.
  // Implementing in-place search will require either implementing a custom
  // sorting code, or custom iterators that abstract away the internal
  // representation.
  std::vector<std::pair<Index, Fractional>> entries;
  entries.reserve(num_entries_.value());
  for (EntryIndex i(0); i < num_entries_; ++i) {
    entries.emplace_back(GetIndex(i), GetCoefficient(i));
  }
  std::stable_sort(
      entries.begin(), entries.end(),
      [](const std::pair<Index, Fractional>& a,
         const std::pair<Index, Fractional>& b) { return a.first < b.first; });

  EntryIndex new_size(0);
  for (int i = 0; i < num_entries_; ++i) {
    const std::pair<Index, Fractional> entry = entries[i];
    if (entry.second == 0.0) continue;
    if (i + 1 == num_entries_ || entry.first != entries[i + 1].first) {
      MutableIndex(new_size) = entry.first;
      MutableCoefficient(new_size) = entry.second;
      ++new_size;
    }
  }
  ResizeDown(new_size);
  may_contain_duplicates_ = false;
}

template <typename IndexType, typename IteratorType>
bool SparseVector<IndexType, IteratorType>::IsCleanedUp() const {
  Index previous_index(-1);
  for (const EntryIndex i : AllEntryIndices()) {
    const Index index = GetIndex(i);
    if (index <= previous_index || GetCoefficient(i) == 0.0) return false;
    previous_index = index;
  }
  may_contain_duplicates_ = false;
  return true;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::PopulateFromSparseVector(
    const SparseVector& sparse_vector) {
  // Clear the sparse vector before reserving the new capacity. If we didn't do
  // this, Reserve would have to copy the current contents of the vector if it
  // allocated a new buffer. This would be wasteful, since we overwrite it in
  // the next step anyway.
  Clear();
  Reserve(sparse_vector.capacity_);
  // If there are no entries, then sparse_vector.index_ or .coefficient_
  // may be nullptr or invalid, and accessing them in memmove is UB,
  // even if the moved size is zero.
  if (sparse_vector.num_entries_ > 0) {
    // NOTE(user): Using a single memmove would be slightly faster, but it
    // would not work correctly if this already had a greater capacity than
    // sparse_vector, because the coefficient_ pointer would be positioned
    // incorrectly.
    std::memmove(index_, sparse_vector.index_,
                 sizeof(Index) * sparse_vector.num_entries_.value());
    std::memmove(coefficient_, sparse_vector.coefficient_,
                 sizeof(Fractional) * sparse_vector.num_entries_.value());
  }
  num_entries_ = sparse_vector.num_entries_;
  may_contain_duplicates_ = sparse_vector.may_contain_duplicates_;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::PopulateFromDenseVector(
    const DenseVector& dense_vector) {
  Clear();
  const Index num_indices(dense_vector.size());
  for (Index index(0); index < num_indices; ++index) {
    if (dense_vector[index] != 0.0) {
      SetCoefficient(index, dense_vector[index]);
    }
  }
  may_contain_duplicates_ = false;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::AppendEntriesWithOffset(
    const SparseVector& sparse_vector, Index offset) {
  for (const EntryIndex i : sparse_vector.AllEntryIndices()) {
    const Index new_index = offset + sparse_vector.GetIndex(i);
    DCHECK_GE(new_index, 0);
    AddEntry(new_index, sparse_vector.GetCoefficient(i));
  }
  may_contain_duplicates_ = true;
}

template <typename IndexType, typename IteratorType>
bool SparseVector<IndexType, IteratorType>::CheckNoDuplicates(
    StrictITIVector<IndexType, bool>* boolean_vector) const {
  RETURN_VALUE_IF_NULL(boolean_vector, false);
  // Note(user): Using num_entries() or any function that call
  // CheckNoDuplicates() again will cause an infinite loop!
  if (!may_contain_duplicates_ || num_entries_ <= 1) return true;

  // Update size if needed.
  const Index max_index =
      *std::max_element(index_, index_ + num_entries_.value());
  if (boolean_vector->size() <= max_index) {
    boolean_vector->resize(max_index + 1, false);
  }

  may_contain_duplicates_ = false;
  for (const EntryIndex i : AllEntryIndices()) {
    const Index index = GetIndex(i);
    if ((*boolean_vector)[index]) {
      may_contain_duplicates_ = true;
      break;
    }
    (*boolean_vector)[index] = true;
  }

  // Reset boolean_vector to false.
  for (const EntryIndex i : AllEntryIndices()) {
    (*boolean_vector)[GetIndex(i)] = false;
  }
  return !may_contain_duplicates_;
}

template <typename IndexType, typename IteratorType>
bool SparseVector<IndexType, IteratorType>::CheckNoDuplicates() const {
  // Using num_entries() or any function in that will call CheckNoDuplicates()
  // again will cause an infinite loop!
  if (!may_contain_duplicates_ || num_entries_ <= 1) return true;
  StrictITIVector<Index, bool> boolean_vector;
  return CheckNoDuplicates(&boolean_vector);
}

// Do not filter out zero values, as a zero value can be added to reset a
// previous value. Zero values and duplicates will be removed by CleanUp.
template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::SetCoefficient(Index index,
                                                           Fractional value) {
  AddEntry(index, value);
  may_contain_duplicates_ = true;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::DeleteEntry(Index index) {
  DCHECK(CheckNoDuplicates());
  EntryIndex i(0);
  const EntryIndex end(num_entries());
  while (i < end && GetIndex(i) != index) {
    ++i;
  }
  if (i == end) return;
  const int num_moved_entries = (num_entries_ - i).value() - 1;
  std::memmove(index_ + i.value(), index_ + i.value() + 1,
               sizeof(Index) * num_moved_entries);
  std::memmove(coefficient_ + i.value(), coefficient_ + i.value() + 1,
               sizeof(Fractional) * num_moved_entries);
  --num_entries_;
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::RemoveNearZeroEntries(
    Fractional threshold) {
  DCHECK(CheckNoDuplicates());
  EntryIndex new_index(0);
  for (const EntryIndex i : AllEntryIndices()) {
    const Fractional magnitude = fabs(GetCoefficient(i));
    if (magnitude > threshold) {
      MutableIndex(new_index) = GetIndex(i);
      MutableCoefficient(new_index) = GetCoefficient(i);
      ++new_index;
    }
  }
  ResizeDown(new_index);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::RemoveNearZeroEntriesWithWeights(
    Fractional threshold, const DenseVector& weights) {
  DCHECK(CheckNoDuplicates());
  EntryIndex new_index(0);
  for (const EntryIndex i : AllEntryIndices()) {
    if (fabs(GetCoefficient(i)) * weights[GetIndex(i)] > threshold) {
      MutableIndex(new_index) = GetIndex(i);
      MutableCoefficient(new_index) = GetCoefficient(i);
      ++new_index;
    }
  }
  ResizeDown(new_index);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::MoveEntryToFirstPosition(
    Index index) {
  DCHECK(CheckNoDuplicates());
  for (const EntryIndex i : AllEntryIndices()) {
    if (GetIndex(i) == index) {
      std::swap(MutableIndex(EntryIndex(0)), MutableIndex(i));
      std::swap(MutableCoefficient(EntryIndex(0)), MutableCoefficient(i));
      return;
    }
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::MoveEntryToLastPosition(
    Index index) {
  DCHECK(CheckNoDuplicates());
  const EntryIndex last_entry = num_entries() - 1;
  for (const EntryIndex i : AllEntryIndices()) {
    if (GetIndex(i) == index) {
      std::swap(MutableIndex(last_entry), MutableIndex(i));
      std::swap(MutableCoefficient(last_entry), MutableCoefficient(i));
      return;
    }
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::MultiplyByConstant(
    Fractional factor) {
  for (const EntryIndex i : AllEntryIndices()) {
    MutableCoefficient(i) *= factor;
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::ComponentWiseMultiply(
    const DenseVector& factors) {
  for (const EntryIndex i : AllEntryIndices()) {
    MutableCoefficient(i) *= factors[GetIndex(i)];
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::DivideByConstant(
    Fractional factor) {
  for (const EntryIndex i : AllEntryIndices()) {
    MutableCoefficient(i) /= factor;
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::ComponentWiseDivide(
    const DenseVector& factors) {
  for (const EntryIndex i : AllEntryIndices()) {
    MutableCoefficient(i) /= factors[GetIndex(i)];
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::CopyToDenseVector(
    Index num_indices, DenseVector* dense_vector) const {
  RETURN_IF_NULL(dense_vector);
  dense_vector->AssignToZero(num_indices);
  for (const EntryIndex i : AllEntryIndices()) {
    (*dense_vector)[GetIndex(i)] = GetCoefficient(i);
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::PermutedCopyToDenseVector(
    const IndexPermutation& index_perm, Index num_indices,
    DenseVector* dense_vector) const {
  RETURN_IF_NULL(dense_vector);
  dense_vector->AssignToZero(num_indices);
  for (const EntryIndex i : AllEntryIndices()) {
    (*dense_vector)[index_perm[GetIndex(i)]] = GetCoefficient(i);
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::AddMultipleToDenseVector(
    Fractional multiplier, DenseVector* dense_vector) const {
  RETURN_IF_NULL(dense_vector);
  if (multiplier == 0.0) return;
  for (const EntryIndex i : AllEntryIndices()) {
    (*dense_vector)[GetIndex(i)] += multiplier * GetCoefficient(i);
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::
    AddMultipleToSparseVectorAndDeleteCommonIndex(
        Fractional multiplier, Index removed_common_index,
        Fractional drop_tolerance, SparseVector* accumulator_vector) const {
  AddMultipleToSparseVectorInternal(true, multiplier, removed_common_index,
                                    drop_tolerance, accumulator_vector);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::
    AddMultipleToSparseVectorAndIgnoreCommonIndex(
        Fractional multiplier, Index removed_common_index,
        Fractional drop_tolerance, SparseVector* accumulator_vector) const {
  AddMultipleToSparseVectorInternal(false, multiplier, removed_common_index,
                                    drop_tolerance, accumulator_vector);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::AddMultipleToSparseVectorInternal(
    bool delete_common_index, Fractional multiplier, Index common_index,
    Fractional drop_tolerance, SparseVector* accumulator_vector) const {
  // DCHECK that the input is correct.
  DCHECK(IsCleanedUp());
  DCHECK(accumulator_vector->IsCleanedUp());
  DCHECK(CheckNoDuplicates());
  DCHECK(accumulator_vector->CheckNoDuplicates());
  DCHECK_NE(0.0, LookUpCoefficient(common_index));
  DCHECK_NE(0.0, accumulator_vector->LookUpCoefficient(common_index));

  // Implementation notes: we create a temporary SparseVector "c" to hold the
  // result. We call "a" the first vector (i.e. the current object, which will
  // be multiplied by "multiplier"), and "b" the second vector (which will be
  // swapped with "c" at the end to hold the result).
  // We incrementally build c as: a * multiplier + b.
  const SparseVector& a = *this;
  const SparseVector& b = *accumulator_vector;
  SparseVector c;
  EntryIndex ia(0);  // Index in the vector "a"
  EntryIndex ib(0);  // ... and "b"
  EntryIndex ic(0);  // ... and "c"
  const EntryIndex size_a = a.num_entries();
  const EntryIndex size_b = b.num_entries();
  const int size_adjustment = delete_common_index ? -2 : 0;
  const EntryIndex new_size_upper_bound = size_a + size_b + size_adjustment;
  c.Reserve(new_size_upper_bound);
  c.num_entries_ = new_size_upper_bound;
  while ((ia < size_a) && (ib < size_b)) {
    const Index index_a = a.GetIndex(ia);
    const Index index_b = b.GetIndex(ib);
    // Benchmarks done by fdid@ in 2012 showed that it was faster to put the
    // "if" clauses in that specific order.
    if (index_a == index_b) {
      if (index_a != common_index) {
        const Fractional a_coeff_mul = multiplier * a.GetCoefficient(ia);
        const Fractional b_coeff = b.GetCoefficient(ib);
        const Fractional sum = a_coeff_mul + b_coeff;

        // We do not want to leave near-zero entries.
        // TODO(user): expose the tolerance used here.
        if (std::abs(sum) > drop_tolerance) {
          c.MutableIndex(ic) = index_a;
          c.MutableCoefficient(ic) = sum;
          ++ic;
        }
      } else if (!delete_common_index) {
        c.MutableIndex(ic) = b.GetIndex(ib);
        c.MutableCoefficient(ic) = b.GetCoefficient(ib);
        ++ic;
      }
      ++ia;
      ++ib;
    } else if (index_a < index_b) {
      c.MutableIndex(ic) = index_a;
      c.MutableCoefficient(ic) = multiplier * a.GetCoefficient(ia);
      ++ia;
      ++ic;
    } else {  // index_b < index_a
      c.MutableIndex(ic) = b.GetIndex(ib);
      c.MutableCoefficient(ic) = b.GetCoefficient(ib);
      ++ib;
      ++ic;
    }
  }
  while (ia < size_a) {
    c.MutableIndex(ic) = a.GetIndex(ia);
    c.MutableCoefficient(ic) = multiplier * a.GetCoefficient(ia);
    ++ia;
    ++ic;
  }
  while (ib < size_b) {
    c.MutableIndex(ic) = b.GetIndex(ib);
    c.MutableCoefficient(ic) = b.GetCoefficient(ib);
    ++ib;
    ++ic;
  }
  c.ResizeDown(ic);
  c.may_contain_duplicates_ = false;
  c.Swap(accumulator_vector);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::ApplyIndexPermutation(
    const IndexPermutation& index_perm) {
  for (const EntryIndex i : AllEntryIndices()) {
    MutableIndex(i) = index_perm[GetIndex(i)];
  }
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::ApplyPartialIndexPermutation(
    const IndexPermutation& index_perm) {
  EntryIndex new_index(0);
  for (const EntryIndex i : AllEntryIndices()) {
    const Index index = GetIndex(i);
    if (index_perm[index] >= 0) {
      MutableIndex(new_index) = index_perm[index];
      MutableCoefficient(new_index) = GetCoefficient(i);
      ++new_index;
    }
  }
  ResizeDown(new_index);
}

template <typename IndexType, typename IteratorType>
void SparseVector<IndexType, IteratorType>::MoveTaggedEntriesTo(
    const IndexPermutation& index_perm, SparseVector* output) {
  // Note that this function is called many times, so performance does matter
  // and it is why we optimized the "nothing to do" case.
  const EntryIndex end(num_entries_);
  EntryIndex i(0);
  while (true) {
    if (i >= end) return;  // "nothing to do" case.
    if (index_perm[GetIndex(i)] >= 0) break;
    ++i;
  }
  output->AddEntry(GetIndex(i), GetCoefficient(i));
  for (EntryIndex j(i + 1); j < end; ++j) {
    if (index_perm[GetIndex(j)] < 0) {
      MutableIndex(i) = GetIndex(j);
      MutableCoefficient(i) = GetCoefficient(j);
      ++i;
    } else {
      output->AddEntry(GetIndex(j), GetCoefficient(j));
    }
  }
  ResizeDown(i);

  // TODO(user): In the way we use this function, we know that will not
  // happen, but it is better to be careful so we can check that properly in
  // debug mode.
  output->may_contain_duplicates_ = true;
}

template <typename IndexType, typename IteratorType>
Fractional SparseVector<IndexType, IteratorType>::LookUpCoefficient(
    Index index) const {
  Fractional value(0.0);
  for (const EntryIndex i : AllEntryIndices()) {
    if (GetIndex(i) == index) {
      // Keep in mind the vector may contains several entries with the same
      // index. In such a case the last one is returned.
      // TODO(user): investigate whether an optimized version of
      // LookUpCoefficient for "clean" columns yields speed-ups.
      value = GetCoefficient(i);
    }
  }
  return value;
}

template <typename IndexType, typename IteratorType>
bool SparseVector<IndexType, IteratorType>::IsEqualTo(
    const SparseVector& other) const {
  // We do not take into account the mutable value may_contain_duplicates_.
  if (num_entries() != other.num_entries()) return false;
  for (const EntryIndex i : AllEntryIndices()) {
    if (GetIndex(i) != other.GetIndex(i)) return false;
    if (GetCoefficient(i) != other.GetCoefficient(i)) return false;
  }
  return true;
}

template <typename IndexType, typename IteratorType>
std::string SparseVector<IndexType, IteratorType>::DebugString() const {
  std::string s;
  for (const EntryIndex i : AllEntryIndices()) {
    if (i != 0) s += ", ";
    absl::StrAppendFormat(&s, "[%d]=%g", GetIndex(i).value(),
                          GetCoefficient(i));
  }
  return s;
}

}  // namespace glop
}  // namespace operations_research

#endif  // OR_TOOLS_LP_DATA_SPARSE_VECTOR_H_