#ifndef MATRIX_H
#define MATRIX_H
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
#include <assert.h>
#include <math.h>
#include <iostream>

template<class T>
class Row {
	public:
	Row(int cols=0):row(NULL), numCols(0) {SetRowSize(cols);}
	~Row() {SetRowSize(0); }
	Row(const Row &r):row(NULL), numCols(0) {
		SetRowSize(r.numCols);
		for (int i=0;i<numCols;i++)
			row[i]=r.row[i];
	}
	
	void SetRowSize(int n) {
		if(row && n != numCols) delete [] row;
		if (n>0) {
			if (n != numCols) row=new T[n];
 			memset(row,0,sizeof(T)*n);
		}
		else row=NULL;
		numCols=n;
	}
	void zero() {
		memset(row,0,sizeof(T)*numCols);
	}
	void load(FILE *fp)
	{
		char line[1500];
		char *p;
		int r;
		int i,j;
		int binary=0;
		if (!fp) return;
		fgets(line,1500,fp);
		r = atoi(line);
		SetRowSize(r);
		i=0;
		while (fgets(line,1500,fp))
		{
			row[i] = atof(line);
			i++;
		}
	}
	T& operator[](int column) {
// 		assert(column<numCols);
		return row[column];
	}
	Row& operator=(const Row& r) {
		SetRowSize(r.numCols);
		for (int i=0;i<numCols;i++)
			row[i]=r.row[i];
		return *this;
	}
	Row<T> operator+ (const Row &r) const
	{
		Row<T> s(numCols);
		int i;
		assert(r.numCols == numCols);
		for (i=0;i<numCols;i++) s.row[i] = r.row[i] + row[i];
		return s;
	}
	Row<T> operator- (const Row &r) const
	{
		Row<T> s(numCols);
		int i;
		assert(r.numCols == numCols);
		for (i=0;i<numCols;i++) s.row[i] = row[i] - r.row[i];
		return s;
	}
	Row<T> operator* (double v) const
	{
		Row<T> s(numCols);
		int i;
		for(i=0;i<numCols;i++) s.row[i] = row[i]*v;
		return s;
	}
	void operator*= (double v) const
	{
		int i;
		for(i=0;i<numCols;i++) row[i] = row[i]*v;
	}
	T operator* (const Row &r) const
	{
		int i;
		T d = 0;
		assert(r.numCols == numCols);
		for (i=0;i<numCols;i++) d += row[i] * r.row[i];
		return d;
	}
	void operator+= (const Row &r)
	{
		int i;
		assert(r.numCols == numCols);
		for (i=0;i<numCols;i++) row[i] += r.row[i];
	}
	int size() { return numCols;}
	T rowsum() {
		int i;
		T sum=0;
		for (i=0;i<numCols;i++) sum += fabs((double)row[i]);
		return sum;
	}
	T* row;

	private:
	int numCols;
};

template<class T>
class MatrixRow {
	// this is just a pointer to an array of bytes in a matrix. It has no brain of its own
	public:
		MatrixRow() {
			data = 0;
		}
		void set_start(T *n, int c) {
			data = n;
			numCols = c;
		}
		T& operator[](int index) { return data[index]; }
		T *datapointer() { return data; }
		T rowsum() {
			int i;
			T sum=0;
			for (i=0;i<numCols;i++) sum += fabs(data[i]);
			return sum;
		}
	private:
		T *data;
		int numCols;
};

template<class T>
class Matrix {
	public:
		Matrix(int rows=0, int cols=0): matrix(NULL), data(NULL), hashed(0) {
			numRows=numCols=0;
			SetSize(rows,cols);
		}
		Matrix(const Matrix& m): matrix(NULL), data(NULL) {
			numRows=numCols=0;
			SetSize(m.numRows,m.numCols);
			memcpy(data,m.data,numRows*numCols*sizeof(T));
		}
		~Matrix() {
			// 			
			delete [] matrix;
			delete [] data;
		}
		void SetSize(int rows, int cols) {
			if (rows == numRows && cols == numCols) return;
			if (matrix) delete [] matrix;
			if (data) delete [] data;
			matrix=new MatrixRow<T>[rows];
			data = new T[rows*cols];
			for (int i=0;i<rows;i++)
				matrix[i].set_start(data+i*cols, cols);
			numCols=cols;numRows=rows;
		}
		void identity() {
			if (numCols != numRows) return;
			memset(data,0,sizeof(T)*numRows*numCols);
			for (int i=0;i<numRows;i++) {
				data[i+i*numCols]=1;
			}
		}
		void zero() {
			bzero(data,sizeof(T)*numRows*numCols);
		}
		int Cols() const { return numCols;}
		int Rows() const { return numRows;}
		MatrixRow<T>& operator[](int index) {
		// 			assert(index<numRows);
			return matrix[index];
		}
		T raw(int row, int column) {
			return data[ row * numCols + column ];
		}
		Matrix& operator=(const Matrix& m) {
			SetSize(m.numRows,m.numCols);
			memcpy(data,m.data,numRows*numCols*sizeof(T));
			return *this;
		}
		Matrix operator+(const Matrix& m) {
	// return a new matrix that is the sum of these two matrices
			assert(m.numRows==numRows&&m.numCols==numCols);
			Matrix p(numRows,numCols);
			int i;
			for (i=0;i<numRows*numCols;i++) p.data[i] = data[i] + m.data[i];
			return p;
		}

		Matrix operator-(const Matrix& m) {
	// return a new matrix that is the sum of these two matrices
			assert(m.numRows==numRows&&m.numCols==numCols);
			Matrix p(numRows,numCols);
			int i;
			for (i=0;i<numRows*numCols;i++) p.data[i] = data[i] - m.data[i];
			return p;
		}

		void operator+=(const Matrix &m) {
			int i;
			assert(m.numRows==numRows && m.numCols==numCols);
			for (i=0;i<numRows*numCols;i++) data[i] += m.data[i];
		}
		Row<T> column(int index) const {
			int i;
			assert(index<numCols);
			Row<T> s(numRows);
			for(i=0;i<numRows;i++) s.row[i] = data[index+i*numCols];
			return s;
		}
		T column_product(int index, MatrixRow<T> &r) const {
			int i;
			assert(index<numCols);
			T sum = 0;
			for (i=0;i<numRows;i++) sum += data[index+i*numCols] * r[i];
			return sum;
		}
		
		Matrix<T> interpolate(Matrix &m, double slope)
		{
			Matrix<T> s(numRows,numCols);
			int i;
			for (i=0;i<numRows*numCols;i++) s.data[i] = (data[i] - m.data[i])*slope + m.data[i];
			return s;
		}

		Matrix<T> operator*(const Matrix &m) {
			assert(m.numRows==numCols);
			Matrix<T> s(numRows,m.numCols);
			int i,j;
			for (i=0;i<s.numRows;i++)
				for (j=0;j<s.numCols;j++) 
					s[i][j] = m.column_product(j, matrix[i]);
			return s;
		}
		T norm()
		{
			T max=-1e40;
			int i;
			for (i=0;i<numRows;i++) {
				if (max < matrix[i].rowsum())
					max = matrix[i].rowsum();
			}
			return max;
		}
		Matrix<T> operator*(double v) {
			int i;
			Matrix<T> s(numRows,numCols);
			for(i=0;i<numRows*numCols;i++) { s.data[i] = data[i]*v; }
			return s;
		}
		void operator*=(double v) {
			int i;
			for (i=0;i<numRows*numCols;i++) { data[i] *= v; }
		}
		Matrix<T> taylor(double t)
		{
	// calculates e^Qt
			int i;
			int s = Cols();
			Matrix<T> sum(Rows(),Cols());
			Matrix<T> I(Rows(),Cols());
			Matrix<T> p(Rows(),Cols());

			for (i=0;i<s;i++) { I[i][i] = 1; p[i][i] = 1; }

			sum += I;
			for (i=1;i<t*30;i++)
			{
				p = p * (*this * (t/i));
				sum += p;
			}
			return sum;
		}
		Matrix<T> pade(double t)
		{
			T c=1;
			Matrix<T> D(Rows(),Cols()), N(Rows(),Cols()), X(Rows(),Cols());
			int j = 1 + (int)floor(log(norm())/log(2.0));
			if (j < 0) j = 0;
			int q = 20;
			Matrix<T> A(*this);
			A *= t/pow(2.0,j);
			D.identity();
			N.identity();
			X.identity();
			int k, ik = -1;
			for (k=1;k<q+1;k++)
			{
				c = c*(q - k + 1)/((2*q-k+1)*k);
				X = A * X; N += X * c; D += X * (c * ik);
				ik = ik * -1;
			}
			Matrix<T> F = D.inverse() * N;
			for (k=1;k<=j;k++) {
				F = F * F;
			}
			return F;
		}

		Matrix<T> transpose()
		{
			Matrix<T> trans(Cols(), Rows());

			int i,j;
			for (i=0;i<Rows();i++)
				for(j=0;j<Cols();j++)
					trans[j][i] = matrix[i][j];
			return trans;
		}

		Matrix<T> inverse()
		{
			T tmp;
			int i,j,k,retval;
			int n = Rows();

			retval = 0;
			Matrix<T> inv(Rows(),Cols());
			Matrix<T> orig(*this);
			inv.identity();

			for (k=0;k<n;k++) {
				tmp = orig[k][k];

				for (j=0;j<n;j++) {
			// Don't bother with previous entries
					if (j>k) orig[k][j] /= tmp;
					inv[k][j] /= tmp;
				}
				for (i=0;i<n;i++) {             // Loop over rows
					if (i == k) continue;
					tmp = orig[i][k];
					for (j=0;j<n;j++) {
						if (j>k) orig[i][j] -= orig[k][j]*tmp;
						inv[i][j] -= inv[k][j]*tmp;
					}
				}
			}
	
	// Copy inverse to source matrix
			return inv;
		}


		void print()
		{
			save(stdout);
		}

		void save(FILE *fp)
		{
			int i,j;
			int s=numCols;
			int t=numRows;
			fprintf(fp,"%d %d\n", numCols, numRows);
			for (i=0;i<t;i++) {
				for (j=0;j<s;j++) {
					fprintf(fp,"%e\t",matrix[i][j]);
				}
				fprintf (fp,"\n");
			}
		}
		void bytesave(FILE *fp)
		{
			int i, j;
			int s = numCols;
			int t = numRows;
			double bt;
			fprintf(fp,"b%d %d\n", numCols, numRows);
			for (i=0;i<t;i++) {
				for (j=0;j<s;j++) {
					bt = matrix[i][j];
					fwrite(&bt,sizeof(double),1,fp);
				}
			}
		}
		void byteload(FILE *fp)
		{
			double b;
			int i,j;
			double *block = new double[numRows*numCols];
			fread(block,sizeof(double),numRows*numCols,fp);
			for (i=0;i<numRows;i++)
				for (j=0;j<numCols;j++) {
				matrix[i][j]=block[j+i*numCols];
				}
				delete block;
		}
		void load(FILE *fp)
		{
			char line[150000];
			char *p;
			int c, r;
			int i,j;
			int binary=0;
			if (!fp) return;
			fgets(line,150000,fp);
			p = strtok(line," ");
			if (p[0] == 'b') {
				c = atoi(p+1);
				binary=1;
			} else
				c = atoi(p);
				p = strtok(NULL,"\n");
				r = atoi(p);
				SetSize(r,c);
				if (binary) {
					byteload(fp);
					return;
				}
				i=0;
				while (fgets(line,150000,fp))
				{
					p = strtok(line,"\t");
					for (j=0;j<c;j++) {
						matrix[i][j] = atof(p);
						p = strtok(NULL,"\t");
					}
					i++;
				}
		}
		int hashed;

	private:
		int numCols, numRows;
		MatrixRow<T>* matrix;
		T *data;

};
		
#endif