// matrix/sparse-matrix.cc
  
  // Copyright 2015     Johns Hopkins University (author: Daniel Povey)
  //           2015     Guoguo Chen
  //           2017     Shiyin Kang
  
  // See ../../COPYING for clarification regarding multiple authors
  //
  // 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
  //
  // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
  // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
  // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
  // MERCHANTABLITY OR NON-INFRINGEMENT.
  // See the Apache 2 License for the specific language governing permissions and
  // limitations under the License.
  
  #include <algorithm>
  #include <limits>
  #include <string>
  
  #include "matrix/sparse-matrix.h"
  #include "matrix/kaldi-matrix.h"
  
  namespace kaldi {
  
  template <typename Real>
  std::pair<MatrixIndexT, Real>* SparseVector<Real>::Data() {
    if (pairs_.empty())
      return NULL;
    else
      return &(pairs_[0]);
  }
  
  template <typename Real>
  const std::pair<MatrixIndexT, Real>* SparseVector<Real>::Data() const {
    if (pairs_.empty())
      return NULL;
    else
      return &(pairs_[0]);
  }
  
  template <typename Real>
  Real SparseVector<Real>::Sum() const {
    Real sum = 0;
    for (int32 i = 0; i < pairs_.size(); ++i) {
      sum += pairs_[i].second;
    }
    return sum;
  }
  
  template<typename Real>
  void SparseVector<Real>::Scale(Real alpha) {
    for (int32 i = 0; i < pairs_.size(); ++i)
      pairs_[i].second *= alpha;
  }
  
  template <typename Real>
  template <typename OtherReal>
  void SparseVector<Real>::CopyElementsToVec(VectorBase<OtherReal> *vec) const {
    KALDI_ASSERT(vec->Dim() == this->dim_);
    vec->SetZero();
    OtherReal *other_data = vec->Data();
    typename std::vector<std::pair<MatrixIndexT, Real> >::const_iterator
        iter = pairs_.begin(), end = pairs_.end();
    for (; iter != end; ++iter)
      other_data[iter->first] = iter->second;
  }
  template
  void SparseVector<float>::CopyElementsToVec(VectorBase<float> *vec) const;
  template
  void SparseVector<float>::CopyElementsToVec(VectorBase<double> *vec) const;
  template
  void SparseVector<double>::CopyElementsToVec(VectorBase<float> *vec) const;
  template
  void SparseVector<double>::CopyElementsToVec(VectorBase<double> *vec) const;
  
  template <typename Real>
  template <typename OtherReal>
  void SparseVector<Real>::AddToVec(Real alpha,
                                    VectorBase<OtherReal> *vec) const {
    KALDI_ASSERT(vec->Dim() == dim_);
    OtherReal *other_data = vec->Data();
    typename std::vector<std::pair<MatrixIndexT, Real> >::const_iterator
        iter = pairs_.begin(), end = pairs_.end();
    if (alpha == 1.0) {  // treat alpha==1.0 case specially.
      for (; iter != end; ++iter)
        other_data[iter->first] += iter->second;
    } else {
      for (; iter != end; ++iter)
        other_data[iter->first] += alpha * iter->second;
    }
  }
  
  template
  void SparseVector<float>::AddToVec(float alpha, VectorBase<float> *vec) const;
  template
  void SparseVector<float>::AddToVec(float alpha, VectorBase<double> *vec) const;
  template
  void SparseVector<double>::AddToVec(double alpha, VectorBase<float> *vec) const;
  template
  void SparseVector<double>::AddToVec(double alpha,
                                      VectorBase<double> *vec) const;
  
  template <typename Real>
  template <typename OtherReal>
  void SparseVector<Real>::CopyFromSvec(const SparseVector<OtherReal> &other) {
    dim_ = other.Dim();
    pairs_.clear();
    if (dim_ == 0) return;
    for (int32 i = 0; i < other.NumElements(); ++i) {
      pairs_.push_back(std::make_pair(
          other.GetElement(i).first,
          static_cast<Real>(other.GetElement(i).second)));
    }
  }
  template
  void SparseVector<float>::CopyFromSvec(const SparseVector<float> &svec);
  template
  void SparseVector<float>::CopyFromSvec(const SparseVector<double> &svec);
  template
  void SparseVector<double>::CopyFromSvec(const SparseVector<float> &svec);
  template
  void SparseVector<double>::CopyFromSvec(const SparseVector<double> &svec);
  
  
  template <typename Real>
  SparseVector<Real>& SparseVector<Real>::operator = (
      const SparseVector<Real> &other) {
    this->CopyFromSvec(other);
    dim_ = other.dim_;
    pairs_ = other.pairs_;
    return *this;
  }
  
  template <typename Real>
  void SparseVector<Real>::Swap(SparseVector<Real> *other) {
    pairs_.swap(other->pairs_);
    std::swap(dim_, other->dim_);
  }
  
  template <typename Real>
  void SparseVector<Real>::Write(std::ostream &os, bool binary) const {
    if (binary) {
      WriteToken(os, binary, "SV");
      WriteBasicType(os, binary, dim_);
      MatrixIndexT num_elems = pairs_.size();
      WriteBasicType(os, binary, num_elems);
      typename std::vector<std::pair<MatrixIndexT, Real> >::const_iterator
          iter = pairs_.begin(), end = pairs_.end();
      for (; iter != end; ++iter) {
        WriteBasicType(os, binary, iter->first);
        WriteBasicType(os, binary, iter->second);
      }
    } else {
      // In text-mode, use a human-friendly, script-friendly format;
      // format is "dim=5 [ 0 0.2 3 0.9 ] "
      os << "dim=" << dim_ << " [ ";
      typename std::vector<std::pair<MatrixIndexT, Real> >::const_iterator
          iter = pairs_.begin(), end = pairs_.end();
      for (; iter != end; ++iter)
        os << iter->first << ' ' << iter->second << ' ';
      os << "] ";
    }
  }
  
  
  template <typename Real>
  void SparseVector<Real>::Read(std::istream &is, bool binary) {
    if (binary) {
      ExpectToken(is, binary, "SV");
      ReadBasicType(is, binary, &dim_);
      KALDI_ASSERT(dim_ >= 0);
      int32 num_elems;
      ReadBasicType(is, binary, &num_elems);
      KALDI_ASSERT(num_elems >= 0 && num_elems <= dim_);
      pairs_.resize(num_elems);
      typename std::vector<std::pair<MatrixIndexT, Real> >::iterator
          iter = pairs_.begin(), end = pairs_.end();
      for (; iter != end; ++iter) {
        ReadBasicType(is, binary, &(iter->first));
        ReadBasicType(is, binary, &(iter->second));
      }
    } else {
      // In text-mode, format is "dim=5 [ 0 0.2 3 0.9 ]
      std::string str;
      is >> str;
      if (str.substr(0, 4) != "dim=")
        KALDI_ERR << "Reading sparse vector, expected 'dim=xxx', got " << str;
      std::string dim_str = str.substr(4, std::string::npos);
      std::istringstream dim_istr(dim_str);
      int32 dim = -1;
      dim_istr >> dim;
      if (dim < 0 || dim_istr.fail()) {
        KALDI_ERR << "Reading sparse vector, expected 'dim=[int]', got " << str;
      }
      dim_ = dim;
      is >> std::ws;
      is >> str;
      if (str != "[")
        KALDI_ERR << "Reading sparse vector, expected '[', got " << str;
      pairs_.clear();
      while (1) {
        is >> std::ws;
        if (is.peek() == ']') {
          is.get();
          break;
        }
        MatrixIndexT i;
        BaseFloat p;
        is >> i >> p;
        if (is.fail())
          KALDI_ERR << "Error reading sparse vector, expecting numbers.";
        KALDI_ASSERT(i >= 0 && i < dim
                     && (pairs_.empty() || i > pairs_.back().first));
        pairs_.push_back(std::pair<MatrixIndexT, BaseFloat>(i, p));
      }
    }
  }
  
  
  namespace sparse_vector_utils {
  template <typename Real>
  struct CompareFirst {
    inline bool operator() (const std::pair<MatrixIndexT, Real> &p1,
                             const std::pair<MatrixIndexT, Real> &p2) const {
      return p1.first < p2.first;
    }
  };
  }
  
  template <typename Real>
  SparseVector<Real>::SparseVector(
      MatrixIndexT dim, const std::vector<std::pair<MatrixIndexT, Real> > &pairs):
      dim_(dim),
      pairs_(pairs) {
    std::sort(pairs_.begin(), pairs_.end(),
              sparse_vector_utils::CompareFirst<Real>());
    typename std::vector<std::pair<MatrixIndexT, Real> >::iterator
        out = pairs_.begin(), in = out,  end = pairs_.end();
    // special case: while there is nothing to be changed, skip over
    // initial input (avoids unnecessary copying).
    while (in + 1 < end && in[0].first != in[1].first && in[0].second != 0.0) {
      in++;
      out++;
    }
    while (in < end) {
      // We reach this point only at the first element of
      // each stretch of identical .first elements.
      *out = *in;
      ++in;
      while (in < end && in->first == out->first) {
        out->second += in->second;  // this is the merge operation.
        ++in;
      }
      if (out->second != Real(0.0))  // Don't keep zero elements.
        out++;
    }
    pairs_.erase(out, end);
    if (!pairs_.empty()) {
      // range check.
      KALDI_ASSERT(pairs_.front().first >= 0 && pairs_.back().first < dim_);
    }
  }
  
  template <typename Real>
  void SparseVector<Real>::SetRandn(BaseFloat zero_prob) {
    pairs_.clear();
    KALDI_ASSERT(zero_prob >= 0 && zero_prob <= 1.0);
    for (MatrixIndexT i = 0; i < dim_; i++)
      if (WithProb(1.0 - zero_prob))
        pairs_.push_back(std::pair<MatrixIndexT, Real>(i, RandGauss()));
  }
  
  template <typename Real>
  void SparseVector<Real>::Resize(MatrixIndexT dim,
                                  MatrixResizeType resize_type) {
    if (resize_type != kCopyData || dim == 0)
      pairs_.clear();
    KALDI_ASSERT(dim >= 0);
    if (dim < dim_ && resize_type == kCopyData)
      while (!pairs_.empty() && pairs_.back().first >= dim)
        pairs_.pop_back();
    dim_ = dim;
  }
  
  template <typename Real>
  MatrixIndexT SparseMatrix<Real>::NumRows() const {
    return rows_.size();
  }
  
  template <typename Real>
  MatrixIndexT SparseMatrix<Real>::NumCols() const {
    if (rows_.empty())
      return 0.0;
    else
      return rows_[0].Dim();
  }
  
  template <typename Real>
  MatrixIndexT SparseMatrix<Real>::NumElements() const {
    int32 num_elements = 0;
    for (int32 i = 0; i < rows_.size(); ++i) {
      num_elements += rows_[i].NumElements();
    }
    return num_elements;
  }
  
  template <typename Real>
  SparseVector<Real>* SparseMatrix<Real>::Data() {
    if (rows_.empty())
      return NULL;
    else
      return rows_.data();
  }
  
  template <typename Real>
  const SparseVector<Real>* SparseMatrix<Real>::Data() const {
    if (rows_.empty())
      return NULL;
    else
      return rows_.data();
  }
  
  template <typename Real>
  Real SparseMatrix<Real>::Sum() const {
    Real sum = 0;
    for (int32 i = 0; i < rows_.size(); ++i) {
      sum += rows_[i].Sum();
    }
    return sum;
  }
  
  template<typename Real>
  Real SparseMatrix<Real>::FrobeniusNorm() const {
    Real squared_sum = 0;
    for (int32 i = 0; i < rows_.size(); ++i) {
      const std::pair<MatrixIndexT, Real> *row_data = rows_[i].Data();
      for (int32 j = 0; j < rows_[i].NumElements(); ++j) {
        squared_sum += row_data[j].second * row_data[j].second;
      }
    }
    return std::sqrt(squared_sum);
  }
  
  template <typename Real>
  template <typename OtherReal>
  void SparseMatrix<Real>::CopyToMat(MatrixBase<OtherReal> *other,
                                     MatrixTransposeType trans) const {
    if (trans == kNoTrans) {
      MatrixIndexT num_rows = rows_.size();
      KALDI_ASSERT(other->NumRows() == num_rows);
      for (MatrixIndexT i = 0; i < num_rows; i++) {
        SubVector<OtherReal> vec(*other, i);
        rows_[i].CopyElementsToVec(&vec);
      }
    } else {
      OtherReal *other_col_data = other->Data();
      MatrixIndexT other_stride = other->Stride(),
          num_rows = NumRows(), num_cols = NumCols();
      KALDI_ASSERT(num_rows == other->NumCols() && num_cols == other->NumRows());
      other->SetZero();
      for (MatrixIndexT row = 0; row < num_rows; row++, other_col_data++) {
        const SparseVector<Real> &svec = rows_[row];
        MatrixIndexT num_elems = svec.NumElements();
        const std::pair<MatrixIndexT, Real> *sdata = svec.Data();
        for (MatrixIndexT e = 0; e < num_elems; e++)
          other_col_data[sdata[e].first * other_stride] = sdata[e].second;
      }
    }
  }
  
  template
  void SparseMatrix<float>::CopyToMat(MatrixBase<float> *other,
                                      MatrixTransposeType trans) const;
  template
  void SparseMatrix<float>::CopyToMat(MatrixBase<double> *other,
                                      MatrixTransposeType trans) const;
  template
  void SparseMatrix<double>::CopyToMat(MatrixBase<float> *other,
                                      MatrixTransposeType trans) const;
  template
  void SparseMatrix<double>::CopyToMat(MatrixBase<double> *other,
                                      MatrixTransposeType trans) const;
  
  template <typename Real>
  void SparseMatrix<Real>::CopyElementsToVec(VectorBase<Real> *other) const {
    KALDI_ASSERT(other->Dim() == NumElements());
    Real *dst_data = other->Data();
    int32 dst_index = 0;
    for (int32 i = 0; i < rows_.size(); ++i) {
      for (int32 j = 0; j < rows_[i].NumElements(); ++j) {
        dst_data[dst_index] =
            static_cast<Real>(rows_[i].GetElement(j).second);
        dst_index++;
      }
    }
  }
  
  template<typename Real>
  template<typename OtherReal>
  void SparseMatrix<Real>::CopyFromSmat(const SparseMatrix<OtherReal> &other,
                                        MatrixTransposeType trans) {
    if (trans == kNoTrans) {
      rows_.resize(other.NumRows());
      if (rows_.size() == 0)
        return;
      for (int32 r = 0; r < rows_.size(); ++r) {
        rows_[r].CopyFromSvec(other.Row(r));
      }
    } else {
      std::vector<std::vector<std::pair<MatrixIndexT, Real> > > pairs(
          other.NumCols());
      for (MatrixIndexT i = 0; i < other.NumRows(); ++i) {
        for (int id = 0; id < other.Row(i).NumElements(); ++id) {
          MatrixIndexT j = other.Row(i).GetElement(id).first;
          Real v = static_cast<Real>(other.Row(i).GetElement(id).second);
          pairs[j].push_back( { i, v });
        }
      }
      SparseMatrix<Real> temp(other.NumRows(), pairs);
      Swap(&temp);
    }
  }
  template
  void SparseMatrix<float>::CopyFromSmat(const SparseMatrix<float> &other,
                                         MatrixTransposeType trans);
  template
  void SparseMatrix<float>::CopyFromSmat(const SparseMatrix<double> &other,
                                         MatrixTransposeType trans);
  template
  void SparseMatrix<double>::CopyFromSmat(const SparseMatrix<float> &other,
                                          MatrixTransposeType trans);
  template
  void SparseMatrix<double>::CopyFromSmat(const SparseMatrix<double> &other,
                                          MatrixTransposeType trans);
  
  template <typename Real>
  void SparseMatrix<Real>::Write(std::ostream &os, bool binary) const {
    if (binary) {
      // Note: we can use the same marker for float and double SparseMatrix,
      // because internally we use WriteBasicType and ReadBasicType to read the
      // floats and doubles, and this will automatically take care of type
      // conversion.
      WriteToken(os, binary, "SM");
      int32 num_rows = rows_.size();
      WriteBasicType(os, binary, num_rows);
      for (int32 row = 0; row < num_rows; row++)
        rows_[row].Write(os, binary);
    } else {
      // The format is "rows=10 dim=20 [ 1 0.4  9 1.2 ] dim=20 [ 3 1.7 19 0.6 ] ..
      // not 100% efficient, but easy to work with, and we can re-use the
      // read/write code from SparseVector.
      int32 num_rows = rows_.size();
      os << "rows=" << num_rows << " ";
      for (int32 row = 0; row < num_rows; row++)
        rows_[row].Write(os, binary);
      os << "
  ";  // Might make it a little more readable.
    }
  }
  
  template <typename Real>
  void SparseMatrix<Real>::Read(std::istream &is, bool binary) {
    if (binary) {
      ExpectToken(is, binary, "SM");
      int32 num_rows;
      ReadBasicType(is, binary, &num_rows);
      KALDI_ASSERT(num_rows >= 0 && num_rows < 10000000);
      rows_.resize(num_rows);
      for (int32 row = 0; row < num_rows; row++)
        rows_[row].Read(is, binary);
    } else {
      std::string str;
      is >> str;
      if (str.substr(0, 5) != "rows=")
        KALDI_ERR << "Reading sparse matrix, expected 'rows=xxx', got " << str;
      std::string rows_str = str.substr(5, std::string::npos);
      std::istringstream rows_istr(rows_str);
      int32 num_rows = -1;
      rows_istr >> num_rows;
      if (num_rows < 0 || rows_istr.fail()) {
        KALDI_ERR << "Reading sparse vector, expected 'rows=[int]', got " << str;
      }
      rows_.resize(num_rows);
      for (int32 row = 0; row < num_rows; row++)
        rows_[row].Read(is, binary);
    }
  }
  
  
  template <typename Real>
  void SparseMatrix<Real>::AddToMat(BaseFloat alpha,
                                    MatrixBase<Real> *other,
                                    MatrixTransposeType trans) const {
    if (trans == kNoTrans) {
      MatrixIndexT num_rows = rows_.size();
      KALDI_ASSERT(other->NumRows() == num_rows);
      for (MatrixIndexT i = 0; i < num_rows; i++) {
        SubVector<Real> vec(*other, i);
        rows_[i].AddToVec(alpha, &vec);
      }
    } else {
      Real *other_col_data = other->Data();
      MatrixIndexT other_stride = other->Stride(),
          num_rows = NumRows(), num_cols = NumCols();
      KALDI_ASSERT(num_rows == other->NumCols() && num_cols == other->NumRows());
      for (MatrixIndexT row = 0; row < num_rows; row++, other_col_data++) {
        const SparseVector<Real> &svec = rows_[row];
        MatrixIndexT num_elems = svec.NumElements();
        const std::pair<MatrixIndexT, Real> *sdata = svec.Data();
        for (MatrixIndexT e = 0; e < num_elems; e++)
          other_col_data[sdata[e].first * other_stride] +=
              alpha * sdata[e].second;
      }
    }
  }
  
  template <typename Real>
  Real VecSvec(const VectorBase<Real> &vec,
               const SparseVector<Real> &svec) {
    KALDI_ASSERT(vec.Dim() == svec.Dim());
    MatrixIndexT n = svec.NumElements();
    const std::pair<MatrixIndexT, Real> *sdata = svec.Data();
    const Real *data = vec.Data();
    Real ans = 0.0;
    for (MatrixIndexT i = 0; i < n; i++)
      ans += data[sdata[i].first] * sdata[i].second;
    return ans;
  }
  
  template
  float VecSvec(const VectorBase<float> &vec,
                const SparseVector<float> &svec);
  template
  double VecSvec(const VectorBase<double> &vec,
                const SparseVector<double> &svec);
  
  template <typename Real>
  const SparseVector<Real> &SparseMatrix<Real>::Row(MatrixIndexT r) const {
    KALDI_ASSERT(static_cast<size_t>(r) < rows_.size());
    return rows_[r];
  }
  
  template <typename Real>
  void SparseMatrix<Real>::SetRow(int32 r, const SparseVector<Real> &vec) {
    KALDI_ASSERT(static_cast<size_t>(r) < rows_.size() &&
                 vec.Dim() == rows_[0].Dim());
    rows_[r] = vec;
  }
  
  
  template<typename Real>
  void SparseMatrix<Real>::SelectRows(const std::vector<int32> &row_indexes,
                                      const SparseMatrix<Real> &smat_other) {
    Resize(row_indexes.size(), smat_other.NumCols());
    for (int i = 0; i < row_indexes.size(); ++i) {
      SetRow(i, smat_other.Row(row_indexes[i]));
    }
  }
  
  template<typename Real>
  SparseMatrix<Real>::SparseMatrix(const std::vector<int32> &indexes, int32 dim,
                                   MatrixTransposeType trans) {
    const std::vector<int32>& idx = indexes;
    std::vector<std::vector<std::pair<MatrixIndexT, Real> > > pair(idx.size());
    for (int i = 0; i < idx.size(); ++i) {
      if (idx[i] >= 0) {
        pair[i].push_back( { idx[i], Real(1) });
      }
    }
    SparseMatrix<Real> smat_cpu(dim, pair);
    if (trans == kNoTrans) {
      this->Swap(&smat_cpu);
    } else {
      SparseMatrix<Real> tmp(smat_cpu, kTrans);
      this->Swap(&tmp);
    }
  }
  
  template<typename Real>
  SparseMatrix<Real>::SparseMatrix(const std::vector<int32> &indexes,
                                   const VectorBase<Real> &weights, int32 dim,
                                   MatrixTransposeType trans) {
    const std::vector<int32>& idx = indexes;
    const VectorBase<Real>& w = weights;
    std::vector<std::vector<std::pair<MatrixIndexT, Real> > > pair(idx.size());
    for (int i = 0; i < idx.size(); ++i) {
      if (idx[i] >= 0) {
        pair[i].push_back( { idx[i], w(i) });
      }
    }
    SparseMatrix<Real> smat_cpu(dim, pair);
    if (trans == kNoTrans) {
      this->Swap(&smat_cpu);
    } else {
      SparseMatrix<Real> tmp(smat_cpu, kTrans);
      this->Swap(&tmp);
    }
  }
  
  template <typename Real>
  SparseMatrix<Real>& SparseMatrix<Real>::operator = (
      const SparseMatrix<Real> &other) {
    rows_ = other.rows_;
    return *this;
  }
  
  template <typename Real>
  void SparseMatrix<Real>::Swap(SparseMatrix<Real> *other) {
    rows_.swap(other->rows_);
  }
  
  template<typename Real>
  SparseMatrix<Real>::SparseMatrix(
      MatrixIndexT dim,
      const std::vector<std::vector<std::pair<MatrixIndexT, Real> > > &pairs):
      rows_(pairs.size()) {
    MatrixIndexT num_rows = pairs.size();
    for (MatrixIndexT row = 0; row < num_rows; row++) {
      SparseVector<Real> svec(dim, pairs[row]);
      rows_[row].Swap(&svec);
    }
  }
  
  template <typename Real>
  void SparseMatrix<Real>::SetRandn(BaseFloat zero_prob) {
    MatrixIndexT num_rows = rows_.size();
    for (MatrixIndexT row = 0; row < num_rows; row++)
      rows_[row].SetRandn(zero_prob);
  }
  
  template <typename Real>
  void SparseMatrix<Real>::Resize(MatrixIndexT num_rows,
                                  MatrixIndexT num_cols,
                                  MatrixResizeType resize_type) {
    KALDI_ASSERT(num_rows >= 0 && num_cols >= 0);
    if (resize_type == kSetZero || resize_type == kUndefined) {
      rows_.clear();
      Resize(num_rows, num_cols, kCopyData);
    } else {
      // Assume resize_type == kCopyData from here.
      int32 old_num_rows = rows_.size(), old_num_cols = NumCols();
      SparseVector<Real> initializer(num_cols);
      rows_.resize(num_rows, initializer);
      if (num_cols != old_num_cols)
        for (int32 row = 0; row < old_num_rows; row++)
          rows_[row].Resize(num_cols, kCopyData);
    }
  }
  
  template <typename Real>
  void SparseMatrix<Real>::AppendSparseMatrixRows(
      std::vector<SparseMatrix<Real> > *inputs) {
    rows_.clear();
    size_t num_rows = 0;
    typename std::vector<SparseMatrix<Real> >::iterator
        input_iter = inputs->begin(),
        input_end = inputs->end();
    for (; input_iter != input_end; ++input_iter)
      num_rows += input_iter->rows_.size();
    rows_.resize(num_rows);
    typename std::vector<SparseVector<Real> >::iterator
        row_iter = rows_.begin(),
        row_end = rows_.end();
    for (input_iter = inputs->begin(); input_iter != input_end; ++input_iter) {
      typename std::vector<SparseVector<Real> >::iterator
          input_row_iter = input_iter->rows_.begin(),
          input_row_end = input_iter->rows_.end();
      for (; input_row_iter != input_row_end; ++input_row_iter, ++row_iter)
        row_iter->Swap(&(*input_row_iter));
    }
    KALDI_ASSERT(row_iter == row_end);
    int32 num_cols = NumCols();
    for (row_iter = rows_.begin(); row_iter != row_end; ++row_iter) {
      if (row_iter->Dim() != num_cols)
        KALDI_ERR << "Appending rows with inconsistent dimensions, "
                  << row_iter->Dim() << " vs. " << num_cols;
    }
    inputs->clear();
  }
  
  template<typename Real>
  void SparseMatrix<Real>::Scale(Real alpha) {
    MatrixIndexT num_rows = rows_.size();
    for (MatrixIndexT row = 0; row < num_rows; row++)
      rows_[row].Scale(alpha);
  }
  
  template<typename Real>
  SparseMatrix<Real>::SparseMatrix(const MatrixBase<Real> &mat) {
    MatrixIndexT num_rows = mat.NumRows();
    rows_.resize(num_rows);
    for (int32 row = 0; row < num_rows; row++) {
      SparseVector<Real> this_row(mat.Row(row));
      rows_[row].Swap(&this_row);
    }
  }
  
  template<typename Real>
  Real TraceMatSmat(const MatrixBase<Real> &A,
                    const SparseMatrix<Real> &B,
                    MatrixTransposeType trans) {
    Real sum = 0.0;
    if (trans == kTrans) {
      MatrixIndexT num_rows = A.NumRows();
      KALDI_ASSERT(B.NumRows() == num_rows);
      for (MatrixIndexT r = 0; r < num_rows; r++)
        sum += VecSvec(A.Row(r), B.Row(r));
    } else {
      const Real *A_col_data = A.Data();
      MatrixIndexT Astride = A.Stride(), Acols = A.NumCols(), Arows = A.NumRows();
      KALDI_ASSERT(Arows == B.NumCols() && Acols == B.NumRows());
      sum = 0.0;
      for (MatrixIndexT i = 0; i < Acols; i++, A_col_data++) {
        Real col_sum = 0.0;
        const SparseVector<Real> &svec = B.Row(i);
        MatrixIndexT num_elems = svec.NumElements();
        const std::pair<MatrixIndexT, Real> *sdata = svec.Data();
        for (MatrixIndexT e = 0; e < num_elems; e++)
          col_sum += A_col_data[Astride * sdata[e].first] * sdata[e].second;
        sum += col_sum;
      }
    }
    return sum;
  }
  
  template
  float TraceMatSmat(const MatrixBase<float> &A,
                     const SparseMatrix<float> &B,
                     MatrixTransposeType trans);
  template
  double TraceMatSmat(const MatrixBase<double> &A,
                     const SparseMatrix<double> &B,
                     MatrixTransposeType trans);
  
  void GeneralMatrix::Clear() {
    mat_.Resize(0, 0);
    cmat_.Clear();
    smat_.Resize(0, 0);
  }
  
  GeneralMatrix& GeneralMatrix::operator= (const MatrixBase<BaseFloat> &mat) {
    Clear();
    mat_ = mat;
    return *this;
  }
  
  GeneralMatrix& GeneralMatrix::operator= (const CompressedMatrix &cmat) {
    Clear();
    cmat_ = cmat;
    return *this;
  }
  
  GeneralMatrix& GeneralMatrix::operator= (const SparseMatrix<BaseFloat> &smat) {
    Clear();
    smat_ = smat;
    return *this;
  }
  
  GeneralMatrix& GeneralMatrix::operator= (const GeneralMatrix &gmat) {
    mat_ = gmat.mat_;
    smat_ = gmat.smat_;
    cmat_ = gmat.cmat_;
    return *this;
  }
  
  
  GeneralMatrixType GeneralMatrix::Type() const {
    if (smat_.NumRows() != 0)
      return kSparseMatrix;
    else if (cmat_.NumRows() != 0)
      return kCompressedMatrix;
    else
      return kFullMatrix;
  }
  
  MatrixIndexT GeneralMatrix::NumRows() const {
    MatrixIndexT r = smat_.NumRows();
    if (r != 0)
      return r;
    r = cmat_.NumRows();
    if (r != 0)
      return r;
    return mat_.NumRows();
  }
  
  MatrixIndexT GeneralMatrix::NumCols() const {
    MatrixIndexT r = smat_.NumCols();
    if (r != 0)
      return r;
    r = cmat_.NumCols();
    if (r != 0)
      return r;
    return mat_.NumCols();
  }
  
  
  void GeneralMatrix::Compress() {
    if (mat_.NumRows() != 0) {
      cmat_.CopyFromMat(mat_);
      mat_.Resize(0, 0);
    }
  }
  
  void GeneralMatrix::Uncompress() {
    if (cmat_.NumRows() != 0) {
      mat_.Resize(cmat_.NumRows(), cmat_.NumCols(), kUndefined);
      cmat_.CopyToMat(&mat_);
      cmat_.Clear();
    }
  }
  
  void GeneralMatrix::GetMatrix(Matrix<BaseFloat> *mat) const {
    if (mat_.NumRows() !=0) {
      *mat = mat_;
    } else if (cmat_.NumRows() != 0) {
      mat->Resize(cmat_.NumRows(), cmat_.NumCols(), kUndefined);
      cmat_.CopyToMat(mat);
    } else if (smat_.NumRows() != 0) {
      mat->Resize(smat_.NumRows(), smat_.NumCols(), kUndefined);
      smat_.CopyToMat(mat);
    } else {
      mat->Resize(0, 0);
    }
  }
  
  void GeneralMatrix::CopyToMat(MatrixBase<BaseFloat> *mat,
                                MatrixTransposeType trans) const {
    if (mat_.NumRows() !=0) {
      mat->CopyFromMat(mat_, trans);
    } else if (cmat_.NumRows() != 0) {
      cmat_.CopyToMat(mat, trans);
    } else if (smat_.NumRows() != 0) {
      smat_.CopyToMat(mat, trans);
    } else {
      KALDI_ASSERT(mat->NumRows() == 0);
    }
  }
  
  void GeneralMatrix::Scale(BaseFloat alpha) {
    if (mat_.NumRows() != 0) {
      mat_.Scale(alpha);
    } else if (cmat_.NumRows() != 0) {
      cmat_.Scale(alpha);
    } else if (smat_.NumRows() != 0) {
      smat_.Scale(alpha);
    }
  
  }
  const SparseMatrix<BaseFloat>& GeneralMatrix::GetSparseMatrix() const {
    if (mat_.NumRows() != 0 || cmat_.NumRows() != 0)
      KALDI_ERR << "GetSparseMatrix called on GeneralMatrix of wrong type.";
    return smat_;
  }
  
  void GeneralMatrix::SwapSparseMatrix(SparseMatrix<BaseFloat> *smat) {
    if (mat_.NumRows() != 0 || cmat_.NumRows() != 0)
      KALDI_ERR << "GetSparseMatrix called on GeneralMatrix of wrong type.";
    smat->Swap(&smat_);
  }
  
  void GeneralMatrix::SwapCompressedMatrix(CompressedMatrix *cmat) {
    if (mat_.NumRows() != 0 || smat_.NumRows() != 0)
      KALDI_ERR << "GetSparseMatrix called on GeneralMatrix of wrong type.";
    cmat->Swap(&cmat_);
  }
  
  const CompressedMatrix &GeneralMatrix::GetCompressedMatrix() const {
    if (mat_.NumRows() != 0 || smat_.NumRows() != 0)
      KALDI_ERR << "GetCompressedMatrix called on GeneralMatrix of wrong type.";
    return cmat_;
  }
  
  const Matrix<BaseFloat> &GeneralMatrix::GetFullMatrix() const {
    if (smat_.NumRows() != 0 || cmat_.NumRows() != 0)
      KALDI_ERR << "GetFullMatrix called on GeneralMatrix of wrong type.";
    return mat_;
  }
  
  
  void GeneralMatrix::SwapFullMatrix(Matrix<BaseFloat> *mat) {
    if (cmat_.NumRows() != 0 || smat_.NumRows() != 0)
      KALDI_ERR << "SwapMatrix called on GeneralMatrix of wrong type.";
    mat->Swap(&mat_);
  }
  
  void GeneralMatrix::Write(std::ostream &os, bool binary) const {
    if (smat_.NumRows() != 0) {
      smat_.Write(os, binary);
    } else if (cmat_.NumRows() != 0) {
      cmat_.Write(os, binary);
    } else {
      mat_.Write(os, binary);
    }
  }
  
  void GeneralMatrix::Read(std::istream &is, bool binary) {
    Clear();
    if (binary) {
      int peekval = is.peek();
      if (peekval == 'C') {
        // Token CM for compressed matrix
        cmat_.Read(is, binary);
      } else if (peekval == 'S') {
        // Token SM for sparse matrix
        smat_.Read(is, binary);
      } else {
        mat_.Read(is, binary);
      }
    } else {
      // note: in text mode we will only ever read regular
      // or sparse matrices, because the compressed-matrix format just
      // gets written as a regular matrix in text mode.
      is >> std::ws;  // Eat up white space.
      int peekval = is.peek();
      if (peekval == 'r') {  // sparse format starts rows=[int].
        smat_.Read(is, binary);
      } else {
        mat_.Read(is, binary);
      }
    }
  }
  
  
  void AppendGeneralMatrixRows(const std::vector<const GeneralMatrix *> &src,
                               GeneralMatrix *mat) {
    mat->Clear();
    int32 size = src.size();
    if (size == 0)
      return;
    bool all_sparse = true;
    for (int32 i = 0; i < size; i++) {
      if (src[i]->Type() != kSparseMatrix && src[i]->NumRows() != 0) {
        all_sparse = false;
        break;
      }
    }
    if (all_sparse) {
      std::vector<SparseMatrix<BaseFloat> > sparse_mats(size);
      for (int32 i = 0; i < size; i++)
        sparse_mats[i] = src[i]->GetSparseMatrix();
      SparseMatrix<BaseFloat> appended_mat;
      appended_mat.AppendSparseMatrixRows(&sparse_mats);
      mat->SwapSparseMatrix(&appended_mat);
    } else {
      int32 tot_rows = 0, num_cols = -1;
      for (int32 i = 0; i < size; i++) {
        const GeneralMatrix &src_mat = *(src[i]);
        int32 src_rows = src_mat.NumRows(), src_cols = src_mat.NumCols();
        if (src_rows != 0) {
          tot_rows += src_rows;
          if (num_cols == -1) num_cols = src_cols;
          else if (num_cols != src_cols)
            KALDI_ERR << "Appending rows of matrices with inconsistent num-cols: "
                      << num_cols << " vs. " << src_cols;
        }
      }
      Matrix<BaseFloat> appended_mat(tot_rows, num_cols, kUndefined);
      int32 row_offset = 0;
      for (int32 i = 0; i < size; i++) {
        const GeneralMatrix &src_mat = *(src[i]);
        int32 src_rows = src_mat.NumRows();
        if (src_rows != 0) {
          SubMatrix<BaseFloat> dest_submat(appended_mat, row_offset, src_rows,
                                           0, num_cols);
          src_mat.CopyToMat(&dest_submat);
          row_offset += src_rows;
        }
      }
      KALDI_ASSERT(row_offset == tot_rows);
      mat->SwapFullMatrix(&appended_mat);
    }
  }
  
  void FilterCompressedMatrixRows(const CompressedMatrix &in,
                                  const std::vector<bool> &keep_rows,
                                  Matrix<BaseFloat> *out) {
    KALDI_ASSERT(keep_rows.size() == static_cast<size_t>(in.NumRows()));
    int32 num_kept_rows = 0;
    std::vector<bool>::const_iterator iter = keep_rows.begin(),
                                       end = keep_rows.end();
    for (; iter != end; ++iter)
      if (*iter)
        num_kept_rows++;
    if (num_kept_rows == 0)
      KALDI_ERR << "No kept rows";
    if (num_kept_rows == static_cast<int32>(keep_rows.size())) {
      out->Resize(in.NumRows(), in.NumCols(), kUndefined);
      in.CopyToMat(out);
      return;
    }
    const BaseFloat heuristic = 0.33;
    // should be > 0 and < 1.0.  represents the performance hit we get from
    // iterating row-wise versus column-wise in compressed-matrix uncompression.
  
    if (num_kept_rows > heuristic * in.NumRows()) {
      // if quite a few of the the rows are kept, it may be more efficient
      // to uncompress the entire compressed matrix, since per-column operation
      // is more efficient.
      Matrix<BaseFloat> full_mat(in);
      FilterMatrixRows(full_mat, keep_rows, out);
    } else {
      out->Resize(num_kept_rows, in.NumCols(), kUndefined);
  
      iter = keep_rows.begin();
      int32 out_row = 0;
      for (int32 in_row = 0; iter != end; ++iter, ++in_row) {
        if (*iter) {
          SubVector<BaseFloat> dest(*out, out_row);
          in.CopyRowToVec(in_row, &dest);
          out_row++;
        }
      }
      KALDI_ASSERT(out_row == num_kept_rows);
    }
  }
  
  template <typename Real>
  void FilterMatrixRows(const Matrix<Real> &in,
                        const std::vector<bool> &keep_rows,
                        Matrix<Real> *out) {
    KALDI_ASSERT(keep_rows.size() == static_cast<size_t>(in.NumRows()));
    int32 num_kept_rows = 0;
    std::vector<bool>::const_iterator iter = keep_rows.begin(),
                                       end = keep_rows.end();
    for (; iter != end; ++iter)
      if (*iter)
        num_kept_rows++;
    if (num_kept_rows == 0)
      KALDI_ERR << "No kept rows";
    if (num_kept_rows == static_cast<int32>(keep_rows.size())) {
      *out = in;
      return;
    }
    out->Resize(num_kept_rows, in.NumCols(), kUndefined);
    iter = keep_rows.begin();
    int32 out_row = 0;
    for (int32 in_row = 0; iter != end; ++iter, ++in_row) {
      if (*iter) {
        SubVector<Real> src(in, in_row);
        SubVector<Real> dest(*out, out_row);
        dest.CopyFromVec(src);
        out_row++;
      }
    }
    KALDI_ASSERT(out_row == num_kept_rows);
  }
  
  template
  void FilterMatrixRows(const Matrix<float> &in,
                        const std::vector<bool> &keep_rows,
                        Matrix<float> *out);
  template
  void FilterMatrixRows(const Matrix<double> &in,
                        const std::vector<bool> &keep_rows,
                        Matrix<double> *out);
  
  template <typename Real>
  void FilterSparseMatrixRows(const SparseMatrix<Real> &in,
                              const std::vector<bool> &keep_rows,
                              SparseMatrix<Real> *out) {
    KALDI_ASSERT(keep_rows.size() == static_cast<size_t>(in.NumRows()));
    int32 num_kept_rows = 0;
    std::vector<bool>::const_iterator iter = keep_rows.begin(),
                                       end = keep_rows.end();
    for (; iter != end; ++iter)
      if (*iter)
        num_kept_rows++;
    if (num_kept_rows == 0)
      KALDI_ERR << "No kept rows";
    if (num_kept_rows == static_cast<int32>(keep_rows.size())) {
      *out = in;
      return;
    }
    out->Resize(num_kept_rows, in.NumCols(), kUndefined);
    iter = keep_rows.begin();
    int32 out_row = 0;
    for (int32 in_row = 0; iter != end; ++iter, ++in_row) {
      if (*iter) {
        out->SetRow(out_row, in.Row(in_row));
        out_row++;
      }
    }
    KALDI_ASSERT(out_row == num_kept_rows);
  }
  
  template
  void FilterSparseMatrixRows(const SparseMatrix<float> &in,
                              const std::vector<bool> &keep_rows,
                              SparseMatrix<float> *out);
  template
  void FilterSparseMatrixRows(const SparseMatrix<double> &in,
                              const std::vector<bool> &keep_rows,
                              SparseMatrix<double> *out);
  
  
  void FilterGeneralMatrixRows(const GeneralMatrix &in,
                               const std::vector<bool> &keep_rows,
                               GeneralMatrix *out) {
    out->Clear();
    KALDI_ASSERT(keep_rows.size() == static_cast<size_t>(in.NumRows()));
    int32 num_kept_rows = 0;
    std::vector<bool>::const_iterator iter = keep_rows.begin(),
                                       end = keep_rows.end();
    for (; iter != end; ++iter)
      if (*iter)
        num_kept_rows++;
    if (num_kept_rows == 0)
      KALDI_ERR << "No kept rows";
    if (num_kept_rows == static_cast<int32>(keep_rows.size())) {
      *out = in;
      return;
    }
    switch (in.Type()) {
      case kCompressedMatrix: {
        const CompressedMatrix &cmat = in.GetCompressedMatrix();
        Matrix<BaseFloat> full_mat;
        FilterCompressedMatrixRows(cmat, keep_rows, &full_mat);
        out->SwapFullMatrix(&full_mat);
        return;
      }
      case kSparseMatrix: {
        const SparseMatrix<BaseFloat> &smat = in.GetSparseMatrix();
        SparseMatrix<BaseFloat> smat_out;
        FilterSparseMatrixRows(smat, keep_rows, &smat_out);
        out->SwapSparseMatrix(&smat_out);
        return;
      }
      case kFullMatrix: {
        const Matrix<BaseFloat> &full_mat = in.GetFullMatrix();
        Matrix<BaseFloat> full_mat_out;
        FilterMatrixRows(full_mat, keep_rows, &full_mat_out);
        out->SwapFullMatrix(&full_mat_out);
        return;
      }
      default:
        KALDI_ERR << "Invalid general-matrix type.";
    }
  }
  
  void GeneralMatrix::AddToMat(BaseFloat alpha, MatrixBase<BaseFloat> *mat,
                               MatrixTransposeType trans) const {
    switch (this->Type()) {
      case kFullMatrix: {
        mat->AddMat(alpha, mat_, trans);
        break;
      }
      case kSparseMatrix: {
        smat_.AddToMat(alpha, mat, trans);
        break;
      }
      case kCompressedMatrix: {
        Matrix<BaseFloat> temp_mat(cmat_);
        mat->AddMat(alpha, temp_mat, trans);
        break;
      }
      default:
        KALDI_ERR << "Invalid general-matrix type.";
    }
  }
  
  template <class Real>
  Real SparseVector<Real>::Max(int32 *index_out) const {
    KALDI_ASSERT(dim_ > 0 && pairs_.size() <= static_cast<size_t>(dim_));
    Real ans = -std::numeric_limits<Real>::infinity();
    int32 index = 0;
    typename std::vector<std::pair<MatrixIndexT, Real> >::const_iterator
        iter = pairs_.begin(), end = pairs_.end();
    for (; iter != end; ++iter) {
      if (iter->second > ans) {
        ans = iter->second;
        index = iter->first;
      }
    }
    if (ans >= 0 || pairs_.size() == dim_) {
      // ans >= 0 will be the normal case.
      // if pairs_.size() == dim_ then we need to return
      // even a negative answer as there are no spaces (hence no unlisted zeros).
      *index_out = index;
      return ans;
    }
    // all the stored elements are < 0, but there are unlisted
    // elements -> pick the first unlisted element.
    // Note that this class requires that the indexes are sorted
    // and unique.
    index = 0;  // "index" will always be the next index, that
                // we haven't seen listed yet.
    iter = pairs_.begin();
    for (; iter != end; ++iter) {
      if (iter->first > index) {  // index "index" is not listed.
        *index_out = index;
        return 0.0;
      } else {
        // index is the next potential gap in the indexes.
        index = iter->first + 1;
      }
    }
    // we can reach here if either pairs_.empty(), or
    // pairs_ is nonempty but contains a sequence (0, 1, 2,...).
    if (!pairs_.empty())
      index = pairs_.back().first + 1;
    // else leave index at zero
    KALDI_ASSERT(index < dim_);
    *index_out = index;
    return 0.0;
  }
  
  template <typename Real>
  SparseVector<Real>::SparseVector(const VectorBase<Real> &vec) {
    MatrixIndexT dim = vec.Dim();
    dim_ = dim;
    if (dim == 0)
      return;
    const Real *ptr = vec.Data();
    for (MatrixIndexT i = 0; i < dim; i++) {
      Real val = ptr[i];
      if (val != 0.0)
        pairs_.push_back(std::pair<MatrixIndexT,Real>(i,val));
    }
  }
  
  void GeneralMatrix::Swap(GeneralMatrix *other) {
    mat_.Swap(&(other->mat_));
    cmat_.Swap(&(other->cmat_));
    smat_.Swap(&(other->smat_));
  }
  
  
  void ExtractRowRangeWithPadding(
      const GeneralMatrix &in,
      int32 row_offset,
      int32 num_rows,
      GeneralMatrix *out) {
    // make sure 'out' is empty to start with.
    Matrix<BaseFloat> empty_mat;
    *out = empty_mat;
    if (num_rows == 0) return;
    switch (in.Type()) {
      case kFullMatrix: {
        const Matrix<BaseFloat> &mat_in = in.GetFullMatrix();
        int32 num_rows_in = mat_in.NumRows(), num_cols = mat_in.NumCols();
        KALDI_ASSERT(num_rows_in > 0);  // we can't extract >0 rows from an empty
                                        // matrix.
        Matrix<BaseFloat> mat_out(num_rows, num_cols, kUndefined);
        for (int32 row = 0; row < num_rows; row++) {
          int32 row_in = row + row_offset;
          if (row_in < 0) row_in = 0;
          else if (row_in >= num_rows_in) row_in = num_rows_in - 1;
          SubVector<BaseFloat> vec_in(mat_in, row_in),
              vec_out(mat_out, row);
          vec_out.CopyFromVec(vec_in);
        }
        out->SwapFullMatrix(&mat_out);
        break;
      }
      case kSparseMatrix: {
        const SparseMatrix<BaseFloat> &smat_in = in.GetSparseMatrix();
        int32 num_rows_in = smat_in.NumRows(),
            num_cols = smat_in.NumCols();
        KALDI_ASSERT(num_rows_in > 0);  // we can't extract >0 rows from an empty
                                        // matrix.
        SparseMatrix<BaseFloat> smat_out(num_rows, num_cols);
        for (int32 row = 0; row < num_rows; row++) {
          int32 row_in = row + row_offset;
          if (row_in < 0) row_in = 0;
          else if (row_in >= num_rows_in) row_in = num_rows_in - 1;
          smat_out.SetRow(row, smat_in.Row(row_in));
        }
        out->SwapSparseMatrix(&smat_out);
        break;
      }
      case kCompressedMatrix: {
        const CompressedMatrix &cmat_in = in.GetCompressedMatrix();
        bool allow_padding = true;
        CompressedMatrix cmat_out(cmat_in, row_offset, num_rows,
                                  0, cmat_in.NumCols(), allow_padding);
        out->SwapCompressedMatrix(&cmat_out);
        break;
      }
      default:
        KALDI_ERR << "Bad matrix type.";
    }
  }
  
  
  
  template class SparseVector<float>;
  template class SparseVector<double>;
  template class SparseMatrix<float>;
  template class SparseMatrix<double>;
  
  }  // namespace kaldi