diff --git a/src/main/java/org/apache/sysds/runtime/util/EMAUtils.java b/src/main/java/org/apache/sysds/runtime/util/EMAUtils.java index f6fd1d14294..2e544ee32c9 100644 --- a/src/main/java/org/apache/sysds/runtime/util/EMAUtils.java +++ b/src/main/java/org/apache/sysds/runtime/util/EMAUtils.java @@ -72,7 +72,7 @@ public double coef() { public class EMAUtils { public static FrameBlock exponentialMovingAverageImputation(FrameBlock block, int search_iterations, String mode, - int freq, Double alpha, Double beta, Double gamma) { + int freq, Double alpha, Double beta, Double gamma) { int cols = block.getNumColumns(); int rows = block.getNumRows(); @@ -145,7 +145,32 @@ public Container(Double[] vals, double error) { double rsme; } + public static void validateRange(Double... values) { + if (values == null || values.length == 0){ + throw new IllegalArgumentException("No arguments provided"); + } + for (Double value : values) { + if (value > 1 || value < 0) { + throw new IllegalArgumentException( + "Input is out of range [0,1]" + ); + } + } + + } + + public static void validateFrequency(int freq, Double[] data){ + if (freq <= 0){ + throw new IllegalArgumentException("Value too low"); + } else if (2*freq >= data.length) { + throw new ArrayIndexOutOfBoundsException("Frequency too high"); + } + + } public static Container single_exponential_smoothing(Double[] data, Double alpha) { + + validateRange(alpha); + int n = data.length; Double[] pred = new Double[n]; @@ -186,12 +211,26 @@ public static Container single_exponential_smoothing(Double[] data, Double alpha } public static Container double_exponential_smoothing(Double[] data, Double alpha, Double beta) { + + // Exception handling + validateRange(alpha, beta); + + if (data == null || data.length < 3) { + throw new IllegalArgumentException("At least three observations are needed to calculate RMSE"); + } + + if (data[0] == null || data[1] == null || Double.isNaN(data[0]) || Double.isNaN(data[1])) { + throw new IllegalArgumentException("The first two observations must not be missing"); + } + + // Initializing variables int n = data.length; ArrayList pred = new ArrayList<>(n-1); Double[] s = new Double[n-1]; Double[] b = new Double[n-1]; + // Setting up s_0, b_0, y_0 s[0] = data[1]; b[0] = data[1] - data[0]; pred.add(s[0] + b[0]); @@ -202,7 +241,9 @@ public static Container double_exponential_smoothing(Double[] data, Double alpha ArrayList not_missing_pred = new ArrayList<>(); int n_size = 0; + // Iterative updating of the model for (int i = 1; i < n-1; i++) { + // Imputes missing data if (Double.isNaN(data[i + 1])) { val = pred.get(i - 1); } else { @@ -211,11 +252,15 @@ public static Container double_exponential_smoothing(Double[] data, Double alpha s[i] = alpha * val + (1 - alpha) * (s[i-1] + b[i-1]); b[i] = beta * (s[i] - s[i-1]) + (1 - beta) * b[i-1]; - pred.add(s[i] + b[i]); - } + if (i < n - 2){ + pred.add(s[i] + b[i]); + } + } + pred.add(0, data[1]); pred.add(0, data[0]); + // Initializing the data for calculation of RMSE for (int i = 0; i < data.length; i++) { if (!Double.isNaN(data[i])) { not_missing.add(data[i]); @@ -224,6 +269,7 @@ public static Container double_exponential_smoothing(Double[] data, Double alpha } } + // Calculation of RMSE double sum = .0; for (int i = 0; i < not_missing.size(); i++) { sum += Math.pow(not_missing.get(i) - not_missing_pred.get(i), 2); @@ -235,7 +281,13 @@ public static Container double_exponential_smoothing(Double[] data, Double alpha } public static Container triple_exponential_smoothing(Double[] data, Double alpha, Double beta, - Double gamma, Integer freq) { + Double gamma, Integer freq) { + + // Exception handling + validateRange(alpha, beta, gamma); + validateFrequency(freq, data); + + //Initialization data double l = freq * 2; ArrayList start_data = new ArrayList<>(); @@ -243,11 +295,15 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha start_data.add(data[i]); } + // Building a centered moving average for the inititialization of the level variable ArrayList filt = new ArrayList<>(); ArrayList trend = new ArrayList<>(); double len = freq; + // To set up the centered moving average filter, there is two cases, odd and even frequency. + // For the odd frequency, there is a true center point, so the filter weights are the same. + // For the even frequency, there is no center point, therefore, the edge weights are halved. if (freq % 2 == 0) { len = freq - 1; } @@ -261,9 +317,11 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha filt.add(0.5/freq); } + // Applying the CMA filter to the trend double trend_len = l - filt.size() + 1; - for (int i = 0; i < l - trend_len; i++) { + // Goes through the start data, and sums the average trend over the first to seasons. + for (int i = 0; i < trend_len; i++) { double sum = 0; for (int j = i; j < i + filt.size(); j++) { sum += start_data.get(j) * filt.get(j - i); @@ -272,6 +330,7 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha trend.add(sum); } + // Removing the trend, to reveal the seasonal pattern int cut = (int) (l - trend.size()) / 2; ArrayList season_tmp = new ArrayList<>(); @@ -282,6 +341,7 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha Double[] season = new Double[freq]; + // Average by season position for (int i = 0; i < freq; i++) { double combined = 0; if (i + freq < trend.size()) { @@ -293,6 +353,7 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha season[(i + (freq / 2)) % freq] = combined; } + // Preparing the data for linear regression, in order to initialize level double sum = 0; for (int i = 0; i < season.length; i++) { sum += season[i]; @@ -316,6 +377,8 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha int n = data.length; + + // Initalize arrays and create first forecast double[] s = new double[n - freq]; s[0] = linreg.intercept(); @@ -338,6 +401,7 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha int n_size = 0; + // Update the model per data point for (int i = 1; i < n - freq; i++) { if (Double.isNaN(data[i + freq - 1] )) { val = pred.get(i - 1); @@ -356,6 +420,7 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha pred.add(i, data[i]); } + // Initialize data for the RMSE calculation for (int i = 0; i < data.length; i++) { if (!Double.isNaN(data[i])) { not_missing.add(data[i]); @@ -364,6 +429,7 @@ public static Container triple_exponential_smoothing(Double[] data, Double alpha } } + // Calculate the RMSE sum = .0; for (int i = 0; i < not_missing.size(); i++) { sum += Math.pow(not_missing.get(i) - not_missing_pred.get(i), 2); diff --git a/src/test/java/org/apache/sysds/runtime/util/EMAUtilsTest.java b/src/test/java/org/apache/sysds/runtime/util/EMAUtilsTest.java new file mode 100644 index 00000000000..e0284ffe9f6 --- /dev/null +++ b/src/test/java/org/apache/sysds/runtime/util/EMAUtilsTest.java @@ -0,0 +1,151 @@ +/* + * 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.util; + +import org.apache.sysds.runtime.matrix.data.MatrixValue; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.junit.Assert; +import org.junit.Test; + +import java.util.HashMap; + +import static org.apache.sysds.test.TestUtils.readRMatrixFromFS; +import static org.apache.sysds.runtime.util.EMAUtils.*; + +public class EMAUtilsTest extends AutomatedTestBase { + private static final String TEST_NAME = "exponentialSmoothing"; + private static final String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = TEST_DIR + EMAUtilsTest.class.getSimpleName() + "/"; + + private static final double REFERENCE_EPS = 1e-4; + private static final Double[] TS_DATA = { + 112., 118., 132., 129., 121., 135., 148., 148., 136., 119., 104., 118., + 115., 126., 141., 135., 125., 149., 170., 170., 158., 133., 114., 140., + 145., 150., 178., 163., 172., 178., 199., 199., 184., 162., 146., 166., + 171., 180., 193., 181., 183., 218., 230., 242., 209., 191., 172., 194., + 196., 196., 236., 235., 229., 243., 264., 272., 237., 211., 180., 201. + }; + @Override + public void setUp() { + addTestConfiguration(TEST_NAME, + new TestConfiguration(TEST_CLASS_DIR, TEST_NAME, new String[] { "R_double", "R_triple", "R_triple_odd" })); + } + + @Test + public void testDoubleExponentialSmoothing() { + fullRScriptName = SCRIPT_DIR + TEST_DIR + TEST_NAME + ".R"; + String resultPath = output("R_double"); + + rCmd = getRCmd("double", resultPath); + runRScript(true); + + HashMap rOutput = readRMatrixFromFS(resultPath); + Assert.assertFalse("R output should contain values", rOutput.isEmpty()); + + int numRows = (int) rOutput.keySet().stream().mapToLong(idx -> idx.row).max().orElse(0); + double[] rFittedValues = new double[numRows]; + for (int i = 1; i <= numRows; i++) { + rFittedValues[i - 1] = rOutput.get(new MatrixValue.CellIndex(i, 1)); + } + + + Container javaFittedValues = double_exponential_smoothing(TS_DATA, 0.3, 0.1); + + Assert.assertEquals("Array length mismatch", rFittedValues.length, javaFittedValues.values.length-2); + for (int i = 0; i < rFittedValues.length; i++) { + Assert.assertEquals("Mismatch at index " + i, rFittedValues[i], javaFittedValues.values[i+2], REFERENCE_EPS); + } + } + + @Test + public void testTripleExponentialSmoothing(){ + int freq = 12; + fullRScriptName = SCRIPT_DIR + TEST_DIR + TEST_NAME + ".R"; + String resultPath = output("R_triple"); + + rCmd = getRCmd("triple", resultPath); + runRScript(true); + + HashMap rOutputTriple = readRMatrixFromFS(resultPath); + Assert.assertFalse( + "R output should contain values", + rOutputTriple.isEmpty() + ); + + int numRows = (int) rOutputTriple.keySet().stream() + .mapToLong(idx -> idx.row) + .max() + .orElse(0); + + double[] rFittedValues = new double[numRows]; + + for (int i = 0; i < numRows; i++) { + rFittedValues[i] = rOutputTriple.getOrDefault( + new MatrixValue.CellIndex(i + 1, 1), 0.0); + } + + Container javaResult = triple_exponential_smoothing( + TS_DATA, 0.3, 0.1, 0.2, freq); + + Assert.assertEquals("Array length mismatch", rFittedValues.length, javaResult.values.length - freq); + + for (int i = 0; i < rFittedValues.length; i++) { + Assert.assertEquals("Mismatch at index " + i, rFittedValues[i], javaResult.values[i + freq], REFERENCE_EPS); + } + } + + @Test + public void testOddTripleExponentialSmoothing(){ + int freq = 7; + fullRScriptName = SCRIPT_DIR + TEST_DIR + TEST_NAME + ".R"; + String resultPath = output("R_triple_odd"); + + rCmd = getRCmd("odd_triple", resultPath); + runRScript(true); + + HashMap rOutputTriple = readRMatrixFromFS(resultPath); + Assert.assertFalse( + "R output should contain values", + rOutputTriple.isEmpty() + ); + + int numRows = (int) rOutputTriple.keySet().stream() + .mapToLong(idx -> idx.row) + .max() + .orElse(0); + + double[] rFittedValues = new double[numRows]; + + for (int i = 0; i < numRows; i++) { + rFittedValues[i] = rOutputTriple.getOrDefault( + new MatrixValue.CellIndex(i + 1, 1), 0.0); + } + + Container javaResult = triple_exponential_smoothing( + TS_DATA, 0.3, 0.1, 0.2, freq); + + Assert.assertEquals("Array length mismatch", rFittedValues.length, javaResult.values.length - freq); + + for (int i = 0; i < rFittedValues.length; i++) { + Assert.assertEquals("Mismatch at index " + i, rFittedValues[i], javaResult.values[i + freq], REFERENCE_EPS); + } + } +} \ No newline at end of file diff --git a/src/test/scripts/functions/builtin/exponentialSmoothing.R b/src/test/scripts/functions/builtin/exponentialSmoothing.R new file mode 100644 index 00000000000..50c30e29121 --- /dev/null +++ b/src/test/scripts/functions/builtin/exponentialSmoothing.R @@ -0,0 +1,95 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +args <- commandArgs(trailingOnly = TRUE) + +if (length(args) != 2L) { + stop("Expected arguments: ") +} + +mode <- args[1] +output_path <- args[2] + +ts_data <- c( + 112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118, + 115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140, + 145, 150, 178, 163, 172, 178, 199, 199, 184, 162, 146, 166, + 171, 180, 193, 181, 183, 218, 230, 242, 209, 191, 172, 194, + 196, 196, 236, 235, 229, 243, 264, 272, 237, 211, 180, 201 +) + +if (mode == "double") { + hw_model <- stats::HoltWinters( + ts_data, + alpha=0.3, + beta = 0.1, + gamma = FALSE) + fit_vals <- as.numeric(fitted(hw_model)[, "xhat"]) + + res_matrix <- matrix(0, nrow = length(fit_vals), ncol = 1) + res_matrix[, 1] <- fit_vals + +} else if (mode == "triple") { + ts_obj <- ts(ts_data, frequency = 12) + + hw_model <- stats::HoltWinters( + ts_obj, + alpha = 0.3, + beta = 0.1, + gamma = 0.2, + seasonal = "additive" + ) + + fit_vals <- as.numeric(fitted(hw_model)[, "xhat"]) + + res_matrix <- matrix(0, nrow = length(fit_vals), ncol = 1) + res_matrix[, 1] <- fit_vals + +} else if (mode == "odd_triple") { + ts_obj <- ts(ts_data, frequency = 7) + + hw_model <- stats::HoltWinters( + ts_obj, + alpha = 0.3, + beta = 0.1, + gamma = 0.2, + seasonal = "additive" + ) + + fit_vals <- as.numeric(fitted(hw_model)[, "xhat"]) + + res_matrix <- matrix(0, nrow = length(fit_vals), ncol = 1) + res_matrix[, 1] <- fit_vals + + } else { + stop(paste("Unknown smoothing mode:", mode)) +} + +dir.create( + dirname(output_path), + recursive = TRUE, + showWarnings = FALSE +) + +Matrix::writeMM( + Matrix::Matrix(res_matrix, sparse = TRUE), + file = output_path +) \ No newline at end of file