PLearn 0.1
RowMapSparseMatrix.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PLearn (A C++ Machine Learning Library)
00004 // Copyright (C) 1999-2002 Yoshua Bengio
00005 //
00006 
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 // 
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 // 
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 // 
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 // 
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00037 #ifndef ROWMAPSPARSEMATRIX
00038 #define ROWMAPSPARSEMATRIX
00039 
00040 #include <map>
00041 #include "Mat.h"
00042 #include "TMat.h"
00043 #include "SparseMatrix.h"
00044 
00045 namespace PLearn {
00046 using namespace std;
00047 
00048 
00068 template <class T>
00069 class RowMapSparseMatrix : public PPointable {
00070 public:
00071     mutable vector< map<int,T> > rows;
00072 protected:
00073     int _width;
00074 
00075 public:
00076     bool save_binary;
00077     T null_elem;
00078 
00079     RowMapSparseMatrix(int n_rows=0, int n_columns=0, T nullelem=0) 
00080         : rows(n_rows), _width(n_columns), save_binary(true), null_elem(nullelem) {}
00081 
00082 //    RowMapSparseMatrix(string filename) { load(filename); }
00083 
00084     RowMapSparseMatrix(const Mat& m, bool fill_all=true, T nullelem=0) : rows(m.length()), _width(m.width()), 
00085                                                                          save_binary(true), null_elem(nullelem)
00086     {
00087         if (fill_all)
00088             for (int i=0;i<length();i++)
00089             {
00090                 real* r=m[i];
00091                 map<int,T>& row_i=rows[i];
00092                 for (int j=0;j<width();j++)
00093                     row_i[j]=T(r[j]);
00094             }
00095         else
00096             for (int i=0;i<length();i++)
00097             {
00098                 real* r=m[i];
00099                 map<int,T>& row_i=rows[i];
00100                 for (int j=0;j<width();j++){
00101                     if(T(r[j])!=0)
00102                         row_i[j]=T(r[j]);
00103                 }
00104             }
00105     }
00106 
00108     RowMapSparseMatrix(const SparseMatrix& sm, int n_rows, int n_cols, T nullelem=0) : 
00109         rows(n_rows), _width(n_cols), save_binary(false), null_elem(nullelem) {
00110 
00111         for (int j = 0; j < n_cols; j++) {
00112             int bcol_j = (int)sm.beginRow[j];
00113             int ecol_j = (int)sm.endRow[j];
00114             for (int row_i = bcol_j; row_i <= ecol_j; row_i++) {
00115                 int i = (int)sm.row[row_i];
00116                 (*this)(i, j) = (T)sm.values[row_i];
00117             }
00118         }
00119 
00120     }
00121 
00123     void resize(int n_rows, int n_columns) {
00124         rows.resize(n_rows);
00125         _width=n_columns;
00126         clear();
00127     }
00128       
00130     void clear() {
00131         int s=rows.size();
00132         for (int i=0;i<s;i++)
00133             rows[i].clear();
00134     }
00135 
00136     void clearRow(int i)
00137     {
00138         rows[i].clear();
00139     }
00140 
00141     T& operator()(int i, int j) { 
00142 #ifdef BOUNDCHECK      
00143         if (i<0 || i>=length() || j<0 || j>=width())
00144             PLERROR("RowMapSparseMatrix: out-of-bound access to (%d,%d), dims=(%d,%d)",
00145                     i,j,length(),width());
00146 #endif
00147         map<int,T>& row_i = rows[i];
00148         typename map<int,T>::iterator it = row_i.find(j);
00149         if (it==row_i.end())
00150             return null_elem;
00151         return it->second;
00152     }
00153 
00154 
00155     const T& operator()(int i, int j) const { 
00156 #ifdef BOUNDCHECK      
00157         if (i<0 || i>=length() || j<0 || j>=width())
00158             PLERROR("RowMapSparseMatrix: out-of-bound access to (%d,%d), dims=(%d,%d)",
00159                     i,j,length(),width());
00160 #endif
00161         map<int,T>& row_i = rows[i];
00162         typename map<int,T>::const_iterator it = row_i.find(j);
00163         if (it==row_i.end())
00164             return null_elem;
00165         return it->second;
00166     }
00167 
00168     T get(int i, int j) const { return (*this)(i,j); }
00169 
00170     void set(int i, int j, T v) const { 
00171 #ifdef BOUNDCHECK      
00172         if (i<0 || i>=length() || j<0 || j>=width())
00173             PLERROR("RowMapSparseMatrix: out-of-bound access to (%d,%d), dims=(%d,%d)",
00174                     i,j,length(),width());
00175 #endif
00176         map<int,T>& row_i = rows[i];
00177         row_i[j] = v;
00178     }
00179 
00180     bool exists(int i, int j) const { 
00181 #ifdef BOUNDCHECK      
00182         if (i<0 || i>=length() || j<0 || j>=width())
00183             PLERROR("RowMapSparseMatrix: out-of-bound access to (%d,%d), dims=(%d,%d)",
00184                     i,j,length(),width());
00185 #endif
00186         map<int,T>& row_i = rows[i];
00187         typename map<int,T>::const_iterator it = row_i.find(j);
00188         return (it!=row_i.end());
00189     }
00190 
00203     map<int,T>& operator()(int i) { return rows[i]; }
00204     map<int,T>& getRow(int i) { return rows[i]; }
00205 
00207     int size() const { 
00208         int _size=0;
00209         for (unsigned int i=0;i<rows.size();i++)
00210             _size += rows[i].size();
00211         return _size;
00212     }
00213     int length() const { return rows.size(); }
00214     int width() const { return _width; }
00215 
00216     void copyRowFrom(int i, const map<int,T>& from_row, bool clear_rest=true)
00217     {
00218         map<int,T>& row = rows[i];
00219         if (clear_rest)
00220             clearRow(i);
00221         typename map<int,T>::const_iterator it = from_row.begin();
00222         typename map<int,T>::const_iterator end = from_row.end();
00223         for (;it!=end;++it)
00224             row[it->first]=it->second;
00225     }
00226 
00227 #if 0
00228     void save(string filename) const { 
00229         ofstream out(filename.c_str()); write(out); 
00230     }
00231     void saveAscii(string filename) { 
00232         bool old_mode=save_binary;
00233         save_binary = false; save(filename); 
00234         save_binary=old_mode;
00235     }
00237     void load(string filename) { ifstream in(filename.c_str()); read(in); }
00238     void read(istream& in) {
00239         readHeader(in,"RowMapSparseMatrix");
00240         int n_bytes=0;
00241         readField(in,"sizeof(T)",n_bytes,4);
00242         readField(in,"save_binary",save_binary);
00243         if (n_bytes!=sizeof(T) && save_binary)
00244             PLERROR("RowMapSparseMatrix<T>::read, input matrix has sizeof(T)=%d, expected %d",n_bytes,sizeof(T));
00245         int n_rows;
00246         readField(in,"n_rows",n_rows);
00247         rows.resize(n_rows);
00248         readField(in,"width",_width);
00249         readField(in,"null_elem",null_elem);
00250         for (int i=0;i<n_rows;i++)
00251         {
00252             map<int,T>& row_i = rows[i];
00253             int n_elements;
00254             if (save_binary)
00255             {
00256                 binread(in,n_elements);
00257                 for (int j=0;j<n_elements;j++)
00258                 {
00259                     int col;
00260                     T value;
00261                     binread(in,col);
00262                     binread(in,value);
00263                     row_i[col]=value;
00264                 }
00265             } else
00266             {
00267                 PLearn::read(in,n_elements);
00268                 for (int j=0;j<n_elements;j++)
00269                 {
00270                     int col;
00271                     T value;
00272                     PLearn::read(in,col);
00273                     PLearn::read(in,value);
00274                     row_i[col]=value;
00275                 }
00276             }
00277         }
00278         readFooter(in,"RowMapSparseMatrix");
00279     }
00280 
00282     void write(ostream& out) const {
00283         writeHeader(out,"RowMapSparseMatrix");
00284         int sizeofT = sizeof(T);
00285         writeField(out,"sizeof(T)",sizeofT);
00286         writeField(out,"save_binary",save_binary);
00287         writeField(out,"n_rows",rows.size());
00288         writeField(out,"width",_width);
00289         writeField(out,"null_elem",null_elem);
00290         int s=rows.size();
00291         for (int i=0;i<s;i++)
00292         {
00293             const map<int,T>& row_i = rows[i];
00294             typename map<int,T>::const_iterator it = row_i.begin();
00295             typename map<int,T>::const_iterator end = row_i.end();
00296             int n_elements = row_i.size();
00297             if (save_binary)
00298             {
00299                 binwrite(out,n_elements);
00300                 for (;it!=end;++it)
00301                 {
00302                     binwrite(out,it->first);
00303                     binwrite(out,it->second);
00304                 }
00305             }
00306             else
00307             {
00308                 PLearn::write(out,n_elements);
00309                 for (;it!=end;++it)
00310                 {
00311                     PLearn::write(out,it->first);
00312                     PLearn::write(out,it->second);
00313                 }
00314             }
00315         }
00316         writeFooter(out,"RowMapSparseMatrix");
00317     }
00322     void saveNonZeroElements(string filename) const {
00323         ofstream out(filename.c_str());
00324         out << length() << " " << width() << " " << size() << endl;
00325         int s=rows.size();
00326         for (int i=0;i<s;i++)
00327         {
00328             const map<int,T>& row_i = rows[i];
00329             typename map<int,T>::const_iterator it = row_i.begin();
00330             typename map<int,T>::const_iterator end = row_i.end();
00331             for (;it!=end;++it)
00332                 out << i << " " << it->first << " " << it->second << endl;
00333         }
00334     }
00335 
00336 #endif
00337 
00338 
00339     void product(const Vec& x, Vec& y) {
00340         real* _y=y.data();
00341         real* _x=x.data();
00342 #ifdef BOUNDCHECK
00343         if (y.length()!=length() || x.length()!=width())
00344             PLERROR("RowMapSparseMatrix::product: inconsistent arguments (%d,%d), dims=(%d,%d)",
00345                     x.length(),y.length(),length(),width());
00346         if (null_elem!=0)
00347             PLERROR("RowMapSparseMatrix::product works only when null_elem=0, but was %g\n",null_elem);
00348 #endif
00349         int s=rows.size();
00350         for (int i=0;i<s;i++)
00351         {
00352             real res=0;
00353             map<int,T>& row_i = rows[i];
00354             typename map<int,T>::const_iterator it = row_i.begin();
00355             typename map<int,T>::const_iterator end = row_i.end();
00356             for (;it!=end;++it)
00357                 res += it->second * _x[it->first];
00358             _y[i] = res;
00359         }
00360     }
00361 
00362     Mat toMat()
00363     {
00364         Mat res(length(),width());
00365         if (null_elem!=0) res.fill(null_elem);
00366         for (int i=0;i<length();i++)
00367         {
00368             map<int,T>& row_i = rows[i];
00369             typename map<int,T>::const_iterator it = row_i.begin();
00370             typename map<int,T>::const_iterator end = row_i.end();
00371             real* res_i=res[i];
00372             for (;it!=end;++it)
00373                 res_i[it->first] = it->second;
00374         }
00375         return res;
00376     }
00377 
00378     bool isSymmetric(real tolerance=0)
00379     {
00380         for (int i=0;i<length();i++)
00381         {
00382             map<int,T>& row_i = rows[i];
00383             typename map<int,T>::const_iterator it = row_i.begin();
00384             typename map<int,T>::const_iterator end = row_i.end();
00385             for (;it!=end;++it)
00386                 if (it->first!=i)
00387                 {
00388                     T& other_guy = (*this)(it->first,i);
00389                     if (fabs(other_guy - it->second)>tolerance)
00390                         return false;
00391                 }
00392         }
00393         return true;
00394     }
00395 
00403     bool fillSymmetricPart(real tolerance=0)
00404     {
00405         for (int i=0;i<length();i++)
00406         {
00407             map<int,T>& row_i = rows[i];
00408             typename map<int,T>::const_iterator it = row_i.begin();
00409             typename map<int,T>::const_iterator end = row_i.end();
00410             for (;it!=end;++it)
00411             {
00412                 int j=it->first;
00413                 if (j!=i)
00414                 {
00416                     map<int,T>& row_j= rows[j];
00417                     typename map<int,T>::iterator symm_it =row_j.find(i);
00418                     if (symm_it==row_j.end())
00420                         row_j[i]=it->second;
00421                     else
00423                         if (fabs(symm_it->second - it->second)>tolerance)
00424                             return false;
00425                 }
00426             }
00427         }
00428         return true;
00429     }
00430 
00432     void diag(Vec& d)
00433     {
00434         real* d_=d.data();
00435         for (int i=0;i<length();i++)
00436             d_[i] = (*this)(i,i);
00437     }
00438 
00441     void diagonalOfSquare(Vec& d)
00442     {
00443         real* d_=d.data();
00444         for (int i=0;i<length();i++)
00445         {
00446             real sum2=0;
00447             map<int,T>& row_i = rows[i];
00448             typename map<int,T>::const_iterator it = row_i.begin();
00449             typename map<int,T>::const_iterator end = row_i.end();
00450             for (;it!=end;++it)
00451                 sum2 += it->second * it->second;
00452             d_[i] = sum2;
00453         }
00454     }
00455 
00457     real dotRow(int i, Vec v)
00458     {
00459 #ifdef BOUNDCHECK
00460         if (v.length()!=width())
00461             PLERROR("RowMapSparseMatrix::dotRow(%d,v), v.length_=%d != matrix width=%d",
00462                     i,v.length(),width());
00463         if (null_elem!=0)
00464             PLERROR("RowMapSparseMatrix::dotRow works only when null_elem=0, but was %g\n",null_elem);
00465 #endif
00466         real s = 0;
00467         real* v_=v.data();
00468         map<int,T>& row_i = rows[i];
00469         typename map<int,T>::const_iterator it = row_i.begin();
00470         typename map<int,T>::const_iterator end = row_i.end();
00471         for (;it!=end;++it)
00472             s += it->second * v_[it->first];
00473         return s;
00474     }
00475 
00477     real dotColumn(int j, Vec v)
00478     {
00479 #ifdef BOUNDCHECK
00480         if (v.length()!=length())
00481             PLERROR("RowMapSparseMatrix::dotColumn(%d,v), v.length_=%d != matrix length=%d",
00482                     j,v.length(),length());
00483         if (null_elem!=0)
00484             PLERROR("RowMapSparseMatrix::dotColumn works only when null_elem=0, but was %g\n",null_elem);
00485 #endif
00486         PLWARNING("RowMapSparseMatrix is not appropriate to perform dotColumn operations");
00487         real s=0;
00488         real* v_=v.data();
00489         for (int i=0;i<v.length();i++)
00490             s += (*this)(i,j) * v_[i];
00491         return s;
00492     }
00493 
00495     void transposeProduct(RowMapSparseMatrix& m, bool verbose = false)
00496     {
00497 #ifdef BOUNDCHECK
00498         if (null_elem!=0)
00499             PLERROR("RowMapSparseMatrix::dotColumn works only when null_elem=0, but was %g\n",null_elem);
00500 #endif
00501         RowMapSparseMatrix<T>& self = (*this);
00502         int n = self.length();
00503         RowMapSparseMatrix<T> mt(n, n);
00504         transpose(m, mt);
00505         int nnz = 0;
00506         for (int i = 0; i < n; i++)
00507         {
00508             if (verbose)
00509             {
00510                 if (i % 10 == 0 && i != 0) { 
00511                     cout << "[" << i << "]" << " ";
00512                     if (i % 100 == 0) cout << endl;
00513                     else cout.flush();
00514                 }
00515             }
00516             for (int j = 0; j < n; j++)
00517             {
00518                 T val =  multiplyVecs(mt(i), mt(j));
00519                 if (val != 0)
00520                 {
00521                     nnz++;
00522                     self(i, j) = val;
00523                 }
00524             }
00525         }
00526         if (verbose) 
00527             cout << endl;
00528     }
00529 
00532     void addToRows(Vec row, bool only_on_non_zeros=true)
00533     {
00534 #ifdef BOUNDCHECK
00535         if (null_elem!=0)
00536             PLERROR("RowMapSparseMatrix::add2Rows works only when null_elem=0, but was %g\n",null_elem);
00537 #endif
00538         if (!only_on_non_zeros)
00539             PLERROR("RowMapSparseMatrix::add2Rows(Vec,bool) works only with bool=true");
00540         real* r=row.data();
00541         for (int i=0;i<length();i++)
00542         {
00543             map<int,T>& row_i = rows[i];
00544             typename map<int,T>::iterator it = row_i.begin();
00545             typename map<int,T>::const_iterator end = row_i.end();
00546             for (;it!=end;++it)
00547                 it->second += r[it->first];
00548         }
00549     }
00552     void addToColumns(Vec col, bool only_on_non_zeros=true)
00553     {
00554 #ifdef BOUNDCHECK
00555         if (null_elem!=0)
00556             PLERROR("RowMapSparseMatrix::add2Columns works only when null_elem=0, but was %g\n",null_elem);
00557 #endif
00558         if (!only_on_non_zeros)
00559             PLERROR("RowMapSparseMatrix::add2Columns(Vec,bool) works only with bool=true");
00560         for (int i=0;i<length();i++)
00561         {
00562             real col_i=col[i];
00563             map<int,T>& row_i = rows[i];
00564             typename map<int,T>::iterator it = row_i.begin();
00565             typename map<int,T>::const_iterator end = row_i.end();
00566             for (;it!=end;++it)
00567                 it->second += col_i;
00568         }
00569     }
00572     void add(real scalar, bool only_on_non_zeros=true)
00573     {
00574         if (!only_on_non_zeros)
00575             PLERROR("RowMapSparseMatrix::add(real,bool) works only with bool=true");
00576         for (int i=0;i<length();i++)
00577         {
00578             map<int,T>& row_i = rows[i];
00579             typename map<int,T>::iterator it = row_i.begin();
00580             typename map<int,T>::const_iterator end = row_i.end();
00581             for (;it!=end;++it)
00582                 it->second += scalar;
00583         }
00584     }
00585 
00586 
00592     void averageAcrossRowsAndColumns(Vec avg_across_rows, Vec avg_across_columns, 
00593                                      bool only_on_non_zeros=true)
00594     {
00595         if (!only_on_non_zeros)
00596             PLERROR("RowMapSparseMatrix::averageAcrossRowsAndColumns works only with only_on_non_zeros=true");
00597         avg_across_rows.resize(width());
00598         avg_across_columns.resize(length());
00599         avg_across_rows.clear();
00600         avg_across_columns.clear();
00601         TVec<int> column_counts(width());
00602         for (int i=0;i<length();i++)
00603         {
00604             real& avg_cols_i=avg_across_columns[i];
00605             real* avg_rows = avg_across_rows.data();
00606             map<int,T>& row_i = rows[i];
00607             typename map<int,T>::const_iterator it = row_i.begin();
00608             typename map<int,T>::const_iterator end = row_i.end();
00609             int n=0;
00610             for (;it!=end;++it)
00611             {
00612                 avg_cols_i += it->second;
00613                 int j=it->first;
00614                 avg_rows[j] += it->second;
00615                 n++;
00616                 column_counts[j]++;
00617             }
00618             if (n>0)
00619                 avg_cols_i /= n;
00620         }
00621         for (int j=0;j<width();j++)
00622             if (column_counts[j]>0)
00623                 avg_across_rows[j] /= column_counts[j];
00624     }
00625 
00631     real sumRow(int i)
00632     {
00633         real s=0;
00634         map<int,T>& row_i = rows[i];
00635         typename map<int,T>::const_iterator it = row_i.begin();
00636         typename map<int,T>::const_iterator end = row_i.end();
00637         for (;it!=end;++it)
00638             s += it->second;
00639         return s;
00640     }
00641 
00642 
00644     void operator*=(real scalar)
00645     {
00646 #ifdef BOUNDCHECK
00647         if (null_elem!=0)
00648             PLERROR("RowMapSparseMatrix::operator* works only when null_elem=0, but was %g\n",null_elem);
00649 #endif
00650         for (int i=0;i<length();i++)
00651         {
00652             map<int,T>& row_i = rows[i];
00653             typename map<int,T>::iterator it = row_i.begin();
00654             typename map<int,T>::const_iterator end = row_i.end();
00655             for (;it!=end;++it)
00656                 it->second *= scalar;
00657         }
00658     }
00659       
00661     static void transpose(RowMapSparseMatrix<T>& src, RowMapSparseMatrix<T>& dest)
00662     {
00663         dest.clear();
00664         for (int i = 0; i < src.length(); i++)
00665         {
00666             map<int, T>& row = src(i);
00667             typename map<int, T>::iterator beg = row.begin();
00668             typename map<int, T>::iterator end = row.end();
00669             while (beg != end)
00670             {
00671                 int col = beg->first;
00672                 T val = beg->second;
00673                 dest(col, i) = val;
00674                 ++beg;
00675             }
00676         }
00677     }
00678       
00680     static real multiplyVecs(map<int, T>& map1, map<int, T>& map2) 
00681     {
00682         if (map1.size() == 0 || map2.size() == 0)
00683             return 0;
00684         typename map<int, T>::iterator beg1 = map1.begin();
00685         typename map<int, T>::iterator beg2 = map2.begin();
00686         typename map<int, T>::iterator end1 = map1.end();
00687         typename map<int, T>::iterator end2 = map2.end();
00688         int col1, col2;
00689         T val1, val2, sum = 0;
00690         bool fend1 = (beg1 == end1), fend2 = (beg2 == end2);
00691         int OUT = getMaxColumnIndex(map1, map2) + 1;
00692       
00693         while (!fend1 || !fend2) {
00694             if (!fend1)
00695                 col1 = beg1->first;
00696             else
00697                 col1 = OUT;
00698             if (!fend2)
00699                 col2 = beg2->first;
00700             else
00701                 col2 = OUT;
00702             val1 = beg1->second;
00703             val2 = beg2->second;
00704             if (col1 == col2) {
00705                 sum += (val1 * val2);
00706                 ++beg1;
00707                 if (beg1 == end1) fend1 = true;
00708                 ++beg2;
00709                 if (beg2 == end2) fend2 = true;
00710             } else if (col1 < col2) {
00711                 //sum += (val1 * val2);
00712                 ++beg1;
00713                 if (beg1 == end1) fend1 = true;
00714             } else if (col1 > col2) {
00715                 //sum += (val1 * val2);
00716                 ++beg2;
00717                 if (beg2 == end2) fend2 = true;
00718             }
00719         }
00720         return sum;
00721     }
00722 
00723     static real euclidianDistance(map<int, real>& map1, map<int, real>& map2) 
00724     {
00725         if (map1.size() == 0 || map2.size() == 0)
00726             return 0;
00727         map<int, real>::iterator beg1 = map1.begin();
00728         map<int, real>::iterator beg2 = map2.begin();
00729         map<int, real>::iterator end1 = map1.end();
00730         map<int, real>::iterator end2 = map2.end();
00731         int col1, col2;
00732         real val1, val2, diff, sum = 0;
00733         bool fend1 = (beg1 == end1), fend2 = (beg2 == end2);
00734         int OUT = getMaxColumnIndex(map1, map2) + 1;
00735         
00736         while (!fend1 || !fend2) 
00737         {
00738             if (!fend1)
00739                 col1 = beg1->first;
00740             else
00741                 col1 = OUT;
00742             if (!fend2)
00743                 col2 = beg2->first;
00744             else
00745                 col2 = OUT;
00746             val1 = beg1->second;
00747             val2 = beg2->second;
00748             if (col1 == col2) 
00749             {
00750                 diff = val1 - val2;
00751                 sum += (diff * diff);
00752                 beg1++;
00753                 if (beg1 == end1) fend1 = true;
00754                 beg2++;
00755                 if (beg2 == end2) fend2 = true;
00756             } else if (col1 < col2) 
00757             {
00758                 diff = val1;
00759                 sum += (diff * diff);
00760                 beg1++;
00761                 if (beg1 == end1) fend1 = true;
00762             } else if (col1 > col2) 
00763             {
00764                 diff = val2;
00765                 sum += (diff * diff);
00766                 beg2++;
00767                 if (beg2 == end2) fend2 = true;
00768             }
00769         }
00770         //return sqrt(sum);
00772         return sum;
00773     }
00774 
00776     static void substractVecs(map<int, real>& map1, map<int, real>& map2, Vec& dest) 
00777     {
00778         map<int, real>::iterator beg1 = map1.begin();
00779         map<int, real>::iterator beg2 = map2.begin();
00780         map<int, real>::iterator end1 = map1.end();
00781         map<int, real>::iterator end2 = map2.end();
00782         int col1, col2;
00783         real val1, val2;
00784         bool fend1 = (beg1 == end1), fend2 = (beg2 == end2);
00785         int OUT = getMaxColumnIndex(map1, map2) + 1;
00786       
00787         while (!fend1 || !fend2) 
00788         {
00789             if (!fend1)
00790                 col1 = beg1->first;
00791             else
00792                 col1 = OUT;
00793             if (!fend2)
00794                 col2 = beg2->first;
00795             else
00796                 col2 = OUT;
00797             val1 = beg1->second;
00798             val2 = beg2->second;
00799             if (col1 == col2) 
00800             {
00801                 dest[col1] = val1 - val2;
00802                 beg1++;
00803                 if (beg1 == end1) fend1 = true;
00804                 beg2++;
00805                 if (beg2 == end2) fend2 = true;
00806             } else if (col1 < col2) 
00807             {
00808                 dest[col1] = val1;
00809                 beg1++;
00810                 if (beg1 == end1) fend1 = true;
00811             } else if (col1 > col2) 
00812             {
00813                 dest[col2] = 0 - val2;
00814                 beg2++;
00815                 if (beg2 == end2) fend2 = true;
00816             }
00817         }
00818     } 
00819     
00822     static int getMaxColumnIndex(map<int, T>& map1, map<int, T>& map2) 
00823     {
00824         map<int, real>::iterator end1 = map1.end();
00825         map<int, real>::iterator end2 = map2.end();
00826         --end1;
00827         --end2;
00828         return MAX(end1->first, end2->first);
00829     }
00830 
00835     void exportToMatlabReadableFormat(string filename)
00836     {
00837         ofstream out(filename.c_str());
00838 #ifdef USEDOUBLE
00839         out.precision(20);
00840 #else
00841         out.precision(6);
00842 #endif
00843         for (unsigned int i = 0; i < rows.size(); i++)
00844         {
00845             typename map<int, T>::iterator beg = rows[i].begin();
00846             typename map<int, T>::iterator end = rows[i].end();
00847             while (beg != end)
00848             {
00849                 out << i + 1 << " " << beg->first + 1 << " " << beg->second << endl;
00850                 beg++;
00851             }
00852         }
00855         int l = length();
00856         int w = width();
00857         l--; w--;
00858         map<int, T>& row_l=rows[l];
00859         if (row_l.find(w) == row_l.end())
00860             out << l + 1 << " " << w + 1 << " " << 0 << endl;
00861     }
00862 
00865     real density()
00866     {
00867         return size() / real(length() * width());
00868     }
00869 
00870 };
00871 
00872 template <class T>
00873 void product(RowMapSparseMatrix<T>& M, const Vec& x, Vec& y) { M.product(x,y); }
00874 
00875 // result[j] = sum_i mat[i,j]
00876 template <class T>
00877 void columnSum(RowMapSparseMatrix<T> mat, TVec<T>& result)
00878 {
00879     int n = mat.length();
00880     result.resize(n);
00881     result.clear();
00882     for (int i=0;i<n;i++)
00883     {
00884         map<int,T>& row_i = mat.rows[i];
00885         typename map<int,T>::iterator Rit = row_i.begin();
00886         typename map<int,T>::const_iterator Rend = row_i.end();
00887         for (;Rit!=Rend;++Rit)
00888             result[Rit->first] += Rit->second;
00889     }
00890 }
00891 
00892 // res[i,j] = scale*(mat[i,j] - avg[i] - avg[j] + mean(avg))
00893 template <class T>
00894 void doubleCentering(RowMapSparseMatrix<T>& mat, TVec<T>& avg, RowMapSparseMatrix<T>& res, T scale=1)
00895 {
00896     T moy = mean(avg);
00897     int n=avg.length();
00898     T* a = avg.data();
00899     if (scale==T(1))
00900     {
00901         if (&mat != &res)
00902             for (int i=0;i<n;i++)
00903             {
00904                 map<int,T>& Mi = mat.rows[i];
00905                 map<int,T>& Ri = res.rows[i];
00906                 T term = moy-a[i];
00907                 typename map<int,T>::iterator Mit = Mi.begin();
00908                 typename map<int,T>::const_iterator Mend = Mi.end();
00909                 typename map<int,T>::iterator Rit = Ri.begin();
00910                 typename map<int,T>::const_iterator Rend = Ri.end();
00911                 for (;Mit!=Mend && Rit!=Rend;)
00912                 {
00913                     if (Mit->first==Rit->first)
00914                         Rit->second = Mit->second - a[Rit->first] + term;
00915                     else if (Mit->first<Rit->first)
00916                         ++Mit;
00917                     else
00918                         ++Rit;
00919                 }
00920             }
00921         else
00922             for (int i=0;i<n;i++)
00923             {
00924                 map<int,T>& Ri = res.rows[i];
00925                 typename map<int,T>::iterator Rit = Ri.begin();
00926                 typename map<int,T>::const_iterator Rend = Ri.end();
00927                 T term = moy-a[i];
00928                 for (;Rit!=Rend;++Rit)
00929                     Rit->second += term - a[Rit->first];
00930 
00931             }
00932     }
00933     else
00934     {
00935         if (&mat != &res)
00936             for (int i=0;i<n;i++)
00937             {
00938                 map<int,T>& Mi = mat.rows[i];
00939                 map<int,T>& Ri = res.rows[i];
00940                 T term = moy-a[i];
00941                 typename map<int,T>::iterator Mit = Mi.begin();
00942                 typename map<int,T>::const_iterator Mend = Mi.end();
00943                 typename map<int,T>::iterator Rit = Ri.begin();
00944                 typename map<int,T>::const_iterator Rend = Ri.end();
00945                 for (;Mit!=Mend && Rit!=Rend;)
00946                 {
00947                     if (Mit->first==Rit->first)
00948                         Rit->second = scale*(Mit->second - a[Rit->first] + term);
00949                     else if (Mit->first<Rit->first)
00950                         ++Mit;
00951                     else
00952                         ++Rit;
00953                 }
00954             }
00955         else
00956             for (int i=0;i<n;i++)
00957             {
00958                 map<int,T>& Ri = res.rows[i];
00959                 typename map<int,T>::iterator Rit = Ri.begin();
00960                 typename map<int,T>::const_iterator Rend = Ri.end();
00961                 T term = moy-a[i];
00962                 for (;Rit!=Rend;++Rit)
00963                     Rit->second = scale*(Rit->second + term - a[Rit->first]);
00964 
00965             }
00966     }
00967 }
00968 
00969 
00970 
00971 template <class T>
00972 void averageAcrossRowsAndColumns(RowMapSparseMatrix<T> mat, 
00973                                  Vec avg_across_rows, Vec avg_across_columns, 
00974                                  bool only_on_non_zeros=true)
00975 {
00976     mat.averageAcrossRowsAndColumns(avg_across_rows, avg_across_columns, only_on_non_zeros);
00977 }
00978 
00979 template <class T>
00980 void addToRows(RowMapSparseMatrix<T> mat, Vec row, bool only_on_non_zeros=true)
00981 {
00982     mat.addToRows(row,only_on_non_zeros);
00983 }
00984 
00985 template <class T>
00986 void addToColumns(RowMapSparseMatrix<T> mat, Vec row, bool only_on_non_zeros=true)
00987 {
00988     mat.addToColumns(row,only_on_non_zeros);
00989 }
00990 
00991 } // end of namespace PLearn
00992 
00993 #endif
00994 
00995 
00996 /*
00997   Local Variables:
00998   mode:c++
00999   c-basic-offset:4
01000   c-file-style:"stroustrup"
01001   c-file-offsets:((innamespace . 0)(inline-open . 0))
01002   indent-tabs-mode:nil
01003   fill-column:79
01004   End:
01005 */
01006 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines