diff --git a/src/main/java/org/apache/sysds/runtime/instructions/cp/QuantilePickCPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/cp/QuantilePickCPInstruction.java index a7bfbf5a16c..12bfc249ffe 100644 --- a/src/main/java/org/apache/sysds/runtime/instructions/cp/QuantilePickCPInstruction.java +++ b/src/main/java/org/apache/sysds/runtime/instructions/cp/QuantilePickCPInstruction.java @@ -90,16 +90,12 @@ public void processInstruction(ExecutionContext ec) { if ( input2.getDataType() == DataType.SCALAR ) { ScalarObject quantile = ec.getScalarInput(input2); - //pick value w/ explicit averaging for even-length arrays - double picked = matBlock.pickValue( - quantile.getDoubleValue(), matBlock.getLength()%2==0); + double picked = matBlock.pickValue(quantile.getDoubleValue()); ec.setScalarOutput(output.getName(), new DoubleObject(picked)); - } + } else { MatrixBlock quantiles = ec.getMatrixInput(input2.getName()); - //pick value w/ explicit averaging for even-length arrays - MatrixBlock resultBlock = matBlock.pickValues( - quantiles, new MatrixBlock(), matBlock.getLength()%2==0); + MatrixBlock resultBlock = matBlock.pickValues(quantiles, new MatrixBlock()); quantiles = null; ec.releaseMatrixInput(input2.getName()); ec.setMatrixOutput(output.getName(), resultBlock); diff --git a/src/main/java/org/apache/sysds/runtime/instructions/fed/QuantilePickFEDInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/fed/QuantilePickFEDInstruction.java index d16e1e16ff1..6d1367e472a 100644 --- a/src/main/java/org/apache/sysds/runtime/instructions/fed/QuantilePickFEDInstruction.java +++ b/src/main/java/org/apache/sysds/runtime/instructions/fed/QuantilePickFEDInstruction.java @@ -6,9 +6,9 @@ * 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 @@ -26,8 +26,8 @@ import java.util.List; import java.util.Map; import java.util.Set; +import java.util.TreeSet; import java.util.stream.Collectors; -import java.util.stream.Stream; import org.apache.commons.lang3.ArrayUtils; import org.apache.commons.lang3.tuple.ImmutablePair; @@ -63,6 +63,8 @@ @SuppressWarnings("unchecked") public class QuantilePickFEDInstruction extends BinaryFEDInstruction { + private static final int NUM_BUCKETS = 256; + private final OperationTypes _type; public QuantilePickFEDInstruction(Operator op, CPOperand in, CPOperand out, OperationTypes type, boolean inmem, @@ -85,7 +87,7 @@ public static QuantilePickFEDInstruction parseInstruction(QuantilePickCPInstruct return new QuantilePickFEDInstruction(instr.getOperator(), instr.input1, instr.input2, instr.output, instr.getOperationType(), instr.isInMem(), instr.getOpcode(), instr.getInstructionString()); } - + public static QuantilePickFEDInstruction parseInstruction(QuantilePickSPInstruction instr) { return new QuantilePickFEDInstruction(instr.getOperator(), instr.input1, instr.input2, instr.output, instr.getOperationType(), false, instr.getOpcode(), instr.getInstructionString()); @@ -136,7 +138,7 @@ public void processInstruction(ExecutionContext ec) { processRowQPick(ec); } - public MatrixBlock getEquiHeightBins(ExecutionContext ec, int colID, double[] quantiles) { + public MatrixBlock getEquiHeightBins(ExecutionContext ec, int colID, double[] quantiles) { FrameObject inFrame = ec.getFrameObject(input1); FederationMap frameFedMap = inFrame.getFedMapping(); @@ -207,43 +209,54 @@ public MatrixBlock getEquiHeightBins(ExecutionContext ec, int colID, double[ }); // Find weights sum, min and max - double globalMin = Double.MAX_VALUE, globalMax = Double.MIN_VALUE, vectorLength = inFrame.getNumColumns() == 2 ? 0 : inFrame.getNumRows(); + double globalMin = Double.MAX_VALUE, globalMax = Double.MIN_VALUE; + int vectorLength = inFrame.getNumColumns() == 2 ? 0 : (int) inFrame.getNumRows(); for(double[] values : minMax) { globalMin = Math.min(globalMin, values[0]); globalMax = Math.max(globalMax, values[1]); } - // If multiple quantiles take first histogram and reuse bins, otherwise recursively get bin with result - int numBuckets = 256; // (int) Math.round(in.getNumRows() / 2.0); - - T ret = createHistogram(in, (int) vectorLength, globalMin, globalMax, numBuckets, -1, false); - - // Compute and set results - MatrixBlock quantileValues = computeMultipleQuantiles(ec, in, (int[]) ret, quantiles, (int) vectorLength, varID, (globalMax-globalMin) / numBuckets, globalMin, _type, true); + // Equi-height bin boundaries use plain ceil-based ranks (no R type 7 interpolation). + final double N = vectorLength; + int[] ranks = Arrays.stream(quantiles).mapToInt(q -> (int) Math.round(N * q)).toArray(); + Map rankToValue = pickMultipleRanks(in, ranks, vectorLength, varID, globalMin, globalMax); ec.removeVariable(String.valueOf(varID)); - // Add min to the result - MatrixBlock res = new MatrixBlock(quantileValues.getNumRows() + 1, 1, false); - res.set(0,0, globalMin); - res.copy(1, quantileValues.getNumRows(), 0, 0, quantileValues,false); - + // Result: [globalMin, v_r1, v_r2, ...] as a column vector. + MatrixBlock res = new MatrixBlock(quantiles.length + 1, 1, false); + res.set(0, 0, globalMin); + for(int i = 0; i < quantiles.length; i++) + res.set(i + 1, 0, rankToValue.get(ranks[i])); return res; } - public void processRowQPick(ExecutionContext ec) { + public void processRowQPick(ExecutionContext ec) { MatrixObject in = ec.getMatrixObject(input1); FederationMap fedMap = in.getFedMapping(); - boolean average = _type == OperationTypes.MEDIAN || _type == OperationTypes.VALUEPICK; - double[] quantiles = input2 != null ? (input2.isMatrix() ? ec.getMatrixInput(input2).getDenseBlockValues() : - input2.isScalar() ? new double[] {ec.getScalarInput(input2).getDoubleValue()} : null) : - (average ? new double[] {0.5} : _type == OperationTypes.IQM ? new double[] {0.25, 0.75} : null); + // Resolve requested probabilities. + final double[] quantiles; + if(input2 != null) { + if(input2.isMatrix()) + quantiles = ec.getMatrixInput(input2).getDenseBlockValues(); + else if(input2.isScalar()) + quantiles = new double[] {ec.getScalarInput(input2).getDoubleValue()}; + else + throw new DMLRuntimeException( + "QuantilePickFEDInstruction: unsupported input2 data type " + input2.getDataType()); + } + else if(_type == OperationTypes.MEDIAN) + quantiles = new double[] {0.5}; + else if(_type == OperationTypes.IQM) + quantiles = new double[] {0.25, 0.75}; + else + throw new DMLRuntimeException("QuantilePickFEDInstruction: " + _type + " requires a probability input"); - if (input2 != null && input2.isMatrix()) + if(input2 != null && input2.isMatrix()) ec.releaseMatrixInput(input2.getName()); - // Find min and max + // Fetch per-worker min/max/weights. long varID = FederationUtils.getNextFedDataID(); List minMax = new ArrayList<>(); fedMap.mapParallel(varID, (range, data) -> { @@ -253,9 +266,7 @@ public void processRowQPick(ExecutionContext ec) { new QuantilePickFEDInstruction.MinMax(data.getVarID()))).get(); if(!response.isSuccessful()) response.throwExceptionFromResponse(); - double[] rangeMinMax = (double[]) response.getData()[0]; - minMax.add(rangeMinMax); - + minMax.add((double[]) response.getData()[0]); return null; } catch(Exception e) { @@ -263,174 +274,253 @@ public void processRowQPick(ExecutionContext ec) { } }); - // Find weights sum, min and max - double globalMin = Double.MAX_VALUE, globalMax = Double.MIN_VALUE, vectorLength = in.getNumColumns() == 2 ? 0 : in.getNumRows(), sumWeights = 0.0; + double globalMin = Double.MAX_VALUE, globalMax = Double.MIN_VALUE; + int vectorLength = in.getNumColumns() == 2 ? 0 : (int) in.getNumRows(); + double sumWeights = 0.0; for(double[] values : minMax) { globalMin = Math.min(globalMin, values[0]); globalMax = Math.max(globalMax, values[1]); if(in.getNumColumns() == 2) - vectorLength += values[2]; + vectorLength += (int) values[2]; sumWeights += values[3]; } - // Average for median - average = average && (in.getNumColumns() == 2 ? sumWeights : in.getNumRows()) % 2 == 0; + if(_type == OperationTypes.IQM) { + computeIqm(ec, in, fedMap, varID, vectorLength, globalMin, globalMax); + return; + } - // If multiple quantiles take first histogram and reuse bins, otherwise recursively get bin with result - int numBuckets = 256; // (int) Math.round(in.getNumRows() / 2.0); - int quantileIndex = quantiles != null && quantiles.length == 1 ? (int) Math.round(vectorLength * quantiles[0]) : -1; + // VALUEPICK / MEDIAN: R quantile type 7 — for each probability derive rank pair (lo, hi, g), + // look up the deduplicated ranks once through the shared multi-rank pipeline, then interpolate. + final long N = in.getNumColumns() == 2 ? Math.round(sumWeights) : vectorLength; + final int[] los = new int[quantiles.length]; + final int[] his = new int[quantiles.length]; + final double[] gs = new double[quantiles.length]; + final Set rankSet = new TreeSet<>(); + for(int i = 0; i < quantiles.length; i++) { + final double[] r = MatrixBlock.computeType7Rank(N, quantiles[i]); + los[i] = (int) r[0]; + his[i] = (int) r[1]; + gs[i] = r[2]; + rankSet.add(los[i]); + if(gs[i] > 0.0 && his[i] != los[i]) + rankSet.add(his[i]); + } + final int[] ranks = rankSet.stream().mapToInt(Integer::intValue).toArray(); - T ret = createHistogram(in, (int) vectorLength, globalMin, globalMax, numBuckets, quantileIndex, average); + final Map rankToValue = pickMultipleRanks(in, ranks, vectorLength, varID, globalMin, + globalMax); - // Compute and set results - if(quantiles != null && quantiles.length > 1) { - double finalVectorLength = vectorLength; - quantiles = Arrays.stream(quantiles).map(val -> (int) Math.round(finalVectorLength * val)).toArray(); - computeMultipleQuantiles(ec, in, (int[]) ret, quantiles, (int) vectorLength, varID, (globalMax-globalMin) / numBuckets, globalMin, _type, false); + if(quantiles.length == 1) { + ec.setScalarOutput(output.getName(), + new DoubleObject(interpolateType7(los[0], his[0], gs[0], rankToValue))); + } + else { + MatrixBlock out = new MatrixBlock(quantiles.length, 1, false); + for(int i = 0; i < quantiles.length; i++) + out.set(i, 0, interpolateType7(los[i], his[i], gs[i], rankToValue)); + ec.setMatrixOutput(output.getName(), out); } - else - getSingleQuantileResult(ret, ec, fedMap, varID, average, false, (int) vectorLength, null); } - private MatrixBlock computeMultipleQuantiles(ExecutionContext ec, MatrixObject in, int[] bucketsFrequencies, double[] quantiles, - int vectorLength, long varID, double bucketRange, double min, OperationTypes type, boolean returnOutput) { - MatrixBlock out = new MatrixBlock(quantiles.length, 1, false); - ImmutableTriple>[] bucketsWithIndex = new ImmutableTriple[quantiles.length]; + private static double interpolateType7(int lo, int hi, double g, Map rankToValue) { + final double loVal = rankToValue.get(lo); + return (g == 0.0 || hi == lo) ? loVal : (1.0 - g) * loVal + g * rankToValue.get(hi); + } - // Find bins with each quantile for first histogram - int sizeBeforeTmp = 0, sizeBefore = 0, countFoundBins = 0; - for(int j = 0; j < bucketsFrequencies.length; j++) { - sizeBeforeTmp += bucketsFrequencies[j]; + // IQM is a trimmed weighted mean, not an R type-7 pick. Uses raw ceil-based q25/q75 boundaries and the closed-form + // boundary correction — kept structurally identical to the pre-3953 IQM math so IQMTest stays a green guardrail. + private void computeIqm(ExecutionContext ec, MatrixObject in, FederationMap fedMap, long varID, int vectorLength, + double globalMin, double globalMax) { + final int q25Rank = (int) Math.ceil(0.25 * vectorLength); + final int q75Rank = (int) Math.ceil(0.75 * vectorLength); - for(int i = 0; i < quantiles.length; i++) { + final double bucketRange = (globalMax - globalMin) / NUM_BUCKETS; + final int[] bucketsFrequencies = createHistogram(in, vectorLength, globalMin, globalMax, NUM_BUCKETS, -1); - ImmutablePair bucketWithQ; + final int[] ranks = new int[] {q25Rank, q75Rank}; + final ImmutableTriple>[] bucketsWithIndex = locateInitialBuckets( + bucketsFrequencies, ranks, globalMin, bucketRange); - if(quantiles[i] > sizeBefore && quantiles[i] <= sizeBeforeTmp) { - bucketWithQ = new ImmutablePair<>(min + (j * bucketRange), min + ((j+1) * bucketRange)); - bucketsWithIndex[i] = new ImmutableTriple<>(quantiles[i] == 1 ? 1 : - (int) quantiles[i] - sizeBefore, bucketsFrequencies[j], bucketWithQ); - countFoundBins++; - } + double q25Left = 0, q25Right = 0, q75Left = 0, q75Right = 0; + for(int i = 0; i < ranks.length; i++) { + final Object hist = refineBucket(in, vectorLength, bucketsWithIndex[i]); + final double left = hist instanceof ImmutablePair ? ((ImmutablePair) hist).left : (Double) hist; + final double right = hist instanceof ImmutablePair ? ((ImmutablePair) hist).right : (Double) hist; + if(i == 0) { + q25Left = left; + q25Right = right; + } + else { + q75Left = left; + q75Right = right; } - - sizeBefore = sizeBeforeTmp; - if(countFoundBins == quantiles.length) - break; } - // Find each quantile bin recursively - Map retBuckets = new HashMap<>(); + if(q25Right == q75Right) { + ec.setScalarOutput(output.getName(), new DoubleObject(q25Left)); + return; + } - double q25Left = 0, q25Right = 0, q75Left = 0, q75Right = 0; - for(int i = 0; i < bucketsWithIndex.length; i++) { - int nextNumBuckets = bucketsWithIndex[i].middle < 100 ? bucketsWithIndex[i].middle * 2 : (int) Math.round(bucketsWithIndex[i].middle / 2.0); - T hist = createHistogram(in, vectorLength, bucketsWithIndex[i].right.left, bucketsWithIndex[i].right.right, nextNumBuckets, bucketsWithIndex[i].left, false); - - if(_type == OperationTypes.IQM) { - q25Right = i == 0 ? hist instanceof ImmutablePair ? ((ImmutablePair)hist).right : (Double) hist : q25Right; - q25Left = i == 0 ? hist instanceof ImmutablePair ? ((ImmutablePair)hist).left : (Double) hist : q25Left; - q75Right = i == 1 ? hist instanceof ImmutablePair ? ((ImmutablePair)hist).right : (Double) hist : q75Right; - q75Left = i == 1 ? hist instanceof ImmutablePair ? ((ImmutablePair)hist).left : (Double) hist : q75Left; - } else { - if(hist instanceof ImmutablePair) - retBuckets.put(i, hist); // set value if returned double instead of bin - else - out.set(i, 0, (Double) hist); + final ImmutablePair iqmRange = new ImmutablePair<>(q25Right, q75Right); + final ImmutablePair bounds = new ImmutablePair<>(q25Left, q75Left); + final List perWorker = new ArrayList<>(); + fedMap.mapParallel(varID, (range, data) -> { + try { + FederatedResponse response = data + .executeFederatedOperation(new FederatedRequest(FederatedRequest.RequestType.EXEC_UDF, -1, + new QuantilePickFEDInstruction.GetValuesInRange(data.getVarID(), iqmRange, true, bounds))) + .get(); + if(!response.isSuccessful()) + response.throwExceptionFromResponse(); + perWorker.add((double[]) response.getData()[0]); + return null; } - } + catch(Exception e) { + throw new DMLRuntimeException(e); + } + }); + + double sum = 0, q25Part = 0, q25Val = 0, q75Part = 0, q75Val = 0; + for(double[] vals : perWorker) { + sum += vals[0]; + q25Part += vals[1]; + q25Val += vals[2]; + q75Part += vals[3]; + q75Val += vals[4]; + } + q25Part -= (0.25 * vectorLength); + q75Part -= (0.75 * vectorLength); + final double result = (sum + q25Part * q25Val - q75Part * q75Val) / (vectorLength * 0.5); + ec.setScalarOutput(output.getName(), new DoubleObject(result)); + } - if(type == OperationTypes.IQM) { - ImmutablePair IQMRange = new ImmutablePair<>(q25Right, q75Right); - if(q25Right == q75Right) - ec.setScalarOutput(output.getName(), new DoubleObject(q25Left)); + // Look up the value at each requested (deduplicated, sorted) 1-based rank. Builds the coarse histogram once and + // refines per rank — the multi-rank pipeline the type-7 pair lookup rides on. + private Map pickMultipleRanks(MatrixObject in, int[] ranks, int vectorLength, long varID, + double globalMin, double globalMax) { + final Map result = new HashMap<>(); + if(ranks.length == 0) + return result; + + // Single rank: skip the shared-histogram scaffolding and let createHistogram do its own initial build + refine. + if(ranks.length == 1) { + final Object hist = createHistogram(in, vectorLength, globalMin, globalMax, NUM_BUCKETS, ranks[0]); + if(hist instanceof ImmutablePair) + result.put(ranks[0], + fetchValueInRange(in.getFedMapping(), varID, (ImmutablePair) hist)); else - getSingleQuantileResult(IQMRange, ec, in.getFedMapping(), varID, false, true, vectorLength, new ImmutablePair<>(q25Left, q75Left)); + result.put(ranks[0], (Double) hist); + return result; } - else { - if(!retBuckets.isEmpty()) { - // Search for values within bucket range where it as returned - in.getFedMapping().mapParallel(varID, (range, data) -> { - try { - FederatedResponse response = data.executeFederatedOperation(new FederatedRequest( - FederatedRequest.RequestType.EXEC_UDF, - -1, - new QuantilePickFEDInstruction.GetValuesInRanges(data.getVarID(), quantiles.length, (HashMap>) retBuckets))).get(); - if(!response.isSuccessful()) - response.throwExceptionFromResponse(); - - // Add results by row - MatrixBlock tmp = (MatrixBlock) response.getData()[0]; - synchronized(out) { - out.binaryOperationsInPlace(InstructionUtils.parseBinaryOperator(Opcodes.PLUS.toString()), tmp); - } - return null; - } - catch(Exception e) { - throw new DMLRuntimeException(e); - } - }); - } - if(returnOutput) - return out; + + final double bucketRange = (globalMax - globalMin) / NUM_BUCKETS; + final int[] bucketsFrequencies = createHistogram(in, vectorLength, globalMin, globalMax, NUM_BUCKETS, -1); + final ImmutableTriple>[] bucketsWithIndex = locateInitialBuckets( + bucketsFrequencies, ranks, globalMin, bucketRange); + + final HashMap> retBuckets = new HashMap<>(); + for(int i = 0; i < ranks.length; i++) { + final Object hist = refineBucket(in, vectorLength, bucketsWithIndex[i]); + if(hist instanceof ImmutablePair) + retBuckets.put(i, (ImmutablePair) hist); else - ec.setMatrixOutput(output.getName(), out); + result.put(ranks[i], (Double) hist); } - return null; - } - private void getSingleQuantileResult(T ret, ExecutionContext ec, FederationMap fedMap, long varID, boolean average, boolean isIQM, int vectorLength, ImmutablePair iqmRange) { - double result = 0.0, q25Part = 0, q25Val = 0, q75Val = 0, q75Part = 0; - if(ret instanceof ImmutablePair) { - // Search for values within bucket range - List values = new ArrayList<>(); - List iqmValues = new ArrayList<>(); - fedMap.mapParallel(varID, (range, data) -> { + if(!retBuckets.isEmpty()) { + final MatrixBlock resolved = new MatrixBlock(ranks.length, 1, false); + in.getFedMapping().mapParallel(varID, (range, data) -> { try { - FederatedResponse response = data.executeFederatedOperation(new FederatedRequest(FederatedRequest.RequestType.EXEC_UDF, -1, - new QuantilePickFEDInstruction.GetValuesInRange(data.getVarID(), (ImmutablePair) ret, isIQM, iqmRange))).get(); + FederatedResponse response = data.executeFederatedOperation(new FederatedRequest( + FederatedRequest.RequestType.EXEC_UDF, -1, + new QuantilePickFEDInstruction.GetValuesInRanges(data.getVarID(), ranks.length, retBuckets))) + .get(); if(!response.isSuccessful()) response.throwExceptionFromResponse(); - if(isIQM) - iqmValues.add((double[]) response.getData()[0]); - else - values.add((double) response.getData()[0]); + MatrixBlock tmp = (MatrixBlock) response.getData()[0]; + synchronized(resolved) { + resolved.binaryOperationsInPlace(InstructionUtils.parseBinaryOperator(Opcodes.PLUS.toString()), + tmp); + } return null; } catch(Exception e) { throw new DMLRuntimeException(e); } }); + for(Map.Entry> entry : retBuckets.entrySet()) + result.put(ranks[entry.getKey()], resolved.get(entry.getKey(), 0)); + } + return result; + } - if(isIQM) { - for(double[] vals : iqmValues) { - result += vals[0]; - q25Part += vals[1]; - q25Val += vals[2]; - q75Part += vals[3]; - q75Val += vals[4]; - } - q25Part -= (0.25 * vectorLength); - q75Part -= (0.75 * vectorLength); - } else - result = values.stream().reduce(0.0, Double::sum); - - } else - result = (Double) ret; + // Refine a coarse-histogram bucket into a finer sub-histogram covering just that bucket's range, and recurse + // into it for the given target rank. Shared by computeIqm and pickMultipleRanks so the nextNumBuckets heuristic + // lives in one place. Returns either the final value (Double) or the bucket range (ImmutablePair) — same + // polymorphic shape as createHistogram, callers instanceof-check. + private Object refineBucket(MatrixObject in, int vectorLength, + ImmutableTriple> bucketWithIndex) { + final int nextNumBuckets = bucketWithIndex.middle < 100 ? bucketWithIndex.middle * + 2 : (int) Math.round(bucketWithIndex.middle / 2.0); + return createHistogram(in, vectorLength, bucketWithIndex.right.left, bucketWithIndex.right.right, + nextNumBuckets, bucketWithIndex.left); + } - result = average ? result / 2 : (isIQM ? ((result + q25Part*q25Val - q75Part*q75Val) / (vectorLength * 0.5)) : result); + // Scan the coarse histogram once and record, for each target rank, the bucket range containing it plus the + // rank's offset within that bucket. Extracted from the multi-rank pipeline so IQM can reuse it verbatim. + // Triple layout per rank i: left = rank offset within the bucket (1-based, i.e. how many entries into the + // bucket the target sits), middle = bucket frequency, right = (bucketMin, bucketMax) sub-range to recurse into. + private static ImmutableTriple>[] locateInitialBuckets( + int[] bucketsFrequencies, int[] ranks, double globalMin, double bucketRange) { + final ImmutableTriple>[] bucketsWithIndex = new ImmutableTriple[ranks.length]; + int sizeBeforeTmp = 0, sizeBefore = 0, countFoundBins = 0; + for(int j = 0; j < bucketsFrequencies.length; j++) { + sizeBeforeTmp += bucketsFrequencies[j]; + for(int i = 0; i < ranks.length; i++) { + if(bucketsWithIndex[i] == null && ranks[i] > sizeBefore && ranks[i] <= sizeBeforeTmp) { + ImmutablePair bucketWithR = new ImmutablePair<>(globalMin + (j * bucketRange), + globalMin + ((j + 1) * bucketRange)); + bucketsWithIndex[i] = new ImmutableTriple<>(ranks[i] == 1 ? 1 : ranks[i] - sizeBefore, + bucketsFrequencies[j], bucketWithR); + countFoundBins++; + } + } + sizeBefore = sizeBeforeTmp; + if(countFoundBins == ranks.length) + break; + } + return bucketsWithIndex; + } - ec.setScalarOutput(output.getName(), new DoubleObject(result)); + private double fetchValueInRange(FederationMap fedMap, long varID, ImmutablePair range) { + final List values = new ArrayList<>(); + fedMap.mapParallel(varID, (r, data) -> { + try { + FederatedResponse response = data + .executeFederatedOperation(new FederatedRequest(FederatedRequest.RequestType.EXEC_UDF, -1, + new QuantilePickFEDInstruction.GetValuesInRange(data.getVarID(), range, false, null))) + .get(); + if(!response.isSuccessful()) + response.throwExceptionFromResponse(); + values.add((double) response.getData()[0]); + return null; + } + catch(Exception e) { + throw new DMLRuntimeException(e); + } + }); + return values.stream().reduce(0.0, Double::sum); } - public T createHistogram(CacheableData in, int vectorLength, double globalMin, double globalMax, int numBuckets, int quantileIndex, boolean average) { + public T createHistogram(CacheableData in, int vectorLength, double globalMin, double globalMax, + int numBuckets, int quantileIndex) { FederationMap fedMap = in.getFedMapping(); List hists = new ArrayList<>(); List> distincts = new ArrayList<>(); - double bucketRange = (globalMax-globalMin) / numBuckets; - boolean isEvenNumRows = vectorLength % 2 == 0; + double bucketRange = (globalMax - globalMin) / numBuckets; // Create histograms long varID = FederationUtils.getNextFedDataID(); @@ -462,49 +552,48 @@ public T createHistogram(CacheableData in, int vectorLength, double glob return (T) bucketsFrequencies; // Find bucket with quantile - ImmutableTriple> bucketWithIndex = getBucketWithIndex(bucketsFrequencies, globalMin, quantileIndex, average, isEvenNumRows, bucketRange); + ImmutableTriple> bucketWithIndex = getBucketWithIndex( + bucketsFrequencies, globalMin, quantileIndex, bucketRange); - // Check if can terminate + // Check if we can terminate early using merged per-worker distincts. Set distinctValues = distincts.stream().flatMap(Set::stream).collect(Collectors.toSet()); + if(distinctValues.size() > quantileIndex - 1) + return (T) distinctValues.stream().sorted().toArray()[quantileIndex > 0 ? quantileIndex - 1 : 0]; - if(distinctValues.size() > quantileIndex-1 && !average) - return (T) distinctValues.stream().sorted().toArray()[quantileIndex > 0 ? quantileIndex-1 : 0]; - - if(average && distinctValues.size() > quantileIndex) { - Double[] distinctsSorted = distinctValues.stream().flatMap(Stream::of).sorted().toArray(Double[]::new); - Double medianSum = Double.sum(distinctsSorted[quantileIndex-1], distinctsSorted[quantileIndex]); - return (T) medianSum; - } - - if((average && distinctValues.size() == 2) || (!average && distinctValues.size() == 1)) + if(distinctValues.size() == 1) return (T) distinctValues.stream().reduce(0.0, Double::sum); ImmutablePair finalBucketWithQ = bucketWithIndex.right; - List distinctInNewBucket = distinctValues.stream().filter( e -> e >= finalBucketWithQ.left && e <= finalBucketWithQ.right).collect(Collectors.toList()); - if((distinctInNewBucket.size() == 1 && !average) || (average && distinctInNewBucket.size() == 2)) - return (T) distinctInNewBucket.stream().reduce(0.0, Double::sum); + List distinctInNewBucket = distinctValues.stream() + .filter(e -> e >= finalBucketWithQ.left && e <= finalBucketWithQ.right).collect(Collectors.toList()); + if(distinctInNewBucket.size() == 1) + return (T) distinctInNewBucket.get(0); - if(!average) { - Set distinctsSet = new HashSet<>(distinctInNewBucket); - if(distinctsSet.size() == 1) - return (T) distinctsSet.toArray()[0]; - } + Set distinctsSet = new HashSet<>(distinctInNewBucket); + if(distinctsSet.size() == 1) + return (T) distinctsSet.toArray()[0]; - if(distinctValues.size() == 1 || (bucketWithIndex.middle == 1 && !average) || (bucketWithIndex.middle == 2 && isEvenNumRows && average) || - globalMin == globalMax) + if(bucketWithIndex.middle == 1 || globalMin == globalMax) return (T) bucketWithIndex.right; - int nextNumBuckets = bucketWithIndex.middle < 100 ? bucketWithIndex.middle * 2 : (int) Math.round(bucketWithIndex.middle / 2.0); + int nextNumBuckets = bucketWithIndex.middle < 100 ? bucketWithIndex.middle * + 2 : (int) Math.round(bucketWithIndex.middle / 2.0); // Add more bins to not stuck - if(numBuckets == nextNumBuckets && globalMin == bucketWithIndex.right.left && globalMax == bucketWithIndex.right.right) { + if(numBuckets == nextNumBuckets && globalMin == bucketWithIndex.right.left && + globalMax == bucketWithIndex.right.right) { nextNumBuckets *= 2; } - return createHistogram(in, vectorLength, bucketWithIndex.right.left, bucketWithIndex.right.right, nextNumBuckets, bucketWithIndex.left, average); + return createHistogram(in, vectorLength, bucketWithIndex.right.left, bucketWithIndex.right.right, + nextNumBuckets, bucketWithIndex.left); } - private ImmutableTriple> getBucketWithIndex(int[] bucketFrequencies, double min, int quantileIndex, boolean average, boolean isEvenNumRows, double bucketRange) { + // Locate the single bucket containing quantileIndex in a refinement histogram (called during recursion). + // Triple layout: left = rank offset within the bucket (1-based), middle = bucket frequency, + // right = (bucketMin, bucketMax) sub-range to recurse into next. + private ImmutableTriple> getBucketWithIndex(int[] bucketFrequencies, + double min, int quantileIndex, double bucketRange) { int sizeBeforeTmp = 0, sizeBefore = 0, bucketWithQSize = 0; ImmutablePair bucketWithQ = null; @@ -516,17 +605,7 @@ private ImmutableTriple> getBuck bucketWithQSize = bucketFrequencies[i]; sizeBeforeTmp -= bucketWithQSize; sizeBefore = sizeBeforeTmp; - - if(!average || sizeBefore + bucketWithQSize >= quantileIndex + 1) - break; - } else if(quantileIndex + 1 <= sizeBeforeTmp + bucketWithQSize && isEvenNumRows && average) { - // Add right bin that contains second index - int bucket2Size = bucketFrequencies[i]; - if (bucket2Size != 0) { - bucketWithQ = new ImmutablePair<>(bucketWithQ.left, tmpBinLeft + bucketRange); - bucketWithQSize += bucket2Size; - break; - } + break; } tmpBinLeft += bucketRange; } @@ -705,9 +784,11 @@ public void processColumnQPick(ExecutionContext ec) { new QuantilePickFEDInstruction.ColIQM(data.getVarID()))).get(); break; case MEDIAN: + // MEDIAN is VALUEPICK at p = 0.5 once the kernel handles R type 7 at both parities. response = data .executeFederatedOperation(new FederatedRequest(FederatedRequest.RequestType.EXEC_UDF, -1, - new QuantilePickFEDInstruction.ColMedian(data.getVarID()))).get(); + new QuantilePickFEDInstruction.ValuePick(data.getVarID(), new DoubleObject(0.5)))) + .get(); break; default: throw new DMLRuntimeException("Unsupported qpick operation type: "+_type); @@ -756,10 +837,10 @@ public FederatedResponse execute(ExecutionContext ec, Data... data) { MatrixBlock picked; if (_quantiles.getLength() == 1) { return new FederatedResponse(FederatedResponse.ResponseType.SUCCESS, - new Object[] {mb.pickValue(_quantiles.get(0, 0), mb.getLength() % 2 == 0)}); + new Object[] {mb.pickValue(_quantiles.get(0, 0))}); } else { - picked = mb.pickValues(_quantiles, new MatrixBlock(), mb.getLength() % 2 == 0); + picked = mb.pickValues(_quantiles, new MatrixBlock()); return new FederatedResponse(FederatedResponse.ResponseType.SUCCESS, new Object[] {picked}); } @@ -851,23 +932,4 @@ public Pair getLineageItem(ExecutionContext ec) { return null; } } - - private static class ColMedian extends FederatedUDF { - - private static final long serialVersionUID = -2808597461054603816L; - - protected ColMedian(long input) { - super(new long[] {input}); - } - @Override - public FederatedResponse execute(ExecutionContext ec, Data... data) { - MatrixBlock mb = ((MatrixObject)data[0]).acquireReadAndRelease(); - return new FederatedResponse(FederatedResponse.ResponseType.SUCCESS, - new Object[] {mb.median()}); - } - @Override - public Pair getLineageItem(ExecutionContext ec) { - return null; - } - } } diff --git a/src/main/java/org/apache/sysds/runtime/instructions/spark/QuantilePickSPInstruction.java b/src/main/java/org/apache/sysds/runtime/instructions/spark/QuantilePickSPInstruction.java index 75f84882478..9d6c8704e0a 100644 --- a/src/main/java/org/apache/sysds/runtime/instructions/spark/QuantilePickSPInstruction.java +++ b/src/main/java/org/apache/sysds/runtime/instructions/spark/QuantilePickSPInstruction.java @@ -32,7 +32,6 @@ import org.apache.sysds.runtime.instructions.InstructionUtils; import org.apache.sysds.runtime.instructions.cp.CPOperand; import org.apache.sysds.runtime.instructions.cp.DoubleObject; -import org.apache.sysds.runtime.instructions.cp.ScalarObject; import org.apache.sysds.runtime.instructions.spark.utils.RDDAggregateUtils; import org.apache.sysds.runtime.matrix.data.MatrixBlock; import org.apache.sysds.runtime.matrix.data.MatrixIndexes; @@ -47,7 +46,6 @@ import java.util.Collections; import java.util.Iterator; import java.util.List; -import java.util.stream.IntStream; public class QuantilePickSPInstruction extends BinarySPInstruction { private OperationTypes _type = null; @@ -110,124 +108,159 @@ public void processInstruction(ExecutionContext ec) { //NOTE: no difference between inmem/mr pick (see related cp instruction), but wrt w/ w/o weights //(in contrast to cp instructions, w/o weights does not materializes weights of 1) - switch( _type ) { + switch(_type) { case VALUEPICK: { - if( input2.isScalar() ) { - ScalarObject quantile = ec.getScalarInput(input2); - double[] wt = getWeightedQuantileSummary(in, mc, - new double[] {quantile.getDoubleValue()}, true); - ec.setScalarOutput(output.getName(), new DoubleObject(wt[3])); + if(input2.isScalar()) { + double picked = pickQuantileValues(in, mc, + new double[] {ec.getScalarInput(input2).getDoubleValue()})[0]; + ec.setScalarOutput(output.getName(), new DoubleObject(picked)); } else { - double[] wt = getWeightedQuantileSummary(in, mc, - DataConverter.convertToDoubleVector(ec.getMatrixInput(input2.getName())), true); + double[] values = pickQuantileValues(in, mc, + DataConverter.convertToDoubleVector(ec.getMatrixInput(input2.getName()))); ec.releaseMatrixInput(input2.getName()); - int qlen = wt.length/3; - MatrixBlock out = new MatrixBlock(qlen,1,false); - IntStream.range(0, out.getNumRows()) - .forEach(i -> out.set(i, 0, wt[2*qlen+i+1])); + MatrixBlock out = new MatrixBlock(values.length, 1, false); + for(int i = 0; i < values.length; i++) + out.set(i, 0, values[i]); ec.setMatrixOutput(output.getName(), out); } break; } case MEDIAN: { - double[] wt = getWeightedQuantileSummary(in, mc, new double[] {0.5}, true); - ec.setScalarOutput(output.getName(), new DoubleObject(wt[3])); + double median = pickQuantileValues(in, mc, new double[] {0.5})[0]; + ec.setScalarOutput(output.getName(), new DoubleObject(median)); break; } case IQM: { - double[] wt = getWeightedQuantileSummary(in, mc, new double[] {0.25, 0.75}, false); - long key25 = (long) Math.ceil(wt[1]); - long key75 = (long) Math.ceil(wt[2]); - JavaPairRDD out = in - .filter(new FilterFunction(key25 + 1, key75, mc.getBlocksize())) - .mapToPair(new ExtractAndSumFunction(key25 + 1, key75, mc.getBlocksize())); - double sum = RDDAggregateUtils.sumStable(out).get(0, 0); - double val = MatrixBlock.computeIQMCorrection( - sum, wt[0], wt[3], wt[5], wt[4], wt[6]); + double val = computeIqm(in, mc); ec.setScalarOutput(output.getName(), new DoubleObject(val)); break; } - + default: - throw new DMLRuntimeException("Unsupported qpick operation type: "+_type); + throw new DMLRuntimeException("Unsupported qpick operation type: " + _type); } } - + /** - * Get a summary of weighted quantiles in in the following form: - * sum of weights, (keys of quantiles), (portions of quantiles), (values of quantiles) - * - * @param w rdd containing values and optionally weights, sorted by value - * @param mc matrix characteristics - * @param quantiles one or more quantiles between 0 and 1. - * @return a summary of weighted quantiles + * Pick one R quantile type 7 value per requested probability. Used by VALUEPICK / MEDIAN. Two-column input is a + * weighted sequence — treated as an expanded sorted sequence of length sum(weights) and picked with the same h/lo/ + * hi/g formula against cumulative weights. */ - private static double[] getWeightedQuantileSummary(JavaPairRDD w, - DataCharacteristics mc, double[] quantiles, boolean average) { - double[] ret = new double[3 * quantiles.length + 1]; - if(mc.getCols() == 2) // weighted - { - //sort blocks (values sorted but blocks and partitions are not) - w = w.sortByKey(); - - //compute cumsum weights per partition - //with assumption that partition aggregates fit into memory - List> partWeights = w + private static double[] pickQuantileValues(JavaPairRDD w, DataCharacteristics mc, + double[] quantiles) { + final int blen = mc.getBlocksize(); + final double[] values = new double[quantiles.length]; + if(mc.getCols() == 2) { + final JavaPairRDD sorted = w.sortByKey(); + final List> partWeights = sorted .mapPartitionsWithIndex(new SumWeightsFunction(), false).collect(); - - //compute sum of weights - ret[0] = partWeights.stream().mapToDouble(p -> p._2()).sum(); - - //compute total cumsum and determine partitions - double[] qdKeys = new double[quantiles.length]; - long[] qiKeys = new long[quantiles.length]; - int[] partitionIDs = new int[quantiles.length]; - double[] offsets = new double[quantiles.length]; - for( int i=0; i psum : partWeights ) { - double tmp = cumSum + psum._2(); - for(int i=0; i= qiKeys[i] && partitionIDs[i] == 0 ) { - partitionIDs[i] = psum._1(); - offsets[i] = cumSum; - } - cumSum = tmp; + final long sumWt = Math.round(partWeights.stream().mapToDouble(p -> p._2()).sum()); + // Two keys per quantile (lo, hi) for type-7 interpolation; qdKey == qiKey since posPart is unused here. + final int nk = 2 * quantiles.length; + final double[] qdKeys = new double[nk]; + final long[] qiKeys = new long[nk]; + final double[] gs = new double[quantiles.length]; + for(int i = 0; i < quantiles.length; i++) { + final double[] r = MatrixBlock.computeType7Rank(sumWt, quantiles[i]); + qiKeys[2 * i] = (long) r[0]; + qiKeys[2 * i + 1] = (long) r[1]; + qdKeys[2 * i] = r[0]; + qdKeys[2 * i + 1] = r[1]; + gs[i] = r[2]; } - - //get keys and values for quantile cutoffs - List> qVals = w - .mapPartitionsWithIndex(new ExtractWeightedQuantileFunction( - mc, qdKeys, qiKeys, partitionIDs, offsets), false).collect(); - for( Tuple2 qVal : qVals ) { - ret[qVal._1()+1] = qVal._2()[0]; - ret[qVal._1()+quantiles.length+1] = qVal._2()[1]; - ret[qVal._1()+2*quantiles.length+1] = qVal._2()[2]; + final double[][] triples = extractWeightedTriples(sorted, mc, partWeights, qdKeys, qiKeys); + for(int i = 0; i < quantiles.length; i++) { + final double loVal = triples[2 * i][2]; + // hi == lo covers the p == 1 clamp. + values[i] = (gs[i] == 0.0 || qiKeys[2 * i + 1] == qiKeys[2 * i]) ? loVal : (1.0 - gs[i]) * loVal + + gs[i] * triples[2 * i + 1][2]; } } else { - ret[0] = mc.getRows(); + final long N = mc.getRows(); for(int i = 0; i < quantiles.length; i++) { - ret[i + 1] = quantiles[i] * mc.getRows(); - ret[i + quantiles.length + 1] = Math.ceil(ret[i + 1]) - ret[i + 1]; - long key = (long) Math.ceil(ret[i + 1]); - ret[i + 2 * quantiles.length + 1] = lookupKey(w, key, mc.getBlocksize()); + final double[] r = MatrixBlock.computeType7Rank(N, quantiles[i]); + final long lo = (long) r[0], hi = (long) r[1]; + final double g = r[2]; + final double loVal = lookupKey(w, lo, blen); + values[i] = (g == 0.0 || hi == lo) ? loVal : (1.0 - g) * loVal + g * lookupKey(w, hi, blen); + } + } + return values; + } + + /** + * Compute the interquartile mean: trimmed mean of values between the ceil-based q25 and q75 ranks, with fractional + * boundary corrections applied via {@link MatrixBlock#computeIQMCorrection}. IQM uses raw ceil-based boundaries + * (not R type 7) so the middle-range sum and boundary portions align with the closed-form correction formula. + */ + private static double computeIqm(JavaPairRDD in, DataCharacteristics mc) { + final int blen = mc.getBlocksize(); + final double sumWt, q25Position, q75Position, q25Portion, q75Portion, q25Value, q75Value; + if(mc.getCols() == 2) { + final JavaPairRDD sorted = in.sortByKey(); + final List> partWeights = sorted + .mapPartitionsWithIndex(new SumWeightsFunction(), false).collect(); + sumWt = partWeights.stream().mapToDouble(p -> p._2()).sum(); + final double[] qdKeys = {0.25 * sumWt, 0.75 * sumWt}; + final long[] qiKeys = {(long) Math.ceil(qdKeys[0]), (long) Math.ceil(qdKeys[1])}; + final double[][] triples = extractWeightedTriples(sorted, mc, partWeights, qdKeys, qiKeys); + // For weighted data q25/q75Position is the sorted matrix row index (from extract), not + // 0.25 * sumWt in the expanded sequence — the IQM filter runs on RDD row-block coordinates. + q25Position = triples[0][0]; + q75Position = triples[1][0]; + q25Portion = triples[0][1]; + q75Portion = triples[1][1]; + q25Value = triples[0][2]; + q75Value = triples[1][2]; + } + else { + final long N = mc.getRows(); + sumWt = N; + q25Position = 0.25 * N; + q75Position = 0.75 * N; + q25Portion = Math.ceil(q25Position) - q25Position; + q75Portion = Math.ceil(q75Position) - q75Position; + q25Value = lookupKey(in, (long) Math.ceil(q25Position), blen); + q75Value = lookupKey(in, (long) Math.ceil(q75Position), blen); + } + final long key25 = (long) Math.ceil(q25Position); + final long key75 = (long) Math.ceil(q75Position); + JavaPairRDD out = in.filter(new FilterFunction(key25 + 1, key75, blen)) + .mapToPair(new ExtractAndSumFunction(key25 + 1, key75, blen)); + double sum = RDDAggregateUtils.sumStable(out).get(0, 0); + return MatrixBlock.computeIQMCorrection(sum, sumWt, q25Portion, q25Value, q75Portion, q75Value); + } - // average w/ next value for even-length arrays (mirrors CP QuantilePickCPInstruction) - if(average && mc.getRows() % 2 == 0 && key < (mc.getRows() - 1)) { - ret[i + 2 * quantiles.length + 1] += lookupKey(w, key + 1, mc.getBlocksize()); - ret[i + 2 * quantiles.length + 1] /= 2; + /** + * Locate the partition holding each ceil-based key by scanning cumulative per-partition weights, then invoke + * {@link ExtractWeightedQuantileFunction} to fetch the (position, posPart, value) triples. Returns a nk-length + * array indexed by the caller's key index, one triple per key. + */ + private static double[][] extractWeightedTriples(JavaPairRDD sorted, + DataCharacteristics mc, List> partWeights, double[] qdKeys, long[] qiKeys) { + final int nk = qiKeys.length; + final int[] partitionIDs = new int[nk]; + final double[] offsets = new double[nk]; + double cumSum = 0; + for(Tuple2 psum : partWeights) { + final double tmp = cumSum + psum._2(); + for(int i = 0; i < nk; i++) + if(tmp >= qiKeys[i] && partitionIDs[i] == 0) { + partitionIDs[i] = psum._1(); + offsets[i] = cumSum; } - } + cumSum = tmp; } - - return ret; + final List> qVals = sorted.mapPartitionsWithIndex( + new ExtractWeightedQuantileFunction(mc, qdKeys, qiKeys, partitionIDs, offsets), false).collect(); + final double[][] triples = new double[nk][]; + for(Tuple2 qVal : qVals) + triples[qVal._1()] = qVal._2(); + return triples; } private static double lookupKey(JavaPairRDD in, long key, int blen) { diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/MatrixBlock.java b/src/main/java/org/apache/sysds/runtime/matrix/data/MatrixBlock.java index 7525dab2f7f..87c4e4f7dc6 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/MatrixBlock.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/MatrixBlock.java @@ -4749,142 +4749,117 @@ public double interQuartileMean() { //compute final IQM, incl. correction for q25 and q75 portions return computeIQMCorrection(sum, sum_wt, q25Part, q25Val, q75Part, q75Val); } - - public static double computeIQMCorrection(double sum, double sum_wt, - double q25Part, double q25Val, double q75Part, double q75Val) { - return (sum + q25Part*q25Val - q75Part*q75Val) / (sum_wt*0.5); + + public static double computeIQMCorrection(double sum, double sum_wt, double q25Part, double q25Val, double q75Part, + double q75Val) { + return (sum + q25Part * q25Val - q75Part * q75Val) / (sum_wt * 0.5); } - + + /** + * R quantile type 7 rank triple {lo, hi, g} for a sequence of length n and probability p. lo and hi are 1-based + * order-statistic indices (returned as double for a single primitive-array return; cast to long at the call site); + * g in [0, 1) is the interpolation weight. The picked quantile is (1 - g) * x[lo] + g * x[hi], reducing to x[lo] + * when g == 0 or hi == lo (the p == 1 upper-bound clamp). Used by single-block picking and the distributed and + * federated pick paths. + */ + public static double[] computeType7Rank(long n, double p) { + final double h = (n - 1) * p + 1.0; + final double lo = Math.max(1.0, Math.min(Math.floor(h), (double) n)); + return new double[] {lo, Math.min(lo + 1.0, (double) n), h - Math.floor(h)}; + } + /** * Pick the quantiles out of this matrix. If this matrix contains two columns it is weighted quantile picking. * If a single column it is unweighted. - * + * * Note the values are assumed to be sorted. - * + * * @param quantiles The quantiles to pick * @param ret The result matrix * @return The result matrix */ - public final MatrixBlock pickValues(MatrixValue quantiles, MatrixValue ret) { - return pickValues(quantiles, ret, false); - } - - public MatrixBlock pickValues(MatrixValue quantiles, MatrixValue ret, boolean average) { - MatrixBlock qs=checkType(quantiles); - - if ( qs.clen != 1 ) { + public MatrixBlock pickValues(MatrixValue quantiles, MatrixValue ret) { + MatrixBlock qs = checkType(quantiles); + + if(qs.clen != 1) { throw new DMLRuntimeException("Multiple quantiles can only be computed on a 1D matrix"); } - + MatrixBlock output = checkType(ret); - if(output==null) - output=new MatrixBlock(qs.rlen, qs.clen, false); // resulting matrix is mostly likely be dense + if(output == null) + output = new MatrixBlock(qs.rlen, qs.clen, false); // resulting matrix is mostly likely be dense else output.reset(qs.rlen, qs.clen, false); for(int i = 0; i < qs.rlen; i++) { - // FIXME: include the average parameter here to fix SYSTEMDS-3953 output.set(i, 0, this.pickValue(qs.get(i, 0))); } - + return output; } - + /** * Pick the median value from this matrix. If this matrix has two columns it is weighted picking using the * weight column, otherwise it is unweighted over the single column. - * + * * Note the values are assumed to be sorted. - * + * * @return The median value */ public double median() { - if(getNumColumns() == 1) - return pickValue(0.5, getNumRows() % 2 == 0); - double sum_wt = sumWeightForQuantile(); - return pickValue(0.5, sum_wt%2==0); + return pickValue(0.5); } /** - * Pick a specific quantile from this matrix. If this matrix has two columns it is weighted picking, otherwise it is unweighted. - * - * Note the values are assumed to be sorted. - * - * @param quantile The quantile to pick - * @return The quantile - */ - public final double pickValue(double quantile){ - return pickValue(quantile, false); - } - - /** - * Pick a specific quantile from this matrix. If this matrix has two columns it is weighted picking, otherwise it is unweighted. - * + * Pick a specific quantile from this matrix using R's default (type 7) definition: linear interpolation between the + * two adjacent order statistics. If this matrix has two columns the second is treated as integer weights. + * * Note the values are assumed to be sorted. - * - * @param quantile The quantile to pick - * @param average If the quantile is averaged. + * + * @param quantile The quantile in [0, 1] to pick * @return The quantile */ - public final double pickValue(double quantile, boolean average) { + public final double pickValue(double quantile) { if(this.getNumColumns() == 1) - return pickUnweightedValue(quantile, average); - return pickWeightedValue(quantile, average); + return pickUnweightedValue(quantile); + return pickWeightedValue(quantile); } - private double pickUnweightedValue(double quantile, boolean average) { - // Mirror the weighted convention (pickWeightedValue) with an implicit weight of 1 per value, so a single - // column yields the same quantile as the equivalent two-column (value, weight) representation: take the - // ceil-based rank and only average adjacent order statistics when an even number of values straddles it. - final int rows = getNumRows(); - average = average && (rows % 2 == 0); - final int pos = (int) Math.ceil(quantile * rows); // 1-based rank - final int i = Math.min(Math.max(pos - 1, 0), rows - 1); - if(average && pos > 0 && pos < rows) - return (get(i, 0) + get(i + 1, 0)) / 2; - return get(i, 0); + private double pickUnweightedValue(double quantile) { + final double[] r = computeType7Rank(getNumRows(), quantile); + final long lo = (long) r[0], hi = (long) r[1]; + final double g = r[2]; + final double loVal = get((int) (lo - 1), 0); + return (g == 0.0 || hi == lo) ? loVal : (1.0 - g) * loVal + g * get((int) (hi - 1), 0); } - private double pickWeightedValue(double quantile, boolean average) { - double sum_wt = sumWeightForQuantile(); - - // do averaging only if it is asked for; and sum_wt is even - average = average && (sum_wt%2 == 0); - - int pos = (int) Math.ceil(quantile*sum_wt); - - int t = 0, i=-1; + private double pickWeightedValue(double quantile) { + // R quantile type 7 generalized to integer weights: treat as expanded sorted sequence of length sum_wt. + final double[] r = computeType7Rank(Math.round(sumWeightForQuantile()), quantile); + final long lo = (long) r[0], hi = (long) r[1]; + final double g = r[2]; + final double loVal = valueAtWeightedRank(lo); + return (g == 0.0 || hi == lo) ? loVal : (1.0 - g) * loVal + g * valueAtWeightedRank(hi); + } + + private double valueAtWeightedRank(long rank) { + // Walk cumulative weights until we reach the requested 1-based rank in the expanded sequence. + final int rows = getNumRows(); + long t = 0; + int i = -1; do { i++; - t += get(i,1); - } while(t idx 0 - assertEquals("q=0.2", 10, mb.pickValue(0.2, false), 0); // rank ceil(1.0)=1 -> idx 0 - assertEquals("q=0.5", 30, mb.pickValue(0.5, false), 0); // rank ceil(2.5)=3 -> idx 2 - assertEquals("q=0.75", 40, mb.pickValue(0.75, false), 0); // rank ceil(3.75)=4 -> idx 3 - assertEquals("q=1.0", 50, mb.pickValue(1.0, false), 0); // rank ceil(5.0)=5 -> idx 4 - } - - @Test - public void pickOddLengthAverageSuppressed() { - // Odd number of values -> averaging is suppressed, so average matches no-average. + public void pickOddLength() { + // n=5: h = 4p + 1. At classical quantiles the rank is integer, so no interpolation is needed. MatrixBlock mb = singleColumn(new double[] {10, 20, 30, 40, 50}, false); - assertEquals("q=0.5 avg", 30, mb.pickValue(0.5, true), 0); - assertEquals("q=0.75 avg", 40, mb.pickValue(0.75, true), 0); - } - - @Test - public void pickEvenLengthAverage() { - // Even number of values -> averaging of adjacent order statistics applies. - MatrixBlock mb = singleColumn(new double[] {10, 20, 30, 40}, false); - assertEquals("q=0.25 avg", 15, mb.pickValue(0.25, true), 0); // rank 1 -> (idx0+idx1)/2 - assertEquals("q=0.375 avg", 25, mb.pickValue(0.375, true), 0); // rank ceil(1.5)=2 -> (idx1+idx2)/2 - assertEquals("q=0.5 avg", 25, mb.pickValue(0.5, true), 0); // rank 2 -> (idx1+idx2)/2 - assertEquals("q=0.75 avg", 35, mb.pickValue(0.75, true), 0); // rank 3 -> (idx2+idx3)/2 + assertEquals("q=0.0", 10, mb.pickValue(0.0), EPS); + assertEquals("q=0.2", 18, mb.pickValue(0.2), EPS); // h=1.8, 0.2*10 + 0.8*20 + assertEquals("q=0.25", 20, mb.pickValue(0.25), EPS); // h=2 -> x[2]=20 + assertEquals("q=0.5", 30, mb.pickValue(0.5), EPS); // h=3 -> x[3]=30 + assertEquals("q=0.75", 40, mb.pickValue(0.75), EPS); // h=4 -> x[4]=40 + assertEquals("q=1.0", 50, mb.pickValue(1.0), EPS); } @Test - public void pickEvenLengthNoAverage() { + public void pickEvenLength() { + // n=4: h = 3p + 1. Classical p != 0.5 land between order statistics and interpolate. MatrixBlock mb = singleColumn(new double[] {10, 20, 30, 40}, false); - assertEquals("q=0.25", 10, mb.pickValue(0.25, false), 0); // rank 1 -> idx 0 - assertEquals("q=0.5", 20, mb.pickValue(0.5, false), 0); // rank 2 -> idx 1 - assertEquals("q=0.75", 30, mb.pickValue(0.75, false), 0); // rank 3 -> idx 2 + assertEquals("q=0.25", 17.5, mb.pickValue(0.25), EPS); // h=1.75 -> 0.25*10 + 0.75*20 + assertEquals("q=0.375", 21.25, mb.pickValue(0.375), EPS); // h=2.125 -> 0.875*20 + 0.125*30 + assertEquals("q=0.5", 25, mb.pickValue(0.5), EPS); // h=2.5 -> 0.5*20 + 0.5*30 + assertEquals("q=0.75", 32.5, mb.pickValue(0.75), EPS); // h=3.25 -> 0.75*30 + 0.25*40 } @Test - public void pickAverageClampedAtTop() { - // Top quantile: rank reaches the last element so there is no successor to average with. + public void pickClampedAtTop() { + // Top quantile is clamped so no successor is required for interpolation. MatrixBlock even = singleColumn(new double[] {10, 20, 30, 40}, false); - assertEquals("even q=0.95 avg", 40, even.pickValue(0.95, true), 0); // rank ceil(3.8)=4 -> idx 3, no avg - assertEquals("even q=1.0 avg", 40, even.pickValue(1.0, true), 0); + assertEquals("even q=0.95", 38.5, even.pickValue(0.95), EPS); // h=3.85 -> 0.15*30 + 0.85*40 + assertEquals("even q=1.0", 40, even.pickValue(1.0), EPS); MatrixBlock odd = singleColumn(new double[] {10, 20, 30, 40, 50}, false); - assertEquals("odd q=0.95 avg", 50, odd.pickValue(0.95, true), 0); // odd -> avg suppressed + assertEquals("odd q=0.95", 48, odd.pickValue(0.95), EPS); // h=4.8 -> 0.2*40 + 0.8*50 } @Test public void pickSingleElement() { MatrixBlock mb = singleColumn(new double[] {42}, false); - assertEquals("q=0.0", 42, mb.pickValue(0.0, false), 0); - assertEquals("q=0.5", 42, mb.pickValue(0.5, false), 0); - assertEquals("q=1.0", 42, mb.pickValue(1.0, false), 0); - assertEquals("q=0.5 avg", 42, mb.pickValue(0.5, true), 0); - assertEquals("median", 42, mb.median(), 0); + assertEquals("q=0.0", 42, mb.pickValue(0.0), EPS); + assertEquals("q=0.5", 42, mb.pickValue(0.5), EPS); + assertEquals("q=1.0", 42, mb.pickValue(1.0), EPS); + assertEquals("median", 42, mb.median(), EPS); } @Test public void pickSparseSingleColumnWithZeros() { - // Sorted ascending including leading zeros, stored sparse. + // Sorted ascending including leading zeros, stored sparse. n=5, h = 4p + 1. MatrixBlock mb = singleColumn(new double[] {0, 0, 10, 20, 30}, true); - assertEquals("q=0.0", 0, mb.pickValue(0.0, false), 0); // rank 0 -> idx 0 (zero) - assertEquals("q=0.5", 10, mb.pickValue(0.5, false), 0); // rank ceil(2.5)=3 -> idx 2 - assertEquals("q=0.75", 20, mb.pickValue(0.75, false), 0); // rank ceil(3.75)=4 -> idx 3 - assertEquals("q=1.0", 30, mb.pickValue(1.0, false), 0); // rank 5 -> idx 4 + assertEquals("q=0.0", 0, mb.pickValue(0.0), EPS); + assertEquals("q=0.25", 0, mb.pickValue(0.25), EPS); // h=2 -> x[2]=0 + assertEquals("q=0.5", 10, mb.pickValue(0.5), EPS); // h=3 -> x[3]=10 + assertEquals("q=0.75", 20, mb.pickValue(0.75), EPS); // h=4 -> x[4]=20 + assertEquals("q=1.0", 30, mb.pickValue(1.0), EPS); // h=5 -> x[5]=30 } @Test public void medianSingleColumn() { - // Odd length -> middle element; even length -> average of the two middle elements. - assertEquals("odd median", 30, singleColumn(new double[] {10, 20, 30, 40, 50}, false).median(), 0); - assertEquals("even median", 25, singleColumn(new double[] {10, 20, 30, 40}, false).median(), 0); - assertEquals("sparse median", 10, singleColumn(new double[] {0, 0, 10, 20, 30}, true).median(), 0); + assertEquals("odd median", 30, singleColumn(new double[] {10, 20, 30, 40, 50}, false).median(), EPS); + assertEquals("even median", 25, singleColumn(new double[] {10, 20, 30, 40}, false).median(), EPS); + assertEquals("sparse median", 10, singleColumn(new double[] {0, 0, 10, 20, 30}, true).median(), EPS); } @Test @@ -123,7 +107,6 @@ public void pickSingleColumnMatchesDenseAndSparse() { MatrixBlock dense = singleColumn(v, false); MatrixBlock sparse = singleColumn(v, true); for(double q : new double[] {0.0, 0.1, 0.25, 0.5, 0.75, 0.9, 1.0}) - for(boolean avg : new boolean[] {false, true}) - assertEquals("q=" + q + " avg=" + avg, dense.pickValue(q, avg), sparse.pickValue(q, avg), 0); + assertEquals("q=" + q, dense.pickValue(q), sparse.pickValue(q), EPS); } } diff --git a/src/test/java/org/apache/sysds/test/functions/binary/matrix/QuantileTest.java b/src/test/java/org/apache/sysds/test/functions/binary/matrix/QuantileTest.java index 810e2614e7c..b06cf7be452 100644 --- a/src/test/java/org/apache/sysds/test/functions/binary/matrix/QuantileTest.java +++ b/src/test/java/org/apache/sysds/test/functions/binary/matrix/QuantileTest.java @@ -21,7 +21,6 @@ import java.util.HashMap; -import org.junit.Ignore; import org.junit.Test; import org.apache.sysds.common.Types.ExecMode; import org.apache.sysds.common.Types.ExecType; @@ -187,35 +186,68 @@ public void testMedianBugSP() { } @Test - @Ignore // FIXME: fix SYSTEMDS-3953 public void testQuartileArrayCP() { runQuantileTest(TEST_NAME6, 0, false, ExecType.CP); } @Test - @Ignore // FIXME: fix SYSTEMDS-3953 public void testQuartileArraySP() { runQuantileTest(TEST_NAME6, 0, false, ExecType.SPARK); } - private void runQuantileTest( String TEST_NAME, double p, boolean sparse, ExecType et) - { + // SYSTEMDS-3953: even-length quantile picks must match R type 7 at p != 0.5. The odd-length TEST_NAME1 + // cases above cannot exercise interpolation because at classical p the rank lands on an integer for odd n. + @Test + public void testQuantileEven1CP() { + runQuantileTest(TEST_NAME1, 0.25, false, ExecType.CP, 128); + } + + @Test + public void testQuantileEven2CP() { + runQuantileTest(TEST_NAME1, 0.50, false, ExecType.CP, 128); + } + + @Test + public void testQuantileEven3CP() { + runQuantileTest(TEST_NAME1, 0.75, false, ExecType.CP, 128); + } + + @Test + public void testQuantileEven1SP() { + runQuantileTest(TEST_NAME1, 0.25, false, ExecType.SPARK, 128); + } + + @Test + public void testQuantileEven2SP() { + runQuantileTest(TEST_NAME1, 0.50, false, ExecType.SPARK, 128); + } + + @Test + public void testQuantileEven3SP() { + runQuantileTest(TEST_NAME1, 0.75, false, ExecType.SPARK, 128); + } + + private void runQuantileTest(String TEST_NAME, double p, boolean sparse, ExecType et) { + runQuantileTest(TEST_NAME, p, sparse, et, rows); + } + + private void runQuantileTest(String TEST_NAME, double p, boolean sparse, ExecType et, int rowCount) { ExecMode platformOld = setExecMode(et); - + try { getAndLoadTestConfiguration(TEST_NAME); - + String HOME = SCRIPT_DIR + TEST_DIR; fullDMLScriptName = HOME + TEST_NAME + ".dml"; programArgs = new String[]{"-args", input("A"), Double.toString(p), output("R")}; fullRScriptName = HOME + TEST_NAME + ".R"; rCmd = "Rscript" + " " + fullRScriptName + " " + inputDir() + " " + p + " "+ expectedDir(); - + //generate actual dataset (always dense because values <=0 invalid) if( !TEST_NAME.equals(TEST_NAME4) ) { double sparsitya = sparse ? sparsity2 : sparsity1; - double[][] A = getRandomMatrix(rows, 1, 1, maxVal, sparsitya, 1236); + double[][] A = getRandomMatrix(rowCount, 1, 1, maxVal, sparsitya, 1236); writeInputMatrixWithMTD("A", A, true); }