PLearn 0.1
TMat_impl.h
Go to the documentation of this file.
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines