diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 45b7479b8..25583a654 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -217,12 +217,12 @@ namespace storm { solver->setMatrix(std::move(submatrix)); } - // Potentially show progress. - this->showProgressIterative(iterations); - // Update environment variables. ++iterations; status = updateStatusIfNotConverged(status, x, iterations, dir == storm::OptimizationDirection::Minimize ? SolverGuarantee::GreaterOrEqual : SolverGuarantee::LessOrEqual); + + // Potentially show progress. + this->showProgressIterative(iterations); } while (status == Status::InProgress); reportStatus(status, iterations); @@ -263,6 +263,10 @@ namespace storm { // Start by copying the requirements of the linear equation solver. MinMaxLinearEquationSolverRequirements requirements(this->linearEquationSolverFactory->getRequirements()); + if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + requirements.requireLowerBounds(); + } + // Guide requirements by whether or not we force soundness. if (this->getSettings().getForceSoundness()) { // Only add requirements for value iteration here as the policy iteration requirements are indifferent @@ -277,6 +281,8 @@ namespace storm { } } } + + requirements.requireBounds(); } // 'Regular' requirements (even for non-sound solving techniques). @@ -298,6 +304,8 @@ namespace storm { template typename IterativeMinMaxLinearEquationSolver::ValueIterationResult IterativeMinMaxLinearEquationSolver::performValueIteration(OptimizationDirection dir, std::vector*& currentX, std::vector*& newX, std::vector const& b, ValueType const& precision, bool relative, SolverGuarantee const& guarantee, uint64_t currentIterations) const { + STORM_LOG_ASSERT(currentX != newX, "Vectors must not be aliased."); + // Get handle to linear equation solver. storm::solver::LinearEquationSolver const& linearEquationSolver = *this->linEqSolverA; @@ -320,18 +328,20 @@ namespace storm { linearEquationSolver.multiplyAndReduce(dir, this->A->getRowGroupIndices(), *currentX, &b, *newX); } + std::cout << "position 4: " << (*currentX)[4] << " vs " << (*newX)[4] << std::endl; + // Determine whether the method converged. if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, precision, relative)) { status = Status::Converged; } - // Potentially show progress. - this->showProgressIterative(iterations); - // Update environment variables. std::swap(currentX, newX); ++iterations; status = updateStatusIfNotConverged(status, *currentX, iterations, guarantee); + + // Potentially show progress. + this->showProgressIterative(iterations); } // Swap the pointers so that the output is always in currentX. @@ -570,9 +580,6 @@ namespace storm { } } - // Potentially show progress. - this->showProgressIterative(iterations); - // Update environment variables. ++iterations; doConvergenceCheck = !doConvergenceCheck; @@ -582,6 +589,9 @@ namespace storm { if (upperStep) { status = updateStatusIfNotConverged(status, *upperX, iterations, SolverGuarantee::GreaterOrEqual); } + + // Potentially show progress. + this->showProgressIterative(iterations); } reportStatus(status, iterations); @@ -634,10 +644,10 @@ namespace storm { // If the value does not match the one in the values vector, the given vector is not a solution. if (!comparator.isEqual(groupValue, *valueIt)) { - std::cout << "group value for " << group << " is " << groupValue << " is not " << *valueIt << ", diff is " << (groupValue - *valueIt) << std::endl; + std::cout << "offending group " << group << " : " << groupValue << " vs " << *valueIt << std::endl; return false; - } else { - std::cout << "group value for " << group << " is " << groupValue << " is actually " << *valueIt << std::endl; + } else if (group == 1) { + std::cout << "val for group 1 is " << groupValue << " vs 227630345357/3221225472" << std::endl; } } @@ -649,12 +659,18 @@ namespace storm { template bool IterativeMinMaxLinearEquationSolver::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp) { - for (uint64_t p = 1; p < precision; ++p) { - storm::utility::kwek_mehlhorn::sharpen(precision, x, tmp); + std::cout << "pre sharpen x[0] " << x[0] << ", is smaller? " << (x[0] < storm::utility::convertNumber(std::string("227630345357/3221225472"))) << " with diff " << (x[0] - storm::utility::convertNumber(std::string("227630345357/3221225472"))) << std::endl; + + storm::utility::kwek_mehlhorn::sharpen(precision, x, tmp); - if (IterativeMinMaxLinearEquationSolver::isSolution(dir, A, tmp, b)) { - return true; - } + std::cout << "post sharpen x[0] " << tmp[0] << ", is smaller? " << (tmp[0] < storm::utility::convertNumber(std::string("227630345357/3221225472"))) << " with diff " << (tmp[0] - storm::utility::convertNumber(std::string("227630345357/3221225472"))) << std::endl; + + if (std::is_same::value) { + std::cout << "sharpen smaller? " << (tmp[0] < storm::utility::convertNumber(x[0])) << std::endl; + } + + if (IterativeMinMaxLinearEquationSolver::isSolution(dir, A, tmp, b)) { + return true; } return false; } @@ -769,10 +785,7 @@ namespace storm { auto targetIt = auxiliaryRowGroupVector->begin(); for (auto it = impreciseX.begin(), ite = impreciseX.end(); it != ite; ++it, ++targetIt) { *targetIt = storm::utility::convertNumber(*it); - - std::cout << "converting vector[" << std::distance(impreciseX.begin(), it) << "] from " << *it << " to " << *targetIt << std::endl; } - // Get rid of the superfluous data structures. impreciseX = std::vector(); @@ -810,11 +823,7 @@ namespace storm { template struct TemporaryHelper { static std::vector* getFreeRational(std::vector& rationalX, std::vector*& currentX, std::vector*& newX) { - if (&rationalX == currentX) { - return newX; - } else { - return currentX; - } + return newX; } static void swapSolutions(std::vector& rationalX, std::vector*& rationalSolution, std::vector& x, std::vector*& currentX, std::vector*& newX) { @@ -851,7 +860,7 @@ namespace storm { // At this point, the result of the imprecise value iteration is stored in the (imprecise) current x. ++valueIterationInvocations; - STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << "."); + STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations."); // Count the iterations. overallIterations += result.iterations; @@ -865,7 +874,7 @@ namespace storm { // Sharpen solution and place it in rationalX. bool foundSolution = sharpen(dir, p, rationalA, *currentX, rationalB, *freeRational); - // After sharpen, if a solution was found, it is contained in the current rational x. + // After sharpen, if a solution was found, it is contained in the free rational. if (foundSolution) { status = Status::Converged; @@ -873,7 +882,7 @@ namespace storm { TemporaryHelper::swapSolutions(rationalX, freeRational, x, currentX, newX); } else { // Increase the precision. - precision = precision / 10; + precision = precision / 100; } } diff --git a/src/storm/utility/KwekMehlhorn.h b/src/storm/utility/KwekMehlhorn.h index 30c48cba3..13c976fed 100644 --- a/src/storm/utility/KwekMehlhorn.h +++ b/src/storm/utility/KwekMehlhorn.h @@ -29,7 +29,7 @@ namespace storm { template std::pair::IntegerType, typename NumberTraits::IntegerType> truncateToRational(ImpreciseType const& value, uint64_t precision) { typedef typename NumberTraits::IntegerType IntegerType; - + IntegerType powerOfTen = storm::utility::pow(storm::utility::convertNumber(10ull), precision); IntegerType truncated = storm::utility::trunc(value * powerOfTen); return std::make_pair(truncated, powerOfTen); @@ -38,6 +38,7 @@ namespace storm { template std::pair::IntegerType, typename NumberTraits::IntegerType> truncateToRational(double const& value, uint64_t precision) { STORM_LOG_THROW(precision < 17, storm::exceptions::PrecisionExceededException, "Exceeded precision of double, consider switching to rational numbers."); + double powerOfTen = std::pow(10, precision); double truncated = storm::utility::trunc(value * powerOfTen); return std::make_pair(truncated, powerOfTen); @@ -59,7 +60,13 @@ namespace storm { for (uint64_t index = 0; index < input.size(); ++index) { ImpreciseType integer = storm::utility::floor(input[index]); ImpreciseType fraction = input[index] - integer; - output[index] = storm::utility::convertNumber(integer) + findRational(precision, fraction); + auto rational = findRational(precision, fraction); + output[index] = storm::utility::convertNumber(integer) + rational; + STORM_LOG_ASSERT(storm::utility::isZero(fraction) || !storm::utility::isZero(rational), "Found zero rational for non-zero fraction " << fraction << "."); + STORM_LOG_ASSERT(rational >= storm::utility::zero(), "Expected non-negative rational."); + if (std::is_same::value) { + STORM_LOG_ASSERT(output[index] >= input[index], "Sharpen decreased value from " << input[index] << " to " << output[index] << "."); + } } } diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index 7b0bb1baf..dd9e3cae9 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -106,6 +106,11 @@ namespace storm { return static_cast(number); } + template<> + double convertNumber(std::string const& value){ + return carl::toDouble(carl::parse(value)); + } + template<> storm::storage::sparse::state_type convertNumber(long long const& number){ return static_cast(number); @@ -330,7 +335,7 @@ namespace storm { STORM_LOG_ASSERT(static_cast(number) == number, "Rationalizing failed, because the number is too large."); return carl::rationalize(static_cast(number)); } - + template<> double convertNumber(ClnRationalNumber const& number) { return carl::toDouble(number); @@ -544,8 +549,7 @@ namespace storm { template<> typename NumberTraits::IntegerType trunc(GmpRationalNumber const& number) { - // FIXME: precision issue. - return typename NumberTraits::IntegerType(static_cast(std::trunc(number.get_d()))); + return carl::getNum(number) / carl::getDenom(number); } template<> diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index 0bb468043..e322bb25d 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -797,7 +797,10 @@ namespace storm { } } else { T diff = val1 - val2; - if (storm::utility::abs(diff) > precision) return false; + if (storm::utility::abs(diff) > precision) { + std::cout << "diff " << storm::utility::abs(diff) << " vs precision " << precision << std::endl; + return false; + } } return true; } @@ -815,7 +818,9 @@ namespace storm { } } else { double diff = val1 - val2; - if (storm::utility::abs(diff) > precision) return false; + if (storm::utility::abs(diff) > precision) { + return false; + } } return true; } @@ -833,8 +838,12 @@ namespace storm { bool equalModuloPrecision(std::vector const& vectorLeft, std::vector const& vectorRight, T const& precision, bool relativeError) { STORM_LOG_ASSERT(vectorLeft.size() == vectorRight.size(), "Lengths of vectors does not match."); - for (uint_fast64_t i = 0; i < vectorLeft.size(); ++i) { - if (!equalModuloPrecision(vectorLeft[i], vectorRight[i], precision, relativeError)) { + auto leftIt = vectorLeft.begin(); + auto leftIte = vectorLeft.end(); + auto rightIt = vectorRight.begin(); + for (; leftIt != leftIte; ++leftIt, ++rightIt) { + if (!equalModuloPrecision(*leftIt, *rightIt, precision, relativeError)) { + std::cout << "offending position " << std::distance(vectorLeft.begin(), leftIt) << std::endl; return false; } }