This is an automated email from the ASF dual-hosted git repository. baunsgaard pushed a commit to branch main in repository https://gitbox.apache.org/repos/asf/systemds.git
commit 52639b2e57683cf8c9ab6ecd26c5a961100f2d79 Author: Jessica Eva Sophie Priebe <[email protected]> AuthorDate: Thu Apr 4 17:17:07 2024 +0200 [SYSTEMDS-3685] DML Integration of FFT and IFFT This commits integrate 4 new builtin functions: FFT FFT_LINEARIZED IFFT IFFT_LINEARIZED The functions implement fast fourier transformations and inverces. The linearized functions are for performing the transformations equivalently on each row in a matrix, while the normal ones are able to perform 2-d ffts if given a matrix. The return of a FFC is a normal and complex matrix pair. LDE 23/24 project Co-authored-by: Mufan Wang <[email protected]> Co-authored-by: Frederic Caspar Zoepffel <[email protected]> Co-authored-by: Jessica Eva Sophie Priebe <[email protected]> Closes #1995 --- .../java/org/apache/sysds/common/Builtins.java | 4 + .../java/org/apache/sysds/hops/FunctionOp.java | 36 ++ .../sysds/parser/BuiltinFunctionExpression.java | 141 ++++++ .../org/apache/sysds/parser/DMLTranslator.java | 4 + .../runtime/instructions/CPInstructionParser.java | 10 +- .../runtime/instructions/cp/CPInstruction.java | 2 +- .../cp/MultiReturnBuiltinCPInstruction.java | 34 ++ ...tiReturnComplexMatrixBuiltinCPInstruction.java} | 117 +++-- .../sysds/runtime/matrix/data/LibCommonsMath.java | 167 ++++++- .../runtime/matrix/data/LibMatrixFourier.java | 479 +++++++++++++++++++++ .../sysds/test/component/matrix/FourierTest.java | 344 +++++++++++++++ 11 files changed, 1287 insertions(+), 51 deletions(-) diff --git a/src/main/java/org/apache/sysds/common/Builtins.java b/src/main/java/org/apache/sysds/common/Builtins.java index 4d0e13791f..7e83984e47 100644 --- a/src/main/java/org/apache/sysds/common/Builtins.java +++ b/src/main/java/org/apache/sysds/common/Builtins.java @@ -133,6 +133,8 @@ public enum Builtins { FIT_PIPELINE("fit_pipeline", true), FIX_INVALID_LENGTHS("fixInvalidLengths", true), FIX_INVALID_LENGTHS_APPLY("fixInvalidLengthsApply", true), + FFT("fft", false, ReturnType.MULTI_RETURN), + FFT_LINEARIZED("fft_linearized", false, ReturnType.MULTI_RETURN), FF_TRAIN("ffTrain", true), FF_PREDICT("ffPredict", true), FLOOR("floor", false), @@ -154,6 +156,8 @@ public enum Builtins { HOSPITAL_RESIDENCY_MATCH("hospitalResidencyMatch", true), HYPERBAND("hyperband", true), IFELSE("ifelse", false), + IFFT("ifft", false, ReturnType.MULTI_RETURN), + IFFT_LINEARIZED("ifft_linearized", false, ReturnType.MULTI_RETURN), IMG_MIRROR("img_mirror", true), IMG_MIRROR_LINEARIZED("img_mirror_linearized", true), IMG_BRIGHTNESS("img_brightness", true), diff --git a/src/main/java/org/apache/sysds/hops/FunctionOp.java b/src/main/java/org/apache/sysds/hops/FunctionOp.java index 28cd6eeafb..ffc12c30ee 100644 --- a/src/main/java/org/apache/sysds/hops/FunctionOp.java +++ b/src/main/java/org/apache/sysds/hops/FunctionOp.java @@ -201,6 +201,26 @@ public class FunctionOp extends Hop long outputValues = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(1).getDim1(), 1, 1.0); return outputVectors+outputValues; } + else if ( getFunctionName().equalsIgnoreCase("fft") ) { + long outputRe = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(0).getDim1(), getOutputs().get(0).getDim2(), 1.0); + long outputIm = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(1).getDim1(), getOutputs().get(1).getDim2(), 1.0); + return outputRe+outputIm; + } + else if ( getFunctionName().equalsIgnoreCase("ifft") ) { + long outputRe = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(0).getDim1(), getOutputs().get(0).getDim2(), 1.0); + long outputIm = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(1).getDim1(), getOutputs().get(1).getDim2(), 1.0); + return outputRe+outputIm; + } + else if ( getFunctionName().equalsIgnoreCase("fft_linearized") ) { + long outputRe = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(0).getDim1(), getOutputs().get(0).getDim2(), 1.0); + long outputIm = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(1).getDim1(), getOutputs().get(1).getDim2(), 1.0); + return outputRe+outputIm; + } + else if ( getFunctionName().equalsIgnoreCase("ifft_linearized") ) { + long outputRe = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(0).getDim1(), getOutputs().get(0).getDim2(), 1.0); + long outputIm = OptimizerUtils.estimateSizeExactSparsity(getOutputs().get(1).getDim1(), getOutputs().get(1).getDim2(), 1.0); + return outputRe+outputIm; + } else if ( getFunctionName().equalsIgnoreCase("lstm") || getFunctionName().equalsIgnoreCase("lstm_backward") ) { // TODO: To allow for initial version to always run on the GPU return 0; @@ -250,6 +270,22 @@ public class FunctionOp extends Hop return OptimizerUtils.estimateSizeExactSparsity(getInput().get(0).getDim1(), getInput().get(0).getDim2(), 1.0) + 3*OptimizerUtils.estimateSizeExactSparsity(getInput().get(0).getDim1(), 1, 1.0); } + else if ( getFunctionName().equalsIgnoreCase("fft") ) { + // 2 matrices of size same as the input + return 2*OptimizerUtils.estimateSizeExactSparsity(getInput().get(0).getDim1(), getInput().get(0).getDim2(), 1.0); + } + else if ( getFunctionName().equalsIgnoreCase("ifft") ) { + // 2 matrices of size same as the input + return 2*OptimizerUtils.estimateSizeExactSparsity(getInput().get(0).getDim1(), getInput().get(0).getDim2(), 1.0); + } + else if ( getFunctionName().equalsIgnoreCase("fft_linearized") ) { + // 2 matrices of size same as the input + return 2*OptimizerUtils.estimateSizeExactSparsity(getInput().get(0).getDim1(), getInput().get(0).getDim2(), 1.0); + } + else if ( getFunctionName().equalsIgnoreCase("ifft_linearized") ) { + // 2 matrices of size same as the input + return 2*OptimizerUtils.estimateSizeExactSparsity(getInput().get(0).getDim1(), getInput().get(0).getDim2(), 1.0); + } else if (getFunctionName().equalsIgnoreCase("batch_norm2d") || getFunctionName().equalsIgnoreCase("batch_norm2d_backward") || getFunctionName().equalsIgnoreCase("batch_norm2d_train") || getFunctionName().equalsIgnoreCase("batch_norm2d_test")) { return 0; diff --git a/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java b/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java index b5b27682d5..5e86a2fd8e 100644 --- a/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java +++ b/src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java @@ -380,6 +380,143 @@ public class BuiltinFunctionExpression extends DataIdentifier { break; } + case FFT: { + checkNumParameters(1); + checkMatrixParam(getFirstExpr()); + + // setup output properties + DataIdentifier fftOut1 = (DataIdentifier) getOutputs()[0]; + DataIdentifier fftOut2 = (DataIdentifier) getOutputs()[1]; + + // Output1 - FFT Values + fftOut1.setDataType(DataType.MATRIX); + fftOut1.setValueType(ValueType.FP64); + fftOut1.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + fftOut1.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + // Output2 - FFT Vectors + fftOut2.setDataType(DataType.MATRIX); + fftOut2.setValueType(ValueType.FP64); + fftOut2.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + fftOut2.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + break; + + } + case IFFT: { + Expression expressionTwo = getSecondExpr(); + checkNumParameters(getSecondExpr() != null ? 2 : 1); + checkMatrixParam(getFirstExpr()); + if (expressionTwo != null) + checkMatrixParam(getSecondExpr()); + + // setup output properties + DataIdentifier ifftOut1 = (DataIdentifier) getOutputs()[0]; + DataIdentifier ifftOut2 = (DataIdentifier) getOutputs()[1]; + + // Output1 - ifft Values + ifftOut1.setDataType(DataType.MATRIX); + ifftOut1.setValueType(ValueType.FP64); + ifftOut1.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + ifftOut1.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + // Output2 - ifft Vectors + ifftOut2.setDataType(DataType.MATRIX); + ifftOut2.setValueType(ValueType.FP64); + ifftOut2.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + ifftOut2.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + break; + } + case FFT_LINEARIZED: { + + Expression expressionOne = getFirstExpr(); + Expression expressionTwo = getSecondExpr(); + + if (expressionOne == null) { + raiseValidateError("The first argument to " + _opcode + " cannot be null.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + + else if (expressionOne.getOutput() == null || expressionOne.getOutput().getDim1() == 0 || expressionOne.getOutput().getDim2() == 0) { + raiseValidateError("The first argument to " + _opcode + " cannot be an empty matrix.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + + else if (expressionTwo != null) { + raiseValidateError("Too many arguments. This FFT_LINEARIZED implementation is only defined for real inputs.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + + else if (!isPowerOfTwo(expressionOne.getOutput().getDim2())) { + raiseValidateError("This FFT_LINEARIZED implementation is only defined for matrices with columns that are powers of 2.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + + checkNumParameters(1); + checkMatrixParam(expressionOne); + + DataIdentifier fftOut1 = (DataIdentifier) getOutputs()[0]; + DataIdentifier fftOut2 = (DataIdentifier) getOutputs()[1]; + + fftOut1.setDataType(DataType.MATRIX); + fftOut1.setValueType(ValueType.FP64); + fftOut1.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + fftOut1.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + fftOut2.setDataType(DataType.MATRIX); + fftOut2.setValueType(ValueType.FP64); + fftOut2.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + fftOut2.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + break; + + } + case IFFT_LINEARIZED: { + Expression expressionTwo = getSecondExpr(); + Expression expressionOne = getFirstExpr(); + + if (expressionOne == null) { + raiseValidateError("The first argument to " + _opcode + " cannot be null.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + + else if (expressionOne.getOutput() == null || expressionOne.getOutput().getDim1() == 0 || expressionOne.getOutput().getDim2() == 0) { + raiseValidateError("The first argument to " + _opcode + " cannot be an empty matrix.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + + else if (expressionTwo != null){ + if(expressionTwo.getOutput() == null || expressionTwo.getOutput().getDim1() == 0 || expressionTwo.getOutput().getDim2() == 0) { + raiseValidateError("The second argument to " + _opcode + " cannot be an empty matrix. Provide either only a real matrix or a filled real and imaginary one.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + } + + checkNumParameters(expressionTwo != null ? 2 : 1); + checkMatrixParam(expressionOne); + if(expressionTwo != null && expressionOne != null){ + checkMatrixParam(expressionTwo); + if(expressionOne.getOutput().getDim1() != expressionTwo.getOutput().getDim1() || expressionOne.getOutput().getDim2() != expressionTwo.getOutput().getDim2()) + raiseValidateError("The real and imaginary part of the provided matrix are of different dimensions.", false); + else if (!isPowerOfTwo(expressionTwo.getOutput().getDim2())) { + raiseValidateError("This IFFT_LINEARIZED implementation is only defined for matrices with columns that are powers of 2.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + } + else if(expressionOne != null){ + if (!isPowerOfTwo(expressionOne.getOutput().getDim2())) { + raiseValidateError("This IFFT_LINEARIZED implementation is only defined for matrices with columns that are powers of 2.", false, LanguageErrorCodes.INVALID_PARAMETERS); + } + } + + DataIdentifier ifftOut1 = (DataIdentifier) getOutputs()[0]; + DataIdentifier ifftOut2 = (DataIdentifier) getOutputs()[1]; + + ifftOut1.setDataType(DataType.MATRIX); + ifftOut1.setValueType(ValueType.FP64); + ifftOut1.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + ifftOut1.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + ifftOut2.setDataType(DataType.MATRIX); + ifftOut2.setValueType(ValueType.FP64); + ifftOut2.setDimensions(getFirstExpr().getOutput().getDim1(), getFirstExpr().getOutput().getDim2()); + ifftOut2.setBlocksize(getFirstExpr().getOutput().getBlocksize()); + + break; + } case REMOVE: { checkNumParameters(2); checkListParam(getFirstExpr()); @@ -464,6 +601,10 @@ public class BuiltinFunctionExpression extends DataIdentifier { raiseValidateError("Unknown Builtin Function opcode: " + _opcode, false); } } + + private static boolean isPowerOfTwo(long n) { + return (n > 0) && ((n & (n - 1)) == 0); + } private static void setDimensions(DataIdentifier out, Expression exp) { out.setDataType(DataType.MATRIX); diff --git a/src/main/java/org/apache/sysds/parser/DMLTranslator.java b/src/main/java/org/apache/sysds/parser/DMLTranslator.java index d0232d63bd..97d5523961 100644 --- a/src/main/java/org/apache/sysds/parser/DMLTranslator.java +++ b/src/main/java/org/apache/sysds/parser/DMLTranslator.java @@ -2236,6 +2236,10 @@ public class DMLTranslator case QR: case LU: case EIGEN: + case FFT: + case IFFT: + case FFT_LINEARIZED: + case IFFT_LINEARIZED: case LSTM: case LSTM_BACKWARD: case BATCH_NORM2D: diff --git a/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java b/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java index 52f6257bf9..b78c1e2b49 100644 --- a/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java +++ b/src/main/java/org/apache/sysds/runtime/instructions/CPInstructionParser.java @@ -57,6 +57,7 @@ import org.apache.sysds.runtime.instructions.cp.LocalCPInstruction; import org.apache.sysds.runtime.instructions.cp.MMChainCPInstruction; import org.apache.sysds.runtime.instructions.cp.MMTSJCPInstruction; import org.apache.sysds.runtime.instructions.cp.MultiReturnBuiltinCPInstruction; +import org.apache.sysds.runtime.instructions.cp.MultiReturnComplexMatrixBuiltinCPInstruction; import org.apache.sysds.runtime.instructions.cp.MultiReturnParameterizedBuiltinCPInstruction; import org.apache.sysds.runtime.instructions.cp.PMMJCPInstruction; import org.apache.sysds.runtime.instructions.cp.ParameterizedBuiltinCPInstruction; @@ -291,7 +292,7 @@ public class CPInstructionParser extends InstructionParser { String2CPInstructionType.put( "batch_norm2d_backward", CPType.Dnn); String2CPInstructionType.put( "lstm" , CPType.Dnn); String2CPInstructionType.put( "lstm_backward" , CPType.Dnn); - + // Quaternary instruction opcodes String2CPInstructionType.put( "wsloss" , CPType.Quaternary); String2CPInstructionType.put( "wsigmoid", CPType.Quaternary); @@ -333,6 +334,10 @@ public class CPInstructionParser extends InstructionParser { String2CPInstructionType.put( "qr", CPType.MultiReturnBuiltin); String2CPInstructionType.put( "lu", CPType.MultiReturnBuiltin); String2CPInstructionType.put( "eigen", CPType.MultiReturnBuiltin); + String2CPInstructionType.put( "fft", CPType.MultiReturnBuiltin); + String2CPInstructionType.put( "ifft", CPType.MultiReturnComplexMatrixBuiltin); + String2CPInstructionType.put( "fft_linearized", CPType.MultiReturnBuiltin); + String2CPInstructionType.put( "ifft_linearized", CPType.MultiReturnComplexMatrixBuiltin); String2CPInstructionType.put( "svd", CPType.MultiReturnBuiltin); String2CPInstructionType.put( "partition", CPType.Partition); @@ -424,6 +429,9 @@ public class CPInstructionParser extends InstructionParser { case MultiReturnParameterizedBuiltin: return MultiReturnParameterizedBuiltinCPInstruction.parseInstruction(str); + + case MultiReturnComplexMatrixBuiltin: + return MultiReturnComplexMatrixBuiltinCPInstruction.parseInstruction(str); case MultiReturnBuiltin: return MultiReturnBuiltinCPInstruction.parseInstruction(str); diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java index aec67ce844..f8527276a7 100644 --- a/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java +++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/CPInstruction.java @@ -40,7 +40,7 @@ public abstract class CPInstruction extends Instruction { public enum CPType { AggregateUnary, AggregateBinary, AggregateTernary, Unary, Binary, Ternary, Quaternary, BuiltinNary, Ctable, - MultiReturnParameterizedBuiltin, ParameterizedBuiltin, MultiReturnBuiltin, + MultiReturnParameterizedBuiltin, ParameterizedBuiltin, MultiReturnBuiltin, MultiReturnComplexMatrixBuiltin, Builtin, Reorg, Variable, FCall, Append, Rand, QSort, QPick, Local, MatrixIndexing, MMTSJ, PMMJ, MMChain, Reshape, Partition, Compression, DeCompression, SpoofFused, StringInit, CentralMoment, Covariance, UaggOuterChain, Dnn, Sql, Prefetch, Broadcast, TrigRemote, diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java index 0d28a62f4b..4719948ec0 100644 --- a/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java +++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java @@ -89,6 +89,40 @@ public class MultiReturnBuiltinCPInstruction extends ComputationCPInstruction { return new MultiReturnBuiltinCPInstruction(null, in1, outputs, opcode, str); + } + else if(parts.length == 4 && opcode.equalsIgnoreCase("fft")) { + // one input and two outputs + CPOperand in1 = new CPOperand(parts[1]); + outputs.add(new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX)); + + return new MultiReturnBuiltinCPInstruction(null, in1, outputs, opcode, str); + + } + else if(parts.length == 3 && opcode.equalsIgnoreCase("fft")) { + // one input and two outputs + outputs.add(new CPOperand(parts[1], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX)); + + return new MultiReturnBuiltinCPInstruction(null, null, outputs, opcode, str); + + } + else if(parts.length == 4 && opcode.equalsIgnoreCase("fft_linearized")) { + // one input and two outputs + CPOperand in1 = new CPOperand(parts[1]); + outputs.add(new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX)); + + return new MultiReturnBuiltinCPInstruction(null, in1, outputs, opcode, str); + + } + else if(parts.length == 3 && opcode.equalsIgnoreCase("fft_linearized")) { + // one input and two outputs + outputs.add(new CPOperand(parts[1], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX)); + + return new MultiReturnBuiltinCPInstruction(null, null, outputs, opcode, str); + } else if ( opcode.equalsIgnoreCase("svd") ) { CPOperand in1 = new CPOperand(parts[1]); diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnComplexMatrixBuiltinCPInstruction.java similarity index 51% copy from src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java copy to src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnComplexMatrixBuiltinCPInstruction.java index 0d28a62f4b..f9c92d4858 100644 --- a/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java +++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnComplexMatrixBuiltinCPInstruction.java @@ -34,111 +34,134 @@ import org.apache.sysds.runtime.matrix.data.LibCommonsMath; import org.apache.sysds.runtime.matrix.data.MatrixBlock; import org.apache.sysds.runtime.matrix.operators.Operator; -public class MultiReturnBuiltinCPInstruction extends ComputationCPInstruction { +public class MultiReturnComplexMatrixBuiltinCPInstruction extends ComputationCPInstruction { protected ArrayList<CPOperand> _outputs; - private MultiReturnBuiltinCPInstruction(Operator op, CPOperand input1, ArrayList<CPOperand> outputs, String opcode, - String istr) { + private MultiReturnComplexMatrixBuiltinCPInstruction(Operator op, CPOperand input1, CPOperand input2, + ArrayList<CPOperand> outputs, String opcode, String istr) { + super(CPType.MultiReturnBuiltin, op, input1, input2, outputs.get(0), opcode, istr); + _outputs = outputs; + } + + private MultiReturnComplexMatrixBuiltinCPInstruction(Operator op, CPOperand input1, ArrayList<CPOperand> outputs, + String opcode, String istr) { super(CPType.MultiReturnBuiltin, op, input1, null, outputs.get(0), opcode, istr); _outputs = outputs; } - + public CPOperand getOutput(int i) { return _outputs.get(i); } - public List<CPOperand> getOutputs(){ + public List<CPOperand> getOutputs() { return _outputs; } - public String[] getOutputNames(){ + public String[] getOutputNames() { return _outputs.parallelStream().map(output -> output.getName()).toArray(String[]::new); } - - public static MultiReturnBuiltinCPInstruction parseInstruction ( String str ) { + + public static MultiReturnComplexMatrixBuiltinCPInstruction parseInstruction(String str) { String[] parts = InstructionUtils.getInstructionPartsWithValueType(str); ArrayList<CPOperand> outputs = new ArrayList<>(); // first part is always the opcode String opcode = parts[0]; - - if ( opcode.equalsIgnoreCase("qr") ) { - // one input and two ouputs + + if(parts.length == 5 && opcode.equalsIgnoreCase("ifft")) { + // one input and two outputs CPOperand in1 = new CPOperand(parts[1]); - outputs.add ( new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX) ); - outputs.add ( new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX) ); - - return new MultiReturnBuiltinCPInstruction(null, in1, outputs, opcode, str); + CPOperand in2 = new CPOperand(parts[2]); + outputs.add(new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[4], ValueType.FP64, DataType.MATRIX)); + + return new MultiReturnComplexMatrixBuiltinCPInstruction(null, in1, in2, outputs, opcode, str); } - else if ( opcode.equalsIgnoreCase("lu") ) { + else if(parts.length == 4 && opcode.equalsIgnoreCase("ifft")) { + // one input and two outputs CPOperand in1 = new CPOperand(parts[1]); - - // one input and three outputs - outputs.add ( new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX) ); - outputs.add ( new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX) ); - outputs.add ( new CPOperand(parts[4], ValueType.FP64, DataType.MATRIX) ); - - return new MultiReturnBuiltinCPInstruction(null, in1, outputs, opcode, str); - + outputs.add(new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX)); + + return new MultiReturnComplexMatrixBuiltinCPInstruction(null, in1, outputs, opcode, str); } - else if ( opcode.equalsIgnoreCase("eigen") ) { + else if(parts.length == 5 && opcode.equalsIgnoreCase("ifft_linearized")) { // one input and two outputs CPOperand in1 = new CPOperand(parts[1]); - outputs.add ( new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX) ); - outputs.add ( new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX) ); - - return new MultiReturnBuiltinCPInstruction(null, in1, outputs, opcode, str); - + CPOperand in2 = new CPOperand(parts[2]); + outputs.add(new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[4], ValueType.FP64, DataType.MATRIX)); + + return new MultiReturnComplexMatrixBuiltinCPInstruction(null, in1, in2, outputs, opcode, str); } - else if ( opcode.equalsIgnoreCase("svd") ) { + else if(parts.length == 4 && opcode.equalsIgnoreCase("ifft_linearized")) { + // one input and two outputs CPOperand in1 = new CPOperand(parts[1]); + outputs.add(new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX)); + outputs.add(new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX)); - // one input and three outputs - outputs.add ( new CPOperand(parts[2], ValueType.FP64, DataType.MATRIX) ); - outputs.add ( new CPOperand(parts[3], ValueType.FP64, DataType.MATRIX) ); - outputs.add ( new CPOperand(parts[4], ValueType.FP64, DataType.MATRIX) ); - - return new MultiReturnBuiltinCPInstruction(null, in1, outputs, opcode, str); - + return new MultiReturnComplexMatrixBuiltinCPInstruction(null, in1, outputs, opcode, str); } else { throw new DMLRuntimeException("Invalid opcode in MultiReturnBuiltin instruction: " + opcode); } } - + public int getNumOutputs() { return _outputs.size(); } - @Override + @Override public void processInstruction(ExecutionContext ec) { + if(input2 == null) + processOneInputInstruction(ec); + else + processTwoInputInstruction(ec); + } + + private void processOneInputInstruction(ExecutionContext ec) { if(!LibCommonsMath.isSupportedMultiReturnOperation(getOpcode())) throw new DMLRuntimeException("Invalid opcode in MultiReturnBuiltin instruction: " + getOpcode()); - + MatrixBlock in = ec.getMatrixInput(input1.getName()); MatrixBlock[] out = LibCommonsMath.multiReturnOperations(in, getOpcode()); + ec.releaseMatrixInput(input1.getName()); - for(int i=0; i < _outputs.size(); i++) { + + for(int i = 0; i < _outputs.size(); i++) { + ec.setMatrixOutput(_outputs.get(i).getName(), out[i]); + } + } + + private void processTwoInputInstruction(ExecutionContext ec) { + if(!LibCommonsMath.isSupportedMultiReturnOperation(getOpcode())) + throw new DMLRuntimeException("Invalid opcode in MultiReturnBuiltin instruction: " + getOpcode()); + + MatrixBlock in1 = ec.getMatrixInput(input1.getName()); + MatrixBlock in2 = ec.getMatrixInput(input2.getName()); + MatrixBlock[] out = LibCommonsMath.multiReturnOperations(in1, in2, getOpcode()); + ec.releaseMatrixInput(input1.getName(), input2.getName()); + + for(int i = 0; i < _outputs.size(); i++) { ec.setMatrixOutput(_outputs.get(i).getName(), out[i]); } } - + @Override public boolean hasSingleLineage() { return false; } - @Override @SuppressWarnings("unchecked") public Pair<String, LineageItem>[] getLineageItems(ExecutionContext ec) { LineageItem[] inputLineage = LineageItemUtils.getLineage(ec, input1, input2, input3); - final Pair<String,LineageItem>[] ret = new Pair[_outputs.size()]; - for(int i = 0; i < _outputs.size(); i++){ + final Pair<String, LineageItem>[] ret = new Pair[_outputs.size()]; + for(int i = 0; i < _outputs.size(); i++) { CPOperand out = _outputs.get(i); ret[i] = Pair.of(out.getName(), new LineageItem(getOpcode(), inputLineage)); } - return ret; + return ret; } } diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java index cc399bac79..dc42b15a96 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java @@ -49,6 +49,11 @@ import org.apache.sysds.runtime.matrix.operators.AggregateBinaryOperator; import org.apache.sysds.runtime.matrix.operators.BinaryOperator; import org.apache.sysds.runtime.util.DataConverter; +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft; +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft; +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft_linearized; +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft_linearized; + /** * Library for matrix operations that need invocation of * Apache Commons Math library. @@ -71,7 +76,19 @@ public class LibCommonsMath } public static boolean isSupportedMultiReturnOperation( String opcode ) { - return ( opcode.equals("qr") || opcode.equals("lu") || opcode.equals("eigen") || opcode.equals("svd") ); + + switch (opcode) { + case "qr": + case "lu": + case "eigen": + case "fft": + case "ifft": + case "fft_linearized": + case "ifft_linearized": + case "svd": return true; + default: return false; + } + } public static boolean isSupportedMatrixMatrixOperation( String opcode ) { @@ -91,6 +108,10 @@ public class LibCommonsMath return multiReturnOperations(in, opcode, 1, 1); } + public static MatrixBlock[] multiReturnOperations(MatrixBlock in1, MatrixBlock in2, String opcode) { + return multiReturnOperations(in1, in2, opcode, 1, 1); + } + public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode, int threads, int num_iterations, double tol) { if(opcode.equals("eigen_qr")) return computeEigenQR(in, num_iterations, tol, threads); @@ -113,9 +134,31 @@ public class LibCommonsMath return computeEigenQR(in, threads); else if (opcode.equals("svd")) return computeSvd(in); + else if (opcode.equals("fft")) + return computeFFT(in, threads); + else if (opcode.equals("ifft")) + return computeIFFT(in, threads); + else if (opcode.equals("fft_linearized")) + return computeFFT_LINEARIZED(in, threads); + else if (opcode.equals("ifft_linearized")) + return computeIFFT_LINEARIZED(in, threads); return null; } - + + public static MatrixBlock[] multiReturnOperations(MatrixBlock in1, MatrixBlock in2, String opcode, int threads, + long seed) { + + switch (opcode) { + case "ifft": + return computeIFFT(in1, in2, threads); + case "ifft_linearized": + return computeIFFT_LINEARIZED(in1, in2, threads); + default: + return null; + } + + } + public static MatrixBlock matrixMatrixOperations(MatrixBlock in1, MatrixBlock in2, String opcode) { if(opcode.equals("solve")) { if (in1.getNumRows() != in1.getNumColumns()) @@ -252,6 +295,126 @@ public class LibCommonsMath DataConverter.convertToArray2DRowRealMatrix(in2)); } + /** + * Function to perform FFT on a given matrix. + * + * @param re matrix object + * @param threads number of threads + * @return array of matrix blocks + */ + private static MatrixBlock[] computeFFT(MatrixBlock re, int threads) { + if(re == null) + throw new DMLRuntimeException("Invalid empty block"); + if(re.isEmptyBlock(false)) { + // Return the original matrix as the result + // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0. + return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)}; + } + // run fft + re.sparseToDense(); + return fft(re, threads); + } + + /** + * Function to perform IFFT on a given matrix. + * + * @param re matrix object + * @param im matrix object + * @param threads number of threads + * @return array of matrix blocks + */ + private static MatrixBlock[] computeIFFT(MatrixBlock re, MatrixBlock im, int threads) { + if(re == null) + throw new DMLRuntimeException("Invalid empty block"); + + // run ifft + if(im != null && !im.isEmptyBlock(false)) { + re.sparseToDense(); + im.sparseToDense(); + return ifft(re, im, threads); + } + else { + if(re.isEmptyBlock(false)) { + // Return the original matrix as the result + // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0. + return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)}; + } + re.sparseToDense(); + return ifft(re, threads); + } + } + + /** + * Function to perform IFFT on a given matrix. + * + * @param re matrix object + * @param threads number of threads + * @return array of matrix blocks + */ + private static MatrixBlock[] computeIFFT(MatrixBlock re, int threads) { + return computeIFFT(re, null, threads); + } + + /** + * Function to perform FFT_LINEARIZED on a given matrix. + * + * @param re matrix object + * @param threads number of threads + * @return array of matrix blocks + */ + private static MatrixBlock[] computeFFT_LINEARIZED(MatrixBlock re, int threads) { + if(re == null) + throw new DMLRuntimeException("Invalid empty block"); + if(re.isEmptyBlock(false)) { + // Return the original matrix as the result + // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0. + return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)}; + } + // run fft + re.sparseToDense(); + return fft_linearized(re, threads); + } + + /** + * Function to perform IFFT_LINEARIZED on a given matrix. + * + * @param re matrix object + * @param im matrix object + * @param threads number of threads + * @return array of matrix blocks + */ + private static MatrixBlock[] computeIFFT_LINEARIZED(MatrixBlock re, MatrixBlock im, int threads) { + if(re == null) + throw new DMLRuntimeException("Invalid empty block"); + + // run ifft + if(im != null && !im.isEmptyBlock(false)) { + re.sparseToDense(); + im.sparseToDense(); + return ifft_linearized(re, im, threads); + } + else { + if(re.isEmptyBlock(false)) { + // Return the original matrix as the result + // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0. + return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)}; + } + re.sparseToDense(); + return ifft_linearized(re, threads); + } + } + + /** + * Function to perform IFFT_LINEARIZED on a given matrix + * + * @param re matrix object + * @param threads number of threads + * @return array of matrix blocks + */ + private static MatrixBlock[] computeIFFT_LINEARIZED(MatrixBlock re, int threads) { + return computeIFFT_LINEARIZED(re, null, threads); + } + /** * Performs Singular Value Decomposition. Calls Apache Commons Math SVD. * X = U * Sigma * Vt, where X is the input matrix, diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java new file mode 100644 index 0000000000..f03cb1663f --- /dev/null +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java @@ -0,0 +1,479 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you 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.apache.sysds.runtime.matrix.data; + +import org.apache.commons.math3.util.FastMath; + +import org.apache.sysds.runtime.util.CommonThreadPool; + +import java.lang.ref.SoftReference; +import java.util.List; +import java.util.HashMap; +import java.util.ArrayList; +import java.util.concurrent.Future; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.ExecutionException; + +public class LibMatrixFourier { + + static SoftReference<HashMap<Double, Double>> sinCacheRef = new SoftReference<>(new HashMap<>()); + static SoftReference<HashMap<Double, Double>> cosCacheRef = new SoftReference<>(new HashMap<>()); + + /** + * Function to perform FFT for two given matrices. The first one represents the real values and the second one the + * imaginary values. The output also contains one matrix for the real and one for the imaginary values. + * + * @param re matrix object representing the real values + * @param im matrix object representing the imaginary values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im, int threads) { + + int rows = re.getNumRows(); + int cols = re.getNumColumns(); + if(!isPowerOfTwo(rows) || !isPowerOfTwo(cols)) + throw new RuntimeException("false dimensions"); + + fft(re.getDenseBlockValues(), im.getDenseBlockValues(), rows, cols, threads, true); + + return new MatrixBlock[] {re, im}; + } + + /** + * Function to perform IFFT for two given matrices. The first one represents the real values and the second one the + * imaginary values. The output also contains one matrix for the real and one for the imaginary values. + * + * @param re matrix object representing the real values + * @param im matrix object representing the imaginary values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im, int threads) { + + int rows = re.getNumRows(); + int cols = re.getNumColumns(); + if(!isPowerOfTwo(rows) || !isPowerOfTwo(cols)) + throw new RuntimeException("false dimensions"); + + ifft(re.getDenseBlockValues(), im.getDenseBlockValues(), rows, cols, threads, true); + + return new MatrixBlock[] {re, im}; + } + + /** + * Function to perform FFT for each row of two given matrices. The first one represents the real values and the + * second one the imaginary values. The output also contains one matrix for the real and one for the imaginary + * values. + * + * @param re matrix object representing the real values + * @param im matrix object representing the imaginary values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] fft_linearized(MatrixBlock re, MatrixBlock im, int threads) { + + int rows = re.getNumRows(); + int cols = re.getNumColumns(); + if(!isPowerOfTwo(cols)) + throw new RuntimeException("false dimensions"); + + fft(re.getDenseBlockValues(), im.getDenseBlockValues(), rows, cols, threads, false); + + return new MatrixBlock[] {re, im}; + } + + /** + * Function to perform IFFT for each row of two given matrices. The first one represents the real values and the + * second one the imaginary values. The output also contains one matrix for the real and one for the imaginary + * values. + * + * @param re matrix object representing the real values + * @param im matrix object representing the imaginary values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] ifft_linearized(MatrixBlock re, MatrixBlock im, int threads) { + + int rows = re.getNumRows(); + int cols = re.getNumColumns(); + if(!isPowerOfTwo(cols)) + throw new RuntimeException("false dimensions"); + + ifft(re.getDenseBlockValues(), im.getDenseBlockValues(), rows, cols, threads, false); + + return new MatrixBlock[] {re, im}; + } + + /** + * Function to perform FFT for two given double arrays. The first one represents the real values and the second one + * the imaginary values. Both arrays get updated and contain the result. + * + * @param re array representing the real values + * @param im array representing the imaginary values + * @param rows number of rows + * @param cols number of rows + * @param threads number of threads + * @param inclColCalc if true, fft is also calculated for each column, otherwise only for each row + */ + public static void fft(double[] re, double[] im, int rows, int cols, int threads, boolean inclColCalc) { + + double[] re_inter = new double[rows * cols]; + double[] im_inter = new double[rows * cols]; + + final ExecutorService pool = CommonThreadPool.get(threads); + final List<Future<?>> tasks = new ArrayList<>(); + + try { + final int rBlz = Math.max(rows / threads, 32); + final int cBlz = Math.max(cols / threads, 32); + + for(int i = 0; i < rows; i += rBlz) { + final int start = i; + final int end = Math.min(i + rBlz, rows); + tasks.add(pool.submit(() -> { + for(int j = start; j < end; j++) { + fft_one_dim(re, im, re_inter, im_inter, j * cols, (j + 1) * cols, cols, 1); + } + })); + } + + for(Future<?> f : tasks) + f.get(); + + tasks.clear(); + if(inclColCalc && rows > 1) { + for(int j = 0; j < cols; j += cBlz) { + final int start = j; + final int end = Math.min(j + cBlz, cols); + tasks.add(pool.submit(() -> { + for(int i = start; i < end; i++) { + fft_one_dim(re, im, re_inter, im_inter, i, i + rows * cols, rows, cols); + } + })); + } + + for(Future<?> f : tasks) + f.get(); + } + } + catch(InterruptedException | ExecutionException e) { + throw new RuntimeException(e); + } + finally { + pool.shutdown(); + } + + } + + /** + * Function to perform IFFT for two given double arrays. The first one represents the real values and the second one + * the imaginary values. Both arrays get updated and contain the result. + * + * @param re array representing the real values + * @param im array representing the imaginary values + * @param rows number of rows + * @param cols number of rows + * @param threads number of threads + * @param inclColCalc if true, fft is also calculated for each column, otherwise only for each row + */ + public static void ifft(double[] re, double[] im, int rows, int cols, int threads, boolean inclColCalc) { + + double[] re_inter = new double[rows * cols]; + double[] im_inter = new double[rows * cols]; + + ExecutorService pool = CommonThreadPool.get(threads); + final List<Future<?>> tasks = new ArrayList<>(); + + try { + final int rBlz = Math.max(rows / threads, 32); + final int cBlz = Math.max(cols / threads, 32); + + if(inclColCalc && rows > 1) { + for(int j = 0; j < cols; j += cBlz) { + final int start = j; + final int end = Math.min(j + cBlz, cols); + tasks.add(pool.submit(() -> { + for(int i = start; i < end; i++) { + ifft_one_dim(re, im, re_inter, im_inter, i, i + rows * cols, rows, cols); + } + })); + } + + for(Future<?> f : tasks) + f.get(); + } + + tasks.clear(); + for(int i = 0; i < rows; i += rBlz) { + final int start = i; + final int end = Math.min(i + rBlz, rows); + tasks.add(pool.submit(() -> { + for(int j = start; j < end; j++) { + ifft_one_dim(re, im, re_inter, im_inter, j * cols, (j + 1) * cols, cols, 1); + } + })); + } + + for(Future<?> f : tasks) + f.get(); + } + catch(InterruptedException | ExecutionException e) { + throw new RuntimeException(e); + } + finally { + pool.shutdown(); + } + + } + + /** + * Function to perform one-dimensional FFT for two given double arrays. The first one represents the real values and + * the second one the imaginary values. Both arrays get updated and contain the result. + * + * @param re array representing the real values + * @param im array representing the imaginary values + * @param re_inter array for the real values of intermediate results + * @param im_inter array for the imaginary values of intermediate results + * @param start start index (incl.) + * @param stop stop index (excl.) + * @param num number of values used for fft + * @param minStep step size between values used for fft + */ + public static void fft_one_dim(double[] re, double[] im, double[] re_inter, double[] im_inter, int start, int stop, + int num, int minStep) { + + // start inclusive, stop exclusive + if(num == 1) + return; + + // even indices + for(int step = minStep * (num / 2), subNum = 2; subNum <= num; step /= 2, subNum *= 2) { + double angle = -2 * FastMath.PI / subNum; + fft_one_dim_iter_sub(re, im, re_inter, im_inter, start, stop, num, minStep, subNum, step, angle); + } + } + + /** + * Function to call one-dimensional FFT for all subArrays of two given double arrays. The first one represents the + * real values and the second one the imaginary values. Both arrays get updated and contain the result. + * + * @param re array representing the real values + * @param im array representing the imaginary values + * @param re_inter array for the real values of intermediate results + * @param im_inter array for the imaginary values of intermediate results + * @param start start index (incl.) + * @param stop stop index (excl.) + * @param num number of values used for fft + * @param minStep step size between values used for fft + * @param subNum number of elements in subarray + * @param step step size between elements in subarray + * @param angle angle + */ + private static void fft_one_dim_iter_sub(double[] re, double[] im, double[] re_inter, double[] im_inter, int start, + int stop, int num, int minStep, int subNum, int step, double angle) { + + // use ceil for the main (sub)array + for(int sub = 0; sub < FastMath.ceil(num / (2 * (double) subNum)); sub++) { + + int startSub = start + sub * minStep; + fft_one_dim_sub(re, im, re_inter, im_inter, start, stop, startSub, subNum, step, angle); + + // no odd values for main (sub)array + if(subNum == num) + return; + + startSub = start + sub * minStep + step / 2; + fft_one_dim_sub(re, im, re_inter, im_inter, start, stop, startSub, subNum, step, angle); + } + } + + /** + * Function to perform one-dimensional FFT for subArrays of two given double arrays. The first one represents the + * real values and the second one the imaginary values. Both arrays get updated and contain the result. + * + * @param re array representing the real values + * @param im array representing the imaginary values + * @param re_inter array for the real values of intermediate results + * @param im_inter array for the imaginary values of intermediate results + * @param start start index (incl.) + * @param stop stop index (excl.) + * @param startSub start index of subarray (incl.) + * @param subNum number of elements in subarray + * @param step step size between elements in subarray + * @param angle angle + */ + private static void fft_one_dim_sub(double[] re, double[] im, double[] re_inter, double[] im_inter, int start, + int stop, int startSub, int subNum, int step, double angle) { + + HashMap<Double, Double> sinCache = sinCacheRef.get(); + HashMap<Double, Double> cosCache = cosCacheRef.get(); + + if(sinCache == null) + sinCache = new HashMap<>(); + if(cosCache == null) + cosCache = new HashMap<>(); + + for(int j = startSub, cnt = 0; cnt < subNum / 2; j += 2 * step, cnt++) { + + double omega_pow_re = cos(cnt * angle, cosCache); + double omega_pow_im = sin(cnt * angle, sinCache); + + // calculate m using the result of odd index + double m_re = omega_pow_re * re[j + step] - omega_pow_im * im[j + step]; + double m_im = omega_pow_re * im[j + step] + omega_pow_im * re[j + step]; + + int index = startSub + cnt * step; + re_inter[index] = re[j] + m_re; + re_inter[index + (stop - start) / 2] = re[j] - m_re; + + im_inter[index] = im[j] + m_im; + im_inter[index + (stop - start) / 2] = im[j] - m_im; + } + + for(int j = startSub; j < startSub + (stop - start); j += step) { + re[j] = re_inter[j]; + im[j] = im_inter[j]; + re_inter[j] = 0; + im_inter[j] = 0; + } + } + + /** + * Function to perform one-dimensional IFFT for two given double arrays. The first one represents the real values and + * the second one the imaginary values. Both arrays get updated and contain the result. + * + * @param re array representing the real values + * @param im array representing the imaginary values + * @param re_inter array for the real values of intermediate results + * @param im_inter array for the imaginary values of intermediate results + * @param start start index (incl.) + * @param stop stop index (excl.) + * @param num number of values used for fft + * @param minStep step size between values used for fft + */ + private static void ifft_one_dim(double[] re, double[] im, double[] re_inter, double[] im_inter, int start, int stop, + int num, int minStep) { + + // conjugate input + for(int i = start; i < start + num * minStep; i += minStep) { + im[i] = -im[i]; + } + + // apply fft + fft_one_dim(re, im, re_inter, im_inter, start, stop, num, minStep); + + // conjugate and scale result + for(int i = start; i < start + num * minStep; i += minStep) { + re[i] = re[i] / num; + im[i] = -im[i] / num; + } + } + + /** + * Function for checking whether the given integer is a power of two. + * + * @param n integer to check + * @return true if n is a power of two, false otherwise + */ + public static boolean isPowerOfTwo(int n) { + return (n != 0) && ((n & (n - 1)) == 0); + } + + /** + * Function to perform FFT for a given matrix. The given matrix only represents real values. The output contains one + * matrix for the real and one for the imaginary values. + * + * @param re matrix object representing the real values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] fft(MatrixBlock re, int threads) { + return fft(re, + new MatrixBlock(re.getNumRows(), re.getNumColumns(), new double[re.getNumRows() * re.getNumColumns()]), + threads); + } + + /** + * Function to perform IFFT for a given matrix. The given matrix only represents real values. The output contains one + * matrix for the real and one for the imaginary values. + * + * @param re matrix object representing the real values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] ifft(MatrixBlock re, int threads) { + return ifft(re, + new MatrixBlock(re.getNumRows(), re.getNumColumns(), new double[re.getNumRows() * re.getNumColumns()]), + threads); + } + + /** + * Function to perform FFT for each row of a given matrix. The given matrix only represents real values. The output + * contains one matrix for the real and one for the imaginary values. + * + * @param re matrix object representing the real values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] fft_linearized(MatrixBlock re, int threads) { + return fft_linearized(re, + new MatrixBlock(re.getNumRows(), re.getNumColumns(), new double[re.getNumRows() * re.getNumColumns()]), + threads); + } + + /** + * Function to perform IFFT for each row of a given matrix. The given matrix only represents real values. The output + * contains one matrix for the real and one for the imaginary values. + * + * @param re matrix object representing the real values + * @param threads number of threads + * @return array of two matrix blocks + */ + public static MatrixBlock[] ifft_linearized(MatrixBlock re, int threads) { + return ifft_linearized(re, + new MatrixBlock(re.getNumRows(), re.getNumColumns(), new double[re.getNumRows() * re.getNumColumns()]), + threads); + } + + private static double sin(double angle, HashMap<Double, Double> cache) { + if(!cache.containsKey(angle)) { + cache.put(angle, FastMath.sin(angle)); + limitCache(cache); + } + return cache.get(angle); + } + + private static double cos(double angle, HashMap<Double, Double> cache) { + if(!cache.containsKey(angle)) { + cache.put(angle, FastMath.cos(angle)); + limitCache(cache); + } + return cache.get(angle); + } + + private static void limitCache(HashMap<Double, Double> cache) { + if(cache.size() > 1000) { + // remove oldest key + cache.remove(cache.keySet().iterator().next()); + } + } + +} diff --git a/src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java b/src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java new file mode 100644 index 0000000000..da844ad54a --- /dev/null +++ b/src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java @@ -0,0 +1,344 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you 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.apache.sysds.test.component.matrix; + +import org.apache.sysds.runtime.matrix.data.MatrixBlock; +import org.junit.Test; + +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft; +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft; +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft_linearized; +import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft_linearized; + +import static org.junit.Assert.assertArrayEquals; + +public class FourierTest { + + int threads = Runtime.getRuntime().availableProcessors(); + + @Test + public void test_fft_one_dim() { + + MatrixBlock re = new MatrixBlock(1, 4, new double[] {0, 18, -15, 3}); + MatrixBlock im = new MatrixBlock(1, 4, new double[4]); + + double[] expected_re = {6, 15, -36, 15}; + double[] expected_im = {0, -15, 0, 15}; + + fft(re, im, threads); + + assertArrayEquals(expected_re, re.getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, im.getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_fft_one_dim_2() { + + MatrixBlock re = new MatrixBlock(1, 8, new double[] {0, 18, -15, 3, 5, 10, 5, 9}); + MatrixBlock im = new MatrixBlock(1, 8, new double[8]); + + double[] expected_re = {35, 4.89949, 15, -14.89949, -45, -14.89949, 15, 4.89949}; + double[] expected_im = {0, 18.58579, -16, -21.41421, 0, 21.41421, 16, -18.58579}; + + fft(re, im, threads); + + assertArrayEquals(expected_re, re.getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, im.getDenseBlockValues(), 0.0001); + } + + @Test + public void test_fft_one_dim_matrixBlock() { + + MatrixBlock re = new MatrixBlock(1, 4, new double[] {0, 18, -15, 3}); + MatrixBlock im = new MatrixBlock(1, 4, new double[] {0, 0, 0, 0}); + + double[] expected_re = {6, 15, -36, 15}; + double[] expected_im = {0, -15, 0, 15}; + + MatrixBlock[] res = fft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_ifft_one_dim_matrixBlock_2() { + + double[] in_re = new double[] {1, -2, 3, -4}; + double[] in_im = new double[] {0, 0, 0, 0}; + + MatrixBlock re = new MatrixBlock(1, 4, in_re); + MatrixBlock im = new MatrixBlock(1, 4, in_im); + + MatrixBlock[] inter = fft(re, im, threads); + MatrixBlock[] res = ifft(inter[0], inter[1], threads); + + assertArrayEquals(in_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(in_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_fft_two_dim_matrixBlock() { + + MatrixBlock re = new MatrixBlock(2, 2, new double[] {0, 18, -15, 3}); + MatrixBlock im = new MatrixBlock(2, 2, new double[] {0, 0, 0, 0}); + + double[] expected_re = {6, -36, 30, 0}; + double[] expected_im = {0, 0, 0, 0}; + + MatrixBlock[] res = fft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_ifft_two_dim_matrixBlock() { + + MatrixBlock re = new MatrixBlock(2, 2, new double[] {6, -36, 30, 0}); + MatrixBlock im = new MatrixBlock(2, 2, new double[] {0, 0, 0, 0}); + + double[] expected_re = {0, 18, -15, 3}; + double[] expected_im = {0, 0, 0, 0}; + + MatrixBlock[] res = ifft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_fft_two_dim_matrixBlock_row_1() { + + MatrixBlock re = new MatrixBlock(1, 2, new double[] {0, 18}); + MatrixBlock im = new MatrixBlock(1, 2, new double[] {0, 0}); + + double[] expected_re = {18, -18}; + double[] expected_im = {0, 0}; + + MatrixBlock[] res = fft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_fft_two_dim_matrixBlock_row_2() { + + MatrixBlock re = new MatrixBlock(1, 2, new double[] {-15, 3}); + MatrixBlock im = new MatrixBlock(1, 2, new double[] {0, 0}); + + double[] expected_re = {-12, -18}; + double[] expected_im = {0, 0}; + + MatrixBlock[] res = fft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_ifft_with_complex_numpy_data() { + + // removed 0's at the end, not just real + MatrixBlock re = new MatrixBlock(1, 16, + new double[] {0.5398705320215192, 0.1793355360736929, 0.9065254044489506, 0.45004385530909075, + 0.11128090341119468, 0.11742862805303522, 0.7574827475752481, 0.2778193170985158, 0.9251562928110273, + 0.9414429667551927, 0.45131569795507087, 0.9522067687409731, 0.22491032260636257, 0.6579426733967295, + 0.7021558730366062, 0.7861117825617701,}); + + MatrixBlock im = new MatrixBlock(1, 16, new double[16]); + + double[] expected_re = {0.5613143313659362, -0.020146500453061336, 0.07086545895481336, -0.05003801442765281, + -0.06351635451036074, -0.03346768844048936, 0.07023899089706032, 0.007330763123826495, 0.016022890367311193, + 0.007330763123826495, 0.07023899089706032, -0.03346768844048936, -0.06351635451036074, -0.05003801442765281, + 0.07086545895481336, -0.020146500453061336}; + + double[] expected_im = {0.0, -0.07513090216687965, 0.023854392878864396, -0.018752997939582024, + -0.035626994964481226, -0.07808216015046443, 0.036579082654843526, -0.10605270957897009, 0.0, + 0.10605270957897009, -0.036579082654843526, 0.07808216015046443, 0.035626994964481226, 0.018752997939582024, + -0.023854392878864396, 0.07513090216687965}; + + MatrixBlock[] res = ifft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_ifft_2d_with_generated_data() { + + MatrixBlock re = new MatrixBlock(2, 2, + new double[] {0.6749989259154331, 0.6278845555308362, 0.995916990652601, 0.511472971081564}); + + MatrixBlock im = new MatrixBlock(2, 2, + new double[] {0.8330832079105173, 0.09857986129294982, 0.6808883894146879, 0.28353782431047303}); + + // adjusted the expected output + double[] expected_re = {0.70256836, 0.1328896, -0.05112662, -0.10933241}; + + double[] expected_im = {0.47402232, 0.28296348, -0.00819079, 0.0842882}; + + MatrixBlock[] res = ifft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_ifft_with_real_numpy_data() { + + // not complex + MatrixBlock re = new MatrixBlock(1, 16, + new double[] {0.17768499045697306, 0.3405035491673728, 0.9272906450602005, 0.28247504052271843, + 0.42775517613102865, 0.8338783039818357, 0.9813749624557385, 0.47224489612381737, 0.7936831995784907, + 0.5584182145651306, 0.5296113722056018, 0.6687593295928902, 0.9630598447622862, 0.7130539473424196, + 0.860081483892192, 0.8985058305053549}); + + MatrixBlock im = new MatrixBlock(1, 16, new double[16]); + + // adjusted the expected output + double[] expected_re = {0.6517738, -0.0263837, -0.03631354, -0.01644966, -0.05851095, -0.0849794, -0.01611732, + -0.02618679, 0.05579391, -0.02618679, -0.01611732, -0.0849794, -0.05851095, -0.01644966, -0.03631354, + -0.0263837}; + + double[] expected_im = {0, -0.04125649, -0.07121312, 0.02554502, 0.00774181, -0.08723921, -0.02314382, + -0.02021455, 0, 0.02021455, 0.02314382, 0.08723921, -0.00774181, -0.02554502, 0.07121312, 0.04125649 + + }; + + MatrixBlock[] res = ifft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_fft_two_dim_8_times_8() { + + MatrixBlock re = new MatrixBlock(8, 8, + new double[] {0.8435874964408077, 0.3565209485970835, 0.6221038572251737, 0.05712418097055716, + 0.9301368966310067, 0.7748052735242277, 0.21117129518682443, 0.08407931152930459, 0.5861235649815163, + 0.45860122035396356, 0.6647476180103304, 0.9167930424492593, 0.6310726270028377, 0.11110504251770592, + 0.32369996452324756, 0.5790548902504138, 0.5712310851880162, 0.5967356161025353, 0.6441861776319489, + 0.14402445187596158, 0.22642623625293545, 0.922443731897705, 0.9527667119829785, 0.2250880965427453, + 0.5755375055168817, 0.48898427237526954, 0.24518238824389693, 0.832292384016089, 0.23789083930394805, + 0.5558982102157535, 0.7220016080026206, 0.9666747522359772, 0.20509423975210916, 0.23170117015755587, + 0.7141206714718693, 0.2014158450611332, 0.6486924358372994, 0.9044990419216931, 0.19849364935627056, + 0.23340297110822106, 0.46854050631969246, 0.10134155509558795, 0.5563200388698989, 0.2669820016661475, + 0.8889445005077763, 0.4273462470993935, 0.8269490075576963, 0.044351336481537995, 0.3771564738915597, + 0.11333723996854606, 0.6913138435759023, 0.062431275099310124, 0.8003013976959878, 0.1276686539064662, + 0.975167392001707, 0.44595301043682656, 0.18401328301977316, 0.7158585484384759, 0.3240126702723025, + 0.740836665073052, 0.8890279623888511, 0.8841266040978419, 0.3058930798936259, 0.8987579873722049}); + + MatrixBlock im = new MatrixBlock(1, 16, + new double[] {0.8572457113722648, 0.668182795310341, 0.9739416721141464, 0.8189153345383146, + 0.6425950286263254, 0.3569634253534639, 0.19715070300424575, 0.8915344479242211, 0.39207930659031054, + 0.1625193685179268, 0.2523438052868171, 0.30940628850519547, 0.7461468672112159, 0.7123766750132684, + 0.5261041429273977, 0.867155304805022, 0.7207769261821749, 0.9139070611733158, 0.7638265842242135, + 0.3508092733308539, 0.6075639148195967, 0.9615531048215422, 0.719499617407839, 0.9616615941848492, + 0.2667126256574347, 0.8215093145949468, 0.4240476512138287, 0.5015798652459079, 0.19784651066995873, + 0.42315603332105356, 0.5575575283922164, 0.9051304828282485, 0.30117855478511435, 0.14219967492505514, + 0.32675429179906557, 0.04889894374947912, 0.8338579676700041, 0.370201089804747, 0.06025987717830994, + 0.9407970353033787, 0.9871788482561391, 0.75984297199074, 0.414969247979073, 0.2453785474698862, + 0.06295683447294731, 0.40141192931768566, 0.19520663793867488, 0.3179027928938928, 0.591138083168947, + 0.5318366162549014, 0.04865894304644136, 0.5339043989658795, 0.09892519435896363, 0.31616794516128466, + 0.06702286400447643, 0.8466767273121609, 0.8134875055724791, 0.6232554321597641, 0.21208039111457444, + 0.25629831822305926, 0.7373140896724466, 0.020486629088602437, 0.8666668269441752, 0.20094387974200512}); + + double[] expected_re = {32.51214260297584, -3.4732779237490314, -0.7257760912890102, -1.9627494786611792, + 3.571671446098747, 1.0451692206901078, 0.8970702451384204, -1.3739767803210428, 4.892442103095981, + -1.1656855109832338, -0.5854908742291178, -1.3497699084098418, -1.377003693155216, -2.1698030461214923, + 0.8172129683973663, -0.9259076518379679, -1.1343756245445045, -1.8734967800709579, 1.7367517585478862, + 0.07349671655414491, -1.5933768052439223, 2.7965196291943983, 4.292588673604611, -1.1032899622026413, + -2.4643093702874834, 2.109128987930992, 3.2834030498896456, 0.21371926254596152, -0.3107488550365316, + 0.7293395030253796, -2.542403789759091, -1.8570654162590052, -2.325781245331303, 0.7963395911053484, + -2.351990667205867, -2.4241188304485735, 4.689766636746301, 3.4748121457306116, 0.5539071663846459, + -0.950099313504134, 2.122310975524349, 4.527637759721644, -2.125596093625001, 1.7676565539001592, + 5.748332643019926, 0.860140632830907, 2.9735142186218484, -1.7198774815194848, -0.18418859401548549, + 1.1909629561342188, 0.21710627714418418, -1.5537184277268996, -0.5486540814747869, 0.14807346060743987, + 2.4333154010438087, -3.2930077637380393, -2.3820067665775113, 2.5463581304688057, 2.5927580716559615, + 1.8921802492721915, 0.4957713559465988, -0.4983536537999108, 3.5808175362367676, 0.7530823235547575}; + + double[] expected_im = {32.64565805549281, 5.177639365468945, 1.1792020104097647, -0.4850627423320939, + -1.719468548169175, -3.064146170894837, 3.3226243586118906, 2.3819341640916107, 1.5824429804361388, + -2.192940882164737, 0.5774407122593543, 0.16873948200103983, 4.297014293326352, -3.1712082122026883, + 0.9741291131305898, -2.4929883795121235, 1.111763301820595, -1.4012254390671657, -0.33687898317382636, + 2.324190267133635, -2.8862969254091397, -4.7558982401265135, 1.8244587481290004, 0.5310550630270396, + 2.655726742689745, 2.510014260306531, 0.25589537824783704, 1.8720307201415736, -2.6458046644482884, + 2.1732611302115585, -2.5162250969793227, -0.9103444457919911, 2.2835527482590248, 0.5187392677625127, + -3.335253420903965, 1.4668560670097441, -1.9681585205341436, -2.81914578771063, 4.818094364700921, + -0.877636803126361, 1.803174743823159, 3.1346192487664277, -3.564058675191744, -0.3391381913837902, + -1.2897867384105863, 2.315065426377637, 1.5764817121472765, 2.412894091248795, -2.3182678917385218, + -3.057303547366563, 0.033996764974414395, -0.5825423640666696, 6.395088232674363, -0.5553659624089358, + 1.1079219041153268, 0.1094531803830765, 3.488182265163636, -1.5698242466218544, -0.1013387045518459, + 0.9269290699615746, -0.699890233104248, 3.617209720991753, -0.5565163478425035, 3.502962737763559}; + + MatrixBlock[] res = fft(re, im, threads); + + assertArrayEquals(expected_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, res[1].getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_fft_linearized() { + + MatrixBlock re = new MatrixBlock(2, 4, new double[] {0, 18, -15, 3, 0, 18, -15, 3}); + MatrixBlock im = new MatrixBlock(1, 4, new double[8]); + + double[] expected_re = {6, 15, -36, 15, 6, 15, -36, 15}; + double[] expected_im = {0, -15, 0, 15, 0, -15, 0, 15}; + + fft_linearized(re, im, threads); + + assertArrayEquals(expected_re, re.getDenseBlockValues(), 0.0001); + assertArrayEquals(expected_im, im.getDenseBlockValues(), 0.0001); + + } + + @Test + public void test_ifft_linearized() { + + double[] in_re = new double[] {1, -2, 3, -4, 1, -2, 3, -4}; + double[] in_im = new double[] {0, 0, 0, 0, 0, 0, 0, 0}; + + MatrixBlock re = new MatrixBlock(2, 4, in_re); + MatrixBlock im = new MatrixBlock(2, 4, in_im); + + MatrixBlock[] inter = fft_linearized(re, im, threads); + MatrixBlock[] res = ifft_linearized(inter[0], inter[1], threads); + + assertArrayEquals(in_re, res[0].getDenseBlockValues(), 0.0001); + assertArrayEquals(in_im, res[1].getDenseBlockValues(), 0.0001); + + } + +}
