• Home
• com.actelion.research
• openchemlib
• 2024.9.3
• source code
• MatrixFunctions.java

×

Project download

Many resources are needed to download a project. Please understand that we have to compensate our server costs. Thank you in advance.
Project price only 1 $

You can buy this project and download/modify it how often you want.

com.actelion.research.calc.MatrixFunctions Maven / Gradle / Ivy

/*
 * Copyright (c) 1997 - 2016
 * Actelion Pharmaceuticals Ltd.
 * Gewerbestrasse 16
 * CH-4123 Allschwil, Switzerland
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * 3. Neither the name of the the copyright holder nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

package com.actelion.research.calc;


import com.actelion.research.chem.descriptor.SimilarityCalculatorDoubleArray;
import com.actelion.research.util.DoubleVec;
import com.actelion.research.util.Formatter;
import com.actelion.research.util.datamodel.*;

import java.awt.*;
import java.io.*;
import java.util.*;
import java.util.Base64.Decoder;
import java.util.List;

/**
 * 
Title: MatrixFunctions
 * 
Description: Matrix operations for which a direct access to the matrix
 * member variables is not necessary
 *
 * References
 * http://www.ini.uzh.ch/~fred/java/Matrix.java
 * https://introcs.cs.princeton.edu/java/95linear/Cholesky.java.html
 *
 * @author Modest von Korff
 * @version 1.0
 * 20.11.2003 MvK: Start implementation
 * 10.06.2004 MvK read matrix.
 */

public class MatrixFunctions {

    private static final double TINY = 0.0000001;




    public static Matrix convert2Binary(Matrix A, double thresh) {
        int rows = A.rows();
        int cols = A.cols();
        Matrix B = new Matrix(rows, cols);
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                if(A.get(i,j)>thresh){
                    B.set(i,j,1);
                }
            }
        }
        return B;
    }

    /**
     * Splits the matrix row wise into two matrices.
     * @param A
     * @param row
     * @return
     */
    public static Matrix [] split(Matrix A, int row) {

        int rows = A.rows();
        int cols = A.cols();
        int rowsC = rows-row;

        Matrix B = new Matrix(row, cols);
        Matrix C = new Matrix(rowsC, cols);

        for (int i = 0; i < row; i++) {
            for (int j = 0; j < cols; j++) {
                B.set(i, j, A.get(i,j));
            }
        }

        for (int i = 0; i < rowsC; i++) {
            for (int j = 0; j < cols; j++) {
                C.set(i, j, A.get(i+row,j));
            }
        }

        Matrix [] arr = new Matrix[2];

        arr[0]=B;
        arr[1]=C;

        return arr;
    }

    public static Matrix [] splitCol(Matrix A, int col) {

        int rows = A.rows();
        int cols = A.cols();
        int colsC = cols-col;

        Matrix B = new Matrix(rows, col);
        Matrix C = new Matrix(rows, colsC);

        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < col; j++) {
                B.set(i, j, A.get(i,j));
            }
        }

        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < colsC; j++) {
                C.set(i, j, A.get(i, j+col));
            }
        }

        Matrix [] arr = new Matrix[2];

        arr[0]=B;
        arr[1]=C;

        return arr;
    }

    public static boolean isEqualCol(Matrix A, int colA, Matrix B, int colB, double thresh) {

        boolean equal = true;

        int r = A.rows();

        if(r != B.rows()){
            throw new RuntimeException("Row number differs!");
        }

        for (int i = 0; i < r; i++) {

            double abs = Math.abs(A.get(i,colA)-B.get(i,colB));

            if(abs>thresh){
                equal=false;
                break;
            }
        }

        return equal;
    }


    public static boolean equals(double d1, double d2){
        return (Math.abs(d1-d2)1.0, all others are set to 0.0).
     *
     * @param n the size of the matrix
     * @return the identity matrix
     */
    public static double[][] id(int n){
        double[][] A = new double [n][n];
        for(int i=0; i a){
                    a = Math.abs(R[j][i]);
                    double[] tempLine = R[j]; R[j] = R[i] ; R[i] = tempLine;
                    tempLine = basis[j];
                    basis[j] = basis[i] ;
                    basis[i] = tempLine;
                }
            }
            // Elimination of (i+i)th element of each line >i in R
            for(int j = i+1; j -1; i--) {
                double sum = 0.0;
                for (int j = i + 1; j < n; j++)
                    sum += R[i][j] * invertedA[j][k];
                invertedA[i][k] = (basis[i][k] - sum) / R[i][i];
            }
        }
        return new Matrix(invertedA);
    }


    public static boolean isSymmetric(Matrix A) {
        int rows = A.rows();
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < i; j++) {
                if (A.get(i,j) != A.get(j,i)) return false;
            }
        }
        return true;
    }

    public static boolean isSquare(Matrix A) {
        int rows = A.rows();
        int cols = A.cols();
        return (rows==cols)?true:false;
    }


    /**
     *
     * @param A
     * @return Cholesky factor L of psd matrix A = L L^T
     */
    public static Matrix cholesky(Matrix A) {
        if (!isSquare(A)) {
            throw new RuntimeException("Matrix is not square");
        }
        if (!isSymmetric(A)) {
            throw new RuntimeException("Matrix is not symmetric");
        }

        int rows  = A.rows();
        double[][] L = new double[rows][rows];

        for (int i = 0; i < rows; i++)  {
            for (int j = 0; j <= i; j++) {
                double sum = 0.0;
                for (int k = 0; k < j; k++) {
                    sum += L[i][k] * L[j][k];
                }
                if (i == j)
                    L[i][i] = Math.sqrt(A.get(i,i) - sum);
                else
                    L[i][j] = 1.0 / L[j][j] * (A.get(i,j) - sum);
            }
            if (L[i][i] <= 0) {
                throw new RuntimeException("Matrix not positive definite");
            }
        }

        return new Matrix(L);
    }

    public static Matrix calculateSimilarityMatrixRowWise(Matrix X1, Matrix X2){

        SimilarityCalculatorDoubleArray similarityCalculatorDoubleArray = new SimilarityCalculatorDoubleArray();


        int r1 = X1.rows();
        int r2 = X2.rows();

        Matrix maSim = new Matrix(r1, r2);

        for (int i = 0; i < r1; i++) {

            double [] x1 = X1.getRow(i);

            for (int j = 0; j < r2; j++) {

                double [] x2 = X2.getRow(j);

                double v = similarityCalculatorDoubleArray.getSimilarity(x1, x2);

                maSim.set(i,j,v);
            }
        }

        return maSim;
    }


    public static double [] calculateMaxSimilarity(Matrix XTrain, Matrix Xtest){

        double [] arrSim = new double[Xtest.rows()];

        int r = Xtest.rows();

        for (int i = 0; i < r; i++) {

            double [] arrDescriptor = Xtest.getRow(i);

            IntegerDouble maxSim = MatrixFunctions.calculateMaxSimilarity(XTrain, arrDescriptor);

            arrSim[i] = maxSim.getDouble();

        }

        return arrSim;
    }


    public static IntegerDouble[] calculateMaxSimilarity(Matrix X, double [] arr, int nMostSimilar){

        List li = new ArrayList<>();

        int rows = X.rows();

        for (int i = 0; i < rows; i++) {
            double similarity = DoubleVec.getTanimotoSimilarity(X.getRow(i), arr);
            li.add(new IntegerDouble(i, similarity));
        }

        Collections.sort(li, IntegerDouble.getComparatorDouble());
        Collections.reverse(li);

        int n = Math.min(li.size(), nMostSimilar);

        IntegerDouble [] arrMax = new IntegerDouble[n];

        for (int i = 0; i < n; i++) {

            arrMax[i]=li.get(i);

        }

        return arrMax;
    }

    public static IntegerDouble calculateMaxSimilarity(Matrix X, double [] arr){

        int rows = X.rows();

        IntegerDouble maxSim = new IntegerDouble(-1, 0);
        for (int i = 0; i < rows; i++) {
            double similarity = DoubleVec.getTanimotoSimilarity(X.getRow(i), arr);

            if(similarity>maxSim.getDouble()){
                maxSim.setInteger(i);
                maxSim.setDouble(similarity);
            }
        }

        return maxSim;
    }


    /**
     *
     * @param maQuadratic
     * @return upper triangle from a quadratic matrix as an array.
     */
    public static double [] upperTriangle(Matrix maQuadratic){

        int r = maQuadratic.rows();
        int n = ((r * r)-r) / 2;

        double [] arr= new double[n];

        int cc=0;

        for (int i = 0; i < r; i++) {

            for (int j = i+1; j < r; j++) {
                arr[cc++]=maQuadratic.get(i,j);
            }
        }

        return arr;
    }

    public static Matrix appendRows(List liMatrix) {

        if(liMatrix==null || liMatrix.size()==0){
            return null;
        }

        int cols = liMatrix.get(0).cols();

        int rows = 0;

        for (Matrix m : liMatrix) {
            rows+=m.rows();
        }

        Matrix mAll = new Matrix(rows, cols);

        int offsetRow=0;
        for (Matrix m : liMatrix) {
            mAll.copy(offsetRow, m);
            offsetRow+=m.rows();
        }


        return mAll;
    }


    public static Matrix appendRows(Matrix ma0, Matrix ma1) {

        Matrix ma = new Matrix(ma0.rows() + ma1.rows(), ma0.cols());

        for (int i = 0; i < ma0.rows(); i++) {
            for (int j = 0; j < ma0.cols(); j++) {
                ma.set(i,j, ma0.get(i,j));
            }
        }

        int offsetRows = ma0.rows();

        for (int i = 0; i < ma1.rows(); i++) {
            for (int j = 0; j < ma1.cols(); j++) {
                ma.set(i+offsetRows,j, ma1.get(i,j));
            }
        }

        return ma;
    }

    public static Matrix appendCols(Matrix ma0, Matrix ma1) {

        if(ma0.rows() != ma1.rows()){
            throw new RuntimeException("Number rows differ!");
        }

        Matrix ma = new Matrix(ma0.rows(), ma0.cols()+ma1.cols());

        for (int i = 0; i < ma0.rows(); i++) {
            for (int j = 0; j < ma0.cols(); j++) {
                ma.set(i,j, ma0.get(i,j));
            }
        }

        int offsetCols = ma0.cols();

        for (int i = 0; i < ma1.rows(); i++) {
            for (int j = 0; j < ma1.cols(); j++) {
                ma.set(i,j+offsetCols, ma1.get(i,j));
            }
        }

        return ma;
    }

	/**
	 * List is already used as a constructor for Matrix. So we have to place this method here.
	 * @param li
	 * @return
	 */
	public static Matrix create(List li){

		double [][] a = new double[li.size()][];

		for (int i = 0; i < li.size(); i++) {
			a[i]=li.get(i);
		}

		return new Matrix(a);

	}

    public static int countFieldsBiggerThan(Matrix ma, int row, double thresh) {
    	int cc = 0;
    	for (int i = 0; i < ma.getColDim(); i++) {
			if(ma.get(row,i) > thresh) {
				cc++;
			}
		}
    	return cc;
    }
    
    public static int countFieldsBiggerThanThreshColWise(Matrix ma, int col, double thresh) {
    	int cc = 0;
    	for (int i = 0; i < ma.rows(); i++) {
			if(ma.get(i,col) > thresh) {
				cc++;
			}
		}
    	return cc;
    }
    
    public static Matrix countFieldsBiggerThanThreshColWise(Matrix ma, double thresh) {
    	Matrix maCounts = new Matrix(1, ma.cols());
    	
    	for (int i = 0; i < ma.cols(); i++) {
			int cc = countFieldsBiggerThanThreshColWise(ma, i, thresh);
			
			maCounts.set(0, i, cc);
		}
    	
    	return maCounts;
    }
    
    /**
     * Calculates the inverse Tanimoto coefficient from row wise comparison of the two input matrices.
     * 
     * @param ma1
     * @param ma2
     * @return complete distance matrix calculated between all rows from  the two input matrices.
     */
    public static Matrix getDistanceMatrixTanimotoInv(Matrix ma1, Matrix ma2) {
        Matrix maDist = new Matrix(ma1.getRowDim(), ma2.getRowDim());
        for (int i = 0; i < ma1.getRowDim(); i++) {
            for (int j = 0; j < ma2.getRowDim(); j++) {
                double dist = getDistanceTanimotoInv(ma1, i, ma2, j);
                maDist.set(i,j, dist);
            }
        }
        return maDist;
    }
    
    public static Matrix getDistanceMatrix(List li) {
        Matrix maDist = new Matrix(li.size(), li.size());
        for (int i = 0; i  li1 = new ArrayList();
        List li2 = new ArrayList();
        for (int ii = 0; ii < ma1.getColDim(); ii++) {
            if((ma1.get(row1, ii) != 0) || (ma2.get(row2, ii) != 0)) {
                li1.add(new Double(ma1.get(row1, ii)));
                li2.add(new Double(ma2.get(row2, ii)));
            }
        }

        Matrix maRow1 = new Matrix(true, li1);
        Matrix maRow2 = new Matrix(true, li2);

        double dAtB = maRow1.multiply(0, maRow2, 0);
        double dAtA = maRow1.multiply(0, maRow1, 0);
        double dBtB = maRow2.multiply(0, maRow2, 0);
        dist = dAtB / (dAtA + dBtB - dAtB);
        return 1-dist;
    }


    /**
     * A square-shaped neighborhood that can be used to define a set of cells surrounding a given point.
     * 
     * @param p
     * @param ma
     * @return
     */
    public static List getMooreNeighborhood(Point p, Matrix ma) {
    	List li = new ArrayList();
    	
    	int startX = Math.max(0, p.x - 1);
    	int endX = Math.min(ma.cols(), p.x+2);
    	
    	int startY = Math.max(0, p.y - 1);
    	int endY = Math.min(ma.rows(), p.y+2);
    	
    	for (int i = startY; i < endY; i++) {
			for (int j = startX; j < endX; j++) {
				if(i!=p.y || j!=p.x) {
					if(ma.get(i,j)>0) {
						li.add(new Point(j,i));
					}
				}
			}
		}
    	return li;
    }
    
    public static List getMooreNeighborhood(Point p, int r, Matrix ma) {
    	List li = new ArrayList();
    	
    	int startX = Math.max(0, p.x - r);
    	int endX = Math.min(ma.cols(), p.x+r+1);
    	
    	int startY = Math.max(0, p.y - r);
    	int endY = Math.min(ma.rows(), p.y+r+1);
    	
    	for (int i = startY; i < endY; i++) {
			for (int j = startX; j < endX; j++) {
				if(i!=p.y || j!=p.x) {
					if(ma.get(i,j)>0) {
						li.add(new Point(j,i));
					}
				}
			}
		}
    	return li;
    }


    /**
	 * 
	 * @param ma
	 *            Matrix with ma.rows and 1 col. Containing the min value from
	 *            each row.
	 * @return
	 */
    public static Matrix getRowMinUnique(Matrix ma) {
    	
    	Matrix maTmp = new Matrix(ma);
    	
    	Matrix maMin = new Matrix(maTmp.getRowDim(), 1);
    	
    	for (int i = 0; i < maTmp.getRowDim(); i++) {
        	ScorePoint td = maTmp.getMinPos();
        	maMin.set(td.y,0,td.getScore());
        	maTmp.setRow(td.y, Double.MAX_VALUE);
        	maTmp.setCol(td.x, Double.MAX_VALUE);
		}
    	return maMin;
    }
    
    /**
     * 
     * @param ma
     * @return list with indices Point(row,col) for values>0.
     */
	public static List getPoints(Matrix ma){
		List li = new ArrayList();
		
		for (int i = 0; i < ma.rows(); i++) {
			for (int j = 0; j < ma.cols(); j++) {
				if(ma.get(i, j)>0){
					li.add(new Point(j,i));
				}
			}
		}
		
		return li;
	}
	
	public static List getIndicesUniqueMaxRowWise(Matrix maIn){
		List li = new ArrayList();
		
		Matrix ma = new Matrix(maIn);
		
		for (int i = 0; i < ma.rows(); i++) {
			
			Point p = ma.getMaxIndex();
			
			li.add(p);
			
			int row = p.y;
			
			for (int j = 0; j < ma.cols(); j++) {
				ma.set(row, j, -Double.MAX_VALUE);
			}
		}
		
		return li;
	}
	
	public static List getIndicesUniqueMaxColumnWise(Matrix maIn){
		List li = new ArrayList();
		
		Matrix ma = new Matrix(maIn);
		
		for (int col = 0; col < ma.cols(); col++) {
			
			Point p = ma.getMaxIndex();
			
			li.add(p);
			
			int colMax = p.x;
			
			for (int row = 0; row < ma.rows(); row++) {
				ma.set(row, colMax, -Double.MAX_VALUE);
			}
		}
		
		return li;
	}

    public static Matrix getScaledByFactor(Matrix ma, double fac) {
    	int cols = (int)(ma.cols()*fac+0.5);
    	int rows = (int)(ma.rows()*fac+0.5);
    	
    	Matrix maSc = new Matrix(rows, cols);
    	
    	for (int i = 0; i < rows; i++) {
			for (int j = 0; j < cols; j++) {
				double v = ma.get((int)(i/fac), (int)(j/fac));
				maSc.set(i,j,v);
			}
		}
    	
    	return maSc;
    	
    }

    /**
     * Centers and divides by standard deviation.
     * @param ma
     * @return
     */
    public static Matrix getScaled(Matrix ma) {

	    Matrix maCentered = ma.getCenteredMatrix();

	    Matrix maSD = ma.getStandardDeviationCols();

	    int rows = ma.rows();
	    int cols = ma.cols();

        Matrix maScale = new Matrix(rows, cols);

        for (int i = 0; i < cols; i++) {

            double sdv = maSD.get(0,i);

            for (int j = 0; j < rows; j++) {
                maScale.set(j,i, maCentered.get(j,i) / sdv);
            }
        }

        return maScale;
    }



    /**
     * 03.10.04 MvK
     * @param ma matrix with objects in rows
     * @param k number of desired cluster
     * @return Matrix of cluster centers, k rows and cols equal ma.cols.
     */
    public static Matrix getKMeanClusters(Matrix ma, int k) {

        Matrix maCenters = new Matrix(k,ma.getColDim());

        int iMaxIterations = 100;
        // Array with indices
        // The index specified the corresponding mean cluster
        int[] arrIndex = new int[ma.getRowDim()];

        // Generate the first mean clusters by random
        Random rnd = new Random();
        for (int ii = 0; ii < k; ii++) {
            int rndIndex = rnd.nextInt(ma.getRowDim());
            maCenters.assignRow(ii, ma.getRow(rndIndex));
        }
        // For test
        // maCenters.assignRow(0, ma.getRow(0));
        // maCenters.assignRow(1, ma.getRow(6));
        // maCenters.assignRow(2, ma.getRow(14));

        maCenters = maCenters.getSorted();
        Matrix maCenters2 = new Matrix(maCenters);
        int counter = 0;
        do{
            maCenters = new Matrix(maCenters2);

            for (int ii = 0; ii < arrIndex.length; ii++)
                arrIndex[ii] = -1;

            // Find the next mean cluster for each object in the matrix.
            // Each col represents a mean cluster, each row represents an
            // object.
            Matrix maDist = getDistTanimotoInvReduced(ma, maCenters);
            for (int ii = 0; ii < maDist.getRowDim(); ii++) {
                int index = maDist.getMinRowIndexRND(ii);
                arrIndex[ii] = index;
            }

            // System.out.println("maDist\n" + maDist);

            // Calculate the new mean clusters.
            double[] arrNumObjects = new double[k];
            maCenters2.set(0.0);
            for (int ii = 0; ii < ma.getRowDim(); ii++) {
                int index = arrIndex[ii];
                maCenters2.add2Row(index, ma, ii);
                arrNumObjects[index]++;
            }
            // boolean bEmptyCenter = false;
            for (int ii = 0; ii < maCenters2.getRowDim(); ii++) {
                if(arrNumObjects[ii] > 0)
                    maCenters2.devideRow(ii, arrNumObjects[ii]);
                else {
                    // maCenters2.setRow(ii, -1);
                    maCenters2.assignRow(ii, maCenters.getRow(ii));
                    // bEmptyCenter = true;
                }
            }
/*
            if(bEmptyCenter) {
                String str = "Break because of empty center.\n";
                System.err.print(str);
                break;
            }
*/
            maCenters2 = maCenters2.getSorted();

            // System.out.println("maCenters2\n" + maCenters2);
            counter++;
            if(counter > iMaxIterations) {
                System.err.print("Max num iterations reached.\n");
                // (new Exception()).printStackTrace();
                break;
            }
        } while(!maCenters.equal(maCenters2));

        System.out.println("Number iterations: " + counter);

        return maCenters;
    }
    
    public static final double getCorrPearson(Matrix A, Matrix B) {
        
        final double [] a = A.toArray();
        
        final double [] aCent = ArrayUtilsCalc.getCentered(a);
        
        final double [] aCentNorm = ArrayUtilsCalc.getNormalized(aCent);
        
        final double [] b = B.toArray();
        
        final double [] bCent = ArrayUtilsCalc.getCentered(b);
        
        final double [] bCentNorm = ArrayUtilsCalc.getNormalized(bCent);
        
        final double val = ArrayUtilsCalc.getCorrPearsonStandardized(aCentNorm,bCentNorm);

        return val;
    }

    public static final double getCorrSpearman(Matrix A, Matrix B) {

        final double [] a = A.toArray();
        final double [] b = B.toArray();

        DoubleArray daA = new DoubleArray(a);
        DoubleArray daB = new DoubleArray(b);

        CorrelationCalculator correlationCalculator = new CorrelationCalculator();

        double c = correlationCalculator.calculateCorrelation(daA, daB, CorrelationCalculator.TYPE_SPEARMAN);

        return c;
    }

    public static double getCorrPearson(Matrix A, int col1, int col2) {
        double val = 0;
        Matrix Anorm = A.getCol(col1);
        Anorm = Anorm.getCenteredMatrix();
        Anorm = Anorm.getNormalizedMatrix();

        Matrix Bnorm = A.getCol(col2);
        Bnorm = Bnorm.getCenteredMatrix();
        Bnorm = Bnorm.getNormalizedMatrix();
        val = MatrixFunctions.getCorrPearsonStandardized(Anorm,Bnorm);
        return val;
    }


    private static double getCorrPearsonStandardized(Matrix A, Matrix B) {
        double val = 0;

        double covXY = getCovarianceCentered(A,B);
        double varA = A.getVarianceCentered();
        double varB = B.getVarianceCentered();
        val = covXY / (varA * varB);
        return val;
    }

    public static double getCovariance(Matrix A, Matrix B) {
        double covXY = 0;
        double Amean = A.getMean();
        double Bmean = B.getMean();
        double sum = 0;
        for (int ii = 0; ii < A.getRowDim(); ii++) {
            for (int jj = 0; jj < A.getColDim(); jj++) {
                sum += (A.get(ii,jj) - Amean) * (B.get(ii,jj) - Bmean);
            }
        }
        covXY = sum / (A.getNumElements() - 1);

        return covXY;
    }

    public static double getCovarianceCentered(Matrix A, Matrix B) {
    	
        double covXY = 0;
        
        double sum = 0;
        
    	final int cols = A.cols();
    	
    	final int rows = A.rows();
        
        for (int i = 0; i < rows; i++) {
        	
        	final double [] a = A.getRow(i);
        	
        	final double [] b = B.getRow(i);
        	
            for (int j = 0; j < cols; j++) {
            	
                sum += (a[j] * b[j]);
            }
        }
        
        covXY = sum / (A.getNumElements() - 1);

        return covXY;
    }
    /**
     * generates a matrix with double values between 0 (inclusive) and 1
     * (exclusive).
     * @param rows rows
     * @param cols columns
     * @return matrix
     */
    public static Matrix getRandomMatrix(int rows, int cols) {
        Matrix ma = new Matrix(rows, cols);

        Random rnd = new Random();

        for (int i = 0; i < ma.getRowDim(); i++) {
            for (int j = 0; j < ma.getColDim(); j++) {
                ma.set(i,j, rnd.nextDouble());
            }
        }

        return ma;
    }


    /**
     *
     * @param n
     * @param mean
     * @param sd
     * @return random one column matrix with given mean and standard deviation.
     */
    public static Matrix rnorm(int n, double mean, double sd) {

        Matrix x = new Matrix(n, 1);

        Random random = new Random();

        for (int i = 0; i < n; i++) {
            double v = random.nextGaussian();

            v = v*sd+mean;

            x.set(i,0,v);
        }
        return x;
    }

    /**
     *
     * @param X
     * @return bit wise.
     */
    public static IntVec getNonZeroCols(Matrix X) {

        int s = X.cols() / Integer.SIZE;

        if(X.cols() % Integer.SIZE>0){
            s++;
        }

        IntVec ivMask = new IntVec(s);

        int r = X.rows();
        int c = X.cols();

        for (int i = 0; i < c; i++) {

            boolean zero=true;
            for (int j = 0; j < r; j++) {

                if(Math.abs(X.get(j,i)) > Matrix.TINY16){
                    zero=false;
                    break;
                }
            }
            if(!zero) {
                ivMask.setBit(i);
            }
        }

        return ivMask;
    }

    /**
     * Converts a vector of vectors into doubles, each vector results in a row in
     * the matrix. All vectors have to be of equal length or a runtime  exception
     * is thrown.
     * @param vecvec vector on vectors, has to be converted into doubles
     * @param ma resulting matrix
     */
    public static void vecvec2Matrix(Vector> vecvec, Matrix ma) {

      Iterator> it = vecvec.iterator();
      int iLenCol0 = ((Vector) it.next()).size();

      for( ; it.hasNext(); ) {
        int iLen = ((Vector) it.next()).size();
        if(iLen != iLenCol0) {
          throw new RuntimeException("All vectors must have the same length.");
        }
      }
      int iRows = vecvec.size();

      ma.resize(iRows, iLenCol0);

      it = vecvec.iterator();
      int iRow = 0;
      for( ; it.hasNext(); ) {
        Vector vec = new Vector(((Vector) it.next()));
        for(int ii = 0; ii < vec.size(); ii++) {
          ma.set(iRow, ii, ((Double) vec.get(ii)).doubleValue());
        }
        iRow++;
      }
    }

    public static Matrix readCSV(File fi) throws IOException {

        List li = new ArrayList();

        BufferedReader br = new BufferedReader(new FileReader(fi));

        int cols = -1;

        String l = null;
        while ((l=br.readLine())!=null){

            String [] arr = l.split(",");

            if(cols == -1) {
                cols=arr.length;
            } else {
                if(arr.length!=cols){
                    throw new RuntimeException("Number of columns differ!");
                }
            }

            li.add(arr);
        }

        Matrix ma = new Matrix(li.size(), cols);

        for (int i = 0; i < li.size(); i++) {

            String [] arr = li.get(i);

            for (int j = 0; j < arr.length; j++) {

                double v = Double.parseDouble(arr[j]);

                ma.set(i,j,v);

            }
        }

        return ma;
    }

    public static Matrix read(File fiMatrix) throws IOException {
        return read(fiMatrix, false);
    }

    public static Matrix read(File fiMatrix, boolean skipFirstLine) throws IOException {

        FileInputStream fis = new FileInputStream(fiMatrix);

        Matrix A = read(fis, skipFirstLine);

        fis.close();

        return A;
    }

    public static Matrix read(String s) throws IOException {

        InputStream is = new ByteArrayInputStream(s.getBytes());

        Matrix A = read(is, false);

        is.close();

        return A;
    }

    public static Matrix readAsLineBase64Encoded(String s) throws IOException {

        Decoder decoder = Base64.getDecoder();

        byte [] b = decoder.decode(s);

        InputStream is = new ByteArrayInputStream(b);

        Matrix A = read(is, false);

        is.close();

        return A;
    }


    public static Matrix read(InputStream is) {
        return read(is, false);
    }

    public static Matrix read(InputStream is, boolean skipFirstLine) {

        List li = new ArrayList();

        Scanner scannerLine = new Scanner (is);

        if(skipFirstLine){
            scannerLine.nextLine();
        }

        int rows = 0;
        int cols = 0;
        while(scannerLine.hasNextLine()) {
            ++rows;
            Scanner scannerValue = new Scanner(scannerLine.nextLine());

            DoubleArray arr = null;

            if(cols==0) {
                arr = new DoubleArray();
                while(scannerValue.hasNextDouble()) {
                    double d = scannerValue.nextDouble();
                    arr.add(d);
                    cols++;
                }
            } else {
                arr = new DoubleArray(cols);
                while(scannerValue.hasNextDouble()) {
                    double d = scannerValue.nextDouble();
                    arr.add(d);
                }
            }

            li.add(arr);
        }

        scannerLine.close();

        double[][] a = new double[rows][cols];

        for (int i = 0; i < li.size(); i++) {
            DoubleArray arr = li.get(i);
            for (int j = 0; j < cols; j++) {
                a[i][j] = arr.get(j);
            }
        }

        return new Matrix(a);
    }


    public static void writeQuadraticSymmetricMatrixPairwise(Matrix ma, File fiTxt) throws IOException {

        if(ma.rows() != ma.cols()){
            throw new RuntimeException("Not a quadratic matrix");
        }

        BufferedWriter bw = new BufferedWriter(new FileWriter(fiTxt));

        int n = ma.rows();

        int rows2Write = ((n * n) - n)/2;

        int cc=0;
        for (int i = 0; i < ma.rows(); i++) {

            for (int j = i+1; j < ma.cols(); j++) {

                double v0 = ma.get(i,j);
                double v1 = ma.get(j,i);

                if(v0!=v1){
                    throw new RuntimeException("Matrix not symmetric");
                }

                StringBuilder sb = new StringBuilder();

                sb.append(i);
                sb.append("\t");
                sb.append(j);
                sb.append("\t");
                sb.append(v0);

                if(cc < rows2Write-1){
                    sb.append("\n");
                }
                cc++;

                bw.write(sb.toString());
            }
        }
        bw.close();
    }


    public static void columnIntoDoubleArray(Matrix A, int col, DoubleArray da) {

        da.clear();

        int r = A.rows();

        for (int i = 0; i < r; i++) {
            da.add(A.get(i, col));
        }
    }


    public static List> createIdentifiedObject(Matrix A) {

        int rows = A.rows();

        List> liIdentifiedObject = new ArrayList<>(rows);

        for (int i = 0; i < rows; i++) {

            double [] a = A.getRow(i);

            IdentifiedObject io = new IdentifiedObject<>(a, i);

            liIdentifiedObject.add(io);
        }

        return liIdentifiedObject;
    }
}