diff --git a/build.gradle b/build.gradle index 97dabe968..d0bfb7bd5 100644 --- a/build.gradle +++ b/build.gradle @@ -83,6 +83,7 @@ subprojects { compile 'com.google.code.findbugs:jsr305:3.0.2' // @Nullable testImplementation( 'org.junit.jupiter:junit-jupiter-api:5.4.0' ) + testImplementation( 'org.junit.jupiter:junit-jupiter-params:5.4.0' ) testRuntimeOnly( 'org.junit.jupiter:junit-jupiter-engine:5.4.0' ) ['core','generator-annprocess'].each { String a-> diff --git a/main/autocode/src/org/ejml/GenerateJavaCode32.java b/main/autocode/src/org/ejml/GenerateJavaCode32.java index 41c29cc40..9b50f306f 100644 --- a/main/autocode/src/org/ejml/GenerateJavaCode32.java +++ b/main/autocode/src/org/ejml/GenerateJavaCode32.java @@ -54,6 +54,20 @@ public GenerateJavaCode32() { prefix32.add("FUnary"); prefix64.add("DBinary"); prefix32.add("FBinary"); + prefix64.add("DMonoid"); + prefix32.add("FMonoid"); + prefix64.add("DSemiRing"); + prefix32.add("FSemiRing"); + prefix64.add("DMonoids"); + prefix32.add("FMonoids"); + prefix64.add("DSemiRings"); + prefix32.add("FSemiRings"); + prefix64.add("PrimitiveDMask"); + prefix32.add("PrimitiveFMask"); + prefix64.add("SparseDMask"); + prefix32.add("SparseFMask"); + prefix64.add("DMasks"); + prefix32.add("FMasks"); prefix64.add("DScalar"); prefix32.add("FScalar"); prefix64.add("DMatrix"); @@ -93,6 +107,11 @@ public GenerateJavaCode32() { converter.replacePattern("DScalar", "FScalar"); converter.replacePattern("DUnary", "FUnary"); converter.replacePattern("DBinary", "FBinary"); + converter.replacePattern("DMonoid", "FMonoid"); + converter.replacePattern("DSemiRing", "FSemiRing"); + converter.replacePattern("SparseDMask", "SparseFMask"); + converter.replacePattern("PrimitiveDMask", "PrimitiveFMask"); + converter.replacePattern("DMasks", "FMasks"); converter.replacePattern("ConvertD", "ConvertF"); converter.replacePattern("DGrowArray", "FGrowArray"); converter.replacePattern("DMatrix", "FMatrix"); @@ -133,6 +152,8 @@ public static void main(String[] args ) { "main/ejml-core/test/org/ejml/data", "main/ejml-core/src/org/ejml/ops", "main/ejml-core/test/org/ejml/ops", + "main/ejml-core/src/org/ejml/masks", + "main/ejml-core/test/org/ejml/masks", "main/ejml-experimental/src/org/ejml/dense/row/decomposition/bidiagonal/" }; diff --git a/main/ejml-core/src/org/ejml/data/DMatrixRMaj.java b/main/ejml-core/src/org/ejml/data/DMatrixRMaj.java index 4ee6c9d15..c0891fc00 100644 --- a/main/ejml-core/src/org/ejml/data/DMatrixRMaj.java +++ b/main/ejml-core/src/org/ejml/data/DMatrixRMaj.java @@ -333,6 +333,13 @@ public void zero() { Arrays.fill(data, 0, getNumElements(), 0.0); } + /** + * Sets all elements equal to the specified value. + */ + public void fill(double value) { + Arrays.fill(data, 0, getNumElements(), value); + } + /** * Creates and returns a matrix which is idential to this one. * diff --git a/main/ejml-core/src/org/ejml/masks/DMasks.java b/main/ejml-core/src/org/ejml/masks/DMasks.java new file mode 100644 index 000000000..329dc13b1 --- /dev/null +++ b/main/ejml-core/src/org/ejml/masks/DMasks.java @@ -0,0 +1,43 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +import org.ejml.data.*; + +/** + * Helper class to get the corresponding mask builder based on a matrix or primitive array + */ +public class DMasks { + public static PrimitiveDMask.Builder builder(double[] values) { + return new PrimitiveDMask.Builder(values); + } + + public static PrimitiveDMask.Builder builder(DMatrixD1 matrix) { + return new PrimitiveDMask.Builder(matrix.data).withNumCols(matrix.numCols); + } + + public static MaskBuilder builder(DMatrixSparseCSC matrix, boolean structural){ + if (structural) { + return new SparseStructuralMask.Builder(matrix); + } + else { + return new SparseDMask.Builder(matrix); + } + } +} diff --git a/main/ejml-core/src/org/ejml/masks/Mask.java b/main/ejml-core/src/org/ejml/masks/Mask.java new file mode 100644 index 000000000..2d70fefbe --- /dev/null +++ b/main/ejml-core/src/org/ejml/masks/Mask.java @@ -0,0 +1,75 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + + +import org.ejml.MatrixDimensionException; +import org.ejml.data.Matrix; + +/** + * Mask used for specifying which matrix entries should be computed + */ +public abstract class Mask { + // useful for sparse matrices, as actual negation would be costly and result in dense masks + public final boolean negated; + + // whether the result should replace the initial result (otherwise merge them) + // same as the "GrB_REPLACE" descriptor + public final boolean replace; + + public Mask(boolean negated, boolean replace) { + this.negated = negated; + this.replace = replace; + } + + public abstract boolean isSet(int row, int col); + + public abstract int getNumCols(); + + public abstract int getNumRows(); + + public void print() { + String result = ""; + for (int row = 0; row < getNumRows(); row++) { + for (int col = 0; col < getNumCols(); col++) { + result += isSet(row, col) ? "+ " : "- "; + } + result += System.lineSeparator(); + } + + System.out.println(result); + }; + + /** + * Checks whether the dimensions of the mask and matrix match + * @param matrix the mask is applied to + */ + public void compatible(Matrix matrix) { + if (matrix.getNumCols() != getNumCols() || matrix.getNumRows() != getNumRows()) { + throw new MatrixDimensionException(String.format( + "Mask of (%d, %d) cannot be applied for matrix (%d, %d)", + getNumRows(), getNumCols(), matrix.getNumCols(), matrix.getNumCols() + )); + } + } + + // TODO: use an Iterator as it should be faster as stepping can be used -> no need to call for each entry? + // Problem .. dense matrices are row-based, whereas existing sparse matrix format is column based + //public abstract Iterator +} diff --git a/main/ejml-core/src/org/ejml/masks/MaskBuilder.java b/main/ejml-core/src/org/ejml/masks/MaskBuilder.java new file mode 100644 index 000000000..535cf4ce6 --- /dev/null +++ b/main/ejml-core/src/org/ejml/masks/MaskBuilder.java @@ -0,0 +1,36 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +public abstract class MaskBuilder { + protected boolean negated = false; + protected boolean replace = true; + + public MaskBuilder withNegated(boolean negated) { + this.negated = negated; + return this; + } + + public MaskBuilder withReplace(boolean replace) { + this.replace = replace; + return this; + } + + public abstract MASK build(); +} diff --git a/main/ejml-core/src/org/ejml/masks/PrimitiveDMask.java b/main/ejml-core/src/org/ejml/masks/PrimitiveDMask.java new file mode 100644 index 000000000..57fa5a232 --- /dev/null +++ b/main/ejml-core/src/org/ejml/masks/PrimitiveDMask.java @@ -0,0 +1,93 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +import org.ejml.MatrixDimensionException; + +public class PrimitiveDMask extends Mask { + // TODO: Check if creating dense boolean mask is worth it (currently converting to boolean on the fly) + private final double[] values; + // alternatively declare extra separate masks for primitive and dense DMatrix + private final int numCols; + // + private final double zeroElement; + + public PrimitiveDMask(double[] values, int numCols, boolean negated, boolean replace, double zeroElement) { + // for dense structures they cannot be used for structural masks + super(negated, replace); + this.values = values; + this.numCols = numCols; + this.zeroElement = zeroElement; + } + + @Override + public boolean isSet(int row, int col) { + // XOR as negated flips the mask flag + return negated ^ (values[row * numCols + col] != zeroElement); + } + + @Override + public int getNumCols() { + return numCols; + } + + @Override + public int getNumRows() { + return values.length / numCols; + } + + public boolean isSet(int index) { + return negated ^ (values[index] != zeroElement); + } + + public void compatible(double[] vector) { + if (vector.length != values.length) { + throw new MatrixDimensionException(String.format( + "Mask of length %d cannot be applied to vector of lenght %d", + values.length, vector.length + )); + } + } + + public static class Builder extends MaskBuilder { + private double[] values; + private int numCols = 1; + private double zeroElement = 0; + + public Builder(double[] values) { + this.values = values; + } + + public Builder withNumCols(int numCols) { + this.numCols = numCols; + return this; + } + + public Builder withZeroElement(double zeroElement) { + this.zeroElement = zeroElement; + return this; + } + + @Override + public PrimitiveDMask build() { + return new PrimitiveDMask(values, numCols, negated, replace, zeroElement); + } + } +} + diff --git a/main/ejml-core/src/org/ejml/masks/SparseDMask.java b/main/ejml-core/src/org/ejml/masks/SparseDMask.java new file mode 100644 index 000000000..3dfb06c3d --- /dev/null +++ b/main/ejml-core/src/org/ejml/masks/SparseDMask.java @@ -0,0 +1,66 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +import org.ejml.data.DMatrixSparseCSC; + +public class SparseDMask extends Mask { + protected final DMatrixSparseCSC matrix; + protected final double zeroElement; + + public SparseDMask(DMatrixSparseCSC matrix, boolean negated, boolean replace, double zeroElement) { + super(negated, replace); + this.matrix = matrix; + this.zeroElement = zeroElement; + } + + @Override + public boolean isSet(int row, int col) { + return negated ^ (matrix.unsafe_get(row, col) != zeroElement); + } + + @Override + public int getNumCols() { + return matrix.numCols; + } + + @Override + public int getNumRows() { + return matrix.numRows; + } + + public static class Builder extends MaskBuilder { + private DMatrixSparseCSC matrix; + private double zeroElement = 0; + + public Builder(DMatrixSparseCSC matrix) { + this.matrix = matrix; + } + + public Builder withZeroElement(double zeroElement) { + this.zeroElement = zeroElement; + return this; + } + + @Override + public SparseDMask build() { + return new SparseDMask(matrix, negated, replace, zeroElement); + } + } +} diff --git a/main/ejml-core/src/org/ejml/masks/SparseStructuralMask.java b/main/ejml-core/src/org/ejml/masks/SparseStructuralMask.java new file mode 100644 index 000000000..6c08ab274 --- /dev/null +++ b/main/ejml-core/src/org/ejml/masks/SparseStructuralMask.java @@ -0,0 +1,62 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +import org.ejml.data.MatrixSparse; + +/** + * only looking if the entry is assigned in the source(disregarding the actual stored value) + * ! it does not copy the input matrix -> changing the matrix structure will also affect the mask + */ +public class SparseStructuralMask extends Mask { + private final MatrixSparse matrix; + + public SparseStructuralMask(MatrixSparse matrix, boolean negated, boolean replace) { + super(negated, replace); + this.matrix = matrix; + } + + @Override + public boolean isSet(int row, int col) { + return negated ^ matrix.isAssigned(row, col); + } + + @Override + public int getNumCols() { + return matrix.getNumCols(); + } + + @Override + public int getNumRows() { + return matrix.getNumRows(); + } + + public static class Builder extends MaskBuilder { + private MatrixSparse matrix; + + public Builder(MatrixSparse matrix) { + this.matrix = matrix; + } + + @Override + public SparseStructuralMask build() { + return new SparseStructuralMask(matrix, negated, replace); + } + } +} diff --git a/main/ejml-core/src/org/ejml/ops/DMonoid.java b/main/ejml-core/src/org/ejml/ops/DMonoid.java new file mode 100644 index 000000000..417111626 --- /dev/null +++ b/main/ejml-core/src/org/ejml/ops/DMonoid.java @@ -0,0 +1,34 @@ +/* + * Copyright (c) 2009-2020, 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.ops; + +public class DMonoid { + // possible re-interpreting 0 + public final double id; + public final DBinaryOperator func; + + DMonoid(double id, DBinaryOperator func) { + this.id = id; + this.func = func; + } + + DMonoid(DBinaryOperator func) { + this(0, func); + } +} \ No newline at end of file diff --git a/main/ejml-core/src/org/ejml/ops/DMonoids.java b/main/ejml-core/src/org/ejml/ops/DMonoids.java new file mode 100644 index 000000000..16c841136 --- /dev/null +++ b/main/ejml-core/src/org/ejml/ops/DMonoids.java @@ -0,0 +1,36 @@ +/* + * Copyright (c) 2009-2020, 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.ops; + +/** + * as defined in the graphblas c-api (https://people.eecs.berkeley.edu/~aydin/GraphBLAS_API_C_v13.pdf) + * p. 26 + */ +public final class DMonoids { + public static final DMonoid AND = new DMonoid(1, (a, b) -> (a == 0 || b == 0) ? 0 : 1); + public static final DMonoid OR = new DMonoid(0, (a, b) -> (a != 0 || b != 0) ? 1 : 0); + public static final DMonoid XOR = new DMonoid(0, (a, b) -> ((a == 0 && b == 0) || (a != 0 && b != 0)) ? 0 : 1); + public static final DMonoid XNOR = new DMonoid(0, (a, b) -> ((a == 0 && b == 0) || (a != 0 && b != 0)) ? 1 : 0); + + public static final DMonoid PLUS = new DMonoid(0, Double::sum); + public static final DMonoid TIMES = new DMonoid(1, (a, b) -> a * b); + + public final static DMonoid MIN = new DMonoid(Double.MAX_VALUE, (a, b) -> (a <= b) ? a : b); + public final static DMonoid MAX = new DMonoid(Double.MIN_VALUE, (a, b) -> (a >= b) ? a : b); +} \ No newline at end of file diff --git a/main/ejml-core/src/org/ejml/ops/DSemiRing.java b/main/ejml-core/src/org/ejml/ops/DSemiRing.java new file mode 100644 index 000000000..2720b7fd7 --- /dev/null +++ b/main/ejml-core/src/org/ejml/ops/DSemiRing.java @@ -0,0 +1,29 @@ +/* + * Copyright (c) 2009-2020, 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.ops; + +public class DSemiRing { + public final DMonoid add; + public final DMonoid mult; + + DSemiRing(DMonoid add, DMonoid mult) { + this.add = add; + this.mult = mult; + } +} \ No newline at end of file diff --git a/main/ejml-core/src/org/ejml/ops/DSemiRings.java b/main/ejml-core/src/org/ejml/ops/DSemiRings.java new file mode 100644 index 000000000..d5d29ff75 --- /dev/null +++ b/main/ejml-core/src/org/ejml/ops/DSemiRings.java @@ -0,0 +1,52 @@ +/* + * Copyright (c) 2009-2020, 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.ops; + +import static org.ejml.ops.DMonoids.*; + +/** + * as defined in the graphblas c-api (https://people.eecs.berkeley.edu/~aydin/GraphBLAS_API_C_v13.pdf) + * p. 27-28 + * Note: some don’t have a multiplicative annihilator (thus no "true" semi-rings") + */ +public final class DSemiRings { + public static final DSemiRing PLUS_TIMES = new DSemiRing(PLUS, TIMES); + public static final DSemiRing MIN_PLUS = new DSemiRing(MIN, PLUS); + public static final DSemiRing MAX_PLUS = new DSemiRing(MAX, PLUS); + public static final DSemiRing MIN_TIMES = new DSemiRing(MIN, TIMES); + public static final DSemiRing MIN_MAX = new DSemiRing(MIN, MAX); + public static final DSemiRing MAX_MIN = new DSemiRing(MAX, MIN); + public static final DSemiRing MAX_TIMES = new DSemiRing(MAX, TIMES); + public static final DSemiRing PLUS_MIN = new DSemiRing(PLUS, MIN); + + public static final DSemiRing OR_AND = new DSemiRing(OR, AND); + public static final DSemiRing AND_OR = new DSemiRing(AND, OR); + public static final DSemiRing XOR_AND = new DSemiRing(XOR, AND); + public static final DSemiRing XNOR_OR = new DSemiRing(XNOR, OR); + + // only private as they have no identity element, hence can only be used for add + private static final DMonoid FIRST = new DMonoid(Double.NaN, (x, y) -> x); + private static final DMonoid SECOND = new DMonoid(Double.NaN, (x, y) -> y); + + // semi-rings with no multiplicative annihilator + public static final DSemiRing MIN_FIRST = new DSemiRing(MIN, FIRST); + public static final DSemiRing MIN_SECOND = new DSemiRing(MIN, SECOND); + public static final DSemiRing MAX_FIRST = new DSemiRing(MAX, FIRST); + public static final DSemiRing MAX_SECOND = new DSemiRing(MAX, SECOND); +} \ No newline at end of file diff --git a/main/ejml-core/test/org/ejml/masks/TestPrimitiveDMasks.java b/main/ejml-core/test/org/ejml/masks/TestPrimitiveDMasks.java new file mode 100644 index 000000000..5307e6a3b --- /dev/null +++ b/main/ejml-core/test/org/ejml/masks/TestPrimitiveDMasks.java @@ -0,0 +1,65 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +import org.ejml.data.DMatrixRMaj; +import org.ejml.dense.row.RandomMatrices_DDRM; +import org.junit.jupiter.api.Test; + +import java.util.Random; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +public class TestPrimitiveDMasks { + + @Test + void testBasedOnPrimitiveArrays() { + double[] values = {2, 0, 4, 0, 0, -1}; + + PrimitiveDMask.Builder maskBuilder = new PrimitiveDMask.Builder(values); + PrimitiveDMask mask = maskBuilder.withNegated(false).build(); + PrimitiveDMask negated_mask = maskBuilder.withNegated(true).build(); + + boolean[] expected = {true, false, true, false, false, true}; + + for (int i = 0; i < values.length; i++) { + assertEquals(mask.isSet(i), expected[i]); + assertEquals(negated_mask.isSet(i), !expected[i]); + } + } + + @Test + void testBasedOnDenseMatrix() { + int dim = 20; + DMatrixRMaj matrix = RandomMatrices_DDRM.rectangle(dim, dim, new Random(42)); + + PrimitiveDMask.Builder maskBuilder = DMasks.builder(matrix); + PrimitiveDMask mask = maskBuilder.withNegated(false).build(); + PrimitiveDMask negated_mask = maskBuilder.withNegated(true).build(); + + + for (int row = 0; row < dim; row++) { + for (int col = 0; col < dim; col++) { + boolean expected = (matrix.get(row, col) != 0); + assertEquals(mask.isSet(row, col), expected); + assertEquals(negated_mask.isSet(row, col), !expected); + } + } + } +} diff --git a/main/ejml-core/test/org/ejml/masks/TestSparseDMasks.java b/main/ejml-core/test/org/ejml/masks/TestSparseDMasks.java new file mode 100644 index 000000000..f0494a83a --- /dev/null +++ b/main/ejml-core/test/org/ejml/masks/TestSparseDMasks.java @@ -0,0 +1,50 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.sparse.csc.RandomMatrices_DSCC; +import org.junit.jupiter.api.Test; + +import java.util.Random; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +public class TestSparseDMasks { + + @Test + void testSparseMasks() { + int dim = 20; + DMatrixSparseCSC matrix = RandomMatrices_DSCC.rectangle(dim, dim, 5, new Random(42)); + + SparseStructuralMask.Builder builder = new SparseStructuralMask.Builder(matrix); + Mask mask = builder.withNegated(false).build(); + Mask negated_mask = builder.withNegated(true).build(); + + var it = matrix.createCoordinateIterator(); + + for (int row = 0; row < dim; row++) { + for (int col = 0; col < dim; col++) { + boolean expected = (matrix.get(row, col) != 0); + assertEquals(mask.isSet(row, col), expected); + assertEquals(negated_mask.isSet(row, col), !expected); + } + } + } +} diff --git a/main/ejml-core/test/org/ejml/masks/TestSparseStructuralMasks.java b/main/ejml-core/test/org/ejml/masks/TestSparseStructuralMasks.java new file mode 100644 index 000000000..b47ebb324 --- /dev/null +++ b/main/ejml-core/test/org/ejml/masks/TestSparseStructuralMasks.java @@ -0,0 +1,50 @@ +/* + * Copyright (c) 2009-2020, 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.masks; + +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.sparse.csc.RandomMatrices_DSCC; +import org.junit.jupiter.api.Test; + +import java.util.Random; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +public class TestSparseStructuralMasks { + + @Test + void testSparseStructuralMasks() { + int dim = 20; + DMatrixSparseCSC matrix = RandomMatrices_DSCC.rectangle(dim, dim, 5, new Random(42)); + + SparseStructuralMask.Builder builder = new SparseStructuralMask.Builder(matrix); + Mask mask = builder.withNegated(false).build(); + Mask negated_mask = builder.withNegated(true).build(); + + var it = matrix.createCoordinateIterator(); + + for (int row = 0; row < dim; row++) { + for (int col = 0; col < dim; col++) { + boolean expected = matrix.isAssigned(row, col); + assertEquals(mask.isSet(row, col), expected); + assertEquals(negated_mask.isSet(row, col), !expected); + } + } + } +} diff --git a/main/ejml-core/test/org/ejml/ops/TestAlgebraicStructures.java b/main/ejml-core/test/org/ejml/ops/TestAlgebraicStructures.java new file mode 100644 index 000000000..6f8957b8e --- /dev/null +++ b/main/ejml-core/test/org/ejml/ops/TestAlgebraicStructures.java @@ -0,0 +1,57 @@ +/* + * Copyright (c) 2009-2020, 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.ops; + + +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +public class TestAlgebraicStructures { + + @Test + public void testPreDefinedMonoids() { + assertEquals(0, DMonoids.AND.func.apply(20, 0)); + assertEquals(0, DMonoids.AND.func.apply(0, 20)); + assertEquals(0, DMonoids.AND.func.apply(0, 0)); + assertEquals(1, DMonoids.AND.func.apply(43, 20)); + + assertEquals(1, DMonoids.OR.func.apply(43, 20)); + assertEquals(1, DMonoids.OR.func.apply(0, 20)); + assertEquals(1, DMonoids.OR.func.apply(43, 0)); + assertEquals(0, DMonoids.OR.func.apply(0, 0)); + + assertEquals(0, DMonoids.XOR.func.apply(43, 20)); + assertEquals(1, DMonoids.XOR.func.apply(0, 20)); + assertEquals(1, DMonoids.XOR.func.apply(43, 0)); + assertEquals(0, DMonoids.XOR.func.apply(0, 0)); + + assertEquals(1, DMonoids.XNOR.func.apply(43, 20)); + assertEquals(0, DMonoids.XNOR.func.apply(0, 20)); + assertEquals(0, DMonoids.XNOR.func.apply(43, 0)); + assertEquals(1, DMonoids.XNOR.func.apply(0, 0)); + + assertEquals(43, DMonoids.MAX.func.apply(43, 20)); + assertEquals(20, DMonoids.MIN.func.apply(43, 20)); + + + assertEquals(63, DMonoids.PLUS.func.apply(43, 20)); + assertEquals(860, DMonoids.TIMES.func.apply(43, 20)); + } +} \ No newline at end of file diff --git a/main/ejml-dsparse/src/org/ejml/sparse/csc/CommonOpsWithSemiRing_DSCC.java b/main/ejml-dsparse/src/org/ejml/sparse/csc/CommonOpsWithSemiRing_DSCC.java new file mode 100644 index 000000000..290b1e286 --- /dev/null +++ b/main/ejml-dsparse/src/org/ejml/sparse/csc/CommonOpsWithSemiRing_DSCC.java @@ -0,0 +1,296 @@ +/* + * Copyright (c) 2009-2020, 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 org.ejml.MatrixDimensionException; +import org.ejml.data.DGrowArray; +import org.ejml.data.DMatrixRMaj; +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.data.IGrowArray; +import org.ejml.masks.Mask; +import org.ejml.ops.DBinaryOperator; +import org.ejml.ops.DSemiRing; +import org.ejml.sparse.csc.misc.ImplCommonOpsWithSemiRing_DSCC; +import org.ejml.sparse.csc.mult.ImplSparseSparseMultWithSemiRing_DSCC; + +import javax.annotation.Nullable; + +import static org.ejml.UtilEjml.reshapeOrDeclare; +import static org.ejml.UtilEjml.stringShapes; +import static org.ejml.sparse.csc.MaskUtil_DSCC.combineOutputs; + + +public class CommonOpsWithSemiRing_DSCC { + + public static DMatrixSparseCSC mult(DMatrixSparseCSC A, DMatrixSparseCSC B, @Nullable DMatrixSparseCSC output, DSemiRing semiRing, + @Nullable Mask mask, @Nullable DBinaryOperator accumulator) { + return mult(A, B, output, semiRing, mask, accumulator, null, null); + } + + /** + * Performs matrix multiplication. output = A*B + * + * @param A (Input) Matrix. Not modified. + * @param B (Input) Matrix. Not modified. + * @param output (Output) Storage for results. Data length is increased if increased if insufficient. + * @param mask Mask for specifying which entries should be overwritten + * @param accumulator Operator to combine result with existing entries in output matrix + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static DMatrixSparseCSC mult(DMatrixSparseCSC A, DMatrixSparseCSC B, @Nullable DMatrixSparseCSC output, DSemiRing semiRing, + @Nullable Mask mask, @Nullable DBinaryOperator accumulator, @Nullable IGrowArray gw, + @Nullable DGrowArray gx) { + if (A.numCols != B.numRows) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + + // !! important to do before reshape + DMatrixSparseCSC initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accumulator, output) ? output.copy() : null; + + output = reshapeOrDeclare(output, A, A.numRows, B.numCols); + if (mask != null) { + mask.compatible(output); + } + + ImplSparseSparseMultWithSemiRing_DSCC.mult(A, B, output, semiRing, mask, gw, gx); + + return combineOutputs(output, accumulator, initialOutput); + } + + public static DMatrixSparseCSC multTransA(DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC output, DSemiRing semiRing, + @Nullable Mask mask, @Nullable DBinaryOperator accumulator, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + if (A.numRows != B.numRows) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + // !! important to do before reshape + DMatrixSparseCSC initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accumulator, output) ? output.copy() : null; + output = reshapeOrDeclare(output, A, A.numCols, B.numCols); + if (mask != null) { + mask.compatible(output); + } + + ImplSparseSparseMultWithSemiRing_DSCC.multTransA(A, B, output, semiRing, mask, gw, gx); + + return combineOutputs(output, accumulator, initialOutput); + } + + /** + * Performs matrix multiplication. output = A*BT. B needs to be sorted and will be sorted if it + * has not already been sorted. + * + * @param A (Input) Matrix. Not modified. + * @param B (Input) Matrix. Value not modified but indicies will be sorted if not sorted already. + * @param output (Output) Storage for results. Data length is increased if increased if insufficient. + * @param mask Mask for specifying which entries should be overwritten + * @param accumulator Operator to combine result with existing entries in output matrix + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static DMatrixSparseCSC multTransB(DMatrixSparseCSC A, DMatrixSparseCSC B, @Nullable DMatrixSparseCSC output, DSemiRing semiRing, + @Nullable Mask mask, @Nullable DBinaryOperator accumulator, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + if (A.numCols != B.numCols) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + DMatrixSparseCSC initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accumulator, output) ? output.copy() : null; + output = reshapeOrDeclare(output, A, A.numRows, B.numRows); + if (mask != null) { + mask.compatible(output); + } + if (!B.isIndicesSorted()) B.sortIndices(null); + ImplSparseSparseMultWithSemiRing_DSCC.multTransB(A, B, output, semiRing, mask, gw, gx); + return combineOutputs(output, accumulator, initialOutput); + } + + + // sparse-dense variations (does not support masked yet. look at primitive versions for vector variants) + + /** + * Performs matrix multiplication. output = A*B + * + * @param A Matrix + * @param B Dense Matrix + * @param output Dense Matrix + */ + public static DMatrixRMaj mult(DMatrixSparseCSC A, DMatrixRMaj B, @Nullable DMatrixRMaj output, DSemiRing semiRing) { + if (A.numCols != B.numRows) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + + output = reshapeOrDeclare(output, A.numRows, B.numCols); + + ImplSparseSparseMultWithSemiRing_DSCC.mult(A, B, output, semiRing); + + return output; + } + + /** + *

output = output + A*B

+ */ + public static void multAdd(DMatrixSparseCSC A, DMatrixRMaj B, @Nullable DMatrixRMaj output, DSemiRing semiRing) { + if (A.numRows != output.numRows || B.numCols != output.numCols) + throw new IllegalArgumentException("Inconsistent matrix shapes. " + stringShapes(A, B, output)); + + ImplSparseSparseMultWithSemiRing_DSCC.multAdd(A, B, output, semiRing); + } + + /** + * Performs matrix multiplication. output = AT*B + * + * @param A Matrix + * @param B Dense Matrix + * @param output Dense Matrix + */ + public static DMatrixRMaj multTransA(DMatrixSparseCSC A, DMatrixRMaj B, @Nullable DMatrixRMaj output, DSemiRing semiRing) { + if (A.numRows != B.numRows) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + + output = reshapeOrDeclare(output, A.numCols, B.numCols); + + ImplSparseSparseMultWithSemiRing_DSCC.multTransA(A, B, output, semiRing); + + return output; + } + + /** + *

output = output + AT*B

+ */ + public static void multAddTransA(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj output, DSemiRing semiRing) { + if (A.numCols != output.numRows || B.numCols != output.numCols) + throw new IllegalArgumentException("Inconsistent matrix shapes. " + stringShapes(A, B, output)); + + ImplSparseSparseMultWithSemiRing_DSCC.multAddTransA(A, B, output, semiRing); + } + + /** + * Performs matrix multiplication. output = A*BT + * + * @param A Matrix + * @param B Dense Matrix + * @param output Dense Matrix + */ + public static DMatrixRMaj multTransB(DMatrixSparseCSC A, DMatrixRMaj B, @Nullable DMatrixRMaj output, DSemiRing semiRing) { + // todo: combine with multAdd as only difference is that output is filled with zero before + if (A.numCols != B.numCols) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + output = reshapeOrDeclare(output, A.numRows, B.numRows); + + ImplSparseSparseMultWithSemiRing_DSCC.multTransB(A, B, output, semiRing); + + return output; + } + + /** + *

output = output + A*BT

+ */ + public static void multAddTransB(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj output, DSemiRing semiRing) { + if (A.numRows != output.numRows || B.numRows != output.numCols) + throw new IllegalArgumentException("Inconsistent matrix shapes. " + stringShapes(A, B, output)); + + // TODO: ? this is basically the equivalent of graphblas mult with specified accumulator op + ImplSparseSparseMultWithSemiRing_DSCC.multAddTransB(A, B, output, semiRing); + } + + /** + * Performs matrix multiplication. output = AT*BT + * + * @param A Matrix + * @param B Dense Matrix + * @param output Dense Matrix + */ + public static DMatrixRMaj multTransAB(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj output, DSemiRing semiRing) { + if (A.numRows != B.numCols) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + output = reshapeOrDeclare(output, A.numCols, B.numRows); + + ImplSparseSparseMultWithSemiRing_DSCC.multTransAB(A, B, output, semiRing); + + return output; + } + + + /** + *

C = C + AT*BT

+ */ + public static void multAddTransAB(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + // TODO: this is basically the equivalent of graphblas mult with specified accumulator op + if (A.numCols != C.numRows || B.numRows != C.numCols) + throw new IllegalArgumentException("Inconsistent matrix shapes. " + stringShapes(A, B, C)); + + ImplSparseSparseMultWithSemiRing_DSCC.multAddTransAB(A, B, C, semiRing); + } + + // sparse versions again + + /** + * Performs matrix addition:
+ * output = αA + βB + * + * @param alpha scalar value multiplied against A + * @param A Matrix + * @param beta scalar value multiplied against B + * @param B Matrix + * @param output (Optional) Output matrix. + * @param mask Mask for specifying which entries should be overwritten + * @param accumulator Operator to combine result with existing entries in output matrix + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static DMatrixSparseCSC add(double alpha, DMatrixSparseCSC A, double beta, DMatrixSparseCSC B, @Nullable DMatrixSparseCSC output, DSemiRing semiRing, + @Nullable Mask mask, @Nullable DBinaryOperator accumulator, @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + if (A.numRows != B.numRows || A.numCols != B.numCols) + throw new MatrixDimensionException("Inconsistent matrix shapes. " + stringShapes(A, B)); + DMatrixSparseCSC initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accumulator, output) ? output.copy() : null; + output = reshapeOrDeclare(output, A, A.numRows, A.numCols); + if (mask != null) { + mask.compatible(output); + } + + ImplCommonOpsWithSemiRing_DSCC.add(alpha, A, beta, B, output, semiRing, mask, gw, gx); + + return combineOutputs(output, accumulator, initialOutput); + } + + /** + * 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 mask Mask for specifying which entries should be overwritten + * @param accumulator Operator to combine result with existing entries in output matrix + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static DMatrixSparseCSC elementMult(DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC output, DSemiRing semiRing, + @Nullable Mask mask, @Nullable DBinaryOperator accumulator, @Nullable IGrowArray gw, @Nullable DGrowArray 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)); + DMatrixSparseCSC initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accumulator, output) ? output.copy() : null; + output = reshapeOrDeclare(output, A, A.numRows, A.numCols); + if (mask != null) { + mask.compatible(output); + } + + ImplCommonOpsWithSemiRing_DSCC.elementMult(A, B, output, semiRing, mask, gw, gx); + + return combineOutputs(output, accumulator, initialOutput); + } +} + diff --git a/main/ejml-dsparse/src/org/ejml/sparse/csc/CommonOps_DSCC.java b/main/ejml-dsparse/src/org/ejml/sparse/csc/CommonOps_DSCC.java index 26cf44fc5..4844a3954 100644 --- a/main/ejml-dsparse/src/org/ejml/sparse/csc/CommonOps_DSCC.java +++ b/main/ejml-dsparse/src/org/ejml/sparse/csc/CommonOps_DSCC.java @@ -26,6 +26,7 @@ import org.ejml.dense.row.CommonOps_DDRM; import org.ejml.interfaces.decomposition.LUSparseDecomposition_F64; import org.ejml.interfaces.linsol.LinearSolverSparse; +import org.ejml.masks.Mask; import org.ejml.ops.DBinaryOperator; import org.ejml.ops.DUnaryOperator; import org.ejml.sparse.FillReducing; @@ -38,6 +39,7 @@ import java.util.Arrays; import static org.ejml.UtilEjml.*; +import static org.ejml.sparse.csc.MaskUtil_DSCC.combineOutputs; /** * @author Peter Abeles @@ -352,24 +354,26 @@ public static void symmLowerToFull( DMatrixSparseCSC A , DMatrixSparseCSC B , @N /** * Performs matrix addition:
- * C = αA + βB - * - * @param alpha scalar value multiplied against A + * output = αA + βB + * @param alpha scalar value multiplied against A * @param A Matrix * @param beta scalar value multiplied against B * @param B Matrix - * @param C Output matrix. + * @param output Output matrix. * @param gw (Optional) Storage for internal workspace. Can be null. * @param gx (Optional) Storage for internal workspace. Can be null. + * @return */ - public static void add(double alpha, DMatrixSparseCSC A, double beta, DMatrixSparseCSC B, DMatrixSparseCSC C, - @Nullable IGrowArray gw, @Nullable DGrowArray gx) + public static DMatrixSparseCSC add(double alpha, DMatrixSparseCSC A, double beta, DMatrixSparseCSC B, @Nullable DMatrixSparseCSC output, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { if( A.numRows != B.numRows || A.numCols != B.numCols ) throw new MatrixDimensionException("Inconsistent matrix shapes. "+stringShapes(A,B)); - C.reshape(A.numRows,A.numCols); + output = reshapeOrDeclare(output, A, A.numRows,A.numCols); + + ImplCommonOps_DSCC.add(alpha,A,beta,B,output, gw, gx); - ImplCommonOps_DSCC.add(alpha,A,beta,B,C, gw, gx); + return output; } public static DMatrixSparseCSC identity(int length ) { @@ -578,22 +582,23 @@ public static double elementSum( DMatrixSparseCSC A ) { /** * Performs an element-wise multiplication.
- * C[i,j] = A[i,j]*B[i,j]
+ * output[i,j] = A[i,j]*B[i,j]
* All matrices must have the same shape. - * - * @param A (Input) Matrix. + * @param A (Input) Matrix. * @param B (Input) Matrix - * @param C (Output) Matrix. data array is grown to min(A.nz_length,B.nz_length), resulting a in a large speed boost. + * @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 void elementMult( DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC C , - @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + public static DMatrixSparseCSC elementMult(DMatrixSparseCSC A, DMatrixSparseCSC B, @Nullable DMatrixSparseCSC output , + @Nullable IGrowArray gw, @Nullable DGrowArray 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)); - C.reshape(A.numRows,A.numCols); + output = reshapeOrDeclare(output, A, A.numRows,A.numCols); + + ImplCommonOps_DSCC.elementMult(A,B,output,gw,gx); - ImplCommonOps_DSCC.elementMult(A,B,C,gw,gx); + return output; } /** @@ -1889,6 +1894,30 @@ public static DMatrixSparseCSC apply(DMatrixSparseCSC input, DUnaryOperator func return apply(input, func, input); } + // !! FIXME: apply also calculates unneeded fields .. this can produce wrong results with another accumulator than "first" + public static DMatrixSparseCSC apply(DMatrixSparseCSC input, DUnaryOperator func, + @Nullable DMatrixSparseCSC output, @Nullable Mask mask, @Nullable DBinaryOperator accum) { + DMatrixSparseCSC initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accum, output) ? output.copy() : null; + // set correct structure + if (output == null) { + output = input.createLike(); + } else if (input != output) { + output.copyStructure(input); + } + if (mask != null) { + mask.compatible(output); + } + + + // dont check for the mask here as the computation is too simple to get any improvements when skipping + // would more likely prohibit vectorisation + for (int i = 0; i < input.nz_length; i++) { + output.nz_values[i] = func.apply(input.nz_values[i]); + } + + return combineOutputs(output, initialOutput, mask, accum); + } + /** * This accumulates the matrix values to a scalar value. * @@ -1932,11 +1961,12 @@ public static double reduceScalar(DMatrixSparseCSC input, DBinaryOperator func) * @param output output (Output) Vector, where result can be stored in * @return a column-vector, where v[i] == values of column i reduced to scalar based on `func` */ - public static DMatrixRMaj reduceColumnWise(DMatrixSparseCSC input, double initValue, DBinaryOperator func, @Nullable DMatrixRMaj output) { - if (output == null) { - output = new DMatrixRMaj(1, input.numCols); - } else { - output.reshape(1, input.numCols); + public static DMatrixRMaj reduceColumnWise(DMatrixSparseCSC input, double initValue, DBinaryOperator func, + @Nullable DMatrixRMaj output, @Nullable Mask mask, @Nullable DBinaryOperator accum) { + DMatrixRMaj initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accum, output) ? output.copy() : null; + output = reshapeOrDeclare(output, 1, input.numCols); + if (mask != null) { + mask.compatible(output); } for (int col = 0; col < input.numCols; col++) { @@ -1945,14 +1975,18 @@ public static DMatrixRMaj reduceColumnWise(DMatrixSparseCSC input, double initVa double acc = initValue; for (int i = start; i < end; i++) { + // (assumption) not checking for mask here as computation should be faster than check acc = func.apply(acc, input.nz_values[i]); } - // TODO: allow optional resultAccumulator function (use tmp_result array to save reduce result and than combine arrays f.i. 2nd func) output.data[col] = acc; } - return output; + return combineOutputs(output, initialOutput, mask, accum); + } + + public static DMatrixRMaj reduceColumnWise(DMatrixSparseCSC input, double initValue, DBinaryOperator func, @Nullable DMatrixRMaj output) { + return reduceColumnWise(input, initValue, func, output, null, null); } /** @@ -1971,25 +2005,34 @@ public static DMatrixRMaj reduceColumnWise(DMatrixSparseCSC input, double initVa * @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 DMatrixRMaj reduceRowWise(DMatrixSparseCSC input, double initValue, DBinaryOperator func, @Nullable DMatrixRMaj output) { - if (output == null) { - output = new DMatrixRMaj(1, input.numRows); - } else { - output.reshape(1, input.numCols); + public static DMatrixRMaj reduceRowWise(DMatrixSparseCSC input, double initValue, DBinaryOperator func, + @Nullable DMatrixRMaj output, @Nullable Mask mask, @Nullable DBinaryOperator accum) { + DMatrixRMaj initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accum, output) ? output.copy() : null; + output = reshapeOrDeclare(output, input.numRows, 1); + if (mask != null) { + mask.compatible(output); } - // TODO: allow optional resultAccumulator function (use tmp_result array to save reduce result and than combine arrays f.i. 2nd func) + Arrays.fill(output.data, initValue); + // TODO here it might be cheaper to have a bitset for the mask as isSet() gets called multiple times for each entry + //boolean[] isAssigned = mask.createBitSet(); + 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++) { + // (assumption) not checking for mask here as computation should be faster than check output.data[input.nz_rows[i]] = func.apply(output.data[input.nz_rows[i]], input.nz_values[i]); } } - return output; + return combineOutputs(output, initialOutput, mask, accum); + } + + public static DMatrixRMaj reduceRowWise(DMatrixSparseCSC input, double initValue, DBinaryOperator func, @Nullable DMatrixRMaj output) { + return reduceRowWise(input, initValue, func, output, null, null); } } diff --git a/main/ejml-dsparse/src/org/ejml/sparse/csc/MaskUtil_DSCC.java b/main/ejml-dsparse/src/org/ejml/sparse/csc/MaskUtil_DSCC.java new file mode 100644 index 000000000..3ec8d864d --- /dev/null +++ b/main/ejml-dsparse/src/org/ejml/sparse/csc/MaskUtil_DSCC.java @@ -0,0 +1,190 @@ +/* + * Copyright (c) 2009-2020, 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 org.ejml.data.*; +import org.ejml.masks.Mask; +import org.ejml.masks.PrimitiveDMask; +import org.ejml.ops.DBinaryOperator; + +import javax.annotation.Nullable; + +import static org.ejml.UtilEjml.adjust; +import static org.ejml.UtilEjml.checkSameShape; + +public class MaskUtil_DSCC { + private static final DBinaryOperator SECOND = (x, y) -> y; + + // for applying mask and accumulator (output gets overwritten) + // ! assumes that the mask is already applied to the output .. e.g. unset fields not even computed (and not assigned) + // ! the mask is only needed for `apply` there the mask is not applied to the result structure before + public static DMatrixSparseCSC combineOutputs(DMatrixSparseCSC output, DMatrixSparseCSC initialOutput, @Nullable Mask mask, @Nullable DBinaryOperator accum) { + if (initialOutput != null) { + if(accum == null) { + // e.g. just take the newly computed value + accum = SECOND; + } + + // memory overhead .. maybe also can reuse something? + DMatrixSparseCSC combinedOutput = output.createLike(); + // instead of "semiRing.add" this could be a dedicated accumulator + add(initialOutput, output, combinedOutput, mask, accum, null, null); + // is the previous result of C gc-able? (should be) + output = combinedOutput; + } + return output; + } + + static DMatrixSparseCSC combineOutputs(DMatrixSparseCSC output, DBinaryOperator accum, DMatrixSparseCSC initialOutput) { + return combineOutputs(output, initialOutput, null, accum); + } + + static DMatrixRMaj combineOutputs(DMatrixRMaj output, DMatrixRMaj initialOutput, Mask mask, @Nullable DBinaryOperator accum) { + if (initialOutput != null) { + checkSameShape(initialOutput, output, true); + + if(accum == null) { + // e.g. just take the newly computed value + accum = SECOND; + } + + // TODO: operate on a bitset/boolean[] here -> also just one for-loop needed + for (int col = 0; col < output.getNumCols(); col++) { + for (int row = 0; row < output.numRows; row++) { + if (mask.isSet(row, col)) { + output.unsafe_set(row, col, accum.apply(initialOutput.get(row, col), output.get(row, col))); + } else { + // just use previous value as it shouldnt be computed in the first place + output.unsafe_set(row, col, initialOutput.get(row, col)); + } + } + } + } + return output; + } + + public static double[] combineOutputs(@Nullable double[] initialOutput, double[] output, @Nullable PrimitiveDMask mask, @Nullable DBinaryOperator accum) { + if (initialOutput != null) { + if(accum == null) { + // e.g. just take the newly computed value + accum = SECOND; + } + + for (int i = 0; i < output.length; i++) { + if (mask == null || mask.isSet(i)) { + output[i] = accum.apply(initialOutput[i], output[i]); + } + else { + output[i] = initialOutput[i]; + } + } + } + + return output; + } + + /** + * Performs matrix addition:
+ * C = A + B + * + * .. f.i. to combine intialOutput and computedOutput (needed for masks and accumulators in GraphHBLAS) + * @param A Matrix + * @param B Matrix + * @param C Output matrix. + * @param accum accumulator + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static void add(DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC C, @Nullable Mask mask, DBinaryOperator accum, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + double[] x = adjust(gx, A.numRows); + int[] w = adjust(gw, A.numRows, A.numRows); + + C.indicesSorted = false; + C.nz_length = 0; + + for (int col = 0; col < A.numCols; col++) { + C.col_idx[col] = C.nz_length; + + // always take the values of A + // second as x[row] would be the first argument + multAddColA(A, col, C, col + 1, null, (a, b) -> b, x, w); + // accum values of A and B (if entry is set in mask) + multAddColA(B, col, C, col + 1, mask, accum, x, w); + + // take the values in the dense vector 'x' and put them into 'C' + int idxC0 = C.col_idx[col]; + int idxC1 = C.col_idx[col + 1]; + + for (int i = idxC0; i < idxC1; i++) { + C.nz_values[i] = x[C.nz_rows[i]]; + } + } + C.col_idx[A.numCols] = C.nz_length; + } + + /** + * Performs the performing operation x = x + A(:,i) + * for applying a accumulator + */ + public static void multAddColA(DMatrixSparseCSC A, int colA, + DMatrixSparseCSC C, int mark, + @Nullable Mask mask, + DBinaryOperator accum, + double x[], int w[]) { + int idxA0 = A.col_idx[colA]; + int idxA1 = A.col_idx[colA + 1]; + + for (int j = idxA0; j < idxA1; j++) { + int row = A.nz_rows[j]; + if (mask == null || mask.isSet(row, mark - 1)) { + if (w[row] < mark) { + if (C.nz_length >= C.nz_rows.length) { + C.growMaxLength(C.nz_length * 2 + 1, true); + } + + w[row] = mark; + C.nz_rows[C.nz_length] = row; + C.col_idx[mark] = ++C.nz_length; + x[row] = A.nz_values[j]; + } else if (accum != null) { + // if it is null .. x[row] can just stay the same + x[row] = accum.apply(x[row], A.nz_values[j]); + } + } + } + } + + // TODO: remove replace flag and just decide based on accumulator and intialOutput? (replace flag seems redundant) + // drawback: user has to specify accumulator everytime (e.g. not default to SECOND) + + /** + * Check if initialOutput matrix needs to be cached for later merge with actual result + * + * ! mask.replace takes precedence before existing accumlator + */ + public static boolean useInitialOutput(Mask mask, DBinaryOperator accumulator, Matrix initialOutput) { + return initialOutput != null && ((mask != null && !mask.replace) || (accumulator != null && mask == null)); + } + + // primitive version + public static boolean useInitialOutput(Mask mask, DBinaryOperator accumulator, double[] initialOutput) { + return initialOutput != null && ((mask != null && !mask.replace) || (accumulator != null && mask == null)); + } +} diff --git a/main/ejml-dsparse/src/org/ejml/sparse/csc/graphAlgos/Bfs_DSCC.java b/main/ejml-dsparse/src/org/ejml/sparse/csc/graphAlgos/Bfs_DSCC.java new file mode 100644 index 000000000..0cb1dafb9 --- /dev/null +++ b/main/ejml-dsparse/src/org/ejml/sparse/csc/graphAlgos/Bfs_DSCC.java @@ -0,0 +1,269 @@ +/* + * Copyright (c) 2009-2020, 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.graphAlgos; + +import org.ejml.data.DGrowArray; +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.data.IGrowArray; +import org.ejml.masks.DMasks; +import org.ejml.masks.Mask; +import org.ejml.masks.PrimitiveDMask; +import org.ejml.ops.DSemiRing; +import org.ejml.ops.DSemiRings; +import org.ejml.sparse.csc.CommonOpsWithSemiRing_DSCC; +import org.ejml.sparse.csc.CommonOps_DSCC; +import org.ejml.sparse.csc.MaskUtil_DSCC; +import org.ejml.sparse.csc.mult.MatrixVectorMultWithSemiRing_DSCC; + +import java.util.Arrays; + +// variants: boolean/parents/level/multi-bfs + sparse/dense result vector +public class Bfs_DSCC { + + private static DSemiRing getSemiRing(BfsVariation variation) { + return variation == BfsVariation.PARENTS ? DSemiRings.MIN_FIRST : DSemiRings.OR_AND; + } + // TODO: try to reuse working arrays/matrices e.g. gw, gx (and potentially initialOutput vector/array) + // TODO: also return needed iterations + // TODO: is the tmp iterationResult really necessary? (just the inputVector could be enough) + + + // TODO: use dense matrix instead of pure primitive array (f.i. to use `apply` and to support MSBFS) + public BfsDenseResult computeDense(DMatrixSparseCSC adjacencyMatrix, BfsVariation bfsVariation, int startNode, int maxIterations) { + DSemiRing semiRing = getSemiRing(bfsVariation); + double[] result = new double[adjacencyMatrix.numCols]; + Arrays.fill(result, semiRing.add.id); + + if (bfsVariation == BfsVariation.PARENTS) { + result[startNode] = startNode + 1; + } else { + result[startNode] = 1; + } + + + // or use dense matrix and reduceScalar to count non-zero elements + double[] iterationResult = new double[adjacencyMatrix.numCols]; + + int visitedNodes = 1; + int prevVisitedNodes; + + double[] inputVector = result.clone(); + boolean isFixPoint = false; + int iteration = 1; + + for (; (iteration <= maxIterations) && !isFixPoint; iteration++) { + // negated -> dont compute values for visited nodes + // replace -> iterationResult is basically the new inputVector + PrimitiveDMask mask = DMasks.builder(result).withZeroElement(semiRing.add.id).withNegated(true).withReplace(true).build(); + // clear iterationsResult to only contain newly discovered nodes + Arrays.fill(iterationResult, semiRing.add.id); + iterationResult = MatrixVectorMultWithSemiRing_DSCC.mult(inputVector, adjacencyMatrix, iterationResult, semiRing, mask, null); + + prevVisitedNodes = visitedNodes; + + // add newly visited nodes + for (double v : iterationResult) { + if (v != semiRing.add.id) { + visitedNodes++; + } + } + + inputVector = iterationResult.clone(); + + if (bfsVariation == BfsVariation.LEVEL) { + // poor mans apply on primitive array + for (int i = 0; i < iterationResult.length; i++) { + if (iterationResult[i] != semiRing.add.id) { + iterationResult[i] = iteration + 1; + } + } + } + + if (bfsVariation == BfsVariation.PARENTS) { + for (int i = 0; i < inputVector.length; i++) { + if (inputVector[i] != semiRing.add.id) { + inputVector[i] = i + 1; + } + } + } + + // FIXME: avoid cloning .. have a combine method that writes into the intialResult + result = MaskUtil_DSCC.combineOutputs(result, iterationResult.clone(), mask, null); + + isFixPoint = (visitedNodes == prevVisitedNodes) || (visitedNodes == adjacencyMatrix.numCols); + } + + return new BfsDenseResult(result, iteration - 1, semiRing.add.id); + } + + public BfsSparseResult computeSparse(DMatrixSparseCSC adjacencyMatrix, BfsVariation bfsVariation, int[] startNodes, int maxIterations) { + // TODO: use transposed result matrix as startNodes.length << adjacencyMatrix.length + // need to transpose result of VxM before combining + DMatrixSparseCSC result = new DMatrixSparseCSC(startNodes.length, adjacencyMatrix.numCols); + // DMatrixSparseCSC iterationResult = result.createLike(); + + // init result vector + for (int i = 0; i < startNodes.length; i++) { + if (bfsVariation == BfsVariation.PARENTS) { + result.set(0, startNodes[i], i + 1); + } else { + result.set(0, startNodes[i], 1); + } + } + + DMatrixSparseCSC inputVector = result.copy(); + DMatrixSparseCSC iterationResult = null; + + // for reusing memory + IGrowArray gw = new IGrowArray(); + DGrowArray gx = new DGrowArray(); + + int visitedNodes = startNodes.length; + int prevVisitedNodes; + + DSemiRing semiRing = getSemiRing(bfsVariation); + + boolean isFixPoint = false; + int iteration = 1; + + for (; (iteration <= maxIterations) && !isFixPoint; iteration++) { + // negated -> dont compute values for visited nodes + // replace -> iterationResult is basically the new inputVector + Mask mask = DMasks.builder(result, true).withNegated(true).withReplace(true).build(); + iterationResult = CommonOpsWithSemiRing_DSCC.mult(inputVector, adjacencyMatrix, iterationResult, semiRing, mask, null, gw, gx); + prevVisitedNodes = visitedNodes; + + if (mask.replace) { + visitedNodes += iterationResult.nz_length; + } else { + visitedNodes = iterationResult.nz_length; + } + + // set inputVector based on newly discovered nodes + // TODO: do this via an `assign` that supports a mask + inputVector = iterationResult.copy(); + + if (bfsVariation == BfsVariation.LEVEL) { + int currentIteration = iteration + 1; + CommonOps_DSCC.apply(iterationResult, x -> currentIteration); + } + + if (bfsVariation == BfsVariation.PARENTS) { + // TODO: generalize to apply for a (row, col, value) -> newValue + // set value to its own id + for (int col = 0; col < inputVector.numCols; col++) { + int idx = inputVector.col_idx[col]; + int endIdx = inputVector.col_idx[col + 1]; + + for (; idx < endIdx; idx++) { + inputVector.nz_values[idx] = col + 1; + } + } + } + + // combine iterationResult and result (basically poor mans `mask.replace=false`) + result = MaskUtil_DSCC.combineOutputs(result, iterationResult, null, null); + + // TODO dont combine result if its known to be a fixPoint? + isFixPoint = (visitedNodes == prevVisitedNodes) || (visitedNodes == adjacencyMatrix.numCols); + } + + // expect the result to be a row vector / row vectors + return new BfsSparseResult(result, iteration - 1, semiRing.add.id); + } + + + public enum BfsVariation { + BOOLEAN, PARENTS, LEVEL + } + + public interface BfsResult { + int iterations(); + int nodesVisited(); + double get(int nodeId); + } + + public class BfsSparseResult implements BfsResult { + private final DMatrixSparseCSC result; + private final double notFoundValue; + private final int iterations; + + public BfsSparseResult(DMatrixSparseCSC result, int iterations, double fallBackValue) { + this.result = result; + this.iterations = iterations; + notFoundValue = fallBackValue; + } + + @Override + public int iterations() { + return this.iterations; + } + + @Override + public int nodesVisited() { + return this.result.getNonZeroLength(); + } + + @Override + public double get(int nodeId) { + return result.get(0, nodeId, notFoundValue); + } + + public DMatrixSparseCSC result() { + return this.result; + } + } + + public class BfsDenseResult implements BfsResult { + private final double[] result; + private final double notFoundValue; + private final int iterations; + + public BfsDenseResult(double[] result, int iterations, double notFoundValue) { + this.result = result; + this.iterations = iterations; + this.notFoundValue = notFoundValue; + } + + @Override + public int iterations() { + return this.iterations; + } + + @Override + public int nodesVisited() { + int visited = 0; + + for (double v : result) { + if (v != notFoundValue) visited++; + } + + return visited; + } + + @Override + public double get(int nodeId) { + return result[nodeId]; + } + + public double[] result() { + return this.result; + } + } +} diff --git a/main/ejml-dsparse/src/org/ejml/sparse/csc/misc/ImplCommonOpsWithSemiRing_DSCC.java b/main/ejml-dsparse/src/org/ejml/sparse/csc/misc/ImplCommonOpsWithSemiRing_DSCC.java new file mode 100644 index 000000000..cc55c15c8 --- /dev/null +++ b/main/ejml-dsparse/src/org/ejml/sparse/csc/misc/ImplCommonOpsWithSemiRing_DSCC.java @@ -0,0 +1,179 @@ +/* + * Copyright (c) 2009-2020, 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.misc; + +import org.ejml.data.DGrowArray; +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.data.IGrowArray; +import org.ejml.masks.Mask; +import org.ejml.ops.DSemiRing; + +import javax.annotation.Nullable; +import java.util.Arrays; + +import static org.ejml.UtilEjml.adjust; +import static org.ejml.sparse.csc.mult.ImplSparseSparseMultWithSemiRing_DSCC.multAddColA; + +/** + * based on ImplCommonOps_DSCC + */ +public class ImplCommonOpsWithSemiRing_DSCC { + + /** + * 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 C Output matrix. + * @param mask Mask for specifying which entries should be overwritten + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static void add(double alpha, DMatrixSparseCSC A, double beta, DMatrixSparseCSC B, DMatrixSparseCSC C, DSemiRing semiRing, @Nullable Mask mask, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + double[] x = adjust(gx, A.numRows); + int[] w = adjust(gw, A.numRows, A.numRows); + + C.indicesSorted = false; + C.nz_length = 0; + + for (int col = 0; col < A.numCols; col++) { + C.col_idx[col] = C.nz_length; + + multAddColA(A, col, alpha, C, col + 1, semiRing, mask, x, w); + multAddColA(B, col, beta, C, col + 1, semiRing, mask, x, w); + + // take the values in the dense vector 'x' and put them into 'C' + int idxC0 = C.col_idx[col]; + int idxC1 = C.col_idx[col + 1]; + + for (int i = idxC0; i < idxC1; i++) { + C.nz_values[i] = x[C.nz_rows[i]]; + } + } + C.col_idx[A.numCols] = C.nz_length; + } + + /** + * Adds the results of adding a column in A and B as a new column in C.
+ * C(:,end+1) = A(:,colA) + B(:,colB) + * + * @param A matrix + * @param colA column in A + * @param B matrix + * @param colB column in B + * @param C Column in C + * @param gw workspace + */ + public static void addColAppend(DMatrixSparseCSC A, int colA, DMatrixSparseCSC B, int colB, + DMatrixSparseCSC C, DSemiRing semiRing, @Nullable IGrowArray gw) { + if (A.numRows != B.numRows || A.numRows != C.numRows) + throw new IllegalArgumentException("Number of rows in A, B, and C do not match"); + + int idxA0 = A.col_idx[colA]; + int idxA1 = A.col_idx[colA + 1]; + int idxB0 = B.col_idx[colB]; + int idxB1 = B.col_idx[colB + 1]; + + C.growMaxColumns(++C.numCols, true); + C.growMaxLength(C.nz_length + idxA1 - idxA0 + idxB1 - idxB0, true); + + int[] w = adjust(gw, A.numRows); + Arrays.fill(w, 0, A.numRows, -1); + + for (int i = idxA0; i < idxA1; i++) { + int row = A.nz_rows[i]; + C.nz_rows[C.nz_length] = row; + C.nz_values[C.nz_length] = A.nz_values[i]; + w[row] = C.nz_length++; + } + + for (int i = idxB0; i < idxB1; i++) { + int row = B.nz_rows[i]; + if (w[row] != -1) { + C.nz_values[w[row]] = semiRing.add.func.apply(C.nz_values[w[row]], B.nz_values[i]); + } else { + C.nz_values[C.nz_length] = B.nz_values[i]; + C.nz_rows[C.nz_length++] = row; + } + } + C.col_idx[C.numCols] = C.nz_length; + } + + /** + * Performs element-wise multiplication:
+ * C_ij = A_ij * B_ij + * + * @param A (Input) Matrix + * @param B (Input) Matrix + * @param C (Output) Matrix. + * @param mask Mask for specifying which entries should be overwritten + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static void elementMult(DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC C, DSemiRing semiRing, @Nullable Mask mask, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + double[] x = adjust(gx, A.numRows); + int[] w = adjust(gw, A.numRows); + Arrays.fill(w, 0, A.numRows, -1); // fill with -1. This will be a value less than column + + C.growMaxLength(Math.min(A.nz_length, B.nz_length), false); + C.indicesSorted = false; // Hmm I think if B is storted then C will be sorted... + C.nz_length = 0; + + for (int col = 0; col < A.numCols; col++) { + int idxA0 = A.col_idx[col]; + int idxA1 = A.col_idx[col + 1]; + int idxB0 = B.col_idx[col]; + int idxB1 = B.col_idx[col + 1]; + + // compute the maximum number of elements that there can be in this row + int maxInRow = Math.min(idxA1 - idxA0, idxB1 - idxB0); + + // make sure there are enough non-zero elements in C + if (C.nz_length + maxInRow > C.nz_values.length) + C.growMaxLength(C.nz_values.length + maxInRow, true); + + // update the structure of C + C.col_idx[col] = C.nz_length; + + // mark the rows that appear in A and save their value + for (int i = idxA0; i < idxA1; i++) { + int row = A.nz_rows[i]; + w[row] = col; + x[row] = A.nz_values[i]; + } + + // If a row appears in A and B, multiply and set as an element in C + for (int i = idxB0; i < idxB1; i++) { + int row = B.nz_rows[i]; + if (w[row] == col) { + if (mask == null || mask.isSet(row, col)) { + C.nz_values[C.nz_length] = semiRing.mult.func.apply(x[row], B.nz_values[i]); + C.nz_rows[C.nz_length++] = row; + } + } + } + } + C.col_idx[C.numCols] = C.nz_length; + } +} \ No newline at end of file diff --git a/main/ejml-dsparse/src/org/ejml/sparse/csc/mult/ImplSparseSparseMultWithSemiRing_DSCC.java b/main/ejml-dsparse/src/org/ejml/sparse/csc/mult/ImplSparseSparseMultWithSemiRing_DSCC.java new file mode 100644 index 000000000..538eaf35e --- /dev/null +++ b/main/ejml-dsparse/src/org/ejml/sparse/csc/mult/ImplSparseSparseMultWithSemiRing_DSCC.java @@ -0,0 +1,365 @@ +/* + * Copyright (c) 2009-2019, 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.mult; + +import org.ejml.data.DGrowArray; +import org.ejml.data.DMatrixRMaj; +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.data.IGrowArray; +import org.ejml.masks.Mask; +import org.ejml.ops.DSemiRing; +import org.ejml.sparse.csc.CommonOps_DSCC; + +import javax.annotation.Nullable; + +import static org.ejml.UtilEjml.adjust; + +/** + * based on ImplSparseSparseGraphMult_DSCC + */ +public class ImplSparseSparseMultWithSemiRing_DSCC { + + /** + * Performs matrix multiplication. C = A*B + * + * @param A Matrix + * @param B Matrix + * @param C Storage for results. Data length is increased if increased if insufficient. + * @param mask Mask for specifying which entries should be overwritten + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static void mult(DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC C, DSemiRing semiRing, @Nullable Mask mask, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + double[] x = adjust(gx, A.numRows); + int[] w = adjust(gw, A.numRows, A.numRows); + + C.growMaxLength(A.nz_length + B.nz_length, false); + C.indicesSorted = false; + C.nz_length = 0; + + // C(i,j) = sum_k A(i,k) * B(k,j) + int idx0 = B.col_idx[0]; + for (int bj = 1; bj <= B.numCols; bj++) { + int colB = bj - 1; + int idx1 = B.col_idx[bj]; + C.col_idx[bj] = C.nz_length; + + if (idx0 == idx1) { + continue; + } + + // C(:,j) = sum_k A(:,k)*B(k,j) + for (int bi = idx0; bi < idx1; bi++) { + int rowB = B.nz_rows[bi]; + double valB = B.nz_values[bi]; // B(k,j) k=rowB j=colB + + multAddColA(A, rowB, valB, C, colB + 1, semiRing, mask, x, w); + } + + // take the values in the dense vector 'x' and put them into 'C' + int idxC0 = C.col_idx[colB]; + int idxC1 = C.col_idx[colB + 1]; + + for (int i = idxC0; i < idxC1; i++) { + C.nz_values[i] = x[C.nz_rows[i]]; + } + + idx0 = idx1; + } + } + + /** + * Performs matrix multiplication. C = AT*B + * + * @param A Matrix + * @param B Matrix + * @param C Storage for results. Data length is increased if increased if insufficient. + * @param mask Mask for specifying which entries should be overwritten + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static void multTransA(DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC C, DSemiRing semiRing, @Nullable Mask mask, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + double[] x = adjust(gx, A.numRows); + int[] w = adjust(gw, A.numRows, A.numRows); + + C.growMaxLength(A.nz_length + B.nz_length, false); + C.indicesSorted = true; + C.nz_length = 0; + C.col_idx[0] = 0; + + int idxB0 = B.col_idx[0]; + for (int bj = 1; bj <= B.numCols; bj++) { + int idxB1 = B.col_idx[bj]; + C.col_idx[bj] = C.nz_length; + + if (idxB0 == idxB1) { + continue; + } + + // convert the column of B into a dense format and mark which rows are used + for (int bi = idxB0; bi < idxB1; bi++) { + int rowB = B.nz_rows[bi]; + x[rowB] = B.nz_values[bi]; + // TODO: wouldn't a BitSet also work instead of an IGrowArray? + w[rowB] = bj; + } + + // C(colA,colB) = A(:,colA)*B(:,colB) + for (int colA = 0; colA < A.numCols; colA++) { + if (mask == null || mask.isSet(colA, bj - 1)) { + int idxA0 = A.col_idx[colA]; + int idxA1 = A.col_idx[colA + 1]; + + double sum = semiRing.add.id; + for (int ai = idxA0; ai < idxA1; ai++) { + int rowA = A.nz_rows[ai]; + if (w[rowA] == bj) { + sum = semiRing.add.func.apply(sum, semiRing.mult.func.apply(x[rowA], A.nz_values[ai])); + } + } + + if (sum != semiRing.add.id) { + if (C.nz_length == C.nz_values.length) { + C.growMaxLength(C.nz_length * 2 + 1, true); + } + C.nz_values[C.nz_length] = sum; + C.nz_rows[C.nz_length++] = colA; + } + } + } + C.col_idx[bj] = C.nz_length; + idxB0 = idxB1; + } + } + + /** + * Performs matrix multiplication. C = A*BT + * + * @param A Matrix + * @param B Matrix + * @param C Storage for results. Data length is increased if increased if insufficient. + * @param mask Mask for specifying which entries should be overwritten + * @param gw (Optional) Storage for internal workspace. Can be null. + * @param gx (Optional) Storage for internal workspace. Can be null. + */ + public static void multTransB(DMatrixSparseCSC A, DMatrixSparseCSC B, DMatrixSparseCSC C, DSemiRing semiRing, @Nullable Mask mask, + @Nullable IGrowArray gw, @Nullable DGrowArray gx) { + if (!B.isIndicesSorted()) + throw new IllegalArgumentException("B must have its indices sorted."); + else if (!CommonOps_DSCC.checkIndicesSorted(B)) { + throw new IllegalArgumentException("Crap. Not really sorted"); + } + + double[] x = adjust(gx, A.numRows); + int[] w = adjust(gw, A.numRows + B.numCols, A.numRows); + + C.growMaxLength(A.nz_length + B.nz_length, false); + C.indicesSorted = false; + C.nz_length = 0; + C.col_idx[0] = 0; + + // initialize w is the first index in each column of B + int locationB = A.numRows; + System.arraycopy(B.col_idx, 0, w, locationB, B.numCols); + + for (int colC = 0; colC < B.numRows; colC++) { + C.col_idx[colC + 1] = C.nz_length; // needs a value of B has nothing in the row + + // find the column in the transposed B + int mark = colC + 1; + for (int colB = 0; colB < B.numCols; colB++) { + int bi = w[locationB + colB]; + if (bi < B.col_idx[colB + 1]) { + int row = B.nz_rows[bi]; + if (row == colC) { + multAddColA(A, colB, B.nz_values[bi], C, mark, semiRing, mask, x, w); + w[locationB + colB]++; + } + } + } + + // take the values in the dense vector 'x' and put them into 'C' + int idxC0 = C.col_idx[colC]; + int idxC1 = C.col_idx[colC + 1]; + + for (int i = idxC0; i < idxC1; i++) { + C.nz_values[i] = x[C.nz_rows[i]]; + } + } + } + + /** + * Performs the performing operation x = x + A(:,i)*alpha + * + *

NOTE: This is the same as cs_scatter() in csparse.

+ */ + public static void multAddColA(DMatrixSparseCSC A, int colA, + double alpha, + DMatrixSparseCSC C, int mark, + DSemiRing semiRing, @Nullable Mask mask, + double x[], int w[]) { + int idxA0 = A.col_idx[colA]; + int idxA1 = A.col_idx[colA + 1]; + + for (int j = idxA0; j < idxA1; j++) { + int row = A.nz_rows[j]; + // TODO: only use iterator over a single column mask values here + // mark - 1 is the actual target column + if (mask == null || mask.isSet(row, mark - 1)) { + if (w[row] < mark) { + if (C.nz_length >= C.nz_rows.length) { + C.growMaxLength(C.nz_length * 2 + 1, true); + } + + w[row] = mark; + C.nz_rows[C.nz_length] = row; + C.col_idx[mark] = ++C.nz_length; + x[row] = semiRing.mult.func.apply(A.nz_values[j], alpha); + } else { + x[row] = semiRing.add.func.apply(x[row], semiRing.mult.func.apply(A.nz_values[j], alpha)); + } + } + } + } + + // sparse-dense variants + + public static void mult(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + C.fill(semiRing.add.id); + multAdd(A, B, C, semiRing); + } + + public static void multAdd(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + // C(i,j) = sum_k A(i,k) * B(k,j) + for (int k = 0; k < A.numCols; k++) { + int idx0 = A.col_idx[k]; + int idx1 = A.col_idx[k + 1]; + + for (int indexA = idx0; indexA < idx1; indexA++) { + int i = A.nz_rows[indexA]; + double valueA = A.nz_values[indexA]; + + int indexB = k * B.numCols; + int indexC = i * C.numCols; + int end = indexB + B.numCols; + +// for (int j = 0; j < B.numCols; j++) { + while (indexB < end) { + C.data[indexC++] = semiRing.add.func.apply(C.data[indexC++], semiRing.mult.func.apply(valueA, B.data[indexB++])); + } + } + } + } + + public static void multTransA(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + + // C(i,j) = sum_k A(k,i) * B(k,j) + for (int j = 0; j < B.numCols; j++) { + + for (int i = 0; i < A.numCols; i++) { + int idx0 = A.col_idx[i]; + int idx1 = A.col_idx[i + 1]; + + double sum = semiRing.add.id; + for (int indexA = idx0; indexA < idx1; indexA++) { + int rowK = A.nz_rows[indexA]; + sum = semiRing.add.func.apply(sum, semiRing.mult.func.apply(A.nz_values[indexA], B.data[rowK * B.numCols + j])); + } + + C.data[i * C.numCols + j] = sum; + } + } + } + + public static void multAddTransA(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + // C(i,j) = sum_k A(k,i) * B(k,j) + for (int j = 0; j < B.numCols; j++) { + + for (int i = 0; i < A.numCols; i++) { + int idx0 = A.col_idx[i]; + int idx1 = A.col_idx[i + 1]; + + double sum = semiRing.add.id; + for (int indexA = idx0; indexA < idx1; indexA++) { + int rowK = A.nz_rows[indexA]; + sum = semiRing.add.func.apply(sum, semiRing.mult.func.apply(A.nz_values[indexA], B.data[rowK * B.numCols + j])); + } + + C.data[i * C.numCols + j] = semiRing.add.func.apply(C.data[i * C.numCols + j], sum); + } + } + } + + public static void multTransB(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + C.zero(); + multAddTransB(A, B, C, semiRing); + } + + public static void multAddTransB(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + + // C(i,j) = sum_k A(i,k) * B(j,k) + for (int k = 0; k < A.numCols; k++) { + int idx0 = A.col_idx[k]; + int idx1 = A.col_idx[k + 1]; + for (int indexA = idx0; indexA < idx1; indexA++) { + for (int j = 0; j < B.numRows; j++) { + int i = A.nz_rows[indexA]; + C.data[i * C.numCols + j] = semiRing.add.func.apply( + C.data[i * C.numCols + j], + semiRing.mult.func.apply( + A.nz_values[indexA], + B.data[j * B.numCols + k] + ) + ); + } + } + } + } + + public static void multTransAB(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + C.zero(); + multAddTransAB(A, B, C, semiRing); + } + + public static void multAddTransAB(DMatrixSparseCSC A, DMatrixRMaj B, DMatrixRMaj C, DSemiRing semiRing) { + // C(i,j) = sum_k A(k,i) * B(j,K) + for (int i = 0; i < A.numCols; i++) { + int idx0 = A.col_idx[i]; + int idx1 = A.col_idx[i + 1]; + + for (int indexA = idx0; indexA < idx1; indexA++) { + for (int j = 0; j < B.numRows; j++) { + int indexB = j * B.numCols; + + int k = A.nz_rows[indexA]; + + C.data[i * C.numCols + j] = semiRing.add.func.apply( + C.data[i * C.numCols + j], + semiRing.mult.func.apply( + A.nz_values[indexA], + B.data[indexB + k] + ) + ); + } + } + } + } +} \ No newline at end of file diff --git a/main/ejml-dsparse/src/org/ejml/sparse/csc/mult/MatrixVectorMultWithSemiRing_DSCC.java b/main/ejml-dsparse/src/org/ejml/sparse/csc/mult/MatrixVectorMultWithSemiRing_DSCC.java new file mode 100644 index 000000000..b5a91e3ef --- /dev/null +++ b/main/ejml-dsparse/src/org/ejml/sparse/csc/mult/MatrixVectorMultWithSemiRing_DSCC.java @@ -0,0 +1,156 @@ +/* + * Copyright (c) 2009-2018, 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.mult; + +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.masks.PrimitiveDMask; +import org.ejml.ops.DBinaryOperator; +import org.ejml.ops.DSemiRing; +import org.ejml.sparse.csc.MaskUtil_DSCC; + +import javax.annotation.Nullable; +import java.util.Arrays; + +/** + * based on MartrixVectorMult_DSCC + */ +public class MatrixVectorMultWithSemiRing_DSCC { + /** + * output = A*b + * + * @param A (Input) Matrix + * @param b (Input) vector + * @param output (Output) vector + * @param mask Mask for specifying which entries should be overwritten + */ + public static double[] mult(DMatrixSparseCSC A, double b[], double output[], DSemiRing semiRing, + @Nullable PrimitiveDMask mask, @Nullable DBinaryOperator accumulator) { + double[] initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accumulator, output) ? output.clone() : null; + if (mask != null) { + mask.compatible(output); + } + + Arrays.fill(output, semiRing.add.id); + + return multAdd(A, b, output, initialOutput, semiRing, mask, accumulator); + } + + public static double[] mult(DMatrixSparseCSC A, double b[], double output[], DSemiRing semiRing) { + return mult(A, b, output, semiRing, null, null); + } + + /** + * output = initialOutput + A*b + * ! difference to normal `multAdd` without semiRing: apply + on result of A*b + * + * @param A (Input) Matrix + * @param b (Input) vector + * @param output (Output) vector + * @param mask Mask for specifying which entries should be overwritten + * @param semiRing + * @param accumulator (Optional) accumulator for output + (A*b), elso use `add` from the semiRing + */ + public static double[] multAdd(DMatrixSparseCSC A, double[] b, double[] output, @Nullable double[] initialOutput, + DSemiRing semiRing, @Nullable PrimitiveDMask mask, @Nullable DBinaryOperator accumulator) { + for (int k = 0; k < A.numCols; k++) { + int idx0 = A.col_idx[k]; + int idx1 = A.col_idx[k + 1]; + + for (int indexA = idx0; indexA < idx1; indexA++) { + if (mask == null || mask.isSet(A.nz_rows[indexA])) { + output[A.nz_rows[indexA]] = semiRing.add.func.apply( + output[A.nz_rows[indexA]], + semiRing.mult.func.apply(A.nz_values[indexA], b[k])); + } + } + } + + // initialOutput + output + return MaskUtil_DSCC.combineOutputs(initialOutput, output, mask, accumulator); + } + + /** + * output = aT*B + * + * @param a (Input) vector + * @param B (Input) Matrix + * @param output (Output) vector + * @param mask Mask for specifying which entries should be overwritten + */ + public static double[] mult(double a[], DMatrixSparseCSC B, double output[], DSemiRing semiRing, + @Nullable PrimitiveDMask mask, @Nullable DBinaryOperator accumulator) { + double[] initialOutput = MaskUtil_DSCC.useInitialOutput(mask, accumulator, output) ? output.clone() : null; + if (mask != null) { + mask.compatible(output); + } + + for (int k = 0; k < B.numCols; k++) { + if (mask == null || mask.isSet(k)) { + int idx0 = B.col_idx[k]; + int idx1 = B.col_idx[k + 1]; + + + double sum = semiRing.add.id; + for (int indexB = idx0; indexB < idx1; indexB++) { + sum = semiRing.add.func.apply(sum, semiRing.mult.func.apply(a[B.nz_rows[indexB]], B.nz_values[indexB])); + } + output[k] = sum; + } + } + + return MaskUtil_DSCC.combineOutputs(initialOutput, output, mask, accumulator); + } + + public static double[] mult(double a[], DMatrixSparseCSC B, double c[], DSemiRing semiRing) { + return mult(a, B, c, semiRing, null, null); + } + + /** + * scalar = AT*B*C + * + * @param a (Input) vector + * @param offsetA Input) first index in vector a + * @param B (Input) Matrix + * @param c (Output) vector + * @param offsetC (Output) first index in vector c + */ + public static double innerProduct(double a[], int offsetA, + DMatrixSparseCSC B, + double c[], int offsetC, DSemiRing semiRing) { + if (a.length - offsetA < B.numRows) + throw new IllegalArgumentException("Length of 'a' isn't long enough"); + if (c.length - offsetC < B.numCols) + throw new IllegalArgumentException("Length of 'c' isn't long enough"); + + double output = 0; + + for (int k = 0; k < B.numCols; k++) { + int idx0 = B.col_idx[k]; + int idx1 = B.col_idx[k + 1]; + + double sum = 0; + for (int indexB = idx0; indexB < idx1; indexB++) { + sum = semiRing.add.func.apply(sum, semiRing.mult.func.apply(a[offsetA + B.nz_rows[indexB]], B.nz_values[indexB])); + } + output = semiRing.add.func.apply(output, semiRing.mult.func.apply(sum, c[offsetC + k])); + } + + return output; + } +} diff --git a/main/ejml-dsparse/test/org/ejml/sparse/csc/graphAlgos/BfsTest_DSCC.java b/main/ejml-dsparse/test/org/ejml/sparse/csc/graphAlgos/BfsTest_DSCC.java new file mode 100644 index 000000000..d87122fbe --- /dev/null +++ b/main/ejml-dsparse/test/org/ejml/sparse/csc/graphAlgos/BfsTest_DSCC.java @@ -0,0 +1,104 @@ +/* + * Copyright (c) 2009-2020, 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.graphAlgos; + +import org.ejml.EjmlUnitTests; +import org.ejml.data.DMatrixRMaj; +import org.ejml.data.DMatrixSparseCSC; +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.api.Test; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.Arguments; +import org.junit.jupiter.params.provider.MethodSource; + +import java.util.Arrays; +import java.util.stream.Stream; + +import static org.ejml.sparse.csc.graphAlgos.Bfs_DSCC.*; +import static org.junit.jupiter.api.Assertions.*; + +public class BfsTest_DSCC { + DMatrixSparseCSC inputMatrix; + + Bfs_DSCC bfs = new Bfs_DSCC(); + + private static Stream bfsVariantSource() { + + return Stream.of( + Arguments.of(BfsVariation.BOOLEAN, new double[] {1,1,1,1,1,1,1}), + Arguments.of(BfsVariation.LEVEL, new double[] {1,2,3,2,3,4,3}), + Arguments.of(BfsVariation.PARENTS, new double[] {1,1,4,1,2,3,2}) + + ); + } + + @BeforeEach + public void setUp() { + // based on example in http://mit.bme.hu/~szarnyas/grb/graphblas-introduction.pdf + inputMatrix = new DMatrixSparseCSC(7, 7, 12); + inputMatrix.set(0, 1, 1); + inputMatrix.set(0, 3, 1); + inputMatrix.set(1, 4, 1); + inputMatrix.set(1, 6, 1); + inputMatrix.set(2, 5, 1); + inputMatrix.set(3, 0, 0.2); + inputMatrix.set(3, 2, 0.4); + inputMatrix.set(4, 5, 1); + inputMatrix.set(5, 2, 0.5); + inputMatrix.set(6, 2, 1); + inputMatrix.set(6, 3, 1); + inputMatrix.set(6, 4, 1); + } + + @ParameterizedTest + @MethodSource("bfsVariantSource") + public void testSparseVariations(BfsVariation variation, double[] expected) { + int[] startNodes = {0}; + int maxIterations = 20; + Bfs_DSCC.BfsSparseResult result = bfs.computeSparse(inputMatrix, variation, startNodes, maxIterations); + + DMatrixRMaj expectedMatrix = new DMatrixRMaj(1, inputMatrix.numCols, true, expected); + + assertEquals(result.iterations(), 3); + EjmlUnitTests.assertEquals(expectedMatrix, result.result()); + } + + @ParameterizedTest + @MethodSource("bfsVariantSource") + public void testDenseVariations(BfsVariation variation, double[] expected) { + int startNode = 0; + int maxIterations = 20; + + Bfs_DSCC.BfsDenseResult result = bfs.computeDense(inputMatrix, variation, startNode, maxIterations); + + assertEquals(result.iterations(), 3); + assertTrue(Arrays.equals(expected, result.result())); + } + + @Test + public void testEmptyResult() { + BfsResult denseIt = bfs.computeDense(new DMatrixSparseCSC(2, 2), BfsVariation.LEVEL, 0, 5); + + int[] startNodes = {0}; + BfsResult sparseIt = bfs.computeSparse(new DMatrixSparseCSC(2, 2), BfsVariation.LEVEL, startNodes, 5); + + assertEquals(denseIt.iterations(), sparseIt.iterations()); + assertEquals(denseIt.nodesVisited(), sparseIt.nodesVisited()); + } +} diff --git a/main/ejml-dsparse/test/org/ejml/sparse/csc/misc/TestImplCommonOpsWithSemiRing_DSCC.java b/main/ejml-dsparse/test/org/ejml/sparse/csc/misc/TestImplCommonOpsWithSemiRing_DSCC.java new file mode 100644 index 000000000..4301b63e3 --- /dev/null +++ b/main/ejml-dsparse/test/org/ejml/sparse/csc/misc/TestImplCommonOpsWithSemiRing_DSCC.java @@ -0,0 +1,98 @@ +/* + * Copyright (c) 2009-2020, 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.misc; + +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.ops.DSemiRing; +import org.ejml.ops.DSemiRings; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.Arguments; +import org.junit.jupiter.params.provider.MethodSource; + +import java.util.Arrays; +import java.util.stream.Stream; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertTrue; + +public class TestImplCommonOpsWithSemiRing_DSCC { + + @ParameterizedTest + @MethodSource("elementWiseAddSource") + public void add(DSemiRing semiRing, double[] expected) { + // == graph unions + DMatrixSparseCSC a = new DMatrixSparseCSC(3, 3); + DMatrixSparseCSC b = a.copy(); + DMatrixSparseCSC c = a.copy(); + + a.set(1, 1, 2); + + b.set(1, 1, 3); + b.set(0, 0, 4); + + ImplCommonOpsWithSemiRing_DSCC.add(1, a, 1, b, c, semiRing, null, null, null); + + double[] found = new double[]{c.get(0, 0), c.get(1, 1)}; + assertTrue(c.getNumElements() == 2); + + assertTrue(Arrays.equals(expected, found)); + } + + @ParameterizedTest + @MethodSource("elementWiseMultSource") + public void testelementWiseMult(DSemiRing semiRing, double[] expected) { + // == graph intersection + DMatrixSparseCSC matrix = new DMatrixSparseCSC(3, 3, 4); + matrix.set(1, 1, 4); + matrix.set(1, 2, -2); + + DMatrixSparseCSC otherMatrix = matrix.copy(); + otherMatrix.set(1, 1, 3); + otherMatrix.set(1, 2, 1); + + + matrix.set(0, 2, 1); + otherMatrix.set(2, 0, 1); + + DMatrixSparseCSC result = new DMatrixSparseCSC(3, 3, 0); + + ImplCommonOpsWithSemiRing_DSCC.elementMult(matrix, otherMatrix, result, semiRing, null, null, null); + + assertEquals(2, result.getNumElements()); + assertTrue(expected[0] == result.get(1, 1)); + assertTrue(expected[1] == result.get(1, 2)); + } + + private static Stream elementWiseAddSource() { + return Stream.of( + Arguments.of(DSemiRings.PLUS_TIMES, new double[]{4, 5}), + Arguments.of(DSemiRings.MIN_MAX, new double[]{4, 2}), + Arguments.of(DSemiRings.OR_AND, new double[]{1, 1}) + ); + } + + private static Stream elementWiseMultSource() { + return Stream.of( + Arguments.of(DSemiRings.PLUS_TIMES, new double[]{12, -2}), + Arguments.of(DSemiRings.PLUS_MIN, new double[]{3, -2}), + Arguments.of(DSemiRings.MIN_MAX, new double[]{4, 1}), + Arguments.of(DSemiRings.OR_AND, new double[]{1, 1}) + ); + } +} diff --git a/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/BaseTestMatrixMatrixOpsWithSemiRing_DSCC.java b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/BaseTestMatrixMatrixOpsWithSemiRing_DSCC.java new file mode 100644 index 000000000..8c458278a --- /dev/null +++ b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/BaseTestMatrixMatrixOpsWithSemiRing_DSCC.java @@ -0,0 +1,45 @@ +/* + * Copyright (c) 2009-2020, 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.mult; + +import org.ejml.data.DMatrixSparseCSC; +import org.junit.jupiter.api.BeforeEach; + +public class BaseTestMatrixMatrixOpsWithSemiRing_DSCC { + DMatrixSparseCSC inputMatrix; + + @BeforeEach + public void setUp() { + // based on example in http://mit.bme.hu/~szarnyas/grb/graphblas-introduction.pdf + inputMatrix = new DMatrixSparseCSC(7, 7, 12); + inputMatrix.set(0, 1, 1); + inputMatrix.set(0, 3, 1); + inputMatrix.set(1, 4, 1); + inputMatrix.set(1, 6, 1); + inputMatrix.set(2, 5, 1); + inputMatrix.set(3, 0, 0.2); + inputMatrix.set(3, 2, 0.4); + inputMatrix.set(4, 5, 1); + inputMatrix.set(5, 2, 0.5); + inputMatrix.set(6, 2, 1); + inputMatrix.set(6, 3, 1); + inputMatrix.set(6, 4, 1); + + } +} \ No newline at end of file diff --git a/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMaskedOperators_DSCC.java b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMaskedOperators_DSCC.java new file mode 100644 index 000000000..9e6e09728 --- /dev/null +++ b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMaskedOperators_DSCC.java @@ -0,0 +1,279 @@ +/* + * Copyright (c) 2009-2020, 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.mult; + +import org.ejml.data.*; +import org.ejml.masks.DMasks; +import org.ejml.masks.Mask; +import org.ejml.masks.PrimitiveDMask; +import org.ejml.ops.DBinaryOperator; +import org.ejml.ops.DSemiRing; +import org.ejml.ops.DSemiRings; +import org.ejml.sparse.csc.CommonOpsWithSemiRing_DSCC; +import org.ejml.sparse.csc.CommonOps_DSCC; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.Arguments; +import org.junit.jupiter.params.provider.MethodSource; +import org.junit.jupiter.params.provider.ValueSource; + +import java.util.Arrays; +import java.util.Iterator; +import java.util.stream.Stream; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +public class TestMaskedOperators_DSCC extends BaseTestMatrixMatrixOpsWithSemiRing_DSCC { + // TODO test with replace = false and generate all possible masks (not hard-coded) + // test accumulator + static DSemiRing semiRing = DSemiRings.PLUS_TIMES; + + private static Stream primitiveVectorSource() { + double[] inputVector = new double[7]; + Arrays.fill(inputVector, semiRing.add.id); + inputVector[3] = 0.5; + inputVector[0] = 0.6; + + double[] prevResult = new double[7]; + prevResult[2] = 99; + prevResult[0] = 42; + + PrimitiveDMask.Builder maskBuilder = new PrimitiveDMask.Builder(prevResult); + + return Stream.of( + Arguments.of(inputVector, prevResult, maskBuilder.withNegated(true).withReplace(false).build()), + Arguments.of(inputVector, prevResult, maskBuilder.withNegated(false).withReplace(false).build()) + ); + } + + private static Stream sparseVectorSource() { + DMatrixSparseCSC otherMatrix = new DMatrixSparseCSC(1, 7); + otherMatrix.set(0, 3, 0.5); + otherMatrix.set(0, 4, 0.6); + + DMatrixSparseCSC prevResult = otherMatrix.copy(); + prevResult.set(0, 2, 99); + + return Stream.of( + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, false).withNegated(true).withReplace(false).build()), + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, true).withNegated(false).withReplace(false).build()), + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, false).withNegated(false).withReplace(false).build()), + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, true).withNegated(true).withReplace(false).build()) + ); + } + + private static Stream sparseMatrixSource() { + DMatrixSparseCSC otherMatrix = new DMatrixSparseCSC(7, 7); + otherMatrix.set(0, 3, 0.5); + otherMatrix.set(0, 5, 0.6); + + DMatrixSparseCSC prevResult = new DMatrixSparseCSC(7, 7); + prevResult.set(0, 0, 99); + prevResult.set(0, 3, 42); + + return Stream.of( + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, false).withNegated(true).withReplace(false).build()), + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, true).withNegated(false).withReplace(false).build()), + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, false).withNegated(false).withReplace(false).build()), + Arguments.of(otherMatrix, prevResult, DMasks.builder(prevResult, true).withNegated(true).withReplace(false).build()) + ); + } + + @ParameterizedTest + @MethodSource("primitiveVectorSource") + public void mult_v_A(double[] inputVector, double[] prevResult, PrimitiveDMask mask) { + double[] found = prevResult.clone(); + double[] foundWithMask = prevResult.clone(); + + found = MatrixVectorMultWithSemiRing_DSCC.mult(inputVector, inputMatrix, found, semiRing); + foundWithMask = MatrixVectorMultWithSemiRing_DSCC.mult(inputVector, inputMatrix, foundWithMask, semiRing, mask, null); + + assertMaskedResult(prevResult, found, foundWithMask, mask); + } + + @ParameterizedTest + @MethodSource("primitiveVectorSource") + public void mult_A_v(double[] inputVector, double[] prevResult, PrimitiveDMask mask) { + double[] found = prevResult.clone(); + double[] foundWithMask = prevResult.clone(); + + found = MatrixVectorMultWithSemiRing_DSCC.mult(inputMatrix, inputVector, found, semiRing); + foundWithMask = MatrixVectorMultWithSemiRing_DSCC.mult(inputMatrix, inputVector, foundWithMask, semiRing, mask, null); + + assertMaskedResult(prevResult, found, foundWithMask, mask); + } + + // matrix, matrix ops + + @ParameterizedTest + @MethodSource("sparseVectorSource") + public void mult_A_B(DMatrixSparseCSC sparseVector, DMatrixSparseCSC prevResult, Mask mask) { + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.mult(sparseVector, inputMatrix, prevResult.copy(), semiRing, null, null); + DMatrixSparseCSC foundWithMask = CommonOpsWithSemiRing_DSCC.mult(sparseVector, inputMatrix, prevResult.copy(), semiRing, mask, null); + + assertMaskedResult(prevResult, found, foundWithMask, mask); + } + + @ParameterizedTest + @MethodSource("sparseVectorSource") + public void mult_T_A_B(DMatrixSparseCSC sparseVector, DMatrixSparseCSC prevResult, Mask mask) { + DMatrixSparseCSC transposed_vector = CommonOps_DSCC.transpose(sparseVector, null, null); + + DMatrixSparseCSC foundTA = CommonOpsWithSemiRing_DSCC.multTransA( + transposed_vector, inputMatrix, prevResult.copy(), semiRing, null, null, null, null); + DMatrixSparseCSC foundTAWithMask = CommonOpsWithSemiRing_DSCC.multTransA( + transposed_vector, inputMatrix, prevResult.copy(), semiRing, mask, null, null, null); + + assertMaskedResult(prevResult, foundTA, foundTAWithMask, mask); + } + + @ParameterizedTest + @MethodSource("sparseVectorSource") + public void mult_A_T_B(DMatrixSparseCSC sparseVector, DMatrixSparseCSC prevResult, Mask mask) { + DMatrixSparseCSC transposed_matrix = CommonOps_DSCC.transpose(inputMatrix, null, null); + + DMatrixSparseCSC foundTB = CommonOpsWithSemiRing_DSCC.multTransB( + sparseVector, transposed_matrix, prevResult.copy(), semiRing, null, null, null, null); + DMatrixSparseCSC foundTBWithMask = CommonOpsWithSemiRing_DSCC.multTransB( + sparseVector, transposed_matrix, prevResult.copy(), semiRing, mask, null, null, null); + + assertMaskedResult(prevResult, foundTB, foundTBWithMask, mask); + } + + + @ParameterizedTest + @MethodSource("sparseMatrixSource") + public void add_A_B(DMatrixSparseCSC otherMatrix, DMatrixSparseCSC prevResult, Mask mask) { + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.add( + 1, otherMatrix, 1, inputMatrix, prevResult.copy(), semiRing, null, null, null, null); + DMatrixSparseCSC foundWithMask = CommonOpsWithSemiRing_DSCC.add( + 1, otherMatrix, 1, inputMatrix, prevResult.copy(), semiRing, mask, null, null, null); + + assertMaskedResult(prevResult, found, foundWithMask, mask); + } + + @ParameterizedTest + @MethodSource("sparseMatrixSource") + public void elementWiseMult(DMatrixSparseCSC otherMatrix, DMatrixSparseCSC prevResult, Mask mask) { + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.elementMult( + otherMatrix, inputMatrix, prevResult.copy(), semiRing, null, null, null, null); + DMatrixSparseCSC foundWithMask = CommonOpsWithSemiRing_DSCC.elementMult( + otherMatrix, inputMatrix, prevResult.copy(), semiRing, mask, null, null, null); + + assertMaskedResult(prevResult, found, foundWithMask, mask); + } + + @ParameterizedTest + @MethodSource("sparseVectorSource") + public void apply(DMatrixSparseCSC vector, DMatrixSparseCSC prevResult, Mask mask) { + DBinaryOperator second = (x, y) -> y; + DMatrixSparseCSC result = CommonOps_DSCC.apply(vector, a -> a * 2, prevResult.copy(), null, null); + DMatrixSparseCSC resultWithMask = CommonOps_DSCC.apply(vector, a -> a * 2, prevResult.copy(), mask, second); + + assertMaskedResult(prevResult, result, resultWithMask, mask); + } + + @ParameterizedTest + @ValueSource(strings = {"true", "false"}) + public void reduceRowWise(boolean negatedMask) { + DMatrixSparseCSC matrix = new DMatrixSparseCSC(7, 7); + matrix.set(0, 3, 0.5); + matrix.set(0, 4, 0.6); + matrix.set(4, 0, 0.1); + matrix.set(3, 0, 0.2); + + double[] prevPrimitiveResult = new double[7]; + prevPrimitiveResult[2] = 42; + prevPrimitiveResult[0] = 99; + + DMatrixRMaj prevResult = DMatrixRMaj.wrap(matrix.numRows, 1, prevPrimitiveResult); + Mask mask = DMasks.builder(prevResult).withNegated(negatedMask).withReplace(false).build(); + + DMatrixRMaj result = CommonOps_DSCC.reduceRowWise(matrix, 0, Double::sum, prevResult.copy(), null, null); + DMatrixRMaj resultWithMask = CommonOps_DSCC.reduceRowWise(matrix, 0, Double::sum, prevResult.copy(), mask, null); + + assertMaskedResult(prevResult, result, resultWithMask, mask); + } + + @ParameterizedTest + @ValueSource(strings = {"true", "false"}) + public void reduceColumnWise(boolean negatedMask) { + DMatrixSparseCSC matrix = new DMatrixSparseCSC(7, 7); + matrix.set(0, 3, 0.5); + matrix.set(0, 4, 0.6); + matrix.set(4, 0, 0.1); + matrix.set(3, 0, 0.2); + + double[] prevPrimitiveResult = new double[7]; + prevPrimitiveResult[2] = 42; + prevPrimitiveResult[0] = 99; + + DMatrixRMaj prevResult = DMatrixRMaj.wrap(1, matrix.numCols, prevPrimitiveResult); + Mask mask = DMasks.builder(prevResult).withNegated(negatedMask).withReplace(false).build(); + + DMatrixRMaj result = CommonOps_DSCC.reduceColumnWise(matrix, 0, Double::sum, prevResult.copy(), null, null); + DMatrixRMaj resultWithMask = CommonOps_DSCC.reduceColumnWise(matrix, 0, Double::sum, prevResult.copy(), mask, null); + + assertMaskedResult(prevResult, result, resultWithMask, mask); + } + + + private void assertMaskedResult(DMatrixSparseCSC prevResult, DMatrixSparseCSC found, DMatrixSparseCSC foundWithMask, Mask mask) { + Iterator it = found.createCoordinateIterator(); + // check that existing result were not overwritten + it.forEachRemaining(value -> { + if (mask.isSet(value.row, value.col)) { + assertEquals(found.get(value.row, value.col), foundWithMask.get(value.row, value.col), "Field should have been computed"); + } else if(!mask.replace) { + assertEquals(prevResult.get(value.row, value.col), foundWithMask.get(value.row, value.col), "Field from initial result was overwritten"); + } + }); + + // checking that untouched cells are still present + prevResult.createCoordinateIterator().forEachRemaining(value -> { + if (!mask.isSet(value.row, value.col) && !mask.replace) { + assertEquals(prevResult.get(value.row, value.col), foundWithMask.get(value.row, value.col), "Field from initial result was deleted"); + } + else if (mask.replace) { + assertEquals(found.isAssigned(value.row, value.col), foundWithMask.isAssigned(value.row, value.col)); + } + }); + } + + private void assertMaskedResult(DMatrixD1 prevResult, DMatrixD1 found, DMatrixD1 foundWithMask, Mask mask) { + for (int row = 0; row < found.getNumRows(); row++) { + for (int col = 0; col < found.getNumCols(); col++) { + if (mask.isSet(row, col) || mask.replace) { + assertEquals(found.get(row, col), foundWithMask.get(row, col), "Field should have been computed/Initial result cleared"); + } else { + assertEquals(prevResult.get(row, col), foundWithMask.get(row, col), "Field from initial result was overwritten"); + } + } + } + } + + private void assertMaskedResult(double[] prevResult, double[] found, double[] foundWithMask, PrimitiveDMask mask) { + for (int i = 0; i < found.length; i++) { + if (mask.isSet(i) || mask.replace) { + assertEquals(foundWithMask[i], found[i], "Computation differs"); + } else { + assertEquals(prevResult[i], foundWithMask[i], "Initial result was overwritten"); + } + } + } +} \ No newline at end of file diff --git a/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMatrixMatrixMultWithSemiRing_DSCC.java b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMatrixMatrixMultWithSemiRing_DSCC.java new file mode 100644 index 000000000..3004d634d --- /dev/null +++ b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMatrixMatrixMultWithSemiRing_DSCC.java @@ -0,0 +1,120 @@ +/* + * Copyright (c) 2009-2020, 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.mult; + +import org.ejml.EjmlUnitTests; +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.ops.DSemiRing; +import org.ejml.ops.DSemiRings; +import org.ejml.sparse.csc.CommonOpsWithSemiRing_DSCC; +import org.ejml.sparse.csc.CommonOps_DSCC; +import org.ejml.sparse.csc.RandomMatrices_DSCC; +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.Arguments; +import org.junit.jupiter.params.provider.MethodSource; + +import java.util.Random; +import java.util.stream.Stream; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +public class TestMatrixMatrixMultWithSemiRing_DSCC extends BaseTestMatrixMatrixOpsWithSemiRing_DSCC { + @ParameterizedTest(name = "{0}") + @MethodSource("sparseVectorMatrixMultSources") + public void mult_v_A(String desc, DSemiRing semiRing, double[] expected) { + // graphblas == following outgoing edges of source nodes + DMatrixSparseCSC vector = new DMatrixSparseCSC(1, 7); + vector.set(0, 3, 0.5); + vector.set(0, 5, 0.6); + + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.mult(vector, inputMatrix, null, semiRing, null, null); + + assertEquals(expected[0], found.get(0, 0)); + assertEquals(expected[1], found.get(0, 2)); + } + + @ParameterizedTest(name = "{0}") + @MethodSource("sparseMatrixSources") + public void mult_AT_B(String desc, DMatrixSparseCSC matrix, DMatrixSparseCSC otherMatrix) { + DSemiRing semiRing = DSemiRings.PLUS_TIMES; + + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.multTransA(matrix, otherMatrix, null, semiRing, null, null, null, null); + DMatrixSparseCSC expected = CommonOps_DSCC.multTransA(matrix, otherMatrix, null, null, null); + + EjmlUnitTests.assertEquals(expected, found); + } + + @ParameterizedTest(name = "{0}") + @MethodSource("sparseMatrixSources") + public void mult_A_BT(String desc, DMatrixSparseCSC matrix, DMatrixSparseCSC otherMatrix) { + DSemiRing semiRing = DSemiRings.PLUS_TIMES; + + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.multTransB(matrix, otherMatrix, null, semiRing, null, null, null, null); + DMatrixSparseCSC expected = CommonOps_DSCC.multTransB(matrix, otherMatrix, null, null, null); + + EjmlUnitTests.assertEquals(expected, found); + } + + @ParameterizedTest(name = "{0}") + @MethodSource("sparseMatrixSources") + public void elementMult(String desc, DMatrixSparseCSC matrix, DMatrixSparseCSC otherMatrix) { + DSemiRing semiRing = DSemiRings.PLUS_TIMES; + + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.elementMult(matrix, otherMatrix, null, semiRing, null, null, null, null); + DMatrixSparseCSC expected = CommonOps_DSCC.elementMult(matrix, otherMatrix, null, null, null); + + EjmlUnitTests.assertEquals(expected, found); + } + + @ParameterizedTest(name = "{0}") + @MethodSource("sparseMatrixSources") + public void add(String desc, DMatrixSparseCSC matrix, DMatrixSparseCSC otherMatrix) { + DSemiRing semiRing = DSemiRings.PLUS_TIMES; + + DMatrixSparseCSC found = CommonOpsWithSemiRing_DSCC.add(1, matrix, 1, otherMatrix, null, semiRing, null, null, null, null); + DMatrixSparseCSC expected = CommonOps_DSCC.add(1, matrix, 1, otherMatrix, null, null, null); + + EjmlUnitTests.assertEquals(expected, found); + } + + private static Stream sparseVectorMatrixMultSources() { + return Stream.of( + // expected entries for (0, 0) and (0, 2) + Arguments.of("Plus, Times", DSemiRings.PLUS_TIMES, new double[]{0.1, 0.5}), + Arguments.of("OR, AND", DSemiRings.OR_AND, new double[]{1, 1}), + Arguments.of("MIN, PLUS", DSemiRings.MIN_PLUS, new double[]{0.7, 0.9}), + Arguments.of("MAX, PLUS", DSemiRings.MAX_PLUS, new double[]{0.7, 1.1}), + Arguments.of("MIN, TIMES", DSemiRings.MIN_TIMES, new double[]{0.1, 0.2}), + Arguments.of("MAX, MIN", DSemiRings.MAX_MIN, new double[]{0.2, 0.5}) + ); + } + + private static Stream sparseMatrixSources() { + Random rand = new Random(42); + DMatrixSparseCSC sparseMatrix = RandomMatrices_DSCC.rectangle(10, 10, 15, rand); + DMatrixSparseCSC denseMatrix = RandomMatrices_DSCC.rectangle(10, 10, 90, rand); + return Stream.of( + Arguments.of("Both really sparse", sparseMatrix, sparseMatrix), + Arguments.of("Sparse, dense", sparseMatrix, denseMatrix), + Arguments.of("dense, sparse", denseMatrix, sparseMatrix), + Arguments.of("dense, dense", denseMatrix, denseMatrix) + ); + } +} \ No newline at end of file diff --git a/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMatrixVectorMultWithSemiRing_DSCC.java b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMatrixVectorMultWithSemiRing_DSCC.java new file mode 100644 index 000000000..e4bd1ec95 --- /dev/null +++ b/main/ejml-dsparse/test/org/ejml/sparse/csc/mult/TestMatrixVectorMultWithSemiRing_DSCC.java @@ -0,0 +1,124 @@ +/* + * Copyright (c) 2009-2020, 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.mult; + +import org.ejml.data.DMatrixSparseCSC; +import org.ejml.ops.DSemiRing; +import org.ejml.ops.DSemiRings; +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.Arguments; +import org.junit.jupiter.params.provider.MethodSource; + +import java.util.Arrays; +import java.util.stream.Stream; + +import static org.junit.jupiter.api.Assertions.assertTrue; + +public class TestMatrixVectorMultWithSemiRing_DSCC { + DMatrixSparseCSC inputMatrix; + + @BeforeEach + public void setUp() { + // based on example in http://mit.bme.hu/~szarnyas/grb/graphblas-introduction.pdf + inputMatrix = new DMatrixSparseCSC(7, 7, 12); + inputMatrix.set(0, 1, 1); + inputMatrix.set(0, 3, 1); + inputMatrix.set(1, 4, 1); + inputMatrix.set(1, 6, 1); + inputMatrix.set(2, 5, 1); + inputMatrix.set(3, 0, 0.2); + inputMatrix.set(3, 2, 0.4); + inputMatrix.set(4, 5, 1); + inputMatrix.set(5, 2, 0.5); + inputMatrix.set(6, 2, 1); + inputMatrix.set(6, 3, 1); + inputMatrix.set(6, 4, 1); + + } + + @ParameterizedTest(name = "{0}") + @MethodSource("vectorMatrixMultSources") + public void mult_v_A(String desc, DSemiRing semiRing, double[] expected) { + // graphblas == following outgoing edges of source nodes + double v[] = new double[7]; + Arrays.fill(v, semiRing.add.id); + v[3] = 0.5; + v[5] = 0.6; + + System.out.println("input vector = " + Arrays.toString(v)); + + //input.print(); + + double found[] = new double[7]; + + MatrixVectorMultWithSemiRing_DSCC.mult(v, inputMatrix, found, semiRing); + + System.out.println("found = " + Arrays.toString(found)); + System.out.println("expected = " + Arrays.toString(expected)); + assertTrue(Arrays.equals(found, expected)); + } + + @ParameterizedTest(name = "{0}") + @MethodSource("matrixVectorMultSources") + public void mult_A_v(String desc, DSemiRing semiRing, double[] expected) { + // graphblas == following incoming edges of source nodes + double v[] = new double[7]; + Arrays.fill(v, semiRing.add.id); + v[3] = 0.5; + v[4] = 0.6; + + System.out.println("input vector = " + Arrays.toString(v)); + + //input.print(); + + double found[] = new double[7]; + + MatrixVectorMultWithSemiRing_DSCC.mult(inputMatrix, v, found, semiRing); + + System.out.println("found = " + Arrays.toString(found)); + System.out.println("expected = " + Arrays.toString(expected)); + assertTrue(Arrays.equals(found, expected)); + } + + private static Stream vectorMatrixMultSources() { + return Stream.of( + Arguments.of("Plus, Times", DSemiRings.PLUS_TIMES, new double[]{0.1, 0, 0.5, 0, 0, 0, 0}), + Arguments.of("OR, AND", DSemiRings.OR_AND, new double[]{1, 0, 1, 0, 0, 0, 0}), + Arguments.of("MIN, PLUS", DSemiRings.MIN_PLUS, + new double[]{0.7, Double.MAX_VALUE, 0.9, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE}), +// This only works on sparse input vectors (here resulting in 1.0 instead of Double.MIN_Value as max(1.0, Double.MIN_VALUE) = 1.0) +// Arguments.of("MAX, PLUS", PreDefinedDoubleSemiRings.MAX_PLUS, +// new double[]{0.7, Double.MIN_VALUE, 1.1, Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_VALUE}), + Arguments.of("MIN, TIMES", DSemiRings.MIN_TIMES, + new double[]{0.1, Double.MAX_VALUE, 0.2, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE}), + Arguments.of("MAX, MIN", DSemiRings.MAX_MIN, + new double[]{0.2, Double.MIN_VALUE, 0.5, Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_VALUE}) + ); + } + + private static Stream matrixVectorMultSources() { + return Stream.of( + Arguments.of("Plus, Times", DSemiRings.PLUS_TIMES, new double[]{0.5, 0.6, 0, 0, 0, 0, 1.1}), + Arguments.of("OR, AND", DSemiRings.OR_AND, new double[]{1, 1, 0, 0, 0, 0, 1}), + Arguments.of("MIN, PLUS", DSemiRings.MIN_PLUS, + new double[]{1.5, 1.6, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, 1.5}) + ); + } +} \ No newline at end of file