Marsyas  0.5.0-beta1
/Users/jleben/code/marsyas/src/marsyas/realvec.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2010 George Tzanetakis <gtzan@cs.uvic.ca>
00003 **
00004 ** This program is free software; you can redistribute it and/or modify
00005 ** it under the terms of the GNU General Public License as published by
00006 ** the Free Software Foundation; either version 2 of the License, or
00007 ** (at your option) any later version.
00008 **
00009 ** This program is distributed in the hope that it will be useful,
00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 ** GNU General Public License for more details.
00013 **
00014 ** You should have received a copy of the GNU General Public License
00015 ** along with this program; if not, write to the Free Software
00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 */
00018 
00019 #include <marsyas/realvec.h>
00020 #include <marsyas/common_source.h>
00021 #include <marsyas/NumericLib.h>
00022 #include "Communicator.h"
00023 
00024 #include <algorithm>
00025 #include <limits>
00026 
00027 #ifdef MARSYAS_MATLAB
00028 //#define _MATLAB_REALVEC_
00029 #endif
00030 
00031 using namespace std;
00032 
00033 namespace Marsyas
00034 {
00035 
00036 
00037 
00039 realvec::realvec()
00040   : size_(0),
00041     allocatedSize_(0),
00042     data_(NULL),
00043     rows_(1),
00044     cols_(0)
00045 {
00046   allocateData(size_);
00047 }
00048 
00049 realvec::~realvec()
00050 {
00051   delete[] data_;
00052   data_ = NULL;
00053 }
00054 
00055 realvec::realvec(mrs_natural size)
00056   : size_(size),
00057     allocatedSize_(0),
00058     data_(NULL),
00059     rows_(1),
00060     cols_(size)
00061 {
00062   allocateData(size_);
00063 }
00064 
00065 realvec::realvec(mrs_natural rows, mrs_natural cols, mrs_real value):
00066   size_(rows * cols),
00067   allocatedSize_(rows * cols),
00068   data_(0),
00069   rows_(rows),
00070   cols_(cols)
00071 {
00072   if (size_ > 0)
00073   {
00074     data_ = new mrs_real[size_];
00075     for (mrs_natural i=0; i<size_; ++i)
00076       data_[i] = value;
00077   }
00078 }
00079 
00080 realvec::realvec(const realvec& a)
00081   : size_(a.size_),
00082     allocatedSize_(0),
00083     data_(NULL),
00084     rows_(a.rows_),
00085     cols_(a.cols_)
00086 {
00087   allocateData(size_);
00088 
00089   for (mrs_natural i=0; i<size_; ++i)
00090     data_[i] = a.data_[i];
00091 
00092 }
00093 
00094 realvec&
00095 realvec::operator=(const realvec& a)
00096 {
00097   if (this != &a)
00098   {
00099     size_ = a.size_;
00100     rows_ = a.rows_;
00101     cols_ = a.cols_;
00102 
00103     //check if we need to allocate more memory
00104     if (allocatedSize_ < size_)
00105       allocateData(size_);
00106 
00107     //copy data
00108     // for (mrs_natural i=0; i < size_; ++i)
00109     // 'data_[i] = a.data_[i];
00110 
00111     memcpy (data_, a.data_, sizeof(mrs_real)*size_);
00112 
00113   }
00114 
00115   return *this;
00116 }
00117 
00118 mrs_real *
00119 realvec::getData() const
00120 {
00121   return data_;
00122 }
00123 
00124 void
00125 realvec::appendRealvec(const realvec newValues)
00126 {
00127   mrs_natural origSize = size_;
00128 
00129   stretch(origSize + newValues.getSize());
00130 
00131   for (mrs_natural i=0; i<newValues.getSize(); ++i)
00132     data_[origSize + i] = newValues.data_[i];
00133 }
00134 
00135 realvec
00136 realvec::getSubVector(mrs_natural startPos, mrs_natural length) const
00137 {
00138   realvec subVector(length);
00139   for (mrs_natural i=0; i<length; ++i)
00140     subVector.data_[i] = data_[startPos + i];
00141   return subVector;
00142 }
00143 
00144 void
00145 realvec::transpose()
00146 {
00147   mrs_real *tmp_ = new mrs_real[size_];
00148 
00149   for (mrs_natural i=0; i < rows_; ++i)
00150     for (mrs_natural j=0; j < cols_; j++)
00151       tmp_[i * cols_ + j] = data_[j * rows_ + i];
00152 
00153   mrs_natural tmp = rows_;
00154   rows_ = cols_;
00155   cols_ = tmp;
00156 
00157   delete [] data_;
00158   data_ = tmp_;
00159 }
00160 
00161 void
00162 realvec::fliplr()
00163 {
00164   mrs_natural i,j,jtmp;
00165 
00166   for (i=0; i < rows_; i++)
00167     for (j=0, jtmp = cols_-1; j < cols_/2; j++, jtmp--)
00168     {
00169       mrs_real tmp  = (*this)(i,j);
00170       (*this)(i,j)  = (*this)(i,jtmp);
00171       (*this)(i,jtmp)   = tmp;
00172     }
00173 }
00174 
00175 void
00176 realvec::flipud()
00177 {
00178   mrs_natural i,j,itmp;
00179 
00180   for (i=0, itmp = rows_-1; i < rows_/2; i++, itmp--)
00181     for (j=0; j < cols_; j++)
00182     {
00183       mrs_real tmp  = (*this)(i,j);
00184       (*this)(i,j)  = (*this)(itmp,j);
00185       (*this)(itmp,j)   = tmp;
00186     }
00187 }
00188 
00189 mrs_real
00190 realvec::median() const
00191 {
00192   if (size_ == 0)
00193     return 0;
00194   realvec tmp(*this);
00195   mrs_real *tmpData = tmp.data_;
00196   std::sort(tmpData, tmpData+size_);
00197   return tmpData[size_/2];
00198 }
00199 
00200 mrs_real
00201 realvec::mean() const
00202 {
00203   mrs_real sum = 0.0;
00204   for (mrs_natural i=0; i < size_; ++i)
00205   {
00206     sum += data_[i];
00207   }
00208   if (sum != 0.0) sum /= size_;
00209   return sum;
00210 }
00211 
00212 void
00213 realvec::sort() //non-descending sort - assumes one dimensional vector only!
00214 {
00215   std::sort(data_, data_+size_);
00216 }
00217 
00218 mrs_real
00219 realvec::sum() const
00220 {
00221   mrs_real sum = 0.0;
00222   for (mrs_natural i=0; i < size_; ++i)
00223   {
00224     sum += data_[i];
00225   }
00226   return sum;
00227 }
00228 
00229 
00230 mrs_real
00231 realvec::var() const
00232 {
00233   mrs_real sum = 0.0;
00234   mrs_real sum_sq = 0.0;
00235   mrs_real val;
00236   mrs_real var;
00237 
00238   for (mrs_natural i=0; i < size_; ++i)
00239   {
00240     val = data_[i];
00241     sum += val;
00242     sum_sq += (val * val);
00243   }
00244   if (sum != 0.0) sum /= size_;
00245   if (sum_sq != 0.0) sum_sq /= size_;
00246 
00247   var = sum_sq - sum * sum;
00248 
00249   if (var < 0.0)
00250     var = 0.0;
00251   return var;
00252 }
00253 
00254 mrs_real
00255 realvec::std() const
00256 {
00257   mrs_real vr = var();
00258   if (vr != 0)
00259     return sqrt(vr);
00260   else
00261     return 0.0;
00262 }
00263 
00264 mrs_natural
00265 realvec::getRows() const
00266 {
00267   return rows_;
00268 }
00269 
00270 mrs_natural
00271 realvec::getCols() const
00272 {
00273   return cols_;
00274 }
00275 
00276 mrs_natural
00277 realvec::getSize() const
00278 {
00279   return size_;
00280 }
00281 
00282 void
00283 realvec::debug_info()
00284 {
00285   MRSERR("realvec information");
00286   MRSERR("size = " << size_);
00287 }
00288 
00289 /* keep the old data and possibly extend */
00290 void
00291 realvec::stretch(mrs_natural size)
00292 {
00293   if (size_ == size)
00294     return; //no need for more memory allocation
00295 
00296   if (size < allocatedSize_)
00297   {
00298     size_ = size;
00299     rows_ = 1;
00300     cols_ = size_;
00301     return; //no need for more memory allocation
00302   }
00303 
00304   // we need more memory allocation!
00305   // size should be always >= 1 at this point
00306   mrs_real* ndata = new mrs_real[size];
00307 
00308   mrs_natural i=0;
00309   for (; i < size_; ++i)
00310     ndata[i] = data_[i]; //copy existing data
00311   for(; i < size; ++i)
00312     ndata[i] = 0.0; //zero new data
00313 
00314   //deallocate existing memory...
00315   delete [] data_;
00316   //...and point to new data memory (if any - it can be NULL, when size == 0)
00317   data_ = ndata;
00318 
00319   //update internal parameters
00320   size_ = size;
00321   allocatedSize_ = size;
00322   rows_ = 1;
00323   cols_ = size_;
00324 }
00325 
00326 /* keep the old data and possibly extend */
00327 void
00328 realvec::stretch(mrs_natural rows, mrs_natural cols)
00329 {
00330   mrs_natural size = rows * cols;
00331 
00332   if ((rows == rows_)&&(cols == cols_))
00333     return;
00334 
00335   if(size == 0) {
00336     size_ = size;
00337     rows_ = rows;
00338     cols_ = cols;
00339     return;
00340   }
00341 
00342     /*lgmartins: this messes up the internal data of the realvec!!
00343     //FIXME: for this optimization to work on 2D realvecs, data has to be
00344     //reorganized inside the realvec...
00345   if (size < allocatedSize_)
00346   {
00347     size_ = size;
00348     rows_ = rows;
00349     cols_ = cols;
00350     return; //no need for more memory allocation
00351   }
00352      */
00353 
00354   mrs_real *ndata = new mrs_real[size];
00355 
00356   // copy and zero new data
00357   for (mrs_natural r=0; r < rows; ++r)
00358     for (mrs_natural c = 0; c < cols; ++c)
00359     {
00360       if ((r < rows_)&&(c < cols_))
00361         ndata[c * rows + r] = data_[c * rows_ + r];
00362       else
00363         ndata[c * rows + r] = 0.0;
00364     }
00365 
00366   delete [] data_;
00367   data_ = ndata;
00368 
00369   size_ = size;
00370   allocatedSize_ = size;
00371   rows_ = rows;
00372   cols_ = cols;
00373 }
00374 
00375 void
00376 realvec::stretchWrite(const mrs_natural pos, const mrs_real val)
00377 {
00378   // don't forget the zero-indexing position!
00379   //   i.e.  pos=0  is the first value to store
00380   mrs_natural wantSize = pos+1;
00381   if (wantSize > size_) {
00382     if ( wantSize < 2*size_ )
00383       // grow exponentially with sequential access
00384       stretch( 2*size_ );
00385     else
00386       // if we have a sudden large value, don't double it
00387       stretch( wantSize );
00388   }
00389   // FIXME: add a MRSASSERT here for debugging.
00390   // cout<<"List stretched to "<<size_<<endl;
00391   data_[pos] = val;
00392 }
00393 
00394 void
00395 realvec::stretchWrite(const mrs_natural r, const mrs_natural c, const mrs_real val)
00396 {
00397   mrs_natural nextR = rows_;
00398   mrs_natural nextC = cols_;
00399   mrs_natural wantR = r+1;
00400   mrs_natural wantC = c+1;
00401   if ( (wantR >= rows_) || (wantC >= cols_) )
00402   {
00403     if ( wantR >= rows_ ) {
00404       if ( wantR < 2*rows_ ) {
00405         nextR = 2*rows_;
00406       } else {
00407         nextR = wantR;
00408       }
00409     }
00410 
00411     if ( wantC >= cols_ ) {
00412       if ( wantC < 2*cols_ ) {
00413         nextC = 2*cols_;
00414       } else {
00415         nextC = wantC;
00416       }
00417     }
00418 
00419     stretch( nextR, nextC );
00420     // FIXME: add a MRSASSERT here for debugging.
00421     // cout<<"List stretched to "<<rows_<<", "<<cols_<<endl;
00422   }
00423   data_[c * rows_ + r] = val;
00424 }
00425 
00426 void
00427 realvec::create(mrs_natural size)
00428 {
00429   size_ = size;
00430   allocateData(size_);
00431   rows_ = 1;
00432   cols_ = size_;
00433 }
00434 
00435 void
00436 realvec::create(mrs_natural rows, mrs_natural cols)
00437 {
00438   size_ = rows * cols;
00439   rows_ = rows;
00440   cols_ = cols;
00441   allocateData(size_);
00442 }
00443 
00444 
00445 void
00446 realvec::create(mrs_real val, mrs_natural rows, mrs_natural cols)
00447 {
00448   size_ = rows * cols;
00449   rows_ = rows;
00450   cols_ = cols;
00451   delete [] data_;
00452   data_ = NULL;
00453   if (size_ > 0)
00454     data_ = new mrs_real[size_];
00455   for (mrs_natural i=0; i<size_; ++i)
00456     data_[i] = val;
00457   allocatedSize_ = size_;
00458 }
00459 
00460 
00461 
00462 void
00463 realvec::allocate(mrs_natural size)
00464 {
00465   delete [] data_;
00466   data_ = NULL;
00467   size_ = size;
00468   cols_ = size_;
00469   rows_ = 1;
00470   allocatedSize_ = size;
00471   if (size_ > 0)
00472     data_ = new mrs_real[size_];
00473 }
00474 
00475 void
00476 realvec::allocate(mrs_natural rows, mrs_natural cols)
00477 {
00478   delete [] data_;
00479   data_ = NULL;
00480   size_ = rows*cols;
00481   cols_ = cols;
00482   rows_ = rows;
00483   allocatedSize_ = size_;
00484   if (size_ > 0)
00485     data_ = new mrs_real[size_];
00486 }
00487 
00488 
00489 void
00490 realvec::setval(mrs_natural start, mrs_natural end, mrs_real val)
00491 {
00492   MRSASSERT(start >= (mrs_natural)0);
00493   MRSASSERT(start < (mrs_natural)size_);
00494   MRSASSERT(end < (mrs_natural)size_);
00495 
00496   for (mrs_natural i=start; i<end; ++i)
00497   {
00498     data_[i] = val;
00499   }
00500 }
00501 
00502 void
00503 realvec::apply(mrs_real (*func) (mrs_real))
00504 {
00505   for (mrs_natural i=0; i<size_; ++i)
00506   {
00507     data_[i] = func(data_[i]);
00508   }
00509 }
00510 
00511 void
00512 realvec::setval(mrs_real val)
00513 {
00514   for (mrs_natural i=0; i<size_; ++i)
00515   {
00516     data_[i] = val;
00517   }
00518 }
00519 
00520 void
00521 realvec::abs()
00522 {
00523   for (mrs_natural i=0; i<size_; ++i)
00524   {
00525     data_[i] = fabs(data_[i]);
00526   }
00527 }
00528 
00529 void
00530 realvec::pow (mrs_real exp)
00531 {
00532   for (mrs_natural i=0; i<size_; i++)
00533   {
00534     data_[i] = ::pow (data_[i], exp);
00535   }
00536 }
00537 
00538 void
00539 realvec::norm()
00540 {
00541   mrs_real mean = this->mean();
00542   mrs_real std = this->std();
00543   for (mrs_natural i=0; i < size_; ++i)
00544   {
00545     data_[i] = (data_[i] - mean) / std;
00546   }
00547 }
00548 
00549 
00550 void
00551 realvec::normMaxMin()
00552 {
00553   mrs_real max = DBL_MIN;
00554   mrs_real min = DBL_MAX;
00555 
00556   for (mrs_natural i=0; i < size_; ++i)
00557   {
00558     if (data_[i] > max)
00559       max = data_[i];
00560     if (data_[i] < min)
00561       min = data_[i];
00562   }
00563 
00564 
00565   for (mrs_natural i=0; i < size_; ++i)
00566   {
00567     data_[i] = (data_[i] - min) / (max - min);
00568   }
00569 
00570 
00571 }
00572 
00573 
00574 void
00575 realvec::norm(mrs_real mean, mrs_real std)
00576 {
00577   for (mrs_natural i=0; i < size_; ++i)
00578   {
00579     data_[i] = (data_[i] - mean) / std;
00580   }
00581 }
00582 
00583 void
00584 realvec::renorm(mrs_real old_mean, mrs_real old_std, mrs_real new_mean, mrs_real new_std)
00585 {
00586   for (mrs_natural i=0; i < size_; ++i)
00587   {
00588     data_[i] = (data_[i] - old_mean) / old_std;
00589     data_[i] *= new_std;
00590     data_[i] += new_mean;
00591   }
00592 }
00593 
00594 mrs_natural
00595 realvec::invert(realvec& res)
00596 {
00597   if (rows_ != cols_)
00598   {
00599     MRSERR("realvec::invert() - matrix should be square!");
00600     res.create(0);
00601     return -1;
00602   }
00603   if (this != &res)
00604   {
00605     mrs_natural rank;
00606     mrs_natural r,c,i;
00607     mrs_real temp;
00608 
00609     res.create(rows_, cols_);
00610 
00611     rank = 0;
00612     for (r = 0; r < rows_; ++r)
00613       for (c=0; c < cols_; ++c)
00614       {
00615         if (r == c)
00616           res(r,c) = 1.0;
00617         else
00618           res(r,c) = 0.0;
00619       }
00620     for (i = 0; i < rows_; ++i)
00621     {
00622       if ((*this)(i,i) == 0)
00623       {
00624         for (r = i; r < rows_; ++r)
00625           for (c = 0; c < cols_; ++c)
00626           {
00627             (*this)(i,c) += (*this)(r,c);
00628             res(i,c) += res(r,c);
00629           }
00630       }
00631       for (r = i; r < rows_; ++r)
00632       {
00633         temp = (*this)(r,i);
00634         if (temp != 0)
00635           for (c =0; c < cols_; ++c)
00636           {
00637             (*this)(r,c) /= temp;
00638             res(r,c) /= temp;
00639           }
00640       }
00641       if (i != rows_-1)
00642       {
00643         for (r = i+1; r < rows_; ++r)
00644         {
00645           temp = (*this)(r,i);
00646           if (temp != 0.0)
00647             for (c=0; c < cols_; ++c)
00648             {
00649               (*this)(r,c) -= (*this)(i,c);
00650               res(r,c) -= res(i,c);
00651             }
00652         }
00653       }
00654     }
00655     for (i=1; i < rows_; ++i)
00656       for (r=0; r < i; ++r)
00657       {
00658         temp = (*this)(r,i);
00659         for (c=0; c < cols_; ++c)
00660         {
00661           (*this)(r,c) -= (temp * (*this)(i,c));
00662           res(r,c) -= (temp * res(i,c));
00663         }
00664       }
00665     for (r =0; r < rows_; ++r)
00666       for (c=0; c < cols_; ++c)
00667         (*this)(r,c) = res(r,c);
00668     return rank;
00669   }
00670   else
00671   {
00672     res.create(0);
00673     MRSERR("realvec::invert() - inPlace operation not supported - returning empty result vector!");
00674     return -1;
00675   }
00676 }
00677 
00678 void
00679 realvec::sqr()
00680 {
00681   for (mrs_natural i=0; i<size_; ++i)
00682   {
00683     data_[i] *= data_[i];
00684   }
00685 }
00686 
00687 mrs_natural
00688 realvec::search(mrs_real val)
00689 {
00690   mrs_real minDiff = MAXREAL;
00691   mrs_natural index=-1;
00692   for (mrs_natural i=0; i<size_; ++i)
00693     if (fabs(data_[i]-val)< minDiff)
00694     {
00695       minDiff = fabs(data_[i]-val);
00696       index=i;
00697     }
00698   return index;
00699 }
00700 
00701 void
00702 realvec::sqroot()
00703 {
00704   for (mrs_natural i=0; i<size_; ++i)
00705   {
00706     data_[i] = sqrt(data_[i]);
00707   }
00708 }
00709 
00710 void
00711 realvec::send(Communicator *com)
00712 {
00713   static char *buf = new char[256];
00714   mrs_string message;
00715   sprintf(buf, "%i\n", (int)size_);
00716   message = buf;
00717   com->send_message(message);
00718   for (mrs_natural i=0; i<size_; ++i)
00719   {
00720     sprintf(buf, "%f\n", data_[i]);
00721     message = buf;
00722     com->send_message(message);
00723   }
00724   MRSERR("realvec::send numbers were sent to the client");
00725 
00726 }
00727 
00728 bool
00729 realvec::read(mrs_string filename)
00730 {
00731   ifstream from(filename.c_str());
00732   if (from.is_open())
00733   {
00734     from >> (*this);
00735     return true;
00736   }
00737   else
00738   {
00739     MRSERR("realvec::read: failed to open file: " << filename);
00740     return false;
00741   }
00742 }
00743 
00744 
00745 void
00746 realvec::shuffle()
00747 {
00748   // Use a Fisher-Yates shuffle : http://en.wikipedia.org/wiki/Fisher-Yates_shuffle
00749   unsigned int n = cols_;
00750   while (n > 1)
00751   {
00752     // Generate a random index in the range [0, n).
00753     unsigned int k = (unsigned int)((mrs_real)n * (mrs_real)rand() / (mrs_real)(RAND_MAX + 1.0));
00754 
00755     n--;
00756 
00757     // Swap column n and column k.
00758     if (k != n)
00759     {
00760       for (int r = 0; r < rows_; ++r)
00761         swap(data_[k * rows_ + r], data_[n * rows_ + r]);
00762     }
00763   }
00764 }
00765 
00766 
00767 bool
00768 realvec::write(mrs_string filename) const
00769 {
00770   ofstream os(filename.c_str());
00771   if (os.is_open())
00772   {
00773     os << (*this) << endl;
00774     return true;
00775   }
00776   else
00777   {
00778     MRSERR("realvec::write: failed to open file to write: filename");
00779     return false;
00780   }
00781 }
00782 
00783 void
00784 realvec::dump()
00785 {
00786   for (mrs_natural i =0 ; i< size_ ; ++i)
00787     MRSMSG(data_[i] << " ") ;
00788   MRSMSG(endl);
00789 }
00790 
00791 bool
00792 realvec::readText(mrs_string filename)
00793 {
00794   ifstream infile(filename.c_str());
00795   if (infile.is_open())
00796   {
00797     int i = 0;
00798     if (size_ == 0)
00799       create(1);
00800 
00801     mrs_real value;
00802     while (infile >> value)
00803     {
00804       // slower but safer
00805       stretchWrite(i,value);
00806       ++i;
00807     }
00808     stretch(i-1);
00809     infile.close();
00810     return true;
00811   }
00812   else
00813   {
00814     MRSERR("realvec::readText: failed to open file: " << filename);
00815     return false;
00816   }
00817 }
00818 
00819 bool
00820 realvec::writeText(mrs_string filename)
00821 {
00822   if (size_ == 0)
00823     return true; //[?]
00824 
00825   ofstream outfile(filename.c_str());
00826   if (outfile.is_open())
00827   {
00828     for (mrs_natural i=0; i<size_; ++i)
00829     {
00830       outfile << data_[i] <<endl;
00831     }
00832     outfile.close();
00833     return true;
00834   }
00835   else
00836   {
00837     MRSERR("realvec::writeText: failed to open file: " << filename);
00838     return false;
00839   }
00840 }
00841 
00846 void
00847 realvec::dumpDataOnly(std::ostream& o, std::string columnSep, std::string rowSep) const
00848 {
00849   for (mrs_natural r = 0; r < this->rows_; ++r)
00850   {
00851     for (mrs_natural c = 0; c < this->cols_; ++c)
00852     {
00853       o << this->data_[c * this->rows_ + r];
00854       if (c < this->cols_ - 1)
00855       {
00856         o << columnSep ;
00857       }
00858     }
00859     if (r < this->rows_ - 1)
00860     {
00861       o << rowSep;
00862     }
00863   }
00864 }
00865 
00866 ostream&
00867 operator<< (ostream& o, const realvec& vec)
00868 {
00869   o << "# MARSYAS mrs_realvec" << endl;
00870   o << "# Size = " << vec.size_ << endl << endl;
00871   o << endl;
00872 
00873   o << "# type: matrix" << endl;
00874   o << "# rows: " << vec.rows_ << endl;
00875   o << "# columns: " << vec.cols_ << endl;
00876 
00877   vec.dumpDataOnly(o, " ", "\n");
00878   o << endl;
00879 
00880   // for (mrs_natural i=0; i<vec.size_; ++i)
00881   // o << vec.data_[i] << endl;
00882   o << endl;
00883   o << "# Size = " << vec.size_ << endl;
00884   o << "# MARSYAS mrs_realvec" << endl;
00885   return o;
00886 }
00887 
00888 istream&
00889 operator>>(istream& is, realvec& vec)
00890 {
00891   // [WTF] ... is this necessary?  doesn't allocate() delete the data array, and reallocate? (Jen)
00892   //  if (vec.size_ != 0)
00893   //    {
00894   //      MRSERR("realvec::operator>>: realvec already allocated cannot read from istream");
00895   //      MRSERR("vec.size_ = " << vec.size_);
00896   //
00897   //      return is;
00898   //    }
00899   mrs_string str0,str1,str2;
00900   mrs_natural size;
00901   mrs_natural i;
00902 
00903   is >> str0;
00904   is >> str1;
00905   is >> str2;
00906 
00907 
00908   if ((str0 != "#") || (str1 != "MARSYAS") || (str2 != "mrs_realvec"))
00909   {
00910     MRSERR("realvec::operator>>: Problem1 reading realvec object from istream");
00911     MRSERR("-str0 = " << str0 << endl);
00912     MRSERR("-str1 = " << str1 << endl);
00913     MRSERR("-str2 = " << str2 << endl);
00914     return is;
00915   }
00916   is >> str0;
00917   is >> str1;
00918   is >> str2;
00919 
00920 
00921   if ((str0 != "#") || (str1 != "Size") || (str2 != "="))
00922   {
00923     MRSERR("realvec::operator>>: Problem2 reading realvec object from istream");
00924     MRSERR("-str0 = " << str0 << endl);
00925     MRSERR("-str1 = " << str1 << endl);
00926     MRSERR("-str2 = " << str2 << endl);
00927     return is;
00928   }
00929   is >> size;
00930 
00931   vec.create(size);
00932   for (i=0; i<3; ++i)
00933   {
00934     is >> str0;
00935   }
00936   is >> str0 >> str1 >> vec.rows_;
00937   is >> str0 >> str1 >> vec.cols_;
00938 
00939   for (mrs_natural r = 0; r < vec.rows_; ++r)
00940     for (mrs_natural c = 0; c < vec.cols_; ++c)
00941     {
00942       is >> str0;
00943       if (str0 == "nan") {
00944         vec.data_[c*vec.rows_ + r] = std::numeric_limits<double>::quiet_NaN();
00945       } else {
00946         std::stringstream s(str0);
00947         s >> vec.data_[c * vec.rows_ + r];
00948       }
00949     }
00950 
00951   is >> str0;
00952   is >> str1;
00953   is >> str2;
00954   if ((str0 != "#") || (str1 != "Size") || (str2 != "="))
00955   {
00956     MRSERR("realvec::operator>>: Problem3 reading realvec object from istream");
00957     MRSERR("-str0 = " << str0 << endl);
00958     MRSERR("-str1 = " << str1 << endl);
00959     MRSERR("-str2 = " << str2 << endl);
00960     is >> str0;
00961     is >> str1;
00962     is >> str2;
00963     MRSERR("-str0 = " << str0 << endl);
00964     MRSERR("-str1 = " << str1 << endl);
00965     MRSERR("-str2 = " << str2 << endl);
00966     return is;
00967   }
00968   is >> size;
00969   is >> str0;
00970   is >> str1;
00971   is >> str2;
00972 
00973   if ((str0 != "#") || (str1 != "MARSYAS") || (str2 != "mrs_realvec"))
00974   {
00975     MRSERR("realvec::operator>>: Problem4 reading realvec object from istream");
00976     MRSERR("-str0 = " << str0 << endl);
00977     MRSERR("-str1 = " << str1 << endl);
00978     MRSERR("-str2 = " << str2 << endl);
00979     return is;
00980   }
00981   return is;
00982 }
00983 
00984 realvec
00985 realvec::operator()(std::string r, std::string c)
00986 {
00987   mrs_string::size_type r_l = r.length();
00988   mrs_string::size_type c_l = c.length();
00989 
00990   mrs_string::size_type r_c = r.find(":");
00991   mrs_string::size_type c_c = c.find(":");
00992 
00993   mrs_string::size_type r_a;
00994   mrs_string::size_type r_b;
00995 
00996   mrs_string::size_type c_a;
00997   mrs_string::size_type c_b;
00998 
00999   char *endptr;
01000 
01001   MRSASSERT( (r_c == 0 && r_l == 1) || (r_c == mrs_string::npos) || (r_c>0 && r_l-r_c>1) );
01002   MRSASSERT( (c_c == 0 && c_l == 1) || (c_c == mrs_string::npos) || (c_c>0 && c_l-c_c>1) );
01003 
01004   if ( r_c != mrs_string::npos && r_l > 1 )
01005   {
01006     r_a = strtol( r.substr(0,r_c).c_str() , &endptr , 10  );
01007     MRSASSERT( *endptr == '\0' );
01008     r_b = strtol( r.substr(r_c+1,r_l-r_c-1).c_str() , &endptr , 10  );
01009     MRSASSERT( *endptr == '\0' );
01010   }
01011   else if ( r_c == mrs_string::npos )
01012   {
01013     r_a = r_b = strtol( r.c_str() , &endptr , 10 );
01014     MRSASSERT( *endptr == '\0' );
01015   }
01016   else
01017   {
01018     r_a = 0;
01019     r_b = rows_-1;
01020   }
01021 
01022   MRSASSERT( (mrs_natural)r_b < rows_ );
01023 
01024   if ( c_c != mrs_string::npos && c_l > 1 )
01025   {
01026     c_a = strtol( c.substr(0,c_c).c_str() , &endptr , 10  );
01027     MRSASSERT( *endptr == '\0' );
01028     c_b = strtol( c.substr(c_c+1,c_l-c_c-1).c_str() , &endptr , 10 );
01029     MRSASSERT( *endptr == '\0' );
01030   }
01031   else if ( c_c == mrs_string::npos )
01032   {
01033     c_a = c_b = strtol( c.c_str() , &endptr , 10 );
01034     MRSASSERT( *endptr == '\0' );
01035   }
01036   else
01037   {
01038     c_a = 0;
01039     c_b = cols_-1;
01040   }
01041 
01042   MRSASSERT( (mrs_natural)c_b < cols_ );
01043 
01044   r_l = r_b - r_a + 1;
01045   c_l = c_b - c_a + 1;
01046 
01047   realvec matrix;
01048 
01049   matrix.create( (mrs_natural) r_l , (mrs_natural) c_l );
01050 
01051   for ( c_c = c_a ; c_c <= c_b ; c_c++ )
01052   {
01053     for ( r_c = r_a ; r_c <= r_b ; r_c++ )
01054     {
01055       matrix.data_[(c_c-c_a) * r_l + (r_c-r_a)] = data_[c_c * rows_ + r_c];
01056     }
01057   }
01058 
01059   return matrix;
01060 }
01061 
01062 realvec
01063 realvec::operator()(std::string c)
01064 {
01065   mrs_string::size_type c_l = c.length();
01066   mrs_string::size_type c_c = c.find(":");
01067   mrs_string::size_type c_a;
01068   mrs_string::size_type c_b;
01069   char *endptr;
01070 
01071   MRSASSERT( (c_c == 0 && c_l == 1) || (c_c == mrs_string::npos) || (c_c>0 && c_l-c_c>1) );
01072 
01073   if ( c_c != mrs_string::npos && c_l > 1 )
01074   {
01075     c_a = strtol( c.substr(0,c_c).c_str() , &endptr , 10  );
01076     MRSASSERT( *endptr == '\0' );
01077     c_b = strtol( c.substr(c_c+1,c_l-c_c-1).c_str() , &endptr , 10  );
01078     MRSASSERT( *endptr == '\0' );
01079   }
01080   else if ( c_c == mrs_string::npos )
01081   {
01082     c_a = c_b = strtol( c.c_str() , &endptr , 10 );
01083     MRSASSERT( *endptr == '\0' );
01084   }
01085   else
01086   {
01087     c_a = 0;
01088     c_b = (rows_*cols_)-1;
01089   }
01090 
01091   MRSASSERT( (mrs_natural)c_b < rows_*cols_ );
01092   c_l = c_b - c_a + 1;
01093 
01094   realvec matrix;
01095   matrix.create( (mrs_natural) c_l );
01096 
01097   for ( c_c = c_a ; c_c <= c_b ; c_c++ )
01098   {
01099     matrix.data_[(c_c-c_a)] = data_[c_c];
01100   }
01101   return matrix;
01102 }
01103 
01104 void
01105 realvec::getRow(const mrs_natural r, realvec& res) const
01106 {
01107   if (this != &res)
01108   {
01109     if (r >= rows_ )
01110     {
01111       MRSERR("realvec::getRow() - row index greater than realvec number of rows! Returning empty result vector.");
01112       res.create(0);
01113       return;
01114     }
01115     res.stretch(cols_);
01116     for (mrs_natural c=0; c < cols_; ++c)
01117     {
01118       res(c) = (*this)(r,c);
01119     }
01120   }
01121   else
01122   {
01123     res.create(0);
01124     MRSERR("realvec::getRow() - inPlace operation not supported - returning empty result vector!");
01125   }
01126 }
01127 
01128 void
01129 realvec::getCol(const mrs_natural c, realvec& res) const
01130 {
01131   if (this != &res)
01132   {
01133     if (c >= cols_)
01134     {
01135       MRSERR("realvec::getCol() - row index greater than realvec number of rows! Returning empty result vector.");
01136       res.create(0);
01137       return;
01138     }
01139     res.stretch(rows_,1);
01140     for (mrs_natural r=0; r < rows_; ++r)
01141     {
01142       res(r) = (*this)(r,c);
01143     }
01144   }
01145   else
01146   {
01147     res.create(0);
01148     MRSERR("realvec::getCol() - inPlace operation not supported - returning empty result vector!");
01149   }
01150 }
01151 void
01152 realvec::getSubMatrix (const mrs_natural r, const mrs_natural c, realvec& res)
01153 {
01154   if (this != &res)
01155   {
01156     mrs_natural numRows = res.getRows (),
01157                 numCols = res.getCols ();
01158     //if (c + numCols > cols_ || r + numRows > rows_) this might also be a reasonable check...
01159     if (c >= cols_ || r >= rows_)
01160     {
01161       MRSERR("realvec::getSubMatrix() - index larger than realvec number of rows/cols! Returning empty result vector.");
01162       res.create(0);
01163       return;
01164     }
01165     mrs_natural m,n,i,j;
01166     mrs_natural endRow  = min(rows_, r + numRows),
01167                  endCol = min(cols_, c + numCols);
01168     for (m=r, i=0; m < endRow; m++,++i)
01169     {
01170       for (n=c, j=0; n < endCol; n++, j++)
01171         res(i,j) = (*this)(m,n);
01172     }
01173 
01174     // if there are remaining elements, fill up with zeros (or should we throw something? see MRSERR check above)
01175     for (i=endRow-r; i < numRows; ++i)
01176       for (j=0; j < numCols; j++)
01177         res(i,j)    = 0;
01178     for (j=endCol-c; j < numCols; j++)
01179       for (i=0; i < numRows; ++i)
01180         res(i,j)    = 0;
01181   }
01182   else
01183   {
01184     res.create(0);
01185     MRSERR("realvec::getSubMatrix() - inPlace operation not supported - returning empty result vector!");
01186   }
01187 }
01188 
01189 void
01190 realvec::setRow (const mrs_natural r, const realvec src)
01191 {
01192   setSubMatrix (r,0, src);
01193 }
01194 void
01195 realvec::setCol (const mrs_natural c, const realvec src)
01196 {
01197   setSubMatrix (0,c,src);
01198 }
01199 void
01200 realvec::setSubMatrix (const mrs_natural r, const mrs_natural c, const realvec src)
01201 {
01202   mrs_natural   m,n;
01203   mrs_natural   numRows = src.getRows (),
01204               numCols   = src.getCols ();
01205   if (c+numCols > cols_ || r+numRows > rows_)
01206   {
01207     MRSERR("realvec::setSubMatrix() - dimension mismatch! Abort.");
01208     return;
01209   }
01210 
01211   mrs_natural   endRow  = min(rows_, r + numRows),
01212                endCol   = min(cols_, c + numCols);
01213   for (m=r; m < endRow; m++)
01214   {
01215     for (n=c; n < endCol; n++)
01216       (*this)(m,n)  = src(m-r,n-c);
01217   }
01218 }
01219 
01220 
01221 mrs_real
01222 realvec::maxval(mrs_natural* index) const
01223 {
01224   mrs_real max = numeric_limits<mrs_real>::max() * -1.0;
01225   mrs_natural ind = 0;
01226   for (mrs_natural i=0; i < size_; ++i)
01227   {
01228     if (data_[i] > max)
01229     {
01230       max = data_[i];
01231       ind = i;
01232     }
01233   }
01234   if (index)
01235     *index = ind;
01236   return max;
01237 }
01238 
01239 mrs_real
01240 realvec::minval() const
01241 {
01242   mrs_real min = numeric_limits<mrs_real>::max();
01243   for (mrs_natural i=0; i < size_; ++i)
01244   {
01245     if (data_[i] < min)
01246       min = data_[i];
01247   }
01248   return min;
01249 }
01250 
01251 void
01252 realvec::meanObs(realvec& res) const
01253 {
01254   if (this != &res)
01255   {
01256     realvec obsrow(cols_); //row vector //[TODO]
01257     res.stretch(rows_, 1); //column vector
01258 
01259     for (mrs_natural r=0; r<rows_; ++r)
01260     {
01261       //obsrow = getRow(r);
01262       for (mrs_natural c=0; c<cols_; ++c)
01263       {
01264         obsrow(c) = (*this)(r,c); //get observation row
01265       }
01266       res(r,0) = obsrow.mean();
01267     }
01268   }
01269   else
01270   {
01271     res.create(0);
01272     MRSERR("realvec::meanObs() - inPlace operation not supported - returning empty result vector!");
01273   }
01274 }
01275 
01276 void
01277 realvec::varObs(realvec& res) const
01278 {
01279   if (this != &res)
01280   {
01281     res.create(rows_, 1); //column vector
01282     realvec obsrow(cols_); //row vector //[TODO]
01283 
01284     for (mrs_natural r=0; r<rows_; ++r)
01285     {
01286       //obsrow = getRow(r);
01287       for (mrs_natural c=0; c<cols_; ++c)
01288       {
01289         obsrow(c) = (*this)(r,c); //get observation row
01290       }
01291       res(r,0) = obsrow.var();
01292     }
01293   }
01294   else
01295   {
01296     res.create(0);
01297     MRSERR("realvec::varObs() - inPlace operation not supported - returning empty result vector!");
01298   }
01299 }
01300 
01301 void
01302 realvec::stdObs(realvec& res) const
01303 {
01304   if (this != &res)
01305   {
01306     realvec obsrow(cols_); //row vector //[TODO]
01307     res.stretch(rows_, 1); //column vector
01308     for (mrs_natural r=0; r<rows_; ++r)
01309     {
01310       //obsrow = getRow(r);
01311       for (mrs_natural c=0; c < cols_; ++c)
01312       {
01313         obsrow(c) = (*this)(r,c); //get observation row
01314       }
01315       res(r,0) = obsrow.std();
01316     }
01317   }
01318   else
01319   {
01320     res.create(0);
01321     MRSERR("realvec::stdObs() - inPlace operation not supported - returning empty result vector!");
01322   }
01323 }
01324 
01325 void
01326 realvec::normObs()
01327 {
01328   //normalize (aka standardize) matrix
01329   //using observations means and standard deviations
01330   realvec obsrow(cols_); //[TODO]
01331   mrs_real mean, std;
01332 
01333   for (mrs_natural r=0; r < rows_; ++r)
01334   {
01335     for (mrs_natural c=0; c < cols_; ++c)
01336     {
01337       obsrow(c) = (*this)(r,c); //get observation row
01338     }
01339     //obsrow = getRow(r);
01340     mean = obsrow.mean();
01341     std = obsrow.std();
01342     for (mrs_natural c=0; c < cols_; ++c)
01343     {
01344       (*this)(r,c) -= mean;
01345       (*this)(r,c) /= std;
01346     }
01347   }
01348 }
01349 
01350 void
01351 realvec::normObsMinMax()
01352 {
01353   //normalize (aka standardize) matrix
01354   //using observations min and max values
01355 
01356   realvec obsrow(cols_);
01357   mrs_real min, max, dif;
01358 
01359   for (mrs_natural r=0; r < rows_; ++r)
01360   {
01361     //for (mrs_natural c=0; c < cols_; ++c)
01362     //{
01363     //  obsrow(c) = (*this)(r,c); //get observation row
01364     //}
01365     getRow(r, obsrow); //[TODO]
01366     min = obsrow.minval();
01367     max = obsrow.maxval();
01368     dif = max-min;
01369     if (dif ==0)
01370       dif=1.0;
01371     for (mrs_natural c=0; c < cols_; ++c)
01372     {
01373       (*this)(r,c) -= min;
01374       (*this)(r,c) /= dif;
01375     }
01376   }
01377 }
01378 
01379 void
01380 realvec::normSpl(mrs_natural index)
01381 {
01382   //normalize (aka standardize) matrix
01383   //using ???       [?] [TODO]
01384   mrs_real mean;
01385   mrs_real std;
01386 
01387   realvec colVec;//[TODO]
01388 
01389   if (!index)
01390     index=cols_;
01391   for (mrs_natural r=0; r < index; ++r)
01392   {
01393     //for (mrs_natural c=0; c < cols_; ++c)
01394     //{
01395     //  obsrow(c) = (*this)(r,c); //get observation row
01396     //}
01397 
01398     getCol(r, colVec);//[TODO]
01399     mean = colVec.mean();
01400     std = colVec.std();
01401 
01402     if (std)
01403       for (mrs_natural c=0; c < rows_; ++c)
01404       {
01405         (*this)(c, r) -= mean;
01406         (*this)(c, r) /= std;
01407       }
01408   }
01409 }
01410 
01411 void
01412 realvec::normSplMinMax(mrs_natural index)
01413 {
01414   //normalize (aka standardize) matrix
01415   //using ???       [?] [TODO]
01416   mrs_real min;
01417   mrs_real max, dif;
01418   realvec colVec; //[TODO]
01419 
01420   if (!index)
01421     index=cols_;
01422   for (mrs_natural r=0; r < index; ++r)
01423   {
01424     //for (mrs_natural c=0; c < cols_; ++c)
01425     //{
01426     //  obsrow(c) = (*this)(r,c); //get observation row
01427     //}
01428     getCol(r, colVec);
01429     min = colVec.minval(); //[TODO]
01430     max = colVec.maxval(); //[TODO]
01431     dif = max-min;
01432     if (dif ==0)
01433       dif=1.0;
01434     if (max)
01435       for (mrs_natural c=0; c < rows_; ++c)
01436       {
01437         (*this)(c, r) -= min;
01438         (*this)(c, r) /= dif;
01439       }
01440   }
01441 }
01442 
01443 void
01444 realvec::correlation(realvec& res) const
01445 {
01446   //correlation as computed in Marsyas0.1
01447   //computes the correlation between observations
01448   //Assumes data points (i.e. examples) in columns and features (i.e. observations) in rows (as usual in Marsyas0.2).
01449   if (size_ == 0)
01450   {
01451     MRSERR("realvec::correlation() : empty input matrix! returning empty correlation matrix!");
01452     res.create(0);
01453     return;
01454   }
01455   if (this != &res)
01456   {
01457     res.stretch(rows_, rows_);//correlation matrix
01458     // normalize observations (i.e subtract obs. mean, divide by obs. standard dev)
01459     realvec temp = (*this); //[TODO]
01460     temp.normObs();
01461 
01462     mrs_real sum = 0.0;
01463     for (mrs_natural r1=0; r1< rows_; ++r1)
01464     {
01465       for (mrs_natural r2=0; r2 < rows_; ++r2)
01466       {
01467         sum = 0.0;
01468 
01469         for (mrs_natural c=0; c < cols_; ++c)
01470           sum += (temp(r1,c) * temp(r2,c));
01471 
01472         sum /= cols_;
01473         res(r1,r2) = sum;
01474       }
01475     }
01476   }
01477   else
01478   {
01479     res.create(0);
01480     MRSERR("realvec::correlation() - inPlace operation not supported - returning empty result vector!");
01481   }
01482 }
01483 
01484 void
01485 realvec::covariance(realvec& res) const
01486 {
01487   //computes the (unbiased estimate) covariance between observations (as in MATLAB cov()).
01488   //Assumes data points (i.e. examples) in columns and features (i.e. observations) in rows (as usual in Marsyas0.2).
01489   //This method assumes non-standardized data (typical covariance calculation).
01490   if (size_ == 0)
01491   {
01492     MRSERR("realvec::covariance(): empty input matrix! returning empty covariance matrix!");
01493     res.create(0);
01494     return;
01495   }
01496 
01497   if (this != &res)
01498   {
01499     res.stretch(rows_, rows_); //covariance matrix
01500     //check if there are sufficient data points for a good covariance estimation...
01501     if (cols_ < (rows_ + 1))
01502     {
01503       MRSWARN("realvec::covariance() : nr. data points < nr. observations + 1 => covariance matrix is SINGULAR!");
01504     }
01505     if ( (mrs_real)cols_ < ((mrs_real)rows_*(mrs_real)(rows_-1.0)/2.0))
01506     {
01507       MRSWARN("realvec::covariance() : too few data points => ill-calculation of covariance matrix!");
01508     }
01509 
01510     realvec meanobs;
01511     this->meanObs(meanobs);//observation means //[TODO]
01512 
01513     mrs_real sum = 0.0;
01514     for (mrs_natural r1=0; r1< rows_; ++r1)
01515     {
01516       for (mrs_natural r2=0; r2 < rows_; ++r2)
01517       {
01518         sum = 0.0;
01519         for (mrs_natural c=0; c < cols_; ++c)
01520           sum += ((data_[c * rows_ + r1] - meanobs(r1)) * (data_[c * rows_ + r2]- meanobs(r2)));
01521 
01522         if (cols_ > 1)
01523           sum /= (cols_ - 1);//unbiased covariance estimate
01524         else
01525           sum /= cols_; //biased covariance estimate
01526 
01527         res(r1,r2) = sum;
01528       }
01529     }
01530   }
01531   else
01532   {
01533     res.create(0);
01534     MRSERR("realvec::covariance() - inPlace operation not supported - returning empty result vector!");
01535   }
01536 }
01537 
01538 void
01539 realvec::covariance2(realvec& res) const
01540 {
01541   //covariance matrix as computed in Marsyas0.1
01542   //computes the covariance between observations.
01543   //Assumes data points (i.e. examples) in columns and features (i.e. observations) in rows (as usual in Marsyas0.2).
01544   //This calculation assumes standardized data at its input...
01545   if (size_ == 0)
01546   {
01547     MRSERR("realvec::covariance() : empty input matrix! returning empty and invalid covariance matrix!");
01548     res.create(0);
01549     return;
01550   }
01551   if (this != &res)
01552   {
01553     res.stretch(rows_, rows_); //covariance matrix
01554     //check if there are sufficient data points for a good covariance estimation...
01555     if (cols_ < (rows_ + 1))
01556     {
01557       MRSWARN("realvec::covariance() : nr. data points < nr. observations + 1 => covariance matrix is SINGULAR!");
01558     }
01559     if ( (mrs_real)cols_ < ((mrs_real)rows_*(mrs_real)(rows_-1.0)/2.0))
01560     {
01561       MRSWARN("realvec::covariance() : too few data points => ill-calculation of covariance matrix!");
01562     }
01563 
01564     for (mrs_natural r1=0; r1< rows_; ++r1)
01565     {
01566       for (mrs_natural r2=0; r2 < rows_; ++r2)
01567       {
01568         mrs_real sum(0);
01569 
01570         for (mrs_natural c=0; c < cols_; ++c)
01571           sum += (data_[c * rows_ + r1] * data_[c * rows_ + r2]);
01572 
01573         sum /= cols_;
01574         res(r1,r2) = sum;
01575       }
01576     }
01577   }
01578   else
01579   {
01580     res.create(0);
01581     MRSERR("realvec::covariance2() - inPlace operation not supported - returning empty result vector!");
01582   }
01583 }
01584 
01585 mrs_real
01586 realvec::trace() const
01587 {
01588   if (cols_ != rows_)
01589   {
01590     MRSWARN("realvec::trace() - matrix is not square!");
01591   }
01592 
01593   mrs_real res(0);
01594 
01595   for (mrs_natural i = 0; i < size_;)
01596   {
01597     res += data_[i];
01598     i += cols_+1;
01599   }
01600 
01601   return res;
01602 }
01603 
01604 mrs_real
01605 realvec::det() const
01606 {
01607   NumericLib numlib;
01608   return numlib.determinant(*this);
01609 }
01610 
01611 
01612 
01613 // allocate data and initialize to zero
01614 // minimum of 1 element is allocated to prevent null pointers
01615 // not that the size_ of the object will stay 0 just like cols_
01616 // Allocating data doen't have anything to do with the size count   of the object.
01617 void realvec::allocateData(mrs_natural size)
01618 {
01619   // if data is not null delete it
01620   delete[] data_;
01621   data_ = NULL;
01622 
01623   allocatedSize_ = size > 0 ? size : 1; // minimum of 1
01624   data_ = new mrs_real[allocatedSize_];
01625 
01626   for (mrs_natural i=0; i<allocatedSize_; ++i)
01627     data_[i] = 0.0;
01628 }
01629 
01630 
01632 
01634 
01635 
01636 
01637 
01638 
01639 
01640 bool realvec::operator==(const realvec &v1) const
01641 {
01642   //different size -> not equal
01643   if (v1.getRows()!=rows_) return false;
01644   if (v1.getCols()!=cols_) return false;
01645 
01646   for(int r=0; r<v1.getRows(); ++r)
01647   {
01648     for(int c=0; c<v1.getCols(); ++c)
01649     {
01650       if(v1(r,c)!=data_[c * rows_ + r])return false;
01651     }
01652   }
01653   return true;
01654 }
01655 
01656 
01657 bool realvec::operator!=(const realvec &v1)  const
01658 {
01659   return !(*this == v1);
01660 }
01661 
01662 realvec&
01663 realvec::operator+=(const realvec& rhs)
01664 {
01665   if (size_ != rhs.size_)
01666     throw std::runtime_error("realvec: Trying to sum matrices of incompatible size.");
01667 
01668   for (mrs_natural i=0; i<size_; ++i)
01669     data_[i] += rhs.data_[i];
01670   return *this;
01671 }
01672 
01673 
01674 realvec&
01675 realvec::operator-=(const realvec& rhs)
01676 {
01677   if (size_ != rhs.size_)
01678     throw std::runtime_error("realvec: Trying to subtract matrices of incompatible size.");
01679 
01680   for (mrs_natural i=0; i<size_; ++i)
01681     data_[i] -= rhs.data_[i];
01682   return *this;
01683 }
01684 
01685 
01686 realvec&
01687 realvec::operator*=(const realvec& rhs)
01688 {
01689   if (size_ != rhs.size_)
01690     throw std::runtime_error("realvec: Trying to multiply matrices of incompatible size.");
01691 
01692   for (mrs_natural i=0; i<size_; ++i)
01693     data_[i] *= rhs.data_[i];
01694   return *this;
01695 }
01696 
01697 
01698 realvec&
01699 realvec::operator/=(const realvec& rhs)
01700 {
01701   if (size_ != rhs.size_)
01702     throw std::runtime_error("realvec: Trying to divide matrices of incompatible size.");
01703 
01704   for (mrs_natural i=0; i<size_; ++i)
01705   {
01706     MRSASSERT(rhs.data_[i] != 0);
01707     data_[i] /= rhs.data_[i];
01708   }
01709   return *this;
01710 }
01711 
01722 void realvec::matrixMulti(const mrs_realvec& a,const mrs_realvec& b,mrs_realvec& out)
01723 {
01724   //naive Matrix multiplication
01725 
01726   MRSASSERT(a.getCols()==b.getRows());
01727   MRSASSERT(out.getRows()==a.getRows());
01728   MRSASSERT(out.getCols()==b.getCols());
01729 
01730   out.setval (0.);
01731 
01732   for (mrs_natural r=0; r<out.getRows(); ++r)
01733   {
01734     for (mrs_natural c=0; c<out.getCols(); ++c)
01735     {
01736       for (mrs_natural i=0; i<a.getCols(); ++i)
01737       {
01738         out(r,c)+=a(r,i)*b(i,c);
01739       }
01740     }
01741   }
01742 }
01743 
01744 
01745 
01754 mrs_real realvec::getValueFenced(const mrs_natural i) const
01755 {
01756 
01757   if (i < 0 || i >= (mrs_natural)size_) {
01758     // TODO: use a Marsyas branded exception here?
01759     throw std::out_of_range("realvec indexing out of bounds.");
01760   }
01761   return data_[i];
01762 }
01763 
01773 mrs_real realvec::getValueFenced(const mrs_natural r, const mrs_natural c) const
01774 {
01775   if (r < 0 || r >= rows_ || c < 0 || c >= cols_) {
01776     // TODO: use a Marsyas branded exception here?
01777     throw std::out_of_range("realvec indexing out of bounds.");
01778   }
01779   return data_[c * rows_ + r];
01780 }
01781 
01792 mrs_real& realvec::getValueFenced(const mrs_natural i)
01793 {
01794   if (i < 0 || i >= (mrs_natural)size_) {
01795     // TODO: use a Marsyas branded exception here?
01796     throw std::out_of_range("realvec indexing out of bounds.");
01797   }
01798   return data_[i];
01799 }
01800 
01812 mrs_real& realvec::getValueFenced(const mrs_natural r, const mrs_natural c)
01813 {
01814   if (r < 0 || r >= rows_ || c < 0 || c >= cols_) {
01815     // TODO: use a Marsyas branded exception here?
01816     throw std::out_of_range("realvec indexing out of bounds.");
01817   }
01818   return data_[c * rows_ + r];
01819 }
01820 
01821 
01822 }