org.ejml.sparse.csc.CommonOps_FSCC Maven / Gradle / Ivy
Show all versions of ejml-fsparse Show documentation
/*
* Copyright (c) 2022, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.ejml.sparse.csc;
import javax.annotation.Generated;
import org.ejml.MatrixDimensionException;
import org.ejml.data.FGrowArray;
import org.ejml.data.FMatrixRMaj;
import org.ejml.data.FMatrixSparseCSC;
import org.ejml.data.IGrowArray;
import org.ejml.dense.row.CommonOps_FDRM;
import org.ejml.interfaces.decomposition.LUSparseDecomposition_F32;
import org.ejml.interfaces.linsol.LinearSolverSparse;
import org.ejml.masks.Mask;
import org.ejml.ops.FOperatorBinary;
import org.ejml.ops.FOperatorBinaryIdx;
import org.ejml.ops.FOperatorUnary;
import org.ejml.ops.IPredicateBinary;
import org.ejml.sparse.FillReducing;
import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC;
import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC;
import org.ejml.sparse.csc.misc.ImplCommonOps_FSCC;
import org.ejml.sparse.csc.mult.ImplMultiplication_FSCC;
import org.jetbrains.annotations.Nullable;
import java.util.Arrays;
import static org.ejml.UtilEjml.*;
/**
* Most common operations on {@link FMatrixSparseCSC}. For example, addition, matrix multiplication, ... etc
*
* @author Peter Abeles
*/
@Generated("org.ejml.sparse.csc.CommonOps_DSCC")
public class CommonOps_FSCC {
private CommonOps_FSCC(){}
/**
* Checks to see if row indicies are sorted into ascending order. O(N)
*
* @return true if sorted and false if not
*/
public static boolean checkIndicesSorted( FMatrixSparseCSC A ) {
for (int j = 0; j < A.numCols; j++) {
int idx0 = A.col_idx[j];
int idx1 = A.col_idx[j + 1];
if (idx0 != idx1 && A.nz_rows[idx0] >= A.numRows)
return false;
for (int i = idx0 + 1; i < idx1; i++) {
int row = A.nz_rows[i];
if (A.nz_rows[i - 1] >= row)
return false;
if (row >= A.numRows)
return false;
}
}
return true;
}
public static boolean checkStructure( FMatrixSparseCSC A ) {
if (A.col_idx.length < A.numCols + 1)
return false;
if (A.col_idx[A.numCols] != A.nz_length)
return false;
if (A.nz_rows.length < A.nz_length)
return false;
if (A.nz_values.length < A.nz_length)
return false;
if (A.col_idx[0] != 0)
return false;
for (int i = 0; i < A.numCols; i++) {
if (A.col_idx[i] > A.col_idx[i + 1]) {
return false;
}
if (A.col_idx[i + 1] - A.col_idx[i] > A.numRows)
return false;
}
if (!checkSortedFlag(A))
return false;
if (checkDuplicateElements(A))
return false;
return true;
}
public static boolean checkSortedFlag( FMatrixSparseCSC A ) {
if (A.indicesSorted)
return checkIndicesSorted(A);
return true;
}
/**
* Checks for duplicate elements. A is sorted
*
* @param A Matrix to be tested.
* @return true if duplicates or false if false duplicates
*/
public static boolean checkDuplicateElements( FMatrixSparseCSC A ) {
A = A.copy(); // create a copy so that it doesn't modify A
A.sortIndices(null);
return !checkSortedFlag(A);
}
/**
* Perform matrix transpose
*
* @param A Input matrix. Not modified
* @param A_t Storage for transpose of 'a'. Must be correct shape. data length might be adjusted.
* @param gw (Optional) Storage for internal workspace. Can be null.
* @return The transposed matrix
*/
public static FMatrixSparseCSC transpose( FMatrixSparseCSC A, @Nullable FMatrixSparseCSC A_t, @Nullable IGrowArray gw ) {
A_t = reshapeOrDeclare(A_t, A.numCols, A.numRows, A.nz_length);
ImplCommonOps_FSCC.transpose(A, A_t, gw);
return A_t;
}
public static FMatrixSparseCSC mult( FMatrixSparseCSC A, FMatrixSparseCSC B,
@Nullable FMatrixSparseCSC outputC ) {
return mult(A, B, outputC, null, null);
}
/**
* Performs matrix multiplication. C = A*B
*
* @param A (Input) Matrix. Not modified.
* @param B (Input) Matrix. Not modified.
* @param outputC (Output) Storage for results. Data length is increased if insufficient.
* @param gw (Optional) Storage for internal workspace. Can be null.
* @param gx (Optional) Storage for internal workspace. Can be null.
*/
public static FMatrixSparseCSC mult( FMatrixSparseCSC A, FMatrixSparseCSC B,
@Nullable FMatrixSparseCSC outputC,
@Nullable IGrowArray gw, @Nullable FGrowArray gx ) {
if (A.numCols != B.numRows)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
outputC = reshapeOrDeclare(outputC, A, A.numRows, B.numCols);
ImplMultiplication_FSCC.mult(A, B, outputC, gw, gx);
return outputC;
}
/**
* Performs matrix multiplication. C = A*B
*
* @param A Matrix
* @param B Dense Matrix
* @param outputC Dense Matrix
*/
public static FMatrixRMaj mult( FMatrixSparseCSC A, FMatrixRMaj B, @Nullable FMatrixRMaj outputC ) {
if (A.numCols != B.numRows)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
outputC = reshapeOrDeclare(outputC, A.numRows, B.numCols);
ImplMultiplication_FSCC.mult(A, B, outputC);
return outputC;
}
/**
* C = C + A*B
*/
public static void multAdd( FMatrixSparseCSC A, FMatrixRMaj B, FMatrixRMaj outputC ) {
if (A.numCols != B.numRows)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
if (A.numRows != outputC.numRows || B.numCols != outputC.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B, outputC));
ImplMultiplication_FSCC.multAdd(A, B, outputC);
}
/**
* Performs matrix multiplication. C = AT*B
*
* @param A Matrix
* @param B Dense Matrix
* @param outputC Dense Matrix
*/
public static FMatrixRMaj multTransA( FMatrixSparseCSC A, FMatrixRMaj B, @Nullable FMatrixRMaj outputC,
@Nullable FGrowArray work ) {
if (A.numRows != B.numRows)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
if (work == null)
work = new FGrowArray();
outputC = reshapeOrDeclare(outputC, A.numCols, B.numCols);
ImplMultiplication_FSCC.multTransA(A, B, outputC, work);
return outputC;
}
/**
* C = C + AT*B
*/
public static void multAddTransA( FMatrixSparseCSC A, FMatrixRMaj B, FMatrixRMaj outputC,
@Nullable FGrowArray work ) {
if (A.numRows != B.numRows)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
if (A.numCols != outputC.numRows || B.numCols != outputC.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B, outputC));
if (work == null)
work = new FGrowArray();
ImplMultiplication_FSCC.multAddTransA(A, B, outputC, work);
}
/**
* Performs matrix multiplication. C = A*BT
*
* @param A Matrix
* @param B Dense Matrix
* @param outputC Dense Matrix
*/
public static FMatrixRMaj multTransB( FMatrixSparseCSC A, FMatrixRMaj B, @Nullable FMatrixRMaj outputC,
@Nullable FGrowArray work ) {
if (A.numCols != B.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
outputC = reshapeOrDeclare(outputC, A.numRows, B.numRows);
if (work == null)
work = new FGrowArray();
ImplMultiplication_FSCC.multTransB(A, B, outputC, work);
return outputC;
}
/**
* C = C + A*BT
*/
public static void multAddTransB( FMatrixSparseCSC A, FMatrixRMaj B, FMatrixRMaj outputC,
@Nullable FGrowArray work ) {
if (A.numCols != B.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
if (A.numRows != outputC.numRows || B.numRows != outputC.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B, outputC));
if (work == null)
work = new FGrowArray();
ImplMultiplication_FSCC.multAddTransB(A, B, outputC, work);
}
/**
* Performs matrix multiplication. C = AT*BT
*
* @param A Matrix
* @param B Dense Matrix
* @param outputC Dense Matrix
*/
public static FMatrixRMaj multTransAB( FMatrixSparseCSC A, FMatrixRMaj B, FMatrixRMaj outputC ) {
if (A.numRows != B.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
outputC = reshapeOrDeclare(outputC, A.numCols, B.numRows);
ImplMultiplication_FSCC.multTransAB(A, B, outputC);
return outputC;
}
/**
* C = C + AT*BT
*/
public static void multAddTransAB( FMatrixSparseCSC A, FMatrixRMaj B, FMatrixRMaj outputC ) {
if (A.numRows != B.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
if (A.numCols != outputC.numRows || B.numRows != outputC.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B, outputC));
ImplMultiplication_FSCC.multAddTransAB(A, B, outputC);
}
/**
* Given a symmetric matrix, which is represented by a lower triangular matrix, convert it back into
* a full symmetric matrix
*
* @param A (Input) Lower triangular matrix
* @param outputB (Output) Symmetric matrix.
* @param gw (Optional) Workspace. Can be null.
*/
public static void symmLowerToFull( FMatrixSparseCSC A, FMatrixSparseCSC outputB, @Nullable IGrowArray gw ) {
ImplCommonOps_FSCC.symmLowerToFull(A, outputB, gw);
}
/**
* Performs matrix addition:
* C = αA + βB
*
* @param alpha scalar value multiplied against A
* @param A Matrix
* @param beta scalar value multiplied against B
* @param B Matrix
* @param outputC Output matrix.
* @param gw (Optional) Storage for internal workspace. Can be null.
* @param gx (Optional) Storage for internal workspace. Can be null.
*/
public static FMatrixSparseCSC add( float alpha, FMatrixSparseCSC A, float beta, FMatrixSparseCSC B,
@Nullable FMatrixSparseCSC outputC,
@Nullable IGrowArray gw, @Nullable FGrowArray gx ) {
if (A.numRows != B.numRows || A.numCols != B.numCols)
throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B));
outputC = reshapeOrDeclare(outputC, A, A.numRows, A.numCols);
ImplCommonOps_FSCC.add(alpha, A, beta, B, outputC, gw, gx);
return outputC;
}
public static FMatrixSparseCSC identity( int length ) {
return identity(length, length);
}
public static FMatrixSparseCSC identity( int numRows, int numCols ) {
int min = Math.min(numRows, numCols);
FMatrixSparseCSC A = new FMatrixSparseCSC(numRows, numCols, min);
setIdentity(A);
return A;
}
public static void setIdentity( FMatrixSparseCSC A ) {
int min = Math.min(A.numRows, A.numCols);
A.growMaxLength(min, false);
A.nz_length = min;
Arrays.fill(A.nz_values, 0, min, 1);
for (int i = 1; i <= min; i++) {
A.col_idx[i] = i;
A.nz_rows[i - 1] = i - 1;
}
for (int i = min + 1; i <= A.numCols; i++) {
A.col_idx[i] = min;
}
}
/**
* B = scalar*A. A and B can be the same instance.
*
* @param scalar (Input) Scalar value
* @param A (Input) Matrix. Not modified.
* @param outputB (Output) Matrix. Modified.
*/
public static void scale( float scalar, FMatrixSparseCSC A, FMatrixSparseCSC outputB ) {
if (A != outputB) {
outputB.copyStructure(A);
for (int i = 0; i < A.nz_length; i++) {
outputB.nz_values[i] = A.nz_values[i]*scalar;
}
} else {
for (int i = 0; i < A.nz_length; i++) {
outputB.nz_values[i] *= scalar;
}
}
}
/**
* B = A/scalar. A and B can be the same instance.
*
* @param scalar (Input) Scalar value
* @param A (Input) Matrix. Not modified.
* @param outputB (Output) Matrix. Modified.
*/
public static void divide( FMatrixSparseCSC A, float scalar, FMatrixSparseCSC outputB ) {
if (A != outputB) {
outputB.copyStructure(A);
for (int i = 0; i < A.nz_length; i++) {
outputB.nz_values[i] = A.nz_values[i]/scalar;
}
} else {
for (int i = 0; i < A.nz_length; i++) {
A.nz_values[i] /= scalar;
}
}
}
/**
* B = scalar/A. A and B can be the same instance. Only non-zero values are affected
*
* @param A (Input) Matrix. Not modified.
* @param scalar (Input) Scalar value
* @param outputB (Output) Matrix. Modified.
*/
public static void divide( float scalar, FMatrixSparseCSC A, FMatrixSparseCSC outputB ) {
if (A != outputB) {
outputB.copyStructure(A);
}
for (int i = 0; i < A.nz_length; i++) {
outputB.nz_values[i] = scalar/A.nz_values[i];
}
}
/**
* B = -A. Changes the sign of elements in A and stores it in B. A and B can be the same instance.
*
* @param A (Input) Matrix. Not modified.
* @param outputB (Output) Matrix. Modified.
*/
public static void changeSign( FMatrixSparseCSC A, FMatrixSparseCSC outputB ) {
if (A != outputB) {
outputB.copyStructure(A);
}
for (int i = 0; i < A.nz_length; i++) {
outputB.nz_values[i] = -A.nz_values[i];
}
}
/**
* Returns the value of the element with the smallest abs()
*
* @param A (Input) Matrix. Not modified.
* @return scalar
*/
public static float elementMinAbs( FMatrixSparseCSC A ) {
if (A.nz_length == 0)
return 0;
float min = A.isFull() ? Math.abs(A.nz_values[0]) : 0;
for (int i = 0; i < A.nz_length; i++) {
float val = Math.abs(A.nz_values[i]);
if (val < min) {
min = val;
}
}
return min;
}
/**
* Returns the value of the element with the largest abs()
*
* @param A (Input) Matrix. Not modified.
* @return scalar
*/
public static float elementMaxAbs( FMatrixSparseCSC A ) {
if (A.nz_length == 0)
return 0;
float max = A.isFull() ? Math.abs(A.nz_values[0]) : 0;
for (int i = 0; i < A.nz_length; i++) {
float val = Math.abs(A.nz_values[i]);
if (val > max) {
max = val;
}
}
return max;
}
/**
* Returns the value of the element with the minimum value
*
* @param A (Input) Matrix. Not modified.
* @return scalar
*/
public static float elementMin( FMatrixSparseCSC A ) {
if (A.nz_length == 0)
return 0;
// if every element is assigned a value then the first element can be a minimum.
// Otherwise zero needs to be considered
float min = A.isFull() ? A.nz_values[0] : 0;
for (int i = 0; i < A.nz_length; i++) {
float val = A.nz_values[i];
if (val < min) {
min = val;
}
}
return min;
}
/**
* Returns the value of the element with the largest value
*
* @param A (Input) Matrix. Not modified.
* @return scalar
*/
public static float elementMax( FMatrixSparseCSC A ) {
if (A.nz_length == 0)
return 0;
// if every element is assigned a value then the first element can be a max.
// Otherwise zero needs to be considered
float max = A.isFull() ? A.nz_values[0] : 0;
for (int i = 0; i < A.nz_length; i++) {
float val = A.nz_values[i];
if (val > max) {
max = val;
}
}
return max;
}
/**
* Sum of all elements
*
* @param A (Input) Matrix. Not modified.
* @return scalar
*/
public static float elementSum( FMatrixSparseCSC A ) {
if (A.nz_length == 0)
return 0;
float sum = 0;
for (int i = 0; i < A.nz_length; i++) {
sum += A.nz_values[i];
}
return sum;
}
/**
* Performs an element-wise multiplication.
* output[i,j] = A[i,j]*B[i,j]
* All matrices must have the same shape.
*
* @param A (Input) Matrix.
* @param B (Input) Matrix
* @param output (Output) Matrix. data array is grown to min(A.nz_length,B.nz_length), resulting a in a large speed boost.
* @param gw (Optional) Storage for internal workspace. Can be null.
* @param gx (Optional) Storage for internal workspace. Can be null.
*/
public static FMatrixSparseCSC elementMult( FMatrixSparseCSC A, FMatrixSparseCSC B, @Nullable FMatrixSparseCSC output,
@Nullable IGrowArray gw, @Nullable FGrowArray gx ) {
if (A.numCols != B.numCols || A.numRows != B.numRows)
throw new MatrixDimensionException("All inputs must have the same number of rows and columns. " + stringShapes(A, B));
output = reshapeOrDeclare(output, A, A.numRows, A.numCols);
ImplCommonOps_FSCC.elementMult(A, B, output, gw, gx);
return output;
}
/**
* Finds the maximum abs in each column of A and stores it into values
*
* @param A (Input) Matrix
* @param outputB (Output) storage for column max abs
*/
public static void maxAbsCols( FMatrixSparseCSC A, @Nullable FMatrixRMaj outputB ) {
outputB = reshapeOrDeclare(outputB, 1, A.numCols);
for (int i = 0; i < A.numCols; i++) {
int idx0 = A.col_idx[i];
int idx1 = A.col_idx[i + 1];
float maxabs = 0;
for (int j = idx0; j < idx1; j++) {
float v = Math.abs(A.nz_values[j]);
if (v > maxabs)
maxabs = v;
}
outputB.data[i] = maxabs;
}
}
/**
* Multiply all elements of column 'i' by value[i]. A[:,i] *= values[i].
* Equivalent to A = A*diag(values)
*
* @param A (Input/Output) Matrix. Modified.
* @param values (Input) multiplication factor for each column
* @param offset (Input) first index in values to start at
*/
public static void multColumns( FMatrixSparseCSC A, float[] values, int offset ) {
if (values.length + offset < A.numCols)
throw new IllegalArgumentException("Array is too small. " + values.length + " < " + A.numCols);
for (int i = 0; i < A.numCols; i++) {
int idx0 = A.col_idx[i];
int idx1 = A.col_idx[i + 1];
float v = values[offset + i];
for (int j = idx0; j < idx1; j++) {
A.nz_values[j] *= v;
}
}
}
/**
* Divides all elements of column 'i' by values[i]. A[:,i] /= values[i].
* Equivalent to A = A*inv(diag(values))
*
* @param A (Input/Output) Matrix. Modified.
* @param values (Input) multiplication factor for each column
* @param offset (Input) first index in values to start at
*/
public static void divideColumns( FMatrixSparseCSC A, float[] values, int offset ) {
if (values.length + offset < A.numCols)
throw new IllegalArgumentException("Array is too small. " + values.length + " < " + A.numCols);
for (int i = 0; i < A.numCols; i++) {
int idx0 = A.col_idx[i];
int idx1 = A.col_idx[i + 1];
float v = values[offset + i];
for (int j = idx0; j < idx1; j++) {
A.nz_values[j] /= v;
}
}
}
/**
* Equivalent to multiplying a matrix B by two diagonal matrices.
* B = A*B*C, where A=diag(a) and C=diag(c).
*
* @param diagA Array of length offsteA + B.numRows
* @param offsetA First index in A
* @param B Rectangular matrix
* @param diagC Array of length indexC + B.numCols
* @param offsetC First index in C
*/
public static void multRowsCols( float[] diagA, int offsetA,
FMatrixSparseCSC B,
float[] diagC, int offsetC ) {
if (diagA.length + offsetA < B.numRows)
throw new IllegalArgumentException("diagA is too small.");
if (diagC.length + offsetC < B.numCols)
throw new IllegalArgumentException("diagA is too small.");
for (int i = 0; i < B.numCols; i++) {
int idx0 = B.col_idx[i];
int idx1 = B.col_idx[i + 1];
float c = diagC[offsetC + i];
for (int j = idx0; j < idx1; j++) {
B.nz_values[j] *= c*diagA[offsetA + B.nz_rows[j]];
}
}
}
/**
* Equivalent to multiplying a matrix B by the inverse of two diagonal matrices.
* B = inv(A)*B*inv(C), where A=diag(a) and C=diag(c).
*
* @param diagA Array of length offsteA + B.numRows
* @param offsetA First index in A
* @param B Rectangular matrix
* @param diagC Array of length indexC + B.numCols
* @param offsetC First index in C
*/
public static void divideRowsCols( float[] diagA, int offsetA,
FMatrixSparseCSC B,
float[] diagC, int offsetC ) {
if (diagA.length + offsetA < B.numRows)
throw new IllegalArgumentException("diagA is too small.");
if (diagC.length + offsetC < B.numCols)
throw new IllegalArgumentException("diagA is too small.");
for (int i = 0; i < B.numCols; i++) {
int idx0 = B.col_idx[i];
int idx1 = B.col_idx[i + 1];
float c = diagC[offsetC + i];
for (int j = idx0; j < idx1; j++) {
B.nz_values[j] /= c*diagA[offsetA + B.nz_rows[j]];
}
}
}
/**
* Returns a diagonal matrix with the specified diagonal elements.
*
* @param values values of diagonal elements
* @return A diagonal matrix
*/
public static FMatrixSparseCSC diag( float... values ) {
int N = values.length;
return diag(new FMatrixSparseCSC(N, N, N), values, 0, N);
}
/**
* Creates a diagonal matrix from an array. Elements in the array can be offset.
*
* @param A (Optional) Storage for diagonal matrix. If null a new one will be declared.
* @param values First index in the diagonal matirx
* @param length Length of the diagonal matrix
* @param offset First index in values
* @return The diagonal matrix
*/
public static FMatrixSparseCSC diag( @Nullable FMatrixSparseCSC A,
float[] values, int offset, int length ) {
int N = length;
if (A == null)
A = new FMatrixSparseCSC(N, N, N);
else
A.reshape(N, N, N);
A.nz_length = N;
for (int i = 0; i < N; i++) {
A.col_idx[i + 1] = i + 1;
A.nz_rows[i] = i;
A.nz_values[i] = values[offset + i];
}
return A;
}
/**
*
* Extracts the diagonal elements 'src' write it to the 'dst' vector. 'dst'
* can either be a row or column vector.
*
*
* @param A Matrix whose diagonal elements are being extracted. Not modified.
* @param outputB A vector the results will be written into. Modified.
*/
public static void extractDiag( FMatrixSparseCSC A, FMatrixSparseCSC outputB ) {
int N = Math.min(A.numRows, A.numCols);
if (!MatrixFeatures_FSCC.isVector(outputB)) {
outputB.reshape(N, 1, N);
} else if (outputB.numRows*outputB.numCols != N) {
outputB.reshape(N, 1, N);
} else {
outputB.growMaxLength(N, false);
}
outputB.nz_length = N;
outputB.indicesSorted = true;
if (outputB.numRows != 1) {
outputB.col_idx[0] = 0;
outputB.col_idx[1] = N;
for (int i = 0; i < N; i++) {
outputB.nz_values[i] = A.unsafe_get(i, i);
outputB.nz_rows[i] = i;
}
} else {
outputB.col_idx[0] = 0;
for (int i = 0; i < N; i++) {
outputB.nz_values[i] = A.unsafe_get(i, i);
outputB.nz_rows[i] = 0;
outputB.col_idx[i + 1] = i + 1;
}
}
}
/**
*
* Extracts the diagonal elements 'src' write it to the 'dst' vector. 'dst'
* can either be a row or column vector.
*
*
* @param A Matrix whose diagonal elements are being extracted. Not modified.
* @param outputB A vector the results will be written into. Modified.
*/
public static void extractDiag( FMatrixSparseCSC A, FMatrixRMaj outputB ) {
int N = Math.min(A.numRows, A.numCols);
if (outputB.getNumElements() != N || !(outputB.numRows == 1 || outputB.numCols == 1)) {
outputB.reshape(N, 1);
}
for (int i = 0; i < N; i++) {
outputB.data[i] = A.unsafe_get(i, i);
}
}
/**
* Converts the permutation vector into a matrix. B = P*A. B[p[i],:] = A[i,:]
*
* @param p (Input) Permutation vector
* @param inverse (Input) If it is the inverse. B[i,:] = A[p[i],:)
* @param P (Output) Permutation matrix
*/
public static FMatrixSparseCSC permutationMatrix( int[] p, boolean inverse, int N,
@Nullable FMatrixSparseCSC P ) {
if (P == null)
P = new FMatrixSparseCSC(N, N, N);
else
P.reshape(N, N, N);
P.indicesSorted = true;
P.nz_length = N;
// each column should have one element inside of it
if (!inverse) {
for (int i = 0; i < N; i++) {
P.col_idx[i + 1] = i + 1;
P.nz_rows[p[i]] = i;
P.nz_values[i] = 1;
}
} else {
for (int i = 0; i < N; i++) {
P.col_idx[i + 1] = i + 1;
P.nz_rows[i] = p[i];
P.nz_values[i] = 1;
}
}
return P;
}
/**
* Converts the permutation matrix into a vector
*
* @param P (Input) Permutation matrix
* @param vector (Output) Permutation vector
*/
public static void permutationVector( FMatrixSparseCSC P, int[] vector ) {
if (P.numCols != P.numRows) {
throw new MatrixDimensionException("Expected a square matrix");
} else if (P.nz_length != P.numCols) {
throw new IllegalArgumentException("Expected N non-zero elements in permutation matrix");
} else if (vector.length < P.numCols) {
throw new IllegalArgumentException("vector is too short");
}
int M = P.numCols;
for (int i = 0; i < M; i++) {
if (P.col_idx[i + 1] != i + 1)
throw new IllegalArgumentException("Unexpected number of elements in a column");
vector[P.nz_rows[i]] = i;
}
}
/**
* Computes the inverse permutation vector
*
* @param original Original permutation vector
* @param inverse It's inverse
*/
public static void permutationInverse( int[] original, int[] inverse, int length ) {
for (int i = 0; i < length; i++) {
inverse[original[i]] = i;
}
}
public static int[] permutationInverse( int[] original, int length ) {
int[] inverse = new int[length];
permutationInverse(original, inverse, length);
return inverse;
}
/**
* Applies the row permutation specified by the vector to the input matrix and save the results
* in the output matrix. output[perm[j],:] = input[j,:]
*
* @param permInv (Input) Inverse permutation vector. Specifies new order of the rows.
* @param input (Input) Matrix which is to be permuted
* @param output (Output) Matrix which has the permutation stored in it. Is reshaped.
*/
public static void permuteRowInv( int[] permInv, FMatrixSparseCSC input, FMatrixSparseCSC output ) {
if (input.numRows > permInv.length)
throw new IllegalArgumentException("permutation vector must have at least as many elements as input has rows");
output.reshape(input.numRows, input.numCols, input.nz_length);
output.nz_length = input.nz_length;
output.indicesSorted = false;
System.arraycopy(input.nz_values, 0, output.nz_values, 0, input.nz_length);
System.arraycopy(input.col_idx, 0, output.col_idx, 0, input.numCols + 1);
int idx0 = 0;
for (int i = 0; i < input.numCols; i++) {
int idx1 = output.col_idx[i + 1];
for (int j = idx0; j < idx1; j++) {
output.nz_rows[j] = permInv[input.nz_rows[j]];
}
idx0 = idx1;
}
}
/**
* Applies the forward column and inverse row permutation specified by the two vector to the input matrix
* and save the results in the output matrix. output[permRow[j],permCol[i]] = input[j,i]
*
* @param permRowInv (Input) Inverse row permutation vector. Null is the same as passing in identity.
* @param input (Input) Matrix which is to be permuted
* @param permCol (Input) Column permutation vector. Null is the same as passing in identity.
* @param output (Output) Matrix which has the permutation stored in it. Is reshaped.
*/
public static void permute( @Nullable int[] permRowInv, FMatrixSparseCSC input, @Nullable int[] permCol,
FMatrixSparseCSC output ) {
if (permRowInv != null && input.numRows > permRowInv.length)
throw new IllegalArgumentException("rowInv permutation vector must have at least as many elements as input has columns");
if (permCol != null && input.numCols > permCol.length)
throw new IllegalArgumentException("permCol permutation vector must have at least as many elements as input has rows");
output.reshape(input.numRows, input.numCols, input.nz_length);
output.indicesSorted = false;
output.nz_length = input.nz_length;
int N = input.numCols;
// traverse through in order for the output columns
int outputNZ = 0;
for (int i = 0; i < N; i++) {
int inputCol = permCol != null ? permCol[i] : i; // column of input to source from
int inputNZ = input.col_idx[inputCol];
int total = input.col_idx[inputCol + 1] - inputNZ; // total nz in this column
output.col_idx[i + 1] = output.col_idx[i] + total;
for (int j = 0; j < total; j++) {
int row = input.nz_rows[inputNZ];
output.nz_rows[outputNZ] = permRowInv != null ? permRowInv[row] : row;
output.nz_values[outputNZ++] = input.nz_values[inputNZ++];
}
}
}
/**
* Permutes a vector. output[i] = input[perm[i]]
*
* @param perm (Input) permutation vector
* @param input (Input) Vector which is to be permuted
* @param output (Output) Where the permuted vector is stored.
* @param N Number of elements in the vector.
*/
public static void permute( int[] perm, float[] input, float[] output, int N ) {
for (int k = 0; k < N; k++) {
output[k] = input[perm[k]];
}
}
/**
* Permutes a vector in the inverse. output[perm[k]] = input[k]
*
* @param perm (Input) permutation vector
* @param input (Input) Vector which is to be permuted
* @param output (Output) Where the permuted vector is stored.
* @param N Number of elements in the vector.
*/
public static void permuteInv( int[] perm, float[] input, float[] output, int N ) {
for (int k = 0; k < N; k++) {
output[perm[k]] = input[k];
}
}
/**
* Applies the permutation to upper triangular symmetric matrices. Typically a symmetric matrix only stores the
* upper triangular part, so normal permutation will have undesirable results, e.g. the zeros will get mixed
* in and will no longer be symmetric. This algorithm will handle the implicit lower triangular and construct
* new upper triangular matrix.
*
*
See page cs_symperm() on Page 22 of "Direct Methods for Sparse Linear Systems"
*
* @param input (Input) Upper triangular symmetric matrix which is to be permuted.
* Entries below the diagonal are ignored.
* @param permInv (Input) Inverse permutation vector. Specifies new order of the rows and columns.
* @param output (Output) Upper triangular symmetric matrix which has the permutation stored in it. Reshaped.
* @param gw (Optional) Storage for internal workspace. Can be null.
*/
public static void permuteSymmetric( FMatrixSparseCSC input, int[] permInv, FMatrixSparseCSC output,
@Nullable IGrowArray gw ) {
if (input.numRows != input.numCols)
throw new MatrixDimensionException("Input must be a square matrix. " + stringShapes(input, output));
if (input.numRows != permInv.length)
throw new MatrixDimensionException("Number of column in input must match length of permInv");
int N = input.numCols;
int[] w = adjustClear(gw, N); // histogram with column counts
output.reshape(N, N, 0);
output.indicesSorted = false;
output.col_idx[0] = 0;
// determine column counts for output
for (int j = 0; j < N; j++) {
int j2 = permInv[j];
int idx0 = input.col_idx[j];
int idx1 = input.col_idx[j + 1];
for (int p = idx0; p < idx1; p++) {
int i = input.nz_rows[p];
if (i > j) // ignore the lower triangular portion
continue;
int i2 = permInv[i];
w[i2 > j2 ? i2 : j2]++;
}
}
// update structure of output
output.histogramToStructure(w);
System.arraycopy(output.col_idx, 0, w, 0, output.numCols);
for (int j = 0; j < N; j++) {
// column j of Input is row j2 of Output
int j2 = permInv[j];
int idx0 = input.col_idx[j];
int idx1 = input.col_idx[j + 1];
for (int p = idx0; p < idx1; p++) {
int i = input.nz_rows[p];
if (i > j) // ignore the lower triangular portion
continue;
int i2 = permInv[i];
// row i of Input is row i2 of Output
int q = w[i2 > j2 ? i2 : j2]++;
output.nz_rows[q] = i2 < j2 ? i2 : j2;
output.nz_values[q] = input.nz_values[p];
}
}
}
/**
* Concats two matrices along their rows (vertical).
*
* @param top Matrix on the top
* @param bottom Matrix on the bototm
* @param out (Output) (Optional) Storage for combined matrix. Resized.
* @return Combination of the two matrices
*/
public static FMatrixSparseCSC concatRows( FMatrixSparseCSC top, FMatrixSparseCSC bottom,
@Nullable FMatrixSparseCSC out ) {
if (top.numCols != bottom.numCols)
throw new MatrixDimensionException("Number of columns must match. " + stringShapes(top, bottom));
if (out == null)
out = new FMatrixSparseCSC(0, 0, 0);
out.reshape(top.numRows + bottom.numRows, top.numCols, top.nz_length + bottom.nz_length);
out.nz_length = top.nz_length + bottom.nz_length;
int index = 0;
for (int i = 0; i < top.numCols; i++) {
int top0 = top.col_idx[i];
int top1 = top.col_idx[i + 1];
int bot0 = bottom.col_idx[i];
int bot1 = bottom.col_idx[i + 1];
int out0 = out.col_idx[i];
int out1 = out0 + top1 - top0 + bot1 - bot0;
out.col_idx[i + 1] = out1;
for (int j = top0; j < top1; j++, index++) {
out.nz_values[index] = top.nz_values[j];
out.nz_rows[index] = top.nz_rows[j];
}
for (int j = bot0; j < bot1; j++, index++) {
out.nz_values[index] = bottom.nz_values[j];
out.nz_rows[index] = top.numRows + bottom.nz_rows[j];
}
}
out.indicesSorted = false;
return out;
}
/**
* Concats two matrices along their columns (horizontal).
*
* @param left Matrix on the left
* @param right Matrix on the right
* @param out (Output) (Optional) Storage for combined matrix. Resized.
* @return Combination of the two matrices
*/
public static FMatrixSparseCSC concatColumns( FMatrixSparseCSC left, FMatrixSparseCSC right,
@Nullable FMatrixSparseCSC out ) {
if (left.numRows != right.numRows)
throw new MatrixDimensionException("Number of rows must match. " + stringShapes(left, right));
if (out == null)
out = new FMatrixSparseCSC(0, 0, 0);
out.reshape(left.numRows, left.numCols + right.numCols, left.nz_length + right.nz_length);
out.nz_length = left.nz_length + right.nz_length;
System.arraycopy(left.col_idx, 0, out.col_idx, 0, left.numCols + 1);
System.arraycopy(left.nz_rows, 0, out.nz_rows, 0, left.nz_length);
System.arraycopy(left.nz_values, 0, out.nz_values, 0, left.nz_length);
int index = left.nz_length;
for (int i = 0; i < right.numCols; i++) {
int r0 = right.col_idx[i];
int r1 = right.col_idx[i + 1];
out.col_idx[left.numCols + i] = index;
out.col_idx[left.numCols + i + 1] = index + (r1 - r0);
for (int j = r0; j < r1; j++, index++) {
out.nz_rows[index] = right.nz_rows[j];
out.nz_values[index] = right.nz_values[j];
}
}
out.indicesSorted = left.indicesSorted && right.indicesSorted;
return out;
}
/**
* Extracts a column from A and stores it into out.
*
* @param A (Input) Source matrix. not modified.
* @param column The column in A
* @param out (Output, Optional) Storage for column vector
* @return The column of A.
*/
public static FMatrixSparseCSC extractColumn( FMatrixSparseCSC A, int column, @Nullable FMatrixSparseCSC out ) {
if (out == null)
out = new FMatrixSparseCSC(1, 1, 1);
int idx0 = A.col_idx[column];
int idx1 = A.col_idx[column + 1];
out.reshape(A.numRows, 1, idx1 - idx0);
out.nz_length = idx1 - idx0;
out.col_idx[0] = 0;
out.col_idx[1] = out.nz_length;
System.arraycopy(A.nz_values, idx0, out.nz_values, 0, out.nz_length);
System.arraycopy(A.nz_rows, idx0, out.nz_rows, 0, out.nz_length);
return out;
}
/**
* Creates a submatrix by extracting the specified rows from A. rows = {row0 %le; i %le; row1}.
*
* @param A (Input) matrix
* @param row0 First row. Inclusive
* @param row1 Last row+1.
* @param out (Output, Option) Storage for output matrix
* @return The submatrix
*/
public static FMatrixSparseCSC extractRows( FMatrixSparseCSC A, int row0, int row1,
@Nullable FMatrixSparseCSC out ) {
if (out == null)
out = new FMatrixSparseCSC(1, 1, 1);
out.reshape(row1 - row0, A.numCols, A.nz_length);
// out.col_idx[0] = 0;
// out.nz_length = 0;
for (int col = 0; col < A.numCols; col++) {
int idx0 = A.col_idx[col];
int idx1 = A.col_idx[col + 1];
for (int i = idx0; i < idx1; i++) {
int row = A.nz_rows[i];
if (row >= row0 && row < row1) {
out.nz_values[out.nz_length] = A.nz_values[i];
out.nz_rows[out.nz_length++] = row - row0;
}
}
out.col_idx[col + 1] = out.nz_length;
}
return out;
}
/**
*
* Extracts a submatrix from 'src' and inserts it in a submatrix in 'dst'.
*
*
* si-y0 , j-x0 = oij for all y0 ≤ i < y1 and x0 ≤ j < x1
*
* where 'sij' is an element in the submatrix and 'oij' is an element in the
* original matrix.
*
*
* WARNING: This is a very slow operation for sparse matrices. The current implementation is simple but
* involves excessive memory copying.
*
* @param src The original matrix which is to be copied. Not modified.
* @param srcX0 Start column.
* @param srcX1 Stop column+1.
* @param srcY0 Start row.
* @param srcY1 Stop row+1.
* @param dst Where the submatrix are stored. Modified.
* @param dstY0 Start row in dst.
* @param dstX0 start column in dst.
*/
public static void extract( FMatrixSparseCSC src, int srcY0, int srcY1, int srcX0, int srcX1,
FMatrixSparseCSC dst, int dstY0, int dstX0 ) {
if (srcY1 < srcY0 || srcY0 < 0 || srcY1 > src.getNumRows())
throw new MatrixDimensionException("srcY1 < srcY0 || srcY0 < 0 || srcY1 > src.numRows. " + stringShapes(src, dst));
if (srcX1 < srcX0 || srcX0 < 0 || srcX1 > src.getNumCols())
throw new MatrixDimensionException("srcX1 < srcX0 || srcX0 < 0 || srcX1 > src.numCols. " + stringShapes(src, dst));
int w = srcX1 - srcX0;
int h = srcY1 - srcY0;
if (dstY0 + h > dst.getNumRows())
throw new IllegalArgumentException("dst is too small in rows. " + dst.getNumRows() + " < " + (dstY0 + h));
if (dstX0 + w > dst.getNumCols())
throw new IllegalArgumentException("dst is too small in columns. " + dst.getNumCols() + " < " + (dstX0 + w));
zero(dst, dstY0, dstY0 + h, dstX0, dstX0 + w);
// NOTE: One possible optimization would be to determine the non-zero pattern in dst after the change is
// applied, modify it's structure, then copy the values in. That way you aren't shifting memory constantly.
//
// NOTE: Another optimization would be to sort the src so that it doesn't need to go through every row
for (int colSrc = srcX0; colSrc < srcX1; colSrc++) {
int idxS0 = src.col_idx[colSrc];
int idxS1 = src.col_idx[colSrc + 1];
for (int i = idxS0; i < idxS1; i++) {
int row = src.nz_rows[i];
if (row >= srcY0 && row < srcY1) {
dst.set(row - srcY0 + dstY0, colSrc - srcX0 + dstX0, src.nz_values[i]);
}
}
}
}
/**
* Select entries from A and save them in C.
* A more generic but probably also slower version of `extract`
*
* Can be used to select f.i. the lower or upper triangle of a matrix
*
* Simplified version of: GxB_Select
*
* @param A (Input) Matrix. Not modified.
* @param selector Function to decide whether an entry gets selected
* @param output (Optional/Output) Matrix to use for the output. Can be the same as A
* @return Matrix storing the selected entries of A
*/
public static FMatrixSparseCSC select( FMatrixSparseCSC A, IPredicateBinary selector, @Nullable FMatrixSparseCSC output) {
if (output != A) {
output = reshapeOrDeclare(output, A);
}
ImplCommonOps_FSCC.select(A, output, selector);
return output;
}
/**
*
* Sets every element in the matrix to the specified value. This can require a very large amount of
* memory and might exceed the maximum array size
*
* Aij = value
*
*
* @param A A matrix whose elements are about to be set. Modified.
* @param value The value each element will have.
*/
public static void fill( FMatrixSparseCSC A, float value ) {
int N = A.numCols*A.numRows;
A.growMaxLength(N, false);
A.col_idx[0] = 0;
for (int col = 0; col < A.numCols; col++) {
int idx0 = A.col_idx[col];
int idx1 = A.col_idx[col + 1] = idx0 + A.numRows;
for (int i = idx0; i < idx1; i++) {
A.nz_rows[i] = i - idx0;
A.nz_values[i] = value;
}
}
A.nz_length = N;
A.indicesSorted = true;
}
/**
*
* Computes the sum of each column in the input matrix and returns the results in a vector:
*
* bj = sum(i=1:m ; aij)
*
*
* @param input Input matrix
* @param output Optional storage for output. Reshaped into a row vector. Modified.
* @return Vector containing the sum of each column
*/
public static FMatrixRMaj sumCols( FMatrixSparseCSC input, @Nullable FMatrixRMaj output ) {
if (output == null) {
output = new FMatrixRMaj(1, input.numCols);
} else {
output.reshape(1, input.numCols);
}
for (int col = 0; col < input.numCols; col++) {
int idx0 = input.col_idx[col];
int idx1 = input.col_idx[col + 1];
float sum = 0;
for (int i = idx0; i < idx1; i++) {
sum += input.nz_values[i];
}
output.data[col] = sum;
}
return output;
}
/**
*
* Computes the minimum of each column in the input matrix and returns the results in a vector:
*
* bj = min(i=1:m ; aij)
*
*
* @param input Input matrix
* @param output Optional storage for output. Reshaped into a row vector. Modified.
* @return Vector containing the minimums of each column
*/
public static FMatrixRMaj minCols( FMatrixSparseCSC input, @Nullable FMatrixRMaj output ) {
output = reshapeOrDeclare(output, 1, input.numCols);
for (int col = 0; col < input.numCols; col++) {
int idx0 = input.col_idx[col];
int idx1 = input.col_idx[col + 1];
// take in account the zeros
float min = idx1 - idx0 == input.numRows ? Float.MAX_VALUE : 0;
for (int i = idx0; i < idx1; i++) {
float v = input.nz_values[i];
if (min > v) {
min = v;
}
}
output.data[col] = min;
}
return output;
}
/**
*
* Computes the maximums of each column in the input matrix and returns the results in a vector:
*
* bj = max(i=1:m ; aij)
*
*
* @param input Input matrix
* @param output Optional storage for output. Reshaped into a row vector. Modified.
* @return Vector containing the maximums of each column
*/
public static FMatrixRMaj maxCols( FMatrixSparseCSC input, @Nullable FMatrixRMaj output ) {
output = reshapeOrDeclare(output, 1, input.numCols);
for (int col = 0; col < input.numCols; col++) {
int idx0 = input.col_idx[col];
int idx1 = input.col_idx[col + 1];
// take in account the zeros
float max = idx1 - idx0 == input.numRows ? -Float.MAX_VALUE : 0;
for (int i = idx0; i < idx1; i++) {
float v = input.nz_values[i];
if (max < v) {
max = v;
}
}
output.data[col] = max;
}
return output;
}
/**
*
* Computes the sum of each row in the input matrix and returns the results in a vector:
*
* bj = sum(i=1:n ; aji)
*
*
* @param input Input matrix
* @param output Optional storage for output. Reshaped into a column vector. Modified.
* @return Vector containing the sum of each row
*/
public static FMatrixRMaj sumRows( FMatrixSparseCSC input, @Nullable FMatrixRMaj output ) {
output = reshapeOrDeclare(output, input.numRows, 1);
Arrays.fill(output.data, 0, input.numRows, 0);
for (int col = 0; col < input.numCols; col++) {
int idx0 = input.col_idx[col];
int idx1 = input.col_idx[col + 1];
for (int i = idx0; i < idx1; i++) {
output.data[input.nz_rows[i]] += input.nz_values[i];
}
}
return output;
}
/**
*
* Computes the minimum of each row in the input matrix and returns the results in a vector:
*
* bj = min(i=1:n ; aji)
*
*
* @param input Input matrix
* @param output Optional storage for output. Reshaped into a column vector. Modified.
* @param gw work space
* @return Vector containing the minimum of each row
*/
public static FMatrixRMaj minRows( FMatrixSparseCSC input, @Nullable FMatrixRMaj output, @Nullable IGrowArray gw ) {
output = reshapeOrDeclare(output, input.numRows, 1);
int[] w = adjust(gw, input.numRows, input.numRows);
Arrays.fill(output.data, 0, input.numRows, Float.MAX_VALUE);
for (int col = 0; col < input.numCols; col++) {
int idx0 = input.col_idx[col];
int idx1 = input.col_idx[col + 1];
for (int i = idx0; i < idx1; i++) {
int row = input.nz_rows[i];
float v = input.nz_values[i];
if (output.data[row] > v) {
output.data[row] = v;
}
w[row]++;
}
}
for (int row = 0; row < input.numRows; row++) {
// consider the zeros now if a row wasn't filled in all the way
if (w[row] != input.numCols) {
if (output.data[row] > 0) {
output.data[row] = 0;
}
}
}
return output;
}
/**
*
* Computes the maximum of each row in the input matrix and returns the results in a vector:
*
* bj = max(i=1:n ; aji)
*
*
* @param input Input matrix
* @param output Optional storage for output. Reshaped into a column vector. Modified.
* @param gw work space
* @return Vector containing the maximum of each row
*/
public static FMatrixRMaj maxRows( FMatrixSparseCSC input, @Nullable FMatrixRMaj output, @Nullable IGrowArray gw ) {
output = reshapeOrDeclare(output, input.numRows, 1);
int[] w = adjust(gw, input.numRows, input.numRows);
Arrays.fill(output.data, 0, input.numRows, -Float.MAX_VALUE);
for (int col = 0; col < input.numCols; col++) {
int idx0 = input.col_idx[col];
int idx1 = input.col_idx[col + 1];
for (int i = idx0; i < idx1; i++) {
int row = input.nz_rows[i];
float v = input.nz_values[i];
if (output.data[row] < v) {
output.data[row] = v;
}
w[row]++;
}
}
for (int row = 0; row < input.numRows; row++) {
// consider the zeros now if a row wasn't filled in all the way
if (w[row] != input.numCols) {
if (output.data[row] < 0) {
output.data[row] = 0;
}
}
}
return output;
}
/**
* Zeros an inner rectangle inside the matrix.
*
* @param A Matrix that is to be modified.
* @param row0 Start row.
* @param row1 Stop row+1.
* @param col0 Start column.
* @param col1 Stop column+1.
*/
public static void zero( FMatrixSparseCSC A, int row0, int row1, int col0, int col1 ) {
for (int col = col1 - 1; col >= col0; col--) {
int numRemoved = 0;
int idx0 = A.col_idx[col], idx1 = A.col_idx[col + 1];
for (int i = idx0; i < idx1; i++) {
int row = A.nz_rows[i];
// if sorted a faster technique could be used
if (row >= row0 && row < row1) {
numRemoved++;
} else if (numRemoved > 0) {
A.nz_rows[i - numRemoved] = row;
A.nz_values[i - numRemoved] = A.nz_values[i];
}
}
if (numRemoved > 0) {
// this could be done more intelligently. Each time a column is adjusted all the columns are adjusted
// after it. Maybe accumulate the changes in each column and do it in one pass? Need an array to store
// those results though
for (int i = idx1; i < A.nz_length; i++) {
A.nz_rows[i - numRemoved] = A.nz_rows[i];
A.nz_values[i - numRemoved] = A.nz_values[i];
}
A.nz_length -= numRemoved;
for (int i = col + 1; i <= A.numCols; i++) {
A.col_idx[i] -= numRemoved;
}
}
}
}
/**
* Computes the inner product of two column vectors taken from the input matrices.
*
* dot = A(:,colA)'*B(:,colB)
*
* @param A Matrix
* @param colA Column in A
* @param B Matrix
* @param colB Column in B
* @return Dot product
*/
public static float dotInnerColumns( FMatrixSparseCSC A, int colA, FMatrixSparseCSC B, int colB,
@Nullable IGrowArray gw, @Nullable FGrowArray gx ) {
return ImplMultiplication_FSCC.dotInnerColumns(A, colA, B, colB, gw, gx);
}
/**
*
* Solves for x in the following equation:
*
* A*x = b
*
*
*
* If the system could not be solved then false is returned. If it returns true
* that just means the algorithm finished operating, but the results could still be bad
* because 'A' is singular or nearly singular.
*
*
*
* If repeat calls to solve are being made then one should consider using {@link LinearSolverFactory_FSCC}
* instead.
*
*
*
* It is ok for 'b' and 'x' to be the same matrix.
*
*
* @param a (Input) A matrix that is m by n. Not modified.
* @param b (Input) A matrix that is n by k. Not modified.
* @param x (Output) A matrix that is m by k. Modified.
* @return true if it could invert the matrix false if it could not.
*/
public static boolean solve( FMatrixSparseCSC a,
FMatrixRMaj b,
FMatrixRMaj x ) {
x.reshape(a.numCols, b.numCols);
LinearSolverSparse solver;
if (a.numRows > a.numCols) {
solver = LinearSolverFactory_FSCC.qr(FillReducing.NONE);// todo specify a filling that makes sense
} else {
solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE);
}
// Ensure that the input isn't modified
if (solver.modifiesA())
a = a.copy();
if (solver.modifiesB())
b = b.copy();
// decompose then solve the matrix
if (!solver.setA(a))
return false;
solver.solve(b, x);
return true;
}
/**
*
* Solves for x in the following equation:
*
* A*x = b
*
*
*
* If the system could not be solved then false is returned. If it returns true
* that just means the algorithm finished operating, but the results could still be bad
* because 'A' is singular or nearly singular.
*
*
*
* If repeat calls to solve are being made then one should consider using {@link LinearSolverFactory_FSCC}
* instead.
*
*
*
* It is ok for 'b' and 'x' to be the same matrix.
*
*
* @param a (Input) A matrix that is m by n. Not modified.
* @param b (Input) A matrix that is n by k. Not modified.
* @param x (Output) A matrix that is m by k. Modified.
* @return true if it could invert the matrix false if it could not.
*/
public static boolean solve( FMatrixSparseCSC a,
FMatrixSparseCSC b,
FMatrixSparseCSC x ) {
x.reshape(a.numCols, b.numCols);
LinearSolverSparse solver;
if (a.numRows > a.numCols) {
solver = LinearSolverFactory_FSCC.qr(FillReducing.NONE);// todo specify a filling that makes sense
} else {
solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE);
}
// Ensure that the input isn't modified
if (solver.modifiesA())
a = a.copy();
if (solver.modifiesB())
b = b.copy();
// decompose then solve the matrix
if (!solver.setA(a))
return false;
solver.solveSparse(b, x);
return true;
}
/**
*
* Performs a matrix inversion operation that does not modify the original
* and stores the results in another matrix. The two matrices must have the
* same dimension.
*
* B = A-1
*
*
*
* If the algorithm could not invert the matrix then false is returned. If it returns true
* that just means the algorithm finished. The results could still be bad
* because the matrix is singular or nearly singular.
*
*
*
* For medium to large matrices there might be a slight performance boost to using
* {@link LinearSolverFactory_FSCC} instead.
*
*
* @param A (Input) The matrix that is to be inverted. Not modified.
* @param inverse (Output) Where the inverse matrix is stored. Modified.
* @return true if it could invert the matrix false if it could not.
*/
public static boolean invert( FMatrixSparseCSC A, FMatrixRMaj inverse ) {
if (A.numRows != A.numCols)
throw new IllegalArgumentException("A must be a square matrix");
inverse.reshape(A.numRows, A.numCols);
LinearSolverSparse solver;
solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE);
// Ensure that the input isn't modified
if (solver.modifiesA())
A = A.copy();
FMatrixRMaj I = CommonOps_FDRM.identity(A.numRows);
// decompose then solve the matrix
if (!solver.setA(A))
return false;
solver.solve(I, inverse);
return true;
}
/**
* Returns the determinant of the matrix. If the inverse of the matrix is also
* needed, then using {@link org.ejml.interfaces.decomposition.LUDecomposition_F32} directly (or any
* similar algorithm) can be more efficient.
*
* @param A The matrix whose determinant is to be computed. Not modified.
* @return The determinant.
*/
public static float det( FMatrixSparseCSC A ) {
LUSparseDecomposition_F32 alg = DecompositionFactory_FSCC.lu(FillReducing.NONE);
if (alg.inputModified()) {
A = A.copy();
}
if (!alg.decompose(A))
return 0.0f;
return alg.computeDeterminant().real;
}
/**
* Copies all elements from input into output which are > tol.
*
* @param input (Input) input matrix. Not modified.
* @param output (Output) Output matrix. Modified and shaped to match input.
* @param tol Tolerance for defining zero
*/
public static void removeZeros( FMatrixSparseCSC input, FMatrixSparseCSC output, float tol ) {
ImplCommonOps_FSCC.removeZeros(input, output, tol);
}
/**
* Removes all elements from the matrix that are > tol. The modification is done in place and no temporary
* storage is declared.
*
* @param A (Input/Output) input matrix. Modified.
* @param tol Tolerance for defining zero
*/
public static void removeZeros( FMatrixSparseCSC A, float tol ) {
ImplCommonOps_FSCC.removeZeros(A, tol);
}
/**
* For each duplicate element in the matrix it will remove the duplicates and replace them with a single element
* that is the sum of all the duplicates. The result will be a valid matrix. This is done inplace and does not
* require the matrix to be initially sorted.
*
* @param A (Input/Output) input matrix. Modified.
* @param work Nullable. Internal workspace array.
*/
public static void duplicatesAdd( FMatrixSparseCSC A, @Nullable IGrowArray work ) {
ImplCommonOps_FSCC.duplicatesAdd(A, work);
}
/**
* Multiply all elements of row 'i' by value[i]. A[i,:] *= values[i]
*
* @param diag (Input) multiplication factors
* @param offset (Input) First index in values
* @param A (Input/Output) Matrix. Modified.
*/
public static void multRows( float[] diag, int offset, FMatrixSparseCSC A ) {
if (diag.length < A.numRows)
throw new IllegalArgumentException("Array is too small. " + diag.length + " < " + A.numCols);
for (int i = 0; i < A.nz_length; i++) {
A.nz_values[i] *= diag[A.nz_rows[i + offset]];
}
}
/**
* Divides all elements of row 'i' by value[i]. A[i,:] /= values[i]
*
* @param diag (Input) division factors
* @param offset (Input) First index in values
* @param A (Input/Output) Matrix. Modified.
*/
public static void divideRows( float[] diag, int offset, FMatrixSparseCSC A ) {
if (diag.length < A.numRows)
throw new IllegalArgumentException("Array is too small. " + diag.length + " < " + A.numCols);
for (int i = 0; i < A.nz_length; i++) {
A.nz_values[i] /= diag[A.nz_rows[i + offset]];
}
}
/**
*
* This computes the trace of the matrix:
*
* trace = ∑i=1:n { aii }
* where n = min(numRows,numCols)
*
*
* @param A (Input) Matrix. Not modified.
*/
public static float trace( FMatrixSparseCSC A ) {
float output = 0;
int o = Math.min(A.numCols, A.numRows);
for (int col = 0; col < o; col++) {
int idx0 = A.col_idx[col];
int idx1 = A.col_idx[col + 1];
for (int i = idx0; i < idx1; i++) {
if (A.nz_rows[i] == col) {
output += A.nz_values[i];
break;
}
}
}
return output;
}
/**
* This applies a given unary function on every value stored in the matrix
*
* B = f(A). A and B can be the same instance.
*
* @param input (Input) input matrix. Not modified
* @param func Unary function accepting a float
* @param output (Output) Matrix. Modified.
* @return The output matrix
*/
public static FMatrixSparseCSC apply( FMatrixSparseCSC input, FOperatorUnary func, @Nullable FMatrixSparseCSC output ) {
if (output == null) {
output = input.createLike();
}
if (input != output) {
output.copyStructure(input);
}
for (int i = 0; i < input.nz_length; i++) {
output.nz_values[i] = func.apply(input.nz_values[i]);
}
return output;
}
public static FMatrixSparseCSC apply( FMatrixSparseCSC input, FOperatorUnary func ) {
return apply(input, func, input);
}
/**
* This applies a given unary function on every nz row,value stored in the matrix
*
* B = f(A). A and B can be the same instance.
*
* @param input (Input) input matrix. Not modified
* @param func Binary function accepting (row-index, value)
* @param output (Output) Matrix. Modified.
* @return The output matrix
*/
public static FMatrixSparseCSC applyRowIdx( FMatrixSparseCSC input, FOperatorBinaryIdx func, @Nullable FMatrixSparseCSC output ) {
if (output == null) {
output = input.createLike();
}
if (input != output) {
output.copyStructure(input);
}
for (int i = 0; i < input.nz_length; i++) {
output.nz_values[i] = func.apply(input.nz_rows[i], input.nz_values[i]);
}
return output;
}
/**
* This applies a given unary function on every non-zero column, value stored in the matrix
*
* B = f(A). A and B can be the same instance.
*
* @param input (Input) input matrix. Not modified
* @param func Binary function accepting (column-index, value)
* @param output (Output) Matrix. Modified.
* @return The output matrix
*/
public static FMatrixSparseCSC applyColumnIdx( FMatrixSparseCSC input, FOperatorBinaryIdx func, @Nullable FMatrixSparseCSC output ) {
if (output == null) {
output = input.createLike();
}
if (input != output) {
output.copyStructure(input);
}
for (int col = 0; col < input.numCols; col++) {
for (int i = input.col_idx[col]; i < input.col_idx[col + 1]; i++) {
output.nz_values[i] = func.apply(col, input.nz_values[i]);
}
}
return output;
}
/**
* This accumulates the matrix values to a scalar value.
*
*
* result = initialValue
* for all (i,j) result = func( result, A[i,j] )
*
*
* @param input (Input) input matrix. Not modified
* @param initValue initial value for accumulator
* @param func Accumulator function defining "+" for accumulator += cellValue
* @return accumulated value
*/
public static float reduceScalar( FMatrixSparseCSC input, float initValue, FOperatorBinary func ) {
float result = initValue;
for (int i = 0; i < input.nz_length; i++) {
result = func.apply(result, input.nz_values[i]);
}
return result;
}
public static float reduceScalar( FMatrixSparseCSC input, FOperatorBinary func ) {
return reduceScalar(input, 0, func);
}
/**
* This accumulates the values per column to a scalar value
*
*
* for each column `j`,
* result = initialValue
* for each row 'i' result = func( result, A[i,j] )
* output[j] = result
*
*
* @param input (Input) input matrix. Not modified
* @param initValue initial value for accumulator
* @param func Accumulator function defining "+" for accumulator += cellValue
* @param output output (Output) Vector, where result can be stored in
* @param mask (Optional) Mask for specifying which entries should be overwritten
* @return a column-vector, where v[i] == values of column i reduced to scalar based on `func`
*/
public static FMatrixRMaj reduceColumnWise( FMatrixSparseCSC input, float initValue, FOperatorBinary func,
@Nullable FMatrixRMaj output, @Nullable Mask mask ) {
if (output == null) {
output = new FMatrixRMaj(1, input.numCols);
} else {
output.reshape(1, input.numCols);
}
for (int col = 0; col < input.numCols; col++) {
int start = input.col_idx[col];
int end = input.col_idx[col + 1];
float acc = initValue;
if (mask == null || mask.isSet(col)) {
for (int i = start; i < end; i++) {
acc = func.apply(acc, input.nz_values[i]);
}
}
output.data[col] = acc;
}
return output;
}
/**
* This accumulates the values per row to a scalar value
*
*
* for each column `j`,
* result = initialValue
* for each row 'i' result = func( result, A[i,j] )
* output[j] = result
*
*
* @param input (Input) input matrix. Not modified
* @param initValue initial value for accumulator
* @param func Accumulator function defining "+" for accumulator += cellValue
* @param output output (Output) Vector, where result can be stored in
* @return a row-vector, where v[i] == values of row i reduced to scalar based on `func`
*/
public static FMatrixRMaj reduceRowWise( FMatrixSparseCSC input, float initValue, FOperatorBinary func,
@Nullable FMatrixRMaj output ) {
if (output == null) {
output = new FMatrixRMaj(1, input.numRows);
} else {
output.reshape(1, input.numCols);
}
Arrays.fill(output.data, initValue);
for (int col = 0; col < input.numCols; col++) {
int start = input.col_idx[col];
int end = input.col_idx[col + 1];
for (int i = start; i < end; i++) {
output.data[input.nz_rows[i]] = func.apply(output.data[input.nz_rows[i]], input.nz_values[i]);
}
}
return output;
}
}