psychlops/core/math/psychlops_m_matrix.cpp

   1 /*
   2  *  psychlops_m_matrix.cpp
   3  *  Psychlops Standard Library (Universal)
   4  *
   5  *  Last Modified 2006/03/30 by Kenchi HOSOKAWA
   6  *  (C) 2005 Kenchi HOSOKAWA, Kazushi MARUYA, Takao SATO
   7  */
   8
   9
  10
  11 #include <math.h>
  12 #include <string.h>
  13 #include <stdlib.h>
  14 #include <iostream>
  15 #include <fstream>
  16 #include <vector>
  17
  18 #include "../ApplicationInterfaces/psychlops_code.h"
  19 #include "../devices/psychlops_io_file.h"
  20 #include "../data/psychlops_d_misc.h"
  21 #include "psychlops_m_interval.h"
  22 #include "psychlops_m_function.h"
  23 #include "psychlops_m_matrix.h"
  24
  25 namespace Psychlops {
  26
  27
  28
  29         struct Matrix::elem_ptr {
  30                 friend class Matrix;
  31                 public:
  32                 double *origin, *begin, *row, *col, *rowend, *colend;
  33                 int rowstep, colwidth;
  34                 elem_ptr(const Matrix &mtx);
  35         };
  36         Matrix::elem_ptr::elem_ptr(const Matrix &mtx) {
  37                 mtx.get_element_ptr( origin, begin, rowstep);
  38                 rowend = begin + (rowstep*mtx.rows_);
  39         }
  40
  41
  42
  43         Matrix::Matrix(const Matrix &mtx) : rows_(0), cols_(0), has_instance_(false), element_(NULL), op_(INSTANCE) {
  44                 Matrix tmp;
  45                 if(typeid(mtx) == typeid(MatrixExpression)) {
  46                         tmp = dynamic_cast<const MatrixExpression &>(mtx);
  47                         set_roughly(mtx.rows_, mtx.cols_);
  48                         substitute(tmp);
  49                 } else {
  50                         set_roughly(mtx.rows_, mtx.cols_);
  51                         substitute(mtx);
  52                 }
  53         }
  54         Matrix::Matrix() : rows_(0), cols_(0), has_instance_(false), element_(NULL), op_(INSTANCE) {
  55         }
  56         Matrix::Matrix(int rows, int cols) : rows_(0), cols_(0), has_instance_(false), element_(NULL), op_(INSTANCE) {
  57                 set(rows, cols);
  58         }
  59         Matrix & Matrix::set(int rows, int cols) {
  60                 if(!has_instance_) {
  61                         set_roughly(rows, cols);
  62                         substitute(0.0);
  63                 }
  64                 return *this;
  65         }
  66         Matrix & Matrix::from(double *array)
  67         {
  68                 if(has_instance_) {
  69                         memcpy(element_, array, rows_*cols_*sizeof(double));
  70                 }
  71                 return *this;
  72         }
  73         Matrix & Matrix::from(double *array, int rows, int cols)
  74         {
  75                 if(!has_instance_ || rows*cols != rows_*cols_) {\r
  76                         release();
  77                         set_roughly(rows, cols);
  78                 }
  79                 memcpy(element_, array, rows*cols*sizeof(double));
  80                 return *this;
  81         }
  82         Matrix & Matrix::set(Range row, double interval_rows, Range col, double interval_cols) {
  83                 int rows = (int)( (row.end-row.begin) / interval_rows );
  84                 int cols = (int)( (col.end-col.begin) / interval_cols );
  85                 double rown = row.begin - interval_rows, coln = col.begin;
  86                 if(!has_instance_) {
  87                         set_roughly(rows, cols);
  88                         for(int i=1; i<=rows; i++) {
  89                                 rown += interval_rows;
  90                                 coln = col.begin - interval_cols;
  91                                 for(int j=1; j<=cols; j++) {
  92                                         coln += interval_cols;
  93                                         operator ()(i,j) = rown + coln;
  94                                 }
  95                         }
  96                 }
  97                 return *this;
  98         }
  99         Matrix & Matrix::set_roughly(int rows, int cols) {
 100                 if(rows<=0) throw Exception(typeid(*this), "FORMAT ERROR", "Number of rows must be positive integer.");
 101                 if(cols<=0) throw Exception(typeid(*this), "FORMAT ERROR", "Number of columns must be positive integer.");
 102                 if(!has_instance_) {
 103                         rows_ = rows;
 104                         cols_ = cols;
 105                         elementSize_ = rows_ * cols_;
 106                         element_ = new double[elementSize_];
 107                         has_instance_ = true;
 108                 }
 109                 return *this;
 110         }
 111         Matrix::~Matrix() {
 112                 release();
 113         }
 114         void Matrix::destroy() {
 115                 release();
 116         }
 117         void Matrix::release() {
 118                 if(has_instance_) {
 119                         delete [] element_;
 120                         element_ = NULL;
 121                         rows_ = cols_ = elementSize_ = 0;
 122                         has_instance_ = false;
 123                 }
 124         }
 125         void Matrix::swap(Matrix &rhs) {
 126                 int rows, cols;
 127                 int elementSize;
 128                 double *element;
 129                 Matrix mtx;
 130                 if(has_instance_ && rhs.has_instance_) {
 131                         rows = rows_;
 132                         cols = cols_;
 133                         elementSize = elementSize_;
 134                         element = element_;
 135
 136                         rows_ = rhs.rows_;
 137                         cols_ = rhs.cols_;
 138                         elementSize_ = rhs.elementSize_;
 139                         element_ = rhs.element_;
 140
 141                         rhs.rows_ = rows;
 142                         rhs.cols_ = cols;
 143                         rhs.elementSize_ = elementSize;
 144                         rhs.element_ = element;
 145                 } else if(has_instance() && rhs.has_instance()) {
 146                         mtx.set(rows_, cols_);
 147                         mtx = *this;
 148                         *this = rhs;
 149                         rhs = mtx;
 150                 } else {
 151                         throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 152                 }
 153         }
 154
 155         void Matrix::track(int n) const {
 156                 for(int i=0; i<n; i++) std::cout << " ";
 157                 std::cout << "Matrix " << rows_ << "*" << cols_ << ":" << this << std::endl;
 158         }
 159         bool Matrix::has_instance() const {
 160                 return has_instance_;
 161         }
 162         Matrix* Matrix::instance_ptr() const {
 163                 if(has_instance_) return const_cast<Matrix *>(this);
 164                 else return NULL;
 165         }
 166         void Matrix::get_element_ptr(double *&ptr_origin, double *&ptr_start, int &row_step) const {
 167                 if(has_instance_) {
 168                         ptr_origin = element_;
 169                         ptr_start = element_;
 170                         row_step = cols_;
 171                         //col_width = cols_;
 172                         return;
 173                 }
 174                 else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 175         }
 176         void Matrix::track_and_calculate(Matrix& target, bool add_posi, bool mul_posi) {
 177                 throw Exception(typeid(*this), "SYSTEM EXPRESSION ERROR", "Please notify this error to developers");
 178         }
 179         bool Matrix::track_self(const Matrix* self, int &cnt_in_add, int &cnt_in_mul, bool in_mul) {
 180                 if(self==this) {
 181                         if(in_mul) cnt_in_add++;
 182                         else cnt_in_mul++;
 183                         return true;
 184                 } else {
 185                         return false;
 186                 }
 187         }
 188
 189         Matrix& Matrix::operator =(double rhs) {
 190                 if(has_instance_) substitute(rhs);
 191                 else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 192                 return *this;
 193         }
 194         Matrix& Matrix::operator =(const Matrix &rhs) {
 195                 if(!has_instance_) set_roughly(rhs.rows_, rhs.cols_);
 196                 if(check_size_equal(rhs)) {
 197                         substitute(rhs);
 198                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
 199                 return *this;
 200         }
 201         Matrix& Matrix::operator =(const MatrixExpression &rhs) {
 202                 if(!has_instance_) set_roughly(rhs.rows_, rhs.cols_);
 203                 if(check_size_equal(rhs)) {
 204                         if(rhs.has_instance()) substitute(rhs);
 205                         else const_cast<MatrixExpression &>(rhs).evaluate(*this);
 206                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
 207                 return *this;
 208         }
 209         void Matrix::substitute(double rhs) {
 210                 elem_ptr t(*this);
 211                 for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
 212                         t.colend=t.row+this->cols_;
 213                         for(t.col=t.row; t.col<t.colend; t.col++)
 214                                 *(t.col)=rhs;
 215                 }
 216         }
 217         void Matrix::substitute(const Matrix &rhs) {
 218                 elem_ptr t(*this), r(rhs);
 219                 for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
 220                         t.colend=t.row+this->cols_;
 221                         for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
 222                                 *(t.col)=*(r.col);
 223                 }
 224         }
 225
 226         double& Matrix::xy(int x, int y) const {
 227                 return element_[y*cols_+x];
 228         }
 229         double& Matrix::operator ()(int row, int col) const {
 230                 if(element_==NULL) throw Exception(typeid(*this), "OPERAND ERROR", "This operated matrix has no instance");
 231                 return element_[(row-1)*cols_+(col-1)];
 232 //              elem_ptr elm(*this);
 233 //              if(elm.begin==NULL) throw Exception(typeid(*this), "OPERAND ERROR", "This operated matrix has no instance");
 234 //              return *(elm.begin + ((row-1)*elm.rowstep) + (col-1));
 235         }
 236
 237         bool Matrix::check_size_equal(const Matrix& rhs) const {
 238                 if((rows_==rhs.rows_) && (cols_==rhs.cols_) ) return true;
 239                 else return false;
 240         }
 241         bool Matrix::check_size_mul(const Matrix& rhs) const {
 242                 if(cols_==rhs.rows_) return true;
 243                 else return false;
 244         }
 245
 246
 247         //      add or sub rhs to *this
 248         void Matrix::add(double rhs, bool add_posi) {
 249                 elem_ptr t(*this);
 250                 if(add_posi) {
 251                         for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
 252                                 t.colend=t.row+this->cols_;
 253                                 for(t.col=t.row; t.col<t.colend; t.col++)
 254                                         *(t.col)+=rhs;
 255                         }
 256                 } else {
 257                         for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
 258                                 t.colend=t.row+this->cols_;
 259                                 for(t.col=t.row; t.col<t.colend; t.col++)
 260                                         *(t.col)-=rhs;
 261                         }
 262                 }
 263         }
 264         void Matrix::add(const Matrix &rhs, bool add_posi) {
 265                 elem_ptr t(*this), r(rhs);
 266                 if(add_posi) {
 267                         for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
 268                                 t.colend=t.row+this->cols_;
 269                                 for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
 270                                         *(t.col)+=*(r.col);
 271                         }
 272                 } else {
 273                         for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
 274                                 t.colend=t.row+this->cols_;
 275                                 for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
 276                                         *(t.col)-=*(r.col);
 277                         }
 278                 }
 279         }
 280         //      add lhs multified by rhs to *this
 281         void Matrix::mul(const Matrix &lhs, double rhs, bool add_posi) {
 282                 elem_ptr t(*this), l(lhs);
 283                 if(add_posi) {
 284                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 285                                 t.colend=t.row+this->cols_;
 286                                 for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
 287                                         *(t.col)+=*(l.col)*rhs;
 288                         }
 289                 } else {
 290                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 291                                 t.colend=t.row+this->cols_;
 292                                 for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
 293                                         *(t.col)-=*(l.col)*rhs;
 294                         }
 295                 }
 296         }
 297         void Matrix::mul(const Matrix &lhs, const Matrix &rhs, bool add_posi) {
 298                 elem_ptr t(*this), l(lhs), r(rhs);
 299                 if(add_posi) {
 300                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 301                                 t.colend=t.row+this->cols_;
 302                                 l.colend=l.row+lhs.cols_;
 303                                 for(t.col=t.row,r.col=r.begin; t.col<t.colend; t.col++,r.col++) {
 304                                         for(l.col=l.row,r.row=r.col; l.col<l.colend; l.col++,r.row+=r.rowstep) {
 305                                                 (*t.col)+=(*l.col)*(*r.row);
 306                                         }
 307                                 }
 308                         }
 309                 } else {
 310                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 311                                 t.colend=t.row+this->cols_;
 312                                 l.colend=l.row+lhs.cols_;
 313                                 for(t.col=t.row,r.col=r.begin; t.col<t.colend; t.col++,r.col++) {
 314                                         for(l.col=l.row,r.row=r.col; l.col<l.colend; l.col++,r.row+=r.rowstep) {
 315                                                 (*t.col)-=(*l.col)*(*r.row);
 316                                         }
 317                                 }
 318                         }
 319                 }
 320         }
 321         void Matrix::mul_mask(const Matrix &lhs, const Matrix &rhs, bool add_posi) {
 322                 elem_ptr t(*this), l(lhs), r(rhs);
 323                 if(add_posi) {
 324                         for(t.row=t.begin,l.row=l.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep,r.row+=r.rowstep) {
 325                                 t.colend=t.row+this->cols_;
 326                                 for(t.col=t.row,l.col=l.row,r.col=r.row; t.col<t.colend; t.col++,l.col++,r.col++)
 327                                         *(t.col)+=(*(l.col))*(*(r.col));
 328                         }
 329                 } else {
 330                         for(t.row=t.begin,l.row=l.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep,r.row+=r.rowstep) {
 331                                 t.colend=t.row+this->cols_;
 332                                 for(t.col=t.row,l.col=l.row,r.col=r.row; t.col<t.colend; t.col++,l.col++,r.col++)
 333                                         *(t.col)-=(*(l.col))*(*(r.col));
 334                         }
 335                 }
 336         }
 337         //      add lhs multified by rhs to *this
 338         void Matrix::pow(const Matrix &lhs, double rhs, bool add_posi) {
 339                 int n = (int)floor(rhs);
 340                 if(n != rhs) throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix must powered by int.");
 341                 if(n < 0) throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix must powered by int.");
 342                 n--;
 343                 Matrix tmp(lhs);
 344                 for(int i=0; i<n; i++) tmp = tmp * lhs;
 345         }
 346         void Matrix::pow_mask(const Matrix &lhs, double rhs, bool add_posi) {
 347                 elem_ptr t(*this), l(lhs);
 348                 if(add_posi) {
 349                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 350                                 t.colend=t.row+this->cols_;
 351                                 for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
 352                                         *(t.col)+=::pow(*(l.col), rhs);
 353                         }
 354                 } else {
 355                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 356                                 t.colend=t.row+this->cols_;
 357                                 for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
 358                                         *(t.col)-=::pow(*(l.col), rhs);
 359                         }
 360                 }
 361         }
 362         //      add lhs divided by rhs to *this
 363         void Matrix::div(const Matrix &lhs, double rhs, bool add_posi) {
 364                 elem_ptr t(*this), l(lhs);
 365                 if(add_posi) {
 366                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 367                                 t.colend=t.row+this->cols_;
 368                                 for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
 369                                         *(t.col)+=*(l.col)/rhs;
 370                         }
 371                 } else {
 372                         for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
 373                                 t.colend=t.row+this->cols_;
 374                                 for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
 375                                         *(t.col)-=*(l.col)/rhs;
 376                         }
 377                 }
 378         }
 379
 380         Matrix & Matrix::slide(int drow, int dcol) {
 381                 Matrix after(rows_, cols_);
 382                 elem_ptr t(*this), a(after);
 383                 int d_row = (drow>0) ? drow % rows_ : rows_ - (abs(drow) % rows_),
 384                     d_col = (dcol>0) ? dcol % cols_ : cols_ - (abs(dcol) % cols_);
 385                 size_t size_double = sizeof(double);
 386                 for(int i=0; i<rows_; i++) {
 387                         memcpy(reinterpret_cast<void *>( a.begin+(a.rowstep*((d_row+i)%rows_))+d_col ), reinterpret_cast<void *>( t.begin+t.rowstep*i ), size_double*(cols_-d_col));
 388                         memcpy(reinterpret_cast<void *>( a.begin+(a.rowstep*((d_row+i)%rows_)) ), reinterpret_cast<void *>( t.begin+t.rowstep*i+(cols_-d_col) ), size_double*(d_col));
 389                 }
 390                 swap(after);
 391                 return *this;
 392         }
 393         Matrix & Matrix::transpose() {
 394                 Matrix after(cols_, rows_);
 395                 elem_ptr t(*this), l(after);
 396                 for(t.row=t.begin,l.col=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.col++) {
 397                         t.colend=t.row+this->cols_;
 398                         for(t.col=t.row,l.row=l.col; t.col<t.colend; t.col++,l.row+=l.rowstep)
 399                                 *(l.row)=*(t.col);
 400                 }
 401                 swap(after);
 402                 return *this;
 403         }
 404         Matrix & Matrix::rot90(int times) {
 405                 Matrix after;
 406                 int time = times>=0 ? times%4 : 4-(-times)%4;
 407
 408                 switch(time) {
 409                         case 0:
 410                         case 4:
 411                                 return *this;
 412                         case 1:
 413                         case 3:
 414                                 after.set(cols_, rows_);
 415                                 break;
 416                         case 2:
 417                                 after.set(rows_, cols_);
 418                                 break;
 419                 }
 420                 elem_ptr t(*this), l(after);
 421                 switch(time) {
 422                         case 0:
 423                         case 4:
 424                                 break;
 425                         case 1:
 426                                 for(t.row=t.begin,l.col=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.col++) {
 427                                         t.colend=t.row+this->cols_;
 428                                         for(t.col=t.row,l.row=l.col+l.rowstep*(after.rows_-1); t.col<t.colend; t.col++,l.row-=l.rowstep)
 429                                                 *(l.row)=*(t.col);
 430                                 }
 431                                 break;
 432                         case 3:
 433                                 for(t.row=t.begin,l.col=l.begin+(after.cols_-1); t.row<t.rowend; t.row+=t.rowstep,l.col--) {
 434                                         t.colend=t.row+this->cols_;
 435                                         for(t.col=t.row,l.row=l.col; t.col<t.colend; t.col++,l.row+=l.rowstep)
 436                                                 *(l.row)=*(t.col);
 437                                 }
 438                                 break;
 439                         case 2:
 440                                 for(t.row=t.begin,l.row=l.rowend-l.rowstep; t.row<t.rowend; t.row+=t.rowstep,l.row-=l.rowstep) {
 441                                         t.colend=t.row+this->cols_;
 442                                         for(t.col=t.row,l.col=l.row+(after.cols_-1); t.col<t.colend; t.col++,l.col--)
 443                                                 *(l.col)=*(t.col);
 444                                 }
 445                                 break;
 446                         default:
 447                                 throw Exception(typeid(*this), "ARGUMENT ERROR", "Notify this error to developers.");
 448                 }
 449                 swap(after);
 450                 return *this;
 451         }
 452         Matrix & Matrix::reshape(int rows, int cols) {
 453                 if(rows*cols!=rows*cols) throw Exception(typeid(*this), "EXPRESSION ERROR", "Rows*Cols must be equal before and after the reshape.");
 454                 rows_ = rows;
 455                 cols_ = cols;
 456                 return *this;
 457         }
 458         Matrix & Matrix::catRow(Matrix &rhs) {
 459                 Matrix tmp(rows_+rhs.rows_, cols_);
 460                 Range rng1,rng2;
 461                 if(has_instance_ && rhs.has_instance_ && (cols_==rhs.cols_) ) {
 462                         tmp(1<=rng1<=rows_, 1<=rng2<=cols_) = *this;
 463                         tmp(rows_+1<=rng1<=tmp.rows_, 1<=rng2<=cols_) = rhs;
 464                         swap(tmp);
 465                 } else {
 466                         throw Exception(typeid(*this), "EXPRESSION ERROR", "Both cols must be same when cat Matrix in row.");
 467                 }
 468                 return *this;
 469         }
 470
 471
 472         double Matrix::max() {
 473                 double tmpmax=operator ()(1,1), tmp;
 474                 for(int i=1; i<=rows_; i++) {
 475                         for(int j=1; j<=cols_; j++) {
 476                                 tmp = operator ()(i,j);
 477                                 tmpmax = ( tmpmax<=tmp ? tmp : tmpmax );
 478                         }
 479                 }
 480                 return tmpmax;
 481         }
 482         double Matrix::min() {
 483                 double tmpmin=operator ()(1,1), tmp;
 484                 for(int i=1; i<=rows_; i++) {
 485                         for(int j=1; j<=cols_; j++) {
 486                                 tmp = operator ()(i,j);
 487                                 tmpmin = ( tmpmin>=tmp ? tmp : tmpmin );
 488                         }
 489                 }
 490                 return tmpmin;
 491         }
 492
 493         double Matrix::trace() {
 494                 if(rows_!=cols_) throw Exception(typeid(Matrix), "EXPRESSION ERROR", "Matrix size did not matched");
 495                 double tmptrace = 0.0;
 496                 for(int i=1; i<=rows_; i++) {
 497                         tmptrace += operator ()(i,i);
 498                 }
 499                 return tmptrace;
 500         }
 501
 502
 503 /*      Matrix& Matrix::convolute(const Matrix &rhs) {
 504                 elem_ptr t(*this), r(rhs);
 505                 for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
 506                         t.colend=t.row+this->cols_;
 507                         for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
 508                                 *(t.col) *= *(r.col);
 509                 }
 510                 return *this;
 511         }
 512 */
 513
 514         int Matrix::getRows() const { return rows_; }
 515         int Matrix::getCols() const { return cols_; }
 516
 517
 518         Matrix & Matrix::each(double (*f)(double)) {    // this method only for test. Don't use for your code.
 519                 elem_ptr t(*this);
 520                 for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
 521                         t.colend=t.row+this->cols_;
 522                         for(t.col=t.row; t.col<t.colend; t.col++)
 523                                 *(t.col)=f(*(t.col));
 524                 }
 525                 return *this;
 526         }
 527         Matrix & Matrix::each(double (*f)(double, double), double second) {     // this method only for test. Don't use for your code.
 528                 elem_ptr t(*this);
 529                 for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
 530                         t.colend=t.row+this->cols_;
 531                         for(t.col=t.row; t.col<t.colend; t.col++)
 532                                 *(t.col)=f(*(t.col), second);
 533                 }
 534                 return *this;
 535         }
 536
 537 /*
 538         Matrix Matrix::mesh(const Wave &roww, int rows, const Wave &colw, int cols) {
 539                 Range row, col;
 540                 1<=row<=rows;
 541                 1<=col<=cols;
 542                 return mesh(roww, row, colw, col);
 543         }
 544         Matrix Matrix::mesh(const Wave &roww, const Range rowrng, const Wave &colw, const Range colrng) {
 545                 int rowf = rowrng.int_floor(), rowc = rowrng.int_ceil(), colf = colrng.int_floor(), colc = colrng.int_ceil();
 546                 Matrix mtx(rowc-rowf+1, colc-colf+1);
 547                 elem_ptr p(mtx);
 548                 int row, col;
 549                 for(row=rowf, p.row=p.begin; row<=rowc; row++, p.row+=p.rowstep) {
 550                         p.colend=p.row+mtx.cols_;
 551                         for(col=colf, p.col=p.row; col<=colc; col++, p.col++) {
 552                                 *(p.col) = roww(row)+colw(col);
 553                         }
 554                 }
 555                 return mtx;
 556         }
 557         Matrix Matrix::mesh(const Range rowrng, const Range colrng) {
 558                 return mesh(Waveform::LINEAR(1), rowrng, Waveform::LINEAR(1), colrng);
 559         }
 560 //      MatrixExpression Matrix::mesh(Wave &roww, Wave &colw) {
 561 //              return MatrixExpression(SIZE_UNDEFINED, SIZE_UNDEFINED, roww, colw, MatrixExpression::MESH);
 562 //      }
 563 */
 564         void Matrix::mesh(const Range rowrng, Matrix &rowmtx, const Range colrng, Matrix &colmtx) {
 565                 int rowf = rowrng.int_floor(), rowc = rowrng.int_ceil(), colf = colrng.int_floor(), colc = colrng.int_ceil();
 566                 rowmtx = mesh(rowrng, colc-colf+1);
 567                 colmtx = mesh(rowc-rowf+1, colrng);
 568         }
 569         Matrix Matrix::mesh(int rows, const Range colrng) {
 570                 int colf = colrng.int_floor(), colc = colrng.int_ceil();
 571                 Matrix mtx(rows, colc-colf+1);
 572                 elem_ptr p(mtx);
 573                 int row, col;
 574                 for(row=0, p.row=p.begin; row<rows; row++, p.row+=p.rowstep) {
 575                         p.colend=p.row+mtx.cols_;
 576                         for(col=colf, p.col=p.row; col<=colc; col++, p.col++) {
 577                                 *(p.col) = col;
 578                         }
 579                 }
 580                 return mtx;
 581         }
 582         Matrix Matrix::mesh(const Range rowrng, int cols) {
 583                 int rowf = rowrng.int_floor(), rowc = rowrng.int_ceil();
 584                 Matrix mtx(rowc-rowf+1, cols);
 585                 elem_ptr p(mtx);
 586                 int row, col;
 587                 for(row=rowf, p.row=p.begin; row<=rowc; row++, p.row+=p.rowstep) {
 588                         p.colend=p.row+mtx.cols_;
 589                         for(col=0, p.col=p.row; col<cols; col++, p.col++) {
 590                                 *(p.col) = row;
 591                         }
 592                 }
 593                 return mtx;
 594         }
 595
 596         double* Matrix::getElementPtr()
 597         {
 598                 if(op_==INSTANCE) return element_;
 599                 else return 0;
 600         }
 601
 602         void Matrix::copy(double *target, Interval init_col, int row)
 603         {
 604                 int top = init_col.int_ceil(cols_), bottom = init_col.int_floor(0), num = top-bottom;
 605                 for(int i=0; i<num; i++) {
 606                         target[i] = (*this)(row, i+bottom);
 607                 }
 608         }
 609
 610         void Matrix::load(std::string filename)
 611         {
 612                 std::vector<std::vector<std::string> > csv;
 613                 Data::loadCSVraw( File::decodePath(filename).c_str(), csv );
 614
 615                 int max_row = csv.size();
 616                 int max_col = 0, tmp_col = 0;
 617                 for(int row=0; row<max_row; row++) {
 618                         tmp_col = csv[row].size();
 619                         if(tmp_col > max_col) max_col = tmp_col;
 620                 }
 621
 622                 set(max_row, max_col);
 623                 double num;
 624                 for(int row=0; row<max_row; row++) {
 625                         for(int col=0; col<csv[row].size(); col++) {
 626                                 num = atof( csv[row][col].c_str() );
 627                                 (*this)(row+1, col+1) = num;
 628                         }
 629                 }
 630
 631         }
 632         void Matrix::save(const std::string filename)
 633         {
 634                 std::ofstream of( File::decodePath(filename).c_str());
 635                 for(int row=0; row<getRows(); row++) {
 636                         for(int col=0; col<getCols()-1; col++) {
 637                                 of << (*this)(row+1, col+1) << ", ";
 638                         }
 639                         of << (*this)(row+1, getCols()) << std::endl;
 640                 }
 641         }
 642
 643
 644
 645
 646         char MatrixExpression::OPERAND_SYMBOL[10] = {'N','S','+','-','*','/','+','-','*','/'};
 647
 648         MatrixExpression Matrix::operator ()(int row, Range col) {
 649                 Range rng;
 650                 return MatrixExpression(1, col.int_ceil(cols_)-col.int_floor(1)+1, *this, row<=rng<=row, col, MatrixExpression::SLICE);
 651         }
 652         MatrixExpression Matrix::operator ()(Range row, int col) {
 653                 Range rng;
 654                 return MatrixExpression(row.int_ceil(rows_)-row.int_floor(1)+1, 1, *this, row, col<=rng<=col, MatrixExpression::SLICE);
 655         }
 656         MatrixExpression Matrix::operator ()(Range row, Range col) {
 657                 return MatrixExpression(row.int_ceil(rows_)-row.int_floor(1)+1, col.int_ceil(cols_)-col.int_floor(1)+1, *this, row, col, MatrixExpression::SLICE);
 658         }
 659         MatrixExpression operator +(const Matrix &lhs, double rhs) {
 660                 return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::ADD_SC);
 661         }
 662         MatrixExpression operator -(const Matrix &lhs, double rhs) {
 663                 return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::SUB_SC);
 664         }
 665         MatrixExpression operator *(const Matrix &lhs, double rhs) {
 666                 return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::MUL_SC);
 667         }
 668         MatrixExpression operator /(const Matrix &lhs, double rhs) {
 669                 return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::DIV_SC);
 670         }
 671         MatrixExpression operator ^(const Matrix &lhs, double rhs) {
 672                 if(lhs.op_!=Matrix::MASK) {
 673                         if(lhs.rows_==lhs.cols_) return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::POW);
 674                         else throw Exception("Argument Error ^: This operator is under testing.");
 675                 } else {
 676                         return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::POW_MASK);
 677                 }
 678         }
 679         MatrixExpression operator +(const Matrix &lhs, const Matrix &rhs) {
 680                 if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::ADD);
 681                 else throw Exception("Argument Error +: Matrix size was unmatched.");
 682         }
 683         MatrixExpression operator -(const Matrix &lhs, const Matrix &rhs) {
 684                 if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::SUB);
 685                 else throw Exception("Argument Error -: Matrix size was unmatched.");
 686         }
 687         MatrixExpression operator *(const Matrix &lhs, const Matrix &rhs) {
 688                 if(lhs.op_!=Matrix::MASK && rhs.op_!=Matrix::MASK) {
 689                         if(lhs.check_size_mul(rhs)) return MatrixExpression(lhs.rows_, rhs.cols_, lhs, rhs, MatrixExpression::MUL);
 690                         else throw Exception("Argument Error *: Matrix size was unmatched.");
 691                 } else {
 692                         if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, rhs.cols_, lhs, rhs, MatrixExpression::MUL_MASK);
 693                         else throw Exception("Argument Error *~: Matrix size was unmatched.");
 694                 }
 695         }
 696         MatrixExpression operator ~(const Matrix &lhs) {
 697                 Range rng;
 698                 return MatrixExpression(lhs.rows_, lhs.cols_, lhs, 1<=rng<=lhs.rows_, 1<=rng<=lhs.cols_, MatrixExpression::MASK);
 699         }
 700 //      MatrixExpression operator /(const Matrix &lhs, const Matrix &rhs) {
 701 //              if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, rhs.cols_, lhs, rhs, MatrixExpression::DIV);
 702 //              else throw Exception();
 703 //      }
 704
 705         MatrixExpression::MatrixExpression() : lhs_(NULL), rhs_(NULL) {}
 706         MatrixExpression::MatrixExpression(int rows, int cols, const Matrix &f, const Matrix &f2, OPERAND op)
 707                 : lhs_(const_cast<Matrix*>(&f)), rhs_(const_cast<Matrix*>(&f2)) {
 708                 if(op==ADD || op==SUB || op==MUL || op==DIV || op==MUL_MASK) {
 709                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
 710                 op_ = op;
 711                 rows_ = rows;
 712                 cols_ = cols;
 713                 element_ = NULL;
 714         }
 715         MatrixExpression::MatrixExpression(int rows, int cols, const Matrix &f, double f2, OPERAND op)
 716                 : lhs_(const_cast<Matrix*>(&f)), rhs_(NULL), scholar_(f2) {
 717                 if(op==ADD_SC || op==SUB_SC || op==MUL_SC || op==DIV_SC || op==POW || op==POW_MASK) {
 718                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
 719                 op_ = op;
 720                 rows_ = rows;
 721                 cols_ = cols;
 722                 element_ = NULL;
 723         }
 724         MatrixExpression::MatrixExpression(int rows, int cols, const Matrix& f, Range& rowrng, Range& colrng, OPERAND op)
 725                 : lhs_(const_cast<Matrix*>(&f)), rhs_(NULL), rowrng_(rowrng), colrng_(colrng) {
 726                 if(op==SLICE || op==MASK) {
 727                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
 728                 op_ = op;
 729                 rows_ = rows;
 730                 cols_ = cols;
 731                 element_ = NULL;
 732         }
 733         MatrixExpression::MatrixExpression(int rows, int cols, Wave &roww, Wave &colw, OPERAND op)
 734                 : lhs_(NULL), rhs_(NULL), meshrow_(&roww), meshcol_(&colw) {
 735                 if(op==MESH) {
 736                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
 737                 op_ = op;
 738                 rows_ = rows;
 739                 cols_ = cols;
 740                 element_ = NULL;
 741         }
 742         MatrixExpression::MatrixExpression(const MatrixExpression &mtxc) : lhs_(mtxc.rhs_), rhs_(mtxc.rhs_) {
 743                 op_ = mtxc.op_;
 744         }
 745
 746
 747         void MatrixExpression::swap(Matrix &rhs) {
 748                 Matrix mtx;
 749                 if(has_instance() && rhs.has_instance()) {
 750                         mtx.set(rows_, cols_);
 751                         mtx = *this;
 752                         *this = rhs;
 753                         rhs = mtx;
 754                 }
 755                 else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 756         }
 757         void MatrixExpression::track(int n) const {
 758                 for(int i=0; i<n; i++) std::cout << " ";
 759                 std::cout << "Matrix op" << MatrixExpression::OPERAND_SYMBOL[op_] << ":" << this << std::endl;
 760                 if(lhs_!=(Matrix*)this) lhs_->track(++n);
 761                 if(rhs_!=(Matrix*)this) rhs_->track(n);
 762         }
 763         MatrixExpression::~MatrixExpression() {
 764                 //std::cout << "DEST-MatrixExpression op" << MatrixExpression::OPERAND_SYMBOL[op_] << " " << rows_ << "*" << cols_ << ":" << this << std::endl;
 765         }
 766
 767
 768         bool MatrixExpression::has_instance() const {
 769                 if(op_==SLICE || op_==MASK) return lhs_->has_instance();
 770                 else return false;
 771         }
 772         Matrix* MatrixExpression::instance_ptr() const {
 773                 if(op_==SLICE || op_==MASK) return lhs_->instance_ptr();
 774                 else return NULL;
 775         }
 776         void MatrixExpression::get_element_ptr(double *&ptr_origin, double *&ptr_start, int &row_step) const {
 777                 if(op_==SLICE || op_==MASK) {
 778                         lhs_->get_element_ptr(ptr_origin, ptr_start, row_step);
 779                         ptr_start = ptr_origin+(colrng_.int_floor(1)-1)+((rowrng_.int_floor(1)-1)*row_step);
 780                         //col_width = cols_;
 781                         return;
 782                 } else throw Exception(typeid(*this), "MEMORY ERROR", "No instance.");
 783         }
 784
 785         //void MatrixExpression::substitute(double rhs) {
 786         //      Matrix::substitute(rhs);
 787         //}
 788         //void MatrixExpression::substitute(const Matrix &mtx) {
 789         //      const_cast<Matrix &>(mtx).track_and_calculate(*this, true, true);
 790         //}
 791
 792         void MatrixExpression::evaluate(Matrix& target) {
 793                 Matrix *newself = NULL, *actualtarget = target.instance_ptr();
 794                 int cnt_in_add=0, cnt_in_mul=0;
 795                 if( track_self(actualtarget, cnt_in_add, cnt_in_mul, false) ) {
 796                         newself = new Matrix(actualtarget->rows_, actualtarget->cols_);
 797                         track_and_calculate(*newself, true, true);
 798                         actualtarget->swap(*newself);
 799                         delete newself;
 800                 } else {
 801                         target.substitute(0.0);
 802                         track_and_calculate(target, true, true);
 803                 }
 804         }
 805         void MatrixExpression::evaluate(MatrixExpression& target) {
 806                 Matrix *newself = NULL, *actualtarget = target.instance_ptr();
 807                 int cnt_in_add=0, cnt_in_mul=0;
 808                 if( track_self(actualtarget, cnt_in_add, cnt_in_mul, false) ) {
 809                         newself = new Matrix(*actualtarget);
 810                         actualtarget = target.lhs_;
 811                         target.lhs_ = newself;
 812                         target = 0;
 813                         track_and_calculate(target, true, true);
 814                         actualtarget->swap(*newself);
 815                         delete newself;
 816                 } else {
 817                         target.substitute(0.0);
 818                         track_and_calculate(target, true, true);
 819                 }
 820         }
 821         bool MatrixExpression::track_self(const Matrix* self, int &cnt_in_add, int &cnt_in_mul, bool in_mul) {
 822                 bool found = false, in_mult = in_mul;
 823                 switch(op_) {
 824 //                      case ADD:
 825 //                      case SUB:
 826 //                      case ADD_SC:
 827 //                      case SUB_SC:
 828 //                              break;
 829                         case MUL:
 830                         case DIV:
 831                         case MUL_SC:
 832                                 in_mult = true;
 833                                 break;
 834                 }
 835                 if(lhs_!=NULL) found = found || lhs_->track_self(self,cnt_in_add,cnt_in_mul,in_mult);
 836                 if(rhs_!=NULL) found = found || rhs_->track_self(self,cnt_in_add,cnt_in_mul,in_mult);
 837                 return found;
 838         }
 839         void MatrixExpression::track_and_calculate(Matrix& target, bool add_posi, bool mul_posi) {
 840                 Matrix *buf1, *buf2;
 841                 bool made_buf1=false, made_buf2=false;
 842                 if(op_==ADD) {
 843                         if(lhs_->has_instance()) target.add(*lhs_, add_posi);
 844                         else lhs_->track_and_calculate(target, add_posi, mul_posi);
 845                         if(rhs_->has_instance()) target.add(*rhs_, add_posi);
 846                         else rhs_->track_and_calculate(target, add_posi, mul_posi);
 847                 } else if(op_==SUB) {
 848                         if(lhs_->has_instance()) target.add(*lhs_, add_posi);
 849                         else lhs_->track_and_calculate(target, add_posi, mul_posi);
 850                         if(rhs_->has_instance()) target.add(*rhs_, !add_posi);
 851                         else rhs_->track_and_calculate(target, !add_posi, mul_posi);
 852                 } else if(op_==MUL) {
 853                         if(lhs_->has_instance()) buf1=lhs_;
 854                         else {
 855                                 buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
 856                                 made_buf1 = true;
 857                                 lhs_->track_and_calculate(*buf1, true, true);
 858                         }
 859                         if(rhs_->has_instance()) buf2=rhs_;
 860                         else {
 861                                 buf2 = new Matrix(rhs_->rows_, rhs_->cols_);
 862                                 made_buf2 = true;
 863                                 rhs_->track_and_calculate(*buf2, true, true);
 864                         }
 865                         target.mul(*buf1, *buf2, add_posi);
 866                         if(made_buf1) delete buf1;
 867                         if(made_buf2) delete buf2;
 868                 } else if(op_==ADD_SC) {
 869                         if(lhs_->has_instance()) target.add(*lhs_, add_posi);
 870                         else lhs_->track_and_calculate(target, add_posi, mul_posi);
 871                         target.add(scholar_, add_posi);
 872                 } else if(op_==SUB_SC) {
 873                         if(lhs_->has_instance()) target.add(*lhs_, add_posi);
 874                         else lhs_->track_and_calculate(target, add_posi, mul_posi);
 875                         target.add(scholar_, !add_posi);
 876                 } else if(op_==MUL_SC) {
 877                         if(lhs_->has_instance()) buf1=lhs_;
 878                         else {
 879                                 buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
 880                                 made_buf1 = true;
 881                                 lhs_->track_and_calculate(*buf1, true, true);
 882                         }
 883                         target.mul(*buf1, scholar_, add_posi);
 884                         if(made_buf1) delete buf1;
 885                 } else if(op_==DIV_SC) {
 886                         if(lhs_->has_instance()) buf1=lhs_;
 887                         else {
 888                                 buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
 889                                 made_buf1 = true;
 890                                 lhs_->track_and_calculate(*buf1, true, true);
 891                         }
 892                         target.div(*buf1, scholar_, add_posi);
 893                         if(made_buf1) delete buf1;
 894                 } else if(op_==MUL_MASK) {
 895                         if(lhs_->has_instance()) buf1=lhs_;
 896                         else {
 897                                 buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
 898                                 made_buf1 = true;
 899                                 lhs_->track_and_calculate(*buf1, true, true);
 900                         }
 901                         if(rhs_->has_instance()) buf2=rhs_;
 902                         else {
 903                                 buf2 = new Matrix(rhs_->rows_, rhs_->cols_);
 904                                 made_buf2 = true;
 905                                 rhs_->track_and_calculate(*buf2, true, true);
 906                         }
 907                         target.mul_mask(*buf1, *buf2, add_posi);
 908                         if(made_buf1) delete buf1;
 909                         if(made_buf2) delete buf2;
 910                 } else if(op_==SLICE) {
 911                         throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 912 //                      if(lhs_->has_instance()) buf1=lhs_;
 913 //                      else {
 914 //                              buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
 915 //                              made_buf1 = true;
 916 //                              lhs_->track_and_calculate(*buf1, true, true);
 917 //                      }
 918 //                      target.div(*buf1, scholar_, add_posi);
 919 //                      if(made_buf1) delete buf1;
 920                 }
 921         }
 922
 923         double& MatrixExpression::xy(int x, int y) const {
 924                 elem_ptr elm(*this);
 925                 if(op_==SLICE) {
 926                         if(elm.begin==NULL) throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 927                         return *(elm.begin + (y*elm.rowstep) + x);
 928                 } else {
 929                         throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 930                 }
 931         }
 932         double& MatrixExpression::operator ()(int row, int col) const {
 933                 elem_ptr elm(*this);
 934                 if(op_==SLICE) {
 935                         if(elm.begin==NULL) throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 936                         return *(elm.begin + ((row-1)*elm.rowstep) + (col-1));
 937                 } else {
 938                         throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 939                 }
 940         }
 941         MatrixExpression MatrixExpression::operator ()(int row, Range col) { return (*this).Matrix::operator ()(row, col); }
 942         MatrixExpression MatrixExpression::operator ()(Range row, int col) { return (*this).Matrix::operator ()(row, col); }
 943         MatrixExpression MatrixExpression::operator ()(Range row, Range col) { return (*this).Matrix::operator ()(row, col); }
 944
 945         Matrix& MatrixExpression::operator =(double rhs) {
 946                 if(has_instance()) substitute(rhs);
 947                 else throw Exception();
 948                 return *this;
 949         }
 950         Matrix& MatrixExpression::operator =(const Matrix &rhs) {
 951                 if(check_size_equal(rhs)) {
 952                         if(has_instance()) {
 953                                 substitute(rhs);
 954                         } else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 955                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
 956                 return *this;
 957         }
 958         Matrix& MatrixExpression::operator =(const MatrixExpression &rhs) {
 959                 if(check_size_equal(rhs)) {
 960                         if(has_instance()) {
 961                                 if(rhs.has_instance()) substitute(rhs);
 962                                 else const_cast<MatrixExpression &>(rhs).evaluate(*this);
 963                         } else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 964                 } else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
 965                 return *this;
 966         }
 967
 968
 969         Matrix & MatrixExpression::slide(int drow, int dcol) {
 970                 throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 971                 return *this;
 972         }
 973         Matrix & MatrixExpression::transpose() {
 974                 throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
 975                 return *this;
 976         }
 977
 978         std::ostream& operator<<(std::ostream& ost, const Matrix &mtx) {
 979                 if(mtx.has_instance()) {
 980                         for(int row=1; row<=mtx.rows_; row++) {
 981                                 for(int col=1; col<=mtx.cols_; col++) {
 982                                         ost << mtx(row,col) << "\t";
 983                                 }
 984                                 ost << std::endl;
 985                         }
 986                 }
 987                 return ost;
 988         }
 989
 990
 991 namespace MatrixOverload {
 992 Matrix sin(const Matrix &m) {
 993         Matrix mtx(m);
 994         return mtx.each(&::sin);
 995 }
 996 Matrix cos(const Matrix &m) {
 997         Matrix mtx(m);
 998         return mtx.each(&::cos);
 999 }
1000 Matrix tan(const Matrix &m) {
1001         Matrix mtx(m);
1002         return mtx.each(&::tan);
1003 }
1004 Matrix asin(const Matrix &m) {
1005         Matrix mtx(m);
1006         return mtx.each(&::asin);
1007 }
1008 Matrix acos(const Matrix &m) {
1009         Matrix mtx(m);
1010         return mtx.each(&::acos);
1011 }
1012 Matrix atan(const Matrix &m) {
1013         Matrix mtx(m);
1014         return mtx.each(&::atan);
1015 }
1016 Matrix exp(const Matrix &m) {
1017         Matrix mtx(m);
1018         return mtx.each(&::exp);
1019 }
1020 Matrix log(const Matrix &m) {
1021         Matrix mtx(m);
1022         return mtx.each(&::log);
1023 }
1024 Matrix log10(const Matrix &m) {
1025         Matrix mtx(m);
1026         return mtx.each(&::log10);
1027 }
1028 Matrix sqrt(const Matrix &m) {
1029         Matrix mtx(m);
1030         return mtx.each(&::sqrt);
1031 }
1032
1033 Matrix abs(const Matrix &m) {
1034         Matrix mtx(m);
1035         return mtx.each(&::fabs);
1036 }
1037 Matrix pow(const Matrix &m, double ex) {
1038         Matrix mtx(m);
1039         return mtx.each(&::pow, ex);
1040 }
1041
1042
1043 Matrix max(const Matrix& a, const Matrix b) {
1044         if(!a.check_size_equal(b)) throw Exception(typeid(Matrix), "EXPRESSION ERROR", "Matrix size did not matched");
1045         Matrix tmp(a.getRows(), a.getCols());
1046         for(int col=0; col<a.getCols(); col++) {
1047                 for(int row=0; row<a.getRows(); row++) {
1048                         tmp(row, col) = ( a(row, col) > b(row, col) ? a(row, col) : b(row, col) );
1049                 }
1050         }
1051         return tmp;
1052 }
1053
1054 Matrix min(const Matrix& a, const Matrix b) {
1055         if(!a.check_size_equal(b)) throw Exception(typeid(Matrix), "EXPRESSION ERROR", "Matrix size did not matched");
1056         Matrix tmp(a.getRows(), a.getCols());
1057         for(int col=0; col<a.getCols(); col++) {
1058                 for(int row=0; row<a.getRows(); row++) {
1059                         tmp(row, col) = ( a(row, col) < b(row, col) ? a(row, col) : b(row, col) );
1060                 }
1061         }
1062         return tmp;
1063 }
1064 }
1065
1066 }       /*      <- namespace Psycholops         */
1067
1068