matrix.h

/* MLPACK 0.2
 *
 * Copyright (c) 2008, 2009 Alexander Gray,
 *                          Garry Boyer,
 *                          Ryan Riegel,
 *                          Nikolaos Vasiloglou,
 *                          Dongryeol Lee,
 *                          Chip Mappus, 
 *                          Nishant Mehta,
 *                          Hua Ouyang,
 *                          Parikshit Ram,
 *                          Long Tran,
 *                          Wee Chin Wong
 *
 * Copyright (c) 2008, 2009 Georgia Institute of Technology
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301, USA.
 */
#ifndef LA_MATRIX_H
#define LA_MATRIX_H

#include "fastlib/base/base.h"
#ifndef DISABLE_DISK_MATRIX
#include "fastlib/mmanager/memory_manager.h"
#endif

#include <stdlib.h>
#include <string.h>
#include <math.h>

template<typename T>
class GenVector {
 private:
  T *ptr_;
  index_t length_;
  bool should_free_;
  
  OBJECT_TRAVERSAL_ONLY(GenVector) {
    OT_OBJ(length_);
    OT_ALLOC(ptr_, length_);
  }
  OT_REFILL_TRANSIENTS(GenVector) {
    should_free_ = true;
  }
  
 public:
  GenVector() {
    DEBUG_ONLY(Uninitialize_());
  }
  
  GenVector(const GenVector& other) {
    DEBUG_ONLY(Uninitialize_());
    Copy(other);
  }
  ASSIGN_VIA_COPY_CONSTRUCTION(GenVector);
  
  ~GenVector() {
    Destruct();
  }
  
  void Destruct() {
    DEBUG_ASSERT_MSG(ptr_ != BIG_BAD_POINTER(T),
       "You forgot to initialize a Vector before it got automatically freed.");
    
    /* mark slow case as "unlikely" even if it might be the likely case */
    if (unlikely(should_free_)) {
      mem::DebugPoison(ptr_, length_);
      mem::Free(ptr_);
    }
    
    DEBUG_ONLY(Uninitialize_());
  }

  void Init(index_t in_length) {
    ptr_ = mem::Alloc<T>(in_length);
    length_ = in_length;
    should_free_ = true;
  }
  
  void StaticInit(index_t in_length) {
#ifndef DISABLE_DISK_MATRIX
    ptr_ = mmapmm::MemoryManager<false>::allocator_->Alloc<T>(in_length);
    length_ = in_length;
    should_free_ = false;
#else
    Init(in_length);
#endif
  }

  void SetAll(T d) {
    mem::RepeatConstruct(ptr_, d, length_);
  }
  
  void SetZero() {
    // TODO: if IEEE is used, this can be done efficiently with memset
    SetAll(0);
  }

  void Copy(const GenVector& other) {
    Copy(other.ptr(), other.length());
  }

  void Copy(const T *doubles, index_t in_length) {
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = mem::AllocCopy(doubles, in_length);
    length_ = in_length;
    should_free_ = true;
  }

  void StaticCopy(const GenVector& other) {
    StaticCopy(other.ptr(), other.length());
  }

  void StaticCopy(const T *doubles, index_t in_length) {
#ifndef DISABLE_DISK_MATRIX
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = mmapmm::MemoryManager<false>::allocator_->Alloc<T>(in_length);
    mem::Copy<T, T, T>(ptr_, doubles, in_length);
    length_ = in_length;
    should_free_ = false;
#else
    Copy(doubles, in_length);
#endif
  }
  
  void Alias(T *in_ptr, index_t in_length) {
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = in_ptr;
    length_ = in_length;
    should_free_ = false;
  }
  
  void WeakCopy(const GenVector& other) {
    Alias(other);
  }
  
  void Alias(const GenVector& other) {
    // we trust in good faith that a const vector won't be abused
    Alias(other.ptr_, other.length());
  }
  
  void Own(GenVector* other) {
    Own(other->ptr_, other->length());
    
    DEBUG_ASSERT(other->should_free_);
    other->should_free_ = false;
  }
  
  void Own(T *in_ptr, index_t in_length) {
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = in_ptr;
    length_ = in_length;
    should_free_ = true;
  }
  
  void StaticOwn(GenVector* other) {
#ifndef DISABLE_DISK_MATRIX
    StaticOwn(other->ptr_, other->length());
#else
    Own(other);
#endif
  }
  
  void StaticOwn(T *in_ptr, index_t in_length) {
#ifndef DISABLE_DISK_MATRIX
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = in_ptr;
    length_ = in_length;
    should_free_ = false;
#else
    Own(in_ptr, in_length);
#endif
  }
  
  void MakeSubvector(index_t start_index, index_t len, GenVector* dest) {
    DEBUG_BOUNDS(start_index, length_);
    DEBUG_BOUNDS(start_index + len, length_ + 1);
    
    dest->Alias(ptr_ + start_index, len);
  }
  
  void SwapValues(GenVector* other) {
    DEBUG_ASSERT(length() == other->length());
    mem::Swap(ptr_, other->ptr_, length_);
  }
  
  void CopyValues(const GenVector& other) {
    DEBUG_ASSERT(length() == other.length());
    mem::Copy(ptr_, other.ptr_, length_);
  }

  void CopyValues(const T *src_ptr) {
    mem::Copy(ptr_, src_ptr, length_);
  }
  
  void PrintDebug(const char *name = "", FILE *stream = stderr) const {
    fprintf(stream, "----- VECTOR %s ------\n", name);
    for (index_t i = 0; i < length(); i++) {
      fprintf(stream, "%+3.3f ", get(i));
    }
    fprintf(stream, "\n");
  }
  
 public:
  index_t length() const {
    return length_;
  }
  
  T *ptr() {
    return ptr_;
  }
  
  const T *ptr() const {
    return ptr_;
  }
  
  T operator [] (index_t i) const {
    DEBUG_BOUNDS(i, length_);
    return ptr_[i];
  }
  
  T &operator [] (index_t i) {
    DEBUG_BOUNDS(i, length_);
    return ptr_[i];
  }
  
  T get(index_t i) const {
    DEBUG_BOUNDS(i, length_);
    return ptr_[i];
  }
  
 private:
  void AssertUninitialized_() const {
    DEBUG_ASSERT_MSG(length_ == BIG_BAD_NUMBER, "Cannot re-init vectors.");
  }
  
  void Uninitialize_() {
    DEBUG_ONLY(ptr_ = BIG_BAD_POINTER(T));
    DEBUG_ONLY(length_ = BIG_BAD_NUMBER);
  }
  
  void AssertInitialized_() {
    DEBUG_ASSERT_MSG(ptr_ != BIG_BAD_POINTER(T),
        "Vector was not initialized.");
  }
};

typedef GenVector<double> Vector;

template<typename T>
class GenMatrix {
 private:
  T *ptr_;
  index_t n_rows_;
  index_t n_cols_;
  bool should_free_;

  OBJECT_TRAVERSAL_ONLY(GenMatrix) {
    OT_OBJ(n_rows_);
    OT_OBJ(n_cols_);
    OT_ALLOC(ptr_, n_elements());
  }
  OT_REFILL_TRANSIENTS(GenMatrix) {
    should_free_ = false;
  }

 public:
  GenMatrix(index_t in_rows, index_t in_cols) {
    DEBUG_ONLY(Uninitialize_());
    Init(in_rows, in_cols);
  }

  GenMatrix(const GenMatrix<T>& other) {
    DEBUG_ONLY(Uninitialize_());
    Copy(other);
  }
  ASSIGN_VIA_COPY_CONSTRUCTION(GenMatrix);

  GenMatrix() {
    DEBUG_ONLY(Uninitialize_());
  }

  ~GenMatrix() {
    Destruct();
  }
  
  void Destruct() {
    DEBUG_ASSERT_MSG(ptr_ != BIG_BAD_POINTER(T),
       "You forgot to initialize a Matrix before it got automatically freed.");
    if (unlikely(should_free_)) {
      mem::DebugPoison(ptr_, n_rows_ * n_cols_);
      mem::Free(ptr_);
      DEBUG_ONLY(Uninitialize_());
    }
    DEBUG_POISON_PTR(ptr_);
    DEBUG_ONLY(n_rows_ = BIG_BAD_NUMBER);
    DEBUG_ONLY(n_cols_ = BIG_BAD_NUMBER);
  }

  void Init(index_t in_rows, index_t in_cols) {
    DEBUG_ONLY(AssertUninitialized_());
    ptr_ = mem::Alloc<T>(in_rows * in_cols);
    n_rows_ = in_rows;
    n_cols_ = in_cols;
    should_free_ = true;
  }

  void InitDiagonal(const GenVector<T>& v) {
    Init(v.length(), v.length());
    SetDiagonal(v);
  }

  void StaticInit(index_t in_rows, index_t in_cols) {
#ifndef DISABLE_DISK_MATRIX
    DEBUG_ONLY(AssertUninitialized_());
    ptr_ = mmapmm::MemoryManager<false>::allocator_->Alloc<T>
      (in_rows * in_cols);
    n_rows_ = in_rows;
    n_cols_ = in_cols;
    should_free_ = false;
#else
    Init(in_rows, in_cols);
#endif
  }

  void StaticInitDiagonal(const GenVector<T>& v) {
    StaticInit(v.length(), v.length());
    SetDiagonal(v);
  }

  void SetAll(T d) {
    mem::RepeatConstruct(ptr_, d, n_elements());
  }

  void SetZero() {
    // TODO: If IEEE floating point is used, this can just be a memset to
    // zero
    SetAll(0);
  }

  void SetDiagonal(const GenVector<T>& v) {
    DEBUG_ASSERT(n_rows() == v.length());
    DEBUG_ASSERT(n_cols() == v.length());
    SetZero();
    index_t n = v.length();
    for (index_t i = 0; i < n; i++) {
      set(i, i, v[i]);
    }
  }

  void Copy(const GenMatrix& other) {
    Copy(other.ptr(), other.n_rows(), other.n_cols());    
  }

  void Copy(const T *ptr_in, index_t n_rows_in, index_t n_cols_in) {
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = mem::AllocCopy(ptr_in, n_rows_in * n_cols_in);
    n_rows_ = n_rows_in;
    n_cols_ = n_cols_in;
    should_free_ = true;
  }

  void StaticCopy(const GenMatrix& other) {
    StaticCopy(other.ptr(), other.n_rows(), other.n_cols());    
  }

  void StaticCopy(const T *ptr_in, index_t n_rows_in, index_t n_cols_in) {
#ifndef DISABLE_DISK_MATRIX
    DEBUG_ONLY(AssertUninitialized_());

    ptr_ = mmapmm::MemoryManager<false>::allocator_->Alloc<T>
      (n_rows_in * n_cols_in);
    mem::Copy<T, T, T>(ptr_, ptr_in, n_rows_in * n_cols_in);

    n_rows_ = n_rows_in;
    n_cols_ = n_cols_in;
    should_free_ = false;
#else
    Copy(ptr_in, n_rows_in, n_cols_in);
#endif
  }
 
  void Alias(const GenMatrix& other) {
    // we trust in good faith that const-ness won't be abused
    Alias(other.ptr_, other.n_rows(), other.n_cols());
  }
  
  void Alias(T *ptr_in, index_t n_rows_in, index_t n_cols_in) {
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = ptr_in;
    n_rows_ = n_rows_in;
    n_cols_ = n_cols_in;
    should_free_ = false;
  }
  
  void AliasRowVector(const GenVector<T>& row_vector) {
    Alias(const_cast<T*>(row_vector.ptr()), 1, row_vector.length());
  }
  
  void AliasColVector(const GenVector<T>& col_vector) {
    Alias(const_cast<T*>(col_vector.ptr()), col_vector.length(), 1);
  }
  
  void WeakCopy(const GenMatrix& other) {
    Alias(other);
  }
  
  void Own(GenMatrix* other) {
    Own(other->ptr(), other->n_rows(), other->n_cols());
    
    DEBUG_ASSERT(other->should_free_);
    other->should_free_ = false;
  }
  
  void Own(T *ptr_in, index_t n_rows_in, index_t n_cols_in) {
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = ptr_in;
    n_rows_ = n_rows_in;
    n_cols_ = n_cols_in;
    should_free_ = true;
  }

  void StaticOwn(GenMatrix* other) {
#ifndef DISABLE_DISK_MATRIX
    StaticOwn(other->ptr(), other->n_rows(), other->n_cols());
#else
    Own(other);
#endif
  }
  
  void StaticOwn(T *ptr_in, index_t n_rows_in, index_t n_cols_in) {
#ifndef DISABLE_DISK_MATRIX
    DEBUG_ONLY(AssertUninitialized_());
    
    ptr_ = ptr_in;
    n_rows_ = n_rows_in;
    n_cols_ = n_cols_in;
    should_free_ = false;
#else
    Own(ptr_in, n_rows_in, n_cols_in);
#endif
  }

  void MakeColumnSlice(index_t start_col, index_t n_cols_new,
      GenMatrix *dest) const {
    DEBUG_BOUNDS(start_col, n_cols_);
    DEBUG_BOUNDS(start_col + n_cols_new, n_cols_ + 1);
    dest->Alias(ptr_ + start_col * n_rows_,
        n_rows_, n_cols_new);
  }
  
  void MakeReshaped(index_t n_rows_in, index_t n_cols_in,
      GenMatrix *dest) const {
    DEBUG_ASSERT(n_rows_in * n_cols_in == n_rows() * n_cols());
    dest->Alias(ptr_, n_rows_in, n_cols_in);
  }
  
  void MakeColumnVector(index_t col, GenVector<T> *dest) const {
    DEBUG_BOUNDS(col, n_cols_);
    dest->Alias(n_rows_ * col + ptr_, n_rows_);
  }
  
  void MakeColumnSubvector(index_t col, index_t start_row, index_t n_rows_new,
      GenVector<T> *dest) const {
    DEBUG_BOUNDS(col, n_cols_);
    DEBUG_BOUNDS(start_row, n_rows_);
    DEBUG_BOUNDS(start_row + n_rows_new, n_rows_ + 1);
    dest->Alias(n_rows_ * col + start_row + ptr_, n_rows_new);
  }
  
  T *GetColumnPtr(index_t col) {
    DEBUG_BOUNDS(col, n_cols_);
    return n_rows_ * col + ptr_;
  }
  
  const T *GetColumnPtr(index_t col) const {
    DEBUG_BOUNDS(col, n_cols_);
    return n_rows_ * col + ptr_;
  }
  
   void CopyColumnFromMat(index_t col1, index_t col2, GenMatrix<T> &mat) {
     DEBUG_BOUNDS(col1, n_cols_);
     DEBUG_BOUNDS(col2, mat.n_cols());
     DEBUG_ASSERT(n_rows_==mat.n_rows());
     memcpy(ptr_ + n_rows_ * col1, mat.GetColumnPtr(col2), n_rows_*sizeof(T));
   }
   void CopyColumnFromMat(index_t col1, index_t col2, index_t ncols, GenMatrix<T> &mat) {
     DEBUG_BOUNDS(col1, n_cols_);
     DEBUG_BOUNDS(col2, mat.n_cols());
     DEBUG_BOUNDS(col1+ncols-1, n_cols_);
     DEBUG_BOUNDS(col2+ncols-1, mat.n_cols());
     DEBUG_ASSERT(n_rows_==mat.n_rows());
     memcpy(ptr_ + n_rows_ * col1, mat.GetColumnPtr(col2), ncols*n_rows_*sizeof(T));
   }

   void CopyVectorToColumn(index_t col, GenVector<T> &vec) {
     DEBUG_BOUNDS(col, n_cols_);
     memcpy(ptr_ + n_rows_ * col, vec.ptr(), n_rows_*sizeof(T));
   }
 
  void ResizeNoalias(index_t new_n_cols) {
    DEBUG_ASSERT(should_free_); // the best assert we can do
    n_cols_ = new_n_cols;
    ptr_ = mem::Realloc(ptr_, n_elements());
  }
  
  void SwapValues(GenMatrix* other) {
    DEBUG_ASSERT(n_cols() == other->n_cols());
    DEBUG_ASSERT(n_rows() == other->n_rows());
    mem::Swap(ptr_, other->ptr_, n_elements());
  }
  
  void CopyValues(const GenMatrix& other) {
    DEBUG_ASSERT(n_rows() == other.n_rows());
    DEBUG_ASSERT(n_cols() == other.n_cols());
    mem::Copy(ptr_, other.ptr_, n_elements());
  }

  void PrintDebug(const char *name = "", FILE *stream = stderr) const {
    fprintf(stream, "----- MATRIX %s ------\n", name);
    for (index_t r = 0; r < n_rows(); r++) {
      for (index_t c = 0; c < n_cols(); c++) {
        fprintf(stream, "%+3.3f ", get(r, c));
      }
      fprintf(stream, "\n");
    }
  }
  
 public:
  const T *ptr() const {
    return ptr_;
  }
  
  T *ptr() {
    return ptr_;
  }
  
  T get(index_t r, index_t c) const {
    DEBUG_BOUNDS(r, n_rows_);
    DEBUG_BOUNDS(c, n_cols_);
    return ptr_[c * n_rows_ + r];
  }
 
  void set(index_t r, index_t c, T v) {
    DEBUG_BOUNDS(r, n_rows_);
    DEBUG_BOUNDS(c, n_cols_);
    ptr_[c * n_rows_ + r] = v;
  }
  
  T &ref(index_t r, index_t c) {
    DEBUG_BOUNDS(r, n_rows_);
    DEBUG_BOUNDS(c, n_cols_);
    return ptr_[c * n_rows_ + r];
  }
  
  index_t n_cols() const {
    return n_cols_;
  }
  
  index_t n_rows() const {
    return n_rows_;
  }
  
  size_t n_elements() const {
    // TODO: putting the size_t on the outside may be faster (32-bit
    // versus 64-bit multiplication in cases) but is more likely to result
    // in bugs
    return size_t(n_rows_) * size_t(n_cols_);
  }
  
 private:
  void AssertUninitialized_() const {
    DEBUG_ASSERT_MSG(n_rows_ == BIG_BAD_NUMBER, "Cannot re-init matrices.");
  }
  
  void Uninitialize_() {
    DEBUG_POISON_PTR(ptr_);
    DEBUG_ONLY(n_rows_ = BIG_BAD_NUMBER);
    DEBUG_ONLY(n_cols_ = BIG_BAD_NUMBER);
  } 

};

template<int t_length>
class SmallVector : public Vector {
 private:
  double array_[t_length];
  
 public:
  SmallVector() {
    Alias(array_, t_length);
  }
  ~SmallVector() {}
  
 public:
  index_t length() const {
    return t_length;
  }

  double *ptr() {
    return array_;
  }
  
  const double *ptr() const {
    return array_;
  }
  
  double operator [] (index_t i) const {
    DEBUG_BOUNDS(i, t_length);
    return array_[i];
  }
  
  double &operator [] (index_t i) {
    DEBUG_BOUNDS(i, t_length);
    return array_[i];
  }
  
  double get(index_t i) const {
    DEBUG_BOUNDS(i, t_length);
    return array_[i];
  }
};

typedef GenMatrix<double> Matrix;

template<int t_rows, int t_cols>
class SmallMatrix : public Matrix {
 private:
  double array_[t_cols][t_rows];

 public:
  SmallMatrix() {
    Alias(array_[0], t_rows, t_cols);
  }
  ~SmallMatrix() {}

 public:
  const double *ptr() const {
    return array_[0];
  }

  double *ptr() {
    return array_[0];
  }

  double get(index_t r, index_t c) const {
    DEBUG_BOUNDS(r, t_rows);
    DEBUG_BOUNDS(c, t_cols);
    return array_[c][r];
  }

  void set(index_t r, index_t c, double v) {
    DEBUG_BOUNDS(r, t_rows);
    DEBUG_BOUNDS(c, t_cols);
    array_[c][r] = v;
  }

  double &ref(index_t r, index_t c) {
    DEBUG_BOUNDS(r, t_rows);
    DEBUG_BOUNDS(c, t_cols);
    return array_[c][r];
  }

  index_t n_cols() const {
    return t_cols;
  }

  index_t n_rows() const {
    return t_rows;
  }

  size_t n_elements() const {
    // TODO: putting the size_t on the outside may be faster (32-bit
    // versus 64-bit multiplication in cases) but is more likely to result
    // in bugs
    return size_t(t_rows) * size_t(t_cols);
  }

  double *GetColumnPtr(index_t col) {
    DEBUG_BOUNDS(col, t_cols);
    return array_[col];
  }

  const double *GetColumnPtr(index_t col) const {
    DEBUG_BOUNDS(col, t_cols);
    return array_[col];
  }
};

#endif