// 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