// matrix/compressed-matrix.cc

// Copyright 2012    Johns Hopkins University (author: Daniel Povey)
//                   Frantisek Skala, Wei Shi
//           2015    Tom Ko

// 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 "matrix/compressed-matrix.h"
#include <algorithm>

namespace kaldi {

//static
MatrixIndexT CompressedMatrix::DataSize(const GlobalHeader &header) {
  // Returns size in bytes of the data.
  DataFormat format = static_cast<DataFormat>(header.format);
  if (format == kOneByteWithColHeaders) {
    return sizeof(GlobalHeader) +
        header.num_cols * (sizeof(PerColHeader) + header.num_rows);
  } else if (format == kTwoByte) {
    return sizeof(GlobalHeader) +
        2 * header.num_rows * header.num_cols;
  } else {
    KALDI_ASSERT(format == kOneByte);
    return sizeof(GlobalHeader) +
        header.num_rows * header.num_cols;
  }
}

// scale all element of matrix by scaling floats
// in GlobalHeader with alpha.
void CompressedMatrix::Scale(float alpha) {
  if (data_ != NULL) {
    GlobalHeader *h = reinterpret_cast<GlobalHeader*>(data_);
    // scale the floating point values in each PerColHolder
    // and leave all integers the same.
    h->min_value *= alpha;
    h->range *= alpha;
  }
}

template<typename Real>  // static inline
void CompressedMatrix::ComputeGlobalHeader(
    const MatrixBase<Real> &mat, CompressionMethod method,
    GlobalHeader *header) {
  if (method == kAutomaticMethod) {
    if (mat.NumRows() > 8) method = kSpeechFeature;
    else method = kTwoByteAuto;
  }

  switch (method) {
    case kSpeechFeature:
      header->format = static_cast<int32>(kOneByteWithColHeaders);  // 1.
      break;
    case kTwoByteAuto: case kTwoByteSignedInteger:
      header->format = static_cast<int32>(kTwoByte);  // 2.
      break;
    case kOneByteAuto: case kOneByteUnsignedInteger: case kOneByteZeroOne:
      header->format = static_cast<int32>(kOneByte);  // 3.
      break;
    default:
      KALDI_ERR << "Invalid compression type: "
                << static_cast<int32>(method);
  }

  header->num_rows = mat.NumRows();
  header->num_cols = mat.NumCols();

  // Now compute 'min_value' and 'range'.
  switch (method) {
    case kSpeechFeature: case kTwoByteAuto: case kOneByteAuto: {
      float min_value = mat.Min(), max_value = mat.Max();
      // ensure that max_value is strictly greater than min_value, even if matrix is
      // constant; this avoids crashes in ComputeColHeader when compressing speech
      // featupres.
      if (max_value == min_value)
        max_value = min_value + (1.0 + fabs(min_value));
      KALDI_ASSERT(min_value - min_value == 0 &&
                   max_value - max_value == 0 &&
                   "Cannot compress a matrix with Nan's or Inf's");

      header->min_value = min_value;
      header->range = max_value - min_value;

      // we previously checked that max_value != min_value, so their
      // difference should be nonzero.
      KALDI_ASSERT(header->range > 0.0);
      break;
    }
    case kTwoByteSignedInteger: {
      header->min_value = -32768.0;
      header->range = 65535.0;
      break;
    }
    case kOneByteUnsignedInteger: {
      header->min_value = 0.0;
      header->range = 255.0;
      break;
    }
    case kOneByteZeroOne: {
      header->min_value = 0.0;
      header->range = 1.0;
      break;
    }
    default:
      KALDI_ERR << "Unknown compression method = "
                << static_cast<int32>(method);
  }
  KALDI_COMPILE_TIME_ASSERT(sizeof(*header) == 20);  // otherwise
  // something weird is happening and our code probably won't work or
  // won't be robust across platforms.
}

template<typename Real>
void CompressedMatrix::CopyFromMat(
    const MatrixBase<Real> &mat, CompressionMethod method) {
  if (data_ != NULL) {
    delete [] static_cast<float*>(data_);  // call delete [] because was allocated with new float[]
    data_ = NULL;
  }
  if (mat.NumRows() == 0) { return; }  // Zero-size matrix stored as zero pointer.


  GlobalHeader global_header;
  ComputeGlobalHeader(mat, method, &global_header);

  int32 data_size = DataSize(global_header);

  data_ = AllocateData(data_size);

  *(reinterpret_cast<GlobalHeader*>(data_)) = global_header;

  DataFormat format = static_cast<DataFormat>(global_header.format);
  if (format == kOneByteWithColHeaders) {
    PerColHeader *header_data =
        reinterpret_cast<PerColHeader*>(static_cast<char*>(data_) +
                                        sizeof(GlobalHeader));
    uint8 *byte_data =
        reinterpret_cast<uint8*>(header_data + global_header.num_cols);

    const Real *matrix_data = mat.Data();

    for (int32 col = 0; col < global_header.num_cols; col++) {
      CompressColumn(global_header,
                     matrix_data + col, mat.Stride(),
                     global_header.num_rows,
                     header_data, byte_data);
      header_data++;
      byte_data += global_header.num_rows;
    }
  } else if (format == kTwoByte) {
    uint16 *data = reinterpret_cast<uint16*>(static_cast<char*>(data_) +
                                             sizeof(GlobalHeader));
    int32 num_rows = mat.NumRows(), num_cols = mat.NumCols();
    for (int32 r = 0; r < num_rows; r++) {
      const Real *row_data = mat.RowData(r);
      for (int32 c = 0; c < num_cols; c++)
        data[c] = FloatToUint16(global_header, row_data[c]);
      data += num_cols;
    }
  } else {
    KALDI_ASSERT(format == kOneByte);
    uint8 *data = reinterpret_cast<uint8*>(static_cast<char*>(data_) +
                                           sizeof(GlobalHeader));
    int32 num_rows = mat.NumRows(), num_cols = mat.NumCols();
    for (int32 r = 0; r < num_rows; r++) {
      const Real *row_data = mat.RowData(r);
      for (int32 c = 0; c < num_cols; c++)
        data[c] = FloatToUint8(global_header, row_data[c]);
      data += num_cols;
    }
  }
}

// Instantiate the template for float and double.
template
void CompressedMatrix::CopyFromMat(const MatrixBase<float> &mat,
                                   CompressionMethod method);

template
void CompressedMatrix::CopyFromMat(const MatrixBase<double> &mat,
                                   CompressionMethod method);


CompressedMatrix::CompressedMatrix(
    const CompressedMatrix &cmat,
    const MatrixIndexT row_offset,
    const MatrixIndexT num_rows,
    const MatrixIndexT col_offset,
    const MatrixIndexT num_cols,
    bool allow_padding): data_(NULL) {
  int32 old_num_rows = cmat.NumRows(), old_num_cols = cmat.NumCols();

  if (old_num_rows == 0) {
    KALDI_ASSERT(num_rows == 0 && num_cols == 0);
    // The empty matrix is stored as a zero pointer.
    return;
  }

  KALDI_ASSERT(row_offset < old_num_rows);
  KALDI_ASSERT(col_offset < old_num_cols);
  KALDI_ASSERT(row_offset >= 0 || allow_padding);
  KALDI_ASSERT(col_offset >= 0);
  KALDI_ASSERT(row_offset + num_rows <= old_num_rows || allow_padding);
  KALDI_ASSERT(col_offset + num_cols <= old_num_cols);

  if (num_rows == 0 || num_cols == 0) { return; }

  bool padding_is_used = (row_offset < 0 ||
                          row_offset + num_rows > old_num_rows);

  GlobalHeader new_global_header;
  KALDI_COMPILE_TIME_ASSERT(sizeof(new_global_header) == 20);

  GlobalHeader *old_global_header = reinterpret_cast<GlobalHeader*>(cmat.Data());

  new_global_header = *old_global_header;
  new_global_header.num_cols = num_cols;
  new_global_header.num_rows = num_rows;

  // We don't switch format from 1 -> 2 (in case of size reduction) yet; if this
  // is needed, we will do this below by creating a temporary Matrix.
  new_global_header.format = old_global_header->format;

  data_ = AllocateData(DataSize(new_global_header));  // allocate memory
  *(reinterpret_cast<GlobalHeader*>(data_)) = new_global_header;


  DataFormat format = static_cast<DataFormat>(old_global_header->format);
  if (format == kOneByteWithColHeaders) {
    PerColHeader *old_per_col_header =
        reinterpret_cast<PerColHeader*>(old_global_header + 1);
    uint8 *old_byte_data =
        reinterpret_cast<uint8*>(old_per_col_header +
                                 old_global_header->num_cols);
    PerColHeader *new_per_col_header =
        reinterpret_cast<PerColHeader*>(
            reinterpret_cast<GlobalHeader*>(data_) + 1);

    memcpy(new_per_col_header, old_per_col_header + col_offset,
           sizeof(PerColHeader) * num_cols);

    uint8 *new_byte_data =
        reinterpret_cast<uint8*>(new_per_col_header + num_cols);
    if (!padding_is_used) {
      uint8 *old_start_of_subcol =
          old_byte_data + row_offset + (col_offset * old_num_rows),
          *new_start_of_col = new_byte_data;
      for (int32 i = 0; i < num_cols; i++) {
        memcpy(new_start_of_col, old_start_of_subcol, num_rows);
        new_start_of_col += num_rows;
        old_start_of_subcol += old_num_rows;
      }
    } else {
      uint8 *old_start_of_col =
          old_byte_data + (col_offset * old_num_rows),
          *new_start_of_col = new_byte_data;
      for (int32 i = 0; i < num_cols; i++) {

        for (int32 j = 0; j < num_rows; j++) {
          int32 old_j = j + row_offset;
          if (old_j < 0) old_j = 0;
          else if (old_j >= old_num_rows) old_j = old_num_rows - 1;
          new_start_of_col[j] = old_start_of_col[old_j];
        }
        new_start_of_col += num_rows;
        old_start_of_col += old_num_rows;
      }
    }
  } else if (format == kTwoByte) {
    const uint16 *old_data =
        reinterpret_cast<const uint16*>(old_global_header + 1);
    uint16 *new_row_data =
        reinterpret_cast<uint16*>(reinterpret_cast<GlobalHeader*>(data_) + 1);

    for (int32 row = 0; row < num_rows; row++) {
      int32 old_row = row + row_offset;
      // The next two lines are only relevant if padding_is_used.
      if (old_row < 0) old_row = 0;
      else if (old_row >= old_num_rows) old_row = old_num_rows - 1;
      const uint16 *old_row_data =
          old_data + col_offset + (old_num_cols * old_row);
      memcpy(new_row_data, old_row_data, sizeof(uint16) * num_cols);
      new_row_data += num_cols;
    }
  } else {
    KALDI_ASSERT(format == kOneByte);
    const uint8 *old_data =
        reinterpret_cast<const uint8*>(old_global_header + 1);
    uint8 *new_row_data =
        reinterpret_cast<uint8*>(reinterpret_cast<GlobalHeader*>(data_) + 1);

    for (int32 row = 0; row < num_rows; row++) {
      int32 old_row = row + row_offset;
      // The next two lines are only relevant if padding_is_used.
      if (old_row < 0) old_row = 0;
      else if (old_row >= old_num_rows) old_row = old_num_rows - 1;
      const uint8 *old_row_data =
          old_data + col_offset + (old_num_cols * old_row);
      memcpy(new_row_data, old_row_data, sizeof(uint8) * num_cols);
      new_row_data += num_cols;
    }
  }

  if (num_rows < 8 && format == kOneByteWithColHeaders) {
    // format was 1 but we want it to be 2 -> create a temporary
    // Matrix (uncompress), re-compress, and swap.
    // This gives us almost exact reconstruction while saving
    // memory (the elements take more space but there will be
    // no per-column headers).
    Matrix<float> temp(this->NumRows(), this->NumCols(),
                       kUndefined);
    this->CopyToMat(&temp);
    CompressedMatrix temp_cmat(temp, kTwoByteAuto);
    this->Swap(&temp_cmat);
  }
}


template<typename Real>
CompressedMatrix &CompressedMatrix::operator =(const MatrixBase<Real> &mat) {
  this->CopyFromMat(mat);
  return *this;
}

// Instantiate the template for float and double.
template
CompressedMatrix& CompressedMatrix::operator =(const MatrixBase<float> &mat);

template
CompressedMatrix& CompressedMatrix::operator =(const MatrixBase<double> &mat);

inline uint16 CompressedMatrix::FloatToUint16(
    const GlobalHeader &global_header,
    float value) {
  float f = (value - global_header.min_value) /
      global_header.range;
  if (f > 1.0) f = 1.0;  // Note: this should not happen.
  if (f < 0.0) f = 0.0;  // Note: this should not happen.
  return static_cast<int>(f * 65535 + 0.499);  // + 0.499 is to
  // round to closest int; avoids bias.
}


inline uint8 CompressedMatrix::FloatToUint8(
    const GlobalHeader &global_header,
    float value) {
  float f = (value - global_header.min_value) /
      global_header.range;
  if (f > 1.0) f = 1.0;  // Note: this should not happen.
  if (f < 0.0) f = 0.0;  // Note: this should not happen.
  return static_cast<int>(f * 255 + 0.499);  // + 0.499 is to
  // round to closest int; avoids bias.
}


inline float CompressedMatrix::Uint16ToFloat(
    const GlobalHeader &global_header,
    uint16 value) {
  // the constant 1.52590218966964e-05 is 1/65535.
  return global_header.min_value
      + global_header.range * 1.52590218966964e-05F * value;
}

template<typename Real>  // static
void CompressedMatrix::ComputeColHeader(
    const GlobalHeader &global_header,
    const Real *data, MatrixIndexT stride,
    int32 num_rows, CompressedMatrix::PerColHeader *header) {
  KALDI_ASSERT(num_rows > 0);
  std::vector<Real> sdata(num_rows); // the sorted data.
  for (size_t i = 0, size = sdata.size(); i < size; i++)
    sdata[i] = data[i*stride];

  if (num_rows >= 5) {
    int quarter_nr = num_rows/4;
    // std::sort(sdata.begin(), sdata.end());
    // The elements at positions 0, quarter_nr,
    // 3*quarter_nr, and num_rows-1 need to be in sorted order.
    std::nth_element(sdata.begin(), sdata.begin() + quarter_nr, sdata.end());
    // Now, sdata.begin() + quarter_nr contains the element that would appear
    // in sorted order, in that position.
    std::nth_element(sdata.begin(), sdata.begin(), sdata.begin() + quarter_nr);
    // Now, sdata.begin() and sdata.begin() + quarter_nr contain the elements
    // that would appear at those positions in sorted order.
    std::nth_element(sdata.begin() + quarter_nr + 1,
                     sdata.begin() + (3*quarter_nr), sdata.end());
    // Now, sdata.begin(), sdata.begin() + quarter_nr, and sdata.begin() +
    // 3*quarter_nr, contain the elements that would appear at those positions
    // in sorted order.
    std::nth_element(sdata.begin() + (3*quarter_nr) + 1, sdata.end() - 1,
                     sdata.end());
    // Now, sdata.begin(), sdata.begin() + quarter_nr, and sdata.begin() +
    // 3*quarter_nr, and sdata.end() - 1, contain the elements that would appear
    // at those positions in sorted order.

    header->percentile_0 =
        std::min<uint16>(FloatToUint16(global_header, sdata[0]), 65532);
    header->percentile_25 =
        std::min<uint16>(
            std::max<uint16>(
                FloatToUint16(global_header, sdata[quarter_nr]),
                header->percentile_0 + static_cast<uint16>(1)), 65533);
    header->percentile_75 =
        std::min<uint16>(
            std::max<uint16>(
                FloatToUint16(global_header, sdata[3*quarter_nr]),
                header->percentile_25 + static_cast<uint16>(1)), 65534);
    header->percentile_100 = std::max<uint16>(
        FloatToUint16(global_header, sdata[num_rows-1]),
        header->percentile_75 + static_cast<uint16>(1));

  } else {  // handle this pathological case.
    std::sort(sdata.begin(), sdata.end());
    // Note: we know num_rows is at least 1.
    header->percentile_0 =
        std::min<uint16>(FloatToUint16(global_header, sdata[0]),
                         65532);
    if (num_rows > 1)
      header->percentile_25 =
          std::min<uint16>(
              std::max<uint16>(FloatToUint16(global_header, sdata[1]),
                               header->percentile_0 + 1), 65533);
    else
      header->percentile_25 = header->percentile_0 + 1;
    if (num_rows > 2)
      header->percentile_75 =
          std::min<uint16>(
              std::max<uint16>(FloatToUint16(global_header, sdata[2]),
                               header->percentile_25 + 1), 65534);
    else
      header->percentile_75 = header->percentile_25 + 1;
    if (num_rows > 3)
      header->percentile_100 =
          std::max<uint16>(FloatToUint16(global_header, sdata[3]),
                           header->percentile_75 + 1);
    else
      header->percentile_100 = header->percentile_75 + 1;
  }
}

// static
inline uint8 CompressedMatrix::FloatToChar(
    float p0, float p25, float p75, float p100,
    float value) {
  int ans;
  if (value < p25) {  // range [ p0, p25 ) covered by
    // characters 0 .. 64.  We round to the closest int.
    float f = (value - p0) / (p25 - p0);
    ans = static_cast<int>(f * 64 + 0.5);
    // Note: the checks on the next two lines
    // are necessary in pathological cases when all the elements in a row
    // are the same and the percentile_* values are separated by one.
    if (ans < 0) ans = 0;
    if (ans > 64) ans = 64;
  } else if (value < p75) {  // range [ p25, p75 )covered
    // by characters 64 .. 192.  We round to the closest int.
    float f = (value - p25) / (p75 - p25);
    ans = 64 + static_cast<int>(f * 128 + 0.5);
    if (ans < 64) ans = 64;
    if (ans > 192) ans = 192;
  } else {  // range [ p75, p100 ] covered by
    // characters 192 .. 255.  Note: this last range
    // has fewer characters than the left range, because
    // we go up to 255, not 256.
    float f = (value - p75) / (p100 - p75);
    ans = 192 + static_cast<int>(f * 63 + 0.5);
    if (ans < 192) ans = 192;
    if (ans > 255) ans = 255;
  }
  return static_cast<uint8>(ans);
}


// static
inline float CompressedMatrix::CharToFloat(
    float p0, float p25, float p75, float p100,
    uint8 value) {
  if (value <= 64) {
    return p0 + (p25 - p0) * value * (1/64.0);
  } else if (value <= 192) {
    return p25 + (p75 - p25) * (value - 64) * (1/128.0);
  } else {
    return p75 + (p100 - p75) * (value - 192) * (1/63.0);
  }
}


template<typename Real>  // static
void CompressedMatrix::CompressColumn(
    const GlobalHeader &global_header,
    const Real *data, MatrixIndexT stride,
    int32 num_rows, CompressedMatrix::PerColHeader *header,
    uint8 *byte_data) {
  ComputeColHeader(global_header, data, stride,
                   num_rows, header);

  float p0 = Uint16ToFloat(global_header, header->percentile_0),
      p25 = Uint16ToFloat(global_header, header->percentile_25),
      p75 = Uint16ToFloat(global_header, header->percentile_75),
      p100 = Uint16ToFloat(global_header, header->percentile_100);

  for (int32 i = 0; i < num_rows; i++) {
    Real this_data = data[i * stride];
    byte_data[i] = FloatToChar(p0, p25, p75, p100, this_data);
  }
}

// static
void* CompressedMatrix::AllocateData(int32 num_bytes) {
  KALDI_ASSERT(num_bytes > 0);
  KALDI_COMPILE_TIME_ASSERT(sizeof(float) == 4);
  // round size up to nearest number of floats.
  return reinterpret_cast<void*>(new float[(num_bytes/3) + 4]);
}

void CompressedMatrix::Write(std::ostream &os, bool binary) const {
  if (binary) {  // Binary-mode write:
    if (data_ != NULL) {
      GlobalHeader &h = *reinterpret_cast<GlobalHeader*>(data_);
      DataFormat format = static_cast<DataFormat>(h.format);
      if (format == kOneByteWithColHeaders) {
        WriteToken(os, binary, "CM");
      } else if (format == kTwoByte) {
        WriteToken(os, binary, "CM2");
      } else if (format == kOneByte) {
        WriteToken(os, binary, "CM3");
      }
      MatrixIndexT size = DataSize(h);  // total size of data in data_
      // We don't write out the "int32 format", hence the + 4, - 4.
      os.write(reinterpret_cast<const char*>(data_) + 4, size - 4);
    } else {  // special case: where data_ == NULL, we treat it as an empty
      // matrix.
      WriteToken(os, binary, "CM");
      GlobalHeader h;
      h.range = h.min_value = 0.0;
      h.num_rows = h.num_cols = 0;
      os.write(reinterpret_cast<const char*>(&h), sizeof(h));
    }
  } else {
    // In text mode, just use the same format as a regular matrix.
    // This is not compressed.
    Matrix<BaseFloat> temp_mat(this->NumRows(), this->NumCols(),
                               kUndefined);
    this->CopyToMat(&temp_mat);
    temp_mat.Write(os, binary);
  }
  if (os.fail())
    KALDI_ERR << "Error writing compressed matrix to stream.";
}

void CompressedMatrix::Read(std::istream &is, bool binary) {
  if (data_ != NULL) {
    delete [] (static_cast<float*>(data_));
    data_ = NULL;
  }
  if (binary) {
    int peekval = Peek(is, binary);
    if (peekval == 'C') {
      std::string tok; // Should be CM (format 1) or CM2 (format 2)
      ReadToken(is, binary, &tok);
      GlobalHeader h;
      if (tok == "CM") { h.format = 1; } //  kOneByteWithColHeaders
      else if (tok == "CM2") { h.format = 2; }  // kTwoByte
      else if (tok == "CM3") { h.format = 3; }  // kOneByte
      else {
        KALDI_ERR << "Unexpected token " << tok << ", expecting CM, CM2 or CM3";
      }
      // don't read the "format" -> hence + 4, - 4.
      is.read(reinterpret_cast<char*>(&h) + 4, sizeof(h) - 4);
      if (is.fail())
        KALDI_ERR << "Failed to read header";
      if (h.num_cols == 0) // empty matrix.
        return;
      int32 size = DataSize(h), remaining_size = size - sizeof(GlobalHeader);
      data_ = AllocateData(size);
      *(reinterpret_cast<GlobalHeader*>(data_)) = h;
      is.read(reinterpret_cast<char*>(data_) + sizeof(GlobalHeader),
              remaining_size);
    } else {
      // Assume that what we're reading is a regular Matrix.  This might be the
      // case if you changed your code, making a Matrix into a CompressedMatrix,
      // and you want back-compatibility for reading.
      Matrix<BaseFloat> M;
      M.Read(is, binary); // This will crash if it was not a Matrix.
      this->CopyFromMat(M);
    }
  } else {  // Text-mode read.  In this case you don't get to
    // choose the compression type.  Anyway this branch would only
    // be taken when debugging.
    Matrix<BaseFloat> temp;
    temp.Read(is, binary);
    this->CopyFromMat(temp);
  }
  if (is.fail())
    KALDI_ERR << "Failed to read data.";
}

template<typename Real>
void CompressedMatrix::CopyToMat(MatrixBase<Real> *mat,
                                 MatrixTransposeType trans) const {
  if (trans == kTrans) {
    Matrix<Real> temp(this->NumCols(), this->NumRows());
    CopyToMat(&temp, kNoTrans);
    mat->CopyFromMat(temp, kTrans);
    return;
  }

  if (data_ == NULL) {
    KALDI_ASSERT(mat->NumRows() == 0);
    KALDI_ASSERT(mat->NumCols() == 0);
    return;
  }
  GlobalHeader *h = reinterpret_cast<GlobalHeader*>(data_);
  int32 num_cols = h->num_cols, num_rows = h->num_rows;
  KALDI_ASSERT(mat->NumRows() == num_rows);
  KALDI_ASSERT(mat->NumCols() == num_cols);

  DataFormat format = static_cast<DataFormat>(h->format);
  if (format == kOneByteWithColHeaders) {
    PerColHeader *per_col_header = reinterpret_cast<PerColHeader*>(h+1);
    uint8 *byte_data = reinterpret_cast<uint8*>(per_col_header +
                                                h->num_cols);
    for (int32 i = 0; i < num_cols; i++, per_col_header++) {
      float p0 = Uint16ToFloat(*h, per_col_header->percentile_0),
          p25 = Uint16ToFloat(*h, per_col_header->percentile_25),
          p75 = Uint16ToFloat(*h, per_col_header->percentile_75),
          p100 = Uint16ToFloat(*h, per_col_header->percentile_100);
      for (int32 j = 0; j < num_rows; j++, byte_data++) {
        float f = CharToFloat(p0, p25, p75, p100, *byte_data);
        (*mat)(j, i) = f;
      }
    }
  } else if (format == kTwoByte) {
    const uint16 *data = reinterpret_cast<const uint16*>(h + 1);
    float min_value = h->min_value,
        increment = h->range * (1.0 / 65535.0);
    for (int32 i = 0; i < num_rows; i++) {
      Real *row_data = mat->RowData(i);
      for (int32 j = 0; j < num_cols; j++)
        row_data[j] = min_value + data[j] * increment;
      data += num_cols;
    }
  } else {
    KALDI_ASSERT(format == kOneByte);
    float min_value = h->min_value, increment = h->range * (1.0 / 255.0);

    const uint8 *data = reinterpret_cast<const uint8*>(h + 1);
    for (int32 i = 0; i < num_rows; i++) {
      Real *row_data = mat->RowData(i);
      for (int32 j = 0; j < num_cols; j++)
        row_data[j] = min_value + data[j] * increment;
      data += num_cols;
    }
  }
}

// Instantiate the template for float and double.
template
void CompressedMatrix::CopyToMat(MatrixBase<float> *mat,
                                 MatrixTransposeType trans) const;
template
void CompressedMatrix::CopyToMat(MatrixBase<double> *mat,
                                 MatrixTransposeType trans) const;

template<typename Real>
void CompressedMatrix::CopyRowToVec(MatrixIndexT row,
                                    VectorBase<Real> *v) const {
  KALDI_ASSERT(row < this->NumRows());
  KALDI_ASSERT(row >= 0);
  KALDI_ASSERT(v->Dim() == this->NumCols());

  GlobalHeader *h = reinterpret_cast<GlobalHeader*>(data_);
  DataFormat format = static_cast<DataFormat>(h->format);
  if (format == kOneByteWithColHeaders) {
    PerColHeader *per_col_header = reinterpret_cast<PerColHeader*>(h+1);
    uint8 *byte_data = reinterpret_cast<uint8*>(per_col_header +
                                                h->num_cols);
    byte_data += row;  // point to first value we are interested in
    for (int32 i = 0; i < h->num_cols;
         i++, per_col_header++, byte_data += h->num_rows) {
      float p0 = Uint16ToFloat(*h, per_col_header->percentile_0),
          p25 = Uint16ToFloat(*h, per_col_header->percentile_25),
          p75 = Uint16ToFloat(*h, per_col_header->percentile_75),
          p100 = Uint16ToFloat(*h, per_col_header->percentile_100);
      float f = CharToFloat(p0, p25, p75, p100, *byte_data);
      (*v)(i) = f;
    }
  } else if (format == kTwoByte) {
    int32 num_cols = h->num_cols;
    float min_value = h->min_value,
        increment = h->range * (1.0 / 65535.0);
    const uint16 *row_data = reinterpret_cast<uint16*>(h + 1) + (num_cols * row);
    Real *v_data = v->Data();
    for (int32 c = 0; c < num_cols; c++)
      v_data[c] = min_value + row_data[c] * increment;
  } else {
    KALDI_ASSERT(format == kOneByte);
    int32 num_cols = h->num_cols;
    float min_value = h->min_value,
        increment = h->range * (1.0 / 255.0);
    const uint8 *row_data = reinterpret_cast<uint8*>(h + 1) + (num_cols * row);
    Real *v_data = v->Data();
    for (int32 c = 0; c < num_cols; c++)
      v_data[c] = min_value + row_data[c] * increment;
  }
}

template<typename Real>
void CompressedMatrix::CopyColToVec(MatrixIndexT col,
                                    VectorBase<Real> *v) const {
  KALDI_ASSERT(col < this->NumCols());
  KALDI_ASSERT(col >= 0);
  KALDI_ASSERT(v->Dim() == this->NumRows());

  GlobalHeader *h = reinterpret_cast<GlobalHeader*>(data_);

  DataFormat format = static_cast<DataFormat>(h->format);
  if (format == kOneByteWithColHeaders) {
    PerColHeader *per_col_header = reinterpret_cast<PerColHeader*>(h+1);
    uint8 *byte_data = reinterpret_cast<uint8*>(per_col_header +
                                                h->num_cols);
    byte_data += col*h->num_rows;  // point to first value in the column we want
    per_col_header += col;
    float p0 = Uint16ToFloat(*h, per_col_header->percentile_0),
        p25 = Uint16ToFloat(*h, per_col_header->percentile_25),
        p75 = Uint16ToFloat(*h, per_col_header->percentile_75),
        p100 = Uint16ToFloat(*h, per_col_header->percentile_100);
    for (int32 i = 0; i < h->num_rows; i++, byte_data++) {
      float f = CharToFloat(p0, p25, p75, p100, *byte_data);
      (*v)(i) = f;
    }
  } else if (format == kTwoByte) {
    int32 num_rows = h->num_rows, num_cols = h->num_cols;
    float min_value = h->min_value,
        increment = h->range * (1.0 / 65535.0);
    const uint16 *col_data = reinterpret_cast<uint16*>(h + 1) + col;
    Real *v_data = v->Data();
    for (int32 r = 0; r < num_rows; r++)
      v_data[r] = min_value + increment * col_data[r * num_cols];
  } else {
    KALDI_ASSERT(format == kOneByte);
    int32 num_rows = h->num_rows, num_cols = h->num_cols;
    float min_value = h->min_value,
        increment = h->range * (1.0 / 255.0);
    const uint8 *col_data = reinterpret_cast<uint8*>(h + 1) + col;
    Real *v_data = v->Data();
    for (int32 r = 0; r < num_rows; r++)
      v_data[r] = min_value + increment * col_data[r * num_cols];
  }
}

// instantiate the templates.
template void
CompressedMatrix::CopyColToVec(MatrixIndexT, VectorBase<double> *) const;
template void
CompressedMatrix::CopyColToVec(MatrixIndexT, VectorBase<float> *) const;
template void
CompressedMatrix::CopyRowToVec(MatrixIndexT, VectorBase<double> *) const;
template void
CompressedMatrix::CopyRowToVec(MatrixIndexT, VectorBase<float> *) const;

template<typename Real>
void CompressedMatrix::CopyToMat(int32 row_offset,
                                 int32 col_offset,
                                 MatrixBase<Real> *dest) const {
  KALDI_PARANOID_ASSERT(row_offset < this->NumRows());
  KALDI_PARANOID_ASSERT(col_offset < this->NumCols());
  KALDI_PARANOID_ASSERT(row_offset >= 0);
  KALDI_PARANOID_ASSERT(col_offset >= 0);
  KALDI_ASSERT(row_offset+dest->NumRows() <= this->NumRows());
  KALDI_ASSERT(col_offset+dest->NumCols() <= this->NumCols());
  // everything is OK
  GlobalHeader *h = reinterpret_cast<GlobalHeader*>(data_);
  int32 num_rows = h->num_rows, num_cols = h->num_cols,
      tgt_cols = dest->NumCols(), tgt_rows = dest->NumRows();

  DataFormat format = static_cast<DataFormat>(h->format);
  if (format == kOneByteWithColHeaders) {
    PerColHeader *per_col_header = reinterpret_cast<PerColHeader*>(h+1);
    uint8 *byte_data = reinterpret_cast<uint8*>(per_col_header +
                                                h->num_cols);

    uint8 *start_of_subcol = byte_data+row_offset;  // skip appropriate
    // number of columns
    start_of_subcol += col_offset*num_rows;  // skip appropriate number of rows

    per_col_header += col_offset;  // skip the appropriate number of headers

    for (int32 i = 0;
         i < tgt_cols;
         i++, per_col_header++, start_of_subcol+=num_rows) {
      byte_data = start_of_subcol;
      float p0 = Uint16ToFloat(*h, per_col_header->percentile_0),
          p25 = Uint16ToFloat(*h, per_col_header->percentile_25),
          p75 = Uint16ToFloat(*h, per_col_header->percentile_75),
          p100 = Uint16ToFloat(*h, per_col_header->percentile_100);
      for (int32 j = 0; j < tgt_rows; j++, byte_data++) {
        float f = CharToFloat(p0, p25, p75, p100, *byte_data);
        (*dest)(j, i) = f;
      }
    }
  } else if (format == kTwoByte) {
    const uint16 *data = reinterpret_cast<const uint16*>(h+1) + col_offset +
        (num_cols * row_offset);
    float min_value = h->min_value,
        increment = h->range * (1.0 / 65535.0);

    for (int32 row = 0; row < tgt_rows; row++) {
      Real *dest_row = dest->RowData(row);
      for (int32 col = 0; col < tgt_cols; col++)
        dest_row[col] = min_value + increment * data[col];
      data += num_cols;
    }
  } else {
    KALDI_ASSERT(format == kOneByte);
    const uint8 *data = reinterpret_cast<const uint8*>(h+1) + col_offset +
        (num_cols * row_offset);
    float min_value = h->min_value,
        increment = h->range * (1.0 / 255.0);
    for (int32 row = 0; row < tgt_rows; row++) {
      Real *dest_row = dest->RowData(row);
      for (int32 col = 0; col < tgt_cols; col++)
        dest_row[col] = min_value + increment * data[col];
      data += num_cols;
    }
  }
}

// instantiate the templates.
template void CompressedMatrix::CopyToMat(int32,
                                          int32,
                                          MatrixBase<float> *dest) const;
template void CompressedMatrix::CopyToMat(int32,
                                          int32,
                                          MatrixBase<double> *dest) const;

void CompressedMatrix::Clear() {
  if (data_ != NULL) {
    delete [] static_cast<float*>(data_);
    data_ = NULL;
  }
}

CompressedMatrix::CompressedMatrix(const CompressedMatrix &mat): data_(NULL) {
  *this = mat; // use assignment operator.
}

CompressedMatrix &CompressedMatrix::operator = (const CompressedMatrix &mat) {
  Clear(); // now this->data_ == NULL.
  if (mat.data_ != NULL) {
    MatrixIndexT data_size = DataSize(*static_cast<GlobalHeader*>(mat.data_));
    data_ = AllocateData(data_size);
    memcpy(static_cast<void*>(data_),
           static_cast<void*>(mat.data_),
           data_size);
  }
  return *this;
}


}  // namespace kaldi