PLearn 0.1
|
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 :