PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // PLearn (A C++ Machine Learning Library) 00004 // Copyright (C) 1998 Pascal Vincent 00005 // Copyright (C) 1999,2000 Pascal Vincent, Yoshua Bengio and University of Montreal 00006 // 00007 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions are met: 00010 // 00011 // 1. Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // 00014 // 2. Redistributions in binary form must reproduce the above copyright 00015 // notice, this list of conditions and the following disclaimer in the 00016 // documentation and/or other materials provided with the distribution. 00017 // 00018 // 3. The name of the authors may not be used to endorse or promote 00019 // products derived from this software without specific prior written 00020 // permission. 00021 // 00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR 00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN 00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00032 // 00033 // This file is part of the PLearn library. For more information on the PLearn 00034 // library, go to the PLearn Web site at www.plearn.org 00035 00036 00037 00038 00039 /* ******************************************************* 00040 * $Id: TMat_impl.h 10000 2009-03-10 14:33:28Z tihocan $ 00041 * AUTHORS: Pascal Vincent & Yoshua Bengio & Rejean Ducharme 00042 * This file is part of the PLearn library. 00043 ******************************************************* */ 00044 00047 #ifndef TMAT_IMPL_H 00048 #define TMAT_IMPL_H 00049 00050 #include "TMat_decl.h" 00051 #include "TMatElementIterator_impl.h" 00052 #include "TMatRowsIterator_impl.h" 00053 #include "TMatRowsAsArraysIterator_impl.h" 00054 #include "TMatColRowsIterator_impl.h" 00055 00056 //#include "algo.h" 00057 00058 namespace PLearn { 00059 using namespace std; 00060 00061 00062 // ************** 00063 // **** TVec **** 00064 // ************** 00065 00066 00070 template <class T> 00071 TVec<T>::TVec(const T& start, const T& stop, const T& step) 00072 :length_(0), offset_(0) 00073 { 00074 // first count the size n 00075 T val; 00076 int n=0; 00077 for(val=start; val<=stop; val+=step) 00078 ++n; 00079 00080 if(n) 00081 { 00082 resize(n); 00083 iterator it = begin(); 00084 iterator itend = end(); 00085 for(val=start; it!=itend; ++it, val+=step) 00086 *it = val; 00087 } 00088 } 00089 00090 00092 template <class T> 00093 TMat<T> TVec<T>::toMat(int newlength, int newwidth) const 00094 { 00095 TMat<T> tm; 00096 tm.offset_ = offset_; 00097 tm.mod_ = newwidth; 00098 tm.width_ = newwidth; 00099 tm.length_ = newlength; 00100 tm.storage = storage; 00101 return tm; 00102 } 00103 00104 00105 template <class T> 00106 void TVec<T>::input(istream& in) const 00107 { 00108 T* v = data(); 00109 for(int i=0; i<length(); i++) 00110 { 00111 if(!(in>>v[i])) 00112 PLERROR("In TVec::input error encountered while reading vector"); 00113 00114 } 00115 } 00116 00117 template <class T> 00118 void TVec<T>::input(PStream& in) const 00119 { 00120 T* v = data(); 00121 int l = length(); 00122 for(int i=0; i<l; i++) 00123 { 00124 in.skipBlanksAndCommentsAndSeparators(); 00125 if(in.peek()==EOF || in.eof()) 00126 PLERROR("In TVec::input encountered EOF before reading the last element of %d",l); 00127 in>>v[i]; 00128 } 00129 } 00130 00131 template <class T> 00132 void TVec<T>::print(ostream& out) const 00133 { 00134 if(storage && 0 < length()) 00135 { 00136 out.setf(ios::fmtflags(0),ios::floatfield); 00137 T* v = data(); 00138 for(int i=0; i<length(); i++) 00139 out << setiosflags(ios::left) << setprecision(7) << setw(11) << v[i] << ' '; 00140 out.flush(); 00141 } 00142 } 00143 00144 template <class T> 00145 void TVec<T>::print(ostream& out, const string& separator) const 00146 { 00147 out.setf(ios::fmtflags(0),ios::floatfield); 00148 T* v = data(); 00149 for(int i=0; i<length()-1; i++) 00150 out << v[i] << separator; 00151 out << v[length()-1]; 00152 out.flush(); 00153 } 00154 00155 template <class T> 00156 void TVec<T>::printcol(ostream& out) const 00157 { 00158 T* v = data(); 00159 for(int i=0; i<length(); i++) 00160 out << v[i] << "\n"; 00161 out.flush(); 00162 } 00163 00164 00165 00166 00167 // *********************** 00168 // * Fonctions pout TVec * 00169 // *********************** 00170 00176 template<class T, class I> 00177 void selectElements(const TVec<T>& source, const TVec<I>& indices, TVec<T>& destination) 00178 { 00179 int ni = indices.length(); 00180 if (ni!=destination.length()) 00181 PLERROR("select(Vec,Vec,Vec): last 2 arguments have lengths %d != %d", 00182 indices.length(),destination.length()); 00183 I* indx = indices.data(); 00184 T* dest = destination.data(); 00185 T* src = source.data(); 00186 #ifdef BOUNDCHECK 00187 int n=source.length(); 00188 #endif 00189 for (int i=0;i<ni;i++) 00190 { 00191 int pos = int(indx[i]); 00192 #ifdef BOUNDCHECK 00193 if (pos<0 || pos>=n) 00194 PLERROR("select(Vec,Vec,Vec) indices[%d]=%d out of bounds (0,%d)", 00195 i,pos,n-1); 00196 #endif 00197 dest[i] = src[pos]; 00198 } 00199 } 00200 00202 template<class T> 00203 void elementsEqualTo(const TVec<T>& source, const T& value, const TVec<T>& destination) 00204 { 00205 #ifdef BOUNDCHECK 00206 if (source.length()!=destination.length()) 00207 PLERROR("elementsEqualTo(Vec(%d),%f,Vec(%d)): incompatible dimensions", 00208 source.length(),value,destination.length()); 00209 #endif 00210 T* src=source.data(); 00211 T* dst=destination.data(); 00212 for (int i=0;i<destination.length();i++) 00213 if (src[i]==value) dst[i]=1.0; 00214 else dst[i]=0.0; 00215 } 00216 00220 template<class T> 00221 TVec<T> removeElement(const TVec<T>& v, int elemnum) 00222 { 00223 if(elemnum==0) 00224 return v.subVec(1,v.length()-1); 00225 else if(elemnum==v.length()-1) 00226 return v.subVec(0,v.length()-1); 00227 else 00228 return concat(v.subVec(0,elemnum), 00229 v.subVec(elemnum+1,v.length()-(elemnum+1))); 00230 } 00231 00232 // Returns an index vector I so that (*this)(I) returns a sorted version 00233 // of this vec in ascending order. 00234 namespace { 00235 template <class T> 00236 struct index_cmp : public binary_function<int, int, bool> 00237 { 00238 const TVec<T>& m_values; 00239 index_cmp(const Vec& values): m_values(values) { } 00240 bool operator()(int x, int y) { return m_values[x] < m_values[y]; } 00241 }; 00242 template <class T> 00243 struct index_missing_cmp : public binary_function<int, int, bool> 00244 { 00245 const TVec<T>& m_values; 00246 index_missing_cmp(const Vec& values): m_values(values) { } 00247 bool operator()(int x, int y) { 00248 const T& v1 = m_values[x]; 00249 const T& v2 = m_values[y]; 00250 if (is_missing(v1)) 00251 return false; 00252 else if (is_missing(v2)) 00253 return true; 00254 else 00255 return v1 < v2; 00256 } 00257 }; 00258 } 00259 // Actual body of the method 00260 template <class T> 00261 TVec<int> TVec<T>::sortingPermutation(bool stable, bool missing) const 00262 { 00263 TVec<int> indices(length_); 00264 for (int i=0; i < length_; i++) indices[i] = i; 00265 if(stable && ! missing) 00266 stable_sort(indices.begin(), indices.end(), index_cmp<T>(*this)); 00267 else if(! stable && !missing) 00268 sort(indices.begin(), indices.end(), index_cmp<T>(*this)); 00269 else if(stable && missing) 00270 stable_sort(indices.begin(), indices.end(),index_missing_cmp<T>(*this)); 00271 else if(!stable && missing) 00272 sort(indices.begin(), indices.end(), index_missing_cmp<T>(*this)); 00273 return indices; 00274 } 00275 00276 00277 // ************** 00278 // **** TMat **** 00279 // ************** 00280 00281 template <class T> 00282 TMat<T>::TMat(int the_length, int the_width, const TVec<T>& v) 00283 : offset_(v.offset()), mod_(the_width), length_(the_length), width_(the_width), storage(v.storage) 00284 { 00285 if(length()*width()!=v.length()) 00286 PLERROR("In Mat constructor from Vec: length()*width() of matrix must be equal to length() of Vec"); 00287 } 00288 00289 00290 template <class T> 00291 TVec<T> TMat<T>::toVecCopy() const 00292 { 00293 TVec<T> v(length()*width()); 00294 v << *this; 00295 return v; 00296 } 00297 00300 00301 template <class T> 00302 TVec<T> TMat<T>::toVec() const 00303 { 00304 if(length()>1 && width()<mod()) 00305 PLERROR("In Mat::toVec internal structure of this Mat makes it impossible to build a Vec that would view exactly the same data. Consider using toVecCopy() instead!"); 00306 00307 TVec<T> v; 00308 v.offset_ = offset_; 00309 v.length_ = length()*width(); 00310 v.storage = storage; 00311 return v; 00312 } 00313 00314 template <class T> 00315 int TMat<T>::findRow(const TVec<T>& row) const 00316 { 00317 for(int i=0; i<length(); i++) 00318 if( (*this)(i)==row ) 00319 return i; 00320 return -1; 00321 } 00322 00323 template <class T> 00324 void TMat<T>::appendRow(const TVec<T>& newrow) 00325 { 00326 #ifdef BOUNDCHECK 00327 if(newrow.length()!=width() && width() > 0) 00328 PLERROR("In TMat::appendRow newrow vector should be as long as the matrix is wide (%d != %d)", newrow.length(), width()); 00329 #endif 00330 if (storage) { 00331 resize(length()+1, newrow.length(), storage->length()); 00332 } else { 00333 // This Mat is empty: it has no storage, so using storage would crash. 00334 resize(length()+1, newrow.length()); 00335 } 00336 (*this)(length()-1) << newrow; 00337 } 00338 00339 00340 // C++ stream output 00341 template <class T> 00342 void TMat<T>::print(ostream& out) const 00343 { 00344 out.flags(ios::left); 00345 for(int i=0; i<length(); i++) 00346 { 00347 const T* m_i = rowdata(i); 00348 for(int j=0; j<width(); j++) 00349 out << setw(11) << m_i[j] << ' '; 00350 out << "\n"; 00351 } 00352 out.flush(); 00353 } 00354 00355 template <class T> 00356 void TMat<T>::input(istream& in) const 00357 { 00358 for(int i=0; i<length(); i++) 00359 { 00360 T* v = rowdata(i); 00361 for (int j=0;j<width();j++) 00362 { 00363 if(!(in>>v[j])) 00364 PLERROR("In TMat<T>::input error encountered while reading matrix"); 00365 } 00366 } 00367 } 00368 00369 template <class T> 00370 void TMat<T>::input(PStream& in) const 00371 { 00372 for(int i=0; i<length(); i++) 00373 { 00374 T* v = rowdata(i); 00375 for (int j=0;j<width();j++) 00376 { 00377 if (!in) 00378 PLERROR("In TMat<T>::input error encountered while reading matrix"); 00379 else 00380 in>>v[j]; 00381 } 00382 } 00383 } 00384 00385 template <class T> 00386 void TMat<T>::resizePreserve(int new_length, int new_width, int extra) 00387 { 00388 int usage = storage->usage(); 00389 int new_size = new_length*MAX(mod(),new_width); 00390 int new_offset = usage>1?offset_:0; 00391 if (new_size>storage->length() || new_width>mod()) 00392 { 00393 int extracols=0, extrarows=0; 00394 if (extra>min(new_width,new_length)) 00395 { 00396 // if width has increased, bet that it will increase again in the future, 00397 // similarly for length, so allocate the extra as extra mod 00398 float l=float(length_), l1=float(new_length), 00399 w=float(width_), w1=float(new_width), 00400 x=float(extra); 00401 // Solve the following equations to apportion the extra 00402 // while keeping the same percentage increase in width and length: 00403 // Solve[{x+w1*l1==w2*l2,(w2/w1 - 1)/(l2/l1 - 1) == (w1/w - 1)/(l1/l - 1)},{w2,l2}] 00404 // This is a quadratic system which has two solutions: {w2a,l2a} and {w2b,l2b}: 00405 float w2a = 00406 w1*(-1 - l1/(l - l1) + w1/w + (l1*w1)/(l*w - l1*w) + 00407 (2*l*(-w + w1)*x)/ 00408 (2*l*l1*w*w1 - l1*l1*w*w1 - l*l1*w1*w1 + 00409 sqrt(square(l1*l1*w*w1 - l*l1*w1*w1) + 00410 4*l*(l - l1)*l1*w*(w - w1)*w1*(l1*w1 + x)))); 00411 float l2a = -(-l1*l1*w*w1 + l*l1*w1*w1 + 00412 sqrt(square(l1*l1*w*w1 - l*l1*w1*w1) + 00413 4*l*(l - l1)*l1*w*(w - w1)*w1*(l1*w1 + x)))/(2*l*(w - w1)*w1); 00414 float w2b =w1*(-1 - l1/(l - l1) + w1/w + (l1*w1)/(l*w - l1*w) - 00415 (2*l*(-w + w1)*x)/ 00416 (-2*l*l1*w*w1 + l1*l1*w*w1 + l*l1*w1*w1 + 00417 sqrt(square(l1*l1*w*w1 - l*l1*w1*w1) + 00418 4*l*(l - l1)*l1*w*(w - w1)*w1*(l1*w1 + x)))); 00419 float l2b = (l1*l1*w*w1 - l*l1*w1*w1 + 00420 sqrt(square(l1*l1*w*w1 - l*l1*w1*w1) + 00421 4*l*(l - l1)*l1*w*(w - w1)*w1*(l1*w1 + x)))/(2*l*(w - w1)*w1); 00422 00423 // pick one that is feasible and maximizes the mod 00424 if (w2b>w2a && w2b>w1 && l2b>l1) { 00425 extracols=int(ceil(w2b-w1)); 00426 extrarows=int(ceil(l2b-l1)); 00427 } 00428 else if (w2a>w1 && l2a>l1) { 00429 extrarows=int(ceil(l2a-l1)); 00430 extracols=int(ceil(w2a-w1)); 00431 } 00432 else { // no valid solution to the system of equation, use a heuristic 00433 extracols = max(0,int(ceil(sqrt(real(extra))/new_length))); 00434 extrarows = max(0,int((extra+l1*w1)/(w1+extracols) - l1)); 00435 } 00436 00437 } 00438 storage->resizeMat(new_length,new_width,extrarows,extracols, 00439 new_offset,mod_,length_,width_,offset_); 00440 mod_ = new_width + extracols; 00441 } 00442 offset_ = new_offset; 00443 } 00444 00445 template <class T> 00446 inline void TMat<T>::resizeBoundCheck(int new_length, int new_width) 00447 { 00448 if(new_length<0 || new_width<0) 00449 PLERROR("IN TMat::resize(int new_length, int new_width)\nInvalid arguments (%d, %d)", new_length, new_width); 00450 } 00451 00452 template <class T> 00453 void TMat<T>::resizeModError() 00454 { 00455 PLERROR("IN TMat::resize(int new_length, int new_width) - For safety " 00456 "reasons, increasing the width() beyond mod()-offset_ modulo " 00457 "mod() is not allowed when the storage is shared with others"); 00458 } 00459 00460 00461 // Deep copying 00462 00463 template<class T> 00464 void TMat<T>::makeDeepCopyFromShallowCopy(CopiesMap& copies) 00465 { 00466 deepCopyField(storage, copies); 00467 } 00468 00469 template<class T> 00470 TMat<T> TMat<T>::deepCopy(CopiesMap& copies) const 00471 { 00472 // First do a shallow copy 00473 TMat<T> deep_copy = *this; 00474 // Transform the shallow copy into a deep copy 00475 deep_copy.makeDeepCopyFromShallowCopy(copies); 00476 // return the completed deep_copy 00477 return deep_copy; 00478 } 00479 00480 // Iterateurs 00481 00482 template<class T> 00483 TMatElementIterator<T> TMat<T>::begin() const 00484 { return TMatElementIterator<T>(data(), width_, mod_); } 00485 00486 template<class T> 00487 TMatElementIterator<T> TMat<T>::end() const 00488 { return TMatElementIterator<T>(data()+length_*mod_, width_, mod_); } 00489 00490 00491 template<class T> 00492 TMatRowsIterator<T> TMat<T>::rows_begin() { 00493 return TMatRowsIterator<T>(data(), width_, mod_); 00494 } 00495 00496 template<class T> 00497 TMatRowsIterator<T> TMat<T>::rows_end() { 00498 return TMatRowsIterator<T>(data()+length_*mod_, width_, mod_); 00499 } 00500 00501 00502 template<class T> 00503 TMatRowsAsArraysIterator<T> TMat<T>::rows_as_arrays_begin() { 00504 return TMatRowsAsArraysIterator<T>(data(), width_, mod_); 00505 } 00506 00507 template<class T> 00508 TMatRowsAsArraysIterator<T> TMat<T>::rows_as_arrays_end() { 00509 return TMatRowsAsArraysIterator<T>(data()+length_*mod_, width_, mod_); 00510 } 00511 00512 template<class T> 00513 TMatColRowsIterator<T> TMat<T>::col_begin(int column) { 00514 return TMatColRowsIterator<T>(data() + column, mod_); 00515 } 00516 00517 template<class T> 00518 TMatColRowsIterator<T> TMat<T>::col_end(int column) { 00519 return TMatColRowsIterator<T>(data()+length_*mod_+column, mod_); 00520 } 00521 00522 template<class T> 00523 bool TMat<T>::operator==(const TMat<T>& other) const 00524 { 00525 if ( length() != other.length() || width() != other.width() ) 00526 return false; 00527 00528 iterator it = begin(); 00529 iterator end_ = end(); 00530 iterator other_it = other.begin(); 00531 00532 for(; it != end_; ++it, ++other_it) 00533 if(*it != *other_it) 00534 return false; 00535 00536 return true; 00537 } 00538 00539 template<class T> 00540 bool TMat<T>::isEqual(const TMat<T>& other, real precision) const 00541 { 00542 if ( length() != other.length() || width() != other.width() ) 00543 return false; 00544 00545 iterator it = begin(); 00546 iterator end_ = end(); 00547 iterator other_it = other.begin(); 00548 00549 for(; it != end_; ++it, ++other_it) 00550 if( !is_equal(*it,*other_it, 10.0, precision, precision) ) 00551 return false; 00552 00553 return true; 00554 } 00555 00556 00557 00558 // ***************************** 00559 // **** Fonctions pour TMat **** 00560 // ***************************** 00561 00562 00563 // select the rows of the source as specified by the 00564 // vector of indices (between 0 and source.length()-1), copied into 00565 // the destination matrix (which must have the same length() 00566 // as the indices vector). 00567 template <class T, class I> 00568 void selectRows(const TMat<T>& source, const TVec<I>& row_indices, TMat<T>& destination) 00569 { 00570 int ni = row_indices.length(); 00571 if (ni!=destination.length()) 00572 PLERROR("selectRows(Mat,Vec,Mat): last 2 arguments have lengths %d != %d", 00573 ni,destination.length()); 00574 00575 if (row_indices.isEmpty()) 00576 // Nothing to select. In addition, 'destination' is empty too since it 00577 // has zero length, according to the test above. Thus there is nothing 00578 // to do. 00579 return; 00580 00581 I* indx = row_indices.data(); 00582 #ifdef BOUNDCHECK 00583 int n=source.length(); 00584 #endif 00585 for (int i=0;i<ni;i++) 00586 { 00587 int pos = int(indx[i]); 00588 #ifdef BOUNDCHECK 00589 if (pos<0 || pos>=n) 00590 PLERROR("selectRows(Mat,Vec,Mat) indices[%d]=%d out of bounds (0,%d)", 00591 i,pos,n-1); 00592 #endif 00593 destination(i) << source(pos); 00594 } 00595 } 00596 00597 // select the colums of the source as specified by the 00598 // vector of indices (between 0 and source.width()-1), copied into 00599 // the destination matrix (which must have the same width() 00600 // as the indices vector). 00601 template <class T, class I> 00602 void selectColumns(const TMat<T>& source, const TVec<I>& column_indices, TMat<T>& destination) 00603 { 00604 int ni = column_indices.length(); 00605 if (ni!=destination.width()) 00606 PLERROR("selectColums(Mat,Vec,Mat): last 2 arguments have dimensions %d != %d", 00607 ni,destination.width()); 00608 00609 if (column_indices.isEmpty()) 00610 // Nothing to select. In addition, 'destination' is empty too since it 00611 // has zero width, according to the test above. Thus there is nothing 00612 // to do. 00613 return; 00614 00615 I* indx = column_indices.data(); 00616 #ifdef BOUNDCHECK 00617 int n=source.width(); 00618 #endif 00619 for (int i=0;i<ni;i++) 00620 { 00621 int pos = int(indx[i]); 00622 #ifdef BOUNDCHECK 00623 if (pos<0 || pos>=n) 00624 PLERROR("selectColumns(Mat,Vec,Mat) indices[%d]=%d out of bounds (0,%d)", 00625 i,pos,n-1); 00626 #endif 00627 destination.column(i) << source.column(pos); 00628 } 00629 } 00630 00631 // select a submatrix of specified rows and colums of the source with 00632 // two vectors of indices. The elements that are both in the specified rows 00633 // and columns are copied into the destination matrix (which must have the 00634 // same length() as the row_indices vector, and the same width() as the length() 00635 // of the col_indices vector). 00636 template <class T, class I> 00637 void select(const TMat<T>& source, const TVec<I>& row_indices, const TVec<I>& column_indices, TMat<T>& destination) 00638 { 00639 int rni = row_indices.length(); 00640 int cni = column_indices.length(); 00641 if (rni!=destination.length() || cni!=destination.width()) 00642 PLERROR("select(Mat(%d,%d),Vec(%d),Vec(%d),Mat(%d,%d)): arguments have incompatible dimensions", 00643 source.length(),source.width(),rni,cni,destination.length(),destination.width()); 00644 I* rindx = row_indices.data(); 00645 I* cindx = column_indices.data(); 00646 #ifdef BOUNDCHECK 00647 int nr=source.length(); 00648 int nc=source.width(); 00649 #endif 00650 for (int i=0;i<rni;i++) 00651 { 00652 int ri=(int)rindx[i]; 00653 #ifdef BOUNDCHECK 00654 if (ri<0 || ri>=nr) 00655 PLERROR("select(Mat,Vec,Vec,Mat) row_indices[%d]=%d out of bounds (0,%d)", 00656 i,ri,nr-1); 00657 #endif 00658 T* dest_row = destination[i]; 00659 T* src_row = source[ri]; 00660 for (int j=0;j<cni;j++) 00661 { 00662 int cj = (int)cindx[j]; 00663 #ifdef BOUNDCHECK 00664 if (cj<0 || cj>=nc) 00665 PLERROR("select(Mat,Vec,Vec,Mat) col_indices[%d]=%d out of bounds (0,%d)", 00666 i,cj,nc-1); 00667 #endif 00668 dest_row[j] = src_row[cj]; 00669 } 00670 } 00671 } 00672 00673 template<class T> 00674 TMat<T> removeRow(const TMat<T>& m, int rownum) 00675 { 00676 if(rownum==0) 00677 return m.subMatRows(1,m.length()-1); 00678 else if(rownum==m.length()-1) 00679 return m.subMatRows(0,m.length()-1); 00680 else 00681 return vconcat(m.subMatRows(0,rownum), 00682 m.subMatRows(rownum+1,m.length()-(rownum+1))); 00683 } 00684 00685 template<class T> 00686 TMat<T> removeColumn(const TMat<T>& m, int colnum) 00687 { 00688 if(colnum==0) 00689 return m.subMatColumns(1,m.width()-1); 00690 else if(colnum==m.width()-1) 00691 return m.subMatColumns(0,m.width()-1); 00692 else 00693 return hconcat(m.subMatColumns(0,colnum), 00694 m.subMatColumns(colnum+1,m.width()-(colnum+1))); 00695 } 00696 00697 template<class T> 00698 TMat<T> diagonalmatrix(const TVec<T>& v) 00699 { 00700 TMat<T> m(v.length(), v.length()); 00701 for(int i=0; i<v.length(); i++) 00702 m(i,i) = v[i]; 00703 return m; 00704 } 00705 00706 00707 00708 // ***************************** 00709 // **** Fonctions pour TMat **** 00710 // ***************************** 00711 00712 template <class T> inline TMat<T> deepCopy(const TMat<T> source) 00713 { 00714 CopiesMap copies; 00715 return deepCopy(source, copies); 00716 } 00717 00718 template <class T> inline TMat<T> 00719 deepCopy(const TMat<T> source, CopiesMap copies) 00720 { return source.deepCopy(copies); } 00721 00722 template <class T> 00723 inline void deepCopyField(TMat<T>& field, CopiesMap& copies) 00724 { 00725 field.makeDeepCopyFromShallowCopy(copies); 00726 } 00727 00728 template<class T> 00729 void clear(const TMat<T>& x) 00730 { 00731 if(x.isCompact()) 00732 { 00733 typename TMat<T>::compact_iterator it = x.compact_begin(); 00734 typename TMat<T>::compact_iterator itend = x.compact_end(); 00735 for(; it!=itend; ++it) 00736 clear(*it); 00737 } 00738 else 00739 { 00740 typename TMat<T>::iterator it = x.begin(); 00741 typename TMat<T>::iterator itend = x.end(); 00742 for(; it!=itend; ++it) 00743 clear(*it); 00744 } 00745 } 00746 00747 template<class T> 00748 void swap( TMat<T>& a, TMat<T>& b) 00749 { swap_ranges(a.begin(), a.end(), b.begin()); } 00750 00752 template<class T> 00753 inline void operator<<(const TMat<T>& m1, const TMat<T>& m2) 00754 { 00755 #ifdef BOUNDCHECK 00756 if(m1.size()!=m2.size()) 00757 PLERROR("In operator<<(m1,m2) the 2 matrices must have the same number of elements\n" 00758 "m1: (%d, %d) && m2: (%d, %d)", m1.length(), m1.width(), m2.length(), m2.width()); 00759 #endif 00760 if (m1.isNotEmpty()) 00761 copy(m2.begin(), m2.end(), m1.begin()); 00762 } 00763 00765 template<class T, class U> 00766 void operator<<(const TMat<T>& m1, const TMat<U>& m2) 00767 { 00768 #ifdef BOUNDCHECK 00769 if(m1.size()!=m2.size()) 00770 PLERROR("In operator<<(m1,m2) the 2 matrices must have the same number of elements"); 00771 #endif 00772 if (m1.isNotEmpty()) 00773 copy_cast(m2.begin(), m2.end(), m1.begin()); 00774 } 00775 00777 template<class T> 00778 inline void operator<<(const TMat<T>& m1, const TVec<T>& m2) 00779 { 00780 #ifdef BOUNDCHECK 00781 if(m1.size()!=m2.size()) 00782 PLERROR("In operator<<(m1,m2) the 2 matrices must have the same number of elements;\t m1.size()= %d;\t m2.size= %d", m1.size(), m2.size()); 00783 #endif 00784 if (m1.isNotEmpty()) 00785 copy(m2.begin(), m2.end(), m1.begin()); 00786 } 00787 00789 template<class T, class U> 00790 inline void operator<<(const TMat<T>& m1, const TVec<U>& m2) 00791 { 00792 #ifdef BOUNDCHECK 00793 if(m1.size()!=m2.size()) 00794 PLERROR("In operator<<(m1,m2) the 2 matrices must have the same number of elements"); 00795 #endif 00796 if (m1.isNotEmpty()) 00797 copy_cast(m2.begin(), m2.end(), m1.begin()); 00798 } 00799 00801 template<class T> 00802 inline void operator<<(const TVec<T>& m1, const TMat<T>& m2) 00803 { 00804 #ifdef BOUNDCHECK 00805 if(m1.size()!=m2.size()) 00806 PLERROR("In operator<<(m1,m2) the 2 matrices must have the same number of elements"); 00807 #endif 00808 if (m1.isNotEmpty()) 00809 copy(m2.begin(), m2.end(), m1.begin()); 00810 } 00811 00813 template<class T, class U> 00814 inline void operator<<(const TVec<T>& m1, const TMat<U>& m2) 00815 { 00816 #ifdef BOUNDCHECK 00817 if(m1.size()!=m2.size()) 00818 PLERROR("In operator<<(m1,m2) the 2 matrices must have the same number of elements"); 00819 #endif 00820 if (m1.isNotEmpty()) 00821 copy_cast(m2.begin(), m2.end(), m1.begin()); 00822 } 00823 00825 template<class T, class U> 00826 inline void operator>>(const TMat<T>& m1, const TMat<U>& m2) 00827 { m2 << m1; } 00828 00830 template<class T, class U> 00831 inline void operator>>(const TVec<T>& m1, const TMat<U>& m2) 00832 { m2 << m1; } 00833 00835 template<class T, class U> 00836 inline void operator>>(const TMat<T>& m1, const TVec<U>& m2) 00837 { m2 << m1; } 00838 00839 00840 00842 template <class T> 00843 inline ostream& operator<<(ostream& out, const TMat<T>& m) 00844 { 00845 m.print(out); 00846 return out; 00847 } 00848 00850 00851 template <class T> 00852 inline istream& operator>>(istream& in, const TMat<T>& m) 00853 { 00854 m.input(in); 00855 return in; 00856 } 00857 00859 template <class T> 00860 inline TMat<T> rowmatrix(const TVec<T>& v) 00861 { return v.toMat(1,v.length()); } 00862 00864 template <class T> 00865 inline TMat<T> columnmatrix(const TVec<T>& v) 00866 { return v.toMat(v.length(),1); } 00867 00868 // select the rows of the source as specified by the 00869 // vector of indices (between 0 and source.length()-1), copied into 00870 // the destination matrix (which must have the same length() 00871 // as the indices vector). 00872 template <class T, class I> 00873 void selectRows(const TMat<T>& source, const TVec<I>& row_indices, TMat<T>& destination); 00874 00875 // select the colums of the source as specified by the 00876 // vector of indices (between 0 and source.length()-1), copied into 00877 // the destination matrix (which must have the same width() 00878 // as the indices vector). 00879 template <class T, class I> 00880 void selectColumns(const TMat<T>& source, const TVec<I>& column_indices, TMat<T>& destination); 00881 00882 // select a submatrix of specified rows and colums of the source with 00883 // two vectors of indices. The elements that are both in the specified rows 00884 // and columns are copied into the destination matrix (which must have the 00885 // same length() as the row_indices vector, and the same width() as the length() 00886 // of the col_indices vector). 00887 template <class T> 00888 void select(const TMat<T>& source, const TVec<T>& row_indices, const TVec<T>& column_indices, TMat<T>& destination); 00889 00894 template<class T> 00895 TMat<T> removeRow(const TMat<T>& m, int rownum); 00896 00901 template<class T> 00902 TMat<T> removeColumn(const TMat<T>& m, int colnum); 00903 00904 00905 template<class T> 00906 TMat<T> diagonalmatrix(const TVec<T>& v); 00907 00908 // old .pmat format 00909 template<class T> 00910 void savePMat(const string& filename, const TMat<T>& mat) 00911 { PLERROR("savePMat only implemented for float and double"); } 00912 00913 template<class T> 00914 void loadPMat(const string& filename, TMat<float>& mat) 00915 { PLERROR("loadPMat only implemented for float and double"); } 00916 00917 inline void deepCopyField(Mat*& field, CopiesMap& copies) 00918 { 00919 if (field) 00920 { 00921 CopiesMap::iterator it = copies.find(field); 00922 if (it != copies.end()) 00923 field = static_cast<Mat*>(it->second); 00924 else 00925 { 00926 // Throw an error. The reason is that: 00927 // - if the Mat* pointer points to a matrix that has already been 00928 // deep-copied, I am unsure whether 'copies' contains the correct 00929 // pointer, thus we may end up here even though we should reuse the 00930 // previous deep copy, 00931 // - if the Mat* pointer points to a matrix that has not been deep 00932 // copied yet, then when that matrix is deep copied it will not 00933 // actually be the same as the matrix we may create here. 00934 PLERROR("In deepCopyField(Mat*& field, CopiesMap& copies) - You " 00935 "cannot deep copy a Mat* directly."); 00936 /* Old code. 00937 Mat* newM = new Mat; 00938 (*newM) = field->deepCopy(copies); 00939 copies[field] = newM; 00940 field = newM; 00941 */ 00942 } 00943 } 00944 } 00945 00946 00949 00953 00954 template <class T> inline PStream & 00955 operator<<(PStream &out, const TMat<T> &m) 00956 { 00957 m.write(out); 00958 return out; 00959 } 00960 00961 template <class T> 00962 PStream & operator>>(PStream &in, TMat<T> &m) 00963 { 00964 m.read(in); 00965 return in; 00966 } 00967 00968 inline string join(const TVec<string>& s, const string& separator) 00969 { 00970 string result; 00971 for(int i=0; i<s.size(); i++) 00972 { 00973 result += s[i]; 00974 if(i<s.size()-1) 00975 result += separator; 00976 } 00977 return result; 00978 } 00979 00980 } // end of namespace PLearn 00981 00982 #endif // TMAT_IMPL_H 00983 00984 00985 /* 00986 Local Variables: 00987 mode:c++ 00988 c-basic-offset:4 00989 c-file-style:"stroustrup" 00990 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00991 indent-tabs-mode:nil 00992 fill-column:79 00993 End: 00994 */ 00995 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :