TMat_decl.h

Go to the documentation of this file.

// -*- C++ -*-

// PLearn (A C++ Machine Learning Library)
// Copyright (C) 1998 Pascal Vincent
// Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal
//

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
// 
//  1. Redistributions of source code must retain the above copyright
//     notice, this list of conditions and the following disclaimer.
// 
//  2. Redistributions in binary form must reproduce the above copyright
//     notice, this list of conditions and the following disclaimer in the
//     documentation and/or other materials provided with the distribution.
// 
//  3. The name of the authors may not be used to endorse or promote
//     products derived from this software without specific prior written
//     permission.
// 
// THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
// NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
// 
// This file is part of the PLearn library. For more information on the PLearn
// library, go to the PLearn Web site at www.plearn.org


 

/* *******************************************************      
 * $Id: TMat_decl.h,v 1.4 2004/07/22 17:19:47 lapalmej Exp $
 * AUTHORS: Pascal Vincent & Yoshua Bengio
 * This file is part of the PLearn library.
 ******************************************************* */


#ifndef TMat_decl_INC
#define TMat_decl_INC

#include "TVec_impl.h"

namespace PLearn {
using namespace std;

// predeclarations
template<class T> class TMatElementIterator;
template<class T> class TMatRowsIterator;
template<class T> class TMatColRowsIterator;
template<class T> class TMatRowsAsArraysIterator;


template <class T>
class TMat
{
  friend class TVec<T>;
  friend class Variable; 
  friend class VarArray; 

protected:
  int offset_; 
  int mod_; 
  int length_; 
  int width_; 
  PP< Storage<T> > storage; 
public:

  int nrows() const { return length_; }
  int ncols() const { return width_; }

public:

  typedef T value_type;
  typedef int size_type;
  typedef TMatElementIterator<T> iterator; // iterator over elements
  typedef TMatElementIterator<T> const_iterator; // iterator over elements
  typedef T* compact_iterator; // super-efficient iterator over elements but works reliably only for compact matrices
  typedef T* rowelements_iterator; // iterator over elements of a
  // particular row
  typedef TMatRowsIterator<T> rows_iterator;
  typedef TMatRowsAsArraysIterator<T> rows_as_arrays_iterator;
  typedef TMatColRowsIterator<T> colrows_iterator;

  TMat<T>()
    :offset_(0), mod_(0), length_(0), width_(0)
    {}

  TMat<T>(int the_length, int the_width)
    :offset_(0), mod_(0), length_(0), width_(0)
    { resize(the_length, the_width); }

  TMat<T>(int the_length, int the_width, const T& init_value)
    :offset_(0), mod_(0), length_(0), width_(0)
    { 
      resize(the_length, the_width); 
      fill(init_value);
    }

  TMat<T>(int the_length, int the_width, T* the_data)
    :offset_(0), mod_(the_width), length_(the_length), width_(the_width), 
     storage(new Storage<T>(the_length*the_width, the_data))
    {}

  TMat<T>(int the_length, int the_width, const TVec<T>& v);

  inline const TMat<T>& operator=(const TMat<T>& other)
    {
      storage = other.storage;
      offset_ = other.offset_;
      mod_ = other.mod_;
      length_ = other.length_;
      width_ = other.width_;
      return *this;
    }

  inline iterator begin() const;
  inline iterator end() const;

  inline compact_iterator compact_begin() const
    { 
#ifdef BOUNDCHECK
      if(mod()!=width()) 
        PLERROR("You cannot use a compact iterator to iterate over the elements of a non compact matrix");
#endif
      return data();
    }

  inline compact_iterator compact_end() const
    { return data()+size(); }

  inline rowelements_iterator rowelements_begin(int rownum) const
    {
#ifdef BOUNDCHECK
      if(rownum<0 || rownum>=length())
        PLERROR("OUT OF RANGE rownum in rowelements_begin");
#endif
      return data()+rownum*mod();
    }

  inline rowelements_iterator rowelements_end(int rownum) const
    { return data()+rownum*mod()+width(); }

  TMatRowsIterator<T> rows_begin();
  TMatRowsIterator<T> rows_end();

  TMatRowsAsArraysIterator<T> rows_as_arrays_begin();
  TMatRowsAsArraysIterator<T> rows_as_arrays_end();


  TMatColRowsIterator<T> col_begin(int column);

  TMatColRowsIterator<T> col_end(int column);
  

  void resize(int newlength, int newwidth, int extrabytes=0)
    {
#ifdef BOUNDCHECK
      if(newlength<0 || newwidth<0)
        PLERROR("IN TMat::resize(int newlength, int newwidth)\nInvalid arguments (<0)");
#endif
      if (newlength==length_ && newwidth==width_) return;
      if(storage.isNull())
      {
        offset_ = 0;
        length_ = newlength;
        width_ = newwidth;
        mod_ = newwidth;
        storage = new Storage<T>(length()*mod());
      }
      else
      {
        if(storage->usage()>1 && newwidth > mod()-offset_%mod())
          PLERROR("IN TMat::resize(int newlength, int newwidth)\nFor safety reasons, increasing the width() beyond mod() - offset_%mod() is not allowed when the storage is shared with others");

        int newsize = offset_+newlength*MAX(mod(),newwidth);
        if(offset_+newsize>storage->length())
          storage->resize(offset_+newsize+extrabytes);
        length_ = newlength;
        width_ = newwidth;
        if(newwidth>mod())
          mod_ = newwidth;
      }
    }

  inline int length() const
    { return length_; }

  inline int width() const
    { return width_; }

  inline int size() const
    { return length_*width_; }

  inline int mod() const
    { return mod_; }

  inline PP< Storage<T> > getStorage() const 
    { return storage; }

  inline bool isSquare() const
    { return length() == width(); }

  bool hasMissing() const
    {
      iterator it = begin();
      iterator itend = end();
      for(; it!=itend; ++it)
        if(is_missing(*it))
          return true;
      return false;
    }

  inline T* data() const
    {
#ifdef BOUNDCHECK
      if(storage.isNull())
        PLERROR("IN TMat::data()\nAttempted to get a pointer to the data of an empty matrix");
#endif
      return storage->data+offset_; 
    }

  inline T* operator[](int rownum) const
    {
#ifdef BOUNDCHECK
      if(rownum<0 || rownum>=length())
        PLERROR("OUT OF BOUND ACCESS IN TMat::operator[](int rownum=%d), length=%d",rownum,length());
#endif
      return storage->data + offset_ + mod()*rownum; 
    }

  inline T* rowdata(int i) const { return (*this)[i]; }

  inline T& operator()(int rownum, int colnum) const
    {
#ifdef BOUNDCHECK
      if(rownum<0 || rownum>=length() || colnum<0 || colnum>=width())
        PLERROR("OUT OF BOUND ACCESS IN TMat::operator()(int rownum, int colnum)"
                " width=%d; length=%d; colnum=%d; rownum=%d;", width(), length(), colnum, rownum);
#endif
      return storage->data[offset_ + mod()*rownum + colnum];
    }

  inline TVec<T> operator()(int rownum) const
  {
#ifdef BOUNDCHECK
    if(rownum<0 || rownum>=length())
      PLERROR("OUT OF BOUND ACCESS IN TMat_impl::operator()(int rownum)");
#endif
    TVec<T> tv;
    tv.length_ = width();
    tv.offset_ = offset_ + mod()*rownum;
    tv.storage = storage;
    return tv;
  }

  void write(PStream& out) const
    {
      T* ptr = 0;
      if(storage)
        ptr = data();

      switch(out.outmode)
      {
      case PStream::raw_ascii:      
      case PStream::pretty_ascii:
        for(int i=0; i<length_; i++, ptr+=mod_)
        {
          for(int j=0; j<width_; j++)
          {
            out << ptr[j];
            out.put('\t');
          }
          out.put('\n');
        }
        break;
        
      case PStream::raw_binary:
        for(int i=0; i<length_; i++, ptr+=mod_)
          binwrite_(out, ptr, width_);
        break;
        
      case PStream::plearn_ascii:
      {
        if(!out.implicit_storage)
        {
          out.write("TMat("); 
          out << length_ << width_ << mod_ << offset_ << storage;
          out.write(")\n");
        }
        else // implicit storage
        {
          out << length_;
          out.put(' ');
          out << width_;
          out.write(" [ \n");
          for(int i=0; i<length_; i++, ptr+=mod_)
          {
            for(int j=0; j<width_; j++)
            {
              out << ptr[j];
              out.put('\t');
            }
            out.put('\n');
          }
          out.write("]\n");
        }
      }
      break;

      case PStream::plearn_binary:
      {
        if(!out.implicit_storage)
        {
          out.write("TMat("); 
          out << length_ << width_ << mod_ << offset_ << storage;
          out.write(")\n");
        }
        else // implicit storage
        {
          unsigned char typecode;
          if(byte_order()==LITTLE_ENDIAN_ORDER)
          {
            out.put(0x14); // 2D little-endian 
            typecode = TypeTraits<T>::little_endian_typecode();
          }
          else
          {
            out.put(0x15); // 2D big-endian
            typecode = TypeTraits<T>::big_endian_typecode();
          }
              
          // write typecode
          out.put(typecode);
              
          // write length and width in raw_binary 
          out.write((char*)&length_, sizeof(length_));
          out.write((char*)&width_, sizeof(width_));
              
          // write the data
          for(int i=0; i<length_; i++, ptr+=mod_)
            binwrite_(out, ptr, width_);
        }
      }
      break;
      
      default:
        PLERROR("In TMat::write(PStream& out)  unknown outmode!!!!!!!!!");
        break;
      }
    }



  void read(PStream& in)
    {
      switch(in.inmode)
      {
      case PStream::raw_ascii:
      case PStream::raw_binary:
      {
        T* ptr = data();
        for(int i=0; i<length_; i++, ptr+=mod_)
          for(int j=0; j<width_; j++)
            in >> ptr[j];
      }
      break;

      case PStream::plearn_ascii:
      case PStream::plearn_binary:
      {
        in.skipBlanksAndComments();
        int c = in.peek();
        if(c=='T') // explicit storage
        {
          char word[6];
          // !!!! BUG: For some reason, this hangs!!!
          // in.read(word,5);

          for(int i=0; i<5; i++)
            in.get(word[i]);

          word[5]='\0';
          if(strcmp(word,"TMat(")!=0)
            PLERROR("In operator>>(PStream&, TMat&) '%s' not a proper header for a TMat!",word);
          // v.storage = 0;
          in >> length_ >> width_ >> mod_ >> offset_;
          in >> storage;
          in.skipBlanksAndCommentsAndSeparators();
          int c = in.get(); // skip ')'
          if(c!=')')
            PLERROR("In operator>>(PStream&, TMat&) expected a closing parenthesis, found '%c'",c);
        }
        else // implicit storage
        {
          if(isdigit(c)) // ascii mode with length and width given  
          {
            int l,w;
            in >> l >> w;
            in.skipBlanksAndComments();
            c = in.get();
            if(c!='[')
              PLERROR("Error in TMat::read(PStream& in), expected '[', read '%c'",c);
            in.skipBlanksAndCommentsAndSeparators();
            resize(l,w);
            T* ptr = (l>0 && w>0)? data():0;
            for(int i=0; i<length_; i++, ptr+=mod_)
              for(int j=0; j<width_; j++)
              {
                in.skipBlanksAndCommentsAndSeparators();
                in >> ptr[j];
              }
            in.skipBlanksAndCommentsAndSeparators();
            c = in.get();
            if(c!=']')
              PLERROR("Error in TMat::read(PStream& in), expected ']', read '%c'",c);
          }
          else if(c==0x14 || c==0x15) // it's a binary 2D sequence
          {
            in.get(); // eat c
            unsigned char typecode = in.get(); 
            int l, w;                  
            in.read((char*)&l,sizeof(l));
            in.read((char*)&w,sizeof(w));
            bool inverted_byte_order = ((c==0x14 && byte_order()==BIG_ENDIAN_ORDER) 
                                        || (c==0x15 && byte_order()==LITTLE_ENDIAN_ORDER) );
            if(inverted_byte_order)
            {
              endianswap(&l);
              endianswap(&w);
            }
            resize(l,w);
            T* ptr = data();
            for(int i=0; i<length_; i++, ptr+=mod_)                    
              binread_(in, ptr, width_, typecode);
          }
          else
            PLERROR("In TMat::read(PStream& in) Char with ascii code %d not a proper first character in the header of a TMat!",c);
        }
      }
      break;
      
      default:
        PLERROR("In TMat<T>::read(PStream& in)  unknown inmode!!!!!!!!!");
        break;
      }
    }

  // The following methods are deprecated, and just call corresponding functions.
  // Please call those functions directly in new code
  //void write(ostream& out) const { PLearn::write(out, *this); }
  //void read(istream& in) { PLearn::read(in, *this); }
  void save(const string& filename) const { savePMat(filename, *this); }
  void load(const string& filename) { loadPMat(filename, *this); }
    
  inline TMat<T> column(int colnum) const
    { return subMatColumns(colnum, 1); }

  inline TMat<T> firstColumn() const
    { return column(0); }

  inline TMat<T> lastColumn() const
    { return column(width()-1); }

  inline TMat<T> row(int row) const
    { return subMatRows(row, 1); }

  inline T& firstElement() const { return *data(); }
  inline T& lastElement() const { return operator()(length-1,width-1); }

  inline TVec<T> firstRow() const { return operator()(0); } 
  inline TVec<T> lastRow() const { return operator()(length_ - 1); }
  inline TVec<T> front() const { return firstRow(); }
  inline TVec<T> back() const { return lastRow(); }
  
  template<class I>
  inline TMat<T> columns(const TVec<I>& columns) const
    {
      TMat<T> result(length(),columns.length());
      selectColumns(*this,columns,result);
      return result;
    }

  template<class I>
  inline TMat<T> rows(const TVec<I>& rows) const
    {
      TMat<T> result(rows.length(),width());
      selectRows(*this,rows,result);
      return result;
    }

  inline bool operator==(const TMat<T>& other) const
    {
      if(length() != other.length() || width() != other.width())
        return false;

      iterator it = begin();
      iterator end = end();
      iterator otherIt = other.begin();
      for(; it != end; ++it)
        if(*it == *otherIt)
          ++otherIt;
        else
          return false;

      return true;
    }

  template<class I>
  inline TMat<T> operator()(const TVec<I>& rows, const TVec<I>& columns) const
    {
      TMat<T> result(rows.length(),columns.length());
      select(*this,rows,columns,result);
      return result;
    }

  inline TMat<T> subMat(int rowstart, int colstart, int newlength, int newwidth) const
    {
#ifdef BOUNDCHECK
      if(rowstart<0 || newlength<0 || rowstart+newlength>length()
         || colstart<0 || newwidth<0 || colstart+newwidth>width())
        PLERROR("Mat::subMat(int rowstart, int colstart, int newlength, int newwidth) OUT OF BOUNDS"
                "  rowstart=%d colstart=%d newlength=%d newwidth=%d length()=%d width()=%d",
                rowstart, colstart, newlength, newwidth, length(), width());
#endif
      TMat<T> subm = *this;
      subm.length_ = newlength;
      subm.width_ = newwidth;
      subm.offset_ += rowstart*mod() + colstart;
      return subm;
    }

  inline TMat<T> subMatRows(int rowstart, int newlength) const
    {
#ifdef BOUNDCHECK
      if(rowstart<0 || newlength<0 || rowstart+newlength>length())
        PLERROR("TMat::subMatRows(int rowstart, int newlength) OUT OF BOUNDS");
#endif
      TMat<T> subm = *this;
      subm.length_ = newlength;
      subm.offset_ += rowstart*mod();
      return subm;
    }

  inline TMat<T> subMatColumns(int colstart, int newwidth) const
    {
#ifdef BOUNDCHECK
      if(colstart<0 || newwidth<0 || colstart+newwidth>width())
        PLERROR("Mat::subMatColumns(int colstart, int newwidth) OUT OF BOUNDS");
#endif
      TMat<T> subm = *this;
      subm.width_ = newwidth;
      subm.offset_ += colstart;
      return subm;
    }

  TMat<T> copy() const
    {
      TMat<T> freshcopy(length(),width());
      freshcopy << *this;
      return freshcopy;
    }

  void copyTo(T* x) const
    {
      T* row = data(); // get data start
      int k=0;
      for(int i=0; i<length(); i++,row+=mod())
        for (int j=0;j<width();j++,k++)
          x[k] = row[j];
    }

  void makeDeepCopyFromShallowCopy(map<const void*, void*>& copies);

  TMat<T> deepCopy(map<const void*, void*>& copies) const;


  TVec<T> toVecCopy() const;

  TVec<T> toVec() const;

  bool isNull() const 
    { return storage.isNull(); }

  bool isNotNull() const
    { return storage.isNotNull(); }

  bool isEmpty() const
    { return length_==0; }

  bool isNotEmpty() const
    { return length_!=0; }

  /*
  inline operator bool() const
  { return isNotEmpty(); }
  */

  inline bool operator!() const
    { return isEmpty(); }

  void fill(const T& value) const
    {
      if (isNotNull()) {
        if(isCompact())
          fill_n(data(),size(),value); 
        else
        {
          int l = length();
          T* ptr = data();
          while(l--)
          {
            fill_n(ptr, width(), value);
            ptr += mod();
          }
        }
      }
    }
  
  inline void operator=(const T& f) const
    { fill(f); }

  inline void clear() const
    { 
      if(isCompact())
        clear_n(data(),size()); 
      else
      {
        int l = length();
        T* ptr = data();
        while(l--)
        {
          clear_n(ptr, width());
          ptr += mod();
        }
      }
    }

  void swapRows(int i, int j) const
    {
      if(i!=j)
      {
        //T* Mi = rowdata(i);
        //T* Mj = rowdata(j);
        T* Mi = (*this)[i];
        T* Mj = (*this)[j];
        for (int k=0;k<width();k++)
        {
          T tmp = Mi[k];
          Mi[k] = Mj[k];
          Mj[k] = tmp;
        }
      }
    }

  int findRow(const TVec<T>& row) const;

  inline void appendRow(const TVec<T>& newrow);

  inline void push_back(const TVec<T>& newrow) { appendRow(newrow); }
  inline void pop_back() { length_ -= 1; }


  bool isCompact() const
    { return mod() == width(); }

  bool isSymmetric() const 
    {
      if (!isSquare())
        return false;

      if (length() == 0)
      {
        PLWARNING("at bool TMat::isSymmetric(), the size of the matrix is 0\n");
        return false;
      }

      for (int i = 0; i < length() - 1 ; i++)
        for (int j = i + 1; j < width(); j++)
          if ( (*this)[i][j] != (*this)[j][i] )
            return false;

      return true;
    }

  void compact()
    {
      if(storage->length() != length()*width())
      {
        if(storage->usage()>1)
          PLERROR("IN Mat::compact() compact operation not allowed when matrix storage is shared (for obvious reasons)");
        operator=(copy());
      }
    }

  void transpose()
    {
      if(mod()!=width())
        PLERROR("In transpose() can transpose in place only compact matrices whose mod()==width() (that is often not the case for submatrices");
      for(int i=0; i<length(); i++)
      {
        //T* rowi = rowdata(i);
        T* rowi = (*this)[i];
        T* colielem = rowi+i+mod();
        for(int j=i+1; j<width(); j++, colielem+=mod())
          swap(rowi[j], *colielem);
      }
    }

  void swapUpsideDown() const
    {
      int half = length()/2;
      for(int i=0; i<half; i++)
        swapRows(i, length()-i-1);
    }

  void print(ostream& out = cout) const;
  void input(istream& in = cin) const;

  // calls print with cerr, usefull with gdb (> call obj.debugprint() )
  void debugPrint(){print(cerr);}


  inline void operator<<(const string& datastring) const
    { 
      istrstream in(datastring.c_str());
      input(in); 
    }

};

typedef TMat<real> Mat;


// Type traits (especially type "names" for displaying optionHelp() )

template<class T>
class TypeTraits< TMat<T> >
{
public:
  static inline string name()
    { return string("TMat< ") + TypeTraits<T>::name()+" >"; }

  static inline unsigned char little_endian_typecode()
    { return 0xFF; }

  static inline unsigned char big_endian_typecode()
    { return 0xFF; }
};

} // end of namespace PLearn
#endif

---