gaturchenko commented on code in PR #2499: URL: https://github.com/apache/systemds/pull/2499#discussion_r3704570644
########## scripts/builtin/powerTransform.dml: ########## @@ -0,0 +1,335 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- Review Comment: Please, add a line break after the license ########## scripts/builtin/powerTransformApply.dml: ########## @@ -0,0 +1,124 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- Review Comment: Please, add a line break after the license ########## scripts/builtin/powerTransform.dml: ########## @@ -0,0 +1,335 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- +# Power transformation using the selected method. +# Reduces feature skewness by estimating and applying an optimal transformation parameter for each column. +# +# INPUT: +# ------------------------------------------------------------------------------------- +# X Input feature matrix of shape n-by-m +# method Power transformation method: "yeo-johnson" (default) or "box-cox" +# standardize Whether to normalize transformed columns to zero mean and unit variance +# ------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------- +# Y Power-transformed matrix of shape n-by-m +# lambdas Estimated lambda parameters of shape 1-by-m, one per column +# means Transformed column means of shape 1-by-m, or an empty matrix when not standardized +# scales Transformed column scales of shape 1-by-m, or an empty matrix when not standardized +# ------------------------------------------------------------------------------------- + +m_powerTransform = function( + Matrix[Double] X, + String method="yeo-johnson", + Boolean standardize=TRUE) + return ( + Matrix[Double] Y, + Matrix[Double] lambdas, + Matrix[Double] means, + Matrix[Double] scales) +{ + if (method != "yeo-johnson" & method != "box-cox") { + stop("powerTransform: unsupported method '" + method + + "'; expected 'yeo-johnson' or 'box-cox'") + } + + if (method == "box-cox" & min(X) <= 0.0) { + stop("powerTransform: Box-Cox requires strictly positive input") + } + + n = nrow(X) + m = ncol(X) + lambdas = matrix(1.0, rows=1, cols=m) # Initialize first, then replace each column with the best lambdas + + # Estimate lambda for each column separately + for (j in 1:m){ + x = X[,j] + + # Yeo-Johnson leaves constant columns unchanged; Box-Cox rejects them + if (max(x) == min(x)) { + if (method == "yeo-johnson") { + lambdas[1,j] = 1.0; + } + else { + stop("powerTransform: Box-Cox does not support constant columns") + } + } + else{ + lambdas[1,j] = ptEstimateLambda(x, method); + } + } + + # Apply the fitted transformation before optional standardization + emptyStats = matrix(0.0, rows=0, cols=0) + Y = powerTransformApply(X, lambdas, emptyStats, emptyStats, method); + + means = matrix(0.0, rows=0, cols=0) + scales = matrix(0.0, rows=0, cols=0) + + if (standardize) { + means = colMeans(Y) + Y = Y - means + scales = sqrt(colSums(Y^2) / n) + scales = replace(target=scales, pattern=NaN, replacement=1.0) + scales = replace(target=scales, pattern=0.0, replacement=1.0) + Y = Y / scales + } +} +ptEstimateLambda = function(Matrix[Double] x, String method) + return (Double lambda) +{ + lower = -2.0; + upper = 2.0; + + lambda = ptBrentSearch(x, lower, upper, method); +} + +# Compute negative log likelihood; lower lambda score is better + +ptNegLogLikelihood = function( + Matrix[Double] x, + Double lambda, + String method) + return (Double negLogLikelihood) +{ + # powerTransformApply needs lambda as a matrix + lambdaMatrix = matrix(lambda, rows=1, cols=1); + + # Apply one transform with this lambda and use the temporary y for scoring + emptyStats = matrix(0.0, rows=0, cols=0) + y = powerTransformApply(x, lambdaMatrix, emptyStats, emptyStats, method); + + # Safety check; give a huge penalty score when variance is below 0 + n = nrow(x); + yMean = mean(y); + yVariance = sum((y - yMean)^2) / n; + + if (yVariance <= 0.0) { + negLogLikelihood = 1e300 + } + + # Start scoring after the safety check passes + # The objective has two parts + + else { + # Variance Term + logLikelihood = -n / 2.0 * log(yVariance); + + # Jacobian term for the selected transformation + if (method == "box-cox") { + jacobian = (lambda - 1.0) * sum(log(x)); + } + else { + jacobian = (lambda - 1.0) * sum(sign(x) * log(abs(x) + 1.0)); Review Comment: Same logic as for the comment above ########## scripts/builtin/powerTransform.dml: ########## @@ -0,0 +1,335 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- +# Power transformation using the selected method. +# Reduces feature skewness by estimating and applying an optimal transformation parameter for each column. +# +# INPUT: +# ------------------------------------------------------------------------------------- +# X Input feature matrix of shape n-by-m +# method Power transformation method: "yeo-johnson" (default) or "box-cox" +# standardize Whether to normalize transformed columns to zero mean and unit variance +# ------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------- +# Y Power-transformed matrix of shape n-by-m +# lambdas Estimated lambda parameters of shape 1-by-m, one per column +# means Transformed column means of shape 1-by-m, or an empty matrix when not standardized +# scales Transformed column scales of shape 1-by-m, or an empty matrix when not standardized +# ------------------------------------------------------------------------------------- + +m_powerTransform = function( + Matrix[Double] X, + String method="yeo-johnson", + Boolean standardize=TRUE) + return ( + Matrix[Double] Y, + Matrix[Double] lambdas, + Matrix[Double] means, + Matrix[Double] scales) +{ + if (method != "yeo-johnson" & method != "box-cox") { + stop("powerTransform: unsupported method '" + method + + "'; expected 'yeo-johnson' or 'box-cox'") + } + + if (method == "box-cox" & min(X) <= 0.0) { + stop("powerTransform: Box-Cox requires strictly positive input") + } + + n = nrow(X) + m = ncol(X) + lambdas = matrix(1.0, rows=1, cols=m) # Initialize first, then replace each column with the best lambdas + + # Estimate lambda for each column separately + for (j in 1:m){ + x = X[,j] + + # Yeo-Johnson leaves constant columns unchanged; Box-Cox rejects them + if (max(x) == min(x)) { + if (method == "yeo-johnson") { + lambdas[1,j] = 1.0; + } + else { + stop("powerTransform: Box-Cox does not support constant columns") + } + } + else{ + lambdas[1,j] = ptEstimateLambda(x, method); + } + } + + # Apply the fitted transformation before optional standardization + emptyStats = matrix(0.0, rows=0, cols=0) + Y = powerTransformApply(X, lambdas, emptyStats, emptyStats, method); + + means = matrix(0.0, rows=0, cols=0) + scales = matrix(0.0, rows=0, cols=0) + + if (standardize) { + means = colMeans(Y) + Y = Y - means + scales = sqrt(colSums(Y^2) / n) + scales = replace(target=scales, pattern=NaN, replacement=1.0) + scales = replace(target=scales, pattern=0.0, replacement=1.0) + Y = Y / scales + } +} +ptEstimateLambda = function(Matrix[Double] x, String method) + return (Double lambda) +{ + lower = -2.0; + upper = 2.0; + + lambda = ptBrentSearch(x, lower, upper, method); +} + +# Compute negative log likelihood; lower lambda score is better + +ptNegLogLikelihood = function( + Matrix[Double] x, + Double lambda, + String method) + return (Double negLogLikelihood) +{ + # powerTransformApply needs lambda as a matrix + lambdaMatrix = matrix(lambda, rows=1, cols=1); + + # Apply one transform with this lambda and use the temporary y for scoring + emptyStats = matrix(0.0, rows=0, cols=0) + y = powerTransformApply(x, lambdaMatrix, emptyStats, emptyStats, method); + + # Safety check; give a huge penalty score when variance is below 0 + n = nrow(x); + yMean = mean(y); + yVariance = sum((y - yMean)^2) / n; + + if (yVariance <= 0.0) { + negLogLikelihood = 1e300 + } + + # Start scoring after the safety check passes + # The objective has two parts + + else { + # Variance Term + logLikelihood = -n / 2.0 * log(yVariance); + + # Jacobian term for the selected transformation + if (method == "box-cox") { + jacobian = (lambda - 1.0) * sum(log(x)); + } + else { + jacobian = (lambda - 1.0) * sum(sign(x) * log(abs(x) + 1.0)); + } + + # Combine + logLikelihood = logLikelihood + jacobian; + # Return the negative value + negLogLikelihood = -logLikelihood; + } +} + +# Minimize the negative log likelihood with Brent optimization +ptBrentSearch = function( + Matrix[Double] x, + Double lower, + Double upper, + String method) + return (Double lambdaOptimal) +{ + # Expand the initial interval until it brackets a minimum + goldenRatio = 1.618034; + maxBracketIterations = 1000; + lowerScore = ptNegLogLikelihood(x=x, lambda=lower, method=method); + upperScore = ptNegLogLikelihood(x=x, lambda=upper, method=method); + + if (lowerScore < upperScore) { + xa = upper; + fa = upperScore; + xb = lower; + fb = lowerScore; + } + else { + xa = lower; + fa = lowerScore; + xb = upper; + fb = upperScore; + } + + initialXc = xb + goldenRatio * (xb - xa); + initialFc = ptNegLogLikelihood(x=x, lambda=initialXc, method=method); + xc = initialXc; + fc = initialFc; + bracketIteration = 0; + while ((fc < fb) & (bracketIteration < maxBracketIterations)) { + nextXc = xc + goldenRatio * (xc - xb); + nextFc = ptNegLogLikelihood(x=x, lambda=nextXc, method=method); + xa = xb; + fa = fb; + xb = xc; + fb = fc; + xc = nextXc; + fc = nextFc; + bracketIteration = bracketIteration + 1; + } + + if ((bracketIteration >= maxBracketIterations) & (fc < fb)) { + stop("powerTransform: failed to bracket lambda minimum") + } + + validBracket = ((fb < fa) & (fb <= fc)) | ((fb <= fa) & (fb < fc)); + if (!validBracket) { + stop("powerTransform: failed to bracket lambda minimum") + } + + a = min(xa, xc); + b = max(xa, xc); + + # Brent combines inverse parabolic interpolation with a golden-section fallback + goldenMean = 0.3819660112501051; + sqrtEpsilon = sqrt(2.2e-16); + tolerance = 1.48e-8; + maxIterations = 500; + + xf = a + goldenMean * (b - a); + nfc = xf; + fulc = xf; + fx = ptNegLogLikelihood(x, xf, method); + fnfc = fx; + ffulc = fx; + + rat = 0.0; + e = 0.0; + midpoint = 0.5 * (a + b); + tol1 = sqrtEpsilon * abs(xf) + tolerance / 3.0; + tol2 = 2.0 * tol1; + + iteration = 0; + while ((abs(xf - midpoint) > (tol2 - 0.5 * (b - a))) & + (iteration < maxIterations)) { + goldenStep = TRUE; + + # Try an inverse parabolic step when the three tracked points support it + if (abs(e) > tol1) { + goldenStep = FALSE; + r = (xf - nfc) * (fx - ffulc); + q = (xf - fulc) * (fx - fnfc); + p = (xf - fulc) * q - (xf - nfc) * r; + q = 2.0 * (q - r); + + if (q > 0.0) { + p = -p; + } + + q = abs(q); + previousE = e; + e = rat; + + # Reject interpolation outside the bracket or larger than the previous step + if ((q > 0.0) & (abs(p) < abs(0.5 * q * previousE)) & + (p > q * (a - xf)) & (p < q * (b - xf))) { + rat = p / q; + candidate = xf + rat; + + if (((candidate - a) < tol2) | ((b - candidate) < tol2)) { + if (midpoint >= xf) { + rat = tol1; + } + else { + rat = -tol1; + } + } + } + else { + goldenStep = TRUE; + } + } + + # Fall back to a golden-section step when interpolation is unreliable + if (goldenStep) { + if (xf >= midpoint) { + e = a - xf; + } + else { + e = b - xf; + } + rat = goldenMean * e; + } + + if (rat >= 0.0) { + candidate = xf + max(abs(rat), tol1); + } + else { + candidate = xf - max(abs(rat), tol1); + } + + fCandidate = ptNegLogLikelihood(x, candidate, method); + + # Update the bracket and retain the three best distinct points + if (fCandidate <= fx) { + if (candidate >= xf) { + a = xf; + } + else { + b = xf; + } + + fulc = nfc; + ffulc = fnfc; + nfc = xf; + fnfc = fx; + xf = candidate; + fx = fCandidate; + } + else { + if (candidate < xf) { + a = candidate; + } + else { + b = candidate; + } + + if ((fCandidate <= fnfc) | (nfc == xf)) { + fulc = nfc; + ffulc = fnfc; + nfc = candidate; + fnfc = fCandidate; + } + else if ((fCandidate <= ffulc) | (fulc == xf) | (fulc == nfc)) { + fulc = candidate; + ffulc = fCandidate; + } + } + + midpoint = 0.5 * (a + b); + tol1 = sqrtEpsilon * abs(xf) + tolerance / 3.0; + tol2 = 2.0 * tol1; + iteration = iteration + 1; + } + + if ((iteration >= maxIterations) & + (abs(xf - midpoint) > (tol2 - 0.5 * (b - a)))) { + stop("powerTransform: Brent optimization failed to converge") Review Comment: Same `NaN` handling issue ########## scripts/builtin/powerTransformApply.dml: ########## @@ -0,0 +1,124 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- +# Applies a fitted power transformation and optional standardization. +# Transforms each feature using its previously estimated lambda and scaling parameters. +# +# INPUT: +# ------------------------------------------------------------------------------------- +# X Input feature matrix of shape n-by-m +# lambdas Precomputed lambda parameters of shape 1-by-m, one per column +# means Transformed column means of shape 1-by-m; empty to skip standardization +# scales Transformed column scales of shape 1-by-m; empty to skip standardization +# method Power transformation method: "yeo-johnson" (default) or "box-cox" +# ------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------- +# Y Power-transformed matrix of shape n-by-m +# ------------------------------------------------------------------------------------- + + +m_powerTransformApply = function( + Matrix[Double] X, + Matrix[Double] lambdas, + Matrix[Double] means, + Matrix[Double] scales, + String method="yeo-johnson") + return (Matrix[Double] Y) +{ + n = nrow(X) + m = ncol(X) + + if (method != "yeo-johnson" & method != "box-cox") { + stop("powerTransformApply: unsupported method '" + method + + "'; expected 'yeo-johnson' or 'box-cox'") + } + + if (method == "box-cox" & min(X) <= 0.0) { Review Comment: This runs for every call from `powerTransform.dml`, so you are re-validating what you already have ########## scripts/builtin/powerTransform.dml: ########## @@ -0,0 +1,335 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- +# Power transformation using the selected method. +# Reduces feature skewness by estimating and applying an optimal transformation parameter for each column. +# +# INPUT: +# ------------------------------------------------------------------------------------- +# X Input feature matrix of shape n-by-m +# method Power transformation method: "yeo-johnson" (default) or "box-cox" +# standardize Whether to normalize transformed columns to zero mean and unit variance +# ------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------- +# Y Power-transformed matrix of shape n-by-m +# lambdas Estimated lambda parameters of shape 1-by-m, one per column +# means Transformed column means of shape 1-by-m, or an empty matrix when not standardized +# scales Transformed column scales of shape 1-by-m, or an empty matrix when not standardized +# ------------------------------------------------------------------------------------- + +m_powerTransform = function( + Matrix[Double] X, + String method="yeo-johnson", + Boolean standardize=TRUE) + return ( + Matrix[Double] Y, + Matrix[Double] lambdas, + Matrix[Double] means, + Matrix[Double] scales) +{ + if (method != "yeo-johnson" & method != "box-cox") { + stop("powerTransform: unsupported method '" + method + + "'; expected 'yeo-johnson' or 'box-cox'") + } + + if (method == "box-cox" & min(X) <= 0.0) { + stop("powerTransform: Box-Cox requires strictly positive input") + } + + n = nrow(X) + m = ncol(X) + lambdas = matrix(1.0, rows=1, cols=m) # Initialize first, then replace each column with the best lambdas + + # Estimate lambda for each column separately + for (j in 1:m){ + x = X[,j] + + # Yeo-Johnson leaves constant columns unchanged; Box-Cox rejects them + if (max(x) == min(x)) { + if (method == "yeo-johnson") { + lambdas[1,j] = 1.0; + } + else { + stop("powerTransform: Box-Cox does not support constant columns") + } + } + else{ + lambdas[1,j] = ptEstimateLambda(x, method); + } + } + + # Apply the fitted transformation before optional standardization + emptyStats = matrix(0.0, rows=0, cols=0) + Y = powerTransformApply(X, lambdas, emptyStats, emptyStats, method); + + means = matrix(0.0, rows=0, cols=0) + scales = matrix(0.0, rows=0, cols=0) + + if (standardize) { + means = colMeans(Y) + Y = Y - means + scales = sqrt(colSums(Y^2) / n) + scales = replace(target=scales, pattern=NaN, replacement=1.0) + scales = replace(target=scales, pattern=0.0, replacement=1.0) + Y = Y / scales + } +} +ptEstimateLambda = function(Matrix[Double] x, String method) + return (Double lambda) +{ + lower = -2.0; + upper = 2.0; + + lambda = ptBrentSearch(x, lower, upper, method); +} + +# Compute negative log likelihood; lower lambda score is better + +ptNegLogLikelihood = function( + Matrix[Double] x, + Double lambda, + String method) + return (Double negLogLikelihood) +{ + # powerTransformApply needs lambda as a matrix + lambdaMatrix = matrix(lambda, rows=1, cols=1); + + # Apply one transform with this lambda and use the temporary y for scoring + emptyStats = matrix(0.0, rows=0, cols=0) + y = powerTransformApply(x, lambdaMatrix, emptyStats, emptyStats, method); + + # Safety check; give a huge penalty score when variance is below 0 + n = nrow(x); + yMean = mean(y); + yVariance = sum((y - yMean)^2) / n; + + if (yVariance <= 0.0) { + negLogLikelihood = 1e300 + } + + # Start scoring after the safety check passes + # The objective has two parts + + else { + # Variance Term + logLikelihood = -n / 2.0 * log(yVariance); + + # Jacobian term for the selected transformation + if (method == "box-cox") { + jacobian = (lambda - 1.0) * sum(log(x)); Review Comment: Only `lambda` changes here, `x` is always the same, so you can compute `sum(log(x))` just once, which can be done with a separate function. Then you can add something like `jacTerm` as another input to `ptNegLogLikelihood` and get `jacobian = (lambda - 1.0) * jacTerm;`. This should improve performance substantially for large matrices ########## scripts/builtin/powerTransform.dml: ########## @@ -0,0 +1,335 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- +# Power transformation using the selected method. +# Reduces feature skewness by estimating and applying an optimal transformation parameter for each column. +# +# INPUT: +# ------------------------------------------------------------------------------------- +# X Input feature matrix of shape n-by-m +# method Power transformation method: "yeo-johnson" (default) or "box-cox" +# standardize Whether to normalize transformed columns to zero mean and unit variance +# ------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------- +# Y Power-transformed matrix of shape n-by-m +# lambdas Estimated lambda parameters of shape 1-by-m, one per column +# means Transformed column means of shape 1-by-m, or an empty matrix when not standardized +# scales Transformed column scales of shape 1-by-m, or an empty matrix when not standardized +# ------------------------------------------------------------------------------------- + +m_powerTransform = function( + Matrix[Double] X, + String method="yeo-johnson", + Boolean standardize=TRUE) + return ( + Matrix[Double] Y, + Matrix[Double] lambdas, + Matrix[Double] means, + Matrix[Double] scales) +{ + if (method != "yeo-johnson" & method != "box-cox") { + stop("powerTransform: unsupported method '" + method + + "'; expected 'yeo-johnson' or 'box-cox'") + } + + if (method == "box-cox" & min(X) <= 0.0) { + stop("powerTransform: Box-Cox requires strictly positive input") + } + + n = nrow(X) + m = ncol(X) + lambdas = matrix(1.0, rows=1, cols=m) # Initialize first, then replace each column with the best lambdas + + # Estimate lambda for each column separately + for (j in 1:m){ + x = X[,j] + + # Yeo-Johnson leaves constant columns unchanged; Box-Cox rejects them + if (max(x) == min(x)) { + if (method == "yeo-johnson") { + lambdas[1,j] = 1.0; + } + else { + stop("powerTransform: Box-Cox does not support constant columns") + } + } + else{ + lambdas[1,j] = ptEstimateLambda(x, method); + } + } + + # Apply the fitted transformation before optional standardization + emptyStats = matrix(0.0, rows=0, cols=0) + Y = powerTransformApply(X, lambdas, emptyStats, emptyStats, method); + + means = matrix(0.0, rows=0, cols=0) + scales = matrix(0.0, rows=0, cols=0) + + if (standardize) { + means = colMeans(Y) + Y = Y - means + scales = sqrt(colSums(Y^2) / n) + scales = replace(target=scales, pattern=NaN, replacement=1.0) + scales = replace(target=scales, pattern=0.0, replacement=1.0) + Y = Y / scales + } +} +ptEstimateLambda = function(Matrix[Double] x, String method) + return (Double lambda) +{ + lower = -2.0; + upper = 2.0; + + lambda = ptBrentSearch(x, lower, upper, method); +} + +# Compute negative log likelihood; lower lambda score is better + +ptNegLogLikelihood = function( + Matrix[Double] x, + Double lambda, + String method) + return (Double negLogLikelihood) +{ + # powerTransformApply needs lambda as a matrix + lambdaMatrix = matrix(lambda, rows=1, cols=1); + + # Apply one transform with this lambda and use the temporary y for scoring + emptyStats = matrix(0.0, rows=0, cols=0) + y = powerTransformApply(x, lambdaMatrix, emptyStats, emptyStats, method); Review Comment: This has a fixed invocation overhead, and for matrices with, say, 1000 rows and 4 columns, this will dominate the runtime, so it should be inlined somehow ########## scripts/builtin/powerTransform.dml: ########## @@ -0,0 +1,335 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- +# Power transformation using the selected method. +# Reduces feature skewness by estimating and applying an optimal transformation parameter for each column. +# +# INPUT: +# ------------------------------------------------------------------------------------- +# X Input feature matrix of shape n-by-m +# method Power transformation method: "yeo-johnson" (default) or "box-cox" +# standardize Whether to normalize transformed columns to zero mean and unit variance +# ------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------- +# Y Power-transformed matrix of shape n-by-m +# lambdas Estimated lambda parameters of shape 1-by-m, one per column +# means Transformed column means of shape 1-by-m, or an empty matrix when not standardized +# scales Transformed column scales of shape 1-by-m, or an empty matrix when not standardized +# ------------------------------------------------------------------------------------- + +m_powerTransform = function( + Matrix[Double] X, + String method="yeo-johnson", + Boolean standardize=TRUE) + return ( + Matrix[Double] Y, + Matrix[Double] lambdas, + Matrix[Double] means, + Matrix[Double] scales) +{ + if (method != "yeo-johnson" & method != "box-cox") { + stop("powerTransform: unsupported method '" + method + + "'; expected 'yeo-johnson' or 'box-cox'") + } + + if (method == "box-cox" & min(X) <= 0.0) { + stop("powerTransform: Box-Cox requires strictly positive input") + } + + n = nrow(X) + m = ncol(X) + lambdas = matrix(1.0, rows=1, cols=m) # Initialize first, then replace each column with the best lambdas + + # Estimate lambda for each column separately + for (j in 1:m){ + x = X[,j] + + # Yeo-Johnson leaves constant columns unchanged; Box-Cox rejects them + if (max(x) == min(x)) { + if (method == "yeo-johnson") { + lambdas[1,j] = 1.0; + } + else { + stop("powerTransform: Box-Cox does not support constant columns") + } + } + else{ + lambdas[1,j] = ptEstimateLambda(x, method); + } + } + + # Apply the fitted transformation before optional standardization + emptyStats = matrix(0.0, rows=0, cols=0) + Y = powerTransformApply(X, lambdas, emptyStats, emptyStats, method); + + means = matrix(0.0, rows=0, cols=0) + scales = matrix(0.0, rows=0, cols=0) + + if (standardize) { + means = colMeans(Y) + Y = Y - means + scales = sqrt(colSums(Y^2) / n) + scales = replace(target=scales, pattern=NaN, replacement=1.0) + scales = replace(target=scales, pattern=0.0, replacement=1.0) + Y = Y / scales + } +} +ptEstimateLambda = function(Matrix[Double] x, String method) + return (Double lambda) +{ + lower = -2.0; + upper = 2.0; + + lambda = ptBrentSearch(x, lower, upper, method); +} + +# Compute negative log likelihood; lower lambda score is better + +ptNegLogLikelihood = function( + Matrix[Double] x, + Double lambda, + String method) + return (Double negLogLikelihood) +{ + # powerTransformApply needs lambda as a matrix + lambdaMatrix = matrix(lambda, rows=1, cols=1); + + # Apply one transform with this lambda and use the temporary y for scoring + emptyStats = matrix(0.0, rows=0, cols=0) + y = powerTransformApply(x, lambdaMatrix, emptyStats, emptyStats, method); + + # Safety check; give a huge penalty score when variance is below 0 + n = nrow(x); + yMean = mean(y); + yVariance = sum((y - yMean)^2) / n; + + if (yVariance <= 0.0) { + negLogLikelihood = 1e300 + } + + # Start scoring after the safety check passes + # The objective has two parts + + else { + # Variance Term + logLikelihood = -n / 2.0 * log(yVariance); + + # Jacobian term for the selected transformation + if (method == "box-cox") { + jacobian = (lambda - 1.0) * sum(log(x)); + } + else { + jacobian = (lambda - 1.0) * sum(sign(x) * log(abs(x) + 1.0)); + } + + # Combine + logLikelihood = logLikelihood + jacobian; + # Return the negative value + negLogLikelihood = -logLikelihood; + } +} + +# Minimize the negative log likelihood with Brent optimization +ptBrentSearch = function( + Matrix[Double] x, + Double lower, + Double upper, + String method) + return (Double lambdaOptimal) +{ + # Expand the initial interval until it brackets a minimum + goldenRatio = 1.618034; + maxBracketIterations = 1000; + lowerScore = ptNegLogLikelihood(x=x, lambda=lower, method=method); + upperScore = ptNegLogLikelihood(x=x, lambda=upper, method=method); + + if (lowerScore < upperScore) { + xa = upper; + fa = upperScore; + xb = lower; + fb = lowerScore; + } + else { + xa = lower; + fa = lowerScore; + xb = upper; + fb = upperScore; + } + + initialXc = xb + goldenRatio * (xb - xa); + initialFc = ptNegLogLikelihood(x=x, lambda=initialXc, method=method); + xc = initialXc; + fc = initialFc; + bracketIteration = 0; + while ((fc < fb) & (bracketIteration < maxBracketIterations)) { + nextXc = xc + goldenRatio * (xc - xb); + nextFc = ptNegLogLikelihood(x=x, lambda=nextXc, method=method); + xa = xb; + fa = fb; + xb = xc; + fb = fc; + xc = nextXc; + fc = nextFc; + bracketIteration = bracketIteration + 1; + } + + if ((bracketIteration >= maxBracketIterations) & (fc < fb)) { + stop("powerTransform: failed to bracket lambda minimum") Review Comment: IMO this is too conservative, a single missing value in a column will abort the entire execution. Can we not just fall back to the best `lambda` round and just explain how we handle `NaN`s same as `scale.dml` does? ########## src/test/java/org/apache/sysds/test/functions/builtin/part2/BuiltinPowerTransformTest.java: ########## @@ -0,0 +1,271 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part2; + +import java.util.HashMap; + +import org.junit.Assert; +import org.junit.Test; + +import org.apache.sysds.common.Types.ExecMode; +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.runtime.DMLScriptException; +import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; + +public class BuiltinPowerTransformTest extends AutomatedTestBase { + private static final String TRANSFORM_TEST_NAME = "powerTransform"; + private static final String APPLY_TEST_NAME = "powerTransformApply"; + private static final String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = + TEST_DIR + BuiltinPowerTransformTest.class.getSimpleName() + "/"; + + private static final double TRANSFORM_EPS = 1e-6; + private static final double APPLY_EPS = 1e-9; + + @Override + public void setUp() { + addTestConfiguration( + TRANSFORM_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + TRANSFORM_TEST_NAME, + new String[] {"Y", "L", "S"} + ) + ); + addTestConfiguration( + APPLY_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + APPLY_TEST_NAME, + new String[] {"Y"} + ) + ); + } + + @Test + public void testPowerTransformYeoJohnsonDefaultDenseCP() { + double[][] input = { + {-2, 1, 5}, + {-1, 1, 5}, + { 0, 2, 5}, + { 1, 3, 5}, + { 2, 6, 5}, + { 4, 12, 5} + }; + runPowerTransformTest("default", true, input, false); + } + + @Test + public void testPowerTransformBoxCoxUnstandardizedDenseCP() { + double[][] input = { + { 0.5, 1}, + { 1.0, 2}, + { 2.0, 3}, + { 4.0, 5}, + { 8.0, 9}, + {16.0, 17} + }; + runPowerTransformTest("box-cox", false, input, false); + } + + @Test + public void testPowerTransformYeoJohnsonLambdaAboveInitialInterval() { + double[][] input = { + {0.00}, + {0.97}, + {0.98}, + {0.99}, + {1.00} + }; + runPowerTransformTest("yeo-johnson", false, input, false); + assertLambdaOutsideInitialInterval(true); + } + + @Test + public void testPowerTransformBoxCoxLambdaBelowInitialInterval() { + double[][] input = { + {1.00}, + {1.01}, + {1.02}, + {1.03}, + {10.0} + }; + runPowerTransformTest("box-cox", false, input, false); + assertLambdaOutsideInitialInterval(false); + } + + @Test + public void testPowerTransformBoxCoxRejectsNonPositiveInput() { + double[][] input = { + {0, 1}, + {1, 2} + }; + runPowerTransformTest("box-cox", false, input, true); + } + + @Test + public void testPowerTransformApplyYeoJohnsonDenseCP() { + double[][] input = { + {-2, -2, -2}, + {-1, -1, -1}, + { 0, 0, 0}, + { 1, 1, 1}, + { 2, 2, 2} + }; + runPowerTransformApplyTest(ExecType.CP, "yeo-johnson", true, input, false); + } + + @Test + public void testPowerTransformApplyBoxCoxDenseCP() { + double[][] input = { + {0.5, 0.5, 0.5}, + {1.0, 1.0, 1.0}, + {2.0, 2.0, 2.0}, + {4.0, 4.0, 4.0}, + {8.0, 8.0, 8.0} + }; + runPowerTransformApplyTest(ExecType.CP, "box-cox", false, input, false); + } + + @Test + public void testPowerTransformApplyBoxCoxRejectsNonPositiveInput() { + double[][] input = { + {0, 1, 2}, + {1, 2, 3} + }; + runPowerTransformApplyTest(ExecType.CP, "box-cox", false, input, true); + } + + private void runPowerTransformTest( + String method, boolean standardize, double[][] input, boolean shouldFail) { + ExecMode oldExecMode = setExecMode(ExecType.CP); + + try { + loadTestConfiguration(getTestConfiguration(TRANSFORM_TEST_NAME)); + + String home = SCRIPT_DIR + TEST_DIR; + fullDMLScriptName = home + TRANSFORM_TEST_NAME + ".dml"; + fullRScriptName = home + TRANSFORM_TEST_NAME + ".R"; + programArgs = new String[] { + "-exec", "singlenode", Review Comment: These are redundant ########## scripts/builtin/powerTransform.dml: ########## @@ -0,0 +1,335 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- +# Power transformation using the selected method. +# Reduces feature skewness by estimating and applying an optimal transformation parameter for each column. +# +# INPUT: +# ------------------------------------------------------------------------------------- +# X Input feature matrix of shape n-by-m +# method Power transformation method: "yeo-johnson" (default) or "box-cox" +# standardize Whether to normalize transformed columns to zero mean and unit variance +# ------------------------------------------------------------------------------------- +# +# OUTPUT: +# ------------------------------------------------------------------------------------- +# Y Power-transformed matrix of shape n-by-m +# lambdas Estimated lambda parameters of shape 1-by-m, one per column +# means Transformed column means of shape 1-by-m, or an empty matrix when not standardized +# scales Transformed column scales of shape 1-by-m, or an empty matrix when not standardized +# ------------------------------------------------------------------------------------- + +m_powerTransform = function( + Matrix[Double] X, + String method="yeo-johnson", + Boolean standardize=TRUE) + return ( + Matrix[Double] Y, + Matrix[Double] lambdas, + Matrix[Double] means, + Matrix[Double] scales) +{ + if (method != "yeo-johnson" & method != "box-cox") { + stop("powerTransform: unsupported method '" + method + + "'; expected 'yeo-johnson' or 'box-cox'") + } + + if (method == "box-cox" & min(X) <= 0.0) { + stop("powerTransform: Box-Cox requires strictly positive input") + } + + n = nrow(X) + m = ncol(X) + lambdas = matrix(1.0, rows=1, cols=m) # Initialize first, then replace each column with the best lambdas + + # Estimate lambda for each column separately + for (j in 1:m){ + x = X[,j] + + # Yeo-Johnson leaves constant columns unchanged; Box-Cox rejects them + if (max(x) == min(x)) { + if (method == "yeo-johnson") { + lambdas[1,j] = 1.0; + } + else { + stop("powerTransform: Box-Cox does not support constant columns") + } + } + else{ + lambdas[1,j] = ptEstimateLambda(x, method); + } + } + + # Apply the fitted transformation before optional standardization + emptyStats = matrix(0.0, rows=0, cols=0) + Y = powerTransformApply(X, lambdas, emptyStats, emptyStats, method); + + means = matrix(0.0, rows=0, cols=0) + scales = matrix(0.0, rows=0, cols=0) + + if (standardize) { + means = colMeans(Y) + Y = Y - means + scales = sqrt(colSums(Y^2) / n) + scales = replace(target=scales, pattern=NaN, replacement=1.0) + scales = replace(target=scales, pattern=0.0, replacement=1.0) + Y = Y / scales + } +} +ptEstimateLambda = function(Matrix[Double] x, String method) + return (Double lambda) +{ + lower = -2.0; + upper = 2.0; + + lambda = ptBrentSearch(x, lower, upper, method); +} + +# Compute negative log likelihood; lower lambda score is better + +ptNegLogLikelihood = function( + Matrix[Double] x, + Double lambda, + String method) + return (Double negLogLikelihood) +{ + # powerTransformApply needs lambda as a matrix + lambdaMatrix = matrix(lambda, rows=1, cols=1); + + # Apply one transform with this lambda and use the temporary y for scoring + emptyStats = matrix(0.0, rows=0, cols=0) + y = powerTransformApply(x, lambdaMatrix, emptyStats, emptyStats, method); + + # Safety check; give a huge penalty score when variance is below 0 + n = nrow(x); + yMean = mean(y); + yVariance = sum((y - yMean)^2) / n; + + if (yVariance <= 0.0) { + negLogLikelihood = 1e300 + } + + # Start scoring after the safety check passes + # The objective has two parts + + else { + # Variance Term + logLikelihood = -n / 2.0 * log(yVariance); + + # Jacobian term for the selected transformation + if (method == "box-cox") { + jacobian = (lambda - 1.0) * sum(log(x)); + } + else { + jacobian = (lambda - 1.0) * sum(sign(x) * log(abs(x) + 1.0)); + } + + # Combine + logLikelihood = logLikelihood + jacobian; + # Return the negative value + negLogLikelihood = -logLikelihood; + } +} + +# Minimize the negative log likelihood with Brent optimization +ptBrentSearch = function( + Matrix[Double] x, + Double lower, + Double upper, + String method) + return (Double lambdaOptimal) +{ + # Expand the initial interval until it brackets a minimum + goldenRatio = 1.618034; + maxBracketIterations = 1000; + lowerScore = ptNegLogLikelihood(x=x, lambda=lower, method=method); + upperScore = ptNegLogLikelihood(x=x, lambda=upper, method=method); + + if (lowerScore < upperScore) { + xa = upper; + fa = upperScore; + xb = lower; + fb = lowerScore; + } + else { + xa = lower; + fa = lowerScore; + xb = upper; + fb = upperScore; + } + + initialXc = xb + goldenRatio * (xb - xa); + initialFc = ptNegLogLikelihood(x=x, lambda=initialXc, method=method); + xc = initialXc; + fc = initialFc; + bracketIteration = 0; + while ((fc < fb) & (bracketIteration < maxBracketIterations)) { + nextXc = xc + goldenRatio * (xc - xb); + nextFc = ptNegLogLikelihood(x=x, lambda=nextXc, method=method); + xa = xb; + fa = fb; + xb = xc; + fb = fc; + xc = nextXc; + fc = nextFc; + bracketIteration = bracketIteration + 1; + } + + if ((bracketIteration >= maxBracketIterations) & (fc < fb)) { + stop("powerTransform: failed to bracket lambda minimum") + } + + validBracket = ((fb < fa) & (fb <= fc)) | ((fb <= fa) & (fb < fc)); + if (!validBracket) { + stop("powerTransform: failed to bracket lambda minimum") Review Comment: Same `NaN` handling issue ########## src/test/java/org/apache/sysds/test/functions/builtin/part2/BuiltinPowerTransformTest.java: ########## @@ -0,0 +1,271 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part2; + +import java.util.HashMap; + +import org.junit.Assert; +import org.junit.Test; + +import org.apache.sysds.common.Types.ExecMode; +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.runtime.DMLScriptException; +import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; + +public class BuiltinPowerTransformTest extends AutomatedTestBase { + private static final String TRANSFORM_TEST_NAME = "powerTransform"; + private static final String APPLY_TEST_NAME = "powerTransformApply"; + private static final String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = + TEST_DIR + BuiltinPowerTransformTest.class.getSimpleName() + "/"; + + private static final double TRANSFORM_EPS = 1e-6; + private static final double APPLY_EPS = 1e-9; + + @Override + public void setUp() { + addTestConfiguration( + TRANSFORM_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + TRANSFORM_TEST_NAME, + new String[] {"Y", "L", "S"} + ) + ); + addTestConfiguration( + APPLY_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + APPLY_TEST_NAME, + new String[] {"Y"} + ) + ); + } + + @Test + public void testPowerTransformYeoJohnsonDefaultDenseCP() { + double[][] input = { + {-2, 1, 5}, + {-1, 1, 5}, + { 0, 2, 5}, + { 1, 3, 5}, + { 2, 6, 5}, + { 4, 12, 5} + }; + runPowerTransformTest("default", true, input, false); + } + + @Test + public void testPowerTransformBoxCoxUnstandardizedDenseCP() { + double[][] input = { + { 0.5, 1}, + { 1.0, 2}, + { 2.0, 3}, + { 4.0, 5}, + { 8.0, 9}, + {16.0, 17} + }; + runPowerTransformTest("box-cox", false, input, false); + } + + @Test + public void testPowerTransformYeoJohnsonLambdaAboveInitialInterval() { + double[][] input = { + {0.00}, + {0.97}, + {0.98}, + {0.99}, + {1.00} + }; + runPowerTransformTest("yeo-johnson", false, input, false); + assertLambdaOutsideInitialInterval(true); + } + + @Test + public void testPowerTransformBoxCoxLambdaBelowInitialInterval() { + double[][] input = { + {1.00}, + {1.01}, + {1.02}, + {1.03}, + {10.0} + }; + runPowerTransformTest("box-cox", false, input, false); + assertLambdaOutsideInitialInterval(false); + } + + @Test + public void testPowerTransformBoxCoxRejectsNonPositiveInput() { + double[][] input = { + {0, 1}, + {1, 2} + }; + runPowerTransformTest("box-cox", false, input, true); + } + + @Test + public void testPowerTransformApplyYeoJohnsonDenseCP() { + double[][] input = { + {-2, -2, -2}, + {-1, -1, -1}, + { 0, 0, 0}, + { 1, 1, 1}, + { 2, 2, 2} + }; + runPowerTransformApplyTest(ExecType.CP, "yeo-johnson", true, input, false); + } + + @Test + public void testPowerTransformApplyBoxCoxDenseCP() { + double[][] input = { + {0.5, 0.5, 0.5}, + {1.0, 1.0, 1.0}, + {2.0, 2.0, 2.0}, + {4.0, 4.0, 4.0}, + {8.0, 8.0, 8.0} + }; + runPowerTransformApplyTest(ExecType.CP, "box-cox", false, input, false); + } + + @Test + public void testPowerTransformApplyBoxCoxRejectsNonPositiveInput() { + double[][] input = { + {0, 1, 2}, + {1, 2, 3} + }; + runPowerTransformApplyTest(ExecType.CP, "box-cox", false, input, true); + } + + private void runPowerTransformTest( + String method, boolean standardize, double[][] input, boolean shouldFail) { + ExecMode oldExecMode = setExecMode(ExecType.CP); + + try { + loadTestConfiguration(getTestConfiguration(TRANSFORM_TEST_NAME)); + + String home = SCRIPT_DIR + TEST_DIR; + fullDMLScriptName = home + TRANSFORM_TEST_NAME + ".dml"; + fullRScriptName = home + TRANSFORM_TEST_NAME + ".R"; + programArgs = new String[] { + "-exec", "singlenode", + "-args", + input("X"), + output("Y"), + output("L"), + output("S"), + method, + Boolean.toString(standardize) + }; + + String referenceMethod = method.equals("default") ? "yeo-johnson" : method; + rCmd = getRCmd( + inputDir(), + expectedDir(), + referenceMethod, + Boolean.toString(standardize) + ); + + writeInputMatrixWithMTD("X", input, true); + runTest(true, shouldFail, shouldFail ? DMLScriptException.class : null, -1); + if (shouldFail) + return; + + runRScript(true); + compareOutput("Y", TRANSFORM_EPS); + compareOutput("L", TRANSFORM_EPS); + compareOutput("S", TRANSFORM_EPS); + } + catch (Exception exception) { + throw new RuntimeException(exception); + } + finally { + resetExecMode(oldExecMode); + } + } + + private void runPowerTransformApplyTest( + ExecType execType, String method, boolean standardize, double[][] input, boolean shouldFail) { + ExecMode oldExecMode = setExecMode(execType); + + try { + loadTestConfiguration(getTestConfiguration(APPLY_TEST_NAME)); + + String home = SCRIPT_DIR + TEST_DIR; + fullDMLScriptName = home + APPLY_TEST_NAME + ".dml"; + fullRScriptName = home + APPLY_TEST_NAME + ".R"; + programArgs = new String[] { + "-exec", "singlenode", Review Comment: These are redundant ########## src/test/scripts/functions/builtin/powerTransform.R: ########## Review Comment: Implementing `powerTransform` like this defeats the purpose of R tests, as we want to compare against an actual existing R implementation. A better way of doing this is to actually find an R package that implements `powerTransform` with your 2 methods and compare against it ########## src/test/java/org/apache/sysds/test/functions/builtin/part2/BuiltinPowerTransformTest.java: ########## @@ -0,0 +1,271 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part2; + +import java.util.HashMap; + +import org.junit.Assert; +import org.junit.Test; + +import org.apache.sysds.common.Types.ExecMode; +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.runtime.DMLScriptException; +import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; + +public class BuiltinPowerTransformTest extends AutomatedTestBase { + private static final String TRANSFORM_TEST_NAME = "powerTransform"; + private static final String APPLY_TEST_NAME = "powerTransformApply"; + private static final String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = + TEST_DIR + BuiltinPowerTransformTest.class.getSimpleName() + "/"; + + private static final double TRANSFORM_EPS = 1e-6; + private static final double APPLY_EPS = 1e-9; + + @Override + public void setUp() { + addTestConfiguration( + TRANSFORM_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + TRANSFORM_TEST_NAME, + new String[] {"Y", "L", "S"} + ) + ); + addTestConfiguration( + APPLY_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + APPLY_TEST_NAME, + new String[] {"Y"} + ) + ); + } + + @Test + public void testPowerTransformYeoJohnsonDefaultDenseCP() { + double[][] input = { + {-2, 1, 5}, + {-1, 1, 5}, + { 0, 2, 5}, + { 1, 3, 5}, + { 2, 6, 5}, + { 4, 12, 5} + }; + runPowerTransformTest("default", true, input, false); + } + + @Test + public void testPowerTransformBoxCoxUnstandardizedDenseCP() { + double[][] input = { + { 0.5, 1}, + { 1.0, 2}, + { 2.0, 3}, + { 4.0, 5}, + { 8.0, 9}, + {16.0, 17} + }; + runPowerTransformTest("box-cox", false, input, false); + } + + @Test + public void testPowerTransformYeoJohnsonLambdaAboveInitialInterval() { Review Comment: This test fails both in CI and for me locally ########## src/test/scripts/functions/builtin/powerTransformApply.R: ########## Review Comment: Same as for the other R file with the only exception that, if an R package does not implement `transformApply`, there is no reason to run an R test ########## src/test/java/org/apache/sysds/test/functions/builtin/part2/BuiltinPowerTransformTest.java: ########## @@ -0,0 +1,271 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part2; + +import java.util.HashMap; + +import org.junit.Assert; +import org.junit.Test; + +import org.apache.sysds.common.Types.ExecMode; +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.runtime.DMLScriptException; +import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.apache.sysds.test.TestUtils; + +public class BuiltinPowerTransformTest extends AutomatedTestBase { + private static final String TRANSFORM_TEST_NAME = "powerTransform"; + private static final String APPLY_TEST_NAME = "powerTransformApply"; + private static final String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = + TEST_DIR + BuiltinPowerTransformTest.class.getSimpleName() + "/"; + + private static final double TRANSFORM_EPS = 1e-6; + private static final double APPLY_EPS = 1e-9; + + @Override + public void setUp() { + addTestConfiguration( + TRANSFORM_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + TRANSFORM_TEST_NAME, + new String[] {"Y", "L", "S"} + ) + ); + addTestConfiguration( + APPLY_TEST_NAME, + new TestConfiguration( + TEST_CLASS_DIR, + APPLY_TEST_NAME, + new String[] {"Y"} + ) + ); + } + + @Test + public void testPowerTransformYeoJohnsonDefaultDenseCP() { + double[][] input = { + {-2, 1, 5}, + {-1, 1, 5}, + { 0, 2, 5}, + { 1, 3, 5}, + { 2, 6, 5}, + { 4, 12, 5} + }; + runPowerTransformTest("default", true, input, false); + } + + @Test + public void testPowerTransformBoxCoxUnstandardizedDenseCP() { + double[][] input = { + { 0.5, 1}, + { 1.0, 2}, + { 2.0, 3}, + { 4.0, 5}, + { 8.0, 9}, + {16.0, 17} + }; + runPowerTransformTest("box-cox", false, input, false); + } + + @Test + public void testPowerTransformYeoJohnsonLambdaAboveInitialInterval() { + double[][] input = { + {0.00}, + {0.97}, + {0.98}, + {0.99}, + {1.00} + }; + runPowerTransformTest("yeo-johnson", false, input, false); + assertLambdaOutsideInitialInterval(true); + } + + @Test + public void testPowerTransformBoxCoxLambdaBelowInitialInterval() { + double[][] input = { + {1.00}, + {1.01}, + {1.02}, + {1.03}, + {10.0} + }; + runPowerTransformTest("box-cox", false, input, false); + assertLambdaOutsideInitialInterval(false); + } + + @Test + public void testPowerTransformBoxCoxRejectsNonPositiveInput() { + double[][] input = { + {0, 1}, + {1, 2} + }; + runPowerTransformTest("box-cox", false, input, true); + } + + @Test + public void testPowerTransformApplyYeoJohnsonDenseCP() { + double[][] input = { + {-2, -2, -2}, + {-1, -1, -1}, + { 0, 0, 0}, + { 1, 1, 1}, + { 2, 2, 2} + }; + runPowerTransformApplyTest(ExecType.CP, "yeo-johnson", true, input, false); + } + + @Test + public void testPowerTransformApplyBoxCoxDenseCP() { + double[][] input = { + {0.5, 0.5, 0.5}, + {1.0, 1.0, 1.0}, + {2.0, 2.0, 2.0}, + {4.0, 4.0, 4.0}, + {8.0, 8.0, 8.0} + }; + runPowerTransformApplyTest(ExecType.CP, "box-cox", false, input, false); + } + + @Test + public void testPowerTransformApplyBoxCoxRejectsNonPositiveInput() { + double[][] input = { + {0, 1, 2}, + {1, 2, 3} + }; + runPowerTransformApplyTest(ExecType.CP, "box-cox", false, input, true); + } + + private void runPowerTransformTest( + String method, boolean standardize, double[][] input, boolean shouldFail) { + ExecMode oldExecMode = setExecMode(ExecType.CP); Review Comment: You only add CP coverage, would be good to have Spark covered as well. Let us know if there are particular issues with Spark for Power Transform however -- 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: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
