matrix.java,v

包含了模式识别中常用的一些分类器设计算法

head	1.6;
access;
symbols;
locks; strict;
comment	@# @;


1.6
date	2005.06.10.18.48.09;	author rirwin;	state Exp;
branches;
next	1.5;

1.5
date	2005.05.24.15.37.02;	author rirwin;	state Exp;
branches;
next	1.4;

1.4
date	2005.03.11.22.04.38;	author patil;	state Exp;
branches;
next	1.3;

1.3
date	2005.01.25.15.31.01;	author patil;	state Exp;
branches;
next	1.2;

1.2
date	2005.01.20.02.41.42;	author patil;	state Exp;
branches;
next	1.1;

1.1
date	2004.12.28.00.04.32;	author patil;	state Exp;
branches;
next	;


desc
@@


1.6
log
@Establishing RCS version.
@
text
@/**
 * file: Matrix.java
 *
 * last edited: Ryan Irwin
 */


// import java packages
//

//These imports are not needed - Phil T. 6-23-03
//import java.io.*;
//import java.lang.*;

import java.util.*;

/**
 * this class represents a matrix object and performs matrix operations 
 * that are needed to compute the eigenvectors and eigenvalues
 *
 */
public class Matrix 
{

    // *********************************************************************
    //
    // declare global variables and components
    //
    // *********************************************************************

    //  enum MTYPE { FULL = 0, DIAGONAL, SYMMETRIC,
    //LOWER_TRIANGULAR, UPPER_TRIANGULAR, SPARSE,
    //	       UNCHANGED, UNKNOWN, DEF_MTYPE = FULL };
    
    public static final int FULL = 0;
    public static final int DIAGONAL = 1;
    public static final int SYMMETRIC = 2;
    public static final int LOWER_TRIANGULAR = 3;
    public static final int UPPER_TRIANGULAR = 4;
    public static final int SPARSE = 5;

    // declare global variables
    //
    protected int row;
    protected int col;
    protected int type_d = FULL;
    protected double Elem[][];

    // *********************************************************************
    //
    // declare class methods
    //
    // *********************************************************************
    
    /**
     * Gets the number of Rows
     *
     * @@return   number of rows in the matrix
     */
    public int getNumRows() 
    {
	return row;
    }
 
    /**
     * Gets the number of Columns
     *
     * @@return   number of columns in the matrix
     *
     */
    public int getNumColumns() 
    {
	return col;
    }
 
    /**
     * this routine computes the inverse of a matrix
     *
     * @@param    A inverse matrix
     */
    public void copyLowerMatrix(Matrix A) 
    {
	
	// the only really clean way to do this is to manually copy all
	// the non-zero elements. for a truly sparse matrix, this won't be
	// prohibitively expensive.
	//
	for (int i = 0; i < A.row; i++) 
	{
	    for (int j = 0; j <= i; j++)
	    {
		
		// since looking up this value is expensive, we want to save it
		//
		Elem[i][j] = A.getValue(i, j);
	    }
	}	
    }
  
    /**
     * this routine computes the sum of a column
     *
     * @@param    col_a column number to be summed
     *
     * @@return   double value of sum
     *
     */
    public double getColSumSquare(int col_a) 
    {
	
	// the only really clean way to do this is to manually copy all
	// the non-zero elements. for a truly sparse matrix, this won't be
	// prohibitively expensive.
	//
	double sum = 0.0;
	for (int i = 0; i < row; i++) 
	{
	    sum += Elem[i][col_a] * Elem[i][col_a];
	}
	return sum;
    }
  
    /**
     * sets the value given row and column index to the double value passed
     *
     * @@param    r number of rows in the matrix
     * @@param    c number of columns in the matrix
     * @@param    value_a input set of points
     *
     */
    public void setValue(int r, int c, double value_a) 
    {
	Elem[r][c] = value_a;
    }
 
    /**
     * gets value at indexed row and column
     *
     * @@param    r number of rows in the matrix
     * @@param    c number of columns in the matrix
     *
     * @@return   double value at index
     */
    public double getValue(int r, int c) 
    {	
	return Elem[r][c];
    }
 
    /**
     * this method initialize the matrix to all r x n elements being the
     * double value passed
     *
     * @@param    r number of rows in the matrix
     * @@param    c number of columns in the matrix
     * @@param    value_a double value to initialize all elements of matrix
     * @@param    type_a integer value of type
     *
     *
     */
    public void initMatrixValue(int r, int c, double value_a, int type_a) 
    {
	row = r;
	col = c;
	type_d = type_a;

	Elem = new double[row][col];

	for(int i = 0; i < row; i++) 
	{
	    for(int j = 0; j < col; j++) 
	    {
		Elem[i][j] = value_a;
	    }
	}
    }
    
    /**
     * this method initializes matrix to be the values of the 2-d array passed
     *
     * @@param    initVal 2-dimensional array to be copyied
     * @@param    r number of rows in the matrix
     * @@param    c number of columns in the matrix
     */
    public void initMatrix(double initVal[][], int r, int c) 
    {
	
	row = r;
	col = c;

	Elem = new double[row][col];

	for(int i = 0; i < row; i++) 
	{
	    for(int j = 0; j < col; j++) 
	    {
		Elem[i][j] = initVal[i][j];
	    }
	}
    }
 
    /**
     * Creates a 1 row matrix of values specified by the input vector
     *
     * @@param    vec_a values to be put in matrix
     *
     */
    public void initMatrix(Vector<Double> vec_a) 
    {
	
	row = 1;
	col = vec_a.size();

	Elem = new double[row][col];

	for(int i = 0; i < row; i++) 
	{
	    for(int j = 0; j < col; j++) 
	    {
		Elem[i][j] = MathUtil.doubleValue(vec_a, j);
	    }
	}
    }
    
    /**
     * Creates a diagonal matrix of values specified by the input vector
     *
     * @@param    vec_a values to be put in matrix
     *
     */
    public void initDiagonalMatrix(Vector<Double> vec_a) 
    {
	
	row = vec_a.size();
	col = row;
	Elem = new double[row][col];

	for(int i = 0; i < row; i++) 
	{
	    Elem[i][i] = MathUtil.doubleValue(vec_a, i);
	}
    }
 
    /**
     * Puts the diagonal of a matrix in the vector passed
     *
     * @@param    vec_a Vector for values to be added from diagonal of a matrix
     *
     */
    public void getDiagonalVector(Vector<Double> vec_a) 
    {
	
	vec_a.clear();

	for(int i = 0; i < row; i++) 
	{
	    vec_a.add(new Double(Elem[i][i]) );
	}
    }
 
    /**
     * Given a column or row matrix the values can be put into a vector
     *
     * @@param    vec_a Vector for values to be added from column or row matrix
     *
     */
    public void toDoubleVector(Vector<Double> vec_a) 
    {

	int size = 0;

	if ( col == 1 )
	{
	    size = row;
	} 
	else if ( row == 1 )
	{
	    size = col;
	}
	else 
	{
	    // System.out.println("Error in toDoubleVector");
	    Exception e = new Exception();
	    e.printStackTrace();
	}

	vec_a.setSize(size);
	for(int j = 0; j < size; j++) 
	{
	    if ( row == 1 )
	    {
		vec_a.set(j, new Double(Elem[0][j]));
	    }
	    else 
	    {
		vec_a.set(j, new Double(Elem[j][0]));
	    }
	}
    }

    /** 
     * This routine takes a matrix as an argument and copies it to the
     * the matrix object that invoked it
     *
     * @@param    A input matrix to be copyied
     */
    public void copyMatrix(Matrix A) 
    {

	row = A.row; 
	col = A.col;
	
	Elem = new double[row][col];
	
	for(int i = 0; i < row; i++ ) 
	{
	    for(int j = 0; j < col; j++) 
	    {
		Elem[i][j] = A.Elem[i][j];
	    }
	}
    }

    /**
     * This routine takes a matrix as an argument and copies its rows to the
     * the matrix object that invoked it only for columns that the flags are
     * true.
     *
     * @@param    A input matrix to be copyied
     * @@param    flags_a 2-dimensional array of boolean variables to tell which
     *                   elements to copy
     */
    public void copyMatrixRows(Matrix A, boolean[] flags_a) 
    {

	col = A.col; 
	row = 0;

	for (int i = 0; i < flags_a.length; i++ )
	{
	    if ( flags_a[i] )
	    {
		row++;
	    }
	}

	int k = 0 ;
	Elem = new double[row][col];
	
	for (int i = 0; i < flags_a.length; i++ )
	{
	    if ( flags_a[i] )
	    {
		for ( int j = 0; j < col; j++ )
		{
		    Elem[k][j] = A.Elem[i][j];
		}
		k++;
	    }
	}
    }
    
    /**
     * This routine makes an identity matrix
     *
     */
    public void identityMatrix() 
    {

	// reset the matrix
	//
	this.resetMatrix();

	// set the diagonal elements to one
	//
	for(int i = 0; i < row; i++ ) 
	{
	    for(int j = 0; j < col; j++) 
	    {
		if (i == j) 
		{
		    Elem[i][j] = 1.0;
		}
	    }
	}
    }

    /**
     * This method solves the linear equation l * l' * x = b, where L is the
     * cholesky decomposition matrix and b is an input vector.
     *
     * @@param  out_vec_a output vector of matrix type
     * @@param  l_a lower triangular matrix cholesky decomposition marix
     * @@param  in_vec_a input vector of matrix type
     *
     * @@return  true
     */
    public boolean choleskySolve(Matrix out_vec_a,
				 Matrix l_a, Matrix in_vec_a) 
    {

	// get the dimensions of input matrix
	//
	int nrows = l_a.getNumRows();
	out_vec_a.initMatrixValue(1, nrows, 0.0, Matrix.FULL);
	
	// System.out.println("l_a nrows: " + nrows);
	// System.out.println("in_vec_a size: " + in_vec_a.getNumColumns());

	// solve L * y = in_vec, result is stored in out_vec
	//
	for (int i = 0; i < nrows; i++) 
	{
	    // declare an accumulator
	    //
	    double sum = in_vec_a.getValue(0, i);

	    // subtract every previous element
	    //
	    for (int j = i - 1; j >= 0; j--) 
	    {
		sum -= (l_a.getValue(i, j) * out_vec_a.getValue(0, j));
	    }
	    
	    // output the result
	    //
	    out_vec_a.setValue(0, i , (sum / l_a.getValue(i,i)));
	}
	
	// solve for the L' * x = y
	//
	for (int i = nrows - 1; i >= 0; i--) 
	{
	    // declare an accumulator
	    //
	    double sum = out_vec_a.getValue(0, i);

	    // subtract every previous element
	    //
	    for (int j = i + 1; j < nrows; j++) 
	    {
		sum -= l_a.getValue(j, i) * out_vec_a.getValue(0, j);
	    }
	    
	    // output the result
	    //
	    out_vec_a.setValue(0, i, (sum / l_a.getValue(i, i)));
	}
	
	// exit gracefully
	//
	return true;
    }

    /**
     * this method constructs the Cholesky decomposition of an input matrix:
     *
     *  W. Press, S. Teukolsky, W. Vetterling, B. Flannery,
     *  Numerical Recipes in C, Second Edition, Cambridge University Press,
     *  New York, New York, USA, pp. 96-97, 1995.
     *
     * upon output, l' * l = this
     * note that this method only operates on symmetric matrices.
     *
     * @@param  l_a output lower triangular matrix
     *
     * @@return true
     */
    public boolean decompositionCholesky(Matrix l_a) {
	
	// condition: non-symmetric
	//
	/*
	if (!this_a.isSymmetric()) {
	    this_a.debug(L"this_a");    
	    return Error::handle(name(), L"decompositionCholesky()", Error::ARG,
				 __FILE__, __LINE__);
	}

	// type: DIAGONAL
	//
	if (this_a.type_d == Integral::DIAGONAL) {

	    // the diagonal elements should be positive or zero
	    //
	    if ((double)this_a.m_d.min() < 0) {
		this_a.debug(L"this_a");      
		return Error::handle(name(), L"decompositionCholesky()", Error::ARG,
				     __FILE__, __LINE__, Error::WARNING);
	    }

	    // assign the input matrix to lower triangular matrix
	    //
	    l_a.assign(this_a);
	    l_a.m_d.sqrt();
	    l_a.changeType(Integral::LOWER_TRIANGULAR);
	}

	// type: FULL, SYMMETRIC, LOWER_TRIANGULAR, UPPER_TRIANGULAR, SPARSE
	//
	else {
	*/

	// get number of rows of current matrix
	//
	int nrows = getNumRows();

	// make a copy of the current matrix
	//
	Matrix a_M = new Matrix();
	a_M.copyMatrix(this);

	// loop over all rows
	//
	for (int i = 0; i < nrows; i++) 
	{

	    // constructs a lower triangular matrix directly iterating
	    // over columns.  remember an upper triangular matrix is
	    // the transpose matrix of a lower triangular matrix in
	    // our storage format.
	    //
	    for (int j = i; j < nrows; j++) 
	    {

		// accumulate the sum
		//
		double sum = a_M.getValue(i, j);
		for (int k = i - 1; k >= 0; k--) 
		{
		    sum -= a_M.getValue(i, k) * a_M.getValue(j, k);
		}

		// handle diagonal elements separately
		//      
		if (i == j) 
		{
		    
		    // if the sum is not greater than zero, the matrix is
		    // not positive definite
		    //
		    if (sum <= 0.0) 
		    {
			// System.out.println("decompositionCholesky()");
			return false;
		    }

		    // assign the value
		    //
		    a_M.setValue(i, i, Math.sqrt(sum));
		}

		// other elements are handled with the normal calculation
		//      
		else 
		{
		    a_M.setValue(j, i, sum / a_M.getValue(i, i));
		}
	    }
	}
	
	// copy the lower triangular part of the result matrix into l_a and
	// set the diagonal elements of l_a matrix to (TIntegral)1
	//
	l_a.initMatrixValue(nrows, nrows, 0.0, LOWER_TRIANGULAR);
	l_a.copyLowerMatrix(a_M);
	
	return true;
	
    }
    
    /**
     * This routine computes the inverse matrix using the gauss jordan method
     * see numerical recipes, the are of scientific computing, second edition,
     * cambridge university press. pages 36 - 41
     *
     * @@param    a   input matrix as 2-d double array
     * @@param    n   number of rows and columns in the input matrix
     * @@param    b   output matrix as 2-d double array
     * @@param    m   number of rows and columns in the output matrix
     */
    public void gaussj(double a[][], int n, double b[][], int m) 
    {

	int indxc[] = null;
	int indxr[] = null;
	int ipiv[] = null;
	int i = 0;
	int icol = 0;
	int irow = 0;
	int j = 0;
	int k = 0;
	int l = 0;
	int ll = 0;
	double big = 0;
	double dum = 0;
	double pivinv = 0;
	double tmp = 0;

	indxc = new int[n];
	indxr = new int[n];
	ipiv = new int[n];

	for (j = 0; j < n; j++) 
	    ipiv[j] = 0;
	for (i = 0; i < n; i++) 
	{
	    big = 0.0;
	    for (j = 0; j < n; j++)
		if (ipiv[j] != 1)
		    for (k = 0; k < n; k++) 
		    {
			if (ipiv[k] == 0) 
			{
			    if (Math.abs(a[j][k]) >= big) 
			    {
				big = Math.abs(a[j][k]);
				irow = j;
				icol = k;
			    }
			} 
			else if (ipiv[k] > 1) 
		       {
			   return;
		       }
		    }
	    ++(ipiv[icol]);
	    if (irow != icol) 
	    {
		for (l = 0; l < n; l++) 
		{
		    tmp = a[irow][l]; 
		    a[irow][l] = a[icol][l]; 
		    a[icol][l] = tmp;
		}
		for (l = 0; l < m; l++) 
		{
		    tmp = b[irow][l]; 
		    b[irow][l] = b[icol][l]; 
		    b[icol][l] = tmp;
		}
	    }
	    indxr[i] = irow;
	    indxc[i] = icol;
	    if (a[icol][icol] == 0.0) 
	    {
		return;
	    }
	    pivinv = 1.0 / a[icol][icol];
	    a[icol][icol] = 1.0;
	    for (l = 0; l < n; l++) 
		a[icol][l] *= pivinv;
	    for (l = 0; l < m; l++) 
		b[icol][l] *= pivinv;
	    for (ll = 0; ll < n; ll++)
		if (ll != icol) 
		{
		    dum = a[ll][icol];
		    a[ll][icol] = 0.0;
		    for (l = 0; l < n; l++) 
			a[ll][l] -= a[icol][l] * dum;
		    for (l = 0; l < m; l++) 
			b[ll][l] -= b[icol][l] * dum;
		}
	}
	for (l = n - 1; l >= 0; l--) 
	{
	    if (indxr[l] != indxc[l]) 
	    {
		for (k = 0; k < n; k++) 
		{
		    tmp = a[k][indxr[l]]; 
		    a[k][indxr[l]] = a[k][indxc[l]]; 
		    a[k][indxc[l]] = tmp;
		}