/*
 *  Draco Econometrics
 *  Copyright (C) 2006-2007 Approximatrix, LLC
 *  http://www.approximatrix.com
 *  
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package com.approximatrix.draco.math;

import org.netlib.lapack.DGECON;
import org.netlib.lapack.DGETRF;
import org.netlib.lapack.DLANGE;
import org.netlib.util.intW;

/*
 * Created on Jan 27, 2005
 *
 * TODO To change the template for this generated file go to
 * Window - Preferences - Java - Code Style - Code Templates
 */

/**
 * @author armstrongj
 *
 * TODO To change the template for this generated type comment go to
 * Window - Preferences - Java - Code Style - Code Templates
 */
public class MatrixSupport {

	static public double[][] Transpose(double[][] matrix, int x, int y) {
		double[][] holding = new double[y][x];
		for(int i=0;i<x;i++) {
			for(int j=0;j<y;j++) {
				holding[j][i] = matrix[i][j];
			}
		}
		return holding;
	}
	
	static public double[][] Identity(int xy) {
		double[][] I = new double[xy][xy];
		for(int i=0;i<xy;i++) {
			for(int j=0;j<xy;j++) {
				if(i==j) {
					I[i][j] = 1.0;
				} else {
					I[i][j] = 0.0;
				}
			}
		}
		return I;
	}
	static public double[][] ScalarTimesMatrix(double alpha, double[][] matrix, int x, int y) {
		double[][] holding = (double[][])matrix.clone();
		for(int i=0;i<x;i++) {
			for(int j=0;j<y;j++) {
				holding[i][j] = alpha*holding[i][j];
			}
		}
		return holding;		
	}
	
	static public double[][] MatrixofZeros(int xy) {
		double[][] Z = new double[xy][xy];
		for(int i=0;i<xy;i++) {
			for(int j=0;j<xy;j++) {
				Z[i][j] = 0.0;
			}
		}
		return Z;
	}
	
	static public double[][] MatrixofZeros(int x, int y) {
		double[][] Z = new double[x][y];
		for(int i=0;i<x;i++) {
			for(int j=0;j<y;j++) {
				Z[i][j] = 0.0;
			}
		}
		return Z;
	}
	
	static public double[][] MatrixWithDiagonal(double[] diagonal)
	{
		double[][] dZ = new double[diagonal.length][diagonal.length];
		for(int i=0;i<diagonal.length;i++) {
			for(int j=0;j<diagonal.length;j++) {
				if(i==j) {
					dZ[i][j] = diagonal[i];
				} else {
					dZ[i][j] = 0.0;
				}
			}
		}
		return dZ;
	}
	
	static public double[][] GetSubMatrix(double[][] matrix,
			                       int startr, int startc,
								   int endr,int endc) {
		double[][] return_val = new double[endr-startr+1][endc-startc+1];
		
		for(int i=startr;i<=endr;i++) {
			for(int j=startc;i<=endc;i++) {
				return_val[i-startr][j-startc] = matrix[i][j]; 
			}
		}
		
		return return_val;
	}
	
	static public double[] GetDiagonal(double[][] matrix) {
		int x = matrix.length;
		int y = matrix[0].length;
		if(x != y) return null;
		
		double[] return_val = new double[x];
		for(int i=0;i<x;i++) {
			return_val[i] = matrix[i][i];
		}
		return return_val;
	}
	
	static public double GetRecipricolCondition(double[][] matrix) {
		
		int matrix_size = matrix.length;
		if(matrix_size != matrix[0].length) return Double.NaN;
		
		int[] ipiv = new int[matrix_size];
		double[] work = new double[matrix_size];
		double[] one_norm_work = new double[matrix_size*4];
		intW result = new intW(0);
		
		// Determines an estimate of the reciprocal of the condition #
		double[][] temp_matrix = CopyMatrix(matrix);
		double one_norm = DLANGE.DLANGE("1",matrix_size,matrix_size,temp_matrix,work);
		
		temp_matrix = CopyMatrix(matrix);
		DGETRF.DGETRF(matrix_size,matrix_size,temp_matrix,ipiv,result);
		
		temp_matrix = CopyMatrix(matrix);
		org.netlib.util.doubleW rcond = new org.netlib.util.doubleW(0.0);
		DGECON.DGECON("1",matrix_size,temp_matrix,one_norm,rcond,one_norm_work,ipiv,result);
		
		if(result.val != 0) {
			System.err.println("Error computing intermediate matrix condition...");
			return 0.0;
		}

		return rcond.val;
	}
	
	static public double[][] MatrixSum(double[][] matrix1, double[][]matrix2, double mat2_scale) {
		double[][] return_val = MatrixofZeros(matrix1.length,matrix1[0].length);
		for(int i=0;i<matrix1.length && i<matrix2.length; i++) {
			for(int j=0;j<matrix1[i].length && j<matrix2[i].length;j++) {
				return_val[i][j] = matrix1[i][j]+mat2_scale*matrix2[i][j];
			}
		}
		return return_val;
	}
	
	static public double[][] MatrixSum(double[][] matrix1, double[][]matrix2) {
		return MatrixSum(matrix1,matrix2,1.0);
	}
	
	static public double[][] MatrixDifference(double[][] matrix1, double[][]matrix2) {
		return MatrixSum(matrix1,matrix2,-1.0);
	}
	
	static public double[][] CopyMatrix(double[][] matrix) {
		double[][] return_val = new double[matrix.length][matrix[0].length];
		
		for(int i=0;i<matrix.length;i++) {
			for(int j=0;j<matrix[i].length;j++) {
				return_val[i][j] = matrix[i][j];
			}
		}
		
		return return_val;
	}
	
	static public double[][] Negate(double[][] matrix) {
		double[][] return_val = new double[matrix.length][matrix[0].length];
		
		for(int i=0;i<matrix.length;i++) {
			for(int j=0;j<matrix[i].length;j++) {
				return_val[i][j] = -1.0*matrix[i][j];
			}
		}
		
		return return_val;
	}
}
