From e719a37c6cfcea3e901b0c662492b4850c242f12 Mon Sep 17 00:00:00 2001 From: dehnert Date: Thu, 21 Sep 2017 14:35:40 +0200 Subject: [PATCH] fixes related to relative termination criterion --- .../3rdparty/sylvan/src/storm_wrapper.cpp | 6 ++++- .../IterativeMinMaxLinearEquationSolver.cpp | 9 +++++-- .../solver/NativeLinearEquationSolver.cpp | 2 +- .../storage/dd/BisimulationDecomposition.cpp | 2 +- .../dd/bisimulation/PartitionRefiner.cpp | 4 ++-- src/storm/utility/constants.cpp | 4 ++++ src/storm/utility/constants.h | 2 ++ src/storm/utility/vector.h | 24 ++++++++++++++++--- 8 files changed, 43 insertions(+), 10 deletions(-) diff --git a/resources/3rdparty/sylvan/src/storm_wrapper.cpp b/resources/3rdparty/sylvan/src/storm_wrapper.cpp index fb92720d0..78d73c0c3 100644 --- a/resources/3rdparty/sylvan/src/storm_wrapper.cpp +++ b/resources/3rdparty/sylvan/src/storm_wrapper.cpp @@ -310,7 +310,11 @@ int storm_rational_number_equal_modulo_precision(int relative, storm_rational_nu storm::RationalNumber const& srn_p = *(storm::RationalNumber const*)precision; if (relative) { - return carl::abs(srn_a - srn_b)/srn_a < srn_p ? 1 : 0; + if (storm::utility::isZero(srn_a)) { + return storm::utility::isZero(srn_b); + } else { + return carl::abs(srn_a - srn_b)/srn_a < srn_p ? 1 : 0; + } } else { return carl::abs(srn_a - srn_b) < srn_p ? 1 : 0; } diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 7595b3812..2fb2c5d57 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -405,7 +405,7 @@ namespace storm { } while (status == Status::InProgress && iterations < this->getSettings().getMaximalNumberOfIterations()) { // In every thousandth iteration, we improve both bounds. - if (iterations % 1000 == 0) { + if (iterations % 1000 == 0 || lowerDiff == upperDiff) { if (useGaussSeidelMultiplication) { lowerDiff = (*lowerX)[0]; this->linEqSolverA->multiplyAndReduceGaussSeidel(dir, this->A->getRowGroupIndices(), *lowerX, &b); @@ -445,7 +445,12 @@ namespace storm { } } } - + STORM_LOG_ASSERT(lowerDiff >= storm::utility::zero(), "Expected non-negative lower diff."); + STORM_LOG_ASSERT(upperDiff >= storm::utility::zero(), "Expected non-negative upper diff."); + if (iterations % 1000 == 0) { + STORM_LOG_TRACE("Iteration " << iterations << ": lower difference: " << lowerDiff << ", upper difference: " << upperDiff << "."); + } + // Determine whether the method converged. if (storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, precision, this->getSettings().getRelativeTerminationCriterion())) { status = Status::Converged; diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 72b4a0b69..b59eb5c2f 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -458,7 +458,7 @@ namespace storm { bool upperStep = false; // In every thousandth iteration, we improve both bounds. - if (iterations % 1000 == 0) { + if (iterations % 1000 == 0 || lowerDiff == upperDiff) { lowerStep = true; upperStep = true; if (useGaussSeidelMultiplication) { diff --git a/src/storm/storage/dd/BisimulationDecomposition.cpp b/src/storm/storage/dd/BisimulationDecomposition.cpp index be19015ec..6efbd7f5a 100644 --- a/src/storm/storage/dd/BisimulationDecomposition.cpp +++ b/src/storm/storage/dd/BisimulationDecomposition.cpp @@ -81,7 +81,7 @@ namespace storm { auto durationSinceLastMessage = std::chrono::duration_cast(now - timeOfLastMessage).count(); if (static_cast(durationSinceLastMessage) >= showProgressDelay) { auto durationSinceStart = std::chrono::duration_cast(now - start).count(); - std::cout << "State partition after " << iterations << " iterations (" << durationSinceStart << "ms) has " << refiner->getStatePartition().getNumberOfBlocks() << " blocks." << std::endl; + STORM_LOG_INFO("State partition after " << iterations << " iterations (" << durationSinceStart << "ms) has " << refiner->getStatePartition().getNumberOfBlocks() << " blocks."); timeOfLastMessage = std::chrono::high_resolution_clock::now(); } } diff --git a/src/storm/storage/dd/bisimulation/PartitionRefiner.cpp b/src/storm/storage/dd/bisimulation/PartitionRefiner.cpp index e06d9db41..8d566d036 100644 --- a/src/storm/storage/dd/bisimulation/PartitionRefiner.cpp +++ b/src/storm/storage/dd/bisimulation/PartitionRefiner.cpp @@ -47,7 +47,7 @@ namespace storm { auto signature = signatureIterator.next(); auto signatureEnd = std::chrono::high_resolution_clock::now(); totalSignatureTime += (signatureEnd - signatureStart); - STORM_LOG_DEBUG("Signature " << refinements << "[" << index << "] DD has " << signature.getSignatureAdd().getNodeCount() << " nodes."); + STORM_LOG_TRACE("Signature " << refinements << "[" << index << "] DD has " << signature.getSignatureAdd().getNodeCount() << " nodes."); auto refinementStart = std::chrono::high_resolution_clock::now(); newPartition = signatureRefiner.refine(statePartition, signature); @@ -65,7 +65,7 @@ namespace storm { auto totalTimeInRefinement = std::chrono::duration_cast(std::chrono::high_resolution_clock::now() - start).count(); ++refinements; - STORM_LOG_DEBUG("Refinement " << refinements << " produced " << newPartition.getNumberOfBlocks() << " blocks and was completed in " << totalTimeInRefinement << "ms (signature: " << signatureTime << "ms, refinement: " << refinementTime << "ms)."); + STORM_LOG_TRACE("Refinement " << refinements << " produced " << newPartition.getNumberOfBlocks() << " blocks and was completed in " << totalTimeInRefinement << "ms (signature: " << signatureTime << "ms, refinement: " << refinementTime << "ms)."); return newPartition; } else { return oldPartition; diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index 67ac0a6b8..2b05ed4e4 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -39,6 +39,10 @@ namespace storm { bool isZero(ValueType const& a) { return a == zero(); } + + bool isAlmostZero(double const& a) { + return a < 1e-15; + } template bool isConstant(ValueType const&) { diff --git a/src/storm/utility/constants.h b/src/storm/utility/constants.h index 69d253e81..7669524ba 100644 --- a/src/storm/utility/constants.h +++ b/src/storm/utility/constants.h @@ -40,6 +40,8 @@ namespace storm { template bool isZero(ValueType const& a); + bool isAlmostZero(double const& a); + template bool isConstant(ValueType const& a); diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index 069c11759..028e3360c 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -792,10 +792,10 @@ namespace storm { template bool equalModuloPrecision(T const& val1, T const& val2, T const& precision, bool relativeError = true) { if (relativeError) { - if (val2 == 0) { - return (storm::utility::abs(val1) <= precision); + if (storm::utility::isZero(val1)) { + return storm::utility::isZero(val2); } - T relDiff = (val1 - val2)/val2; + T relDiff = (val1 - val2)/val1; if (storm::utility::abs(relDiff) > precision) { return false; } @@ -806,6 +806,24 @@ namespace storm { return true; } + // Specializiation for double as the relative check for doubles very close to zero is not meaningful. + template<> + bool equalModuloPrecision(double const& val1, double const& val2, double const& precision, bool relativeError) { + if (relativeError) { + if (storm::utility::isAlmostZero(val2)) { + return storm::utility::isAlmostZero(val1); + } + double relDiff = (val1 - val2)/val1; + if (storm::utility::abs(relDiff) > precision) { + return false; + } + } else { + double diff = val1 - val2; + if (storm::utility::abs(diff) > precision) return false; + } + return true; + } + /*! * Compares the two vectors and determines whether they are equal modulo the provided precision. Depending on whether the * flag is set, the difference between the vectors is computed relative to the value or in absolute terms.