From cebf29ef23282acb31219647c4fb9c619d8f5d91 Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 19 Jan 2018 16:44:45 +0100 Subject: [PATCH] trying something else --- .../IterativeMinMaxLinearEquationSolver.cpp | 204 ++++++++++++++++-- 1 file changed, 184 insertions(+), 20 deletions(-) diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 38313a57f..d45f1478e 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -679,7 +679,7 @@ namespace storm { template void performIterationStep(storm::storage::SparseMatrix const& A, std::vector const& b) { if (!decisionValueBlocks) { - performIterationStepUpdateDecisionValue2(A, b); + performIterationStepUpdateDecisionValue(A, b); } else { assert(decisionValue == getPrimaryBound()); auto xIt = x.rbegin(); @@ -719,6 +719,90 @@ namespace storm { } } + void performIterationStepMax(storm::storage::SparseMatrix const& A, std::vector const& b) { + if (!decisionValueBlocks) { + performIterationStepUpdateDecisionValueMax(A, b); + } else { + assert(decisionValue == upperBound); + auto xIt = x.rbegin(); + auto yIt = y.rbegin(); + auto groupStartIt = A.getRowGroupIndices().rbegin(); + uint64_t groupEnd = *groupStartIt; + ++groupStartIt; + for (auto groupStartIte = A.getRowGroupIndices().rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + // Perform the iteration for the first row in the group + uint64_t row = *groupStartIt; + ValueType xBest, yBest; + multiplyRow(row, A, b[row], xBest, yBest); + ++row; + // Only do more work if there are still rows in this row group + if (row != groupEnd) { + ValueType xi, yi; + ValueType bestValue = xBest + yBest * upperBound; + for (;row < groupEnd; ++row) { + // Get the multiplication results + multiplyRow(row, A, b[row], xi, yi); + ValueType currentValue = xi + yi * upperBound; + // Check if the current row is better then the previously found one + if (currentValue > bestValue) { + xBest = std::move(xi); + yBest = std::move(yi); + bestValue = std::move(currentValue); + } else if (currentValue == bestValue && yBest > yi) { + // If the value for this row is not strictly better, it might still be equal and have a better y value + xBest = std::move(xi); + yBest = std::move(yi); + } + } + } + *xIt = std::move(xBest); + *yIt = std::move(yBest); + } + } + } + + void performIterationStepMin(storm::storage::SparseMatrix const& A, std::vector const& b) { + if (!decisionValueBlocks) { + performIterationStepUpdateDecisionValueMin(A, b); + } else { + assert(decisionValue == lowerBound); + auto xIt = x.rbegin(); + auto yIt = y.rbegin(); + auto groupStartIt = A.getRowGroupIndices().rbegin(); + uint64_t groupEnd = *groupStartIt; + ++groupStartIt; + for (auto groupStartIte = A.getRowGroupIndices().rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + // Perform the iteration for the first row in the group + uint64_t row = *groupStartIt; + ValueType xBest, yBest; + multiplyRow(row, A, b[row], xBest, yBest); + ++row; + // Only do more work if there are still rows in this row group + if (row != groupEnd) { + ValueType xi, yi; + ValueType bestValue = xBest + yBest * lowerBound; + for (;row < groupEnd; ++row) { + // Get the multiplication results + multiplyRow(row, A, b[row], xi, yi); + ValueType currentValue = xi + yi * lowerBound; + // Check if the current row is better then the previously found one + if (currentValue < bestValue) { + xBest = std::move(xi); + yBest = std::move(yi); + bestValue = std::move(currentValue); + } else if (currentValue == bestValue && yBest > yi) { + // If the value for this row is not strictly better, it might still be equal and have a better y value + xBest = std::move(xi); + yBest = std::move(yi); + } + } + } + *xIt = std::move(xBest); + *yIt = std::move(yBest); + } + } + } + template void performIterationStepUpdateDecisionValue(storm::storage::SparseMatrix const& A, std::vector const& b) { auto xIt = x.rbegin(); @@ -800,30 +884,30 @@ namespace storm { } } - template - void performIterationStepUpdateDecisionValue2(storm::storage::SparseMatrix const& A, std::vector const& b) { - auto const& groupIndices = A.getRowGroupIndices(); - uint64_t group = groupIndices.size(); - while (group != 0) { - uint64_t const& groupEnd = groupIndices[group]; - --group; - uint64_t row = groupIndices[group]; + void performIterationStepUpdateDecisionValueMax(storm::storage::SparseMatrix const& A, std::vector const& b) { + auto xIt = x.rbegin(); + auto yIt = y.rbegin(); + auto groupStartIt = A.getRowGroupIndices().rbegin(); + uint64_t groupEnd = *groupStartIt; + ++groupStartIt; + for (auto groupStartIte = A.getRowGroupIndices().rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + // Perform the iteration for the first row in the group + uint64_t row = *groupStartIt; ValueType xBest, yBest; multiplyRow(row, A, b[row], xBest, yBest); ++row; - // Only do more work if there are still rows in this row group if (row != groupEnd) { ValueType xi, yi; uint64_t xyTmpIndex = 0; - if (hasPrimaryBound()) { - ValueType bestValue = xBest + yBest * getPrimaryBound(); + if (hasUpperBound) { + ValueType bestValue = xBest + yBest * upperBound; for (;row < groupEnd; ++row) { // Get the multiplication results multiplyRow(row, A, b[row], xi, yi); - ValueType currentValue = xi + yi * getPrimaryBound(); + ValueType currentValue = xi + yi * upperBound; // Check if the current row is better then the previously found one - if (better(currentValue, bestValue)) { + if (currentValue > bestValue) { if (yBest < yi) { // We need to store the 'old' best value as it might be relevant for the decision value xTmp[xyTmpIndex] = std::move(xBest); @@ -849,7 +933,7 @@ namespace storm { for (;row < groupEnd; ++row) { multiplyRow(row, A, b[row], xi, yi); // Update the best choice - if (yi > yBest || (yi == yBest && better(xi, xBest))) { + if (yi > yBest || (yi == yBest && xi > xBest)) { xTmp[xyTmpIndex] = std::move(xBest); yTmp[xyTmpIndex] = std::move(yBest); ++xyTmpIndex; @@ -868,15 +952,95 @@ namespace storm { ValueType deltaY = yBest - yTmp[i]; if (deltaY > storm::utility::zero()) { ValueType newDecisionValue = (xTmp[i] - xBest) / deltaY; - if (!hasDecisionValue || better(newDecisionValue, decisionValue)) { + if (!hasDecisionValue || newDecisionValue > decisionValue) { decisionValue = std::move(newDecisionValue); hasDecisionValue = true; } } } } - x[group] = std::move(xBest); - y[group] = std::move(yBest); + *xIt = std::move(xBest); + *yIt = std::move(yBest); + } + } + + void performIterationStepUpdateDecisionValueMin(storm::storage::SparseMatrix const& A, std::vector const& b) { + auto xIt = x.rbegin(); + auto yIt = y.rbegin(); + auto groupStartIt = A.getRowGroupIndices().rbegin(); + uint64_t groupEnd = *groupStartIt; + ++groupStartIt; + for (auto groupStartIte = A.getRowGroupIndices().rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + // Perform the iteration for the first row in the group + uint64_t row = *groupStartIt; + ValueType xBest, yBest; + multiplyRow(row, A, b[row], xBest, yBest); + ++row; + // Only do more work if there are still rows in this row group + if (row != groupEnd) { + ValueType xi, yi; + uint64_t xyTmpIndex = 0; + if (hasLowerBound) { + ValueType bestValue = xBest + yBest * lowerBound; + for (;row < groupEnd; ++row) { + // Get the multiplication results + multiplyRow(row, A, b[row], xi, yi); + ValueType currentValue = xi + yi * lowerBound; + // Check if the current row is better then the previously found one + if (currentValue < bestValue) { + if (yBest < yi) { + // We need to store the 'old' best value as it might be relevant for the decision value + xTmp[xyTmpIndex] = std::move(xBest); + yTmp[xyTmpIndex] = std::move(yBest); + ++xyTmpIndex; + } + xBest = std::move(xi); + yBest = std::move(yi); + bestValue = std::move(currentValue); + } else if (yBest > yi) { + // If the value for this row is not strictly better, it might still be equal and have a better y value + if (currentValue == bestValue) { + xBest = std::move(xi); + yBest = std::move(yi); + } else { + xTmp[xyTmpIndex] = std::move(xi); + yTmp[xyTmpIndex] = std::move(yi); + ++xyTmpIndex; + } + } + } + } else { + for (;row < groupEnd; ++row) { + multiplyRow(row, A, b[row], xi, yi); + // Update the best choice + if (yi > yBest || (yi == yBest && xi < xBest)) { + xTmp[xyTmpIndex] = std::move(xBest); + yTmp[xyTmpIndex] = std::move(yBest); + ++xyTmpIndex; + xBest = std::move(xi); + yBest = std::move(yi); + } else { + xTmp[xyTmpIndex] = std::move(xi); + yTmp[xyTmpIndex] = std::move(yi); + ++xyTmpIndex; + } + } + } + + // Update the decision value + for (uint64_t i = 0; i < xyTmpIndex; ++i) { + ValueType deltaY = yBest - yTmp[i]; + if (deltaY > storm::utility::zero()) { + ValueType newDecisionValue = (xTmp[i] - xBest) / deltaY; + if (!hasDecisionValue || newDecisionValue < decisionValue) { + decisionValue = std::move(newDecisionValue); + hasDecisionValue = true; + } + } + } + } + *xIt = std::move(xBest); + *yIt = std::move(yBest); } } @@ -1086,13 +1250,13 @@ namespace storm { while (status == SolverStatus::InProgress && iterations < env.solver().minMax().getMaximalNumberOfIterations()) { if (minimize(dir)) { - helper.template performIterationStep(*this->A, b); + helper.template performIterationStepMin(*this->A, b); if (helper.template checkConvergenceUpdateBounds(iterations, relevantValuesPtr)) { status = SolverStatus::Converged; } } else { assert(maximize(dir)); - helper.template performIterationStep(*this->A, b); + helper.template performIterationStepMax(*this->A, b); if (helper.template checkConvergenceUpdateBounds(iterations, relevantValuesPtr)) { status = SolverStatus::Converged; }