/*
 * Copyright 2011-2013 JOptimizer
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */
package com.joptimizer.util;

import java.util.Random;

import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.CholeskyDecomposition;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.SingularMatrixException;

import cern.colt.matrix.DoubleFactory2D;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.linalg.Algebra;

/**
 * @author alberto trivellato (alberto.trivellato@gmail.com)
 */
public class Utils {
	
	private static Double RELATIVE_MACHINE_PRECISION = Double.NaN;
	private static Log log = LogFactory.getLog(Utils.class);
	
	public static DoubleMatrix2D randomValuesMatrix(int rows, int cols, double min, double max) {
		return randomValuesMatrix(rows, cols, min, max, null);
	}

	public static DoubleMatrix2D randomValuesMatrix(int rows, int cols, double min, double max, Long seed) {
		Random random = (seed != null) ? new Random(seed) : new Random();

		double[][] matrix = new double[rows][cols];
		for (int i = 0; i < rows; i++) {
			for (int j = 0; j < cols; j++) {
				matrix[i][j] = min + random.nextDouble() * (max - min);
			}
		}
		return DoubleFactory2D.dense.make(matrix);
	}
	
	/**
	 * @TODO: check this!!!
	 * @see "http://mathworld.wolfram.com/PositiveDefiniteMatrix.html"
	 */
	public static DoubleMatrix2D randomValuesPositiveMatrix(int rows, int cols, double min, double max, Long seed) {
		DoubleMatrix2D Q = Utils.randomValuesMatrix(rows, cols, min, max, seed);
		DoubleMatrix2D P = Algebra.DEFAULT.mult(Q, Algebra.DEFAULT.transpose(Q.copy()));
    return Algebra.DEFAULT.mult(P, P);
	}
	
	/**
	 * Residual conditions check after resolution of A.x=b.
	 * 
	 * eps := The relative machine precision
	 * N   := matrix dimension
	 * 
     * Checking the residual of the solution. 
     * Inversion pass if scaled residuals are less than 10:
	 * ||Ax-b||_oo/( (||A||_oo . ||x||_oo + ||b||_oo) . N . eps ) < 10.
	 * 
	 * @param A not-null matrix
	 * @param x not-null vector
	 * @param b not-null vector
	 */
//	public static boolean checkScaledResiduals(DoubleMatrix2D A, DoubleMatrix1D x, DoubleMatrix1D b, Algebra ALG) {
//	  //The relative machine precision
//		double eps = RELATIVE_MACHINE_PRECISION;
//		int N = A.rows();//matrix dimension
//		double residual = -Double.MAX_VALUE;
//		if(Double.compare(ALG.normInfinity(x), 0.)==0 && Double.compare(ALG.normInfinity(b), 0.)==0){
//			return true;
//		}else{
//			residual = ALG.normInfinity(ALG.mult(A, x).assign(b,	Functions.minus)) / 
//	          ((ALG.normInfinity(A)*ALG.normInfinity(x) + ALG.normInfinity(b)) * N * eps);
//			log.debug("scaled residual: " + residual);
//			return residual < 10;
//		}
//	}
	
	/**
	 * The smallest positive (epsilon) such that 1.0 + epsilon != 1.0.
	 * @see http://en.wikipedia.org/wiki/Machine_epsilon#Approximation_using_Java
	 */
	public static final double getDoubleMachineEpsilon() {
		
		if(!Double.isNaN(RELATIVE_MACHINE_PRECISION)){
			return RELATIVE_MACHINE_PRECISION;
		}
		
		synchronized(RELATIVE_MACHINE_PRECISION){
			
			if(!Double.isNaN(RELATIVE_MACHINE_PRECISION)){
				return RELATIVE_MACHINE_PRECISION;
			}
			
			double eps = 1.;
			do {
				eps /= 2.;
			} while ((double) (1. + (eps / 2.)) != 1.);
			
			log.debug("Calculated double machine epsilon: " + eps);
			RELATIVE_MACHINE_PRECISION = eps;
		}
		
		return RELATIVE_MACHINE_PRECISION;
	}
	
	/**
	 * @see http://en.wikipedia.org/wiki/Machine_epsilon#Approximation_using_Java
	 */
//	public static final float calculateFloatMachineEpsilon() {
//		float eps = 1.0f;
//		do {
//			eps /= 2.0f;
//		} while ((float) (1.0 + (eps / 2.0)) != 1.0);
//		log.debug("Calculated float machine epsilon: " + eps);
//		return eps;
//	}
	
	public static RealMatrix squareMatrixInverse(RealMatrix M) throws SingularMatrixException{
		if(!M.isSquare()){
			throw new IllegalArgumentException("Not square matrix!");
		}
		
		// try commons-math cholesky
		try {
			CholeskyDecomposition cd = new CholeskyDecomposition(M);
			return cd.getSolver().getInverse();
		} catch (SingularMatrixException e) {
			throw e;
		} catch (Exception e) {
			log.debug(e.getMessage());
		}
		
		// try joptimizer cholesky
		try {
			CholeskyFactorization cd = new CholeskyFactorization(M.getData());
			double[][] MInv = cd.getInverse();
			return MatrixUtils.createRealMatrix(MInv);
		} catch (Exception e) {
			log.debug(e.getMessage());
		}
		
		// try LU
		try{
			LUDecomposition ld = new LUDecomposition(M);
			return ld.getSolver().getInverse();
		} catch (SingularMatrixException e) {
			throw e;
		} catch (Exception e) {
			log.debug(e.getMessage());
		}

		// try QR
		try {
			QRDecomposition qr = new org.apache.commons.math3.linear.QRDecomposition(M);
			return qr.getSolver().getInverse();
		} catch (SingularMatrixException e) {
			throw e;
		} catch (Exception e) {
			log.debug(e.getMessage());
		}
		
		return null;
	}
	 
	/**
	   * Get the index of the maximum entry.
	   */
	  public static int getMaxIndex(double[] v){
	  	int maxIndex = -1;
	  	double maxValue = -Double.MAX_VALUE;
	  	for(int i=0; imaxValue){
	  			maxIndex = i;
	  			maxValue = v[i]; 
	  		}
	  	}
	  	return maxIndex; 
	  } 
	  
	  
	  
	  /**
	   * Get the index of the minimum entry.
	   */
	  public static int getMinIndex(double[] v){
	  	int minIndex = -1;
	  	double minValue = Double.MAX_VALUE;
	  	for(int i=0; i