RowMapSparseMatrix.h

Go to the documentation of this file.

// -*- C++ -*-

// PLearn (A C++ Machine Learning Library)
// Copyright (C) 1999-2002 Yoshua Bengio
//

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

#ifndef ROWMAPSPARSEMATRIX
#define ROWMAPSPARSEMATRIX

#include <map>
#include "Mat.h"
#include "TMat.h"
#include "SparseMatrix.h"

namespace PLearn {
using namespace std;


  template <class T>
  class RowMapSparseMatrix {
    public:
    vector< map<int,T> > rows;
    protected:
    int _width;

    public:
    bool save_binary;

    RowMapSparseMatrix(int n_rows=0, int n_columns=0) 
      : rows(n_rows), _width(n_columns), save_binary(true) {}

    RowMapSparseMatrix(string filename) { load(filename); }

    RowMapSparseMatrix(const Mat& m, bool fill_all=true) : rows(m.length()), _width(m.width()), save_binary(true)
    {
      if (fill_all)
        for (int i=0;i<length();i++)
        {
          real* r=m[i];
          map<int,T>& row_i=rows[i];
          for (int j=0;j<width();j++)
            row_i[j]=T(r[j]);
        }
      else
        for (int i=0;i<length();i++)
        {
          real* r=m[i];
          map<int,T>& row_i=rows[i];
          for (int j=0;j<width();j++){
            if(T(r[j])!=0)
              row_i[j]=T(r[j]);
          }
        }
    }

    RowMapSparseMatrix(const SparseMatrix& sm, int n_rows, int n_cols) : rows(n_rows), _width(n_cols), save_binary(false) {

      for (int j = 0; j < n_cols; j++) {
        int bcol_j = (int)sm.beginRow[j];
        int ecol_j = (int)sm.endRow[j];
        for (int row_i = bcol_j; row_i <= ecol_j; row_i++) {
          int i = (int)sm.row[row_i];
          (*this)(i, j) = (T)sm.values[row_i];
        }
      }

    }

    void resize(int n_rows, int n_columns) {
      rows.resize(n_rows);
      _width=n_columns;
      clear();
    }
      
    void clear() {
      int s=rows.size();
      for (int i=0;i<s;i++)
        rows[i].clear();
    }

    void clearRow(int i)
    {
      rows[i].clear();
    }

    T& operator()(int i, int j) { 
#ifdef BOUNDCHECK      
      if (i<0 || i>=length() && j<0 || j>=width())
        PLERROR("RowMapSparseMatrix: out-of-bound access to (%d,%d), dims=(%d,%d)",
              i,j,length(),width());
#endif
      return rows[i][j]; 
    }

    const T& operator()(int i, int j) const { 
#ifdef BOUNDCHECK      
      if (i<0 || i>=length() && j<0 || j>=width())
        PLERROR("RowMapSparseMatrix: out-of-bound access to (%d,%d), dims=(%d,%d)",
              i,j,length(),width());
#endif
      map<int,T>& row_i = rows[i];
      typename map<int,T>::const_iterator it = row_i.find(j);
      if (it==row_i.end())
        return 0;
      return it->second;
    }

    map<int,T>& operator()(int i) { return rows[i]; }

    int size() const { 
      int _size=0;
      for (unsigned int i=0;i<rows.size();i++)
        _size += rows[i].size();
      return _size;
    }
    int length() const { return rows.size(); }
    int width() const { return _width; }

    void save(string filename) const { 
      ofstream out(filename.c_str()); write(out); 
    }
    void saveAscii(string filename) { 
      bool old_mode=save_binary;
      save_binary = false; save(filename); 
      save_binary=old_mode;
    }
    void load(string filename) { ifstream in(filename.c_str()); read(in); }
    void read(istream& in) {
      readHeader(in,"RowMapSparseMatrix");
      int n_bytes=0;
      readField(in,"sizeof(T)",n_bytes,4);
      readField(in,"save_binary",save_binary);
      if (n_bytes!=sizeof(T) && save_binary)
        PLERROR("RowMapSparseMatrix<T>::read, input matrix has sizeof(T)=%d, expected %d",n_bytes,sizeof(T));
      int n_rows;
      readField(in,"n_rows",n_rows);
      rows.resize(n_rows);
      readField(in,"width",_width);
      for (int i=0;i<n_rows;i++)
      {
        map<int,T>& row_i = rows[i];
        int n_elements;
        if (save_binary)
        {
          binread(in,n_elements);
          for (int j=0;j<n_elements;j++)
          {
            int col;
            T value;
            binread(in,col);
            binread(in,value);
            row_i[col]=value;
          }
        } else
        {
          PLearn::read(in,n_elements);
          for (int j=0;j<n_elements;j++)
          {
            int col;
            T value;
            PLearn::read(in,col);
            PLearn::read(in,value);
            row_i[col]=value;
          }
        }
      }
      readFooter(in,"RowMapSparseMatrix");
    }

    void write(ostream& out) const {
      writeHeader(out,"RowMapSparseMatrix");
      int sizeofT = sizeof(T);
      writeField(out,"sizeof(T)",sizeofT);
      writeField(out,"save_binary",save_binary);
      writeField(out,"n_rows",rows.size());
      writeField(out,"width",_width);
      int s=rows.size();
      for (int i=0;i<s;i++)
      {
        const map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        int n_elements = row_i.size();
        if (save_binary)
        {
          binwrite(out,n_elements);
          for (;it!=end;++it)
          {
            binwrite(out,it->first);
            binwrite(out,it->second);
          }
        }
        else
        {
          PLearn::write(out,n_elements);
          for (;it!=end;++it)
          {
            PLearn::write(out,it->first);
            PLearn::write(out,it->second);
          }
        }
      }
      writeFooter(out,"RowMapSparseMatrix");
    }

    void saveNonZeroElements(string filename) const {
      ofstream out(filename.c_str());
      out << length() << " " << width() << " " << size() << endl;
      int s=rows.size();
      for (int i=0;i<s;i++)
      {
        const map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          out << i << " " << it->first << " " << it->second << endl;
      }
    }

    void product(const Vec& x, Vec& y) {
      real* _y=y.data();
      real* _x=x.data();
#ifdef BOUNDCHECK
      if (y.length()!=length() || x.length()!=width())
        PLERROR("RowMapSparseMatrix::product: inconsistent arguments (%d,%d), dims=(%d,%d)",
              x.length(),y.length(),length(),width());
#endif
      int s=rows.size();
      for (int i=0;i<s;i++)
      {
        real res=0;
        map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          res += it->second * _x[it->first];
        _y[i] = res;
      }
    }

    Mat toMat()
    {
      Mat res(length(),width());
      for (int i=0;i<length();i++)
      {
        map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        real* res_i=res[i];
        for (;it!=end;++it)
          res_i[it->first] = it->second;
      }
      return res;
    }

    bool isSymmetric(real tolerance=0)
    {
      for (int i=0;i<length();i++)
      {
        map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          if (it->first!=i)
          {
            T& other_guy = (*this)(it->first,i);
            if (fabs(other_guy - it->second)>tolerance)
              return false;
          }
      }
      return true;
    }

    bool fillSymmetricPart(real tolerance=0)
    {
      for (int i=0;i<length();i++)
      {
        map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
        {
          int j=it->first;
          if (j!=i)
          {
            map<int,T>& row_j= rows[j];
            typename map<int,T>::iterator symm_it =row_j.find(i);
            if (symm_it==row_j.end())
              row_j[i]=it->second;
            else
              if (fabs(symm_it->second - it->second)>tolerance)
                return false;
          }
        }
      }
      return true;
    }

    void diag(Vec& d)
    {
      real* d_=d.data();
      for (int i=0;i<length();i++)
        d_[i] = (*this)(i,i);
    }

    void diagonalOfSquare(Vec& d)
    {
      real* d_=d.data();
      for (int i=0;i<length();i++)
      {
        real sum2=0;
        map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          sum2 += it->second * it->second;
        d_[i] = sum2;
      }
    }

    real dotRow(int i, Vec v)
    {
#ifdef BOUNDCHECK
      if (v.length()!=width())
        PLERROR("RowMapSparseMatrix::dotRow(%d,v), v.length_=%d != matrix width=%d",
              i,v.length(),width());
#endif
      real s = 0;
      real* v_=v.data();
      map<int,T>& row_i = rows[i];
      typename map<int,T>::const_iterator it = row_i.begin();
      typename map<int,T>::const_iterator end = row_i.end();
      for (;it!=end;++it)
        s += it->second * v_[it->first];
      return s;
    }

    real dotColumn(int j, Vec v)
    {
#ifdef BOUNDCHECK
      if (v.length()!=length())
        PLERROR("RowMapSparseMatrix::dotColumn(%d,v), v.length_=%d != matrix length=%d",
              j,v.length(),length());
#endif
      PLWARNING("RowMapSparseMatrix is not appropriate to perform dotColumn operations");
      real s=0;
      real* v_=v.data();
      for (int i=0;i<v.length();v++)
        s += (*this)(i,j) * v_[i];
      return s;
    }

    void transposeProduct(RowMapSparseMatrix& m, bool verbose = false)
    {
      RowMapSparseMatrix<T>& self = (*this);
      int n = self.length();
      RowMapSparseMatrix<T> mt(n, n);
      transpose(m, mt);
      int nnz = 0;
      for (int i = 0; i < n; i++)
      {
        if (verbose)
        {
          if (i % 10 == 0 && i != 0) { 
            cout << "[" << i << "]" << " ";
            if (i % 100 == 0) cout << endl;
            else cout.flush();
          }
        }
        for (int j = 0; j < n; j++)
        {
          T val =  multiplyVecs(mt(i), mt(j));
          if (val != 0)
          {
            nnz++;
            self(i, j) = val;
          }
        }
      }
      if (verbose) 
        cout << endl;
    }

    void add2Rows(Vec row, bool only_on_non_zeros=true)
    {
      if (!only_on_non_zeros)
        PLERROR("RowMapSparseMatrix::add2Rows(Vec,bool) works only with bool=true");
      real* r=row.data();
      for (int i=0;i<length();i++)
      {
        map<int,T>& row_i = rows[i];
        typename map<int,T>::iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          it->second += r[it->first];
      }
    }
    void add2Columns(Vec col, bool only_on_non_zeros=true)
    {
      if (!only_on_non_zeros)
        PLERROR("RowMapSparseMatrix::add2Columns(Vec,bool) works only with bool=true");
      for (int i=0;i<length();i++)
      {
        real col_i=col[i];
        map<int,T>& row_i = rows[i];
        typename map<int,T>::iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          it->second += col_i;
      }
    }
    void add(real scalar, bool only_on_non_zeros=true)
    {
      if (!only_on_non_zeros)
        PLERROR("RowMapSparseMatrix::add(real,bool) works only with bool=true");
      for (int i=0;i<length();i++)
      {
        map<int,T>& row_i = rows[i];
        typename map<int,T>::iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          it->second += scalar;
      }
    }

    void averageAcrossRowsAndColumns(Vec avg_across_rows, Vec avg_across_columns, 
                                     bool only_on_non_zeros=true)
    {
      if (!only_on_non_zeros)
        PLERROR("RowMapSparseMatrix::averageAcrossRowsAndColumns works only with only_on_non_zeros=true");
      avg_across_rows.resize(width());
      avg_across_columns.resize(length());
      avg_across_rows.clear();
      avg_across_columns.clear();
      TVec<int> column_counts(width());
      for (int i=0;i<length();i++)
      {
        real& avg_cols_i=avg_across_columns[i];
        real* avg_rows = avg_across_rows.data();
        map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        int n=0;
        for (;it!=end;++it)
        {
          avg_cols_i += it->second;
          int j=it->first;
          avg_rows[j] += it->second;
          n++;
          column_counts[j]++;
        }
        avg_cols_i /= n;
      }
      for (int j=0;j<width();j++)
        avg_across_rows[j] /= column_counts[j];
    }

      real sumRow(int i)
      {
        real s=0;
        map<int,T>& row_i = rows[i];
        typename map<int,T>::const_iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          s += it->second;
        return s;
      }


    void operator*=(real scalar)
    {
      for (int i=0;i<length();i++)
      {
        map<int,T>& row_i = rows[i];
        typename map<int,T>::iterator it = row_i.begin();
        typename map<int,T>::const_iterator end = row_i.end();
        for (;it!=end;++it)
          it->second *= scalar;
      }
    }
      
    static void transpose(RowMapSparseMatrix<T>& src, RowMapSparseMatrix<T>& dest)
    {
      dest.clear();
      for (int i = 0; i < src.length(); i++)
      {
        map<int, T>& row = src(i);
        typename map<int, T>::iterator beg = row.begin();
        typename map<int, T>::iterator end = row.end();
        while (beg != end)
        {
          int col = beg->first;
          T val = beg->second;
          dest(col, i) = val;
          ++beg;
        }
      }
    }
      
    static real multiplyVecs(map<int, T>& map1, map<int, T>& map2) 
    {
      if (map1.size() == 0 || map2.size() == 0)
        return 0;
      typename map<int, T>::iterator beg1 = map1.begin();
      typename map<int, T>::iterator beg2 = map2.begin();
      typename map<int, T>::iterator end1 = map1.end();
      typename map<int, T>::iterator end2 = map2.end();
      int col1, col2;
      T val1, val2, sum = 0;
      bool fend1 = (beg1 == end1), fend2 = (beg2 == end2);
      int OUT = getMaxColumnIndex(map1, map2) + 1;
      
      while (!fend1 || !fend2) {
        if (!fend1)
          col1 = beg1->first;
        else
          col1 = OUT;
        if (!fend2)
          col2 = beg2->first;
        else
          col2 = OUT;
        val1 = beg1->second;
        val2 = beg2->second;
        if (col1 == col2) {
          sum += (val1 * val2);
          ++beg1;
          if (beg1 == end1) fend1 = true;
          ++beg2;
          if (beg2 == end2) fend2 = true;
        } else if (col1 < col2) {
          //sum += (val1 * val2);
          ++beg1;
          if (beg1 == end1) fend1 = true;
        } else if (col1 > col2) {
          //sum += (val1 * val2);
          ++beg2;
          if (beg2 == end2) fend2 = true;
        }
      }
      return sum;
    }

    static real euclidianDistance(map<int, real>& map1, map<int, real>& map2) 
    {
      if (map1.size() == 0 || map2.size() == 0)
        return 0;
      map<int, real>::iterator beg1 = map1.begin();
      map<int, real>::iterator beg2 = map2.begin();
      map<int, real>::iterator end1 = map1.end();
      map<int, real>::iterator end2 = map2.end();
      int col1, col2;
      real val1, val2, diff, sum = 0;
      bool fend1 = (beg1 == end1), fend2 = (beg2 == end2);
      int OUT = getMaxColumnIndex(map1, map2) + 1;
        
      while (!fend1 || !fend2) 
      {
        if (!fend1)
          col1 = beg1->first;
        else
          col1 = OUT;
        if (!fend2)
          col2 = beg2->first;
        else
          col2 = OUT;
        val1 = beg1->second;
        val2 = beg2->second;
        if (col1 == col2) 
        {
          diff = val1 - val2;
          sum += (diff * diff);
          beg1++;
          if (beg1 == end1) fend1 = true;
          beg2++;
          if (beg2 == end2) fend2 = true;
        } else if (col1 < col2) 
        {
          diff = val1;
          sum += (diff * diff);
          beg1++;
          if (beg1 == end1) fend1 = true;
        } else if (col1 > col2) 
        {
          diff = val2;
          sum += (diff * diff);
          beg2++;
          if (beg2 == end2) fend2 = true;
        }
      }
      //return sqrt(sum);
      return sum;
    }

    static void substractVecs(map<int, real>& map1, map<int, real>& map2, Vec& dest) 
    {
      map<int, real>::iterator beg1 = map1.begin();
      map<int, real>::iterator beg2 = map2.begin();
      map<int, real>::iterator end1 = map1.end();
      map<int, real>::iterator end2 = map2.end();
      int col1, col2;
      real val1, val2;
      bool fend1 = (beg1 == end1), fend2 = (beg2 == end2);
      int OUT = getMaxColumnIndex(map1, map2) + 1;
      
      while (!fend1 || !fend2) 
      {
        if (!fend1)
          col1 = beg1->first;
        else
          col1 = OUT;
        if (!fend2)
          col2 = beg2->first;
        else
          col2 = OUT;
        val1 = beg1->second;
        val2 = beg2->second;
        if (col1 == col2) 
        {
          dest[col1] = val1 - val2;
          beg1++;
          if (beg1 == end1) fend1 = true;
          beg2++;
          if (beg2 == end2) fend2 = true;
        } else if (col1 < col2) 
        {
          dest[col1] = val1;
          beg1++;
          if (beg1 == end1) fend1 = true;
        } else if (col1 > col2) 
        {
          dest[col2] = 0 - val2;
          beg2++;
          if (beg2 == end2) fend2 = true;
        }
      }
    } 
    
    static int getMaxColumnIndex(map<int, T>& map1, map<int, T>& map2) 
    {
      map<int, real>::iterator end1 = map1.end();
      map<int, real>::iterator end2 = map2.end();
      --end1;
      --end2;
      return MAX(end1->first, end2->first);
    }

    void exportToMatlabReadableFormat(string filename)
    {
      ofstream out(filename.c_str());
#ifdef USEDOUBLE
      out.precision(20);
#else
      out.precision(6);
#endif
      for (unsigned int i = 0; i < rows.size(); i++)
      {
        typename map<int, T>::iterator beg = rows[i].begin();
        typename map<int, T>::iterator end = rows[i].end();
        while (beg != end)
        {
          out << i + 1 << " " << beg->first + 1 << " " << beg->second << endl;
          beg++;
        }
      }
      int l = length();
      int w = width();
      l--; w--;
      map<int, T>& row_l=rows[l];
      if (row_l.find(w) == row_l.end())
        out << l + 1 << " " << w + 1 << " " << 0 << endl;
    }

    real density()
    {
      return size() / real(length() * width());
    }

  };

template <class T>
void product(RowMapSparseMatrix<T>& M, const Vec& x, Vec& y) { M.product(x,y); }

} // end of namespace PLearn

#endif

---