/*
 *  psychlops_m_matrix.cpp
 *  Psychlops Standard Library (Universal)
 *
 *  Last Modified 2006/03/30 by Kenchi HOSOKAWA
 *  (C) 2005 Kenchi HOSOKAWA, Kazushi MARUYA, Takao SATO
 */




#include <iostream>

#include "math.h"
#include "string.h"
#include "../ApplicationInterfaces/psychlops_code.h"
#include "psychlops_m_interval.h"
#include "psychlops_m_function.h"
#include "psychlops_m_matrix.h"

namespace Psychlops {



	struct Matrix::elem_ptr {
		friend class Matrix;
		public:
		double *origin, *begin, *row, *col, *rowend, *colend;
		int rowstep, colwidth;
		elem_ptr(const Matrix &mtx);
	};
	Matrix::elem_ptr::elem_ptr(const Matrix &mtx) {
		mtx.get_element_ptr( origin, begin, rowstep);
		rowend = begin + (rowstep*mtx.rows_);
	}



	Matrix::Matrix(const Matrix &mtx) : rows_(0), cols_(0), has_instance_(false), element_(NULL), op_(INSTANCE) {
		Matrix tmp;
		if(typeid(mtx) == typeid(MatrixExpression)) {
			tmp = dynamic_cast<const MatrixExpression &>(mtx);
			set_roughly(mtx.rows_, mtx.cols_);
			substitute(tmp);
		} else {
			set_roughly(mtx.rows_, mtx.cols_);
			substitute(mtx);
		}
	}
	Matrix::Matrix() : rows_(0), cols_(0), has_instance_(false), element_(NULL), op_(INSTANCE) {
	}
	Matrix::Matrix(int rows, int cols) : rows_(0), cols_(0), has_instance_(false), element_(NULL), op_(INSTANCE) {
		set(rows, cols);
	}
	Matrix & Matrix::set(int rows, int cols) {
		if(!has_instance_) {
			set_roughly(rows, cols);
			substitute(0.0);
		}
		return *this;
	}
	Matrix & Matrix::set(Range row, double interval_rows, Range col, double interval_cols) {
		int rows = (int)( (row.end-row.begin) / interval_rows );
		int cols = (int)( (col.end-col.begin) / interval_cols );
		double rown = row.begin - interval_rows, coln = col.begin;
		if(!has_instance_) {
			set_roughly(rows, cols);
			for(int i=1; i<=rows; i++) {
				rown += interval_rows;
				coln = col.begin - interval_cols;
				for(int j=1; j<=cols; j++) {
					coln += interval_cols;
					operator ()(i,j) = rown + coln;
				}
			}
		}
		return *this;
	}
	Matrix & Matrix::set_roughly(int rows, int cols) {
		if(rows<=0) throw Exception(typeid(*this), "FORMAT ERROR", "Number of rows must be positive integer.");
		if(cols<=0) throw Exception(typeid(*this), "FORMAT ERROR", "Number of columns must be positive integer.");
		if(!has_instance_) {
			rows_ = rows;
			cols_ = cols;
			elementSize_ = rows_ * cols_;
			element_ = new double[elementSize_];
			has_instance_ = true;
		}
		return *this;
	}
	Matrix::~Matrix() {
		release();
	}
	void Matrix::destroy() {
		release();
	}
	void Matrix::release() {
		if(has_instance_) {
			delete [] element_;
			element_ = NULL;
			rows_ = cols_ = elementSize_ = 0;
			has_instance_ = false;
		}
	}
	void Matrix::swap(Matrix &rhs) {
		int rows, cols;
		int elementSize;
		double *element;
		Matrix mtx;
		if(has_instance_ && rhs.has_instance_) {
			rows = rows_;
			cols = cols_;
			elementSize = elementSize_;
			element = element_;

			rows_ = rhs.rows_;
			cols_ = rhs.cols_;
			elementSize_ = rhs.elementSize_;
			element_ = rhs.element_;

			rhs.rows_ = rows;
			rhs.cols_ = cols;
			rhs.elementSize_ = elementSize;
			rhs.element_ = element;
		} else if(has_instance() && rhs.has_instance()) {
			mtx.set(rows_, cols_);
			mtx = *this;
			*this = rhs;
			rhs = mtx;
		} else {
			throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
		}
	}

	void Matrix::track(int n) const {
		for(int i=0; i<n; i++) std::cout << " ";
		std::cout << "Matrix " << rows_ << "*" << cols_ << ":" << this << std::endl;
	}
	bool Matrix::has_instance() const {
		return has_instance_;
	}
	Matrix* Matrix::instance_ptr() const {
		if(has_instance_) return const_cast<Matrix *>(this);
		else return NULL;
	}
	void Matrix::get_element_ptr(double *&ptr_origin, double *&ptr_start, int &row_step) const {
		if(has_instance_) {
			ptr_origin = element_;
			ptr_start = element_;
			row_step = cols_;
			//col_width = cols_;
			return;
		}
		else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
	}
	void Matrix::track_and_calculate(Matrix& target, bool add_posi, bool mul_posi) {
		throw Exception(typeid(*this), "SYSTEM EXPRESSION ERROR", "Please notify this error to developers");
	}
	bool Matrix::track_self(const Matrix* self, int &cnt_in_add, int &cnt_in_mul, bool in_mul) {
		if(self==this) {
			if(in_mul) cnt_in_add++;
			else cnt_in_mul++;
			return true;
		} else {
			return false;
		}
	}

	Matrix& Matrix::operator =(double rhs) {
		if(has_instance_) substitute(rhs);
		else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
		return *this;
	}
	Matrix& Matrix::operator =(const Matrix &rhs) {
		if(!has_instance_) set_roughly(rhs.rows_, rhs.cols_);
		if(check_size_equal(rhs)) {
			substitute(rhs);
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
		return *this;
	}
	Matrix& Matrix::operator =(const MatrixExpression &rhs) {
		if(!has_instance_) set_roughly(rhs.rows_, rhs.cols_);
		if(check_size_equal(rhs)) {
			if(rhs.has_instance()) substitute(rhs);
			else const_cast<MatrixExpression &>(rhs).evaluate(*this);
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
		return *this;
	}
	void Matrix::substitute(double rhs) {
		elem_ptr t(*this);
		for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
			t.colend=t.row+this->cols_;
			for(t.col=t.row; t.col<t.colend; t.col++)
				*(t.col)=rhs;
		}
	}
	void Matrix::substitute(const Matrix &rhs) {
		elem_ptr t(*this), r(rhs);
		for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
			t.colend=t.row+this->cols_;
			for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
				*(t.col)=*(r.col);
		}
	}

	double& Matrix::operator ()(int row, int col) const {
		if(element_==NULL) throw Exception(typeid(*this), "OPERAND ERROR", "This operated matrix has no instance");
		return element_[(row-1)*cols_+(col-1)];
//		elem_ptr elm(*this);
//		if(elm.begin==NULL) throw Exception(typeid(*this), "OPERAND ERROR", "This operated matrix has no instance");
//		return *(elm.begin + ((row-1)*elm.rowstep) + (col-1));
	}

	bool Matrix::check_size_equal(const Matrix& rhs) const {
		if((rows_==rhs.rows_) && (cols_==rhs.cols_) ) return true;
		else return false;
	}
	bool Matrix::check_size_mul(const Matrix& rhs) const {
		if(cols_==rhs.rows_) return true;
		else return false;
	}


	//	add or sub rhs to *this
	void Matrix::add(double rhs, bool add_posi) {
		elem_ptr t(*this);
		if(add_posi) {
			for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row; t.col<t.colend; t.col++)
					*(t.col)+=rhs;
			}
		} else {
			for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row; t.col<t.colend; t.col++)
					*(t.col)-=rhs;
			}
		}
	}
	void Matrix::add(const Matrix &rhs, bool add_posi) {
		elem_ptr t(*this), r(rhs);
		if(add_posi) {
			for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
					*(t.col)+=*(r.col);
			}
		} else {
			for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
					*(t.col)-=*(r.col);
			}
		}
	}
	//	add lhs multified by rhs to *this
	void Matrix::mul(const Matrix &lhs, double rhs, bool add_posi) {
		elem_ptr t(*this), l(lhs);
		if(add_posi) {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
					*(t.col)+=*(l.col)*rhs;
			}
		} else {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
					*(t.col)-=*(l.col)*rhs;
			}
		}
	}
	void Matrix::mul(const Matrix &lhs, const Matrix &rhs, bool add_posi) {
		elem_ptr t(*this), l(lhs), r(rhs);
		if(add_posi) {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				l.colend=l.row+lhs.cols_;
				for(t.col=t.row,r.col=r.begin; t.col<t.colend; t.col++,r.col++) {
					for(l.col=l.row,r.row=r.col; l.col<l.colend; l.col++,r.row+=r.rowstep) {
						(*t.col)+=(*l.col)*(*r.row);
					}
				}
			}
		} else {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				l.colend=l.row+lhs.cols_;
				for(t.col=t.row,r.col=r.begin; t.col<t.colend; t.col++,r.col++) {
					for(l.col=l.row,r.row=r.col; l.col<l.colend; l.col++,r.row+=r.rowstep) {
						(*t.col)-=(*l.col)*(*r.row);
					}
				}
			}
		}
	}
	void Matrix::mul_mask(const Matrix &lhs, const Matrix &rhs, bool add_posi) {
		elem_ptr t(*this), l(lhs), r(rhs);
		if(add_posi) {
			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) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row,r.col=r.row; t.col<t.colend; t.col++,l.col++,r.col++)
					*(t.col)+=(*(l.col))*(*(r.col));
			}
		} else {
			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) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row,r.col=r.row; t.col<t.colend; t.col++,l.col++,r.col++)
					*(t.col)-=(*(l.col))*(*(r.col));
			}
		}
	}
	//	add lhs multified by rhs to *this
	void Matrix::pow(const Matrix &lhs, double rhs, bool add_posi) {
		int n = (int)floor(rhs);
		if(n != rhs) throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix must powered by int.");
		if(n < 0) throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix must powered by int.");
		n--;
		Matrix tmp(lhs);
		for(int i=0; i<n; i++) tmp = tmp * lhs;
	}
	void Matrix::pow_mask(const Matrix &lhs, double rhs, bool add_posi) {
		elem_ptr t(*this), l(lhs);
		if(add_posi) {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
					*(t.col)+=::pow(*(l.col), rhs);
			}
		} else {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
					*(t.col)-=::pow(*(l.col), rhs);
			}
		}
	}
	//	add lhs divided by rhs to *this
	void Matrix::div(const Matrix &lhs, double rhs, bool add_posi) {
		elem_ptr t(*this), l(lhs);
		if(add_posi) {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
					*(t.col)+=*(l.col)/rhs;
			}
		} else {
			for(t.row=t.begin,l.row=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.row+=l.rowstep) {
				t.colend=t.row+this->cols_;
				for(t.col=t.row,l.col=l.row; t.col<t.colend; t.col++,l.col++)
					*(t.col)-=*(l.col)/rhs;
			}
		}
	}

	Matrix & Matrix::slide(int drow, int dcol) {
		Matrix after(rows_, cols_);
		elem_ptr t(*this), a(after);
		int d_row = (drow>0) ? drow % rows_ : rows_ - (abs(drow) % rows_),
		    d_col = (dcol>0) ? dcol % cols_ : cols_ - (abs(dcol) % cols_);
		size_t size_double = sizeof(double);
		for(int i=0; i<rows_; i++) {
			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));
			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));
		}
		swap(after);
		return *this;
	}
	Matrix & Matrix::transpose() {
		Matrix after(cols_, rows_);
		elem_ptr t(*this), l(after);
		for(t.row=t.begin,l.col=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.col++) {
			t.colend=t.row+this->cols_;
			for(t.col=t.row,l.row=l.col; t.col<t.colend; t.col++,l.row+=l.rowstep)
				*(l.row)=*(t.col);
		}
		swap(after);
		return *this;
	}
	Matrix & Matrix::rot90(int times) {
		Matrix after;
		int time = times>=0 ? times%4 : 4-(-times)%4;

		switch(time) {
			case 0:
			case 4:
				return *this;
			case 1:
			case 3:
				after.set(cols_, rows_);
				break;
			case 2:
				after.set(rows_, cols_);
				break;
		}
		elem_ptr t(*this), l(after);
		switch(time) {
			case 0:
			case 4:
				break;
			case 1:
				for(t.row=t.begin,l.col=l.begin; t.row<t.rowend; t.row+=t.rowstep,l.col++) {
					t.colend=t.row+this->cols_;
					for(t.col=t.row,l.row=l.col+l.rowstep*(after.rows_-1); t.col<t.colend; t.col++,l.row-=l.rowstep)
						*(l.row)=*(t.col);
				}
				break;
			case 3:
				for(t.row=t.begin,l.col=l.begin+(after.cols_-1); t.row<t.rowend; t.row+=t.rowstep,l.col--) {
					t.colend=t.row+this->cols_;
					for(t.col=t.row,l.row=l.col; t.col<t.colend; t.col++,l.row+=l.rowstep)
						*(l.row)=*(t.col);
				}
				break;
			case 2:
				for(t.row=t.begin,l.row=l.rowend-l.rowstep; t.row<t.rowend; t.row+=t.rowstep,l.row-=l.rowstep) {
					t.colend=t.row+this->cols_;
					for(t.col=t.row,l.col=l.row+(after.cols_-1); t.col<t.colend; t.col++,l.col--)
						*(l.col)=*(t.col);
				}
				break;
			default:
				throw Exception(typeid(*this), "ARGUMENT ERROR", "Notify this error to developers.");
		}
		swap(after);
		return *this;
	}
	Matrix & Matrix::reshape(int rows, int cols) {
		if(rows*cols!=rows*cols) throw Exception(typeid(*this), "EXPRESSION ERROR", "Rows*Cols must be equal before and after the reshape.");
		rows_ = rows;
		cols_ = cols;
		return *this;
	}
	Matrix & Matrix::catRow(Matrix &rhs) {
		Matrix tmp(rows_+rhs.rows_, cols_);
		Range rng1,rng2;
		if(has_instance_ && rhs.has_instance_ && (cols_==rhs.cols_) ) {
			tmp(1<=rng1<=rows_, 1<=rng2<=cols_) = *this;
			tmp(rows_+1<=rng1<=tmp.rows_, 1<=rng2<=cols_) = rhs;
			swap(tmp);
		} else {
			throw Exception(typeid(*this), "EXPRESSION ERROR", "Both cols must be same when cat Matrix in row.");
		}
		return *this;
	}


	double Matrix::max() {
		double tmpmax=operator ()(1,1), tmp;
		for(int i=1; i<=rows_; i++) {
			for(int j=1; j<=cols_; j++) {
				tmp = operator ()(i,j);
				tmpmax = ( tmpmax<=tmp ? tmp : tmpmax );
			}
		}
		return tmpmax;
	}
	double Matrix::min() {
		double tmpmin=operator ()(1,1), tmp;
		for(int i=1; i<=rows_; i++) {
			for(int j=1; j<=cols_; j++) {
				tmp = operator ()(i,j);
				tmpmin = ( tmpmin>=tmp ? tmp : tmpmin );
			}
		}
		return tmpmin;
	}

	double Matrix::trace() {
		if(rows_!=cols_) throw Exception(typeid(Matrix), "EXPRESSION ERROR", "Matrix size did not matched");
		double tmptrace = 0.0;
		for(int i=1; i<=rows_; i++) {
			tmptrace += operator ()(i,i);
		}
		return tmptrace;
	}


/*	Matrix& Matrix::convolute(const Matrix &rhs) {
		elem_ptr t(*this), r(rhs);
		for(t.row=t.begin,r.row=r.begin; t.row<t.rowend; t.row+=t.rowstep,r.row+=r.rowstep) {
			t.colend=t.row+this->cols_;
			for(t.col=t.row,r.col=r.row; t.col<t.colend; t.col++,r.col++)
				*(t.col) *= *(r.col);
		}
		return *this;
	}
*/

	int Matrix::getRows() const { return rows_; }
	int Matrix::getCols() const { return cols_; }


	Matrix & Matrix::each(double (*f)(double)) {	// this method only for test. Don't use for your code.
		elem_ptr t(*this);
		for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
			t.colend=t.row+this->cols_;
			for(t.col=t.row; t.col<t.colend; t.col++)
				*(t.col)=f(*(t.col));
		}
		return *this;
	}
	Matrix & Matrix::each(double (*f)(double, double), double second) {	// this method only for test. Don't use for your code.
		elem_ptr t(*this);
		for(t.row=t.begin; t.row<t.rowend; t.row+=t.rowstep) {
			t.colend=t.row+this->cols_;
			for(t.col=t.row; t.col<t.colend; t.col++)
				*(t.col)=f(*(t.col), second);
		}
		return *this;
	}

/*
	Matrix Matrix::mesh(const Wave &roww, int rows, const Wave &colw, int cols) {
		Range row, col;
		1<=row<=rows;
		1<=col<=cols;
		return mesh(roww, row, colw, col);
	}
	Matrix Matrix::mesh(const Wave &roww, const Range rowrng, const Wave &colw, const Range colrng) {
		int rowf = rowrng.int_floor(), rowc = rowrng.int_ceil(), colf = colrng.int_floor(), colc = colrng.int_ceil();
		Matrix mtx(rowc-rowf+1, colc-colf+1);
		elem_ptr p(mtx);
		int row, col;
		for(row=rowf, p.row=p.begin; row<=rowc; row++, p.row+=p.rowstep) {
			p.colend=p.row+mtx.cols_;
			for(col=colf, p.col=p.row; col<=colc; col++, p.col++) {
				*(p.col) = roww(row)+colw(col);
			}
		}
		return mtx;
	}
	Matrix Matrix::mesh(const Range rowrng, const Range colrng) {
		return mesh(Waveform::LINEAR(1), rowrng, Waveform::LINEAR(1), colrng);
	}
//	MatrixExpression Matrix::mesh(Wave &roww, Wave &colw) {
//		return MatrixExpression(SIZE_UNDEFINED, SIZE_UNDEFINED, roww, colw, MatrixExpression::MESH);
//	}
*/
	void Matrix::mesh(const Range rowrng, Matrix &rowmtx, const Range colrng, Matrix &colmtx) {
		int rowf = rowrng.int_floor(), rowc = rowrng.int_ceil(), colf = colrng.int_floor(), colc = colrng.int_ceil();
		rowmtx = mesh(rowrng, colc-colf+1);
		colmtx = mesh(rowc-rowf+1, colrng);
	}
	Matrix Matrix::mesh(int rows, const Range colrng) {
		int colf = colrng.int_floor(), colc = colrng.int_ceil();
		Matrix mtx(rows, colc-colf+1);
		elem_ptr p(mtx);
		int row, col;
		for(row=0, p.row=p.begin; row<rows; row++, p.row+=p.rowstep) {
			p.colend=p.row+mtx.cols_;
			for(col=colf, p.col=p.row; col<=colc; col++, p.col++) {
				*(p.col) = col;
			}
		}
		return mtx;
	}
	Matrix Matrix::mesh(const Range rowrng, int cols) {
		int rowf = rowrng.int_floor(), rowc = rowrng.int_ceil();
		Matrix mtx(rowc-rowf+1, cols);
		elem_ptr p(mtx);
		int row, col;
		for(row=rowf, p.row=p.begin; row<=rowc; row++, p.row+=p.rowstep) {
			p.colend=p.row+mtx.cols_;
			for(col=0, p.col=p.row; col<cols; col++, p.col++) {
				*(p.col) = row;
			}
		}
		return mtx;
	}




	char MatrixExpression::OPERAND_SYMBOL[10] = {'N','S','+','-','*','/','+','-','*','/'};

	MatrixExpression Matrix::operator ()(int row, Range col) {
		Range rng;
		return MatrixExpression(1, col.int_ceil(cols_)-col.int_floor(1)+1, *this, row<=rng<=row, col, MatrixExpression::SLICE);
	}
	MatrixExpression Matrix::operator ()(Range row, int col) {
		Range rng;
		return MatrixExpression(row.int_ceil(rows_)-row.int_floor(1)+1, 1, *this, row, col<=rng<=col, MatrixExpression::SLICE);
	}
	MatrixExpression Matrix::operator ()(Range row, Range col) {
		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);
	}
	MatrixExpression operator +(const Matrix &lhs, double rhs) {
		return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::ADD_SC);
	}
	MatrixExpression operator -(const Matrix &lhs, double rhs) {
		return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::SUB_SC);
	}
	MatrixExpression operator *(const Matrix &lhs, double rhs) {
		return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::MUL_SC);
	}
	MatrixExpression operator /(const Matrix &lhs, double rhs) {
		return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::DIV_SC);
	}
	MatrixExpression operator ^(const Matrix &lhs, double rhs) {
		if(lhs.op_!=Matrix::MASK) {
			if(lhs.rows_==lhs.cols_) return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::POW);
			else throw Exception("Argument Error ^: This operator is under testing.");
		} else {
			return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::POW_MASK);
		}
	}
	MatrixExpression operator +(const Matrix &lhs, const Matrix &rhs) {
		if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::ADD);
		else throw Exception("Argument Error +: Matrix size was unmatched.");
	}
	MatrixExpression operator -(const Matrix &lhs, const Matrix &rhs) {
		if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, lhs.cols_, lhs, rhs, MatrixExpression::SUB);
		else throw Exception("Argument Error -: Matrix size was unmatched.");
	}
	MatrixExpression operator *(const Matrix &lhs, const Matrix &rhs) {
		if(lhs.op_!=Matrix::MASK && rhs.op_!=Matrix::MASK) {
			if(lhs.check_size_mul(rhs)) return MatrixExpression(lhs.rows_, rhs.cols_, lhs, rhs, MatrixExpression::MUL);
			else throw Exception("Argument Error *: Matrix size was unmatched.");
		} else {
			if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, rhs.cols_, lhs, rhs, MatrixExpression::MUL_MASK);
			else throw Exception("Argument Error *~: Matrix size was unmatched.");
		}
	}
	MatrixExpression operator ~(const Matrix &lhs) {
		Range rng;
		return MatrixExpression(lhs.rows_, lhs.cols_, lhs, 1<=rng<=lhs.rows_, 1<=rng<=lhs.cols_, MatrixExpression::MASK);
	}
//	MatrixExpression operator /(const Matrix &lhs, const Matrix &rhs) {
//		if(lhs.check_size_equal(rhs)) return MatrixExpression(lhs.rows_, rhs.cols_, lhs, rhs, MatrixExpression::DIV);
//		else throw Exception();
//	}

	MatrixExpression::MatrixExpression() : lhs_(NULL), rhs_(NULL) {}
	MatrixExpression::MatrixExpression(int rows, int cols, const Matrix &f, const Matrix &f2, OPERAND op)
		: lhs_(const_cast<Matrix*>(&f)), rhs_(const_cast<Matrix*>(&f2)) {
		if(op==ADD || op==SUB || op==MUL || op==DIV || op==MUL_MASK) {
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
		op_ = op;
		rows_ = rows;
		cols_ = cols;
		element_ = NULL;
	}
	MatrixExpression::MatrixExpression(int rows, int cols, const Matrix &f, double f2, OPERAND op)
		: lhs_(const_cast<Matrix*>(&f)), rhs_(NULL), scholar_(f2) {
		if(op==ADD_SC || op==SUB_SC || op==MUL_SC || op==DIV_SC || op==POW || op==POW_MASK) {
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
		op_ = op;
		rows_ = rows;
		cols_ = cols;
		element_ = NULL;
	}
	MatrixExpression::MatrixExpression(int rows, int cols, const Matrix& f, Range& rowrng, Range& colrng, OPERAND op)
		: lhs_(const_cast<Matrix*>(&f)), rhs_(NULL), rowrng_(rowrng), colrng_(colrng) {
		if(op==SLICE || op==MASK) {
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
		op_ = op;
		rows_ = rows;
		cols_ = cols;
		element_ = NULL;
	}
	MatrixExpression::MatrixExpression(int rows, int cols, Wave &roww, Wave &colw, OPERAND op)
		: lhs_(NULL), rhs_(NULL), meshrow_(&roww), meshcol_(&colw) {
		if(op==MESH) {
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Operad format was failed.");
		op_ = op;
		rows_ = rows;
		cols_ = cols;
		element_ = NULL;
	}
	MatrixExpression::MatrixExpression(const MatrixExpression &mtxc) : lhs_(mtxc.rhs_), rhs_(mtxc.rhs_) {
		op_ = mtxc.op_;
	}


	void MatrixExpression::swap(Matrix &rhs) {
		Matrix mtx;
		if(has_instance() && rhs.has_instance()) {
			mtx.set(rows_, cols_);
			mtx = *this;
			*this = rhs;
			rhs = mtx;
		}
		else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
	}
	void MatrixExpression::track(int n) const {
		for(int i=0; i<n; i++) std::cout << " ";
		std::cout << "Matrix op" << MatrixExpression::OPERAND_SYMBOL[op_] << ":" << this << std::endl;
		if(lhs_!=(Matrix*)this) lhs_->track(++n);
		if(rhs_!=(Matrix*)this) rhs_->track(n);
	}
	MatrixExpression::~MatrixExpression() {
		//std::cout << "DEST-MatrixExpression op" << MatrixExpression::OPERAND_SYMBOL[op_] << " " << rows_ << "*" << cols_ << ":" << this << std::endl;
	}


	bool MatrixExpression::has_instance() const {
		if(op_==SLICE || op_==MASK) return lhs_->has_instance();
		else return false;
	}
	Matrix* MatrixExpression::instance_ptr() const {
		if(op_==SLICE || op_==MASK) return lhs_->instance_ptr();
		else return NULL;
	}
	void MatrixExpression::get_element_ptr(double *&ptr_origin, double *&ptr_start, int &row_step) const {
		if(op_==SLICE || op_==MASK) {
			lhs_->get_element_ptr(ptr_origin, ptr_start, row_step);
			ptr_start = ptr_origin+(colrng_.int_floor(1)-1)+((rowrng_.int_floor(1)-1)*row_step);
			//col_width = cols_;
			return;
		} else throw Exception(typeid(*this), "MEMORY ERROR", "No instance.");
	}

	//void MatrixExpression::substitute(double rhs) {
	//	Matrix::substitute(rhs);
	//}
	//void MatrixExpression::substitute(const Matrix &mtx) {
	//	const_cast<Matrix &>(mtx).track_and_calculate(*this, true, true);
	//}

	void MatrixExpression::evaluate(Matrix& target) {
		Matrix *newself = NULL, *actualtarget = target.instance_ptr();
		int cnt_in_add=0, cnt_in_mul=0;
		if( track_self(actualtarget, cnt_in_add, cnt_in_mul, false) ) {
			newself = new Matrix(actualtarget->rows_, actualtarget->cols_);
			track_and_calculate(*newself, true, true);
			actualtarget->swap(*newself);
			delete newself;
		} else {
			target.substitute(0.0);
			track_and_calculate(target, true, true);
		}
	}
	void MatrixExpression::evaluate(MatrixExpression& target) {
		Matrix *newself = NULL, *actualtarget = target.instance_ptr();
		int cnt_in_add=0, cnt_in_mul=0;
		if( track_self(actualtarget, cnt_in_add, cnt_in_mul, false) ) {
			newself = new Matrix(*actualtarget);
			actualtarget = target.lhs_;
			target.lhs_ = newself;
			target = 0;
			track_and_calculate(target, true, true);
			actualtarget->swap(*newself);
			delete newself;
		} else {
			target.substitute(0.0);
			track_and_calculate(target, true, true);
		}
	}
	bool MatrixExpression::track_self(const Matrix* self, int &cnt_in_add, int &cnt_in_mul, bool in_mul) {
		bool found = false, in_mult = in_mul;
		switch(op_) {
//			case ADD:
//			case SUB:
//			case ADD_SC:
//			case SUB_SC:
//				break;
			case MUL:
			case DIV:
			case MUL_SC:
				in_mult = true;
				break;
		}
		if(lhs_!=NULL) found = found || lhs_->track_self(self,cnt_in_add,cnt_in_mul,in_mult);
		if(rhs_!=NULL) found = found || rhs_->track_self(self,cnt_in_add,cnt_in_mul,in_mult);
		return found;
	}
	void MatrixExpression::track_and_calculate(Matrix& target, bool add_posi, bool mul_posi) {
		Matrix *buf1, *buf2;
		bool made_buf1=false, made_buf2=false;
		if(op_==ADD) {
			if(lhs_->has_instance()) target.add(*lhs_, add_posi);
			else lhs_->track_and_calculate(target, add_posi, mul_posi);
			if(rhs_->has_instance()) target.add(*rhs_, add_posi);
			else rhs_->track_and_calculate(target, add_posi, mul_posi);
		} else if(op_==SUB) {
			if(lhs_->has_instance()) target.add(*lhs_, add_posi);
			else lhs_->track_and_calculate(target, add_posi, mul_posi);
			if(rhs_->has_instance()) target.add(*rhs_, !add_posi);
			else rhs_->track_and_calculate(target, !add_posi, mul_posi);
		} else if(op_==MUL) {
			if(lhs_->has_instance()) buf1=lhs_;
			else {
				buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
				made_buf1 = true;
				lhs_->track_and_calculate(*buf1, true, true);
			}
			if(rhs_->has_instance()) buf2=rhs_;
			else {
				buf2 = new Matrix(rhs_->rows_, rhs_->cols_);
				made_buf2 = true;
				rhs_->track_and_calculate(*buf2, true, true);
			}
			target.mul(*buf1, *buf2, add_posi);
			if(made_buf1) delete buf1;
			if(made_buf2) delete buf2;
		} else if(op_==ADD_SC) {
			if(lhs_->has_instance()) target.add(*lhs_, add_posi);
			else lhs_->track_and_calculate(target, add_posi, mul_posi);
			target.add(scholar_, add_posi);
		} else if(op_==SUB_SC) {
			if(lhs_->has_instance()) target.add(*lhs_, add_posi);
			else lhs_->track_and_calculate(target, add_posi, mul_posi);
			target.add(scholar_, !add_posi);
		} else if(op_==MUL_SC) {
			if(lhs_->has_instance()) buf1=lhs_;
			else {
				buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
				made_buf1 = true;
				lhs_->track_and_calculate(*buf1, true, true);
			}
			target.mul(*buf1, scholar_, add_posi);
			if(made_buf1) delete buf1;
		} else if(op_==DIV_SC) {
			if(lhs_->has_instance()) buf1=lhs_;
			else {
				buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
				made_buf1 = true;
				lhs_->track_and_calculate(*buf1, true, true);
			}
			target.div(*buf1, scholar_, add_posi);
			if(made_buf1) delete buf1;
		} else if(op_==MUL_MASK) {
			if(lhs_->has_instance()) buf1=lhs_;
			else {
				buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
				made_buf1 = true;
				lhs_->track_and_calculate(*buf1, true, true);
			}
			if(rhs_->has_instance()) buf2=rhs_;
			else {
				buf2 = new Matrix(rhs_->rows_, rhs_->cols_);
				made_buf2 = true;
				rhs_->track_and_calculate(*buf2, true, true);
			}
			target.mul_mask(*buf1, *buf2, add_posi);
			if(made_buf1) delete buf1;
			if(made_buf2) delete buf2;
		} else if(op_==SLICE) {
			throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
//			if(lhs_->has_instance()) buf1=lhs_;
//			else {
//				buf1 = new Matrix(lhs_->rows_, lhs_->cols_);
//				made_buf1 = true;
//				lhs_->track_and_calculate(*buf1, true, true);
//			}
//			target.div(*buf1, scholar_, add_posi);
//			if(made_buf1) delete buf1;
		}
	}

	double& MatrixExpression::operator ()(int row, int col) const {
		elem_ptr elm(*this);
		if(op_==SLICE) {
			if(elm.begin==NULL) throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
			return *(elm.begin + ((row-1)*elm.rowstep) + (col-1));
		} else {
			throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
		}
	}
	MatrixExpression MatrixExpression::operator ()(int row, Range col) { return (*this).Matrix::operator ()(row, col); }
	MatrixExpression MatrixExpression::operator ()(Range row, int col) { return (*this).Matrix::operator ()(row, col); }
	MatrixExpression MatrixExpression::operator ()(Range row, Range col) { return (*this).Matrix::operator ()(row, col); }

	Matrix& MatrixExpression::operator =(double rhs) {
		if(has_instance()) substitute(rhs);
		else throw Exception();
		return *this;
	}
	Matrix& MatrixExpression::operator =(const Matrix &rhs) {
		if(check_size_equal(rhs)) {
			if(has_instance()) {
				substitute(rhs);
			} else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
		return *this;
	}
	Matrix& MatrixExpression::operator =(const MatrixExpression &rhs) {
		if(check_size_equal(rhs)) {
			if(has_instance()) {
				if(rhs.has_instance()) substitute(rhs);
				else const_cast<MatrixExpression &>(rhs).evaluate(*this);
			} else throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
		} else throw Exception(typeid(*this), "EXPRESSION ERROR", "Matrix size did not matched");
		return *this;
	}


	Matrix & MatrixExpression::slide(int drow, int dcol) {
		throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
		return *this;
	}
	Matrix & MatrixExpression::transpose() {
		throw Exception(typeid(*this), "NO MEMORY ERROR", "Memory not allocated");
		return *this;
	}

	std::ostream& operator<<(std::ostream& ost, const Matrix &mtx) {
		if(mtx.has_instance()) {
			for(int row=1; row<=mtx.rows_; row++) {
				for(int col=1; col<=mtx.cols_; col++) {
					ost << mtx(row,col) << "\t";
				}
				ost << std::endl;
			}
		}
		return ost;
	}


namespace MatrixOverload {
Matrix sin(const Matrix &m) {
	Matrix mtx(m);
	return mtx.each(&::sin);
}
Matrix cos(const Matrix &m) {
	Matrix mtx(m);
	return mtx.each(&::cos);
}
Matrix tan(const Matrix &m) {
	Matrix mtx(m);
	return mtx.each(&::tan);
}
Matrix exp(const Matrix &m) {
	Matrix mtx(m);
	return mtx.each(&::exp);
}Matrix sqrt(const Matrix &m) {
	Matrix mtx(m);
	return mtx.each(&::sqrt);
}
Matrix pow(const Matrix &m, double ex) {
	Matrix mtx(m);
	return mtx.each(&::pow, ex);
}


Matrix max(const Matrix& a, const Matrix b) {
	if(!a.check_size_equal(b)) throw Exception(typeid(Matrix), "EXPRESSION ERROR", "Matrix size did not matched");
	Matrix tmp(a.getRows(), a.getCols());
	for(int col=0; col<a.getCols(); col++) {
		for(int row=0; row<a.getRows(); row++) {
			tmp(row, col) = ( a(row, col) > b(row, col) ? a(row, col) : b(row, col) );
		}
	}
	return tmp;
}

Matrix min(const Matrix& a, const Matrix b) {
	if(!a.check_size_equal(b)) throw Exception(typeid(Matrix), "EXPRESSION ERROR", "Matrix size did not matched");
	Matrix tmp(a.getRows(), a.getCols());
	for(int col=0; col<a.getCols(); col++) {
		for(int row=0; row<a.getRows(); row++) {
			tmp(row, col) = ( a(row, col) < b(row, col) ? a(row, col) : b(row, col) );
		}
	}
	return tmp;
}
}

}	/*	<- namespace Psycholops 	*/


