A java class that creates GaussianMixture model. Created for reference from stackoverflow.

GaussianComponent.java

package com.example.javagmm;

import java.util.Random;


/**
 * <b>Gaussian Component</b>
 *
 * <p>Description: </p>
 * A <code>GaussianMixture</code> consists of several gaussian components, which
 * are modeled by this class. A gaussian component is a n-dimensional gaussian
 * distribution.
 *
 * @see comirva.audio.util.gmm.GaussianMixture
 * @author Klaus Seyerlehner
 * @version 1.1
 */
public final class GaussianComponent implements java.io.Serializable
{
	/**
	 * 
	 */
	private static final long serialVersionUID = 4992376536218581556L;
	
	//fields defining the gaussian componente
	private int dimension;
	private Matrix mean;
	private double componentWeight;
	private Matrix covariances;

	//implementation details (optimizations)
	private double coefficient;
	//private Matrix invCovariances;
	private CholeskyDecomposition chol;
	private Matrix invChol;

	private static Random rnd = new Random();
	/**
	 * don't make the threshold too large, or the correction routine cannot
	 * eliminate the points the corrupt component is focusing on.
	 **/
	private static double SINGULARITY_DETECTION_THRESHOLD = 1.0E-20;


	/**
	 * Creates a gaussian component and checks the component settings for
	 * correctness.
	 *
	 * @param componentWeight double the weight of this component
	 * @param mean Matrix the mean vector of this component
	 * @param covariances Matrix the covariance matrix of this component
	 * @throws IllegalArgumentException thrown if any invaild or incorrect
	 *                                  settings were found
	 */
	public GaussianComponent(double componentWeight, Matrix mean, Matrix covariances) throws IllegalArgumentException
	{
		//check parameters
		if(componentWeight < 0 || componentWeight > 1)
			throw new IllegalArgumentException("the weight of the component must be in the intervall [0,1];");

		if(mean == null || covariances == null)
			throw new IllegalArgumentException("mean and covariances must not be null values;");

		//set weight, mean and covarainces
		this.mean = mean;
		this.covariances = covariances;
		this.componentWeight = componentWeight;

		//set dimension of the component
		dimension = mean.getRowDimension();

		//check if the settings
		try
		{
			covariances.times(mean);

			//actualize optimization fields
			actualizeOptimizationFields();
		}
		catch(Exception e)
		{
			throw new IllegalArgumentException("mean and covariance matrix must have compatible shapes and the covariance matrix must not be singular;");
		}
	}


	/**
	 * Computes the components weight or in other words the prior probability of
	 * being generated by component i.<br>
	 * <br>
	 * [P(C = i)]
	 *
	 * @return double the components weight
	 */
	public double getComponentWeight()
	{
		return componentWeight;
	}


	/**
	 * Returns the probaility of drawing the given sample from this n-dimensional
	 * gaussian distribution.<br>
	 * <br>
	 * [p(x | C = i)]
	 *
	 * @param x Matrix a sample
	 * @return double the probability of this sample with respect to this
	 *                distribution
	 */
	/*public double getSampleProbability(Matrix x)
  {
    Matrix diff, result;
    double value;

    diff = x.minus(mean);
    result = diff.transpose();
    result = result.times(invCovariances);
    result = result.times(diff);
    value = -0.5d * result.get(0,0);
    return coefficient * Math.exp(value);
  }*/


	/**
	 * Returns the probability of drawing the given sample from this n-dimensional
	 * gaussian distribution weighted with the prior probability of this
	 * component.<br>
	 * <br>
	 * [p(x | C = i) * P(C = i)]<br>
	 * <br >
	 * So this is the probability of the sample under this gaussian distribution,
	 * which is only a part of the wohl distribution. This is an optimized
	 * implementation using the cholesky decomposition.
	 *
	 * @param x Matrix a sample
	 * @return double the probability of this sample with respect to this
	 *                distribution weighted with the prior probability of this
	 *                component.
	 */
	public double getWeightedSampleProbability(Matrix x)
	{
		Matrix diff, result;
		double value = 0;;
		double[] row;

		diff = x.minus(mean);
		result = diff.transpose();
		result = result.timesTriangular(invChol);

		row = result.getArray()[0];
		for(int i = 0; i < row.length; i++)
			value += row[i]*row[i];

		value *= -0.5d;
		return coefficient * Math.exp(value) * componentWeight;
	}


	/**
	 * This method performs the maximization step for this component given the
	 * sample points and the esimates <span>p_ij = P(C=i | x_j)</span> for sample
	 * j of being generated by this component under the assumption that this
	 * sample has been drawn from this GMM. So i is fixed to the index of this
	 * component.
	 *
	 * @param samplePoints PointList the sample points
	 * @param p_ij double[] the estimates
	 *
	 * @throws CovarianceSingularityException thrown if this component got
	 *                                        singular druing the maximization step
	 */
	protected void maximise(PointList samplePoints, double[] p_ij) throws CovarianceSingularityException
	{
		//compute new component weight [p_i = SUM over all j: P(C=i | x_j)]
		double p_i = 0;
		for(int j = 0; j < samplePoints.size(); j++)
			p_i += p_ij[j];

		//compute new mean
		Matrix mean = new Matrix(dimension, 1);
		for(int j = 0; j < samplePoints.size(); j++)
		{
			Matrix x = (Matrix) samplePoints.get(j);
			x = x.times(p_ij[j]);
			mean.plusEquals(x);
		}
		mean.timesEquals(1/p_i);

		//compute new covarince matrix
		Matrix covariances = new Matrix(dimension, dimension);
		for(int j = 0; j < samplePoints.size(); j++)
		{
			Matrix x = (Matrix) samplePoints.get(j);
			Matrix diff = x.minus(mean);
			diff = diff.times(diff.transpose());
			//weighted by normalized posteriors
			diff.timesEquals(p_ij[j]);
			covariances.plusEquals(diff);
		}
		covariances.timesEquals(1/p_i);

		//set covarince matrix
		this.covariances = covariances;
		//set new mean
		this.mean = mean;
		//set new component weight
		componentWeight = p_i/samplePoints.size();
		//actualize optimization fields
		actualizeOptimizationFields();
	}


	/**
	 * Returns a vector, whose components are n independent standard normal
	 * variates.
	 *
	 * @param dimension int the number of components (dimensionality) of the
	 *                      vector
	 * @return double[] a vector, whose components are standard normal variates
	 */
	public static double[] getStandardNormalVector(int dimension)
	{
		double[] v = new double[dimension];
		for(int i = 0; i < v.length; i++)
			v[i] = rnd.nextGaussian();
		return v;
	}


	/**
	 * Returns the dimension of this n-dimensional gaussian distribution.
	 *
	 * @return int number of dimensions
	 */
	public int getDimension()
	{
		return this.dimension;
	}


	/**
	 * Prints some information about this componet.
	 * This is for debugging purpose only.
	 */
	public void print()
	{
		System.out.println("mean:");
		mean.print(6, 3);
		System.out.println("covariance matrix:");
		covariances.print(6,3);
		System.out.println("component weight:");
		System.out.println(componentWeight);
		System.out.println("---------------------------------------------");
	}


	/**
	 * For testing purpose only.
	 *
	 * @return Matrix the mean vector
	 */
	public Matrix getMean()
	{
		return mean;
	}


	/**
	 * Recomputes the optimization fields. The opimization fields store values,
	 * that stay the same till the configuration of the component changes.
	 * Consequently they don't have to be computed for each sample, but for each
	 * configuration.
	 *
	 * @throws CovarianceSingularityException thrown if this component's
	 *                                        covariance matrix got singular
	 */
	private void actualizeOptimizationFields() throws CovarianceSingularityException
	{
		double detCovariances = 0.0d;

		//compute cholesky decomposition
		chol = covariances.chol();

		//root of the determinante is the product of all diagonal elements of the cholesky decomposition
		double[][] covars = chol.getL().getArray();
		detCovariances = 1.0d;
		for(int i = 0; i < covars.length; i++)
			detCovariances *= covars[i][i];
		detCovariances *= detCovariances;

		//is covariance matrix is singular?
		if(detCovariances < SINGULARITY_DETECTION_THRESHOLD)
			throw new CovarianceSingularityException(null);

		//get inverse of chol decomposition
		invChol = chol.getL().inverse();

		//invCovariances = covariances.inverse();
		coefficient = 1.0d / (Math.pow((2.0d * Math.PI), dimension/2.0) * Math.pow(detCovariances, 0.5d));
	}

}

GaussianMixture.java

package com.example.javagmm;

import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.lang.ref.SoftReference;
import java.util.Random;



/**
 * <b>Gaussian Mixture Model</b>
 *
 * <p>Description</p>
 * This class implements in combination with <code>GaussianComponent</code> a
 * gaussian mixture model. To create a gaussian mixture model you have to
 * specify means, covariances and the weight for each of the gaussian components.
 * <code>GaussianMixture</code> and GaussianComponent do not only model a GMM,
 * but also support training GMMs using the EM algorithm.<br>
 * <br>
 * One noteable aspect regarding the implementation is the fact that if
 * the covariance matrix of any of the components of this GMM gets singular
 * during the training process a <code>CovarianceSingularityException</code> is
 * thrown. The <code>CovarianceSingularityException</code> contains a reduced
 * <code>PointList</code>. All points belonging to the singular component have
 * been removed. So after the reduction one can try to rerun the training
 * algorithem with the reduced <code>PointList</code>.<br>
 * <br>
 * Another aspect of the design of this class was influenced by the limited
 * memory on real world computers. To improve performace a buffer to store some
 * estimations is used. This buffer is <i>static</i> to reduce garbage
 * collection time and all training processes are synchronized on this buffer.
 * Consequently one can only train one GMM instance at a time.<br>
 * <br>
 * <b>New in version 1.1:</b><br>
 * - The cholesky decomposition is used to speed up computations.
 *
 * @see comirva.audio.util.gmm.GaussianComponent
 * @author Klaus Seyerlehner
 * @version 1.1
 */
public final class GaussianMixture implements java.io.Serializable
{
	/**
	 * 
	 */
	private static final long serialVersionUID = -3450643365494865923L;
	
	private int dimension = 0;                                                                    //dimension of the gmm
	private GaussianComponent[] components = new GaussianComponent[0];                            //gaussian components

	private static double[][] p_ij = new double[1][1];                                           //hard reference to the buffer of current estimate
	private static SoftReference<double[][]> p_ij_SoftRef = new SoftReference<double[][]>(p_ij); //soft reference to the buffer of current estimates                                                //defines the maximum number of training iterations
	private static Random rnd = new Random();                                                     //a random number generator


	/**
	 * This constructor creates a GMM and checks the parameters for plausibility.
	 * The weights, means and covariances of every component are passed as arrays
	 * to the constructor. The i-th component therefore is completely defined by
	 * the i-th entries within these arrays.
	 *
	 * @param componentWeights double[] specifies the components weights
	 * @param means Matrix[] specifies the components mean vectors
	 * @param covariances Matrix[] specifies the components covariance matrices
	 *
	 * @throws IllegalArgumentException if any invalid parameter settings are
	 *                                  detected while checking them
	 */
	public GaussianMixture(double[] componentWeights, Matrix[] means, Matrix[] covariances) throws IllegalArgumentException
	{
		//check if all parameters are valid
		if(componentWeights.length != means.length || means.length != covariances.length || componentWeights.length < 1)
			throw new IllegalArgumentException("all arrays must have the same length with size greater than 0;");

		//create component array
		components = new GaussianComponent[componentWeights.length];

		//check and create the components
		double sum = 0;
		for(int i = 0; i < components.length; i++)
		{
			if(means[i] == null || covariances[i] == null)
				throw new IllegalArgumentException("all mean and covarince matrices must not be null values;");

			sum += componentWeights[i];

			components[i] = new GaussianComponent(componentWeights[i], means[i], covariances[i]);
		}

		//check if the component weights are set correctly
		if( sum < 0.99 || sum > 1.01)
			throw new IllegalArgumentException("the sum over all component weights must be in the interval [0.99, 1.01];");

		//set dimension
		this.dimension = components[0].getDimension();

		//check if all the components have the same dimensions
		for(int i = 0; i < components.length; i++)
			if(components[i].getDimension() != dimension)
				throw new IllegalArgumentException("the dimensions of all components must be the same;");

	}

	/**
	 * This constructor creates a GMM and checks the components for compatibility.
	 * The components themselfs have been checked during their construction.
	 *
	 * @param components GaussianComponent[] an array of gaussian components
	 *
	 * @throws IllegalArgumentException if the passed components are not
	 *                                  compatible
	 */
	public GaussianMixture(GaussianComponent[] components) throws IllegalArgumentException
	{
		if(components == null)
			throw new IllegalArgumentException("the component array must not be null;");

		//check the components
		double sum = 0;
		for(int i = 0; i < components.length; i++)
		{
			if(components[i] == null)
				throw new IllegalArgumentException("all components in the array must not be null;");

			sum += components[i].getComponentWeight();
		}

		//check if the component weights are set correctly
		if( sum < 0.99 || sum > 1.01)
			throw new IllegalArgumentException("the sum over all component weights must be in the interval [0.99, 1.01];");

		this.components = components;
		this.dimension = components[0].getDimension();

		//check if all the components have the same dimensions
		for(int i = 0; i < components.length; i++)
			if(components[i].getDimension() != dimension)
				throw new IllegalArgumentException("the dimensions of all components must be the same;");
	}



	/**
	 * Returns the log likelihood of the points stored in the pointlist under the
	 * assumption the these points where sample from this GMM.<br>
	 * <br>
	 * [SUM over all j: log (SUM over all i:(p(x_j | C = i) * P(C = i)))]
	 *
	 * @param points PointList list of sample points to estimate the log
	 *                              likelihood of
	 * @return double the log likelihood of drawing these samples from this gmm
	 */
	public double getLogLikelihood(PointList points)
	{
		double p = 0;
		for (int j = 0; j < points.size(); j++)
			p += Math.log(getProbability((Matrix) points.get(j)));
		return p;
	}


	/**
	 * Returns the probability of a single sample point under the assumption that
	 * it was draw from the distribution represented by this GMM.<br>
	 * <br>
	 * [SUM over all i:(p(x | C = i) * P(C = i))]
	 *
	 * @param x Matrix a sample point
	 * @return double the probability of the given sample
	 */
	public double getProbability(Matrix x)
	{
		double p = 0;

		for(int i = 0; i < components.length; i++)
			p += components[i].getWeightedSampleProbability(x);

		return p;
	}


	/**
	 * Returns the number of dimensions of the GMM.
	 *
	 * @return int number of dimmensions
	 */
	public int getDimension()
	{
		return dimension;
	}


	/**
	 * Prints some information about this gaussian component.
	 * This is for debugging purpose only.
	 */
	public void print()
	{
		for(int i = 0; i < components.length; i++)
		{
			System.out.println("Component " + i + ":");
			components[i].print();
		}
	}


	/**
	 * For testing purpose only.
	 *
	 * @param numberOfComponent int the number of the component
	 * @return Matrix the mean vector
	 */
	public Matrix getMean(int numberOfComponent)
	{
		return components[numberOfComponent].getMean();
	}


	/**
	 * This method returns a reference to a buffer for storing estimates of the
	 * sample points. The buffer will be reused if possible or reallocated, if
	 * it is too small or if the garbage collector already captured the buffer.
	 *
	 * @param nrComponents int the number of components of the gmm to allocate the
	 *                         buffer for
	 * @param nrSamplePoints int the number of sample points of the gmm to
	 *                           allocate the buffer for
	 */
	protected static void getBuffer(int nrComponents, int nrSamplePoints)
	{
		//get the buffer from the soft ref => now hard ref
		p_ij = p_ij_SoftRef.get();

		if(p_ij == null)
		{
			//reallocate since gc collected the buffer
			p_ij = new double[nrComponents][2*nrSamplePoints];
			p_ij_SoftRef = new SoftReference<double[][]>(p_ij);
		}

		//check if buffer is too small
		if (p_ij[0].length >= nrSamplePoints && p_ij.length >= nrComponents)
			return;

		//to prevent gc runs take double of the current buffer size
		if (p_ij[0].length < nrSamplePoints)
			nrSamplePoints += nrSamplePoints;

		//reallocate since buffer was too small
		p_ij = new double[nrComponents][nrSamplePoints];
		p_ij_SoftRef = new SoftReference<double[][]>(p_ij);

		//run gc to collect old buffer
		System.gc();
	}
	
	
	/**
	 * Reads a GaussianMixture object by deserializing it from disk
	 * @author Florian Schulze
	 * @param path Path of the serialized GMM file
	 * @return The deserialized GMM
	 */
	public static GaussianMixture readGMM(String path) {
		GaussianMixture gmm = null;
		ObjectInputStream ois;
		try {
			ois = new ObjectInputStream(new FileInputStream(path));
			gmm = (GaussianMixture) ois.readObject();
			ois.close();
		} catch (IOException e) {
			e.printStackTrace();
		} catch (ClassNotFoundException e) {
			e.printStackTrace();
		}
		return gmm;
	}
	
	
	/**
	 * Stores a GaussianMixture object on the hard disk by serializing it
	 * @param path Where to store the file
	 */
	public void writeGMM(String path) {
		try {
			ObjectOutputStream oos = new ObjectOutputStream(new FileOutputStream(path));
			oos.writeObject(this);
			oos.flush();
			oos.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}
}

Matrix.java

package com.example.javagmm;

import android.util.Log;

import java.io.PrintWriter;
import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.text.NumberFormat;
import java.util.Locale;
import java.util.Vector;


/**
Jama = Java Matrix class.
<P>
The Java Matrix Class provides the fundamental operations of numerical
linear algebra.  Various constructors create Matrices from two dimensional
arrays of double precision floating point numbers.  Various "gets" and
"sets" provide access to submatrices and matrix elements.  Several methods
implement basic matrix arithmetic, including matrix addition and
multiplication, matrix norms, and element-by-element array operations.
Methods for reading and printing matrices are also included.  All the
operations in this version of the Matrix Class involve real matrices.
Complex matrices may be handled in a future version.
<P>
Five fundamental matrix decompositions, which consist of pairs or triples
of matrices, permutation vectors, and the like, produce results in five
decomposition classes.  These decompositions are accessed by the Matrix
class to compute solutions of simultaneous linear equations, determinants,
inverses and other matrix functions.  The five decompositions are:
<P><UL>
<LI>Cholesky Decomposition of symmetric, positive definite matrices.
<LI>LU Decomposition of rectangular matrices.
<LI>QR Decomposition of rectangular matrices.
<LI>Singular Value Decomposition of rectangular matrices.
<LI>Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices.
</UL>
<DL>
<DT><B>Example of use:</B></DT>
<P>
<DD>Solve a linear system A x = b and compute the residual norm, ||b - A x||.
<P><PRE>
   double[][] vals = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
   Matrix A = new Matrix(vals);
   Matrix b = Matrix.random(3,1);
   Matrix x = A.solve(b);
   Matrix r = A.times(x).minus(b);
   double rnorm = r.normInf();
</PRE></DD>
</DL>

@author The MathWorks, Inc. and the National Institute of Standards and Technology.
@version 5 August 1998
 */

public class Matrix implements Cloneable, java.io.Serializable
{
	/**
	 * 
	 */
	private static final long serialVersionUID = -6800942850845397017L;
	private static final String TAG = "MatrixClass";
	/* ------------------------
	   Class variables
	 * ------------------------ */

	   /** Array for internal storage of elements.
	   @serial internal array storage.
	   */
	   public double[][] A;

	   /** Row and column dimensions.
	   @serial row dimension.
	   @serial column dimension.
	   */
	   private int m, n;


	/* ------------------------
Constructors
	 * ------------------------ */

	/** Construct an m-by-n matrix of zeros.
@param m    Number of rows.
@param n    Number of colums.
	 */

	public Matrix (int m, int n) {
		this.m = m;
		this.n = n;
		A = new double[m][n];
	}

	/** Construct an m-by-n constant matrix.
@param m    Number of rows.
@param n    Number of colums.
@param s    Fill the matrix with this scalar value.
	 */

	public Matrix (int m, int n, double s) {
		this.m = m;
		this.n = n;
		A = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = s;
			}
		}
	}

	/** Construct a matrix from a 2-D array.
@param A    Two-dimensional array of doubles.
@exception  IllegalArgumentException All rows must have the same length
@see        #constructWithCopy
	 */

	public Matrix (double[][] A) {
		m = A.length;
		n = A[0].length;
		for (int i = 0; i < m; i++) {
			if (A[i].length != n) {
				throw new IllegalArgumentException("All rows must have the same length.");
			}
		}
		this.A = A;
	}

	/** Construct a matrix quickly without checking arguments.
@param A    Two-dimensional array of doubles.
@param m    Number of rows.
@param n    Number of colums.
	 */

	public Matrix (double[][] A, int m, int n) {
		this.A = A;
		this.m = m;
		this.n = n;
	}

	/** Construct a matrix from a one-dimensional packed array
@param vals One-dimensional array of doubles, packed by columns (ala Fortran).
@param m    Number of rows.
@exception  IllegalArgumentException Array length must be a multiple of m.
	 */

	public Matrix (double vals[], int m) {
		this.m = m;
		n = (m != 0 ? vals.length/m : 0);
		if (m*n != vals.length) {
			throw new IllegalArgumentException("Array length must be a multiple of m.");
		}
		A = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = vals[i+j*m];
			}
		}
	}

	/**
	 * Construct a matrix from a Vector of double[]
	 * @param rows	A Vector<double[]> containing double[] items. Each double[] item
	 * is a row of the matrix to create. All double[] items must be of the same length.
	 */
	public Matrix (Vector<double[]> rows, boolean clone){
		m = rows.size();
		n = rows.get(0).length;
		A = new double[m][n];
		if(clone){
			for(int i=0; i<m; i++)
				A[i] = rows.get(i).clone();
		}
		else{
			for(int i=0; i<m; i++)
				A[i] = rows.get(i);
		}

		// sanity check:
		for(int i=0; i<m; i++)
			if(A[i].length != n)
				(new IllegalArgumentException("Length of row "+i+" is "+ A[i].length+". Should be "+n)).printStackTrace();
	}

	/* ------------------------
Public Methods
	 * ------------------------ */

	/** Construct a matrix from a copy of a 2-D array.
@param A    Two-dimensional array of doubles.
@exception  IllegalArgumentException All rows must have the same length
	 */

	public static Matrix constructWithCopy(double[][] A) {
		int m = A.length;
		int n = A[0].length;
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			if (A[i].length != n) {
				throw new IllegalArgumentException
				("All rows must have the same length.");
			}
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j];
			}
		}
		return X;
	}

	/** Make a deep copy of a matrix
	 */

	public Matrix copy () {
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j];
			}
		}
		return X;
	}

	/** Clone the Matrix object.
	 */

	public Object clone () {
		return this.copy();
	}

	/** Access the internal two-dimensional array.
@return     Pointer to the two-dimensional array of matrix elements.
	 */

	public double[][] getArray () {
		return A;
	}

	/** Copy the internal two-dimensional array.
@return     Two-dimensional array copy of matrix elements.
	 */

	public double[][] getArrayCopy () {
		double[][] C = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j];
			}
		}
		return C;
	}

	/** Make a one-dimensional column packed copy of the internal array.
@return     Matrix elements packed in a one-dimensional array by columns.
	 */

	public double[] getColumnPackedCopy () {
		double[] vals = new double[m*n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				vals[i+j*m] = A[i][j];
			}
		}
		return vals;
	}

	/** Make a one-dimensional row packed copy of the internal array.
@return     Matrix elements packed in a one-dimensional array by rows.
	 */

	public double[] getRowPackedCopy () {
		double[] vals = new double[m*n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				vals[i*n+j] = A[i][j];
			}
		}
		return vals;
	}

	/** Get row dimension.
@return     m, the number of rows.
	 */

	public int getRowDimension () {
		return m;
	}

	/** Get column dimension.
@return     n, the number of columns.
	 */

	public int getColumnDimension () {
		return n;
	}

	/** Get a single element.
@param i    Row index.
@param j    Column index.
@return     A(i,j)
@exception  ArrayIndexOutOfBoundsException
	 */

	public double get (int i, int j) {
		return A[i][j];
	}

	/** Get a submatrix.
@param i0   Initial row index
@param i1   Final row index
@param j0   Initial column index
@param j1   Final column index
@return     A(i0:i1,j0:j1)
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public Matrix getMatrix (int i0, int i1, int j0, int j1) {
		Matrix X = new Matrix(i1-i0+1,j1-j0+1);
		double[][] B = X.getArray();
		try {
			for (int i = i0; i <= i1; i++) {
				for (int j = j0; j <= j1; j++) {
					B[i-i0][j-j0] = A[i][j];
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
		return X;
	}

	/** Get a submatrix.
@param r    Array of row indices.
@param c    Array of column indices.
@return     A(r(:),c(:))
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public Matrix getMatrix (int[] r, int[] c) {
		Matrix X = new Matrix(r.length,c.length);
		double[][] B = X.getArray();
		try {
			for (int i = 0; i < r.length; i++) {
				for (int j = 0; j < c.length; j++) {
					B[i][j] = A[r[i]][c[j]];
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
		return X;
	}

	/** Get a submatrix.
@param i0   Initial row index
@param i1   Final row index
@param c    Array of column indices.
@return     A(i0:i1,c(:))
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public Matrix getMatrix (int i0, int i1, int[] c) {
		Matrix X = new Matrix(i1-i0+1,c.length);
		double[][] B = X.getArray();
		try {
			for (int i = i0; i <= i1; i++) {
				for (int j = 0; j < c.length; j++) {
					B[i-i0][j] = A[i][c[j]];
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
		return X;
	}

	/** Get a submatrix.
@param r    Array of row indices.
@param j0   Initial column index
@param j1   Final column index
@return     A(r(:),j0:j1)
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public Matrix getMatrix (int[] r, int j0, int j1) {
		Matrix X = new Matrix(r.length,j1-j0+1);
		double[][] B = X.getArray();
		try {
			for (int i = 0; i < r.length; i++) {
				for (int j = j0; j <= j1; j++) {
					B[i][j-j0] = A[r[i]][j];
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
		return X;
	}

	/** Set a single element.
@param i    Row index.
@param j    Column index.
@param s    A(i,j).
@exception  ArrayIndexOutOfBoundsException
	 */

	public void set (int i, int j, double s) {
		A[i][j] = s;
	}

	/** Set a submatrix.
@param i0   Initial row index
@param i1   Final row index
@param j0   Initial column index
@param j1   Final column index
@param X    A(i0:i1,j0:j1)
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public void setMatrix (int i0, int i1, int j0, int j1, Matrix X) {
		try {
			for (int i = i0; i <= i1; i++) {
				for (int j = j0; j <= j1; j++) {
					A[i][j] = X.get(i-i0,j-j0);
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
	}

	/** Set a submatrix.
@param r    Array of row indices.
@param c    Array of column indices.
@param X    A(r(:),c(:))
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public void setMatrix (int[] r, int[] c, Matrix X) {
		try {
			for (int i = 0; i < r.length; i++) {
				for (int j = 0; j < c.length; j++) {
					A[r[i]][c[j]] = X.get(i,j);
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
	}

	/** Set a submatrix.
@param r    Array of row indices.
@param j0   Initial column index
@param j1   Final column index
@param X    A(r(:),j0:j1)
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public void setMatrix (int[] r, int j0, int j1, Matrix X) {
		try {
			for (int i = 0; i < r.length; i++) {
				for (int j = j0; j <= j1; j++) {
					A[r[i]][j] = X.get(i,j-j0);
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
	}

	/** Set a submatrix.
@param i0   Initial row index
@param i1   Final row index
@param c    Array of column indices.
@param X    A(i0:i1,c(:))
@exception  ArrayIndexOutOfBoundsException Submatrix indices
	 */

	public void setMatrix (int i0, int i1, int[] c, Matrix X) {
		try {
			for (int i = i0; i <= i1; i++) {
				for (int j = 0; j < c.length; j++) {
					A[i][c[j]] = X.get(i-i0,j);
				}
			}
		} catch(ArrayIndexOutOfBoundsException e) {
			throw new ArrayIndexOutOfBoundsException("Submatrix indices");
		}
	}

	/** Matrix transpose.
@return    A'
	 */

	public Matrix transpose () {
		Matrix X = new Matrix(n,m);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[j][i] = A[i][j];
			}
		}
		return X;
	}

	/** One norm
@return    maximum column sum.
	 */

	public double norm1 () {
		double f = 0;
		for (int j = 0; j < n; j++) {
			double s = 0;
			for (int i = 0; i < m; i++) {
				s += Math.abs(A[i][j]);
			}
			f = Math.max(f,s);
		}
		return f;
	}

	/** Two norm
@return    maximum singular value.
	 */

	public double norm2 () {
		return (new SingularValueDecomposition(this).norm2());
	}

	/** Infinity norm
@return    maximum row sum.
	 */

	public double normInf () {
		double f = 0;
		for (int i = 0; i < m; i++) {
			double s = 0;
			for (int j = 0; j < n; j++) {
				s += Math.abs(A[i][j]);
			}
			f = Math.max(f,s);
		}
		return f;
	}

	/** Frobenius norm
@return    sqrt of sum of squares of all elements.
	 */

	public double normF () {
		double f = 0;
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				f = MathUtils.hypot(f,A[i][j]);
			}
		}
		return f;
	}

	/**  Unary minus
@return    -A
	 */

	public Matrix uminus () {
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = -A[i][j];
			}
		}
		return X;
	}

	/** C = A + B
@param B    another matrix
@return     A + B
	 */

	public Matrix plus (Matrix B) {
		checkMatrixDimensions(B);
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j] + B.A[i][j];
			}
		}
		return X;
	}

	/** A = A + B
@param B    another matrix
@return     A + B
	 */

	public Matrix plusEquals (Matrix B) {
		checkMatrixDimensions(B);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = A[i][j] + B.A[i][j];
			}
		}
		return this;
	}

	/** C = A - B
@param B    another matrix
@return     A - B
	 */

	public Matrix minus (Matrix B) {
		checkMatrixDimensions(B);
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j] - B.A[i][j];
			}
		}
		return X;
	}

	/** A = A - B
@param B    another matrix
@return     A - B
	 */

	public Matrix minusEquals (Matrix B) {
		checkMatrixDimensions(B);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = A[i][j] - B.A[i][j];
			}
		}
		return this;
	}

	/** Element-by-element multiplication, C = A.*B
@param B    another matrix
@return     A.*B
	 */

	public Matrix arrayTimes (Matrix B) {
		checkMatrixDimensions(B);
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j] * B.A[i][j];
			}
		}
		return X;
	}

	/** Element-by-element multiplication in place, A = A.*B
@param B    another matrix
@return     A.*B
	 */

	public Matrix arrayTimesEquals (Matrix B) {
		checkMatrixDimensions(B);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = A[i][j] * B.A[i][j];
			}
		}
		return this;
	}

	/** Element-by-element right division, C = A./B
@param B    another matrix
@return     A./B
	 */

	public Matrix arrayRightDivide (Matrix B) {
		checkMatrixDimensions(B);
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = A[i][j] / B.A[i][j];
			}
		}
		return X;
	}

	/** Element-by-element right division in place, A = A./B
@param B    another matrix
@return     A./B
	 */

	public Matrix arrayRightDivideEquals (Matrix B) {
		checkMatrixDimensions(B);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = A[i][j] / B.A[i][j];
			}
		}
		return this;
	}

	/** Element-by-element left division, C = A.\B
@param B    another matrix
@return     A.\B
	 */

	public Matrix arrayLeftDivide (Matrix B) {
		checkMatrixDimensions(B);
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = B.A[i][j] / A[i][j];
			}
		}
		return X;
	}

	/** Element-by-element left division in place, A = A.\B
@param B    another matrix
@return     A.\B
	 */

	public Matrix arrayLeftDivideEquals (Matrix B) {
		checkMatrixDimensions(B);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = B.A[i][j] / A[i][j];
			}
		}
		return this;
	}

	/** Multiply a matrix by a scalar, C = s*A
@param s    scalar
@return     s*A
	 */

	public Matrix times (double s) {
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[i][j] = s*A[i][j];
			}
		}
		return X;
	}

	/** Multiply a matrix by a scalar in place, A = s*A
@param s    scalar
@return     replace A by s*A
	 */

	public Matrix timesEquals (double s) {
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = s*A[i][j];
			}
		}
		return this;
	}

	/** Linear algebraic matrix multiplication, A * B
@param B    another matrix
@return     Matrix product, A * B
@exception  IllegalArgumentException Matrix inner dimensions must agree.
	 */

	public Matrix times (Matrix B) {
		Log.i(TAG, "covariance.times(mean)");
		if (B.m != n) {
            Log.i(TAG, "covariance.times(mean) inside\nB.m: " + B.m + " n: " + n + "\n");
			throw new IllegalArgumentException("Matrix inner dimensions must agree.");
		}
        Log.i(TAG, "covariance.times(mean)");
		Matrix X = new Matrix(m,B.n);
        Log.i(TAG, "covariance.times(mean)");
		double[][] C = X.getArray();
        Log.i(TAG, "covariance.times(mean)");
		double[] Bcolj = new double[n];


		for (int j = 0; j < B.n; j++) {
			for (int k = 0; k < n; k++) {
				Bcolj[k] = B.A[k][j];
			}
			for (int i = 0; i < m; i++) {
				double[] Arowi = A[i];
				double s = 0;
				for (int k = 0; k < n; k++) {
					s += Arowi[k]*Bcolj[k];
				}
				C[i][j] = s;
			}
		}
		return X;
	}


	/**
	 * Linear algebraic matrix multiplication, A * B
	 * B being a triangular matrix
	 * <b>Note:</b>
	 * Actually the matrix should be a <b>column orienten, upper triangular
	 * matrix</b> but use the <b>row oriented, lower triangular matrix</b>
	 * instead (transposed), because this is faster due to the easyer array
	 * access.
	 *
	 * @param B    another matrix
	 * @return     Matrix product, A * B
	 *
	 * @exception  IllegalArgumentException Matrix inner dimensions must agree.
	 */
	public Matrix timesTriangular (Matrix B)
	{
		if (B.m != n)
			throw new IllegalArgumentException("Matrix inner dimensions must agree.");

		Matrix X = new Matrix(m,B.n);
		double[][] c = X.getArray();
		double[][] b;
		double s = 0;
		double[] Arowi;
		double[] Browj;

		b = B.getArray();
		//multiply with each row of A
		for (int i = 0; i < m; i++)
		{
			Arowi = A[i];

			//for all columns of B
			for (int j = 0; j < B.n; j++)
			{
				s = 0;
				Browj = b[j];
				//since B being triangular, this loop uses k <= j
				for (int k = 0; k <= j; k++)
				{
					s += Arowi[k] * Browj[k];
				}
				c[i][j] = s;
			}
		}
		return X;
	}


	/**
	 * X.diffEquals() calculates differences between adjacent columns of this
	 * matrix. Consequently the size of the matrix is reduced by one. The result
	 * is stored in this matrix object again.
	 */
	public void diffEquals()
	{
		double[] col = null;
		for(int i = 0; i < A.length; i++)
		{
			col = new double[A[i].length - 1];

			for(int j = 1; j < A[i].length; j++)
				col[j-1] = Math.abs(A[i][j] - A[i][j - 1]);

			A[i] = col;
		}
		n--;
	}

	/**
	 * X.logEquals() calculates the natural logarithem of each element of the
	 * matrix. The result is stored in this matrix object again.
	 */
	public void logEquals()
	{
		for(int i = 0; i < A.length; i++)
			for(int j = 0; j < A[i].length; j++)
				A[i][j] = Math.log(A[i][j]);
	}

	/**
	 * X.powEquals() calculates the power of each element of the matrix. The
	 * result is stored in this matrix object again.
	 */
	public void powEquals(double exp)
	{
		for(int i = 0; i < A.length; i++)
			for(int j = 0; j < A[i].length; j++)
				A[i][j] = Math.pow(A[i][j], exp);
	}

	/**
	 * X.powEquals() calculates the power of each element of the matrix.
	 *
	 * @return Matrix
	 */
	public Matrix pow(double exp)
	{
		Matrix X = new Matrix(m,n);
		double[][] C = X.getArray();

		for (int i = 0; i < m; i++)
			for (int j = 0; j < n; j++)
				C[i][j] = Math.pow(A[i][j], exp);

		return X;
	}





	/**
	 * X.thrunkAtLowerBoundariy(). All values smaller than the given one are set
	 * to this lower boundary.
	 */
	public void thrunkAtLowerBoundary(double value)
	{
		for(int i = 0; i < A.length; i++)
			for(int j = 0; j < A[i].length; j++)
			{
				if(A[i][j] < value)
					A[i][j] = value;
			}
	}

	public CholeskyDecomposition chol () {
		return new CholeskyDecomposition(this);
	}

	/** Singular Value Decomposition
@return     SingularValueDecomposition
@see SingularValueDecomposition
	 */

	public SingularValueDecomposition svd () {
		return new SingularValueDecomposition(this);
	}



	/** Matrix rank
@return     effective numerical rank, obtained from SVD.
	 */

	public int rank () {
		return new SingularValueDecomposition(this).rank();
	}

	/** Matrix condition (2 norm)
@return     ratio of largest to smallest singular value.
	 */

	public double cond () {
		return new SingularValueDecomposition(this).cond();
	}

	/** Matrix trace.
@return     sum of the diagonal elements.
	 */

	public double trace () {
		double t = 0;
		for (int i = 0; i < Math.min(m,n); i++) {
			t += A[i][i];
		}
		return t;
	}

	/** Generate matrix with random elements
@param m    Number of rows.
@param n    Number of colums.
@return     An m-by-n matrix with uniformly distributed random elements.
	 */

	public static Matrix random (int m, int n) {
		Matrix A = new Matrix(m,n);
		double[][] X = A.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				X[i][j] = Math.random();
			}
		}
		return A;
	}

	/** Generate identity matrix
@param m    Number of rows.
@param n    Number of colums.
@return     An m-by-n matrix with ones on the diagonal and zeros elsewhere.
	 */

	public static Matrix identity (int m, int n) {
		Matrix A = new Matrix(m,n);
		double[][] X = A.getArray();
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				X[i][j] = (i == j ? 1.0 : 0.0);
			}
		}
		return A;
	}

	public Matrix solve (Matrix B) {
		return (m == n ? (new LUDecomposition(this)).solve(B) :
			(new QRDecomposition(this)).solve(B));
	}

	/** Solve X*A = B, which is also A'*X' = B'
	   @param B    right hand side
	   @return     solution if A is square, least squares solution otherwise.
	 */

	public Matrix solveTranspose (Matrix B) {
		return transpose().solve(B.transpose());
	}

	/** Matrix inverse or pseudoinverse
	   @return     inverse(A) if A is square, pseudoinverse otherwise.
	 */

	public Matrix inverse () {
		return solve(identity(m,m));
	}


	/** Print the matrix to stdout.   Line the elements up in columns
	 * with a Fortran-like 'Fw.d' style format.
@param w    Column width.
@param d    Number of digits after the decimal.
	 */

	public void print (int w, int d) {
		print(new PrintWriter(System.out,true),w,d); }

	/** Print the matrix to the output stream.   Line the elements up in
	 * columns with a Fortran-like 'Fw.d' style format.
@param output Output stream.
@param w      Column width.
@param d      Number of digits after the decimal.
	 */

	public void print (PrintWriter output, int w, int d) {
		DecimalFormat format = new DecimalFormat();
		format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.US));
		format.setMinimumIntegerDigits(1);
		format.setMaximumFractionDigits(d);
		format.setMinimumFractionDigits(d);
		format.setGroupingUsed(false);
		print(output,format,w+2);
	}

	/** Print the matrix to stdout.  Line the elements up in columns.
	 * Use the format object, and right justify within columns of width
	 * characters.
	 * Note that is the matrix is to be read back in, you probably will want
	 * to use a NumberFormat that is set to US Locale.
@param format A  Formatting object for individual elements.
@param width     Field width for each column.
@see java.text.DecimalFormat#setDecimalFormatSymbols
	 */

	public void print (NumberFormat format, int width) {
		print(new PrintWriter(System.out,true),format,width); }

	// DecimalFormat is a little disappointing coming from Fortran or C's printf.
	// Since it doesn't pad on the left, the elements will come out different
	// widths.  Consequently, we'll pass the desired column width in as an
	// argument and do the extra padding ourselves.

	/** Print the matrix to the output stream.  Line the elements up in columns.
	 * Use the format object, and right justify within columns of width
	 * characters.
	 * Note that is the matrix is to be read back in, you probably will want
	 * to use a NumberFormat that is set to US Locale.
@param output the output stream.
@param format A formatting object to format the matrix elements
@param width  Column width.
@see java.text.DecimalFormat#setDecimalFormatSymbols
	 */

	public void print (PrintWriter output, NumberFormat format, int width) {
		output.println();  // start on new line.
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				String s = format.format(A[i][j]); // format the number
				int padding = Math.max(1,width-s.length()); // At _least_ 1 space
				for (int k = 0; k < padding; k++)
					output.print(' ');
				output.print(s);
			}
			output.println();
		}
		output.println();   // end with blank line.
	}




	/* ------------------------
Private Methods
	 * ------------------------ */

	/** Check if size(A) == size(B) **/

	private void checkMatrixDimensions (Matrix B) {
		if (B.m != m || B.n != n) {
			throw new IllegalArgumentException("Matrix dimensions must agree.");
		}
	}


	/**
	 * Returns the mean values along the specified dimension.
	 *
	 * @param dim
	 *            If 1, then the mean of each column is returned in a row
	 *            vector. If 2, then the mean of each row is returned in a
	 *            column vector.
	 * @return A vector containing the mean values along the specified
	 *         dimension.
	 */
	public Matrix mean(int dim) {
		Matrix result;
		switch (dim) {
		case 1:
			result = new Matrix(1, n);
			for (int currN = 0; currN < n; currN++) {
				for (int currM = 0; currM < m; currM++)
					result.A[0][currN] += A[currM][currN];
				result.A[0][currN] /= m;
			}
			return result;
		case 2:
			result = new Matrix(m, 1);
			for (int currM = 0; currM < m; currM++) {
				for (int currN = 0; currN < n; currN++) {
					result.A[currM][0] += A[currM][currN];
				}
				result.A[currM][0] /= n;
			}
			return result;
		default:
			(new IllegalArgumentException(
					"dim must be either 1 or 2, and not: " + dim))
					.printStackTrace();
			return null;
		}
	}

	/**
	 * Calculate the full covariance matrix.
	 * @return the covariance matrix
	 */
	public Matrix cov() {
		Matrix transe = this.transpose();
		Matrix result = new Matrix(transe.m, transe.m);
		Log.i(TAG, "In Matrix cov()");
		for(int currM = 0; currM < transe.m; currM++){
			for(int currN = currM; currN < transe.m; currN++){
				double covMN = cov(transe.A[currM], transe.A[currN]);
				result.A[currM][currN] = covMN;
				result.A[currN][currM] = covMN;
			}
		}
		return result;
	}

	/**
	 * Calculate the covariance between the two vectors.
	 * @param vec1
	 * @param vec2
	 * @return the covariance between the two vectors.
	 */
	private double cov(double[] vec1, double[] vec2){
		double result = 0;
		int dim = vec1.length;
		if(vec2.length != dim)
			(new IllegalArgumentException("vectors are not of same length")).printStackTrace();
		double meanVec1 = mean(vec1), meanVec2 = mean(vec2);
		for(int i=0; i<dim; i++){
			result += (vec1[i]-meanVec1)*(vec2[i]-meanVec2);
		}
		Log.i(TAG, "In double cov()");
		return result / Math.max(1, dim-1);
		//		int dim = vec1.length;
		//		if(vec2.length != dim)
		//			(new IllegalArgumentException("vectors are not of same length")).printStackTrace();
		//		double[] times = new double[dim];
		//		for(int i=0; i<dim; i++)
		//			times[i] += vec1[i]*vec2[i];
		//		return mean(times) - mean(vec1)*mean(vec2);
	}

	/**
	 * the mean of the values in the double array
	 * @param vec	double values
	 * @return	the mean of the values in vec
	 */
	private double mean(double[] vec){
		double result = 0;
		for(int i=0; i<vec.length; i++)
			result += vec[i];
		return result / vec.length;
	}

	/**
	 * Returns the sum of the component of the matrix.
	 * @return the sum
	 */
	public double sum(){
		double result = 0;
		for(double[] dArr : A)
			for(double d : dArr)
				result += d;
		return result;
	}

	/**
	 * returns a new Matrix object, where each value is set to the absolute value
	 * @return a new Matrix with all values being positive
	 */
	public Matrix abs() {
		Matrix result = new Matrix(m, n); // don't use clone(), as the values are assigned in the loop.
		for(int i=0; i<result.A.length; i++){
			for(int j=0; j<result.A[i].length; j++)
				result.A[i][j] = Math.abs(A[i][j]);
		}
		return result;
	}

}