package ml.utils;
import static ml.utils.ArrayOperator.allocate1DArray;
import static ml.utils.ArrayOperator.allocate2DArray;
import static ml.utils.ArrayOperator.allocateVector;
import static ml.utils.ArrayOperator.divideAssign;
import static ml.utils.ArrayOperator.quickSort;
import static ml.utils.InPlaceOperator.assign;
import static ml.utils.InPlaceOperator.clear;
import static ml.utils.InPlaceOperator.minusAssign;
import static ml.utils.InPlaceOperator.timesAssign;
import static ml.utils.Printer.disp;
import static ml.utils.Printer.err;
import static ml.utils.Printer.fprintf;
import static ml.utils.Printer.printMatrix;
import static ml.utils.Time.tic;
import static ml.utils.Time.toc;
import static ml.utils.Utility.exit;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.Map.Entry;
import java.util.Random;
import java.util.Set;
import java.util.TreeMap;
import java.util.TreeSet;
import la.decomposition.EigenValueDecomposition;
import la.decomposition.LUDecomposition;
import la.decomposition.QRDecomposition;
import la.decomposition.SingularValueDecomposition;
import la.matrix.DenseMatrix;
import la.matrix.Matrix;
import la.matrix.SparseMatrix;
import la.vector.DenseVector;
import la.vector.SparseVector;
import la.vector.Vector;
import ml.random.MultivariateGaussianDistribution;
public class Matlab {
/**
* @param args
*/
public static void main(String[] args) {
double[] Vec = new double[] {4, 2, 3, 6, 1, 8, 5, 9, 7};
double start = 0;
start = tic();
disp(max(Vec));
System.out.format("Elapsed time: %.9f seconds.%n", toc(start));
start = tic();
double max = Vec[0];
for (int i = 1; i < Vec.length; i++) {
if (max < Vec[i])
max = Vec[i];
}
disp(max);
System.out.format("Elapsed time: %.9f seconds.%n", toc(start));
double[][] data = { {10d,-5d,0d,3d}, {2d,0d,1d,2d}, {1d,6d,0d,5d}};
Matrix A = new DenseMatrix(data);
disp(sigmoid(A));
tic();
/*
* Allocate 1000 x 1000 matrix costs 0.02 seconds!
* In effect, during the iteration of an algorithm,
* memory allocation for temporary matrix variables
* will cost considerable time. OpenCV is right. Never
* allocate memory in the inner part of a function.
*
* Avoid implicitly use the garbage collection from
* within the Java virtual machine.
*/
for (int i = 0; i < 0; i++) {
A = new DenseMatrix(1000, 1000);
}
System.out.format("Elapsed time: %.9f seconds.%n", toc());
int m = 4;
int n = 3;
A = hilb(m, n);
int[] rIndices = new int[] {0, 1, 3, 1, 2, 2, 3, 2, 3};
int[] cIndices = new int[] {0, 0, 0, 1, 1, 2, 2, 3, 3};
double[] values = new double[] {10, 3.2, 3, 9, 7, 8, 8, 7, 7};
int numRows = 4;
int numColumns = 4;
int nzmax = rIndices.length;
A = new SparseMatrix(rIndices, cIndices, values, numRows, numColumns, nzmax);
fprintf("A:%n");
printMatrix(A);
fprintf("sum(A):%n");
disp(sum(A));
fprintf("sum(A, 2):%n");
disp(sum(A, 2));
disp("mean(A):");
disp(mean(A, 1));
disp("std(A):");
disp(std(A, 0, 1));
fprintf("max(A):%n");
disp(max(A)[0]);
fprintf("max(A, 2):%n");
disp(max(A, 2)[0]);
fprintf("min(A):%n");
disp(min(A)[0]);
fprintf("min(A, 2):%n");
disp(min(A, 2)[0]);
A = full(A);
fprintf("A:%n");
disp(A);
fprintf("sum(A):%n");
disp(sum(A));
fprintf("sum(A, 2):%n");
disp(sum(A, 2));
fprintf("max(A):%n");
disp(max(A)[0]);
fprintf("max(A, 2):%n");
disp(max(A, 2)[0]);
fprintf("min(A):%n");
disp(min(A)[0]);
fprintf("min(A, 2):%n");
disp(min(A, 2)[0]);
fprintf("A'A:%n");
disp(A.transpose().mtimes(A));
Matrix[] VD = EigenValueDecomposition.decompose(A.transpose().mtimes(A));
Matrix V = VD[0];
Matrix D = VD[1];
fprintf("V:%n");
printMatrix(V);
fprintf("D:%n");
printMatrix(D);
fprintf("VDV':%n");
disp(V.mtimes(D).mtimes(V.transpose()));
fprintf("A'A:%n");
printMatrix(A.transpose().mtimes(A));
fprintf("V'V:%n");
printMatrix(V.transpose().mtimes((V)));
fprintf("norm(A, 2):%n");
disp(norm(A, 2));
fprintf("rank(A):%n");
disp(rank(A));
Vector V1 = new SparseVector(3);
V1.set(1, 1);
V1.set(2, -1);
disp("V1:");
disp(V1);
Vector V2 = new SparseVector(3);
V2.set(0, -1);
V2.set(2, 1);
disp("V2:");
disp(V2);
fprintf("max(V1, -1)%n");
disp(max(V1, -1));
fprintf("max(V1, 1)%n");
disp(max(V1, 1));
fprintf("max(V1, 0)%n");
disp(max(V1, 0));
fprintf("max(V1, V2)%n");
disp(max(V1, V2));
fprintf("min(V1, -1)%n");
disp(min(V1, -1));
fprintf("min(V1, 1)%n");
disp(min(V1, 1));
fprintf("min(V1, 0)%n");
disp(min(V1, 0));
fprintf("min(V1, V2)%n");
disp(min(V1, V2));
A = new SparseMatrix(rIndices, cIndices, values, numRows, numColumns, nzmax);
disp("A:");
printMatrix(A);
Vector[] Vs = Matlab.sparseMatrix2SparseRowVectors(A);
Matrix S = Matlab.sparseRowVectors2SparseMatrix(Vs);
disp("S:");
printMatrix(S);
Vs = Matlab.sparseMatrix2SparseColumnVectors(S);
S = Matlab.sparseColumnVectors2SparseMatrix(Vs);
disp("S:");
printMatrix(S);
Matrix B = sparse(new DenseMatrix(size(A), 2).times(A.transpose()));
disp("B:");
printMatrix(B);
disp("max(A, 5)");
printMatrix(max(A, 5d));
disp("max(A, -2)");
printMatrix(max(A, -2d));
disp("max(A, B)");
printMatrix(max(A, B));
disp("A:");
printMatrix(A);
disp("B:");
printMatrix(B);
disp("min(A, 5)");
printMatrix(min(A, 5d));
disp("min(A, -2)");
printMatrix(min(A, -2d));
disp("min(A, B)");
printMatrix(min(A, B));
disp("A:");
printMatrix(A);
Matrix[] sortRes = sort((A), 1, "ascend");
disp("Sorted values:");
printMatrix(sortRes[0]);
disp("Sorted indices:");
disp(sortRes[1]);
Vector V3 = new SparseVector(8);
V3.set(1, 6);
V3.set(2, -2);
V3.set(4, 9);
V3.set(6, 8);
disp("V3:");
disp(V3);
double[] IX = sort(V3, "ascend");
disp("Sorted V3:");
disp(V3);
disp("Sorted indices:");
disp(IX);
/*disp("A:");
A.setEntry(2, 3, 0);
A.setEntry(3, 3, 0);
printMatrix(A);
printMatrix(A.getSubMatrix(1, 3, 1, 3));
printMatrix(getRows(A, 1, 3, 2));
printMatrix(getColumns(A, 1, 3));
printMatrix(getColumns(A, new int[] {3, 2}));*/
disp("A:");
printMatrix(A);
disp("repmat(A, 2, 3):");
printMatrix(repmat(A, 2, 3));
disp("vec(A)");
disp(vec(A));
disp("reshape(vec(A), 4, 4)");
printMatrix(reshape(vec(A), 4, 4));
A = full(A);
disp("full(A)");
printMatrix(A);
disp("repmat(A, 2, 3):");
disp(repmat(A, 2, 3));
disp("vec(A)");
disp(vec(A));
disp("reshape(vec(A), 4, 4)");
disp(reshape(vec(A), 4, 4));
B = new DenseMatrix(new double[][] {{3, 2}, {0, 2}});
disp("sparse(A)");
printMatrix(sparse(A));
disp("sparse(B)");
printMatrix(sparse(B));
printMatrix(kron(full(A), full(B)));
printMatrix(kron(full(A), sparse(B)));
printMatrix(kron(sparse(A), full(B)));
printMatrix(kron(sparse(A), sparse(B)));
/*printMatrix(A);
printMatrix(A.getSubMatrix(1, 3, 1, 3));
printMatrix(getRows(A, 1, 3, 2));
printMatrix(getColumns(A, 1, 3));
printMatrix(getColumns(A, new int[] {3, 2}));*/
}
/**
* A constant holding the smallest positive
* nonzero value of type double, 2-1074.
*/
public static double eps = Double.MIN_VALUE;
/**
* A constant holding the positive infinity of type double.
*/
public static double inf = Double.POSITIVE_INFINITY;
/**
* Compute the mean for rows or columns of a real matrix.
*
* @param A a real matrix
*
* @param dim direction, 1 for column-wise, 2 for row-wise
*
* @return mean(A, dim)
*/
/*public static Vector mean(Matrix A, int dim) {
if (dim != 1 && dim != 2) {
err("dim should be 1 or 2.");
exit(1);
}
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (M == 1) { // A is a row matrix
if (dim == 1) {
return A.getRowVector(0);
} else if (dim == 2) {
double mean = sumAll(A) / N;
return new DenseVector(new double[] {mean});
}
} else if (N == 1) { // A is a column matrix
if (dim == 2) {
return A.getColumnVector(0);
} else if (dim == 1) {
double mean = sumAll(A) / M;
return new DenseVector(new double[] {mean});
}
} else { // A is a normal matrix
Vector res = sum(A, dim);
if (dim == 1)
timesAssign(res, 1.0 / M);
else if (dim == 2)
timesAssign(res, 1.0 / N);
}
}*/
/**
* Compute standard deviation.
*
* @param X a real matrix
*
* @param flag 0: n - 1 in the divisor, 1: n in the divisor
*
* @param dim direction, 1 for column-wise, 2 for row-wise
*
* @return a real dense vector
*/
public static DenseVector std(Matrix X, int flag, int dim) {
if (dim != 1 && dim != 2) {
err("dim should be 1 or 2.");
exit(1);
}
int M = X.getRowDimension();
int N = X.getColumnDimension();
DenseVector mean = mean(X, dim);
Matrix meanMat = null;
Matrix temp = null;
if (dim == 1) {
meanMat = rowVector2RowMatrix(mean);
temp = repmat(meanMat, M, 1);
} else {
meanMat = columnVector2ColumnMatrix(mean);
temp = repmat(meanMat, 1, N);
}
minusAssign(temp, X);
timesAssign(temp, temp);
int num = size(X, dim);
double[] res = sum(temp, dim).getPr();
if (flag == 0) { // flag = 0
if (num == 1)
clear(res);
else
divideAssign(res, num - 1);
} else if (flag == 1) { // flag = 1
if (num != 1)
divideAssign(res, num);
}
for (int k = 0; k < res.length; k++)
res[k] = Math.sqrt(res[k]);
return new DenseVector(res);
}
/**
* Compute standard deviation for columns of a real matrix.
*
* @param X a real matrix
*
* @param flag 0: n - 1 in the divisor, 1: n in the divisor
*
* @return a real dense vector
*/
public static DenseVector std(Matrix X, int flag) {
return std(X, flag, 1);
}
/**
* Compute standard deviation for columns of a real matrix with
* n - 1 in the divisor.
*
* @param X a real matrix
*
* @return a real dense vector
*/
public static DenseVector std(Matrix X) {
return std(X, 0);
}
/**
* Set submatrix of A with selected rows and selected columns by elements of B.
* B should have the same shape to the submatrix of A to be set. It is equivalent
* to the syntax A(selectedRows, selectedColumns) = B.
*
* @param A a matrix whose submatrix is to be set
*
* @param selectedRows {@code int[]} holding indices of selected rows
*
* @param selectedColumns {@code int[]} holding indices of selected columns
*
* @param B a matrix to set the submatrix of A
*
*/
public static void setSubMatrix(Matrix A, int[] selectedRows,
int[] selectedColumns, Matrix B) {
int r, c;
for (int i = 0; i < selectedRows.length; i++) {
for (int j = 0; j < selectedColumns.length; j++) {
r = selectedRows[i];
c = selectedColumns[j];
A.setEntry(r, c, B.getEntry(i, j));
}
}
}
/**
* Reshape a matrix to a new shape specified by a two dimensional
* integer array.
*
* @param A a matrix
*
* @param size a two dimensional integer array describing a new shape
*
* @return a new matrix with a shape specified by size
*
*/
public static Matrix reshape(Matrix A, int[] size) {
if (size.length != 2) {
System.err.println("Input vector should have two elements!");
}
int M = size[0];
int N = size[1];
if (M * N != A.getRowDimension() * A.getColumnDimension()) {
System.err.println("Wrong shape!");
exit(1);
}
Matrix res = null;
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
double[][] resData = new double[M][];
double[] resRow = null;
int r, c;
for (int i = 0, shiftI = 0; i < M; i++, shiftI++) {
resData[i] = new double[N];
resRow = resData[i];
for (int j = 0, shiftJ = shiftI; j < N; j++, shiftJ += M) {
r = shiftJ % A.getRowDimension();
c = shiftJ / A.getRowDimension();
resRow[j] = AData[r][c];
}
}
res = new DenseMatrix(resData);
} else if (A instanceof SparseMatrix) {
int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
double[] pr = ((SparseMatrix) A).getPr();
int nnz = ((SparseMatrix) A).getNNZ();
int[] resIr = new int[nnz];
int[] resJc = new int[N + 1];
double[] resPr = new double[nnz];
System.arraycopy(pr, 0, resPr, 0, nnz);
int lastColIdx = -1;
int currentColIdx = 0;
int idx = 0;
for (int j = 0, shiftJ = 0; j < A.getColumnDimension(); j++, shiftJ += A.getRowDimension()) {
for (int k = jc[j]; k < jc[j + 1]; k++) {
idx = ir[k] + shiftJ;
currentColIdx = idx / M;
resIr[k] = idx % M;
while (lastColIdx < currentColIdx) {
resJc[lastColIdx + 1] = k;
lastColIdx++;
}
}
}
resJc[N] = nnz;
res = SparseMatrix.createSparseMatrixByCSCArrays(resIr, resJc, resPr, M, N, nnz);
}
return res;
}
/**
* Reshape a matrix to a new shape specified number of rows and
* columns.
*
* @param A a matrix
*
* @param M number of rows of the new shape
*
* @param N number of columns of the new shape
*
* @return a new M-by-N matrix whose elements are taken columnwise
* from A
*
*/
public static Matrix reshape(Matrix A, int M, int N) {
return reshape(A, new int[]{M, N});
}
/**
* Reshape a vector to a matrix with a shape specified by a two dimensional
* integer array.
*
* @param V a vector
*
* @param size a two dimensional integer array describing a new shape
*
* @return a new matrix with a shape specified by size
*
*/
public static Matrix reshape(Vector V, int[] size) {
if (size.length != 2) {
System.err.println("Input vector should have two elements!");
exit(1);
}
int dim = V.getDim();
if (size[0] * size[1] != dim) {
System.err.println("Wrong shape!");
}
Matrix res = null;
int M = size[0];
int N = size[1];
if (V instanceof DenseVector) {
double[][] resData = new double[M][];
double[] resRow = null;
double[] pr = ((DenseVector) V).getPr();
for (int i = 0, shiftI = 0; i < M; i++, shiftI++) {
resData[i] = new double[N];
resRow = resData[i];
for (int j = 0, shiftJ = shiftI; j < N; j++, shiftJ += M) {
resRow[j] = pr[shiftJ];
}
}
res = new DenseMatrix(resData);
} else if (V instanceof SparseVector) {
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int[] resIr = new int[nnz];
int[] resJc = new int[N + 1];
double[] resPr = new double[nnz];
System.arraycopy(pr, 0, resPr, 0, nnz);
int lastColIdx = -1;
int currentColIdx = 0;
int idx = 0;
for (int k = 0; k < nnz; k++) {
idx = ir[k];
currentColIdx = idx / M;
resIr[k] = idx % M;
while (lastColIdx < currentColIdx) {
resJc[lastColIdx + 1] = k;
lastColIdx++;
}
}
resJc[N] = nnz;
res = SparseMatrix.createSparseMatrixByCSCArrays(resIr, resJc, resPr, M, N, nnz);
}
return res;
}
/**
* Reshape a vector to a new shape specified number of rows and
* columns.
*
* @param V a vector
*
* @param M number of rows of the new shape
*
* @param N number of columns of the new shape
*
* @return a new M-by-N matrix whose elements are taken from V
*
*/
public static Matrix reshape(Vector V, int M, int N) {
int dim = V.getDim();
if (M * N != dim) {
System.err.println("Wrong shape!");
}
Matrix res = null;
if (V instanceof DenseVector) {
double[][] resData = new double[M][];
double[] resRow = null;
double[] pr = ((DenseVector) V).getPr();
for (int i = 0, shiftI = 0; i < M; i++, shiftI++) {
resData[i] = new double[N];
resRow = resData[i];
for (int j = 0, shiftJ = shiftI; j < N; j++, shiftJ += M) {
resRow[j] = pr[shiftJ];
}
}
res = new DenseMatrix(resData);
} else if (V instanceof SparseVector) {
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int[] resIr = new int[nnz];
int[] resJc = new int[N + 1];
double[] resPr = new double[nnz];
System.arraycopy(pr, 0, resPr, 0, nnz);
int lastColIdx = -1;
int currentColIdx = 0;
int idx = 0;
for (int k = 0; k < nnz; k++) {
idx = ir[k];
currentColIdx = idx / M;
resIr[k] = idx % M;
while (lastColIdx < currentColIdx) {
resJc[lastColIdx + 1] = k;
lastColIdx++;
}
}
resJc[N] = nnz;
res = SparseMatrix.createSparseMatrixByCSCArrays(resIr, resJc, resPr, M, N, nnz);
}
return res;
}
/**
* Vectorize a matrix A.
*
* @param A a matrix
*
* @return Vectorization of a matrix A
*/
public static Matrix vec(Matrix A) {
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (N == 1) {
return A;
}
Matrix res = null;
int dim = M * N;
if (A instanceof DenseMatrix) {
double[][] resData = new double[dim][];
double[][] AData = ((DenseMatrix) A).getData();
for (int j = 0, shift = 0; j < N; j++, shift += M) {
for (int i = 0, shiftI = shift; i < M; i++, shiftI++) {
resData[shiftI] = new double[] {AData[i][j]};
}
}
res = new DenseMatrix(resData);
} else if (A instanceof SparseMatrix) {
int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
double[] pr = ((SparseMatrix) A).getPr();
int nnz = ((SparseMatrix) A).getNNZ();
int[] resIr = new int[nnz];
int[] resJc = new int[] {0, nnz};
double[] resPr = new double[nnz];
System.arraycopy(pr, 0, resPr, 0, nnz);
int cnt = 0;
for (int j = 0, shift = 0; j < N; j++, shift += M) {
for (int k = jc[j]; k <jc[j + 1]; k++) {
resIr[cnt++] = ir[k] + shift;
}
}
res = SparseMatrix.createSparseMatrixByCSCArrays(resIr, resJc, resPr, dim, 1, nnz);
}
return res;
}
/**
* Compute the "economy size" matrix singular value decomposition.
*
* @param A a real matrix
*
* @return a matrix array [U, S, V] where U is left orthonormal matrix, S is a
* a diagonal matrix, and V is the right orthonormal matrix such that
* A = U * S * V'
*
*/
public static Matrix[] svd(Matrix A) {
SingularValueDecomposition svdImpl = new SingularValueDecomposition(A);
Matrix U = svdImpl.getU();
Matrix S = svdImpl.getS();
Matrix V = svdImpl.getV();
Matrix[] res = new Matrix[3];
res[0] = U;
res[1] = S;
res[2] = V;
return res;
}
/**
* Generate random samples chosen from the multivariate Gaussian
* distribution with mean MU and covariance SIGMA.
*
* @param MU 1 x d mean vector
*
* @param SIGMA covariance matrix
*
* @param cases number of d dimensional random samples
*
* @return cases-by-d sample matrix subject to the multivariate
* Gaussian distribution N(MU, SIGMA)
*
*/
public static Matrix mvnrnd(Matrix MU, Matrix SIGMA, int cases) {
return MultivariateGaussianDistribution.mvnrnd(MU, SIGMA, cases);
}
/**
* Generate random samples chosen from the multivariate Gaussian
* distribution with mean MU and covariance SIGMA.
*
* @param MU a 1D {@code double} array holding the mean vector
*
* @param SIGMA a 2D {@code double} array holding the covariance matrix
*
* @param cases number of d dimensional random samples
*
* @return cases-by-d sample matrix subject to the multivariate
* Gaussian distribution N(MU, SIGMA)
*
*/
public static Matrix mvnrnd(double[] MU, double[][] SIGMA, int cases) {
return mvnrnd(new DenseMatrix(MU, 2), new DenseMatrix(SIGMA), cases);
}
/**
* Generate random samples chosen from the multivariate Gaussian
* distribution with mean MU and a diagonal covariance SIGMA.
*
* @param MU a 1D {@code double} array holding the mean vector
*
* @param SIGMA a 1D {@code double} array holding the diagonal elements
* of the covariance matrix
*
* @param cases number of d dimensional random samples
*
* @return cases-by-d sample matrix subject to the multivariate
* Gaussian distribution N(MU, SIGMA)
*
*/
public static Matrix mvnrnd(double[] MU, double[] SIGMA, int cases) {
return mvnrnd(new DenseMatrix(MU, 2), diag(SIGMA), cases);
}
/**
* Replicate and tile an array.
*
* @param A a matrix
*
* @param M number of rows to replicate
*
* @param N number of columns to replicate
*
* @return repmat(A, M, N)
*
*/
public static Matrix repmat(Matrix A, int M, int N) {
Matrix res = null;
int nRow = M * A.getRowDimension();
int nCol = N * A.getColumnDimension();
if (A instanceof DenseMatrix) {
double[][] resData = allocate2DArray(nRow, nCol);
double[] resRow = null;
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
int r;
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
r = i % A.getRowDimension();
ARow = AData[r];
for (int k = 0, shift = 0; k < N; k++, shift += A.getColumnDimension()) {
System.arraycopy(ARow, 0, resRow, shift, A.getColumnDimension());
}
}
res = new DenseMatrix(resData);
} else if (A instanceof SparseMatrix) {
int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
double[] pr = ((SparseMatrix) A).getPr();
int nnz = ((SparseMatrix) A).getNNZ();
int resNNZ = nnz * N * M;
int[] resIr = new int[resNNZ];
int[] resJc = new int[nCol + 1];
double[] resPr = new double[resNNZ];
int[] nnzPerColumn = new int[A.getColumnDimension()];
for (int j = 0; j < A.getColumnDimension(); j++) {
nnzPerColumn[j] = M * (jc[j + 1] - jc[j]);
}
resJc[0] = 0;
for (int c = 0; c < nCol; c++) {
int j = c % A.getColumnDimension();
resJc[c + 1] = resJc[c] + nnzPerColumn[j];
}
for (int j = 0, shiftA = 0; j < A.getColumnDimension(); j++) {
int numNNZACol_j = (jc[j + 1] - jc[j]);
int[] irACol_j = new int[numNNZACol_j];
for (int k = 0, shift = shiftA * M; k < N; k++, shift += nnz * M) {
System.arraycopy(ir, shiftA, irACol_j, 0, numNNZACol_j);
for (int i = 0, shift2 = shift; i < M; i++, shift2 += numNNZACol_j) {
System.arraycopy(irACol_j, 0, resIr, shift2, numNNZACol_j);
if (i < M - 1)
for (int t = 0; t < numNNZACol_j; t++)
irACol_j[t] += A.getRowDimension();
System.arraycopy(pr, shiftA, resPr, shift2, numNNZACol_j);
}
}
shiftA += numNNZACol_j;
}
res = SparseMatrix.createSparseMatrixByCSCArrays(resIr, resJc, resPr, nRow, nCol, resNNZ);
}
return res;
}
/**
* Replicate and tile an array.
*
* @param A a matrix
*
* @param size a int[2] vector [M N]
*
* @return repmat(A, size)
*
*/
public static Matrix repmat(Matrix A, int[] size) {
return repmat(A, size[0], size[1]);
}
/**
* M = mean(A, dim) returns the mean values for elements
* along the dimension of A specified by scalar dim.
* For matrices, mean(A, 2) is a column vector containing
* the mean value of each row.
*
* @param X a real matrix
*
* @param dim dimension order
*
* @return mean(A, dim)
*
*/
public static DenseVector mean(Matrix X, int dim) {
int N = size(X, dim);
double[] S = sum(X, dim).getPr();
divideAssign(S, N);
return new DenseVector(S);
}
/**
* Compute eigenvalues and eigenvectors of a symmetric real matrix.
*
* @param A a symmetric real matrix
*
* @param K number of eigenvalues selected
*
* @param sigma either "lm" (largest magnitude) or "sm" (smallest magnitude)
*
* @return a matrix array [V, D], V is the selected K eigenvectors (normalized
* to 1), and D is a diagonal matrix holding selected K eigenvalues.
*
*/
public static Matrix[] eigs(Matrix A, int K, String sigma) {
EigenValueDecomposition eigImpl = new EigenValueDecomposition(A, 1e-6);
Matrix eigV = eigImpl.getV();
Matrix eigD = eigImpl.getD();
/*disp(eigV);
disp(eigD);*/
int N = A.getRowDimension();
Matrix[] res = new Matrix[2];
Vector eigenValueVector = new DenseVector(K);
Matrix eigenVectors = null;
if (sigma.equals("lm")) {
for (int k = 0; k < K; k++)
eigenValueVector.set(k, eigD.getEntry(k, k));
eigenVectors = eigV.getSubMatrix(0, N - 1, 0, K - 1);
} else if (sigma.equals("sm")) {
for (int k = 0; k < K; k++)
eigenValueVector.set(k, eigD.getEntry(N - 1 - k, N - 1 - k));
// eigenVectors = new DenseMatrix(N, K);
double[][] eigenVectorsData = allocate2DArray(N, K);
double[][] eigVData = ((DenseMatrix) eigV).getData();
double[] eigenVectorsRow = null;
double[] eigVRow = null;
// eigenVectors.setColumnVector(k, eigV.getColumnVector(j));
for (int i = 0; i < N; i++) {
eigenVectorsRow = eigenVectorsData[i];
eigVRow = eigVData[i];
for(int j = N - 1, k = 0; k < K ; j--, k++) {
eigenVectorsRow[k] = eigVRow[j];
}
}
eigenVectors = new DenseMatrix(eigenVectorsData);
} else {
System.err.println("sigma should be either \"lm\" or \"sm\"");
System.exit(-1);
}
res[0] = eigenVectors;
res[1] = diag(eigenValueVector);
return res;
}
/**
* Generate a diagonal matrix with its elements of a 1D {@code double}
* array on its main diagonal.
*
* @param V a 1D {@code double} array holding the diagonal elements
*
* @return diag(V)
*
*/
public static Matrix diag(double[] V) {
int d = V.length;
Matrix res = new SparseMatrix(d, d);
for (int i = 0; i < d; i++) {
res.setEntry(i, i, V[i]);
}
return res;
}
/**
* Calculate the element-wise logarithm of a matrix.
*
* @param A a matrix
*
* @return log(A)
*
*/
public static DenseMatrix log(Matrix A) {
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
double[][] resData = allocate2DArray(nRow, nCol);
double[] resRow = null;
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
ARow = AData[i];
for (int j = 0; j < nCol; j++) {
resRow[j] = Math.log(ARow[j]);
}
}
} else if (A instanceof SparseMatrix) {
int[] ic = ((SparseMatrix) A).getIc();
int[] jr = ((SparseMatrix) A).getJr();
int[] valCSRIndices = ((SparseMatrix) A).getValCSRIndices();
double[] pr = ((SparseMatrix) A).getPr();
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
if (jr[i + 1] == jr[i]) {
assign(resRow, Double.NEGATIVE_INFINITY);
continue;
}
int lastIdx = -1;
int currentIdx = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
currentIdx = ic[k];
for (int j = lastIdx + 1; j < currentIdx; j++) {
resRow[j] = Double.NEGATIVE_INFINITY;
}
resRow[currentIdx] = Math.log(pr[valCSRIndices[k]]);
lastIdx = currentIdx;
}
for (int j = lastIdx + 1; j < nCol; j++) {
resRow[j] = Double.NEGATIVE_INFINITY;
}
}
}
/*for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
res.setEntry(i, j, Math.log(A.getEntry(i, j)));
}
}*/
return new DenseMatrix(resData);
}
/**
* Calculate TFIDF of a doc-term-count matrix, each column
* is a data sample.
*
* @param docTermCountMatrix a matrix, each column is a data example
*
* @return TFIDF of docTermCountMatrix
*
*/
public static Matrix getTFIDF(Matrix docTermCountMatrix) {
final int NTerm = docTermCountMatrix.getRowDimension();
final int NDoc = docTermCountMatrix.getColumnDimension();
// Get TF vector
double[] tfVector = new double[NTerm];
for (int i = 0; i < docTermCountMatrix.getRowDimension(); i++) {
tfVector[i] = 0;
for (int j = 0; j < docTermCountMatrix.getColumnDimension(); j++) {
tfVector[i] += docTermCountMatrix.getEntry(i, j) > 0 ? 1 : 0;
}
}
Matrix res = docTermCountMatrix.copy();
for (int i = 0; i < docTermCountMatrix.getRowDimension(); i++) {
for (int j = 0; j < docTermCountMatrix.getColumnDimension(); j++) {
if (res.getEntry(i, j) > 0) {
res.setEntry(i, j, res.getEntry(i, j) * (tfVector[i] > 0 ? Math.log(NDoc / tfVector[i]) : 0));
}
}
}
return res;
}
/**
* Normalize A by columns.
*
* @param A a matrix
*
* @return a column-wise normalized matrix
*/
public static Matrix normalizeByColumns(Matrix A) {
double[] AA = full(sqrt(sum(A.times(A)))).getPr();
Matrix res = A.copy();
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
resRow[j] /= AA[j];
}
}
} else if (res instanceof SparseMatrix) {
// int[] ir = ((SparseMatrix) res).getIr();
int[] jc = ((SparseMatrix) res).getJc();
double[] pr = ((SparseMatrix) res).getPr();
double v = 0;
for (int j = 0; j < N; j++) {
v = AA[j];
for (int k = jc[j]; k < jc[j + 1]; k++) {
pr[k] /= v;
}
}
}
return res;
}
/**
* Random permutation.
* </br>
* randperm(n) returns a random permutation of the integers 1:n.
*
* @param n an integer
*
* @return randperm(n)
*/
public static int[] randperm(int n) {
int[] res = new int[n];
Set<Integer> leftSet = new TreeSet<Integer>();
for (int i = 0; i < n; i++) {
leftSet.add(i);
}
Random generator = new Random();
for (int i = 0; i < n; i++) {
double[] uniformDist = allocateVector(n - i, 1.0 / (n - i));
double rndRealScalor = generator.nextDouble();
double sum = 0;
for (int j = 0, k = 0; j < n; j++) {
if (!leftSet.contains(j))
continue;
sum += uniformDist[k];
if (rndRealScalor <= sum) {
res[i] = j + 1;
leftSet.remove(j);
break;
} else {
k++;
}
}
}
return res;
}
/**
* Find nonzero elements and return their indices.
*
* @param V a real vector
*
* @return an integer array of indices of nonzero elements of V
*
*/
public static int[] find(Vector V) {
int[] indices = null;
if (V instanceof DenseVector) {
ArrayList<Integer> idxList = new ArrayList<Integer>();
double[] pr = ((DenseVector) V).getPr();
double v = 0;
for (int k = 0; k < V.getDim(); k++) {
v = pr[k];
if (v != 0) {
idxList.add(k);
}
}
int nnz = idxList.size();
indices = new int[nnz];
Iterator<Integer> idxIter = idxList.iterator();
int cnt = 0;
while (idxIter.hasNext()) {
indices[cnt++] = idxIter.next();
}
} else if (V instanceof SparseVector) {
((SparseVector) V).clean();
int nnz = ((SparseVector) V).getNNZ();
int[] ir = ((SparseVector) V).getIr();
indices = new int[nnz];
System.arraycopy(ir, 0, indices, 0, nnz);
}
return indices;
}
/**
* Find nonzero elements and return their value, row and column indices.
*
* @param A a matrix
*
* @return a {@code FindResult} data structure which has three instance
* data members:<br/>
* rows: row indices array for non-zero elements of a matrix<br/>
* cols: column indices array for non-zero elements of a matrix<br/>
* vals: values array for non-zero elements of a matrix<br/>
*
*/
public static FindResult find(Matrix A) {
int[] rows = null;
int[] cols = null;
double[] vals = null;
if (A instanceof SparseMatrix) {
((SparseMatrix) A).clean();
int nnz = ((SparseMatrix) A).getNNZ();
rows = new int[nnz];
cols = new int[nnz];
vals = new double[nnz];
int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
double[] pr = ((SparseMatrix) A).getPr();
int cnt = 0;
for (int j = 0; j < A.getColumnDimension(); j++) {
for (int k = jc[j]; k < jc[j + 1]; k++) {
rows[cnt] = ir[k];
cols[cnt] = j;
vals[cnt] = pr[k];
cnt++;
}
}
} else if (A instanceof DenseMatrix) {
int M = A.getRowDimension();
int N = A.getColumnDimension();
ArrayList<Integer> rowIdxList = new ArrayList<Integer>();
ArrayList<Integer> colIdxList = new ArrayList<Integer>();
ArrayList<Double> valList = new ArrayList<Double>();
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
double v = 0;
for (int i = 0; i < M; i++) {
ARow = AData[i];
for (int j = 0; j < N; j++) {
v = ARow[j];
if (v != 0) {
rowIdxList.add(i);
colIdxList.add(j);
valList.add(v);
}
}
}
int nnz = valList.size();
rows = new int[nnz];
cols = new int[nnz];
vals = new double[nnz];
Iterator<Integer> rowIdxIter = rowIdxList.iterator();
Iterator<Integer> colIdxIter = colIdxList.iterator();
Iterator<Double> valIter = valList.iterator();
int cnt = 0;
while (valIter.hasNext()) {
rows[cnt] = rowIdxIter.next();
cols[cnt] = colIdxIter.next();
vals[cnt] = valIter.next();
cnt++;
}
}
return new FindResult(rows, cols, vals);
}
/**
* Compute the element-wise exponential of a matrix
*
* @param A a matrix
*
* @return exp(A)
*
*/
public static Matrix exp(Matrix A) {
int M = A.getRowDimension();
int N = A.getColumnDimension();
Matrix res = new DenseMatrix(M, N, 1);
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
if (A instanceof DenseMatrix) {
double[][] data = A.getData();
double[] row = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
row = data[i];
for (int j = 0; j < N; j++) {
resRow[j] = Math.exp(row[j]);
}
}
} else if (A instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) A).getPr();
int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
for (int j = 0; j < N; j++) {
for (int k = jc[j]; k < jc[j + 1]; k++) {
resData[ir[k]][j] = Math.exp(pr[k]);
}
}
}
return res;
}
/**
* Get the subMatrix containing the elements of the specified rows.
* Rows are indicated counting from 0.
*
* @param A a real matrix
*
* @param startRow initial row index (inclusive)
*
* @param endRow final row index (inclusive)
*
* @return the subMatrix of A containing the data of the specified rows
*
*/
public static Matrix getRows(Matrix A, int startRow, int endRow) {
/*Matrix res = null;
if (A instanceof DenseMatrix) {
double[][] resData = new double[endRow - startRow + 1][];
double[][] AData = ((DenseMatrix) A).getData();
for (int r = startRow, i = 0; r <= endRow; r++, i++) {
resData[i] = AData[r].clone();
}
res = new DenseMatrix(resData);
} else if (A instanceof SparseMatrix) {
Vector[] vectors = sparseMatrix2SparseRowVectors(A);
Vector[] resVectors = new Vector[endRow - startRow + 1];
for (int r = startRow, i = 0; r <= endRow; r++, i++) {
resVectors[i] = vectors[r];
}
res = sparseRowVectors2SparseMatrix(resVectors);
}
return res;*/
return A.getRows(startRow, endRow);
}
/**
* Get the subMatrix containing the elements of the specified rows.
* Rows are indicated counting from 0.
*
* @param A a real matrix
*
* @param selectedRows indices of selected rows
*
* @return the subMatrix of A containing the data of the specified rows
*
*/
public static Matrix getRows(Matrix A, int... selectedRows) {
/*Matrix res = null;
if (A instanceof DenseMatrix) {
double[][] resData = new double[selectedRows.length][];
double[][] AData = ((DenseMatrix) A).getData();
for (int i = 0; i < selectedRows.length; i++) {
resData[i] = AData[selectedRows[i]].clone();
}
res = new DenseMatrix(resData);
} else if (A instanceof SparseMatrix) {
Vector[] vectors = sparseMatrix2SparseRowVectors(A);
Vector[] resVectors = new Vector[selectedRows.length];
for (int i = 0; i < selectedRows.length; i++) {
resVectors[i] = vectors[selectedRows[i]];
}
res = sparseRowVectors2SparseMatrix(resVectors);
}
return res;*/
return A.getRows(selectedRows);
}
/**
* Get the subMatrix containing the elements of the specified columns.
* Columns are indicated counting from 0.
*
* @param A a real matrix
*
* @param startColumn initial column index (inclusive)
*
* @param endColumn final column index (inclusive)
*
* @return the subMatrix of A containing the data of the specified rows
*
*/
public static Matrix getColumns(Matrix A, int startColumn, int endColumn) {
Matrix res = null;
int nRow = A.getRowDimension();
int nCol = endColumn - startColumn + 1;
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
double[][] resData = new double[nRow][nCol];
double[] resRow = null;
for (int r = 0; r < nRow; r++) {
ARow = AData[r];
resRow = new double[nCol];
for (int c = startColumn, j = 0; c <= endColumn; c++, j++) {
resRow[j] = ARow[c];
}
resData[r] = resRow;
}
res = new DenseMatrix(resData);
} else if (A instanceof SparseMatrix) {
Vector[] vectors = sparseMatrix2SparseColumnVectors(A);
Vector[] resVectors = new Vector[nCol];
for (int c = startColumn, j = 0; c <= endColumn; c++, j++) {
resVectors[j] = vectors[c];
}
res = sparseColumnVectors2SparseMatrix(resVectors);
}
return res;
}
/**
* Get the subMatrix containing the elements of the specified columns.
* Columns are indicated counting from 0.
*
* @param A a real matrix
*
* @param selectedColumns indices of selected columns
*
* @return the subMatrix of A containing the data of the specified columns
*
*/
public static Matrix getColumns(Matrix A, int... selectedColumns) {
Matrix res = null;
int nRow = A.getRowDimension();
int nCol = selectedColumns.length;
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
double[][] resData = new double[nRow][nCol];
double[] resRow = null;
for (int r = 0; r < nRow; r++) {
ARow = AData[r];
resRow = new double[nCol];
for (int j = 0; j < nCol; j++) {
resRow[j] = ARow[selectedColumns[j]];
}
resData[r] = resRow;
}
res = new DenseMatrix(resData);
} else if (A instanceof SparseMatrix) {
Vector[] vectors = sparseMatrix2SparseColumnVectors(A);
Vector[] resVectors = new Vector[nCol];
for (int j = 0; j < nCol; j++) {
resVectors[j] = vectors[selectedColumns[j]];
}
res = sparseColumnVectors2SparseMatrix(resVectors);
}
return res;
}
/**
* Concatenate matrices vertically. All matrices in the argument
* list must have the same number of columns.
*
* @param As matrices to be concatenated vertically
*
* @return [A1; A2; ...]
*
*/
public static Matrix vertcat(final Matrix ... As) {
int nM = As.length;
int nRow = 0;
int nCol = 0;
for (int i = 0; i < nM; i++) {
if (As[i] == null)
continue;
nRow += As[i].getRowDimension();
nCol = As[i].getColumnDimension();
}
for (int i = 1; i < nM; i++) {
if (As[i] != null && nCol != As[i].getColumnDimension())
System.err.println("Any matrix in the argument list should either " +
"be empty matrix or have the same number of columns to the others!");
}
if (nRow == 0 || nCol == 0) {
return null;
}
Matrix res = null;
double[][] resData = new double[nRow][];
double[] resRow = null;
int idx = 0;
for (int i = 0; i < nM; i++) {
if (i > 0 && As[i - 1] != null)
idx += As[i - 1].getRowDimension();
if (As[i] == null)
continue;
if (As[i] instanceof DenseMatrix) {
DenseMatrix A = (DenseMatrix) As[i];
double[][] AData = A.getData();
for (int r = 0; r < A.getRowDimension(); r++) {
// res.setRow(idx + r, As[i].getRow(r));
resData[idx + r] = AData[r].clone();
}
} else if (As[i] instanceof SparseMatrix) {
SparseMatrix A = (SparseMatrix) As[i];
double[] pr = A.getPr();
int[] ic = A.getIc();
int[] jr = A.getJr();
int[] valCSRIndices = A.getValCSRIndices();
for (int r = 0; r < A.getRowDimension(); r++) {
resRow = allocate1DArray(nCol, 0);
for (int k = jr[r]; k < jr[r + 1]; k++) {
resRow[ic[k]] = pr[valCSRIndices[k]];
}
resData[idx + r] = resRow;
}
}
}
res = new DenseMatrix(resData);
return res;
}
/**
* Concatenate matrices horizontally. All matrices in the argument
* list must have the same number of rows.
*
* @param As matrices to be concatenated horizontally
*
* @return [A1 A2 ...]
*
*/
public static Matrix horzcat(final Matrix ... As) {
int nM = As.length;
int nCol = 0;
int nRow = 0;
for (int i = 0; i < nM; i++) {
if (As[i] != null) {
nCol += As[i].getColumnDimension();
nRow = As[i].getRowDimension();
}
}
for (int i = 1; i < nM; i++) {
if (As[i] != null && nRow != As[i].getRowDimension())
System.err.println("Any matrix in the argument list should either " +
"be empty matrix or have the same number of rows to the others!");
}
if (nRow == 0 || nCol == 0) {
return null;
}
Matrix res = new DenseMatrix(nRow, nCol, 0);
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
int idx = 0;
for (int r = 0; r < nRow; r++) {
resRow = resData[r];
idx = 0;
for (int i = 0; i < nM; i++) {
if (i > 0 && As[i - 1] != null)
idx += As[i - 1].getColumnDimension();
if (As[i] == null)
continue;
if (As[i] instanceof DenseMatrix) {
DenseMatrix A = (DenseMatrix) As[i];
System.arraycopy(A.getData()[r], 0, resRow, idx, A.getColumnDimension());
} else if (As[i] instanceof SparseMatrix) {
SparseMatrix A = (SparseMatrix) As[i];
double[] pr = A.getPr();
int[] ic = A.getIc();
int[] jr = A.getJr();
int[] valCSRIndices = A.getValCSRIndices();
for (int k = jr[r]; k < jr[r + 1]; k++) {
resRow[idx + ic[k]] = pr[valCSRIndices[k]];
}
}
}
}
return res;
}
/**
* Concatenate matrices along specified dimension.
*
* @param dim specified dimension, can only be either 1 or 2 currently
*
* @param As matrices to be concatenated
*
* @return a concatenation of all the matrices in the argument list
*
*/
public static Matrix cat(int dim, final Matrix... As) {
Matrix res = null;
if (dim == 1)
res = vertcat(As);
else if (dim == 2)
res = horzcat(As);
else
System.err.println("Specified dimension can only be either 1 or 2 currently!");
return res;
}
/**
* Compute the Kronecker tensor product of A and B
*
* @param A a matrix
*
* @param B a matrix
*
* @return Kronecker product of A and B
*
*/
public static Matrix kron(Matrix A, Matrix B) {
Matrix res = null;
int nRowLeft = A.getRowDimension();
int nColLeft = A.getColumnDimension();
int nRowRight = B.getRowDimension();
int nColRight = B.getColumnDimension();
if (A instanceof DenseMatrix && B instanceof DenseMatrix) {
res = new DenseMatrix(nRowLeft * nRowRight, nColLeft * nColRight, 0);
double[][] resData = res.getData();
double[][] BData = B.getData();
for (int i = 0, rShift = 0; i < nRowLeft; i++, rShift += nRowRight) {
for (int j = 0, cShift = 0; j < nColLeft; j++, cShift += nColRight) {
double A_ij = A.getEntry(i, j);
if (A_ij == 0) {
continue;
}
for (int p = 0; p < nRowRight; p++) {
int r = rShift + p;
double[] BRow = BData[p];
double[] resRow = resData[r];
for (int q = 0; q < nColRight; q++) {
if (BRow[q] == 0) {
continue;
}
int c = cShift + q;
resRow[c] = A_ij * BRow[q];
}
}
}
}
} else if (A instanceof DenseMatrix && B instanceof SparseMatrix) {
TreeMap<Pair<Integer, Integer>, Double> map = new TreeMap<Pair<Integer, Integer>, Double>();
int[] ir2 = ((SparseMatrix) B).getIr();
int[] jc2 = ((SparseMatrix) B).getJc();
double[] pr2 = ((SparseMatrix) B).getPr();
for (int i = 0, rShift = 0; i < nRowLeft; i++, rShift += nRowRight) {
for (int j = 0, cShift = 0; j < nColLeft; j++, cShift += nColRight) {
double A_ij = A.getEntry(i, j);
if (A_ij == 0) {
continue;
}
for (int j2 = 0, c = cShift; j2 < nColRight; j2++, c++) {
for (int k2 = jc2[j2]; k2 < jc2[j2 + 1]; k2++) {
if (pr2[k2] == 0) {
continue;
}
int r = rShift + ir2[k2];
map.put(Pair.of(r, c), A_ij * pr2[k2]);
}
}
}
}
res = SparseMatrix.createSparseMatrix(map, nRowLeft * nRowRight, nColLeft * nColRight);
} else if (A instanceof SparseMatrix && B instanceof DenseMatrix) {
TreeMap<Pair<Integer, Integer>, Double> map = new TreeMap<Pair<Integer, Integer>, Double>();
int[] ir1 = ((SparseMatrix) A).getIr();
int[] jc1 = ((SparseMatrix) A).getJc();
double[] pr1 = ((SparseMatrix) A).getPr();
double[][] BData = B.getData();
for (int j1 = 0, cShift = 0; j1 < nColLeft; j1++, cShift += nColRight) {
for (int k1 = jc1[j1]; k1 < jc1[j1 + 1]; k1++) {
if (pr1[k1] == 0) {
continue;
}
int rShift = ir1[k1] * nRowRight;
for (int i = 0; i < nRowRight; i++) {
double[] BRow = BData[i];
int r = rShift + i;
for (int j = 0; j < nColRight; j++) {
if (BRow[j] == 0) {
continue;
}
int c = cShift + j;
map.put(Pair.of(r, c), pr1[k1] * BRow[j]);
}
}
}
}
res = SparseMatrix.createSparseMatrix(map, nRowLeft * nRowRight, nColLeft * nColRight);
} else if (A instanceof SparseMatrix && B instanceof SparseMatrix) {
TreeMap<Pair<Integer, Integer>, Double> map = new TreeMap<Pair<Integer, Integer>, Double>();
int[] ir1 = ((SparseMatrix) A).getIr();
int[] jc1 = ((SparseMatrix) A).getJc();
double[] pr1 = ((SparseMatrix) A).getPr();
int[] ir2 = ((SparseMatrix) B).getIr();
int[] jc2 = ((SparseMatrix) B).getJc();
double[] pr2 = ((SparseMatrix) B).getPr();
for (int j1 = 0, cShift = 0; j1 < nColLeft; j1++, cShift += nColRight) {
for (int k1 = jc1[j1]; k1 < jc1[j1 + 1]; k1++) {
if (pr1[k1] == 0) {
continue;
}
int rShift = ir1[k1] * nRowRight;
for (int j2 = 0, c = cShift; j2 < nColRight; j2++, c++) {
for (int k2 = jc2[j2]; k2 < jc2[j2 + 1]; k2++) {
if (pr2[k2] == 0) {
continue;
}
int r = rShift + ir2[k2];
map.put(Pair.of(r, c), pr1[k1] * pr2[k2]);
}
}
}
}
res = SparseMatrix.createSparseMatrix(map, nRowLeft * nRowRight, nColLeft * nColRight);
}
return res;
}
/**
* Compute the sum of all elements of a matrix.
*
* @param A a matrix
*
* @return sum(sum(A))
*
*/
public static double sumAll(Matrix A) {
return sum(sum(A));
}
/**
* If A is a 1-row or 1-column matrix, then diag(A) is a
* sparse diagonal matrix with elements of A as its main diagonal,
* else diag(A) is a column matrix holding A's diagonal elements.
*
* @param A a matrix
*
* @return diag(A)
*
*/
public static Matrix diag(Matrix A) {
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
Matrix res = null;
if (nRow == 1) {
res = new SparseMatrix(nCol, nCol);
for (int i = 0; i < nCol; i++) {
res.setEntry(i, i, A.getEntry(0, i));
}
} else if (nCol == 1) {
res = new SparseMatrix(nRow, nRow);
for (int i = 0; i < nRow; i++) {
res.setEntry(i, i, A.getEntry(i, 0));
}
} else if (nRow == nCol) {
res = new DenseMatrix(nRow, 1);
for (int i = 0; i < nRow; i++) {
res.setEntry(i, 0, A.getEntry(i, i));
}
}
return res;
}
/**
* Construct a sparse diagonal matrix from a vector.
*
* @param V a real vector
*
* @return diag(V)
*/
public static SparseMatrix diag(Vector V) {
int dim = V.getDim();
SparseMatrix res = new SparseMatrix(dim, dim);
for (int i = 0; i < dim; i++) {
res.setEntry(i, i, V.get(i));
}
return res;
}
/**
* Right array division.
*
* @param A a matrix
*
* @param v a scalar
*
* @return A ./ v
*/
public static Matrix rdivide(Matrix A, double v) {
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
Matrix res = A.copy();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
for (int j = 0; j < nCol; j++) {
resRow[j] /= v;
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
for (int k = 0; k < pr.length; k++) {
pr[k] /= v;
}
}
return res;
}
/**
* Generate an nRow-by-nCol matrix containing pseudo-random values drawn
* from the standard uniform distribution on the open interval (0,1).
*
* @param nRow number of rows
*
* @param nCol number of columns
*
* @return rand(nRow, nCol)
*
*/
public static Matrix rand(int nRow, int nCol) {
Random generator = new Random();
Matrix res = new DenseMatrix(nRow, nCol);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
res.setEntry(i, j, generator.nextDouble());
}
}
return res;
}
/**
* Generate an n-by-n matrix containing pseudo-random values drawn
* from the standard uniform distribution on the open interval (0,1).
*
* @param n number of rows or columns
*
* @return rand(n, n)
*
*/
public static Matrix rand(int n) {
return rand(n, n);
}
/**
* Generate an nRow-by-nCol matrix containing pseudo-random values drawn
* from the standard normal distribution.
*
* @param nRow number of rows
*
* @param nCol number of columns
*
* @return randn(nRow, nCol)
*
*/
public static Matrix randn(int nRow, int nCol) {
Random generator = new Random();
Matrix res = new DenseMatrix(nRow, nCol);
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
res.setEntry(i, j, generator.nextGaussian());
}
}
return res;
}
/**
* Generate an n-by-n matrix containing pseudo-random values drawn
* from the standard normal distribution.
*
* @param n number of rows or columns
*
* @return randn(n, n)
*
*/
public static Matrix randn(int n) {
return randn(n, n);
}
/**
* Signum function.
* <p>
* For each element of X, SIGN(X) returns 1 if the element
* is greater than zero, 0 if it equals zero and -1 if it is
* less than zero.
* </p>
*
* @param A a real matrix
*
* @return sign(A)
*
*/
public static Matrix sign(Matrix A) {
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
Matrix res = A.copy();
double v = 0;
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
for (int j = 0; j < nCol; j++) {
v = resRow[j];
if (v > 0) {
resRow[j] = 1;
} else if (v < 0) {
resRow[j] = -1;
}
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
for (int k = 0; k < pr.length; k++) {
v = pr[k];
if (v > 0) {
pr[k] = 1;
} else if (v < 0) {
pr[k] = -1;
}
}
}
return res;
}
/**
* Compute the squared l2 distance matrix between column vectors in matrix X
* and column vectors in matrix Y.
*
* @param X
* Data matrix with each column being a feature vector.
*
* @param Y
* Data matrix with each column being a feature vector.
*
* @return an n_x X n_y matrix with its (i, j) entry being the squared l2
* distance between i-th feature vector in X and j-th feature
* vector in Y, i.e., || X(:, i) - Y(:, j) ||_2^2
*
*/
@Deprecated
public static Matrix l2DistanceSquare0(Matrix X, Matrix Y) {
int nX = X.getColumnDimension();
int nY = Y.getColumnDimension();
Matrix dist = null;
Matrix part1 = columnVector2ColumnMatrix(sum(times(X, X), 1)).mtimes(ones(1, nY));
Matrix part2 = ones(nX, 1).mtimes(rowVector2RowMatrix(sum(times(Y, Y), 1)));
Matrix part3 = X.transpose().mtimes(Y).times(2);
dist = part1.plus(part2).minus(part3);
Matrix I = lt(dist, 0);
logicalIndexingAssignment(dist, I, 0);
return dist;
}
/**
* Compute the squared l2 distance matrix between row vectors in matrix X
* and row vectors in matrix Y.
*
* @param X
* Data matrix with each row being a feature vector.
*
* @param Y
* Data matrix with each row being a feature vector.
*
* @return an n_x X n_y matrix with its (i, j) entry being the squared l2
* distance between i-th feature vector in X and j-th feature
* vector in Y, i.e., || X(i, :) - Y(j, :) ||_2^2
*
*/
public static Matrix l2DistanceSquare(Matrix X, Matrix Y) {
int nX = X.getRowDimension();
int nY = Y.getRowDimension();
Matrix dist = null;
double[] XX = sum(times(X, X), 2).getPr();
double[] YY = sum(times(Y, Y), 2).getPr();
dist = full(X.mtimes(Y.transpose()).times(-2));
double[][] resData = ((DenseMatrix) dist).getData();
double[] resRow = null;
double s = 0;
double v = 0;
for (int i = 0; i < nX; i++) {
resRow = resData[i];
s = XX[i];
for (int j = 0; j < nY; j++) {
v = resRow[j] + s + YY[j];
// resRow[j] += s + YY[j];
if (v >= 0)
resRow[j] = v;
else
resRow[j] = 0;
}
}
/*Matrix I = lt(dist, 0);
logicalIndexingAssignment(dist, I, 0);*/
return dist;
}
/**
* Compute the squared l2 distance vector between a vector V
* and row vectors in matrix Y.
*
* @param V
* a feature vector
*
* @param Y
* Data matrix with each row being a feature vector
*
* @return an n_y dimensional vector with its i-th entry being the squared l2
* distance between V and the i-th feature vector in Y, i.e., || V - Y(i, :) ||_2^2
*
*/
public static Vector l2DistanceSquare(Vector V, Matrix Y) {
// int nX = 1;
int nY = Y.getRowDimension();
Vector dist = null;
double XX = sum(times(V, V));
double[] YY = sum(times(Y, Y), 2).getPr();
dist = full(Y.operate(V).times(-2));
double[] pr = ((DenseVector) dist).getPr();
double v = 0;
for (int j = 0; j < nY; j++) {
v = pr[j] + XX + YY[j];
if (v >= 0)
pr[j] = v;
else
pr[j] = 0;
// pr[j] += XX + YY[j];
}
/*Matrix I = lt(dist, 0);
logicalIndexingAssignment(dist, I, 0);*/
return dist;
}
/**
* Compute the squared l2 distance matrix between column vectors in matrix X
* and column vectors in matrix Y.
*
* @param X
* Data matrix with each column being a feature vector.
*
* @param Y
* Data matrix with each column being a feature vector.
*
* @return an n_x X n_y matrix with its (i, j) entry being the squared l2
* distance between i-th feature vector in X and j-th feature
* vector in Y, i.e., || X(:, i) - Y(:, j) ||_2^2
*
*/
public static Matrix l2DistanceSquareByColumns(Matrix X, Matrix Y) {
int nX = X.getColumnDimension();
int nY = Y.getColumnDimension();
Matrix dist = null;
double[] XX = sum(times(X, X)).getPr();
double[] YY = sum(times(Y, Y)).getPr();
dist = full(X.transpose().mtimes(Y).times(-2));
double[][] resData = ((DenseMatrix) dist).getData();
double[] resRow = null;
double s = 0;
double v = 0;
for (int i = 0; i < nX; i++) {
resRow = resData[i];
s = XX[i];
for (int j = 0; j < nY; j++) {
v = resRow[j] + s + YY[j];
// resRow[j] += XX[i] + YY[j];
if (v >= 0)
resRow[j] = v;
else
resRow[j] = 0;
}
}
/*Matrix I = lt(dist, 0);
logicalIndexingAssignment(dist, I, 0);*/
return dist;
}
/**
* Compute the squared l2 distance matrix between two sets of vectors X and Y.
*
* @param X
* Data vectors
*
* @param Y
* Data vectors
*
* @return an n_x X n_y matrix with its (i, j) entry being the squared l2
* distance between i-th feature vector in X and j-th feature
* vector in Y, i.e., || X[i] - Y[j] ||_2^2
*
*/
public static Matrix l2DistanceSquare(Vector[] X, Vector[] Y) {
int nX = X.length;
int nY = Y.length;
Matrix dist = null;
double[] XX = new double[nX];
Vector V = null;
for (int i = 0; i < nX; i++) {
V = X[i];
XX[i] = sum(V.times(V));
}
double[] YY = new double[nY];
for (int i = 0; i < nY; i++) {
V = Y[i];
YY[i] = sum(V.times(V));
}
double[][] resData = allocate2DArray(nX, nY, 0);
double[] resRow = null;
double s = 0;
double v = 0;
for (int i = 0; i < nX; i++) {
resRow = resData[i];
V = X[i];
s = XX[i];
for (int j = 0; j < nY; j++) {
v = s + YY[j] - 2 * innerProduct(V, Y[j]);;
// resRow[j] = s + YY[j] - 2 * innerProduct(V, Y[j]);
if (v >= 0)
resRow[j] = v;
else
resRow[j] = 0;
}
}
/*Matrix I = lt(dist, 0);
logicalIndexingAssignment(dist, I, 0);*/
dist = new DenseMatrix(resData);
return dist;
}
/**
* Compute the l2 distance matrix between column vectors in matrix X
* and column vectors in matrix Y.
*
* @param X
* Data matrix with each column being a feature vector.
*
* @param Y
* Data matrix with each column being a feature vector.
*
* @return an n_x X n_y matrix with its (i, j)th entry being the l2
* distance between i-th feature vector in X and j-th feature
* vector in Y, i.e., || X(:, i) - Y(:, j) ||_2
*
*/
public static Matrix l2DistanceByColumns(Matrix X, Matrix Y) {
return sqrt(l2DistanceSquareByColumns(X, Y));
}
/**
* Compute the l2 distance matrix between row vectors in matrix X
* and row vectors in matrix Y.
*
* @param X
* Data matrix with each row being a feature vector.
*
* @param Y
* Data matrix with each row being a feature vector.
*
* @return an n_x X n_y matrix with its (i, j)th entry being the l2
* distance between i-th feature vector in X and j-th feature
* vector in Y, i.e., || X(i, :) - Y(j, :) ||_2
*
*/
public static Matrix l2Distance(Matrix X, Matrix Y) {
return sqrt(l2DistanceSquare(X, Y));
}
/**
* Compute the l2 distance vector between a vector V
* and row vectors in matrix Y.
*
* @param V
* a feature vector
*
* @param Y
* Data matrix with each row being a feature vector
*
* @return an n_y dimensional vector with its i-th entry being the l2
* distance between V and the i-th feature vector in Y, i.e., || V - Y(i, :) ||_2
*
*/
public static Vector l2Distance(Vector V, Matrix Y) {
return sqrt(l2DistanceSquare(V, Y));
}
/**
* Compute the l2 distance matrix between two sets of vectors X and Y.
*
* @param X
* Data vectors
*
* @param Y
* Data vectors
*
* @return an n_x X n_y matrix with its (i, j)th entry being the l2
* distance between i-th feature vector in X and j-th feature
* vector in Y, i.e., || X[i] - Y[j] ||_2
*
*/
public static Matrix l2Distance(Vector[] X, Vector[] Y) {
return sqrt(l2DistanceSquare(X, Y));
}
/**
* Calculate element by element division between a scalar and a vector.
*
* @param v a real scalar
*
* @param V a real vector
*
* @return v ./ V
*/
public static Vector dotDivide(double v, Vector V) {
int dim = V.getDim();
DenseVector res = (DenseVector) full(V).copy();
double pr[] = ((DenseVector) res).getPr();
for (int k = 0; k < dim; k++) {
pr[k] = v / pr[k];
}
return res;
}
/**
* Calculate square root for all elements of a vector V.
*
* @param V a real vector
*
* @return sqrt(V)
*/
public static Vector sqrt(Vector V) {
int dim = V.getDim();
Vector res = V.copy();
if (res instanceof DenseVector) {
double pr[] = ((DenseVector) res).getPr();
for (int k = 0; k < dim; k++) {
pr[k] = Math.sqrt(pr[k]);
}
} else if (res instanceof SparseVector) {
double pr[] = ((SparseVector) res).getPr();
int nnz = ((SparseVector) res).getNNZ();
for (int k = 0; k < nnz; k++) {
pr[k] = Math.sqrt(pr[k]);
}
}
return res;
}
/**
* Calculate square root for all elements of a matrix A.
*
* @param A a matrix
*
* @return sqrt(A)
*
*/
public static Matrix sqrt(Matrix A) {
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
Matrix res = A.copy();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
for (int j = 0; j < nCol; j++) {
resRow[j] = Math.sqrt(resRow[j]);
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
for (int k = 0; k < pr.length; k++) {
pr[k] = Math.sqrt(pr[k]);
}
}
return res;
}
// public static Matrix repmat(Vector)
public static Matrix rowVector2RowMatrix(Vector V) {
Matrix res = null;
if (V instanceof DenseVector) {
res = denseRowVectors2DenseMatrix(new Vector[] {V});
} else if (V instanceof SparseVector) {
res = sparseRowVectors2SparseMatrix(new Vector[] {V});
}
return res;
}
public static Matrix columnVector2ColumnMatrix(Vector V) {
Matrix res = null;
if (V instanceof DenseVector) {
res = denseColumnVectors2DenseMatrix(new Vector[] {V});
} else if (V instanceof SparseVector) {
res = sparseColumnVectors2SparseMatrix(new Vector[] {V});
}
return res;
}
public static Matrix denseRowVectors2DenseMatrix(Vector[] Vs) {
int M = Vs.length;
// int N = Vs[0].getDim();
double[][] resData = new double[M][];
for (int i = 0; i < M; i++) {
resData[i] = ((DenseVector) Vs[i]).getPr().clone();
}
return new DenseMatrix(resData);
}
public static Matrix denseColumnVectors2DenseMatrix(Vector[] Vs) {
int N = Vs.length;
int M = Vs[0].getDim();
double[][] resData = new double[M][];
for (int i = 0; i < M; i++) {
resData[i] = new double[N];
}
for (int j = 0; j < N; j++) {
double[] column = ((DenseVector) Vs[j]).getPr().clone();
for (int i = 0; i < M; i++) {
resData[i][j] = column[i];
}
}
return new DenseMatrix(resData);
}
public static Vector[] denseMatrix2DenseRowVectors(Matrix A) {
int M = A.getRowDimension();
Vector[] res = new Vector[M];
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
for (int i = 0; i < M; i++) {
res[i] = new DenseVector(AData[i]);
}
} else {
System.err.println("The input matrix should be a dense matrix.");
exit(1);
}
return res;
}
public static Vector[] denseMatrix2DenseColumnVectors(Matrix A) {
int N = A.getColumnDimension();
int M = A.getRowDimension();
Vector[] res = new Vector[N];
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
for (int j = 0; j < N; j++) {
double[] column = new double[M];
for (int i = 0; i < M; i++) {
column[i] = AData[i][j];
}
res[j] = new DenseVector(column);
}
} else {
System.err.println("The input matrix should be a dense matrix.");
exit(1);
}
return res;
}
/**
* Generate nRow by nCol all zero matrix.
*
* @param nRow number of rows
*
* @param nCol number of columns
*
* @return zeros(nRow, nCol)
*
*/
public static Matrix zeros(int nRow, int nCol) {
if (nRow == 0 || nCol == 0) {
return null;
}
return new DenseMatrix(nRow, nCol, 0);
}
/**
* Generate an all zero matrix with its size
* specified by a two dimensional integer array.
*
* @param size a two dimensional integer array
*
* @return an all zero matrix with its shape specified by size
*
*/
public static Matrix zeros(int[] size) {
if (size.length != 2) {
System.err.println("Input vector should have two elements!");
}
return zeros(size[0], size[1]);
}
/**
* Generate an n by n all zero matrix.
*
* @param n number of rows and columns
*
* @return ones(n)
*
*/
public static Matrix zeros(int n) {
return zeros(n, n);
}
/**
* Generate an all one matrix with nRow rows and nCol columns.
*
* @param nRow number of rows
*
* @param nCol number of columns
*
* @return ones(nRow, nCol)
*
*/
public static Matrix ones(int nRow, int nCol) {
if (nRow == 0 || nCol == 0) {
return null;
}
return new DenseMatrix(nRow, nCol, 1);
}
/**
* Generate an all one matrix with its size
* specified by a two dimensional integer array.
*
* @param size a two dimensional integer array
*
* @return an all one matrix with its shape specified by size
*
*/
public static Matrix ones(int[] size) {
if (size.length != 2) {
System.err.println("Input vector should have two elements!");
}
return ones(size[0], size[1]);
}
/**
* Generate an n by n all one matrix.
*
* @param n number of rows and columns
*
* @return ones(n)
*
*/
public static Matrix ones(int n) {
return ones(n, n);
}
/**
* Compute the determinant of a real square matrix.
*
* @param A a real square matrix
*
* @return det(A)
*/
public static double det(Matrix A) {
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (M != N) {
System.err.println("Input should be a square matrix.");
exit(1);
}
return new LUDecomposition(A).det();
}
/**
* Compute the inverse of a real square matrix.
*
* @param A a real square matrix
*
* @return inv(A)
*/
public static Matrix inv(Matrix A) {
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (M != N) {
System.err.println("Input should be a square matrix.");
exit(1);
}
LUDecomposition LUDecomp = new LUDecomposition(A);
if (LUDecomp.det() == 0) {
System.err.println("The input matrix is not invertible.");
exit(1);
}
return LUDecomp.inverse();
}
/**
* Sort elements of a vector V in place with an increasing order.
*
* @param V a real vector, it will be sorted in place.
*
* @return sorted indices represented by a 1D {@code double} array
*/
public static double[] sort(Vector V) {
return sort(V, "ascend");
}
/**
* Sort elements of a vector V in place with a specified order.
*
* @param V a real vector, it will be sorted in place.
*
* @param order a {@code String} variable either "ascend" or "descend"
*
* @return sorted indices represented by a 1D {@code double} array
*/
public static double[] sort(Vector V, String order) {
double[] indices = null;
int dim = V.getDim();
if (V instanceof DenseVector) {
double[] pr = ((DenseVector) V).getPr();
indices = new double[dim];
for (int k = 0; k < dim; k++) {
indices[k] = k;
}
quickSort(pr, indices, 0, dim - 1, order);
} else if (V instanceof SparseVector) {
double[] pr = ((SparseVector) V).getPr();
int[] ir = ((SparseVector) V).getIr();
int nnz = ((SparseVector) V).getNNZ();
indices = new double[dim];
int insertionPos = nnz;
if (order.equalsIgnoreCase("ascend")) {
for (int k = 0; k < nnz; k++) {
if (pr[k] >= 0) {
insertionPos--;
}
}
} else if (order.equalsIgnoreCase("descend")) {
for (int k = 0; k < nnz; k++) {
if (pr[k] <= 0) {
insertionPos--;
}
}
}
int lastIdx = -1;
int currentIdx = 0;
int cnt = insertionPos;
for (int k = 0; k < nnz; k++) {
currentIdx = ir[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
indices[cnt++] = idx;
}
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
indices[cnt++] = idx;
}
quickSort(pr, ir, 0, nnz - 1, order);
for (int k = 0; k < insertionPos; k++) {
indices[k] = ir[k];
}
for (int k = insertionPos; k < nnz; k++) {
indices[k + dim - nnz] = ir[k];
}
for (int k = 0; k < nnz; k++) {
if (k < insertionPos) {
ir[k] = k;
} else {
ir[k] = k + dim - nnz;
}
}
}
return indices;
}
/**
* Sort elements of a vector V in place with a specified order.
*
* @param V a real vector, it will be sorted in place.
*
* @param order a {@code String} variable either "ascend" or "descend"
*
* @return sorted indices represented by a 1D {@code double} array
*/
@Deprecated
public static double[] sort1(Vector V, String order) {
double[] indices = null;
int dim = V.getDim();
if (V instanceof DenseVector) {
double[] pr = ((DenseVector) V).getPr();
indices = new double[dim];
for (int k = 0; k < dim; k++) {
indices[k] = k;
}
quickSort(pr, indices, 0, dim - 1, order);
} else if (V instanceof SparseVector) {
double[] pr = ((SparseVector) V).getPr();
int[] ir = ((SparseVector) V).getIr();
int nnz = ((SparseVector) V).getNNZ();
int[] ir_ori = ir.clone();
quickSort(pr, ir, 0, nnz - 1, order);
int insertionPos = nnz;
if (order.equalsIgnoreCase("ascend")) {
for (int k = 0; k < nnz; k++) {
if (pr[k] >= 0) {
insertionPos = k;
break;
}
}
} else if (order.equalsIgnoreCase("descend")) {
for (int k = 0; k < nnz; k++) {
if (pr[k] <= 0) {
insertionPos = k;
break;
}
}
}
indices = new double[dim];
for (int k = 0; k < insertionPos; k++) {
indices[k] = ir[k];
}
int lastIdx = -1;
int currentIdx = 0;
int cnt = insertionPos;
for (int k = 0; k < nnz; k++) {
currentIdx = ir_ori[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
indices[cnt++] = idx;
}
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
indices[cnt++] = idx;
}
for (int k = insertionPos; k < nnz; k++) {
indices[k + dim - nnz] = ir[k];
}
for (int k = 0; k < nnz; k++) {
if (k < insertionPos) {
ir[k] = k;
} else {
ir[k] = k + dim - nnz;
}
}
}
return indices;
}
/**
* Sort elements of a vector V in place with a specified order.
*
* @param V a real vector, it will be sorted in place.
*
* @param order a {@code String} variable either "ascend" or "descend"
*
* @return sorted indices represented by a 1D {@code double} array
*/
@Deprecated
public static double[] sort0(Vector V, String order) {
double[] indices = null;
int dim = V.getDim();
if (V instanceof DenseVector) {
double[] pr = ((DenseVector) V).getPr();
indices = new double[dim];
for (int k = 0; k < dim; k++) {
indices[k] = k;
}
quickSort(pr, indices, 0, dim - 1, order);
} else if (V instanceof SparseVector) {
double[] pr = ((SparseVector) V).getPr();
int[] ir = ((SparseVector) V).getIr();
int nnz = ((SparseVector) V).getNNZ();
int[] ir_ori = ir.clone();
if (order.equalsIgnoreCase("ascend")) {
quickSort(pr, ir, 0, nnz - 1, order);
int numNegatives = nnz;
for (int k = 0; k < nnz; k++) {
if (pr[k] >= 0) {
numNegatives = k;
break;
}
}
indices = new double[dim];
for (int k = 0; k < numNegatives; k++) {
indices[k] = ir[k];
}
int lastIdx = -1;
int currentIdx = 0;
int cnt = numNegatives;
for (int k = 0; k < nnz; k++) {
currentIdx = ir_ori[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
indices[cnt++] = idx;
}
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
indices[cnt++] = idx;
}
/*for (int k = numNegatives; k < numNegatives + dim - nnz; k++) {
}*/
for (int k = numNegatives; k < nnz; k++) {
indices[k + dim - nnz] = ir[k];
}
for (int k = 0; k < nnz; k++) {
if (k < numNegatives) {
ir[k] = k;
} else {
ir[k] = k + dim - nnz;
}
}
} else if (order.equalsIgnoreCase("descend")) {
quickSort(pr, ir, 0, nnz - 1, order);
int numPositives = nnz;
for (int k = 0; k < nnz; k++) {
if (pr[k] <= 0) {
numPositives = k;
break;
}
}
indices = new double[dim];
for (int k = 0; k < numPositives; k++) {
indices[k] = ir[k];
}
int lastIdx = -1;
int currentIdx = 0;
int cnt = numPositives;
for (int k = 0; k < nnz; k++) {
currentIdx = ir_ori[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
indices[cnt++] = idx;
}
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
indices[cnt++] = idx;
}
/*for (int k = numNegatives; k < numNegatives + dim - nnz; k++) {
}*/
for (int k = numPositives; k < nnz; k++) {
indices[k + dim - nnz] = ir[k];
}
for (int k = 0; k < nnz; k++) {
if (k < numPositives) {
ir[k] = k;
} else {
ir[k] = k + dim - nnz;
}
}
}
}
return indices;
}
/**
* Sort elements of a matrix A on a direction in a specified order.
* A will not be modified.
*
* @param A a matrix to be sorted
*
* @param dim sorting direction, 1 for column-wise, 2 for row-wise
*
* @param order sorting order, either "ascend" or "descend"
*
* @return a {@code Matrix} array:
* res[0] is the sorted matrix
* res[1] is the sorted indices
*
*/
public static Matrix[] sort(Matrix A, int dim, String order) {
if (A == null) {
return null;
}
if (dim != 1 && dim != 2) {
System.err.println("Dimension should be either 1 or 2.");
exit(1);
}
if (!order.equalsIgnoreCase("ascend") && !order.equalsIgnoreCase("descend")) {
System.err.println("Order should be either \"ascend\" or \"descend\".");
exit(1);
}
Matrix[] res = new Matrix[2];
int M = A.getRowDimension();
int N = A.getColumnDimension();
Matrix sortedValues = null;
Matrix sortedIndices = null;
double[][] sortedIndexData = null;
if (A instanceof DenseMatrix) {
sortedValues = A.copy();
sortedIndices = null;
double[][] data = ((DenseMatrix) sortedValues).getData();
if (dim == 2) {
sortedIndexData = new double[M][];
double[] values = null;
double[] indices = null;
for (int i = 0; i < M; i++) {
values = data[i];
indices = new double[N];
for (int j = 0; j < N; j++) {
indices[j] = j;
}
quickSort(values, indices, 0, N - 1, order);
sortedIndexData[i] = indices;
}
sortedIndices = new DenseMatrix(sortedIndexData);
} else if (dim == 1) {
Matrix[] res2 = sort(A.transpose(), 2, order);
sortedValues = res2[0].transpose();
sortedIndices = res2[1].transpose();
}
} else if (A instanceof SparseMatrix) {
if (dim == 2) {
Vector[] rowVectors = sparseMatrix2SparseRowVectors(A);
sortedIndexData = new double[M][];
for (int i = 0; i < M; i++) {
sortedIndexData[i] = sort(rowVectors[i], order);
}
sortedIndices = new DenseMatrix(sortedIndexData);
sortedValues = sparseRowVectors2SparseMatrix(rowVectors);
} else if (dim == 1) {
Matrix[] res2 = sort(A.transpose(), 2, order);
sortedValues = res2[0].transpose();
sortedIndices = res2[1].transpose();
}
}
res[0] = sortedValues;
res[1] = sortedIndices;
return res;
}
/**
* Sort elements of a matrix A on a direction in an increasing order.
*
* @param A a matrix to be sorted
*
* @param dim sorting direction, 1 for column-wise, 2 for row-wise
*
* @return a {@code Matrix} array:
* res[0] is the sorted matrix
* res[1] is the sorted indices
*
*/
public static Matrix[] sort(Matrix A, int dim) {
return sort(A, dim, "ascend");
}
/**
* Sort elements of a matrix A by columns in a specified order.
* A will not be modified.
*
* @param A a matrix to be sorted
*
* @param order sorting order, either "ascend" or "descend"
*
* @return a {@code Matrix} array:
* res[0] is the sorted matrix
* res[1] is the sorted indices
*
*/
public static Matrix[] sort(Matrix A, String order) {
return sort(A, 1, order);
}
/**
* Sort elements of a matrix A by columns in an increasing order.
* A will not be modified.
*
* @param A a matrix to be sorted
*
* @return a {@code Matrix} array:
* res[0] is the sorted matrix
* res[1] is the sorted indices
*
*/
public static Matrix[] sort(Matrix A) {
return sort(A, "ascend");
}
/**
* Get a two dimensional integer array for size of a matrix A.
*
* @param A a matrix
*
* @return size(A)
*
*/
public static int[] size(Matrix A) {
return new int[] { A.getRowDimension(), A.getColumnDimension()};
}
/**
* Get the dimensionality on dim-th dimension for a matrix A.
*
* @param A a matrix
*
* @param dim dimension order
*
* @return size(A, dim)
*
*/
public static int size(Matrix A, int dim) {
if (dim == 1) {
return A.getRowDimension();
} else if (dim == 2) {
return A.getColumnDimension();
} else {
System.err.println("Dim error!");
return 0;
}
}
/**
* Compute maximum between elements of A and a real number and return
* as a matrix with the same shape of A.
*
* @param A a real matrix
*
* @param v a real number
*
* @return max(A, v)
*
*/
public static Matrix max(Matrix A, double v) {
if (A == null)
return null;
Matrix res = null;
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (A instanceof DenseMatrix) {
res = A.copy();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] < v) {
resRow[j] = v;
}
}
}
} else if (A instanceof SparseMatrix) {
if (v != 0) {
res = new DenseMatrix(M, N, v);
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
int[] ic = ((SparseMatrix) A).getIc();
int[] jr = ((SparseMatrix) A).getJr();
int[] valCSRIndices = ((SparseMatrix) A).getValCSRIndices();
double[] pr = ((SparseMatrix) A).getPr();
for (int i = 0; i < M; i++) {
resRow = resData[i];
if (jr[i] == jr[i + 1]) {
if (v < 0)
ArrayOperator.assignVector(resRow, 0);
continue;
}
int lastColumnIdx = -1;
int currentColumnIdx = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
currentColumnIdx = ic[k];
if (v < 0)
for (int c = lastColumnIdx + 1; c < currentColumnIdx; c++) {
resRow[c] = 0;
}
resRow[currentColumnIdx] = Math.max(pr[valCSRIndices[k]], v);
lastColumnIdx = currentColumnIdx;
}
for (int c = lastColumnIdx + 1; c < N; c++) {
if (v < 0)
resRow[c] = 0;
}
}
} else { // v == 0
return max(A, new SparseMatrix(M, N));
}
}
return res;
}
/**
* Compute maximum between elements of A and a real number and return
* as a matrix with the same shape of A.
*
* @param v a real number
*
* @param A a real matrix
*
* @return max(v, A)
*
*/
public static Matrix max(double v, Matrix A) {
return max(A, v);
}
/**
* Compute maximum between two real matrices A and B.
*
* @param A a real matrix
*
* @param B a real matrix
*
* @return max(A, B);
*/
public static Matrix max(Matrix A, Matrix B) {
Matrix res = null;
int M = A.getRowDimension();
int N = A.getColumnDimension();
double v = 0;
if (A instanceof DenseMatrix) {
res = A.copy();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
if (B instanceof DenseMatrix) {
double[][] BData = ((DenseMatrix) B).getData();
double[] BRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
BRow = BData[i];
resRow = resData[i];
for (int j = 0; j < N; j++) {
v = BRow[j];
if (resRow[j] < v)
resRow[j] = v;
}
}
} else if (B instanceof SparseMatrix) {
int[] ic = ((SparseMatrix) B).getIc();
int[] jr = ((SparseMatrix) B).getJr();
int[] valCSRIndices = ((SparseMatrix) B).getValCSRIndices();
double[] pr = ((SparseMatrix) B).getPr();
for (int i = 0; i < M; i++) {
resRow = resData[i];
if (jr[i] == jr[i + 1]) {
for (int j = 0; j < N; j++) {
if (resRow[j] < 0)
resRow[j] = 0;
}
continue;
}
int lastColumnIdx = -1;
int currentColumnIdx = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
currentColumnIdx = ic[k];
for (int c = lastColumnIdx + 1; c < currentColumnIdx; c++) {
if (resRow[c] < 0)
resRow[c] = 0;
}
resRow[currentColumnIdx] = Math.max(pr[valCSRIndices[k]], resRow[currentColumnIdx]);
lastColumnIdx = currentColumnIdx;
}
for (int c = lastColumnIdx + 1; c < N; c++) {
if (resRow[c] < 0)
resRow[c] = 0;
}
}
}
} else if (A instanceof SparseMatrix) {
if (B instanceof DenseMatrix) {
return max(B, A);
} else if (B instanceof SparseMatrix) {
res = new SparseMatrix(M, N);
int[] ir1 = null;
int[] jc1 = null;
double[] pr1 = null;
ir1 = ((SparseMatrix) A).getIr();
jc1 = ((SparseMatrix) A).getJc();
pr1 = ((SparseMatrix) A).getPr();
int[] ir2 = null;
int[] jc2 = null;
double[] pr2 = null;
ir2 = ((SparseMatrix) B).getIr();
jc2 = ((SparseMatrix) B).getJc();
pr2 = ((SparseMatrix) B).getPr();
int k1 = 0;
int k2 = 0;
int r1 = -1;
int r2 = -1;
int i = -1;
// double v = 0;
for (int j = 0; j < N; j++) {
k1 = jc1[j];
k2 = jc2[j];
// Both A and B's j-th columns are empty.
if (k1 == jc1[j + 1] && k2 == jc2[j + 1])
continue;
while (k1 < jc1[j + 1] || k2 < jc2[j + 1]) {
if (k2 == jc2[j + 1]) { // B's j-th column has been processed.
i = ir1[k1];
v = pr1[k1];
if (v < 0)
v = 0;
k1++;
} else if (k1 == jc1[j + 1]) { // A's j-th column has been processed.
i = ir2[k2];
v = pr2[k2];
if (v < 0)
v = 0;
k2++;
} else { // Both A and B's j-th columns have not been fully processed.
r1 = ir1[k1];
r2 = ir2[k2];
if (r1 < r2) {
i = r1;
v = pr1[k1];
if (v < 0)
v = 0;
k1++;
} else if (r1 == r2) {
i = r1;
v = Math.max(pr1[k1], pr2[k2]);
k1++;
k2++;
} else {
i = r2;
v = pr2[k2];
if (v < 0)
v = 0;
k2++;
}
}
if (v != 0)
res.setEntry(i, j, v);
}
}
}
}
return res;
}
/**
* Compute minimum between elements of A and a real number and return
* as a matrix with the same shape of A.
*
* @param A a real matrix
*
* @param v a real number
*
* @return min(A, v)
*
*/
public static Matrix min(Matrix A, double v) {
if (A == null)
return null;
Matrix res = null;
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (A instanceof DenseMatrix) {
res = A.copy();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] > v) {
resRow[j] = v;
}
}
}
} else if (A instanceof SparseMatrix) {
if (v != 0) {
res = new DenseMatrix(M, N, v);
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
int[] ic = ((SparseMatrix) A).getIc();
int[] jr = ((SparseMatrix) A).getJr();
int[] valCSRIndices = ((SparseMatrix) A).getValCSRIndices();
double[] pr = ((SparseMatrix) A).getPr();
for (int i = 0; i < M; i++) {
resRow = resData[i];
if (jr[i] == jr[i + 1]) {
if (v > 0)
ArrayOperator.assignVector(resRow, 0);
continue;
}
int lastColumnIdx = -1;
int currentColumnIdx = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
currentColumnIdx = ic[k];
if (v > 0)
for (int c = lastColumnIdx + 1; c < currentColumnIdx; c++) {
resRow[c] = 0;
}
resRow[currentColumnIdx] = Math.min(pr[valCSRIndices[k]], v);
lastColumnIdx = currentColumnIdx;
}
for (int c = lastColumnIdx + 1; c < N; c++) {
if (v > 0)
resRow[c] = 0;
}
}
} else { // v == 0
return min(A, new SparseMatrix(M, N));
}
}
return res;
}
/**
* Compute minimum between elements of A and a real number and return
* as a matrix with the same shape of A.
*
* @param v a real number
*
* @param A a real matrix
*
* @return min(v, A)
*
*/
public static Matrix min(double v, Matrix A) {
return min(A, v);
}
/**
* Compute minimum between two real matrices A and B.
*
* @param A a real matrix
*
* @param B a real matrix
*
* @return max(A, B);
*/
public static Matrix min(Matrix A, Matrix B) {
Matrix res = null;
int M = A.getRowDimension();
int N = A.getColumnDimension();
double v = 0;
if (A instanceof DenseMatrix) {
res = A.copy();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
if (B instanceof DenseMatrix) {
double[][] BData = ((DenseMatrix) B).getData();
double[] BRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
BRow = BData[i];
resRow = resData[i];
for (int j = 0; j < N; j++) {
v = BRow[j];
if (resRow[j] > v)
resRow[j] = v;
}
}
} else if (B instanceof SparseMatrix) {
int[] ic = ((SparseMatrix) B).getIc();
int[] jr = ((SparseMatrix) B).getJr();
int[] valCSRIndices = ((SparseMatrix) B).getValCSRIndices();
double[] pr = ((SparseMatrix) B).getPr();
for (int i = 0; i < M; i++) {
resRow = resData[i];
if (jr[i] == jr[i + 1]) {
for (int j = 0; j < N; j++) {
if (resRow[j] > 0)
resRow[j] = 0;
}
continue;
}
int lastColumnIdx = -1;
int currentColumnIdx = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
currentColumnIdx = ic[k];
for (int c = lastColumnIdx + 1; c < currentColumnIdx; c++) {
if (resRow[c] > 0)
resRow[c] = 0;
}
resRow[currentColumnIdx] = Math.min(pr[valCSRIndices[k]], resRow[currentColumnIdx]);
lastColumnIdx = currentColumnIdx;
}
for (int c = lastColumnIdx + 1; c < N; c++) {
if (resRow[c] > 0)
resRow[c] = 0;
}
}
}
} else if (A instanceof SparseMatrix) {
if (B instanceof DenseMatrix) {
return min(B, A);
} else if (B instanceof SparseMatrix) {
res = new SparseMatrix(M, N);
int[] ir1 = null;
int[] jc1 = null;
double[] pr1 = null;
ir1 = ((SparseMatrix) A).getIr();
jc1 = ((SparseMatrix) A).getJc();
pr1 = ((SparseMatrix) A).getPr();
int[] ir2 = null;
int[] jc2 = null;
double[] pr2 = null;
ir2 = ((SparseMatrix) B).getIr();
jc2 = ((SparseMatrix) B).getJc();
pr2 = ((SparseMatrix) B).getPr();
int k1 = 0;
int k2 = 0;
int r1 = -1;
int r2 = -1;
int i = -1;
// double v = 0;
for (int j = 0; j < N; j++) {
k1 = jc1[j];
k2 = jc2[j];
// Both A and B's j-th columns are empty.
if (k1 == jc1[j + 1] && k2 == jc2[j + 1])
continue;
while (k1 < jc1[j + 1] || k2 < jc2[j + 1]) {
if (k2 == jc2[j + 1]) { // B's j-th column has been processed.
i = ir1[k1];
v = pr1[k1];
if (v > 0)
v = 0;
k1++;
} else if (k1 == jc1[j + 1]) { // A's j-th column has been processed.
i = ir2[k2];
v = pr2[k2];
if (v > 0)
v = 0;
k2++;
} else { // Both A and B's j-th columns have not been fully processed.
r1 = ir1[k1];
r2 = ir2[k2];
if (r1 < r2) {
i = r1;
v = pr1[k1];
if (v > 0)
v = 0;
k1++;
} else if (r1 == r2) {
i = r1;
v = Math.min(pr1[k1], pr2[k2]);
k1++;
k2++;
} else {
i = r2;
v = pr2[k2];
if (v > 0)
v = 0;
k2++;
}
}
if (v != 0)
res.setEntry(i, j, v);
}
}
}
}
return res;
}
/**
* Compute maximum between elements of V and a real number
* and return as a vector with the same shape of V.
*
* @param V a real vector
*
* @param v a real number
*
* @return max(V, v)
*
*/
public static Vector max(Vector V, double v) {
Vector res = null;
int dim = V.getDim();
if (V instanceof DenseVector) {
res = V.copy();
double[] pr = ((DenseVector) res).getPr();
for (int k = 0; k < dim; k++) {
if (pr[k] < v)
pr[k] = v;
}
} else if (V instanceof SparseVector) {
if (v != 0) {
res = new DenseVector(dim, v);
double[] prRes = ((DenseVector) res).getPr();
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int lastIdx = -1;
int currentIdx = -1;
for (int k = 0; k < nnz; k++) {
currentIdx = ir[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
if (v < 0)
prRes[idx] = 0;
}
if (v < pr[k])
prRes[currentIdx] = pr[k];
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
if (v < 0)
prRes[idx] = 0;
}
} else { // v == 0
res = V.copy();
double[] pr = ((SparseVector) res).getPr();
int nnz = ((SparseVector) res).getNNZ();
for (int k = 0; k < nnz; k++) {
if (pr[k] < 0)
pr[k] = 0;
}
((SparseVector) res).clean();
}
}
return res;
}
/**
* Compute maximum between elements of V and a real number
* and return as a vector with the same shape of V.
*
* @param v a real number
*
* @param V a real vector
*
* @return max(v, V)
*
*/
public static Vector max(double v, Vector V) {
return max(V, v);
}
/**
* Compute maximum between two vectors U and V.
*
* @param U a real vector
*
* @param V a real vector
*
* @return max(U, V)
*
*/
public static Vector max(Vector U, Vector V) {
Vector res = null;
int dim = U.getDim();
double v = 0;
if (U instanceof DenseVector) {
res = U.copy();
double[] prRes = ((DenseVector) res).getPr();
if (V instanceof DenseVector) {
// double[] prU = ((DenseVector) U).getPr();
double[] prV = ((DenseVector) V).getPr();
for (int k = 0; k < dim; k++) {
v = prV[k];
if (prRes[k] < v)
prRes[k] = v;
}
} else if (V instanceof SparseVector) {
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int lastIdx = -1;
int currentIdx = -1;
for (int k = 0; k < nnz; k++) {
currentIdx = ir[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
if (prRes[idx] < 0)
prRes[idx] = 0;
}
v = pr[k];
if (prRes[currentIdx] < v)
prRes[currentIdx] = v;
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
if (prRes[idx] < 0)
prRes[idx] = 0;
}
}
} else if (U instanceof SparseVector) {
if (V instanceof DenseVector) {
return max(V, U);
} else if (V instanceof SparseVector) {
res = new SparseVector(dim);
int[] ir1 = ((SparseVector) V).getIr();
double[] pr1 = ((SparseVector) V).getPr();
int nnz1 = ((SparseVector) V).getNNZ();
int[] ir2 = ((SparseVector) U).getIr();
double[] pr2 = ((SparseVector) U).getPr();
int nnz2 = ((SparseVector) U).getNNZ();
if (!(nnz1 == 0 && nnz2 == 0)) {
int k1 = 0;
int k2 = 0;
int r1 = 0;
int r2 = 0;
// double v = 0;
int i = -1;
while (k1 < nnz1 || k2 < nnz2) {
if (k2 == nnz2) { // V has been processed.
i = ir1[k1];
v = pr1[k1];
if (v < 0)
v = 0;
k1++;
} else if (k1 == nnz1) { // this has been processed.
i = ir2[k2];
v = pr2[k2];
if (v < 0)
v = 0;
k2++;
} else { // Both this and V have not been fully processed.
r1 = ir1[k1];
r2 = ir2[k2];
if (r1 < r2) {
i = r1;
v = pr1[k1];
if (v < 0)
v = 0;
k1++;
} else if (r1 == r2) {
i = r1;
v = Math.max(pr1[k1], pr2[k2]);
k1++;
k2++;
} else {
i = r2;
v = pr2[k2];
if (v < 0)
v = 0;
k2++;
}
}
if (v != 0) {
res.set(i, v);
}
}
}
}
}
return res;
}
/**
* Compute minimum between elements of V and a real number
* and return as a vector with the same shape of V.
*
* @param V a real vector
*
* @param v a real number
*
* @return min(V, v)
*
*/
public static Vector min(Vector V, double v) {
Vector res = null;
int dim = V.getDim();
if (V instanceof DenseVector) {
res = V.copy();
double[] pr = ((DenseVector) res).getPr();
for (int k = 0; k < dim; k++) {
if (pr[k] > v)
pr[k] = v;
}
} else if (V instanceof SparseVector) {
if (v != 0) {
res = new DenseVector(dim, v);
double[] prRes = ((DenseVector) res).getPr();
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int lastIdx = -1;
int currentIdx = -1;
for (int k = 0; k < nnz; k++) {
currentIdx = ir[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
if (v > 0)
prRes[idx] = 0;
}
if (v > pr[k])
prRes[currentIdx] = pr[k];
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
if (v > 0)
prRes[idx] = 0;
}
} else { // v == 0
res = V.copy();
double[] pr = ((SparseVector) res).getPr();
int nnz = ((SparseVector) res).getNNZ();
for (int k = 0; k < nnz; k++) {
if (pr[k] > 0)
pr[k] = 0;
}
((SparseVector) res).clean();
}
}
return res;
}
/**
* Compute minimum between elements of V and a real number
* and return as a vector with the same shape of V.
*
* @param v a real number
*
* @param V a real vector
*
* @return min(v, V)
*
*/
public static Vector min(double v, Vector V) {
return min(V, v);
}
/**
* Compute minimum between two vectors U and V.
*
* @param U a real vector
*
* @param V a real vector
*
* @return min(U, V)
*
*/
public static Vector min(Vector U, Vector V) {
Vector res = null;
int dim = U.getDim();
double v = 0;
if (U instanceof DenseVector) {
res = U.copy();
double[] prRes = ((DenseVector) res).getPr();
if (V instanceof DenseVector) {
// double[] prU = ((DenseVector) U).getPr();
double[] prV = ((DenseVector) V).getPr();
for (int k = 0; k < dim; k++) {
v = prV[k];
if (prRes[k] > v)
prRes[k] = v;
}
} else if (V instanceof SparseVector) {
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int lastIdx = -1;
int currentIdx = -1;
for (int k = 0; k < nnz; k++) {
currentIdx = ir[k];
for (int idx = lastIdx + 1; idx < currentIdx; idx++) {
if (prRes[idx] > 0)
prRes[idx] = 0;
}
v = pr[k];
if (prRes[currentIdx] > v)
prRes[currentIdx] = v;
lastIdx = currentIdx;
}
for (int idx = lastIdx + 1; idx < dim; idx++) {
if (prRes[idx] > 0)
prRes[idx] = 0;
}
}
} else if (U instanceof SparseVector) {
if (V instanceof DenseVector) {
return max(V, U);
} else if (V instanceof SparseVector) {
res = new SparseVector(dim);
int[] ir1 = ((SparseVector) V).getIr();
double[] pr1 = ((SparseVector) V).getPr();
int nnz1 = ((SparseVector) V).getNNZ();
int[] ir2 = ((SparseVector) U).getIr();
double[] pr2 = ((SparseVector) U).getPr();
int nnz2 = ((SparseVector) U).getNNZ();
if (!(nnz1 == 0 && nnz2 == 0)) {
int k1 = 0;
int k2 = 0;
int r1 = 0;
int r2 = 0;
// double v = 0;
int i = -1;
while (k1 < nnz1 || k2 < nnz2) {
if (k2 == nnz2) { // V has been processed.
i = ir1[k1];
v = pr1[k1];
if (v > 0)
v = 0;
k1++;
} else if (k1 == nnz1) { // this has been processed.
i = ir2[k2];
v = pr2[k2];
if (v > 0)
v = 0;
k2++;
} else { // Both this and V have not been fully processed.
r1 = ir1[k1];
r2 = ir2[k2];
if (r1 < r2) {
i = r1;
v = pr1[k1];
if (v > 0)
v = 0;
k1++;
} else if (r1 == r2) {
i = r1;
v = Math.min(pr1[k1], pr2[k2]);
k1++;
k2++;
} else {
i = r2;
v = pr2[k2];
if (v > 0)
v = 0;
k2++;
}
}
if (v != 0) {
res.set(i, v);
}
}
}
}
}
return res;
}
/**
* Logical indexing A by B for the syntax A(B). A logical matrix
* provides a different type of array indexing in MATLAB. While
* most indices are numeric, indicating a certain row or column
* number, logical indices are positional. That is, it is the
* position of each 1 in the logical matrix that determines which
* array element is being referred to.
*
* @param A a real matrix
*
* @param B a logical matrix with elements being either 1 or 0
*
* @return A(B)
*
*/
public static Matrix logicalIndexing(Matrix A, Matrix B) {
int nA = A.getColumnDimension();
int dA = A.getRowDimension();
int nB = B.getColumnDimension();
int dB = B.getRowDimension();
if (nA != nB || dA != dB) {
System.err.println("The input matrices should have same size!");
System.exit(1);
}
ArrayList<Double> vals = new ArrayList<Double>();
double b;
for (int j = 0; j < nA; j++) {
for (int i = 0; i < dA; i++) {
b = B.getEntry(i, j);
if (b == 1)
vals.add(A.getEntry(i, j));
else if (b != 0)
System.err.println("Elements of the logical matrix should be either 1 or 0!");
}
}
Double[] Data = new Double[vals.size()];
vals.toArray(Data);
double[] data = new double[vals.size()];
for (int i = 0; i < vals.size(); i++) {
data[i] = Data[i];
}
if (data.length != 0)
return new DenseMatrix(data, 1);
else
return null;
}
/**
* Linear indexing V by an index array.
*
* @param V an {@code int} array
*
* @param indices an {@code int} array of selected indices
*
* @return V(indices)
*
*/
public static int[] linearIndexing(int[] V, int[] indices) {
if (indices == null || indices.length == 0) {
return null;
}
int[] res = new int[indices.length];
for (int i = 0; i < indices.length; i++) {
res[i] = V[indices[i]];
}
return res;
}
/**
* Linear indexing V by an index array.
*
* @param V an {@code double} array
*
* @param indices an {@code int} array of selected indices
*
* @return V(indices)
*
*/
public static double[] linearIndexing(double[] V, int[] indices) {
if (indices == null || indices.length == 0) {
return null;
}
double[] res = new double[indices.length];
for (int i = 0; i < indices.length; i++) {
res[i] = V[indices[i]];
}
return res;
}
/**
* Linear indexing A by an index array.
*
* @param A a real matrix
*
* @param indices an {@code int} array of selected indices
*
* @return A(indices)
*
*/
public static Matrix linearIndexing(Matrix A, int[] indices) {
if (indices == null || indices.length == 0) {
return null;
}
Matrix res = null;
if (A instanceof DenseMatrix) {
res = new DenseMatrix(indices.length, 1);
} else {
res = new SparseMatrix(indices.length, 1);
}
int nRow = A.getRowDimension();
// int nCol = A.getColumnDimension();
int r = -1;
int c = -1;
int index = -1;
for (int i = 0; i < indices.length; i++) {
index = indices[i];
r = index % nRow;
c = index / nRow;
res.setEntry(i, 0, A.getEntry(r, c));
}
return res;
}
/**
* Matrix assignment by linear indexing for the syntax A(B) = V.
*
* @param A a matrix to be assigned
*
* @param idx a linear index vector
*
* @param V a column matrix to assign A(idx)
*
*/
public static void linearIndexingAssignment(Matrix A, int[] idx, Matrix V) {
if (V == null)
return;
int nV = V.getColumnDimension();
int dV = V.getRowDimension();
if (nV != 1)
System.err.println("Assignment matrix should be a column matrix!");
if (idx.length != dV)
System.err.println("Assignment with different number of elements!");
int nRow = A.getRowDimension();
// int nCol = A.getColumnDimension();
int r = -1;
int c = -1;
int index = -1;
for (int i = 0; i < idx.length; i++) {
index = idx[i];
r = index % nRow;
c = index / nRow;
A.setEntry(r, c, V.getEntry(i, 0));
}
}
/**
* Matrix assignment by linear indexing for the syntax A(B) = v.
*
* @param A a matrix to be assigned
*
* @param idx a linear index vector
*
* @param v a real scalar to assign A(idx)
*
*/
public static void linearIndexingAssignment(Matrix A, int[] idx, double v) {
int nRow = A.getRowDimension();
// int nCol = A.getColumnDimension();
int r = -1;
int c = -1;
int index = -1;
for (int i = 0; i < idx.length; i++) {
index = idx[i];
r = index % nRow;
c = index / nRow;
A.setEntry(r, c, v);
}
}
/**
* Matrix assignment by logical indexing for the syntax A(B) = v.
*
* @param A a matrix to be assigned
*
* @param B a logical matrix where position of each 1 determines
* which array element is being assigned
*
* @param v a real scalar to assign A(B)
*
*/
public static void logicalIndexingAssignment(Matrix A, Matrix B, double v) {
int nA = A.getColumnDimension();
int dA = A.getRowDimension();
int nB = B.getColumnDimension();
int dB = B.getRowDimension();
if (nA != nB || dA != dB) {
System.err.println("The input matrices for logical indexing should have same size!");
System.exit(1);
}
double b;
if (B instanceof SparseMatrix) {
int[] ir = ((SparseMatrix) B).getIr();
int[] jc = ((SparseMatrix) B).getJc();
double[] pr = ((SparseMatrix) B).getPr();
for (int j = 0; j < nB; j++) {
for (int k = jc[j]; k < jc[j + 1]; k++) {
b = pr[k];
if (b == 1)
A.setEntry(ir[k], j, v);
else if (b != 0)
err("Elements of the logical matrix should be either 1 or 0!");
}
}
} else if (B instanceof DenseMatrix) {
double[][] BData = ((DenseMatrix) B).getData();
double[] BRow = null;
for (int i = 0; i < dA; i++) {
BRow = BData[i];
for (int j = 0; j < nA; j++) {
b = BRow[j];
if (b == 1) {
A.setEntry(i, j, v);
}
else if (b != 0)
System.err.println("Elements of the logical matrix should be either 1 or 0!");
}
}
}
}
/**
* Matrix assignment by logical indexing for the syntax A(B) = V.
*
* @param A a matrix to be assigned
*
* @param B a logical matrix where position of each 1 determines
* which array element is being assigned
*
* @param V a column matrix to assign A(B)
*
*/
public static void logicalIndexingAssignment(Matrix A, Matrix B, Matrix V) {
int nA = A.getColumnDimension();
int dA = A.getRowDimension();
int nB = B.getColumnDimension();
int dB = B.getRowDimension();
if (nA != nB || dA != dB) {
System.err.println("The input matrices for logical indexing should have same size!");
System.exit(1);
}
if (V == null)
return;
int nV = V.getColumnDimension();
int dV = V.getRowDimension();
if (nV != 1) {
System.err.println("Assignment matrix should be a column matrix!");
exit(1);
}
/*double b;
int cnt = 0;
for (int j = 0; j < nA; j++) {
for (int i = 0; i < dA; i++) {
b = B.getEntry(i, j);
if (b == 1) {
A.setEntry(i, j, V.getEntry(cnt++, 0));
}
else if (b != 0)
System.err.println("Elements of the logical matrix should be either 1 or 0!");
}
}*/
double b;
int cnt = 0;
if (B instanceof SparseMatrix) {
int[] ir = ((SparseMatrix) B).getIr();
int[] jc = ((SparseMatrix) B).getJc();
double[] pr = ((SparseMatrix) B).getPr();
for (int j = 0; j < nB; j++) {
for (int k = jc[j]; k < jc[j + 1]; k++) {
b = pr[k];
if (b == 1)
A.setEntry(ir[k], j, V.getEntry(cnt++, 0));
else if (b != 0)
err("Elements of the logical matrix should be either 1 or 0!");
}
}
} else if (B instanceof DenseMatrix) {
double[][] BData = ((DenseMatrix) B).getData();
for (int j = 0; j < nA; j++) {
for (int i = 0; i < dA; i++) {
b = BData[i][j];
if (b == 1) {
A.setEntry(i, j, V.getEntry(cnt++, 0));
}
else if (b != 0)
System.err.println("Elements of the logical matrix should be either 1 or 0!");
}
}
}
if (cnt != dV)
System.err.println("Assignment with different number of elements!");
}
/**
* Get the nonnegative part of a matrix A.
*
* @param A a matrix
*
* @return a matrix which is the nonnegative part of a matrix A
*
*/
public static Matrix subplus(Matrix A) {
Matrix res = A.copy();
int M = res.getRowDimension();
int N = res.getColumnDimension();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] < 0)
resRow[j] = 0;
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
int nnz = ((SparseMatrix) res).getNNZ();
for (int k = 0; k < nnz; k++) {
if (pr[k] < 0)
pr[k] = 0;
}
((SparseMatrix) res).clean();
}
return res;
}
/**
* Round towards zero.
*
* @param x a real number
*
* @return fix(x)
*
*/
public static int fix(double x) {
if (x > 0) {
return (int)Math.floor(x);
} else {
return (int)Math.ceil(x);
}
}
/**
* Generates a linearly spaced integer array with distance of
* D between two consecutive numbers. colon(J, D, K) is
* the same as [J, J+D, ..., J+m*D] where m = fix((K-J)/D).
*
* @param begin starting point (inclusive)
*
* @param d distance between two consecutive numbers
*
* @param end ending point (inclusive if possible)
*
* @return indices array for the syntax begin:d:end
*
*//*
public static int[] colon(int begin, int d, int end) {
int m = fix((end - begin) / d);
if (m < 0) {
System.err.println("Difference error!");
System.exit(1);
}
int[] res = new int[m + 1];
for (int i = 0; i <= m; i++) {
res[i] = begin + i * d;
}
return res;
}
*//**
* Same as colon(begin, 1, end).
*
* @param begin starting point (inclusive)
*
* @param end ending point (inclusive)
*
* @return indices array for the syntax begin:end
*
*//*
public static int[] colon(int begin, int end) {
return colon(begin, 1, end);
}*/
/**
* Returns an array that contains 1's where
* the elements of X are NaN's and 0's where they are not.
*
* @param A a matrix
*
* @return a 0-1 matrix: isnan(A)
*
*/
public static Matrix isnan(Matrix A) {
Matrix res = A.copy();
int M = res.getRowDimension();
int N = res.getColumnDimension();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (Double.isNaN(resRow[j]))
resRow[j] = 1;
else
resRow[j] = 0;
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
int nnz = ((SparseMatrix) res).getNNZ();
for (int k = 0; k < nnz; k++) {
if (Double.isNaN(pr[k]))
pr[k] = 1;
else
pr[k] = 0;
}
((SparseMatrix) res).clean();
}
return res;
}
/**
* returns an array that contains 1's where the
* elements of X are +Inf or -Inf and 0's where they are not.
*
* @param A a real matrix
*
* @return a 0-1 matrix: isinf(A)
*
*/
public static Matrix isinf(Matrix A) {
Matrix res = A.copy();
int M = res.getRowDimension();
int N = res.getColumnDimension();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (Double.isInfinite(resRow[j]))
resRow[j] = 1;
else
resRow[j] = 0;
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
int nnz = ((SparseMatrix) res).getNNZ();
for (int k = 0; k < nnz; k++) {
if (Double.isInfinite(pr[k]))
pr[k] = 1;
else
pr[k] = 0;
}
((SparseMatrix) res).clean();
}
return res;
}
/**
* Do element by element comparisons between A and B and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
* A and B must have the same dimensions.
*
* @param A a matrix
*
* @param B a matrix
*
* @return A | B or or(A, B)
*
*/
public static Matrix or(Matrix A, Matrix B) {
Matrix res = A.plus(B);
int M = res.getRowDimension();
int N = res.getColumnDimension();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] > 1)
resRow[j] = 1;
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
int nnz = ((SparseMatrix) res).getNNZ();
for (int k = 0; k < nnz; k++) {
if (pr[k] > 1)
pr[k] = 1;
}
}
return res;
}
/**
* Do element by element comparisons between A and B and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
* A and B must have the same dimensions.
*
* @param A a matrix
*
* @param B a matrix
*
* @return A & B or and(A, B)
*
*/
public static Matrix and(Matrix A, Matrix B) {
Matrix res = A.times(B);
return res;
}
/**
* Performs a logical NOT of input array A, and returns an array
* containing elements set to either 1 (TRUE) or 0 (FALSE). An
* element of the output array is set to 1 if A contains a zero
* value element at that same array location. Otherwise, that
* element is set to 0.
*
* @param A a matrix
*
* @return ~A or not(A)
*
*/
public static Matrix not(Matrix A) {
Matrix res = minus(1, A);
return res;
}
/**
* Do element by element comparisons between X and Y and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param Y a matrix
*
* @return X ~= Y or ne(X, Y)
*
*/
public static Matrix ne(Matrix X, Matrix Y) {
Matrix res = X.minus(Y);
int M = res.getRowDimension();
int N = res.getColumnDimension();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] != 0)
resRow[j] = 1;
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
int nnz = ((SparseMatrix) res).getNNZ();
for (int k = 0; k < nnz; k++) {
if (pr[k] != 0)
pr[k] = 1;
}
}
return res;
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param x a real scalar
*
* @return X ~= x or ne(X, x)
*
*/
public static Matrix ne(Matrix X, double x) {
Matrix res = X.minus(x);
// res must be a dense matrix.
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] != 0)
resRow[j] = 1;
}
}
return res;
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param x a real scalar
*
* @param X a matrix
*
* @return X ~= x or ne(X, x)
*
*/
public static Matrix ne(double x, Matrix X) {
Matrix res = X.minus(x);
// res must be a dense matrix.
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] != 0)
resRow[j] = 1;
}
}
return res;
}
/**
* Do element by element comparisons between X and Y and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param Y a matrix
*
* @return X == Y or eq(X, Y)
*
*/
public static Matrix eq(Matrix X, Matrix Y) {
return minus(1, ne(X, Y));
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param x a real scalar
*
* @return X == x or eq(X, x)
*
*/
public static Matrix eq(Matrix X, double x) {
Matrix res = X.minus(x);
// res must be a dense matrix.
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] != 0)
resRow[j] = 0;
else
resRow[j] = 1;
}
}
return res;
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param x a real scalar
*
* @param X a matrix
*
* @return X == x or eq(X, x)
*
*/
public static Matrix eq(double x, Matrix X) {
return eq(X, x);
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param x a real scalar
*
* @return X >= x or ge(X, x)
*
*/
public static Matrix ge(Matrix X, double x) {
Matrix res = X.minus(x);
// res must be a dense matrix.
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] >= 0)
resRow[j] = 1;
else
resRow[j] = 0;
}
}
return res;
}
/**
* Do element by element comparisons between x and X and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param x a real scalar
*
* @param X a matrix
*
* @return x >= X or ge(x, X)
*/
public static Matrix ge(double x, Matrix X) {
Matrix res = minus(x, X);
// res must be a dense matrix.
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] >= 0)
resRow[j] = 1;
else
resRow[j] = 0;
}
}
return res;
}
/**
* Do element by element comparisons between X and Y and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
* X and Y must have the same dimensions.
*
* @param X a matrix
*
* @param Y a matrix
*
* @return X >= Y or ge(X, Y)
*
*/
public static Matrix ge(Matrix X, Matrix Y) {
Matrix res = full(X.minus(Y));
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] >= 0)
resRow[j] = 1;
else
resRow[j] = 0;
}
}
return res;
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param x a real scalar
*
* @return X <= x or le(X, x)
*
*/
public static Matrix le(Matrix X, double x) {
return ge(x, X);
}
/**
* Do element by element comparisons between x and X and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param x a real scalar
*
* @param X a matrix
*
* @return x <= X or le(x, X)
*
*/
public static Matrix le(double x, Matrix X) {
return ge(X, x);
}
/**
* Do element by element comparisons between X and Y and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
* X and Y must have the same dimensions.
*
* @param X a matrix
*
* @param Y a matrix
*
* @return X <= Y or le(X, Y)
*
*/
public static Matrix le(Matrix X, Matrix Y) {
return ge(Y, X);
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param x a real scalar
*
* @return X > x or gt(X, x)
*
*/
public static Matrix gt(Matrix X, double x) {
Matrix res = X.minus(x);
// res must be a dense matrix.
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] > 0)
resRow[j] = 1;
else
resRow[j] = 0;
}
}
return res;
}
/**
* Do element by element comparisons between x and X and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param x a real scalar
*
* @param X a matrix
*
* @return x > X or gt(x, X)
*/
public static Matrix gt(double x, Matrix X) {
Matrix res = minus(x, X);
// res must be a dense matrix.
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] > 0)
resRow[j] = 1;
else
resRow[j] = 0;
}
}
return res;
}
/**
* Do element by element comparisons between X and Y and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
* X and Y must have the same dimensions.
*
* @param X a matrix
*
* @param Y a matrix
*
* @return X > Y or gt(X, Y)
*
*/
public static Matrix gt(Matrix X, Matrix Y) {
Matrix res = full(X.minus(Y));
int M = res.getRowDimension();
int N = res.getColumnDimension();
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < M; i++) {
resRow = resData[i];
for (int j = 0; j < N; j++) {
if (resRow[j] > 0)
resRow[j] = 1;
else
resRow[j] = 0;
}
}
return res;
}
/**
* Do element by element comparisons between X and x and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param X a matrix
*
* @param x a real scalar
*
* @return X < x or lt(X, x)
*
*/
public static Matrix lt(Matrix X, double x) {
return gt(x, X);
}
/**
* Do element by element comparisons between x and X and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
*
* @param x a real scalar
*
* @param X a matrix
*
* @return x < X or lt(x, X)
*
*/
public static Matrix lt(double x, Matrix X) {
return gt(X, x);
}
/**
* Do element by element comparisons between X and Y and returns
* a matrix of the same size with elements set to 1 where the
* relation is true and elements set to 0 where it is not.
* X and Y must have the same dimensions.
*
* @param X a matrix
*
* @param Y a matrix
*
* @return X < Y or lt(X, Y)
*
*/
public static Matrix lt(Matrix X, Matrix Y) {
return gt(Y, X);
}
/**
* Compute element-wise absolute value of all elements of matrix.
*
* @param A a matrix
*
* @return abs(A)
*/
public static Matrix abs(Matrix A) {
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
Matrix res = A.copy();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
for (int j = 0; j < nCol; j++) {
resRow[j] = Math.abs(resRow[j]);
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
for (int k = 0; k < pr.length; k++) {
pr[k] = Math.abs(pr[k]);
}
}
return res;
}
/**
* Compute element-wise absolute value of all elements of a vector.
*
* @param V a vector
*
* @return abs(V)
*/
public static Vector abs(Vector V) {
Vector res = V.copy();
if (res instanceof DenseVector) {
double[] pr = ((DenseVector) res).getPr();
for (int k = 0; k < res.getDim(); k++) {
pr[k] = Math.abs(pr[k]);
}
} else if (res instanceof SparseVector) {
double[] pr = ((SparseVector) res).getPr();
int nnz = ((SparseVector) res).getNNZ();
for (int k = 0; k < nnz; k++) {
pr[k] = Math.abs(pr[k]);
}
}
return res;
}
/**
* Compute element-wise exponentiation of a vector.
*
* @param A a real matrix
*
* @param p exponent
*
* @return A.^p
*/
public static Matrix pow(Matrix A, double p) {
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
Matrix res = A.copy();
if (res instanceof DenseMatrix) {
double[][] resData = ((DenseMatrix) res).getData();
double[] resRow = null;
for (int i = 0; i < nRow; i++) {
resRow = resData[i];
for (int j = 0; j < nCol; j++) {
resRow[j] = Math.pow(resRow[j], p);
}
}
} else if (res instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) res).getPr();
for (int k = 0; k < pr.length; k++) {
pr[k] = Math.pow(pr[k], p);
}
}
return res;
}
/**
* Compute element-wise exponentiation of a vector.
*
* @param V a real vector
*
* @param p exponent
*
* @return V.^p
*/
public static Vector pow(Vector V, double p) {
Vector res = V.copy();
if (res instanceof DenseVector) {
double[] pr = ((DenseVector) res).getPr();
for (int k = 0; k < res.getDim(); k++) {
pr[k] = Math.pow(pr[k], p);
}
} else if (res instanceof SparseVector) {
double[] pr = ((SparseVector) res).getPr();
int nnz = ((SparseVector) res).getNNZ();
for (int k = 0; k < nnz; k++) {
pr[k] = Math.pow(pr[k], p);
}
}
return res;
}
/**
* Compute the maximal value and the corresponding row
* index of a vector V. The index starts from 0.
*
* @param V a real vector
*
* @return a {@code double} array [M, IX] where M is the
* maximal value, and IX is the corresponding
* index
*/
public static double[] max(Vector V) {
double[] res = new double[2];
double maxVal = Double.NEGATIVE_INFINITY;
double maxIdx = -1;
double v = 0;
if (V instanceof DenseVector) {
double[] pr = ((DenseVector) V).getPr();
for (int k = 0; k < pr.length; k++) {
v = pr[k];
if (maxVal < v) {
maxVal = v;
maxIdx = k;
}
}
} else if (V instanceof SparseVector) {
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int dim = V.getDim();
if (nnz == 0) {
maxVal = 0;
maxIdx = 0;
} else {
int lastIdx = -1;
int currentIdx = 0;
for (int k = 0; k < nnz; k++) {
currentIdx = ir[k];
for (int i = lastIdx + 1; i < currentIdx;) {
if (maxVal < 0) {
maxVal = 0;
maxIdx = i;
}
break;
}
v = pr[k];
if (maxVal < v) {
maxVal = v;
maxIdx = currentIdx;
}
lastIdx = currentIdx;
}
for (int i = lastIdx + 1; i < dim;) {
if (maxVal < 0) {
maxVal = 0;
maxIdx = i;
}
break;
}
}
}
res[0] = maxVal;
res[1] = maxIdx;
return res;
}
/**
* Compute the maximal value for each column and the
* corresponding row indices of a matrix A. The row index
* starts from 0.
*
* @param A a real matrix
*
* @return a {@code Vector} array [M, IX] where M contains
* the maximal values, and IX contains the corresponding
* indices
*/
public static Vector[] max(Matrix A) {
return max(A, 1);
}
/**
* Compute the maximal value for each row or each column
* and the corresponding indices of a matrix A. The row or
* column indices start from 0.
*
* @param A a 2D {@code double} array
*
* @param dim 1: column-wise; 2: row-wise
*
* @return a {@code double[]} array [M, IX] where M contains
* the maximal values, and IX contains the corresponding
* indices
*/
public static double[][] max(double[][] A, int dim) {
double[][] res = new double[2][];
double[][] AData = A;
int M = A.length;
int N = A[0].length;
double maxVal = 0;
int maxIdx = -1;
double v = 0;
double[] maxValues = null;
maxValues = ArrayOperator.allocate1DArray(N, Double.NEGATIVE_INFINITY);
double[] maxIndices = null;
maxIndices = ArrayOperator.allocate1DArray(N, 0);
double[] ARow = null;
if (dim == 1) {
maxValues = ArrayOperator.allocate1DArray(N, Double.NEGATIVE_INFINITY);
maxIndices = ArrayOperator.allocate1DArray(N, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
for (int j = 0; j < N; j++) {
v = ARow[j];
if (maxValues[j] < v) {
maxValues[j] = v;
maxIndices[j] = i;
}
}
}
} else if (dim == 2) {
maxValues = ArrayOperator.allocate1DArray(M, Double.NEGATIVE_INFINITY);
maxIndices = ArrayOperator.allocate1DArray(M, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
maxVal = ARow[0];
maxIdx = 0;
for (int j = 1; j < N; j++) {
v = ARow[j];
if (maxVal < v) {
maxVal = v;
maxIdx = j;
}
}
maxValues[i] = maxVal;
maxIndices[i] = maxIdx;
}
}
res[0] = maxValues;
res[1] = maxIndices;
return res;
}
/**
* Compute the maximal value for each row or each column
* and the corresponding indices of a matrix A. The row or
* column indices start from 0.
*
* @param A a real matrix
*
* @param dim 1: column-wise; 2: row-wise
*
* @return a {@code Vector} array [M, IX] where M contains
* the maximal values, and IX contains the corresponding
* indices
*/
public static Vector[] max(Matrix A, int dim) {
Vector[] res = new DenseVector[2];
int M = A.getRowDimension();
int N = A.getColumnDimension();
double maxVal = 0;
int maxIdx = -1;
double v = 0;
double[] maxValues = null;
maxValues = ArrayOperator.allocate1DArray(N, Double.NEGATIVE_INFINITY);
double[] maxIndices = null;
maxIndices = ArrayOperator.allocate1DArray(N, 0);
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
if (dim == 1) {
maxValues = ArrayOperator.allocate1DArray(N, Double.NEGATIVE_INFINITY);
maxIndices = ArrayOperator.allocate1DArray(N, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
for (int j = 0; j < N; j++) {
v = ARow[j];
if (maxValues[j] < v) {
maxValues[j] = v;
maxIndices[j] = i;
}
}
}
} else if (dim == 2) {
maxValues = ArrayOperator.allocate1DArray(M, Double.NEGATIVE_INFINITY);
maxIndices = ArrayOperator.allocate1DArray(M, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
maxVal = ARow[0];
maxIdx = 0;
for (int j = 1; j < N; j++) {
v = ARow[j];
if (maxVal < v) {
maxVal = v;
maxIdx = j;
}
}
maxValues[i] = maxVal;
maxIndices[i] = maxIdx;
}
}
} else if (A instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) A).getPr();
if (dim == 1) {
maxValues = ArrayOperator.allocate1DArray(N, Double.NEGATIVE_INFINITY);
maxIndices = ArrayOperator.allocate1DArray(N, 0);
int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
for (int j = 0; j < N; j++) {
if (jc[j] == jc[j + 1]) {
maxValues[j] = 0;
maxIndices[j] = 0;
continue;
}
maxVal = Double.NEGATIVE_INFINITY;
maxIdx = -1;
int lastRowIdx = -1;
int currentRowIdx = 0;
for (int k = jc[j]; k < jc[j + 1]; k++) {
currentRowIdx = ir[k];
for (int r = lastRowIdx + 1; r < currentRowIdx;) {
if (maxVal < 0) {
maxVal = 0;
maxIdx = r;
}
break;
}
v = pr[k];
if (maxVal < v) {
maxVal = v;
maxIdx = ir[k];
}
lastRowIdx = currentRowIdx;
}
for (int r = lastRowIdx + 1; r < M;) {
if (maxVal < 0) {
maxVal = 0;
maxIdx = r;
}
break;
}
maxValues[j] = maxVal;
maxIndices[j] = maxIdx;
}
} else if (dim == 2) {
maxValues = ArrayOperator.allocate1DArray(M, Double.NEGATIVE_INFINITY);
maxIndices = ArrayOperator.allocate1DArray(M, 0);
int[] ic = ((SparseMatrix) A).getIc();
int[] jr = ((SparseMatrix) A).getJr();
int[] valCSRIndices = ((SparseMatrix) A).getValCSRIndices();
for (int i = 0; i < M; i++) {
if (jr[i] == jr[i + 1]) {
maxValues[i] = 0;
maxIndices[i] = 0;
continue;
}
maxVal = Double.NEGATIVE_INFINITY;
maxIdx = -1;
int lastColumnIdx = -1;
int currentColumnIdx = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
currentColumnIdx = ic[k];
for (int c = lastColumnIdx + 1; c < currentColumnIdx;) {
if (maxVal < 0) {
maxVal = 0;
maxIdx = c;
}
break;
}
v = pr[valCSRIndices[k]];
if (maxVal < v) {
maxVal = v;
maxIdx = ic[k];
}
lastColumnIdx = currentColumnIdx;
}
for (int c = lastColumnIdx + 1; c < N;) {
if (maxVal < 0) {
maxVal = 0;
maxIdx = c;
}
break;
}
maxValues[i] = maxVal;
maxIndices[i] = maxIdx;
}
}
}
res[0] = new DenseVector(maxValues);
res[1] = new DenseVector(maxIndices);
return res;
}
/**
* Compute the minimal value and the corresponding row
* index of a vector V. The index starts from 0.
*
* @param V a real vector
*
* @return a {@code double} array [M, IX] where M is the
* minimal value, and IX is the corresponding
* index
*/
public static double[] min(Vector V) {
double[] res = new double[2];
double minVal = Double.POSITIVE_INFINITY;
double minIdx = -1;
double v = 0;
if (V instanceof DenseVector) {
double[] pr = ((DenseVector) V).getPr();
for (int k = 0; k < pr.length; k++) {
v = pr[k];
if (minVal > v) {
minVal = v;
minIdx = k;
}
}
} else if (V instanceof SparseVector) {
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
int dim = V.getDim();
if (nnz == 0) {
minVal = 0;
minIdx = 0;
} else {
int lastIdx = -1;
int currentIdx = 0;
for (int k = 0; k < nnz; k++) {
currentIdx = ir[k];
for (int i = lastIdx + 1; i < currentIdx;) {
if (minVal > 0) {
minVal = 0;
minIdx = i;
}
break;
}
v = pr[k];
if (minVal > v) {
minVal = v;
minIdx = currentIdx;
}
lastIdx = currentIdx;
}
for (int i = lastIdx + 1; i < dim;) {
if (minVal > 0) {
minVal = 0;
minIdx = i;
}
break;
}
}
}
res[0] = minVal;
res[1] = minIdx;
return res;
}
/**
* Compute the minimal value for each column and the
* corresponding row indices of a matrix A. The row index
* starts from 0.
*
* @param A a real matrix
*
* @return a {@code Vector} array [M, IX] where M contains
* the minimal values, and IX contains the corresponding
* indices
*/
public static Vector[] min(Matrix A) {
return min(A, 1);
}
/**
* Compute the minimal value for each row or each column
* and the corresponding indices of a matrix A. The row or
* column indices start from 0.
*
* @param A a real matrix
*
* @param dim 1: column-wise; 2: row-wise
*
* @return a {@code Vector} array [M, IX] where M contains
* the minimal values, and IX contains the corresponding
* indices
*/
public static Vector[] min(Matrix A, int dim) {
Vector[] res = new DenseVector[2];
int M = A.getRowDimension();
int N = A.getColumnDimension();
double minVal = 0;
int minIdx = -1;
double v = 0;
double[] minValues = null;
minValues = ArrayOperator.allocate1DArray(N, Double.POSITIVE_INFINITY);
double[] minIndices = null;
minIndices = ArrayOperator.allocate1DArray(N, 0);
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
if (dim == 1) {
minValues = ArrayOperator.allocate1DArray(N, Double.POSITIVE_INFINITY);
minIndices = ArrayOperator.allocate1DArray(N, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
for (int j = 0; j < N; j++) {
v = ARow[j];
if (minValues[j] > v) {
minValues[j] = v;
minIndices[j] = i;
}
}
}
} else if (dim == 2) {
minValues = ArrayOperator.allocate1DArray(M, Double.POSITIVE_INFINITY);
minIndices = ArrayOperator.allocate1DArray(M, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
minVal = ARow[0];
minIdx = 0;
for (int j = 1; j < N; j++) {
v = ARow[j];
if (minVal > v) {
minVal = v;
minIdx = j;
}
}
minValues[i] = minVal;
minIndices[i] = minIdx;
}
}
} else if (A instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) A).getPr();
if (dim == 1) {
minValues = ArrayOperator.allocate1DArray(N, Double.POSITIVE_INFINITY);
minIndices = ArrayOperator.allocate1DArray(N, 0);
int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
for (int j = 0; j < N; j++) {
if (jc[j] == jc[j + 1]) {
minValues[j] = 0;
minIndices[j] = 0;
continue;
}
minVal = Double.POSITIVE_INFINITY;
minIdx = -1;
int lastRowIdx = -1;
int currentRowIdx = 0;
for (int k = jc[j]; k < jc[j + 1]; k++) {
currentRowIdx = ir[k];
for (int r = lastRowIdx + 1; r < currentRowIdx;) {
if (minVal > 0) {
minVal = 0;
minIdx = r;
}
break;
}
v = pr[k];
if (minVal > v) {
minVal = v;
minIdx = ir[k];
}
lastRowIdx = currentRowIdx;
}
for (int r = lastRowIdx + 1; r < M;) {
if (minVal > 0) {
minVal = 0;
minIdx = r;
}
break;
}
minValues[j] = minVal;
minIndices[j] = minIdx;
}
} else if (dim == 2) {
minValues = ArrayOperator.allocate1DArray(M, Double.POSITIVE_INFINITY);
minIndices = ArrayOperator.allocate1DArray(M, 0);
int[] ic = ((SparseMatrix) A).getIc();
int[] jr = ((SparseMatrix) A).getJr();
int[] valCSRIndices = ((SparseMatrix) A).getValCSRIndices();
for (int i = 0; i < M; i++) {
if (jr[i] == jr[i + 1]) {
minValues[i] = 0;
minIndices[i] = 0;
continue;
}
minVal = Double.POSITIVE_INFINITY;
minIdx = -1;
int lastColumnIdx = -1;
int currentColumnIdx = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
currentColumnIdx = ic[k];
for (int c = lastColumnIdx + 1; c < currentColumnIdx;) {
if (minVal > 0) {
minVal = 0;
minIdx = c;
}
break;
}
v = pr[valCSRIndices[k]];
if (minVal > v) {
minVal = v;
minIdx = ic[k];
}
lastColumnIdx = currentColumnIdx;
}
for (int c = lastColumnIdx + 1; c < N;) {
if (minVal > 0) {
minVal = 0;
minIdx = c;
}
break;
}
minValues[i] = minVal;
minIndices[i] = minIdx;
}
}
}
res[0] = new DenseVector(minValues);
res[1] = new DenseVector(minIndices);
return res;
}
/**
* Compute the sum of elements for each row of a real matrix.
*
* @param A a real matrix
*
* @return a dense vector containing the sum of elements of
* each row of A, i.e. sum(A, 1)
*/
public static DenseVector sum(Matrix A) {
return sum(A, 1);
}
/**
* Compute the sum of elements for each row or each column of
* a real matrix.
*
* @param A a real matrix
*
* @param dim 1: column-wise; 2: row-wise
*
* @return a dense vector containing the sum of elements of
* each row or each column of A
*/
public static DenseVector sum(Matrix A, int dim) {
double[] sumValues = null;
int M = A.getRowDimension();
int N = A.getColumnDimension();
double s = 0;
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
double[] ARow = null;
if (dim == 1) {
sumValues = ArrayOperator.allocate1DArray(N, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
for (int j = 0; j < N; j++) {
s = ARow[j];
if (s != 0)
sumValues[j] += s;
}
}
} else if (dim == 2) {
sumValues = ArrayOperator.allocate1DArray(M, 0);
for (int i = 0; i < M; i++) {
ARow = AData[i];
s = 0;
for (int j = 0; j < N; j++) {
s += ARow[j];
}
sumValues[i] = s;
}
}
} else if (A instanceof SparseMatrix) {
double[] pr = ((SparseMatrix) A).getPr();
if (dim == 1) {
sumValues = ArrayOperator.allocate1DArray(N, 0);
// int[] ir = ((SparseMatrix) A).getIr();
int[] jc = ((SparseMatrix) A).getJc();
for (int j = 0; j < N; j++) {
if (jc[j] == jc[j + 1]) {
sumValues[j] = 0;
continue;
}
s = 0;
for (int k = jc[j]; k < jc[j + 1]; k++) {
s += pr[k];
}
sumValues[j] = s;
}
} else if (dim == 2) {
sumValues = ArrayOperator.allocate1DArray(M, 0);
// int[] ic = ((SparseMatrix) A).getIc();
int[] jr = ((SparseMatrix) A).getJr();
int[] valCSRIndices = ((SparseMatrix) A).getValCSRIndices();
for (int i = 0; i < M; i++) {
if (jr[i] == jr[i + 1]) {
sumValues[i] = 0;
continue;
}
s = 0;
for (int k = jr[i]; k < jr[i + 1]; k++) {
s += pr[valCSRIndices[k]];
}
sumValues[i] = s;
}
}
}
return new DenseVector(sumValues);
}
/**
* Compute the sum of a 1D {@code double} array.
*
* @param V a 1D {@code double} array
*
* @return sum(V)
*//*
public static double sum(double[] V) {
double res = 0;
for (int i = 0; i < V.length; i++)
res += V[i];
return res;
}*/
/**
* Compute the sum of all elements of a vector.
*
* @param V a real dense or sparse vector
*
* @return sum(V)
*/
public static double sum(Vector V) {
double res = 0;
if (V instanceof DenseVector) {
double[] pr = ((DenseVector) V).getPr();
for (int k = 0; k < pr.length; k++) {
res += pr[k];
}
} else if (V instanceof SparseVector) {
double[] pr = ((SparseVector) V).getPr();
int nnz = ((SparseVector) V).getNNZ();
for (int k = 0; k < nnz; k++) {
res += pr[k];
}
}
return res;
}
/**
* Compute the Euclidean norm of a matrix.
*
* @param A a real matrix
*
* @return ||A||_2
*/
public static double norm(Matrix A) {
return norm(A, 2);
}
/**
* Compute the induced vector norm of a matrix (row or column
* matrices are allowed).
*
* @param A a matrix or a vector
*
* @param type 1, 2, or, inf for a matrix or a positive real
* number for a row or column matrix
*
* @return ||A||_{type}
*
*/
public static double norm(Matrix A, double type) {
double res = 0;
int nRow = A.getRowDimension();
int nCol = A.getColumnDimension();
if (nRow == 1) {
if (Double.isInfinite(type)) {
return max(max(abs(A), 2)[0])[0];
} else if (type > 0) {
return Math.pow(sum(sum(pow(abs(A), type))), 1.0 / type);
} else {
System.err.printf("Error norm type: %f\n", type);
System.exit(1);
}
}
if (nCol == 1) {
if (Double.isInfinite(type)) {
return max(max(abs(A), 1)[0])[0];
} else if (type > 0) {
return Math.pow(sum(sum(pow(abs(A), type))), 1.0 / type);
} else {
System.err.printf("Error norm type: %f\n", type);
System.exit(1);
}
}
if (type == 2) {
double eigenvalue = EigenValueDecomposition.computeEigenvalues(A.transpose().mtimes(A))[0];
res = eigenvalue <= 0 ? 0 : Math.sqrt(eigenvalue);
// res = SingularValueDecomposition.computeSingularValues(A)[0];
} else if (Double.isInfinite(type)) {
res = max(sum(abs(A), 2))[0];
} else if (type == 1) {
res = max(sum(abs(A), 1))[0];
} else {
System.err.printf("Sorry, %f-norm of a matrix is not supported currently.\n", type);
}
return res;
}
/**
* Compute the norm of a matrix (row or column matrices
* are allowed).
*
* @param A a matrix
*
* @param type 1, 2
*
* @return ||A||_{type}
*
*/
public static double norm(Matrix A, int type) {
return norm(A, (double)type);
}
/**
* Calculate the Frobenius norm of a matrix A.
*
* @param A a matrix
*
* @param type can only be "fro"
*
* @return ||A||_F
*
*/
public static double norm(Matrix A, String type) {
double res = 0;
if (type.compareToIgnoreCase("fro") == 0) {
res = Math.sqrt(innerProduct(A, A));
} else if (type.equals("inf")){
res = norm(A, inf);
} else if (type.equals("nuclear")) {
res = ArrayOperator.sum(SingularValueDecomposition.computeSingularValues(A));
} else {
System.err.println(String.format("Norm %s unimplemented!\n" , type));
}
return res;
}
/**
* Compute the norm of a vector.
*
* @param V a real vector
*
* @param p a positive {@code double} value
*
* @return ||V||_p
*/
public static double norm(Vector V, double p) {
if (p == 1) {
return sum(abs(V));
} else if (p == 2) {
return Math.sqrt(innerProduct(V, V));
} else if (Double.isInfinite(p)) {
return max(abs(V))[0];
} else if (p > 0) {
return Math.pow(sum(pow(abs(V), p)), 1.0 / p);
} else {
System.err.println("Wrong argument for p");
System.exit(1);
}
return -1;
}
/**
* Compute the Euclidean norm of a vector.
*
* @param V a real vector
*
* @return ||V||_2
*/
public static double norm(Vector V) {
return norm(V, 2);
}
/**
* Compute the norm of a vector.
*
* @param V a 1D {@code double} array
*
* @param p a positive {@code double} value
*
* @return ||V||_p
*/
public static double norm(double[] V, double p) {
return norm(new DenseVector(V), p);
}
/**
* Compute the Euclidean norm of a vector.
*
* @param V a 1D {@code double} array
*
* @return ||V||_2
*/
public static double norm(double[] V) {
return norm(V, 2);
}
/**
* Compute the inner product of two vectors, i.e. res = <V1, V2>.
*
* @param V1 the first vector
*
* @param V2 the second vector
*
* @return <V1, V2>
*/
public static double innerProduct(Vector V1, Vector V2) {
if (V1.getDim() != V2.getDim()) {
System.err.println("Dimension doesn't match.");
System.exit(1);
}
double res = 0;
double v = 0;
if (V1 instanceof DenseVector) {
double[] pr1 = ((DenseVector) V1).getPr();
if (V2 instanceof DenseVector) {
if (V1 == V2) {
for (int k = 0; k < pr1.length; k++) {
v = pr1[k];
res += v * v;
}
return res;
}
double[] pr2 = ((DenseVector) V2).getPr();
for (int k = 0; k < pr1.length; k++) {
res += pr1[k] * pr2[k];
}
} else if (V2 instanceof SparseVector) {
int[] ir = ((SparseVector) V2).getIr();
double[] pr2 = ((SparseVector) V2).getPr();
int nnz = ((SparseVector) V2).getNNZ();
for (int k = 0; k < nnz; k++) {
res += pr1[ir[k]] * pr2[k];
}
}
} else if (V1 instanceof SparseVector) {
if (V2 instanceof DenseVector) {
return innerProduct(V2, V1);
} else if (V2 instanceof SparseVector) {
int[] ir1 = ((SparseVector) V1).getIr();
double[] pr1 = ((SparseVector) V1).getPr();
int nnz1 = ((SparseVector) V1).getNNZ();
if (V1 == V2) {
for (int k = 0; k < nnz1; k++) {
v = pr1[k];
res += v * v;
}
return res;
} else {
int[] ir2 = ((SparseVector) V2).getIr();
double[] pr2 = ((SparseVector) V2).getPr();
int nnz2 = ((SparseVector) V2).getNNZ();
int k1 = 0;
int k2 = 0;
int i1 = 0;
int i2 = 0;
while (true) {
if (k1 >= nnz1 || k2 >= nnz2) {
break;
}
i1 = ir1[k1];
i2 = ir2[k2];
if (i1 < i2) {
k1++;
} else if (i1 > i2) {
k2++;
} else {
res += pr1[k1] * pr2[k2];
k1++;
k2++;
}
}
}
}
}
return res;
}
/**
* Compute the inner product of two matrices, i.e. res = <A, B>.
*
* @param A a matrix
*
* @param B a matrix
*
* @return <A, B>
*
*/
public static double innerProduct(Matrix A, Matrix B) {
double s = 0;
int M = A.getRowDimension();
int N = A.getColumnDimension();
if (A instanceof DenseMatrix) {
double[][] AData = ((DenseMatrix) A).getData();
if (B instanceof DenseMatrix) {
double[][] BData = ((DenseMatrix) B).getData();
double[] BRow = null;
double[] ARow = null;
// double[] resRow = null;
for (int i = 0; i < M; i++) {
ARow = AData[i];
BRow = BData[i];
for (int j = 0; j < N; j++) {
s += ARow[j] * BRow[j];
}
}
} else if (B instanceof SparseMatrix) {
int[] ir = null;
int[] jc = null;
double[] pr = null;
ir = ((SparseMatrix) B).getIr();
jc = ((SparseMatrix) B).getJc();
pr = ((SparseMatrix) B).getPr();
int r = -1;
for (int j = 0; j < B.getColumnDimension(); j++) {
for (int k = jc[j]; k < jc[j + 1]; k++) {
r = ir[k];
// A[r][j] = pr[k]
s += AData[r][j] * pr[k];
}
}
}
} else if (A instanceof SparseMatrix) {
if (B instanceof DenseMatrix) {
return innerProduct(B, A);
} else if (B instanceof SparseMatrix) {
int[] ir1 = null;
int[] jc1 = null;
double[] pr1 = null;
ir1 = ((SparseMatrix) A).getIr();
jc1 = ((SparseMatrix) A).getJc();
pr1 = ((SparseMatrix) A).getPr();
int[] ir2 = null;
int[] jc2 = null;
double[] pr2 = null;
ir2 = ((SparseMatrix) B).getIr();
jc2 = ((SparseMatrix) B).getJc();
pr2 = ((SparseMatrix) B).getPr();
int k1 = 0;
int k2 = 0;
int r1 = -1;
int r2 = -1;
// int i = -1;
double v = 0;
for (int j = 0; j < N; j++) {
k1 = jc1[j];
k2 = jc2[j];
// If the j-th column of A or this is empty, we don't need to compute.
if (k1 == jc1[j + 1] || k2 == jc2[j + 1])
continue;
while (k1 < jc1[j + 1] && k2 < jc2[j + 1]) {
r1 = ir1[k1];
r2 = ir2[k2];
if (r1 < r2) {
k1++;
} else if (r1 == r2) {
// i = r1;
v = pr1[k1] * pr2[k2];
k1++;
k2++;
if (v != 0) {
s += v;
}
} else {
k2++;
}
}
}
}
}
return s;
}
/**
* Calculate the sigmoid of a matrix A by rows. Specifically, supposing
* that the input activation matrix is [a11, a12; a21, a22], the output
* value is
* <p>
* [exp(a11) / exp(a11) + exp(a12), exp(a12) / exp(a11) + exp(a12);
* </br>
* exp(a21) / exp(a21) + exp(a22), exp(a22) / exp(a21) + exp(a22)].
*
* @param A a real matrix
*
* @return sigmoid(A)
*/
public static Matrix sigmoid(Matrix A) {
double[][] data = full(A.copy()).getData();
int M = A.getRowDimension();
int N = A.getColumnDimension();
double[] row_i = null;
double old = 0;
double current = 0;
double max = 0;
double sum = 0;
double v = 0;
for (int i = 0; i < M; i++) {
row_i = data[i];
old = row_i[0];
current = 0;
max = old;
for (int j = 1; j < N; j++) {
current = row_i[j];
if (max < current)
max = current;
old = current;
}
sum = 0;
for (int j = 0; j < N; j++) {
v = Math.exp(row_i[j] - max);
sum += v;
row_i[j] = v;
}
for (int j = 0; j < N; j++) {
row_i[j] /= sum;
}
}
return new DenseMatrix(data);
}
/**
* Compute the maximum of elements in a 1D {@code double} array.
*
* @param V a 1D {@code double} array
*
* @return max(V)
*/
public static double max(double[] V) {
// double max = max(V, 0, V.length - 1);
double max = V[0];
for (int i = 1; i < V.length; i++) {
if (max < V[i])
max = V[i];
}
return max;
}
/**
* Compute the maximum of elements in an interval
* of a 1D {@code double} array.
*
* @param V a 1D {@code double} array
*
* @param start start index (inclusive)
*
* @param end end index (inclusive)
*
* @return max(V(start:end))
*/
public static double max(double[] V, int start, int end) {
if (start == end)
return V[start];
if (start == end - 1)
return Math.max(V[start], V[end]);
int middle = (start + end) / 2;
double leftMax = max(V, start, middle);
double rightMax = max(V, middle + 1, end);
return Math.max(leftMax, rightMax);
}
/**
* Calculate the left division of the form A \ B. A \ B is the
* matrix division of A into B, which is roughly the same as
* INV(A)*B , except it is computed in a different way. For
* implementation, we actually solve the system of linear
* equations A * X = B.
*
* @param A divisor
*
* @param B dividend
*
* @return A \ B
*
*/
public static Matrix mldivide(Matrix A, Matrix B) {
return new QRDecomposition(A).solve(B);
}
/**
* Calculate the right division of B into A, i.e., A / B. For
* implementation, we actually solve the system of linear
* equations X * B = A.
* <p>
* Note: X = A / B <=> X * B = A <=> B' * X' = A' <=> X' = B' \ A'
* <=> X = (B' \ A')'
* </p>
*
* @param A dividend
*
* @param B divisor
*
* @return A / B
*
*/
public static Matrix mrdivide(Matrix A, Matrix B) {
return mldivide(B.transpose(), A.transpose()).transpose();
}
/**
* Compute the rank of a matrix. The rank function provides
* an estimate of the number of linearly independent rows or
* columns of a matrix.
*
* @param A a matrix
*
* @return rank of the given matrix
*/
public static int rank(Matrix A) {
return SingularValueDecomposition.rank(A);
}
/**
* Construct an m-by-n dense identity matrix.
*
* @param m number of rows
*
* @param n number of columns
*
* @return an m-by-n dense identity matrix
*
*/
public static DenseMatrix eye(int m, int n) {
double[][] res = ArrayOperator.allocate2DArray(m, n, 0);
int len = m >= n ? n : m;
for (int i = 0; i < len; i++) {
res[i][i] = 1;
}
return new DenseMatrix(res);
}
/**
* Construct an n-by-n dense identity matrix.
*
* @param n number of rows and columns
*
* @return an n-by-n dense identity matrix
*
*/
public static DenseMatrix eye(int n) {
return eye(n, n);
}
/**
* Generate a dense identity matrix with its size
* specified by a two dimensional integer array.
*
* @param size a two dimensional integer array
*
* @return a dense identity matrix with its shape specified by size
*
*/
public static Matrix eye(int... size) {
if (size.length != 2) {
System.err.println("Input size vector should have two elements!");
}
return eye(size[0], size[1]);
}
/**
* Construct an m-by-n sparse identity matrix.
*
* @param m number of rows
*
* @param n number of columns
*
* @return an m-by-n sparse identity matrix
*
*/
public static SparseMatrix speye(int m, int n) {
TreeMap<Pair<Integer, Integer>, Double> map = new TreeMap<Pair<Integer, Integer>, Double>();
int len = m >= n ? n : m;
for (int i = 0; i < len; i++) {
map.put(Pair.of(i, i), 1.0);
}
return SparseMatrix.createSparseMatrix(map, m, n);
}
/**
* Construct an n-by-n sparse identity matrix.
*
* @param n number of rows and columns
*
* @return an n-by-n sparse identity matrix
*
*/
public static SparseMatrix speye(int n) {
/*TreeMap<Pair<Integer, Integer>, Double> map = new TreeMap<Pair<Integer, Integer>, Double>();
for (int i = 0; i < n; i++) {
map.put(Pair.of(i, i), 1.0);
}
return SparseMatrix.createSparseMatrix(map, n, n);*/
return speye(n, n);
}
/**
* Generate a sparse identity matrix with its size
* specified by a two dimensional integer array.
*
* @param size a two dimensional integer array
*
* @return a sparse identity matrix with its shape specified by size
*
*/
public static Matrix speye(int... size) {
if (size.length != 2) {
System.err.println("Input size vector should have two elements!");
}
return speye(size[0], size[1]);
}
/**
* Create a m-by-n Hilbert matrix. The elements of the
* Hilbert matrices are H(i, j) = 1 / (i + j + 1), where
* i and j are zero based indices.
*
* @param m number of rows
*
* @param n number of columns
*
* @return a m-by-n Hilbert matrix
*/
public static Matrix hilb(int m, int n) {
DenseMatrix A = new DenseMatrix(m, n);
double[][] data = A.getData();
double[] A_i = null;
for (int i = 0; i < m; i++) {
A_i = data[i];
for (int j = 0; j < n; j++) {
A_i[j] = 1.0 / (i + j + 1);
}
}
return A;
}
public static Vector times(Vector V1, Vector V2) {
return V1.times(V2);
}
public static Vector plus(Vector V1, Vector V2) {
return V1.plus(V2);
}
public static Vector minus(Vector V1, Vector V2) {
return V1.minus(V2);
}
public static Matrix times(Matrix X, Matrix Y) {
int nX = X.getColumnDimension();
int dX = X.getRowDimension();
int nY = Y.getColumnDimension();
int dY = Y.getRowDimension();
if (dX == 1 && nX == 1) {
return times(X.getEntry(0, 0), Y);
} else if (dY == 1 && nY == 1) {
return times(X, Y.getEntry(0, 0));
}
if (nX != nY || dX != dY) {
System.err.println("The operands for Hadmada product should be of same shapes!");
}
return X.times(Y);
}
public static Matrix times(Matrix A, double v) {
return A.times(v);
}
public static Matrix times(double v, Matrix A) {
return A.times(v);
}
public static Matrix mtimes(Matrix M1, Matrix M2) {
return M1.mtimes(M2);
}
public static Matrix plus(Matrix M1, Matrix M2) {
return M1.plus(M2);
}
public static Matrix plus(Matrix A, double v) {
return A.plus(v);
}
public static Matrix plus(double v, Matrix A) {
return A.plus(v);
}
public static Matrix minus(Matrix M1, Matrix M2) {
return M1.minus(M2);
}
public static Matrix minus(Matrix A, double v) {
return A.minus(v);
}
/**
* res = v - A.
* @param v
* @param A
* @return v - A
*/
public static Matrix minus(double v, Matrix A) {
return uminus(A).plus(v);
}
/**
* Set all the elements of X to be those of Y, this is particularly
* useful when we want to change elements of the object referred by X
* rather than the reference X itself.
*
* @param X a matrix to be set
*
* @param Y a matrix to set X
*
*/
public static void setMatrix(Matrix X, Matrix Y) {
assign(X, Y);
}
/**
* Unary minus.
*
* @param A a matrix
*
* @return -A
*/
public static Matrix uminus(Matrix A) {
if (A == null) {
return null;
} else {
return A.times(-1);
}
}
public static DenseVector full(Vector V) {
if (V instanceof SparseVector) {
int dim = V.getDim();
int[] ir = ((SparseVector) V).getIr();
double[] pr = ((SparseVector) V).getPr();
double[] values = ArrayOperator.allocateVector(dim, 0);
for (int k = 0; k < ((SparseVector) V).getNNZ(); k++) {
values[ir[k]] = pr[k];
}
return new DenseVector(values);
} else {
return (DenseVector) V;
}
}
public static SparseVector sparse(Vector V) {
if (V instanceof DenseVector) {
double[] values = ((DenseVector) V).getPr();
TreeMap<Integer, Double> map = new TreeMap<Integer, Double>();
for (int k = 0; k < values.length; k++) {
if (values[k] != 0) {
map.put(k, values[k]);
}
}
int nnz = map.size();
int[] ir = new int[nnz];
double[] pr = new double[nnz];
int dim = values.length;
int ind = 0;
for (Entry<Integer, Double> entry : map.entrySet()) {
ir[ind] = entry.getKey();
pr[ind] = entry.getValue();
ind++;
}
return new SparseVector(ir, pr, nnz, dim);
} else {
return (SparseVector) V;
}
}
public static DenseMatrix full(Matrix S) {
if (S instanceof SparseMatrix) {
int M = S.getRowDimension();
int N = S.getColumnDimension();
double[][] data = new double[M][];
int[] ic = ((SparseMatrix) S).getIc();
int[] jr = ((SparseMatrix) S).getJr();
int[] valCSRIndices = ((SparseMatrix) S).getValCSRIndices();
double[] pr = ((SparseMatrix) S).getPr();
for (int i = 0; i < M; i++) {
double[] rowData = ArrayOperator.allocateVector(N, 0);
for (int k = jr[i]; k < jr[i + 1]; k++) {
rowData[ic[k]] = pr[valCSRIndices[k]];
}
data[i] = rowData;
}
return new DenseMatrix(data);
} else {
return (DenseMatrix) S;
}
}
public static SparseMatrix sparse(Matrix A) {
if (A instanceof DenseMatrix) {
int rIdx = 0;
int cIdx = 0;
int nzmax = 0;
double value = 0;
TreeMap<Pair<Integer, Integer>, Double> map = new TreeMap<Pair<Integer, Integer>, Double>();
int numRows = A.getRowDimension();
int numColumns = A.getColumnDimension();
double[][] data = ((DenseMatrix) A).getData();
for (int j = 0; j < numColumns; j++) {
cIdx = j;
for (int i = 0; i < numRows; i++) {
rIdx = i;
value = data[i][j];
if (value != 0) {
map.put(Pair.of(cIdx, rIdx), value);
nzmax++;
}
}
}
int[] ir = new int[nzmax];
int[] jc = new int[numColumns + 1];
double[] pr = new double[nzmax];
int k = 0;
jc[0] = 0;
int currentColumn = 0;
for (Entry<Pair<Integer, Integer>, Double> entry : map.entrySet()) {
rIdx = entry.getKey().second;
cIdx = entry.getKey().first;
pr[k] = entry.getValue();
ir[k] = rIdx;
if (currentColumn < cIdx) {
jc[currentColumn + 1] = k;
currentColumn++;
}
k++;
}
while (currentColumn < numColumns) {
jc[currentColumn + 1] = k;
currentColumn++;
}
jc[numColumns] = k;
return SparseMatrix.createSparseMatrixByCSCArrays(ir, jc, pr, numRows, numColumns, nzmax);
} else {
return (SparseMatrix) A;
}
}
public static Vector[] sparseMatrix2SparseRowVectors(Matrix S) {
if (!(S instanceof SparseMatrix)) {
System.err.println("SparseMatrix input is expected.");
System.exit(1);
}
int M = S.getRowDimension();
int N = S.getColumnDimension();
Vector[] Vs = new Vector[M];
int[] ic = ((SparseMatrix) S).getIc();
int[] jr = ((SparseMatrix) S).getJr();
double[] pr = ((SparseMatrix) S).getPr();
int[] valCSRIndices = ((SparseMatrix) S).getValCSRIndices();
int[] indices = null;
double[] values = null;
int nnz = 0;
int dim = N;
for (int r = 0; r < M; r++) {
nnz = jr[r + 1] - jr[r];
indices = new int[nnz];
values = new double[nnz];
int idx = 0;
for (int k = jr[r]; k < jr[r + 1]; k++) {
indices[idx] = ic[k];
values[idx] = pr[valCSRIndices[k]];
idx++;
}
Vs[r] = new SparseVector(indices, values, nnz, dim);
}
return Vs;
}
public static Vector[] sparseMatrix2SparseColumnVectors(Matrix S) {
if (!(S instanceof SparseMatrix)) {
System.err.println("SparseMatrix input is expected.");
System.exit(1);
}
int M = S.getRowDimension();
int N = S.getColumnDimension();
Vector[] Vs = new Vector[N];
int[] ir = ((SparseMatrix) S).getIr();
int[] jc = ((SparseMatrix) S).getJc();
double[] pr = ((SparseMatrix) S).getPr();
int[] indices = null;
double[] values = null;
int nnz = 0;
int dim = M;
for (int c = 0; c < N; c++) {
nnz = jc[c + 1] - jc[c];
indices = new int[nnz];
values = new double[nnz];
int idx = 0;
for (int k = jc[c]; k < jc[c + 1]; k++) {
indices[idx] = ir[k];
values[idx] = pr[k];
idx++;
}
Vs[c] = new SparseVector(indices, values, nnz, dim);
}
return Vs;
}
public static Matrix sparseRowVectors2SparseMatrix(Vector[] Vs) {
// return sparseColumnVectors2SparseMatrix(Vs).transpose();
int nnz = 0;
int numColumns = Vs.length > 0 ? Vs[0].getDim() : 0;
for (int i = 0; i < Vs.length; i++) {
if (!(Vs[i] instanceof SparseVector)) {
fprintf("Vs[%d] should be a sparse vector.%n", i);
System.exit(1);
}
nnz += ((SparseVector) Vs[i]).getNNZ();
if (numColumns != Vs[i].getDim()) {
fprintf("Vs[%d]'s dimension doesn't match.%n", i);
System.exit(1);
}
}
int numRows = Vs.length;
int nzmax = nnz;
int[] ic = new int[nzmax];
int[] jr = new int[numRows + 1];
double[] pr = new double[nzmax];
int rIdx = -1;
int cIdx = -1;
int cnt = 0;
jr[0] = 0;
int currentRow = 0;
for (int i = 0; i < numRows; i++) {
int[] indices = ((SparseVector) Vs[i]).getIr();
double[] values = ((SparseVector) Vs[i]).getPr();
nnz = ((SparseVector) Vs[i]).getNNZ();
for (int k = 0; k < nnz; k++) {
cIdx = indices[k];
rIdx = i;
pr[cnt] = values[k];
ic[cnt] = cIdx;
while (currentRow < rIdx) {
jr[currentRow + 1] = cnt;
currentRow++;
}
cnt++;
}
}
while (currentRow < numRows) {
jr[currentRow + 1] = cnt;
currentRow++;
}
// jr[numColumns] = k;
return SparseMatrix.createSparseMatrixByCSRArrays(ic, jr, pr, numRows, numColumns, nzmax);
}
public static Matrix sparseColumnVectors2SparseMatrix(Vector[] Vs) {
int nnz = 0;
int numRows = Vs.length > 0 ? Vs[0].getDim() : 0;
for (int i = 0; i < Vs.length; i++) {
if (!(Vs[i] instanceof SparseVector)) {
fprintf("Vs[%d] should be a sparse vector.%n", i);
System.exit(1);
}
nnz += ((SparseVector) Vs[i]).getNNZ();
if (numRows != Vs[i].getDim()) {
fprintf("Vs[%d]'s dimension doesn't match.%n", i);
System.exit(1);
}
}
int numColumns = Vs.length;
int nzmax = nnz;
int[] ir = new int[nzmax];
int[] jc = new int[numColumns + 1];
double[] pr = new double[nzmax];
int rIdx = -1;
int cIdx = -1;
int k = 0;
jc[0] = 0;
int currentColumn = 0;
for (int c = 0; c < numColumns; c++) {
int[] indices = ((SparseVector) Vs[c]).getIr();
double[] values = ((SparseVector) Vs[c]).getPr();
nnz = ((SparseVector) Vs[c]).getNNZ();
for (int r = 0; r < nnz; r++) {
rIdx = indices[r];
cIdx = c;
pr[k] = values[r];
ir[k] = rIdx;
while (currentColumn < cIdx) {
jc[currentColumn + 1] = k;
currentColumn++;
}
k++;
}
}
while (currentColumn < numColumns) {
jc[currentColumn + 1] = k;
currentColumn++;
}
// jc[numColumns] = k;
return SparseMatrix.createSparseMatrixByCSCArrays(ir, jc, pr, numRows, numColumns, nzmax);
}
}
