matrix.java

完整的模式识别库

/**
 * 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);
	    }