Marsyas
0.6.0-alpha
|
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 }