package xyz.lejon.utils;

import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;
import java.util.concurrent.ForkJoinPool;
import java.util.concurrent.RecursiveAction;
import java.util.concurrent.ThreadLocalRandom;

import org.ejml.data.DenseMatrix64F;

public class MatrixOps {
	Random rnd = new Random();
	static DecimalFormat mydecimalFormat = new DecimalFormat("00.###E0");
	private static ForkJoinPool pool = new ForkJoinPool();
	private static String DEFAULT_TITLE  = "Vector";
	public static int noDigits = 4;

	public static String arrToStr(int [] arr, String title) {
		String res = "";
		res += title + "[" +  arr.length + "]:[";
		for (int j = 0; j < arr.length; j++) {
			res += arr[j];
			if(j+1<arr.length) 
				res += ", ";
		}
		return res + "]";
	}

	public static String arrToCsv(int [] arr) {
		String res = "";
		for (int j = 0; j < arr.length; j++) {
			res += arr[j];
			if(j+1<arr.length) 
				res += ", ";
		}
		return res;
	}

	public static String arrToCsv(double [] arr) {
		String res = "";
		for (int j = 0; j < arr.length; j++) {
			res += arr[j];
			if(j+1<arr.length) 
				res += ", ";
		}
		return res;
	}

	public static String arrToStr(int [] arr) {
		return arrToStr(arr, "IntVector");
	}

	public static String arrToStr(double [] arr) {
		return arrToStr(arr, DEFAULT_TITLE, Integer.MAX_VALUE);
	}

	public static String arrToStr(double [] arr, int maxLen) {
		return arrToStr(arr, DEFAULT_TITLE, maxLen);
	}
	
	public static String arrToStr(double [] arr, String title) {
		return arrToStr(arr, title, Integer.MAX_VALUE);
	}

	public static String arrToStr(double [] arr, String title, int maxLen) {
		String res = "";
		res += title + "[" +  arr.length + "]:";
		for (int j = 0; j < arr.length && j < maxLen; j++) {
			res += formatDouble(arr[j]) + ", ";
		}
		return res;
	}

	public static String doubleArrayToPrintString(double[][] m) {
		return doubleArrayToPrintString(m, ", ", Integer.MAX_VALUE, m.length, Integer.MAX_VALUE, "\n");
	}

	public static String doubleArrayToPrintString(double[][] m, int maxRows) {
		return doubleArrayToPrintString(m, ", ", maxRows, maxRows, Integer.MAX_VALUE, "\n");
	}

	public static String doubleArrayToPrintString(double[][] m, String colDelimiter) {
		return doubleArrayToPrintString(m, colDelimiter, Integer.MAX_VALUE, -1, Integer.MAX_VALUE, "\n");
	}

	public static String doubleArrayToPrintString(double[][] m, String colDelimiter, int toprowlim) {
		return doubleArrayToPrintString(m, colDelimiter, toprowlim, -1, Integer.MAX_VALUE, "\n");
	}

	public static String doubleArrayToPrintString(double[][] m, int toprowlim, int btmrowlim) {
		return doubleArrayToPrintString(m, ", ", toprowlim, btmrowlim, Integer.MAX_VALUE, "\n");
	}

	public static String doubleArrayToPrintString(double[][] m, int toprowlim, int btmrowlim, int collim) {
		return doubleArrayToPrintString(m, ", ", toprowlim, btmrowlim, collim, "\n");
	}

	public static String doubleArrayToPrintString(double[][] m, String colDelimiter, int toprowlim, int btmrowlim) {
		return doubleArrayToPrintString(m, colDelimiter, toprowlim, btmrowlim, Integer.MAX_VALUE, "\n");
	}

	public static String doubleArrayToPrintString(double[][] m, String colDelimiter, int toprowlim, int btmrowlim, int collim) {
		return doubleArrayToPrintString(m, colDelimiter, toprowlim, btmrowlim, collim, "\n");
	}
	
	public static String doubleArrayToPrintString(double[][] m, String colDelimiter, int toprowlim, int btmrowlim, int collim, String sentenceDelimiter) {
		StringBuffer str = new StringBuffer(m.length * m[0].length);

		str.append("Dim:" + m.length + " x " + m[0].length + "\n");

		int i = 0;
		int lastCol = m[i].length - 1;
		for (; i < m.length && i < toprowlim; i++) {
			String rowPref = i < 1000 ? String.format("%03d", i) : String.format("%04d", i);
			str.append(rowPref+": [");
			for (int j = 0; j < lastCol && j < collim; j++) {
				String formatted = formatDouble(m[i][j]);
				str = str.append(formatted);
				str = str.append(colDelimiter);
			}
			if(lastCol>0)
				str = str.append(formatDouble(m[i][lastCol]));

			if( collim == Integer.MAX_VALUE) { 
				str.append("]");
			} else {
				str.append("...]");
			}
			if (i < m.length - 1) {
				str = str.append(sentenceDelimiter);
			}
		}
		if(btmrowlim<0) return str.toString();
		while(i<(m.length-btmrowlim)) i++;
		if( i < m.length) str.append("\t.\n\t.\n\t.\n");
		for (; i < m.length; i++) {
			String rowPref = i < 1000 ? String.format("%03d", i) : String.format("%04d", i);
			str.append(rowPref+": [");
			for (int j = 0; j < lastCol && j < collim; j++) {
				str = str.append(formatDouble(m[i][j]));
				str = str.append(colDelimiter);
			}
			if(lastCol>0)
				str = str.append(formatDouble(m[i][lastCol]));

			if( collim > m[i].length ) { 
				str.append("]");
			} else {
				str.append(", ...]");
			}
			if (i < m.length - 1) {
				str = str.append(sentenceDelimiter);
			}
		}
		return str.toString();
	}

	public static String formatDouble(double d) {
		if ( d == 0.0 ) return "<0.0>";
		if ( d<0.0001 && d>0 || d > -0.0001 && d < 0) {
			return mydecimalFormat.format(d);
		} else {
			String formatString = "%." + noDigits + "f";
			return String.format(formatString, d);
		}
	}

	public static String doubleArrayToString(double[][] m) {
		return doubleArrayToString(m, ",");
	}

	public static String doubleArrayToString(double[][] m, String colDelimiter) {
		StringBuffer str = new StringBuffer(m.length * m[0].length);
		for (int i = 0; i < m.length; i++) {
			for (int j = 0; j < m[i].length - 1; j++) {
				str = str.append(Double.toString(m[i][j]));
				str = str.append(colDelimiter);
			}
			str = str.append(Double.toString(m[i][m[i].length - 1]));
			str = str.append("\n");
		}
		return str.toString();
	}

	public static double [] rep(double val, int times) {
		double [] res = new double[times];
		for (int i = 0; i < res.length; i++) {
			res[i] = val;
		}
		return res;
	}

	public static double [] asVector(double [][] matrix) {
		boolean isCol = matrix.length != 1;
		int n = matrix.length == 1 ? matrix[0].length : matrix.length;
		if(matrix.length != 1 && matrix[0].length!=1) {
			throw new IllegalArgumentException("Cannot convert non-row or col matrix to vactor! Matrix dim: " 
					+ matrix.length + "x" + matrix[0].length);
		}

		double [] res = new double[n];

		if(isCol) {
			for (int j = 0; j < matrix.length; j++) {				
				res[j] = matrix[j][0]; 
			}
		} else {
			for (int j = 0; j < matrix[0].length; j++) {				
				res[j] = matrix[0][j]; 
			}			
		}

		return res;
	}

	/**
	 * This function returns a new matrix which is centered and scaled, i.e each
	 * the global mean is subtracted from each element in the matrix and divided 
	 * by the global matrix standard deviation   
	 *
	 * @return new matrix which is centered (subtracted mean) and scaled (divided with stddev)
	 */
	public static double [][] centerAndScaleGlobal(double [][] matrix) {
		double [][] res = new double[matrix.length][matrix[0].length]; 
		double mean = mean(matrix);
		double std = stdev(matrix);
		for (int i = 0; i < res.length; i++) {
			for (int j = 0; j < res[i].length; j++) {
				res[i][j] = (matrix[i][j]-mean) / std; 
			}
		}
		
		return res;
	}
	
	/**
	 * This function returns a new matrix which is centered and scaled, i.e each
	 * the column mean is subtracted from each column element in the matrix and 
	 * divided by the respective column matrix standard deviation 
	 *   
	 * @return new matrix which is centered (subtracted mean) and scaled (divided with stddev)
	 */
	public static double [][] centerAndScale(double [][] matrix) {
		double [][] res = new double[matrix.length][matrix[0].length]; 
		double [] means = colMeans(matrix);
		for (int i = 0; i < res.length; i++) {
			for (int j = 0; j < res[i].length; j++) {
				res[i][j] = (matrix[i][j]-means[j]); 
			}
		}
		
		double [] std = colStddev(res);
		for (int i = 0; i < res.length; i++) {
			for (int j = 0; j < res[i].length; j++) {
				res[i][j] = res[i][j] / (std[j] == 0 ? 1 : std[j]); 
			}
		}

		return res;
	}
	
	public static double [][] centerAndScaleSametime(double [][] matrix) {
		double [][] res = new double[matrix.length][matrix[0].length]; 
		double [] means = colMeans(matrix);
		double [] std = colStddev(matrix);
		for (int i = 0; i < res.length; i++) {
			for (int j = 0; j < res[i].length; j++) {
				res[i][j] = (matrix[i][j]-means[j]) / (std[j] == 0 ? 1 : std[j]); 
			}
		}

		return res;
	}
	
	/**
	 * This function returns adds a small amount of noise to each column  
	 * 
	 * @return new matrix with added noise
	 */
	public static double [][] addNoise(double [][] matrix) {
		double [][] res = new double[matrix.length][matrix[0].length]; 
		double [] std = colStddev(matrix);
		for (int i = 0; i < res.length; i++) {
			for (int j = 0; j < res[i].length; j++) {
				double noise = rnorm(0, std[j] == 0.0 ? 0.00001 : std[j]/5);
				res[i][j] = matrix[i][j] + noise; 
			}
		}

		return res;
	}

	// Unit Tested
	double[][] transpose(double[][] matrix) {
		return transpose(matrix, 1000);
	}

	/**
	 * Returns a new matrix which is the transpose of input matrix
	 * @param matrix
	 * @return
	 */
	public static double[][] transposeSerial(double[][] matrix) {
		int rows = matrix.length;
		int cols = matrix[0].length;
		double[][] transpose = new double[cols][rows];
		for (int col = 0; col < cols; col++)
			for (int row = 0; row < rows; row++)
				transpose[col][row] = matrix[row][col];
		return transpose;
	}
	
	public static long[][] transposeSerial(long[][] matrix) {
		int rows = matrix.length;
		int cols = matrix[0].length;
		long[][] transpose = new long[cols][rows];
		for (int col = 0; col < cols; col++)
			for (int row = 0; row < rows; row++)
				transpose[col][row] = matrix[row][col];
		return transpose;
	}

	// Unit Tested
	/**
	 * Returns a new matrix which is the transpose of input matrix
	 * @param matrix
	 * @return
	 */
	public double[][] transpose(double[][] matrix, int ll) {
		int cols = matrix[0].length;
		int rows = matrix.length;
		double[][] transpose = new double[cols][rows];
		if(rows < 100 ) {
			for (int i = 0; i < cols; i++)
				for (int j = 0; j < rows; j++)
					transpose[i][j] = matrix[j][i];
		} else {
			MatrixTransposer process = new MatrixTransposer(matrix, transpose,0,rows,ll);                
			pool.invoke(process);
		}
		return transpose;
	}

	class MatrixTransposer extends RecursiveAction {
		private static final long serialVersionUID = 1L;
		double [][] orig;
		double [][] transpose;
		int startRow = -1;
		int endRow = -1;
		int limit = 1000;

		public MatrixTransposer(double [][] orig, double [][] transpose, int startRow, int endRow, int ll) {
			this.limit = ll;
			this.orig = orig;
			this.transpose = transpose;
			this.startRow = startRow;
			this.endRow = endRow;
		}

		public MatrixTransposer(double [][] orig, double [][] transpose, int startRow, int endRow) {
			this.orig = orig;
			this.transpose = transpose;
			this.startRow = startRow;
			this.endRow = endRow;
		}

		@Override
		protected void compute() {
			try {
				if ( (endRow-startRow) <= limit ) {
					int cols = orig[0].length;
					for (int i = 0; i < cols; i++) {
						for (int j = startRow; j < endRow; j++) {
							transpose[i][j] = orig[j][i];
						}
					}
				}
				else {
					int range = (endRow-startRow);
					int startRow1 = startRow;
					int endRow1 = startRow + (range / 2);
					int startRow2 = endRow1;
					int endRow2 = endRow;
					invokeAll(new MatrixTransposer(orig, transpose, startRow1, endRow1, limit),
							new MatrixTransposer(orig, transpose, startRow2, endRow2, limit));
				}
			}
			catch ( Exception e ) {
				e.printStackTrace();
			}
		}
	}

	void plainTranspose(double[][] matrix, double[][] transpose, int startRow, int endRow) {
		int cols = matrix[0].length;
		for (int i = 0; i < cols; i++)
			for (int j = startRow; j < endRow; j++)
				transpose[i][j] = matrix[j][i];
	}

	// Unit Tested
	/**
	 * Returns a new matrix with values exponentiated
	 * @param matrix
	 * @return new matrix with values exponentiated
	 */
	public static double [][] exp(double [][] m1) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				matrix[i][j] = Math.exp(m1[i][j]);
			}
		}
		return matrix;
	}

	// Unit Tested
	/**
	 * Returns new vetcor with vals sqrt'ed
	 * @param vector
	 * @return new vector with values sqrt'ed
	 */
	public static double [] sqrt(double [] v1) {
		double [] vector = new double[v1.length];
		for (int i = 0; i < vector.length; i++) {
			vector[i] = Math.sqrt(v1[i]);
		}
		return vector;
	}

	// Unit Tested
	/**
	 * @param vector
	 * @return mean of values in vector
	 */
	public static double mean(double [] vector) {
		double sum = 0.0;
		for (int i = 0; i < vector.length; i++) {
			sum +=vector[i];
		}
		return sum/vector.length;
	}

	// Unit Tested
	/**
	 * Returns a new matrix with values that are the log of the input matrix
	 * @param matrix
	 * @return  same matrix with values log'ed
	 */
	public static double [][] log(double [][] m1) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				matrix[i][j] = Math.log(m1[i][j]);
			}
		}
		return matrix;
	}

	/**
	 * Returns a new matrix with values that are taken to the power of the input matrix
	 * @param matrix
	 * @param power
	 * @return  same matrix with values pow'ed
	 */
	public static double [][] pow(double [][] m1, double power) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				matrix[i][j] = Math.pow(m1[i][j], power);
			}
		}
		return matrix;
	}

	/**
	 * Returns a new matrix with values that are taken to the power of the input matrix
	 * @param matrix
	 * @param power
	 * @return  same matrix with values pow'ed
	 */
	public static double [] pow(double [] m1, double power) {
		double[] matrix = new double[m1.length];
		for (int i = 0; i < matrix.length; i++) {
			matrix[i] = Math.pow(m1[i], power);
		}
		return matrix;
	}

	/**
	 * Returns a new matrix with values that are the log of the input matrix
	 * @param matrix
	 * @param infAsZero treat +- Infinity as zero, i.e replaces Infinity with 0.0
	 * if set to true
	 * @return  same matrix with values log'ed
	 */
	public static double [][] log(double [][] m1, boolean infAsZero) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				matrix[i][j] = Math.log(m1[i][j]);
				if(infAsZero && Double.isInfinite(matrix[i][j]))
					matrix[i][j] = 0.0;
			}
		}
		return matrix;
	}

	// Unit Tested
	/**
	 * @param matrix
	 * @return scalar inverse of matrix
	 */
	public static double [][] scalarInverse(double [][] m1) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				matrix[i][j] = 1/m1[i][j];
			}
		}
		return matrix;
	}

	// Unit Tested
	/**
	 * @param vector
	 * @return scalar inverse of vector
	 */
	public static double [] scalarInverse(double [] v1) {
		double [] vector = new double[v1.length];
		for (int i = 0; i < vector.length; i++) {
			vector[i] = 1/v1[i];
		}
		return vector;
	}

	/**
	 * @param m
	 * @param n
	 * @return new 2D matrix with normal random values with mean 0 and std. dev 1
	 */
	public static double[][] rnorm(int m, int n) {
		double[][] array = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < array[i].length; j++) {				
				array[i][j] = rnorm(0.0,1.0);
			}
		}
		return array;
	}

	public static double [] rnorm(int n, double [] mus, double [] sigmas) {
		double [] res = new double[n];
		for (int i = 0; i < res.length; i++) {
			res[i] = mus[i] + (ThreadLocalRandom.current().nextGaussian() * sigmas[i]);
		}
		return res; 
	}

	public static double [] rnorm(int n, double mu, double [] sigmas) {
		double [] res = new double[n];
		for (int i = 0; i < res.length; i++) {
			res[i] = mu + (ThreadLocalRandom.current().nextGaussian() * sigmas[i]);
		}
		return res; 
	}

	public static double rnorm() {
		return ThreadLocalRandom.current().nextGaussian();
	}

	public static double rnorm(double mu, double sigma) {
		return mu + (ThreadLocalRandom.current().nextGaussian() * sigma);
	}

	// Unit Tested
	/**
	 * Returns a new matrix of booleans where true is set if the values to the two matrices are
	 * the same at that index
	 * @param matrix1
	 * @param matrix2
	 * @return new matrix with booelans with values matrix1[i,j] == matrix2[i,j]
	 */
	public static boolean [][] equal(double [][] matrix1, double [][] matrix2) {
		boolean [][] equals = new boolean[matrix1.length][matrix1[0].length];
		if( matrix1.length != matrix2.length) {
			throw new IllegalArgumentException("Dimensions does not match");
		}
		if( matrix1[0].length != matrix2[0].length) {
			throw new IllegalArgumentException("Dimensions does not match");
		}
		for (int i = 0; i < matrix1.length; i++) {
			for (int j = 0; j < matrix1[0].length; j++) {
				equals[i][j] = Double.compare(matrix1[i][j], matrix2[i][j]) == 0;
			}
		}
		return equals;
	}

	/**
	 * Returns a new matrix of booleans where true is set if the values to the two matrices are
	 * the same at that index
	 * @param matrix1
	 * @param matrix2
	 * @return new matrix with booelans with values matrix1[i,j] == matrix2[i,j]
	 */
	public static boolean [][] equal(boolean [][] matrix1, boolean [][] matrix2) {
		boolean [][] equals = new boolean[matrix1.length][matrix1[0].length];
		if( matrix1.length != matrix2.length) {
			throw new IllegalArgumentException("Dimensions does not match");
		}
		if( matrix1[0].length != matrix2[0].length) {
			throw new IllegalArgumentException("Dimensions does not match");
		}
		for (int i = 0; i < matrix1.length; i++) {
			for (int j = 0; j < matrix1[0].length; j++) {
				equals[i][j] = (matrix1[i][j] == matrix2[i][j]);
			}
		}
		return equals;
	}

	/**
	 * Returns a new matrix of booleans where true is set if the value in the matrix is
	 * bigger than value
	 * @param matrix
	 * @param value
	 * @return new matrix with booelans with values matrix1[i,j] == matrix2[i,j]
	 */
	public static boolean [][] biggerThan(double [][] matrix, double value) {
		boolean [][] equals = new boolean[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				equals[i][j] = Double.compare(matrix[i][j], value) == 1;
			}
		}
		return equals;
	}

	/**
	 * @param booleans
	 * @return new matrix with booleans which are the negations of the input
	 */
	public static boolean [][] negate(boolean [][] booleans) {
		boolean [][] negates = new boolean[booleans.length][booleans[0].length];
		for (int i = 0; i < booleans.length; i++) {
			for (int j = 0; j < booleans[0].length; j++) {
				negates[i][j] = !booleans[i][j];
			}
		}
		return negates;
	}

	/**
	 * @param booleans
	 * @return
	 */
	public static double [][] abs(boolean [][] booleans) {
		double [][] absolutes = new double[booleans.length][booleans[0].length];
		for (int i = 0; i < booleans.length; i++) {
			for (int j = 0; j < booleans[0].length; j++) {
				absolutes[i][j] = booleans[i][j] ? 1 : 0;
			}
		}
		return absolutes;
	}

	/**
	 * @param vals
	 * @return
	 */
	public static double [][] abs(double [][] vals) {
		double [][] absolutes = new double[vals.length][vals[0].length];
		for (int i = 0; i < vals.length; i++) {
			for (int j = 0; j < vals[0].length; j++) {
				absolutes[i][j] = Math.abs(vals[i][j]);
			}
		}
		return absolutes;
	}

	/**
	 * @param absolutes
	 * @return
	 */
	public static double [] abs(double [] vals) {
		double [] absolutes = new double[vals.length];
		for (int i = 0; i < vals.length; i++) {
			absolutes[i] = Math.abs(vals[i]);
		}
		return absolutes;
	}

	public static double [][] sign(double [][] matrix) {
		double [][] signs = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				signs[i][j] = matrix[i][j] >= 0 ? 1 : -1;
			}
		}
		return signs;
	}

	public static double mean(double [][] matrix) {
		return mean(matrix,2)[0][0];
	}
	
	// Unit Tested
	public static double [][] mean(double [][] matrix, int axis) {
		// Axis = 0 => sum columns
		// Axis = 1 => sum rows
		// Axis = 2 => global (returns a 1 element array with the result)
		double [][] result;
		if( axis == 0) {
			result = new double[1][matrix[0].length];
			for (int j = 0; j < matrix[0].length; j++) {
				double colsum = 0.0;
				for (int i = 0; i < matrix.length; i++) {
					colsum += matrix[i][j];
				}
				result[0][j] = colsum / matrix.length;
			}
		}   else if (axis == 1) {
			result = new double[matrix.length][1];
			for (int i = 0; i < matrix.length; i++) {
				double rowsum = 0.0;
				for (int j = 0; j < matrix[0].length; j++) {
					rowsum += matrix[i][j];
				}
				result[i][0] = rowsum / matrix[0].length;
			}
		}   else if (axis == 2) {
			result = new double[1][1];
			for (int j = 0; j < matrix[0].length; j++) {
				for (int i = 0; i < matrix.length; i++) {
					result[0][0] += matrix[i][j];
				}
			}
			result[0][0] /=  (matrix[0].length *  matrix.length);
		}else {
			throw  new IllegalArgumentException("Axes other than 0,1,2 is unsupported");
		}
		return result;
	}

	// Unit Tested
	// Should be called dim-sum! :)
	public static double [][] sum(double [][] matrix, int axis) {
		// Axis = 0 => sum columns
		// Axis = 1 => sum rows
		double [][] result;
		if( axis == 0) {
			result = new double[1][matrix[0].length];
			for (int j = 0; j < matrix[0].length; j++) {
				double rowsum = 0.0;
				for (int i = 0; i < matrix.length; i++) {
					rowsum += matrix[i][j];
				}
				result[0][j] = rowsum;
			}
		}   else if (axis == 1) {
			result = new double[matrix.length][1];
			for (int i = 0; i < matrix.length; i++) {
				double colsum = 0.0;
				for (int j = 0; j < matrix[0].length; j++) {
					colsum += matrix[i][j];
				}
				result[i][0] = colsum;
			}
		}   else {
			throw  new IllegalArgumentException("Axes other than 0,1 is unsupported");
		}
		return result;
	}
	
	// Unit Tested
	// Should be called dim-sum! :)
	public static int [][] sum(int [][] matrix, int axis) {
		// Axis = 0 => sum columns
		// Axis = 1 => sum rows
		int [][] result;
		if( axis == 0) {
			result = new int[1][matrix[0].length];
			for (int j = 0; j < matrix[0].length; j++) {
				int rowsum = 0;
				for (int i = 0; i < matrix.length; i++) {
					rowsum += matrix[i][j];
				}
				result[0][j] = rowsum;
			}
		}   else if (axis == 1) {
			result = new int[matrix.length][1];
			for (int i = 0; i < matrix.length; i++) {
				int colsum = 0;
				for (int j = 0; j < matrix[0].length; j++) {
					colsum += matrix[i][j];
				}
				result[i][0] = colsum;
			}
		}   else {
			throw  new IllegalArgumentException("Axes other than 0,1 is unsupported");
		}
		return result;
	}
	
	public static int [] colSum(int [][] matrix) {
		int [] result = new int[matrix[0].length];
		for (int j = 0; j < matrix[0].length; j++) {
			int rowsum = 0;
			for (int i = 0; i < matrix.length; i++) {
				rowsum += matrix[i][j];
			}
			result[j] = rowsum;
		}
		return result;
	}
	
	public static double [] colSum(double [][] matrix) {
		double [] result = new double[matrix[0].length];
		for (int j = 0; j < matrix[0].length; j++) {
			int rowsum = 0;
			for (int i = 0; i < matrix.length; i++) {
				rowsum += matrix[i][j];
			}
			result[j] = rowsum;
		}
		return result;
	}

	// Unit Tested
	/**
	 * Returns a new matrix which is the transpose of input matrix
	 * @param matrix
	 * @return
	 */
	public double sumPar(double[][] matrix) {
		int ll = 100;
		int cols = matrix[0].length;
		int rows = matrix.length;
		double [] sums = new double[rows];
		if(rows < ll ) {
			for (int row = 0; row < rows; row++)
				for (int col = 0; col < cols; col++)
					sums[row] += matrix[row][col];
		} else {
			MatrixSummer process = new MatrixSummer(matrix, sums, 0, rows, ll);                
			pool.invoke(process);
		}
		double sum = 0.0;
		for (int i = 0; i < sums.length; i++) {
			sum += sums[i];
		}
		return sum;
	}

	class MatrixSummer extends RecursiveAction {
		private static final long serialVersionUID = 1L;
		double [][] orig;
		double [] sums;
		int startRow = -1;
		int endRow = -1;
		int limit = 1000;

		public MatrixSummer(double [][] orig, double [] sums, int startRow, int endRow, int ll) {
			this.limit = ll;
			this.orig = orig;
			this.sums = sums;
			this.startRow = startRow;
			this.endRow = endRow;
		}

		public MatrixSummer(double [][] orig, double [] transpose, int startRow, int endRow) {
			this.orig = orig;
			this.sums = transpose;
			this.startRow = startRow;
			this.endRow = endRow;
		}

		@Override
		protected void compute() {
			try {
				if ( (endRow-startRow) <= limit ) {
					int cols = orig[0].length;
					for (int row = startRow; row < endRow; row++) {
						for (int i = 0; i < cols; i++) {
							sums[row] += orig[row][i];
						}
					}
				}
				else {
					int range = (endRow-startRow);
					int startRow1 = startRow;
					int endRow1 = startRow + (range / 2);
					int startRow2 = endRow1;
					int endRow2 = endRow;
					invokeAll(new MatrixSummer(orig, sums, startRow1, endRow1, limit),
							new MatrixSummer(orig, sums, startRow2, endRow2, limit));
				}
			}
			catch ( Exception e ) {
				e.printStackTrace();
			}
		}
	}


	// Unit Tested
	/**
	 * @param matrix
	 * @return sum of all values in the matrix
	 */
	public static double sum(double [][] matrix) {
		double sum = 0.0;
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				sum+=matrix[i][j];
			}
		}
		return sum;
	}

	// Unit Tested
	/**
	 * @param matrix
	 * @return sum of all values in the matrix
	 */
	public static int sum(int [][] matrix) {
		int sum = 0;
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				sum+=matrix[i][j];
			}
		}
		return sum;
	}
	
	// Unit Tested
	/**
	 * @param vector
	 * @return sum of all values in the matrix
	 */
	public static int sum(int [] vector) {
		int sum = 0;
		for (int i = 0; i < vector.length; i++) {
			sum+=vector[i];
		}
		return sum;
	}

	public static double sum(double [] vector) {
		double res = 0.0;
		for (int i = 0; i < vector.length; i++) {
			res += vector[i];
		}
		return res;
	}

	/**
	 * Return a new matrix with the max value of either the value in the matrix 
	 * or maxval otherwise 
	 * @param matrix
	 * @param maxval
	 * @return
	 */
	public static double [][] maximum(double [][] matrix, double maxval) {
		double [][] maxed = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				maxed[i][j] = matrix[i][j] > maxval ? matrix[i][j] : maxval;
			}
		}
		return maxed;
	}

	// Unit Tested
	/**
	 * All values in matrix that is less than <code>lessthan</code> is assigned
	 * the value <code>assign</code>
	 * @param matrix
	 * @param lessthan
	 * @param assign
	 * @return
	 */
	public static void assignAllLessThan(double[][] matrix, double lessthan, double assign) {
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				if( matrix[i][j] < lessthan) {
					matrix[i][j] = assign;
				}
			}
		}
	}

	// Unit Tested
	/**
	 * @param matrix
	 * @return a new matrix with the values of matrix squared
	 */
	public static double [][] square(double [][] matrix) {
		return scalarPow(matrix,2);
	}

	/** 
	 * Replaces NaN's with repl
	 * @param matrix
	 * @param repl
	 * @return
	 */
	public static double [][] replaceNaN(double [][] matrix, double repl) {
		double [][] result = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				if(Double.isNaN(matrix[i][j])) {
					result[i][j] = repl;
				} else {
					result[i][j] = matrix[i][j];
				}
			}
		}
		return result;
	}

	/** 
	 * Replaces Infinity's with repl
	 * @param matrix
	 * @param repl
	 * @return
	 */
	public static double [][] replaceInf(double [][] matrix, double repl) {
		double [][] result = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				if(Double.isInfinite(matrix[i][j])) {
					result[i][j] = repl;
				} else {
					result[i][j] = matrix[i][j];
				}
			}
		}
		return result;
	}

	
	public static double [][] scalarPow(double [][] matrix, double power) {
		double [][] result = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				result[i][j] += Math.pow(matrix[i][j],power);
			}
		}
		return result;
	}

	public static double [][] addColumnVector(double [][] matrix, double [][] colvector) {
		double [][] result = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				result[i][j] = matrix[i][j] + colvector[i][0];
			}
		}
		return result;
	}

	public static double [][] addRowVector(double [][] matrix, double [][] rowvector) {
		double [][] result = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				result[i][j] = matrix[i][j] + rowvector[0][j];
			}
		}
		return result;
	}

	public static double [][] fillWithRowOld(double [][] matrix, int row) {
		double [][] result = new double[matrix.length][matrix[0].length];
		for (int i = 0; i < matrix.length; i++) {
			for (int j = 0; j < matrix[0].length; j++) {
				result[i][j] = matrix[row][j];
			}
		}
		return result;
	}

	public static double [][] fillWithRow(double [][] matrix, int row) {
		int rows = matrix.length;
		int cols = matrix[0].length;
		double [][] result = new double[rows][cols];
		for (int i = 0; i < rows; i++) {
			System.arraycopy(matrix[row], 0, result[i], 0, cols);
		}
		return result;
	}

	// Unit Tested
	public static double [][] tile(double [][] matrix, int rowtimes, int coltimes) {
		double [][] result = new double[matrix.length*rowtimes][matrix[0].length*coltimes];
		for (int i = 0, resultrow = 0; i < rowtimes; i++) {
			for (int j = 0; j < matrix.length; j++) {
				for (int k = 0, resultcol = 0; k < coltimes; k++) {
					for (int l = 0; l < matrix[0].length; l++) {
						result[resultrow][resultcol++] = matrix[j][l];
					}
				}
				resultrow++;
			}
		}

		return result;
	}

	public static double[][] normalize(double[][] x, double[] meanX, double[] stdevX) {
		double[][] y = new double[x.length][x[0].length];
		for (int i = 0; i < y.length; i++)
			for (int j = 0; j < y[i].length; j++)
				y[i][j] = (x[i][j] - meanX[j]) / stdevX[j];
		return y;
	}
	
	public static int [] range(int n) {
		int [] result = new int[n];
		for (int i = 0; i < n; i++) {
			result[i] = i;
		}
		return result;
	}

	public static int [] range(int a, int b) {
		if( b < a ) {
			throw new IllegalArgumentException("b has to be larger than a");
		}
		int val = a;
		int [] result = new int[b-a];
		for (int i = 0; i < (b-a); i++) {
			result[i] = val++;
		}
		return result;
	}

	// Unit Tested
	public static int [] concatenate(int [] v1,int [] v2) {
		int [] result = new int[v1.length+v2.length];
		int index = 0;
		for (int i = 0; i < v1.length; i++, index++) {
			result[index] = v1[index];
		}
		for (int i = 0; i < v2.length; i++, index++) {
			result[index] = v2[i];
		}
		return result;
	}
	
	// Unit Tested
		public static double [] concatenateOLD(double [] v1,double [] v2) {
			double [] result = new double[v1.length+v2.length];
			int index = 0;
			for (int i = 0; i < v1.length; i++, index++) {
				result[index] = v1[index];
			}
			for (int i = 0; i < v2.length; i++, index++) {
				result[index] = v2[i];
			}
			return result;
		}
		
		public static double [] concatenate(double [] v1,double [] v2) {
			double [] result = new double[v1.length+v2.length];
			System.arraycopy( v1, 0, result, 0,  v1.length);
			System.arraycopy( v2, 0, result, v1.length, v2.length);			
			return result;
		}
		

	// Unit Tested
	public static double [][] concatenateOLD(double [][] m1,double[][] m2) {
		if(m1.length!=m2.length) throw new IllegalArgumentException("m1 and m2 must have the same number of rows:" + m1.length + " != " + m2.length);
		double [][] result = new double[m1.length][m1[0].length+m2[0].length];
		int resCol = 0;
		for (int i = 0; i < m1.length; i++) {
			resCol = 0;
			for (int j = 0; j < m1[i].length; j++) {
				result[i][resCol++] = m1[i][j];
			}
			for (int j = 0; j < m2[i].length; j++) {
				result[i][resCol++] = m2[i][j];
			}
		}	
		return result;
	}
	
	// Unit Tested
	// New version, hopefully faster than above
	public static double [][] concatenate(double [][] m1,double[][] m2) {
		if(m1.length!=m2.length) throw new IllegalArgumentException("m1 and m2 must have the same number of rows:" + m1.length + " != " + m2.length);
		double [][] result = new double[m1.length][m1[0].length+m2[0].length];
		for (int i = 0; i < m1.length; i++) {
			System.arraycopy( m1[i], 0, result[i], 0,  m1[i].length);
			System.arraycopy( m2[i], 0, result[i], m1[i].length,  m2[i].length);
		}	
		return result;
	}

	// Unit Tested
	public static double [][] concatenate(double [][] m1,double[] v2) {
		if(m1.length!=v2.length) throw new IllegalArgumentException("m1 and v2 must have the same number of rows:" + m1.length + " != " + v2.length);
		double [][] result = new double[m1.length][m1[0].length+1];
		int resCol = 0;
		for (int i = 0; i < m1.length; i++) {
			resCol = 0;
			for (int j = 0; j < m1[i].length; j++) {
				result[i][resCol++] = m1[i][j];
			}
			result[i][resCol++] = v2[i];
		}	
		return result;
	}

	public double [][] scalarMultiply(double [][] m1,double [][] m2) {
		return parScalarMultiply(m1, m2);
	}
	
	// Unit Tested
	public static double [][] sMultiply(double [][] v1,double [][] v2) {
		if( v1.length != v2.length || v1[0].length != v2[0].length ) {
			throw new IllegalArgumentException("a and b has to be of equal dimensions");
		}
		double [][] result = new double[v1.length][v1[0].length];
		for (int i = 0; i < v1.length; i++) {
			for (int j = 0; j < v1[0].length; j++) {
				result[i][j] = v1[i][j] * v2[i][j];
			}
		}
		return result;
	}
	
	public double[][] parScalarMultiply(double [][] m1,double [][] m2) {
		int ll = 600;
		double [][] result = new double[m1.length][m1[0].length];

		MatrixOperator process = new MatrixOperator(m1,m2,result, multiplyop, 0, m1.length,ll);                
		pool.invoke(process);
		return result;
	}

	public double[][] parScalarMinus(double [][] m1,double [][] m2) {
		int ll = 600;
		double [][] result = new double[m1.length][m1[0].length];
		
		MatrixOperator process = new MatrixOperator(m1,m2,result, minusop, 0, m1.length,ll);                
		pool.invoke(process);
		return result;
	}

	public interface MatrixOp {
		double compute(double op1, double op2);
	}

	MatrixOp multiplyop = new MatrixOp() {
		public double compute(double f1, double f2) {
			return f1 * f2;
		}
	};

	MatrixOp minusop = new MatrixOp() {
		public double compute(double f1, double f2) {
			return f1 - f2;
		}
	};

	class MatrixOperator extends RecursiveAction {
		final static long serialVersionUID = 1L;
		double [][] matrix1;
		double [][] matrix2;
		double [][] resultMatrix;
		int startRow = -1;
		int endRow = -1;
		int limit = 1000;
		MatrixOp op;

		public MatrixOperator(double [][] matrix1, double [][] matrix2, double [][] resultMatrix, 
				MatrixOp op, int startRow, int endRow, int ll) {
			this.op = op;
			this.limit = ll;
			this.matrix1 = matrix1;
			this.matrix2 = matrix2;
			this.resultMatrix = resultMatrix;
			this.startRow = startRow;
			this.endRow = endRow;
		}

		@Override
		protected void compute() {
			try {
				if ( (endRow-startRow) <= limit ) {
					int cols = matrix1[0].length;
					for (int i = startRow; i < endRow; i++) {
						for (int j = 0; j < cols; j++) {
							resultMatrix[i][j] = op.compute(matrix1[i][j], matrix2[i][j]);
						}
					}
				}
				else {
					int range = (endRow-startRow);
					int startRow1 = startRow;
					int endRow1 = startRow + (range / 2);
					int startRow2 = endRow1;
					int endRow2 = endRow;
					invokeAll(new MatrixOperator(matrix1, matrix2, resultMatrix, op, startRow1, endRow1, limit),
							new MatrixOperator(matrix1, matrix2, resultMatrix, op, startRow2, endRow2, limit));
				}
			}
			catch ( Exception e ) {
				e.printStackTrace();
			}
		}
	}


	public static void assignAtIndex(double[][] num, int[] range, int[] range1, double value) {
		for (int j = 0; j < range.length; j++) {
			num[range[j]][range1[j]] = value;
		}
	}

	public static double [][] getValuesFromRow(double[][] matrix, int row, int[] indicies) {
		double [][] values = new double[1][indicies.length];
		for (int j = 0; j < indicies.length; j++) {
			values[0][j] = matrix[row][indicies[j]];
		}
		return values;
	}

	public static void assignValuesToRow(double[][] matrix, int row, int[] indicies, double [] values) {
		if( indicies.length != values.length ) {
			throw new IllegalArgumentException("Length of indicies and values have to be equal");
		}
		for (int j = 0; j < indicies.length; j++) {
			matrix[row][indicies[j]] = values[j];
		}
	}

	public static double stdev(double [][] matrix) {
		double m = mean(matrix);

        double total = 0;

        final int N = matrix.length * matrix[0].length;

        for( int i = 0; i < matrix.length; i++ ) {
        	for (int j = 0; j < matrix[i].length; j++) {				
        		double x = matrix[i][j];
        		total += (x - m)*(x - m);
			}
        }

        return Math.sqrt(total / (N-1));	
	}
	
	public static double[] colStddev(double[][] v) {
		double[] var = variance(v);
		for (int i = 0; i < var.length; i++)
			var[i] = Math.sqrt(var[i]);
		return var;
	}

	public static double[] variance(double[][] v) {
		int m = v.length;
		int n = v[0].length;
		double[] var = new double[n];
		int degrees = (m - 1);
		double c;
		double s;
		for (int j = 0; j < n; j++) {
			c = 0;
			s = 0;
			for (int k = 0; k < m; k++)
				s += v[k][j];
			s = s / m;
			for (int k = 0; k < m; k++)
				c += (v[k][j] - s) * (v[k][j] - s);
			var[j] = c / degrees;
		}
		return var;
	}

	/**
	 * @param matrix
	 * @return a new vector with the column means of matrix
	 */
	public static double[] colMeans(double[][] matrix) {
		int rows = matrix.length;
		int cols = matrix[0].length;
		double[] mean = new double[cols];
		for (int i = 0; i < rows; i++)
			for (int j = 0; j < cols; j++)
				mean[j] += matrix[i][j];
		for (int j = 0; j < cols; j++)
			mean[j] /= (double) rows;
		return mean;
	}
	
	public static double[] colMeans(long[][] matrix) {
		int rows = matrix.length;
		int cols = matrix[0].length;
		double[] mean = new double[cols];
		for (int i = 0; i < rows; i++)
			for (int j = 0; j < cols; j++)
				mean[j] += matrix[i][j];
		for (int j = 0; j < cols; j++)
			mean[j] /= (double) rows;
		return mean;
	}

	public static double[][] copyRows(double[][] input, int... indices) {
		double[][] matrix = new double[indices.length][input[0].length];
		for (int i = 0; i < indices.length; i++)
			System.arraycopy(input[indices[i]], 0, matrix[i], 0, input[indices[i]].length);
		return matrix;
	}

	public static double[][] copyCols(double[][] input, int... indices) {
		double[][] matrix = new double[indices.length][input.length];
		for (int i = 0; i < indices.length; i++)
			for (int j = 0; j < input.length; j++) {
				matrix[i][j] = input[j][indices[i]];
			}
		return matrix;
	}
	
	public static double [][] clone(double [][] arr1) {
		double [][] result = new double [arr1.length][];
		
		for (int i = 0; i < arr1.length; i++) {
			result[i] = arr1[i].clone();
		}
		return result;
	}

	public static double[][] fillMatrix(int rows, int cols, double fillvalue) {
		double[][] matrix = new double[rows][cols];
		for (int i = 0; i < matrix.length; i++)
			for (int j = 0; j < matrix[i].length; j++)
				matrix[i][j] = fillvalue;
		return matrix;
	}

	public static double[][] plus(double[][] m1, double[][] m2) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < m1.length; i++)
			for (int j = 0; j < m1[0].length; j++)
				matrix[i][j] = m1[i][j] + m2[i][j];
		return matrix;
	}

	// Unit Tested
	public static double[][] scalarPlus(double[][] m1, double m2) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < m1.length; i++)
			for (int j = 0; j < m1[0].length; j++)
				matrix[i][j] = m1[i][j] + m2;
		return matrix;
	}

	// Unit Tested
	public static double [] scalarPlus(double[] m1, double m2) {
		double[] matrix = new double[m1.length];
		for (int i = 0; i < m1.length; i++)
			matrix[i] = m1[i] + m2;
		return matrix;
	}

	public double[][] minus(double[][] m1, double[][] m2) {
		return parScalarMinus(m1, m2);
	}

	// Unit Tested
	public static double[][] sMinus(double[][] m1, double[][] m2) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < m1.length; i++)
			for (int j = 0; j < m1[0].length; j++)
				matrix[i][j] = m1[i][j] - m2[i][j];
		return matrix;
	}

	// Unit Tested
	public static double[][] scalarDivide(double[][] numerator, double denom) {
		double[][] matrix = new double[numerator.length][numerator[0].length];
		for (int i = 0; i < numerator.length; i++)
			for (int j = 0; j < numerator[i].length; j++)
				matrix[i][j] = numerator[i][j] / denom;
		return matrix;
	}

	// Unit Tested
	public static double [] scalarDivide(double numerator, double[] denom) {
		double[] vector = new double[denom.length];
		for (int i = 0; i < denom.length; i++)
			vector[i] = numerator / denom[i] ;
		return vector;
	}

	// Unit Tested
	public static double [] scalarDivide(double[] numerator, double denom) {
		double[] vector = new double[numerator.length];
		for (int i = 0; i < numerator.length; i++)
			vector[i] = numerator[i] / denom;
		return vector;
	}

	// Unit Tested
	public static double [] scalarDivide(double [] numerator, double[] denom) {
		double[] vector = new double[denom.length];
		for (int i = 0; i < denom.length; i++)
			vector[i] = numerator[i] / denom[i] ;
		return vector;
	}

	// Unit Tested
	public static double[][] scalarDivide(double[][] numerator, double[][] denom) {
		double[][] matrix = new double[numerator.length][numerator[0].length];
		for (int i = 0; i < numerator.length; i++)
			for (int j = 0; j < numerator[i].length; j++)
				matrix[i][j] = numerator[i][j] / denom[i][j];
		return matrix;
	}
	
	// Unit Tested
	public static double[][] scalarMult(double[][] m1, double mul) {
		double[][] matrix = new double[m1.length][m1[0].length];
		for (int i = 0; i < m1.length; i++)
			for (int j = 0; j < m1[i].length; j++)
				matrix[i][j] = m1[i][j] * mul;
		return matrix;
	}

	public static double [] scalarMultiply(double[] m1, double mul) {
		double[] matrix = new double[m1.length];
		for (int i = 0; i < m1.length; i++)
			matrix[i] = m1[i] * mul;
		return matrix;
	}

	public static double [] scalarMultiply(double[] m1, double [] m2) {
		double[] matrix = new double[m1.length];
		for (int i = 0; i < m1.length; i++)
			matrix[i] = m1[i] * m2[i];
		return matrix;
	}

	public static double[][] diag(double[][] ds) {
		boolean isLong = ds.length > ds[0].length;
		int dim = Math.max(ds.length,ds[0].length);
		double [][] result = new double [dim][dim];
		for (int i = 0; i < result.length; i++) {
			for (int j = 0; j < result.length; j++) {
				if(i==j) {
					if(isLong)
						result[i][j] = ds[i][0];
					else
						result[i][j] = ds[0][i];
				}
			}
		}

		return result;
	}
	
	public static double[][] diag(double[] ds) {
		int dim = ds.length;
		double [][] result = new double [dim][dim];
		for (int i = 0; i < result.length; i++) {
			result[i][i] = ds[i];
		}

		return result;
	}

	public static double[][] diag(int dim) {
		double [][] result = new double [dim][dim];
		for (int i = 0; i < result.length; i++) {
			result[i][i] = 1.0;
		}

		return result;
	}

	public static double[][] diag(int dim, double diagVal) {
		double [][] result = new double [dim][dim];
		for (int i = 0; i < result.length; i++) {
			result[i][i] = diagVal;
		}

		return result;
	}

	public static double [][] dot(double [][] a, double [][] b) {
		if(a[0].length!=b.length) throw new IllegalArgumentException("Dims does not match: " + a[0].length  + "!=" + b.length);
		double [][] res = new double[a.length][b[0].length];
		for (int row = 0; row < a.length; row++) {
			for (int col = 0; col < b[row].length; col++) {
				for (int i = 0; i < a[0].length; i++) {
					res[row][col] = a[row][i] * b[i][col];
				}
			}
		}
		return res;
	}

	public static double dot(double [] a, double [] b) {
		if(a.length!=b.length) {
			throw new IllegalArgumentException("Vectors are not of equal length");
		}
		double res = 0.0;
		for (int i = 0; i < b.length; i++) {
			res += a[i] * b[i];
		}
		return res;
	}
	
	public static double dot2P1(double [] a1, double [] a2, double [] b) {
		if((a1.length+a2.length)!=b.length) {
			throw new IllegalArgumentException("Vectors are not of equal length: a1 (" 
					+ a1.length + ") + a2 (" + a2.length + ") => (a1+a2)=" + (a1.length + a2.length) + " != (b)=" + b.length);
		}
		double res = 0.0;
		int bidx = 0;
		for (int i = 0; i < a1.length; i++, bidx++) {
			res += a1[i] * b[bidx];
		}
		for (int i = 0; i < a2.length; i++, bidx++) {
			res += a2[i] * b[bidx];
		}
		return res;
	}

	public static int maxIdx(double[] probs) {
		int maxIdx = 0;
		double max = probs[maxIdx];
		for (int i = 0; i < probs.length; i++) {
			if(probs[i]>max) {
				max = probs[i];
				maxIdx = i;
			}
		}
		return maxIdx;
	}

	public static double[][] extractCol(int col, double[][] matrix) {
		double [][] res = new double[matrix.length][1];
		for (int row = 0; row < matrix.length; row++) {
			res[row][0] = matrix[row][col];
		}
		return res;
	}

	public static double[] extractColVector(int col, double[][] matrix) {
		double [] res = new double[matrix.length];
		for (int row = 0; row < matrix.length; row++) {
			res[row] = matrix[row][col];
		}
		return res;
	}

	public static double [][] extractDoubleArray(DenseMatrix64F p) {
		int rows = p.getNumRows();
		int cols = p.getNumCols();
		double [][] result = new double[rows][cols];
		for (int i = 0; i < rows; i++) {
			for (int j = 0; j < cols; j++) {
				result[i][j] = p.get(i, j);
			}
		}
		return result;
	}
	
	public static double [] extractDoubleVector(DenseMatrix64F p) {
		int rows = p.getNumRows();
		int cols = p.getNumCols();
		
		if(rows != 1 && cols != 1) {
			throw new IllegalArgumentException("Cannot convert a " + rows + "x" + cols + " matrix to a vector");
		}
		
		double [] result;
		if(cols == 1) {
			result = new double[rows];
			for (int j = 0; j < rows; j++) {
				result[j] = p.get(j,0);
			}
		} else {
			result = new double[cols];
			for (int j = 0; j < cols; j++) {
				result[j] = p.get(0,j);
			}			
		}
		return result;
	}

	/**
	 * Extracts column <code>col</code> from rows in <code>cJIdxs</code>
	 * @param col
	 * @param matrix
	 * @param rows
	 * @return
	 */
	public static double[][] extractRowCols(int col, double[][] matrix, int[] rows) {
		int [] cols = {col};
		return extractSubMatrix(cols, rows, matrix);
	}
	
	/**
	 * Extracts columns in <code>cols</code> from rows in <code>rows</code>
	 * @param cols
	 * @param rows
	 * @param matrix
	 * @return
	 */
	public static double[][] extractSubMatrix(int [] cols, int [] rows, double[][] zs2) {
		double [][] res = new double[Math.max(1, rows.length)][cols.length];
		for (int row = 0; row < rows.length; row++) {
			for (int col = 0; col < cols.length; col++) {
				res[row][col] = zs2[rows[row]][cols[col]];				
			}
		}
		return res;
	}
	
	/**
	 * Extracts columns <code>cols</code> from<code>matrix</code>
	 * @param cols
	 * @param rows
	 * @param matrix
	 * @return
	 */
	public static double[][] extractCols(int [] cols, double[][] matrix) {
		double [][] res = new double[matrix.length][cols.length];
		for (int row = 0; row < matrix.length; row++) {
			int resCol = 0;
			for (int col = 0; col < cols.length; col++) {
				res[row][resCol++] = matrix[row][cols[col]];				
			}
		}
		return res;
	}
	
	/**
	 * Extracts columns from <code>start</code> (inclusive) to <code>end</code> (exclusive (i.e [start,end) ) ) from<code>matrix</code>
	 * @param cols
	 * @param rows
	 * @param matrix
	 * @return
	 */
	public static double[][] extractCols(int start, int end, double[][] matrix) {
		if(end < 0) throw new IllegalArgumentException("End (" + end + ") has to be positive!");
		if(end < start) throw new IllegalArgumentException("End (" + end + ") has to be >= start (" + start + ")!");
		double [][] res = new double[matrix.length][end-start];
		if(end>start) {
			for (int row = 0; row < matrix.length; row++) {
				int resCol = 0;
				for (int col = start; col < end; col++) {
					res[row][resCol++] = matrix[row][col];				
				}
			}
		}
		return res;
	}

	public static Integer [] indicesOf(int classIdx, int [] ys) {
		List<Integer> indices = new ArrayList<>();
		for (int row = 0; row < ys.length; row++) {
			if(ys[row]==classIdx) {
				indices.add(row);
			}
		}
		return indices.toArray(new Integer[0]);
	}

	public static double[][] makeDesignMatrix(double[][] xstmp) {
		double [][] xs = new double[xstmp.length][xstmp[0].length+1];
		for (int row = 0; row < xs.length; row++) {
			for (int col = 0; col < xs[0].length; col++) {
				if(col==0) {
					xs[row][col] = 1.0;	
				} else {
					xs[row][col] = xstmp[row][col-1];
				}
			}
		}
		return xs;
	}

	public static double [][] pca(double [][] matrix, int dims) {
		PrincipalComponentAnalysis pca = new PrincipalComponentAnalysis();
		return pca.pca(matrix, dims);
	}

	public static double[][] addIntercept(double[][] xs) {
		double [][] result = new double [xs.length][xs[0].length+1];
		for (int i = 0; i < result.length; i++) {
			for (int j = 0; j < result[0].length; j++) {
				if(j==0) {
					result[i][j] = 1.0;
				} else {
					result[i][j] = xs[i][j-1];
				}
			}
		}
		return result;
	}

	public static double[] toPrimitive(Double[] ds) {
		double [] result = new double[ds.length];
		for (int i = 0; i < ds.length; i++) {
			result[i] = ds[i];
		}
		return result;
	}

	public static int[] toPrimitive(Integer[] ds) {
		int [] result = new int[ds.length];
		for (int i = 0; i < ds.length; i++) {
			result[i] = ds[i];
		}
		return result;
	}

	public static Integer[] toObject(int[] ds) {
		Integer [] result = new Integer[ds.length];
		for (int i = 0; i < ds.length; i++) {
			result[i] = ds[i];
		}
		return result;
	}

	public static Double[] toObject(Double[] ds) {
		Double [] result = new Double[ds.length];
		for (int i = 0; i < ds.length; i++) {
			result[i] = ds[i];
		}
		return result;
	}
}
