Baunsgaard commented on code in PR #1983:
URL: https://github.com/apache/systemds/pull/1983#discussion_r1452814727
##########
src/main/java/org/apache/sysds/parser/BuiltinFunctionExpression.java:
##########
@@ -355,6 +355,62 @@ public void validateExpression(MultiAssignmentStatement
stmt, HashMap<String, Da
break;
}
+ case FFT:
+ {
+ checkNumParameters(1);
+ checkMatrixParam(getFirstExpr());
+
+ // setup output properties
+ DataIdentifier fftOut1 = (DataIdentifier)
getOutputs()[0];
+ DataIdentifier fftOut2 = (DataIdentifier)
getOutputs()[1];
+
+ // if (getFirstExpr().getOutput().getDim2() != 1 ||
getFirstExpr().getOutput().getDim2() != 2) {
+ // raiseValidateError("Eigen Decomposition can
only be done on a square matrix. Input matrix is rectangular (rows=" +
getFirstExpr().getOutput().getDim1() + ", cols="+
getFirstExpr().getOutput().getDim2() +")", conditional);
+ // }
Review Comment:
initially change this to a TODO: add verification check.
##########
pom.xml:
##########
@@ -875,6 +875,12 @@
</profiles>
<dependencies>
+ <dependency>
+ <groupId>org.json</groupId>
+ <artifactId>json</artifactId>
+ <version>20210307</version>
+ </dependency>
Review Comment:
use Jackson instead if you need to parse JSON, we use the Jackson dependency
already.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
Review Comment:
license missing
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java:
##########
@@ -252,6 +258,25 @@ private static EigenDecomposition
computeEigenRegularized(MatrixBlock in) {
DataConverter.convertToArray2DRowRealMatrix(in2));
}
+ private static MatrixBlock[] computeFFT(MatrixBlock in) {
+ if( in == null || in.isEmptyBlock(false) )
+ throw new DMLRuntimeException("Invalid empty block");
+
+ //run fft
Review Comment:
add a sparseToDense call.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
Review Comment:
this also allocate.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+package org.apache.sysds.test.component.matrix;
Review Comment:
license
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
Review Comment:
java doc require you to define `param` and `returns`, see other java doc
examples.
I suggest adding it for all public methods.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTestData.py:
##########
@@ -0,0 +1,165 @@
+import numpy as np
Review Comment:
add license, if we want to keep this file.
Unfortunately it is complicated for us to run this script before we run the
Fourier test, and therefore
the tests would not parse in the cloud because of missing input.
We need to decided what to do.
1. include the generated files in git. (i can live with this if they are
fairly small, but i would prefer not to)
2. Move the tests to execute fully in our PythonAPI, this would be cool, but
some effort to write new test cases there, and add a shallow wrapper fft
function (3-8 lines of code). If you want to do this i can help out quickly in
a call to fix it.
3. Remove it and use smaller controlled test cases in the java test that
have known outputs that you hard code (we need a few more of these hard coded
cases anyway in my opinion )
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
Review Comment:
allocate
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTestWithFiles.java:
##########
@@ -0,0 +1,377 @@
+package org.apache.sysds.test.component.matrix;
+
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.junit.Test;
+
+import java.io.BufferedReader;
+import java.io.FileReader;
+import java.io.IOException;
+
+import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.*;
+import static org.junit.Assert.assertTrue;
+
+public class FourierTestWithFiles {
+ int progressInterval = 5000;
+
+ // prior to executing the following tests it is necessary to run the Numpy
Script in FourierTestData.py
+ // and add the generated files to the root of the project.
+ @Test
+ public void testFftWithNumpyData() throws IOException {
+ String filename = "fft_data.csv"; // Path to your CSV file
+ BufferedReader reader = new BufferedReader(new FileReader(filename));
+ String line;
+ int lineNumber = 0;
+ long totalTime = 0; // Total time for all FFT computations
+ int numCalculations = 0; // Number of FFT computations
+
+ while ((line = reader.readLine()) != null) {
+ lineNumber++;
+
+ String[] values = line.split(",");
+ int n = values.length / 3;
+ double[][][] input = new double[2][1][n];
+ double[][] expected = new double[2][n]; // First row for real,
second row for imaginary parts
+
+ for (int i = 0; i < n; i++) {
+ input[0][0][i] = Double.parseDouble(values[i]);
+ expected[0][i] = Double.parseDouble(values[n + i]); // Real
part
+ expected[1][i] = Double.parseDouble(values[n * 2 + i]); //
Imaginary part
+ }
+
+ long startTime = System.nanoTime();
+ MatrixBlock[] actualBlocks = fft(input);
+ long endTime = System.nanoTime();
+
+ if(lineNumber > 1000){
+ totalTime += (endTime - startTime);
+ numCalculations++;
+
+ if (numCalculations % progressInterval == 0) {
+ double averageTime = (totalTime / 1e6) / numCalculations;
// Average time in milliseconds
+ System.out.println("fft(double[][][] in): Average
execution time after " + numCalculations + " calculations: " +
String.format("%.8f", averageTime/1000) + " s");
+ }
+ }
+
+ // Validate the FFT results
+ validateFftResults(expected, actualBlocks, lineNumber);
+ }
+
+ reader.close();
+
+ }
+
+ private void validateFftResults(double[][] expected, MatrixBlock[]
actualBlocks, int lineNumber) {
+ int length = expected[0].length;
+ for (int i = 0; i < length; i++) {
+ double realActual = actualBlocks[0].getValueDenseUnsafe(0, i);
+ double imagActual = actualBlocks[1].getValueDenseUnsafe(0, i);
+ assertEquals("Mismatch in real part at index " + i + " in line " +
lineNumber, expected[0][i], realActual, 1e-9);
+ assertEquals("Mismatch in imaginary part at index " + i + " in
line " + lineNumber, expected[1][i], imagActual, 1e-9);
+ }
+ if(lineNumber % progressInterval == 0){
+ System.out.println("fft(double[][][] in): Finished processing line
" + lineNumber);
+ }
+
+ }
+
+ @Test
+ public void testFftExecutionTime() throws IOException {
+ String filename = "fft_data.csv"; // Path to your CSV file
+ BufferedReader reader = new BufferedReader(new FileReader(filename));
+ String line;
+ int lineNumber = 0;
+ long totalTime = 0; // Total time for all FFT computations
+ int numCalculations = 0; // Number of FFT computations
+
+ while ((line = reader.readLine()) != null) {
+ lineNumber++;
+ String[] values = line.split(",");
+ int n = values.length / 3;
+ double[][][] input = new double[2][1][n];
+
+ for (int i = 0; i < n; i++) {
+ input[0][0][i] = Double.parseDouble(values[i]); // Real part
+ input[1][0][i] = Double.parseDouble(values[n + i]); //
Imaginary part
+ }
+
+ long startTime = System.nanoTime();
+ fft(input, false);
+ long endTime = System.nanoTime();
+ if(lineNumber > 1000){
+ totalTime += (endTime - startTime);
+ numCalculations++;
+
+ if (numCalculations % progressInterval == 0) {
+ double averageTime = (totalTime / 1e6) / numCalculations;
// Average time in milliseconds
+ System.out.println("fft(double[][][] in, boolean calcInv)
Average execution time after " + numCalculations + " calculations: " +
String.format("%.8f", averageTime/1000) + " s");
+ }
+ }
+ }
+
+ reader.close();
+ }
+
+ @Test
+ public void testFftExecutionTimeOfOneDimFFT() throws IOException {
+ String filename = "fft_data.csv"; // Path to your CSV file
+ BufferedReader reader = new BufferedReader(new FileReader(filename));
+ String line;
+ int lineNumber = 0;
+ long totalTime = 0; // Total time for all FFT computations
+ int numCalculations = 0; // Number of FFT computations
+
+ while ((line = reader.readLine()) != null) {
+ lineNumber++;
+ String[] values = line.split(",");
+ int n = values.length / 2;
+ double[][] input = new double[2][n]; // First row for real, second
row for imaginary parts
+
+ for (int i = 0; i < n; i++) {
+ input[0][i] = Double.parseDouble(values[i]); // Real part
+ input[1][i] = Double.parseDouble(values[n + i]); // Imaginary
part
+ }
+
+ long startTime = System.nanoTime();
+ fft_one_dim(input);
+ long endTime = System.nanoTime();
+ if(lineNumber > 1000){
+ totalTime += (endTime - startTime);
+ numCalculations++;
+
+ if (numCalculations % progressInterval == 0) {
+ double averageTime = (totalTime / 1e6) / numCalculations;
// Average time in milliseconds
+ System.out.println("fft_one_dim: Average execution time
after " + numCalculations + " calculations: " + String.format("%.8f",
averageTime/1000) + " s ");
+ }
+ }
+ }
+
+ reader.close();
+ }
+
+
+ // prior to executing this test it is necessary to run the Numpy Script in
FourierTestData.py and add the generated file to the root of the project.
+ @Test
+ public void testIfftWithRealNumpyData() throws IOException {
+ String filename = "ifft_data.csv"; // Path to your CSV file
+ BufferedReader reader = new BufferedReader(new FileReader(filename));
+ String line;
+ int lineNumber = 0;
+
+ while ((line = reader.readLine()) != null) {
+ lineNumber++;
+ String[] values = line.split(",");
+ int n = values.length / 3;
+ double[][][] input = new double[2][1][n];
+ double[][] expected = new double[2][n]; // First row for real,
second row for imaginary parts
+
+ for (int i = 0; i < n; i++) {
+ input[0][0][i] = Double.parseDouble(values[i]); // Real part
of input
+ // Imaginary part of input is assumed to be 0
+ expected[0][i] = Double.parseDouble(values[n + i]); // Real
part of expected output
+ expected[1][i] = Double.parseDouble(values[n * 2 + i]); //
Imaginary part of expected output
+ }
+
+ double[][][] actualResult = fft(input, true); // Perform IFFT
+
+ // Validate the IFFT results
+ validateFftResults(expected, actualResult, lineNumber);
+ }
+
+ reader.close();
+ }
+
+ private void validateFftResults(double[][] expected, double[][][]
actualResult, int lineNumber) {
+ int length = expected[0].length;
+ for (int i = 0; i < length; i++) {
+ double realActual = actualResult[0][0][i];
+ double imagActual = actualResult[1][0][i];
+ assertEquals("Mismatch in real part at index " + i + " in line " +
lineNumber, expected[0][i], realActual, 1e-9);
+ assertEquals("Mismatch in imaginary part at index " + i + " in
line " + lineNumber, expected[1][i], imagActual, 1e-9);
+ }
+ if(lineNumber % progressInterval == 0){
+ System.out.println("ifft(real input): Finished processing line " +
lineNumber);
+ }
+
+ }
+
+
+
+ @Test
+public void testIfftWithComplexNumpyData() throws IOException {
Review Comment:
indentation
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
Review Comment:
if possible can we use the already allocated res matrix, and index into
cells not calculated yet?
##########
src/main/java/org/apache/sysds/runtime/instructions/cp/MultiReturnBuiltinCPInstruction.java:
##########
@@ -90,6 +90,23 @@ else if ( opcode.equalsIgnoreCase("eigen") ) {
return new MultiReturnBuiltinCPInstruction(null, in1,
outputs, opcode, str);
}
+ else if ( 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 ( 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);
Review Comment:
these two can probably share code.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
Review Comment:
This is a good change but we need to change a few things.
The convertToArray together with convertToMatrixBlocks is dominating the
execution time,
since the implementation now allocate new double[][] arrays in the size of
input and output, and then convert the output once again. In total the
algorithm allocate ~6x the original input.
Instead of please take the input MatrixBlock and call .getDenseBlockValues()
This return an already allocated linearized double[] array of the values
contained.
The algorithm should only allocate two output real double[] and imaginary
double[] array making the total allocation 2x.
do not use double[][] nor double[][][] in location of the algorithm
implementation.
The change should not take that long, since most of the loops then become
some modification of calculating indexes in the matrices while the FFT
computation is the same.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java:
##########
@@ -71,7 +73,7 @@ public static boolean isSupportedUnaryOperation( String
opcode ) {
}
public static boolean isSupportedMultiReturnOperation( String opcode ) {
- return ( opcode.equals("qr") || opcode.equals("lu") ||
opcode.equals("eigen") || opcode.equals("svd") );
+ return ( opcode.equals("qr") || opcode.equals("lu") ||
opcode.equals("eigen") || opcode.equals("fft") || opcode.equals("ifft") ||
opcode.equals("svd") );
Review Comment:
could you change this entire method to a switch statement thanks.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java:
##########
@@ -49,6 +49,8 @@
import org.apache.sysds.runtime.matrix.operators.BinaryOperator;
import org.apache.sysds.runtime.util.DataConverter;
+import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.*;
Review Comment:
we do not allow wildcard imports. you probably have to configure your IDE
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java:
##########
@@ -252,6 +258,25 @@ private static EigenDecomposition
computeEigenRegularized(MatrixBlock in) {
DataConverter.convertToArray2DRowRealMatrix(in2));
}
+ private static MatrixBlock[] computeFFT(MatrixBlock in) {
+ if( in == null || in.isEmptyBlock(false) )
+ throw new DMLRuntimeException("Invalid empty block");
+
+ //run fft
+ return fft(in);
+ }
+
+ private static MatrixBlock[] computeIFFT(MatrixBlock in) {
+ if( in == null || in.isEmptyBlock(false))
+ throw new DMLRuntimeException("Invalid empty block");
+ int rows = in.getNumRows();
+ int cols = in.getNumColumns();
+
+ MatrixBlock inIm = new MatrixBlock(rows, cols, new
double[cols*rows]);
Review Comment:
is it possible not to allocate this, maybe we need two different kernels?
the allocation of this I part is very expensive.
also you probably want to modify the IFFT to optionally allow two inputs
from DML script level.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
+ double[] flattened_im = new double[rows*cols];
+ if(in.length > 1){
+ flattened_im =
Arrays.stream(in[1]).flatMapToDouble(Arrays::stream).toArray();
+ }
+
+ MatrixBlock re = new MatrixBlock(rows, cols, flattened_re);
+ MatrixBlock im = new MatrixBlock(rows, cols, flattened_im);
+
+ return new MatrixBlock[]{re, im};
+ }
+
+ private static MatrixBlock getZeroMatrixBlock(int rows, int cols){
Review Comment:
remove this method.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
+ double[] flattened_im = new double[rows*cols];
+ if(in.length > 1){
+ flattened_im =
Arrays.stream(in[1]).flatMapToDouble(Arrays::stream).toArray();
+ }
+
+ MatrixBlock re = new MatrixBlock(rows, cols, flattened_re);
+ MatrixBlock im = new MatrixBlock(rows, cols, flattened_im);
+
+ return new MatrixBlock[]{re, im};
+ }
+
+ private static MatrixBlock getZeroMatrixBlock(int rows, int cols){
+
+ return new MatrixBlock(rows, cols, new double[cols*rows]);
+
+ }
+
+ private static double[][] convertToArray(MatrixBlock in){
+
+ int rows = in.getNumRows();
+ int cols = in.getNumColumns();
+
+ double[][] out = new double[rows][cols];
+ for(int i = 0; i < rows; i++){
+ out[i] = Arrays.copyOfRange(in.getDenseBlockValues(), i * cols,
(i+1) * cols);
+ }
+
+ return out;
+ }
+ private static double[][][] convertToArray(MatrixBlock[] in){
Review Comment:
avoid, and remove
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
+ double[] flattened_im = new double[rows*cols];
+ if(in.length > 1){
+ flattened_im =
Arrays.stream(in[1]).flatMapToDouble(Arrays::stream).toArray();
+ }
+
+ MatrixBlock re = new MatrixBlock(rows, cols, flattened_re);
+ MatrixBlock im = new MatrixBlock(rows, cols, flattened_im);
+
+ return new MatrixBlock[]{re, im};
+ }
+
+ private static MatrixBlock getZeroMatrixBlock(int rows, int cols){
+
+ return new MatrixBlock(rows, cols, new double[cols*rows]);
+
+ }
+
+ private static double[][] convertToArray(MatrixBlock in){
Review Comment:
remove, and avoid this
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
+ double[] flattened_im = new double[rows*cols];
+ if(in.length > 1){
+ flattened_im =
Arrays.stream(in[1]).flatMapToDouble(Arrays::stream).toArray();
+ }
+
+ MatrixBlock re = new MatrixBlock(rows, cols, flattened_re);
+ MatrixBlock im = new MatrixBlock(rows, cols, flattened_im);
+
+ return new MatrixBlock[]{re, im};
+ }
+
+ private static MatrixBlock getZeroMatrixBlock(int rows, int cols){
+
+ return new MatrixBlock(rows, cols, new double[cols*rows]);
+
+ }
+
+ private static double[][] convertToArray(MatrixBlock in){
+
+ int rows = in.getNumRows();
+ int cols = in.getNumColumns();
+
+ double[][] out = new double[rows][cols];
+ for(int i = 0; i < rows; i++){
+ out[i] = Arrays.copyOfRange(in.getDenseBlockValues(), i * cols,
(i+1) * cols);
+ }
+
+ return out;
+ }
+ private static double[][][] convertToArray(MatrixBlock[] in){
+
+ int rows = in[0].getNumRows();
+ int cols = in[0].getNumColumns();
+
+ double[][][] out = new double[2][rows][cols];
+ for(int k = 0; k < 2; k++){
+ for(int i = 0; i < rows; i++){
+ out[k][i] = Arrays.copyOfRange(in[k].getDenseBlockValues(), i
* cols, (i+1) * cols);
+ }
+ }
+
+ return out;
+ }
+
+ public static double[][] get_complex_row(double[][][] in, int i){
+
+ int cols = in[0][0].length;
+
+ double[][] row = new double[2][cols];
+ // get row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ row[k][j] = in[k][i][j];
+ }
+ }
+ return row;
+ }
+ public static double[][] get_complex_col(double[][][] in, int j){
+
+ int rows = in[0].length;
+
+ double[][] col = new double[2][rows];
+ // get row
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ col[k][i] = in[k][i][j];
+ }
+ }
+ return col;
+ }
+
+ private static boolean isPowerOfTwo(int n){
+ return ((n != 0) && ((n & (n - 1)) == 0)) || n == 1;
+ }
+
+ public static MatrixBlock[] fft(double[] in){
+ double[][][] arr = new double[2][1][in.length];
+ arr[0][0] = in;
+ return fft(convertToMatrixBlocks(arr));
+ }
+
+ public static MatrixBlock[] fft(double[][][] in){
+ return fft(convertToMatrixBlocks(in));
+ }
+
+ public static MatrixBlock[] fft(MatrixBlock[] in){
+ return fft(in[0], in[1]);
+ }
+
+ public static MatrixBlock[] fft(MatrixBlock re){
+ return fft(re, getZeroMatrixBlock(re.getNumRows(),
re.getNumColumns()));
Review Comment:
a zero matrix block is the default, and we do not have to allocate the
double array underneath, because we can simply use an sparse / empty matrix.
(The only way you would know this is by knowing the system)
simply write `new MatrixBlock(rows,cols, true);`
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
Review Comment:
another allocation
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
+ double[] flattened_im = new double[rows*cols];
+ if(in.length > 1){
+ flattened_im =
Arrays.stream(in[1]).flatMapToDouble(Arrays::stream).toArray();
+ }
+
+ MatrixBlock re = new MatrixBlock(rows, cols, flattened_re);
+ MatrixBlock im = new MatrixBlock(rows, cols, flattened_im);
+
+ return new MatrixBlock[]{re, im};
+ }
+
+ private static MatrixBlock getZeroMatrixBlock(int rows, int cols){
+
+ return new MatrixBlock(rows, cols, new double[cols*rows]);
+
+ }
+
+ private static double[][] convertToArray(MatrixBlock in){
+
+ int rows = in.getNumRows();
+ int cols = in.getNumColumns();
+
+ double[][] out = new double[rows][cols];
+ for(int i = 0; i < rows; i++){
+ out[i] = Arrays.copyOfRange(in.getDenseBlockValues(), i * cols,
(i+1) * cols);
+ }
+
+ return out;
+ }
+ private static double[][][] convertToArray(MatrixBlock[] in){
+
+ int rows = in[0].getNumRows();
+ int cols = in[0].getNumColumns();
+
+ double[][][] out = new double[2][rows][cols];
+ for(int k = 0; k < 2; k++){
+ for(int i = 0; i < rows; i++){
+ out[k][i] = Arrays.copyOfRange(in[k].getDenseBlockValues(), i
* cols, (i+1) * cols);
+ }
+ }
+
+ return out;
+ }
+
+ public static double[][] get_complex_row(double[][][] in, int i){
Review Comment:
we should avoid this as well, since it simply index into the matrix, and
extract the cells.
This only makes the implementation slow, and trade a bit of coding.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+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.*;
+import static org.junit.Assert.assertArrayEquals;
+
+public class FourierTest {
+
+ @Test
Review Comment:
change indentation to tabs
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
+ double[] flattened_im = new double[rows*cols];
+ if(in.length > 1){
+ flattened_im =
Arrays.stream(in[1]).flatMapToDouble(Arrays::stream).toArray();
+ }
+
+ MatrixBlock re = new MatrixBlock(rows, cols, flattened_re);
+ MatrixBlock im = new MatrixBlock(rows, cols, flattened_im);
+
+ return new MatrixBlock[]{re, im};
+ }
+
+ private static MatrixBlock getZeroMatrixBlock(int rows, int cols){
+
+ return new MatrixBlock(rows, cols, new double[cols*rows]);
+
+ }
+
+ private static double[][] convertToArray(MatrixBlock in){
+
+ int rows = in.getNumRows();
+ int cols = in.getNumColumns();
+
+ double[][] out = new double[rows][cols];
+ for(int i = 0; i < rows; i++){
+ out[i] = Arrays.copyOfRange(in.getDenseBlockValues(), i * cols,
(i+1) * cols);
+ }
+
+ return out;
+ }
+ private static double[][][] convertToArray(MatrixBlock[] in){
+
+ int rows = in[0].getNumRows();
+ int cols = in[0].getNumColumns();
+
+ double[][][] out = new double[2][rows][cols];
+ for(int k = 0; k < 2; k++){
+ for(int i = 0; i < rows; i++){
+ out[k][i] = Arrays.copyOfRange(in[k].getDenseBlockValues(), i
* cols, (i+1) * cols);
+ }
+ }
+
+ return out;
+ }
+
+ public static double[][] get_complex_row(double[][][] in, int i){
+
+ int cols = in[0][0].length;
+
+ double[][] row = new double[2][cols];
+ // get row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ row[k][j] = in[k][i][j];
+ }
+ }
+ return row;
+ }
+ public static double[][] get_complex_col(double[][][] in, int j){
Review Comment:
index into the input, and remove as well.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
+ * Function to perform Fast Fourier Transformation
+ */
+
+ public static MatrixBlock[] fft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, false);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static MatrixBlock[] ifft(MatrixBlock re, MatrixBlock im){
+
+ int rows = re.getNumRows();
+ int cols = re.getNumColumns();
+
+ double[][][] in = new double[2][rows][cols];
+ in[0] = convertToArray(re);
+ in[1] = convertToArray(im);
+
+ double[][][] res = fft(in, true);
+
+ return convertToMatrixBlocks(res);
+ }
+
+ public static double[][][] fft(double[][][] in, boolean calcInv){
+
+ int rows = in[0].length;
+ int cols = in[0][0].length;
+
+ double[][][] res = new double[2][rows][cols];
+
+ for(int i = 0; i < rows; i++){
+ // use fft or ifft on each row
+ double[][] res_row = calcInv? ifft_one_dim(get_complex_row(in, i))
: fft_one_dim(get_complex_row(in, i));
+
+ // set res row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_row[k][j];
+ }
+ }
+ }
+
+ if(rows == 1) return res;
+
+ for(int j = 0; j < cols; j++){
+ // use fft on each col
+ double[][] res_col = calcInv? ifft_one_dim(get_complex_col(res,
j)) : fft_one_dim(get_complex_col(res, j));
+
+ // set res col
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ res[k][i][j] = res_col[k][i];
+ }
+ }
+ }
+
+ return res;
+ }
+
+ public static double[][] fft_one_dim(double[][] in){
+ // 1st row real part, 2nd row imaginary part
+ if(in == null || in.length != 2 || in[0].length != in[1].length) throw
new RuntimeException("in false dimensions");
+
+ int cols = in[0].length;
+ if(cols == 1) return in;
+
+ double angle = -2*FastMath.PI/cols;
+
+ // split values depending on index
+ double[][] even = new double[2][cols/2];
+ double[][] odd = new double[2][cols/2];
+
+ for(int i = 0; i < 2; i++){
+ for (int j = 0; j < cols/2; j++){
+ even[i][j] = in[i][j*2];
+ odd[i][j] = in[i][j*2+1];
+ }
+ }
+ double[][] res_even = fft_one_dim(even);
+ double[][] res_odd = fft_one_dim(odd);
+
+ double[][] res = new double[2][cols];
+
+ for(int j=0; j < cols/2; j++){
+ double[] omega_pow = new double[]{FastMath.cos(j*angle),
FastMath.sin(j*angle)};
+
+ // m = omega * res_odd[j]
+ double[] m = new double[]{
+ omega_pow[0] * res_odd[0][j] - omega_pow[1] *
res_odd[1][j],
+ omega_pow[0] * res_odd[1][j] + omega_pow[1] *
res_odd[0][j]};
+
+ // res[j] = res_even + m;
+ // res[j+cols/2] = res_even - m;
+ for(int i = 0; i < 2; i++){
+ res[i][j] = res_even[i][j] + m[i];
+ res[i][j+cols/2] = res_even[i][j] - m[i];
+ }
+ }
+
+ return res;
+
+ }
+
+ public static double[][] ifft_one_dim(double[][] in) {
+
+ // cols[0] is real part, cols[1] is imaginary part
+ int cols = in[0].length;
+
+ // conjugate input
+ in[1] = Arrays.stream(in[1]).map(i -> -i).toArray();
+
+ // apply fft
+ double[][] res = fft_one_dim(in);
+
+ // conjugate and scale result
+ res[0] = Arrays.stream(res[0]).map(i -> i/cols).toArray();
+ res[1] = Arrays.stream(res[1]).map(i -> -i/cols).toArray();
+
+ return res;
+ }
+
+ private static MatrixBlock[] convertToMatrixBlocks(double[][][] in){
+
+ int cols = in[0][0].length;
+ int rows = in[0].length;
+
+ double[] flattened_re =
Arrays.stream(in[0]).flatMapToDouble(Arrays::stream).toArray();
+ double[] flattened_im = new double[rows*cols];
+ if(in.length > 1){
+ flattened_im =
Arrays.stream(in[1]).flatMapToDouble(Arrays::stream).toArray();
+ }
+
+ MatrixBlock re = new MatrixBlock(rows, cols, flattened_re);
+ MatrixBlock im = new MatrixBlock(rows, cols, flattened_im);
+
+ return new MatrixBlock[]{re, im};
+ }
+
+ private static MatrixBlock getZeroMatrixBlock(int rows, int cols){
+
+ return new MatrixBlock(rows, cols, new double[cols*rows]);
+
+ }
+
+ private static double[][] convertToArray(MatrixBlock in){
+
+ int rows = in.getNumRows();
+ int cols = in.getNumColumns();
+
+ double[][] out = new double[rows][cols];
+ for(int i = 0; i < rows; i++){
+ out[i] = Arrays.copyOfRange(in.getDenseBlockValues(), i * cols,
(i+1) * cols);
+ }
+
+ return out;
+ }
+ private static double[][][] convertToArray(MatrixBlock[] in){
+
+ int rows = in[0].getNumRows();
+ int cols = in[0].getNumColumns();
+
+ double[][][] out = new double[2][rows][cols];
+ for(int k = 0; k < 2; k++){
+ for(int i = 0; i < rows; i++){
+ out[k][i] = Arrays.copyOfRange(in[k].getDenseBlockValues(), i
* cols, (i+1) * cols);
+ }
+ }
+
+ return out;
+ }
+
+ public static double[][] get_complex_row(double[][][] in, int i){
+
+ int cols = in[0][0].length;
+
+ double[][] row = new double[2][cols];
+ // get row
+ for (int j = 0; j < cols; j++){
+ for( int k = 0; k < 2; k++){
+ row[k][j] = in[k][i][j];
+ }
+ }
+ return row;
+ }
+ public static double[][] get_complex_col(double[][][] in, int j){
+
+ int rows = in[0].length;
+
+ double[][] col = new double[2][rows];
+ // get row
+ for (int i = 0; i < rows; i++){
+ for( int k = 0; k < 2; k++){
+ col[k][i] = in[k][i][j];
+ }
+ }
+ return col;
+ }
+
+ private static boolean isPowerOfTwo(int n){
+ return ((n != 0) && ((n & (n - 1)) == 0)) || n == 1;
Review Comment:
`(n != 0) && ((n & (n - 1)) == 0)`
is sufficient.
also consider if you ever would input 0.
imagine putting in n = 1, then the second part is 1 & 0 == 0, which is true.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+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.*;
Review Comment:
no wildcard imports
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+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.*;
+import static org.junit.Assert.assertArrayEquals;
+
+public class FourierTest {
+
+ @Test
+ public void simple_test_one_dim() {
+ // 1st row real part, 2nd row imaginary part
+ double[][] in = {{0, 18, -15, 3},{0, 0, 0, 0}};
Review Comment:
we need to change the test inputs to MatrixBlock, or linearized double[]
arrays.
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibMatrixFourier.java:
##########
@@ -0,0 +1,236 @@
+package org.apache.sysds.runtime.matrix.data;
+
+import org.apache.commons.math3.util.FastMath;
+import java.util.Arrays;
+
+public class LibMatrixFourier {
+
+ /**
Review Comment:
change indentation to tabs
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+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.*;
+import static org.junit.Assert.assertArrayEquals;
+
+public class FourierTest {
+
+ @Test
+ public void simple_test_one_dim() {
+ // 1st row real part, 2nd row imaginary part
+ double[][] in = {{0, 18, -15, 3},{0, 0, 0, 0}};
+ double[][] expected = {{6, 15, -36, 15},{0, -15, 0, 15}};
+
+ double[][] res = fft_one_dim(in);
+ for(double[] row : res){
+ for (double elem : row){
+ System.out.print(elem + " ");
+ }
+ System.out.println();
Review Comment:
As a general rule in testing (at least from my point of view) never print.
the only thing you print is if there is an error. and then you print only
the needed information to understand the error.
The reason is we run ~1 million tests 20+ times a day and we do not want to
read and store the output prints of tests from each execution. We only want to
know, what tests are failing, and why.
Therefore make the test construct a string, that you add as argument to the
assertArrayEquals such that if the check fails, it prints the information,
otherwise not.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+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.*;
+import static org.junit.Assert.assertArrayEquals;
+
+public class FourierTest {
+
+ @Test
+ public void simple_test_one_dim() {
+ // 1st row real part, 2nd row imaginary part
+ double[][] in = {{0, 18, -15, 3},{0, 0, 0, 0}};
+ double[][] expected = {{6, 15, -36, 15},{0, -15, 0, 15}};
+
+ double[][] res = fft_one_dim(in);
+ for(double[] row : res){
+ for (double elem : row){
+ System.out.print(elem + " ");
+ }
+ System.out.println();
+ }
+ assertArrayEquals(expected[0], res[0], 0.0001);
+ assertArrayEquals(expected[1], res[1], 0.0001);
+ }
+
+ @Test
+ public void simple_test_two_dim() {
+ // tested with numpy
+ double[][][] in = {{{0, 18},{-15, 3}},{{0, 0},{0, 0}}};
+
+ double[][][] expected = {{{6, -36},{30, 0}},{{0, 0},{0, 0}}};
+
+ double[][][] res = fft(in, false);
+
+ for(double[][] matrix : res){
+ for(double[] row : matrix) {
+ for (double elem : row) {
+ System.out.print(elem + " ");
+ }
+ System.out.println();
+ }
+ System.out.println();
+ }
+
+ for(int k = 0; k < 2 ; k++){
+ for(int i = 0; i < res[0].length; i++) {
+ assertArrayEquals(expected[k][i], res[k][i], 0.0001);
+ }
+ }
+ }
+
+ @Test
+ public void simple_test_one_dim_ifft() {
+
+ double[][] in = {{1, -2, 3, -4},{0, 0, 0, 0}};
+
+ double[][] res_fft = fft_one_dim(in);
+ double[][] res = ifft_one_dim(res_fft);
+
+ assertArrayEquals(in[0], res[0], 0.0001);
+ assertArrayEquals(in[1], res[1], 0.0001);
+ }
+
+ @Test
+ public void matrix_block_one_dim_test(){
+
+ double[] in = {0, 18, -15, 3};
+
+ double[] expected_re = {6,15,-36,15};
+ double[] expected_im = {0,-15,0,15};
+
+ MatrixBlock[] res = fft(in);
+ double[] res_re = res[0].getDenseBlockValues();
+ double[] res_im = res[1].getDenseBlockValues();
+
+ for(double elem : res_re){
+ System.out.print(elem+" ");
+ }
+ System.out.println();
+ for(double elem : res_im){
+ System.out.print(elem+" ");
+ }
+
+ assertArrayEquals(expected_re, res_re, 0.0001);
+ assertArrayEquals(expected_im, res_im, 0.0001);
+ }
+ @Test
Review Comment:
a new line between end of methods and new test declarations.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+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.*;
+import static org.junit.Assert.assertArrayEquals;
+
+public class FourierTest {
+
+ @Test
+ public void simple_test_one_dim() {
+ // 1st row real part, 2nd row imaginary part
+ double[][] in = {{0, 18, -15, 3},{0, 0, 0, 0}};
+ double[][] expected = {{6, 15, -36, 15},{0, -15, 0, 15}};
+
+ double[][] res = fft_one_dim(in);
+ for(double[] row : res){
+ for (double elem : row){
+ System.out.print(elem + " ");
+ }
+ System.out.println();
+ }
+ assertArrayEquals(expected[0], res[0], 0.0001);
+ assertArrayEquals(expected[1], res[1], 0.0001);
+ }
+
+ @Test
+ public void simple_test_two_dim() {
+ // tested with numpy
+ double[][][] in = {{{0, 18},{-15, 3}},{{0, 0},{0, 0}}};
+
+ double[][][] expected = {{{6, -36},{30, 0}},{{0, 0},{0, 0}}};
+
+ double[][][] res = fft(in, false);
+
+ for(double[][] matrix : res){
+ for(double[] row : matrix) {
+ for (double elem : row) {
+ System.out.print(elem + " ");
+ }
+ System.out.println();
+ }
+ System.out.println();
Review Comment:
remove prints.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTestWithFiles.java:
##########
@@ -0,0 +1,377 @@
+package org.apache.sysds.test.component.matrix;
Review Comment:
add licence
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTestData.py:
##########
@@ -0,0 +1,165 @@
+import numpy as np
+import csv
+import time
+
+def generate_inputs(num_inputs, max_power):
+ for _ in range(num_inputs):
+ power = np.random.randint(1, max_power+1)
+ length = 2 ** power
+ yield np.random.rand(length) # generate array of random floats
+
+def generate_complex_inputs(num_inputs, max_power):
+ for _ in range(num_inputs):
+ power = np.random.randint(1, max_power+1)
+ length = 2 ** power
+ real_part = np.random.rand(length)
+ imag_part = np.random.rand(length)
+ complex_array = real_part + 1j * imag_part
+ yield complex_array
+
+
+def compute_fft(inputs):
+ total_time = 0
+ num_calculations = 0
+
+ for input_array in inputs:
+ start_time = time.time()
+ result = np.fft.fft(input_array)
+ end_time = time.time()
+
+ total_time += end_time - start_time
+ num_calculations += 1
+
+ if num_calculations % 1000 == 0:
+ average_time = total_time / num_calculations
+ print(f"Average execution time after {num_calculations}
calculations: {average_time:.6f} seconds")
+
+ yield result
+
+def compute_ifft(inputs):
+ total_time = 0
+ num_calculations = 0
+
+ for input_array in inputs:
+ start_time = time.time()
+ result = np.fft.ifft(input_array)
+ end_time = time.time()
+
+ total_time += end_time - start_time
+ num_calculations += 1
+
+ if num_calculations % 1000 == 0:
+ average_time = total_time / num_calculations
+ print(f"Average execution time after {num_calculations}
calculations: {average_time:.6f} seconds")
+
+ yield result
+
+
+
+def save_to_file(inputs, outputs, filename, mode='a'):
+ with open(filename, mode, newline='') as file:
+ writer = csv.writer(file)
+ for input_array, output_array in zip(inputs, outputs):
+ flattened_data = np.concatenate((input_array, output_array.real,
output_array.imag))
+ writer.writerow(flattened_data)
+
+
+def save_to_file_complex(inputs, outputs, filename, mode='a'):
+ with open(filename, mode, newline='') as file:
+ writer = csv.writer(file)
+ for input_array, output_array in zip(inputs, outputs):
+ flattened_data = np.concatenate((input_array.real,
input_array.imag, output_array.real, output_array.imag))
+ writer.writerow(flattened_data)
+
+def generate_complex_inputs_2d(num_inputs, max_power):
+ for _ in range(num_inputs):
+ power = np.random.randint(1, max_power+1)
+ rows = 2 ** power
+ cols = 2 ** power
+ real_part = np.random.rand(rows, cols)
+ imag_part = np.random.rand(rows, cols)
+ complex_array = real_part + 1j * imag_part
+ yield complex_array
+
+def compute_ifft_2d(inputs):
+
+ total_time = 0
+ num_calculations = 0
+
+ for input_array in inputs:
+ start_time = time.time()
+ result = np.fft.ifft2(input_array)
+ end_time = time.time()
+
+ total_time += end_time - start_time
+ num_calculations += 1
+
+ if num_calculations % 1000 == 0:
+ average_time = total_time / num_calculations
+ print(f"Average execution time after {num_calculations}
calculations: {average_time:.6f} seconds")
+
+ yield result
+
+def save_to_file_complex_2d(inputs, outputs, filename, mode='a'):
+ with open(filename, mode, newline='') as file:
+ writer = csv.writer(file)
+ for input_array, output_array in zip(inputs, outputs):
+ flattened_input = np.concatenate((input_array.real.flatten(),
input_array.imag.flatten()))
+ flattened_output = np.concatenate((output_array.real.flatten(),
output_array.imag.flatten()))
+ writer.writerow(np.concatenate((flattened_input,
flattened_output)))
+
+
+# Parameters
+num_inputs = 100000
Review Comment:
add a
`if __name__ == "__main__":`
and indent the rest of the lines.
We do this to potentially allow other files to reuse some of your methods.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTestWithFiles.java:
##########
@@ -0,0 +1,377 @@
+package org.apache.sysds.test.component.matrix;
+
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.junit.Test;
+
+import java.io.BufferedReader;
+import java.io.FileReader;
+import java.io.IOException;
+
+import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.*;
+import static org.junit.Assert.assertTrue;
+
+public class FourierTestWithFiles {
+ int progressInterval = 5000;
+
+ // prior to executing the following tests it is necessary to run the Numpy
Script in FourierTestData.py
+ // and add the generated files to the root of the project.
+ @Test
+ public void testFftWithNumpyData() throws IOException {
+ String filename = "fft_data.csv"; // Path to your CSV file
+ BufferedReader reader = new BufferedReader(new FileReader(filename));
+ String line;
+ int lineNumber = 0;
+ long totalTime = 0; // Total time for all FFT computations
+ int numCalculations = 0; // Number of FFT computations
+
+ while ((line = reader.readLine()) != null) {
+ lineNumber++;
+
+ String[] values = line.split(",");
+ int n = values.length / 3;
+ double[][][] input = new double[2][1][n];
+ double[][] expected = new double[2][n]; // First row for real,
second row for imaginary parts
+
+ for (int i = 0; i < n; i++) {
+ input[0][0][i] = Double.parseDouble(values[i]);
+ expected[0][i] = Double.parseDouble(values[n + i]); // Real
part
+ expected[1][i] = Double.parseDouble(values[n * 2 + i]); //
Imaginary part
+ }
+
+ long startTime = System.nanoTime();
+ MatrixBlock[] actualBlocks = fft(input);
+ long endTime = System.nanoTime();
+
+ if(lineNumber > 1000){
+ totalTime += (endTime - startTime);
+ numCalculations++;
+
+ if (numCalculations % progressInterval == 0) {
+ double averageTime = (totalTime / 1e6) / numCalculations;
// Average time in milliseconds
+ System.out.println("fft(double[][][] in): Average
execution time after " + numCalculations + " calculations: " +
String.format("%.8f", averageTime/1000) + " s");
+ }
+ }
+
+ // Validate the FFT results
+ validateFftResults(expected, actualBlocks, lineNumber);
+ }
+
+ reader.close();
+
+ }
+
+ private void validateFftResults(double[][] expected, MatrixBlock[]
actualBlocks, int lineNumber) {
+ int length = expected[0].length;
+ for (int i = 0; i < length; i++) {
+ double realActual = actualBlocks[0].getValueDenseUnsafe(0, i);
+ double imagActual = actualBlocks[1].getValueDenseUnsafe(0, i);
+ assertEquals("Mismatch in real part at index " + i + " in line " +
lineNumber, expected[0][i], realActual, 1e-9);
+ assertEquals("Mismatch in imaginary part at index " + i + " in
line " + lineNumber, expected[1][i], imagActual, 1e-9);
+ }
+ if(lineNumber % progressInterval == 0){
+ System.out.println("fft(double[][][] in): Finished processing line
" + lineNumber);
+ }
+
+ }
+
+ @Test
+ public void testFftExecutionTime() throws IOException {
Review Comment:
instead of evaluating time here, add a Performance script in:
`src/test/java/org/apache/sysds/performance/micro`
this allows a bit more control for the experiment if you want it, instead of
relying on the test framework for measurements.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTest.java:
##########
@@ -0,0 +1,112 @@
+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.*;
+import static org.junit.Assert.assertArrayEquals;
+
+public class FourierTest {
+
+ @Test
+ public void simple_test_one_dim() {
+ // 1st row real part, 2nd row imaginary part
+ double[][] in = {{0, 18, -15, 3},{0, 0, 0, 0}};
+ double[][] expected = {{6, 15, -36, 15},{0, -15, 0, 15}};
+
+ double[][] res = fft_one_dim(in);
+ for(double[] row : res){
+ for (double elem : row){
+ System.out.print(elem + " ");
+ }
+ System.out.println();
+ }
+ assertArrayEquals(expected[0], res[0], 0.0001);
+ assertArrayEquals(expected[1], res[1], 0.0001);
+ }
+
+ @Test
+ public void simple_test_two_dim() {
+ // tested with numpy
+ double[][][] in = {{{0, 18},{-15, 3}},{{0, 0},{0, 0}}};
+
+ double[][][] expected = {{{6, -36},{30, 0}},{{0, 0},{0, 0}}};
+
+ double[][][] res = fft(in, false);
+
+ for(double[][] matrix : res){
+ for(double[] row : matrix) {
+ for (double elem : row) {
+ System.out.print(elem + " ");
+ }
+ System.out.println();
+ }
+ System.out.println();
+ }
+
+ for(int k = 0; k < 2 ; k++){
+ for(int i = 0; i < res[0].length; i++) {
+ assertArrayEquals(expected[k][i], res[k][i], 0.0001);
+ }
+ }
+ }
+
+ @Test
+ public void simple_test_one_dim_ifft() {
+
+ double[][] in = {{1, -2, 3, -4},{0, 0, 0, 0}};
+
+ double[][] res_fft = fft_one_dim(in);
+ double[][] res = ifft_one_dim(res_fft);
+
+ assertArrayEquals(in[0], res[0], 0.0001);
+ assertArrayEquals(in[1], res[1], 0.0001);
+ }
+
+ @Test
+ public void matrix_block_one_dim_test(){
+
+ double[] in = {0, 18, -15, 3};
+
+ double[] expected_re = {6,15,-36,15};
+ double[] expected_im = {0,-15,0,15};
+
+ MatrixBlock[] res = fft(in);
+ double[] res_re = res[0].getDenseBlockValues();
+ double[] res_im = res[1].getDenseBlockValues();
+
+ for(double elem : res_re){
+ System.out.print(elem+" ");
+ }
+ System.out.println();
+ for(double elem : res_im){
+ System.out.print(elem+" ");
+ }
+
+ assertArrayEquals(expected_re, res_re, 0.0001);
+ assertArrayEquals(expected_im, res_im, 0.0001);
+ }
+ @Test
+ public void matrix_block_two_dim_test(){
+
+ double[][][] in = {{{0, 18},{-15, 3}}};
+
+ double[] flattened_expected_re = {6,-36, 30,0};
+ double[] flattened_expected_im = {0,0,0,0};
+
+ MatrixBlock[] res = fft(in);
+ double[] res_re = res[0].getDenseBlockValues();
+ double[] res_im = res[1].getDenseBlockValues();
+
+ for(double elem : res_re){
+ System.out.print(elem+" ");
+ }
+ System.out.println();
+ for(double elem : res_im){
+ System.out.print(elem+" ");
+ }
+
+ assertArrayEquals(flattened_expected_re, res_re, 0.0001);
+ assertArrayEquals(flattened_expected_im, res_im, 0.0001);
+ }
+
+}
Review Comment:
add a blank newline in the end of files.
##########
src/test/java/org/apache/sysds/test/component/matrix/FourierTestWithFiles.java:
##########
@@ -0,0 +1,377 @@
+package org.apache.sysds.test.component.matrix;
+
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.junit.Test;
+
+import java.io.BufferedReader;
+import java.io.FileReader;
+import java.io.IOException;
+
+import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.*;
+import static org.junit.Assert.assertTrue;
+
+public class FourierTestWithFiles {
+ int progressInterval = 5000;
+
+ // prior to executing the following tests it is necessary to run the Numpy
Script in FourierTestData.py
Review Comment:
indent with tabs.
--
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
To unsubscribe, e-mail: dev-unsubscr...@systemds.apache.org
For queries about this service, please contact Infrastructure at:
us...@infra.apache.org
