/*  * @(#)Matrix.java	1.3	20/9/96 Akio Utsugi */import java.lang.Math;import java.util.Random;/** * Matrix is a class for matrix objects with basic linear calculation methods. * @version 1.2 17 Sep 1996 * @author <a href="http://www.aist.go.jp/NIBH/~b0616/"> Akio Utsugi </a> */public class Matrix{    final static double small = 1e-8;        int row, col;    double value[][];/** * Construct a square matrix object. * @param n size */    public Matrix(int n)    {	this(n, n);    }    /** * Construct a rectangle matrix object. * @param row number of row * @param col number of column */    public Matrix(int row, int col)    {	int i, j;		this.row = row;	this.col = col;	value = new double[row][col];	for(i=0; i<row; i++){	    for(j=0; j<col; j++){		value[i][j] = 0;	    }	}	    }/* * Creation of special matrices */    /** * Make an identity matrix. * @param n size */    static public Matrix identity(int n)    {	int i;		Matrix ret = new Matrix(n);	for(i=0; i<n; i++){	    ret.value[i][i] = 1;	}	return ret;    }/**     * Make a standard Gaussian random Matrix. * @param row number of row * @param col number of column * @param seed random seed */    static public Matrix random(int row, int col, int seed)    {	int i, j;		Random rnd = new Random(seed);	Matrix ret = new Matrix(row, col);	for(i=0; i<row; i++){	    for(j=0; j<col; j++){		ret.value[i][j] = rnd.nextGaussian();	    }	}	return ret;    }/** * Same as random(row, col, 1). */       static public Matrix random(int row, int col)    {	return random(row, col, 1);    }    /** * Copy the matrix. */    public Matrix copy()    {	int i, j;		Matrix ret = new Matrix(row, col);	for(i=0; i<row; i++){	    for(j=0; j<col; j++){		ret.value[i][j] = value[i][j];	    }	}		return ret;    }/* * Basic Matrix Calculations *//** * Transpose: X'. */    public Matrix transpose()    {        int i, j;                Matrix ret = new Matrix(col, row);        for(i=0; i< row; i++){            for(j=0; j< col; j++){                ret.value[j][i] = value[i][j];            }        }        return ret;    }    /** * Linear conversion: a*X + b*Y. */    public Matrix linearConv(double a, double b, Matrix Y)    {	int i, j;		int Y_row = Y.row;	int Y_col = Y.col;		if(row == Y_row && col == Y_col){	    Matrix ret = new Matrix(row, col);	    	    for(i=0; i < row; i++){		for(j=0; j < col; j++){		    ret.value[i][j] = a*value[i][j] + b*Y.value[i][j];		}	    }	    return ret;	}	else{	    return null;	}    }            /** * Update by Linear conversion: X = a*X + b*Y. */    public void updateLinearConv(double a, double b, Matrix Y)    {	int i, j;		int Y_row = Y.row;	int Y_col = Y.col;		for(i=0; i < row; i++){	   for(j=0; j < col; j++){	       value[i][j] = a*value[i][j] + b*Y.value[i][j];	   }	}    }    /** * Scalar product: aX. */    public Matrix multipliedBy(double a)    {	int i, j;		Matrix ret = new Matrix(row, col);	for(i=0; i< row; i++){	    for(j=0; j< col; j++){		ret.value[i][j] = a*value[i][j];	    }	}	return ret;    }/** * Matrix multiplication: X*Y. */    public Matrix multipliedBy(Matrix Y)    {	int i, j, k;	double s;		int Y_row = Y.row, Y_col = Y.col;		if(col != Y_row){	    return null;	}		Matrix ret = new Matrix(row, Y_col);		for(i=0; i< row; i++){	    for(j=0; j<Y_col; j++){		s = 0;		for(k=0; k<col; k++){		    s += value[i][k] * Y.value[k][j];		}		ret.value[i][j] = s;	    }	}		return ret;    }    /** * Matrix division: X\Y. */    public Matrix dividedBy(Matrix Y)    {	return this.dividedBy(Y, 2*col-1);    }/** * Matrix division with bandwidth: X\Y.  */    public Matrix dividedBy(Matrix Y, int bw)    {    	Matrix L, U, ret;	int n, m, i, j, k;	double w, s, pp;		L = this.copy();	U = Y.copy();		m = L.row;	n = U.col;	bw = bw/2 + 1;		for(k = 0; k < m; k++){	    pp=L.value[k][k];	    	    if(Math.abs(pp)< small){		return null;	    }	    w = 1/pp;	    	    for(j = k+1; j < k+bw && j < m ; j++){		s = w*L.value[k][j];		L.value[k][j] = s;		for(i = k+1; i < k+bw+1 && i < m; i++){		    s = L.value[i][j] - L.value[i][k]*L.value[k][j];		    L.value[i][j] = s;		}	    }	    for(j = 0; j < n; j++){		s=w*U.value[k][j];		U.value[k][j] = s;		for(i = k+1; i < k+bw && i < m; i++){		    s=U.value[i][j] - L.value[i][k] * U.value[k][j];		    U.value[i][j] = s;		}	    }	}		ret = new Matrix(m, n);		for(i = m; --i>=0;){	    for(j = 0; j < n; j++){		w=U.value[i][j];				for(k = i+1; i < k+bw && k < m; k++){		    w -= L.value[i][k] * ret.value[k][j];		}		ret.value[i][j] = w;	    }		}		return ret;    }/* * Special purpose calculations *//** * Cross squared Euclidean distance Matrix: Z_{ij} = ||X_i - Y_j||^2  */    public Matrix crossSqDistance(Matrix Y)    {	Matrix ret;	int i, j, k;	double s, d;	int Y_row = Y.row, Y_col = Y.col;		if(col != Y_col){	    return null;	}		ret = new Matrix(row, Y_row);		for(i=0; i<row; i++){	    for(j=0; j<Y_row; j++){		s = 0;		for(k=0; k<col; k++){		    d = value[i][k] - Y.value[j][k];		    s += d*d;		}		ret.value[i][j] = s;	    }	}		return ret;    }    /** * Update by exponential : exp(X_{ij}) */    public void updateExp()    {	int i, j;		for(i=0; i<row; i++){	    for(j=0; j<col; j++){		value[i][j] = Math.exp(value[i][j]);	    }	}    }    /** * Make a horizontal-summation vector. */    public Matrix horizontalSum()    {	int i, j;	double s;		Matrix ret = new Matrix(row, 1);		for(i=0; i<row; i++){	    s = 0;	    for(j=0; j<col; j++){		s += value[i][j];	    }	    ret.value[i][0] = s;	}	return ret;    }/** * Make a vertical-summation vector. */    public Matrix verticalSum()    {	int i, j;	double s;		Matrix ret = new Matrix(1, col);		for(j=0; j<col; j++){	    s = 0;	    for(i=0; i<row; i++){		s += value[i][j];	    }	    ret.value[0][j] = s;	}	return ret;    }/** * Summation of all entries. */        public double sumEntries()    {        int i, j;        double s = 0;                for(i=0; i<row; i++){            for(j=0; j<col; j++){                s += value[i][j];            }        }        return s;    }   /** * Summation of squared entries. */     public double sumSqrEntries()    {        int i, j;        double s = 0;                        for(i=0; i<row; i++){            for(j=0; j<col; j++){                s += value[i][j]*value[i][j];            }        }        return s;    }/** * Multiply entries. */   public Matrix mulipliedEntriesWith(Matrix Y)   {       Matrix ret;       int i, j;       ret = new Matrix(row, col);       for(i=0; i<row; i++){	   for(j=0; j<col; j++){	       ret.value[i][j] = value[i][j]*Y.value[i][j];	   }       }       return ret;    }                  /**  * X_{ij} = X_{ij}/v_{j}. */    public void updateDivideByRowVector(Matrix vec)    {	int i, j;		for(j=0; j<col; j++){	    for(i=0; i<row; i++){		value[i][j] = value[i][j]/vec.value[0][j];	    }	}    }/** * If the matrix is a row or column vector, make a diagonal matrix from it. * Otherwise extract its diagonal as a row vector. */    public Matrix diagonal()    {	Matrix ret;	int i, n;	if(row == 1){	    ret = new Matrix(col, col);	    for(i=0; i< col; i++){		ret.value[i][i] = value[0][i];	    }	}	else if(col == 1){	    ret = new Matrix(row, row);	    for(i=0; i< row; i++){		ret.value[i][i] = value[i][0];	    }	}	else{	    n = Math.min(row, col);	    ret = new Matrix(1, n);	    for(i=0; i< n; i++){		ret.value[0][i] = value[i][i];	    }	}	return ret;    }/** * Make a row vector including the eigenvalues of the matrix. * @param eps allowable absolute value of off-diagonal entries * @lmax maximum number of repeat */    public Matrix eigenvalues(double eps, int lmax)    {	int n;	int i, j, im = 0, jm = 0, ll;	double pmax, pii, pjj, pij, d, tn;	double a, b, p1, p2, q1, q2, dd;	Matrix px, qx, lmd;	px = this.copy();	n = px.row;	qx = Matrix.identity(n);	ll=0;	do{	    pmax = Math.abs(px.value[0][1]);	    ll++;	    for(i=0; i<n-1; i++){		for(j=i+1; j<n; j++){		    if(pmax<= Math.abs(px.value[i][j])){			pmax= Math.abs(px.value[i][j]);			im=i;			jm=j;		    		    }		}	    }	    if(pmax<=eps)	break;	    pii=px.value[im][im];	    pjj=px.value[jm][jm];	    pij=px.value[im][jm];	    d=pii-pjj;	    dd=d*d+4.*pij*pij;	    tn=2.*pij/(d+Math.sqrt(dd)*(d > 0 ? 1 : -1));	    a=1./Math.sqrt(1.+tn*tn);	    b=a*tn;	    for(i=0; i<n; i++){		p1=px.value[i][im];		p2=px.value[i][jm];		px.value[i][im] = p1*a+p2*b;		px.value[i][jm] =  -p1*b+p2*a;		q1=qx.value[i][im];		q2=qx.value[i][jm];		qx.value[i][im] = q1*a+q2*b;		qx.value[i][jm] = -q1*b+q2*a;	    	    }	    for(i=0; i<n; i++){		p1=px.value[im][i];		p2=px.value[jm][i];		px.value[im][i] = p1*a+p2*b;		px.value[jm][i] = -p1*b+p2*a;	    	    }	}while(ll<lmax);	lmd = px.diagonal();	return lmd;    }}