fixes related to relative termination criterion

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<storm::RationalNumber>(srn_a)) {
+            return storm::utility::isZero<storm::RationalNumber>(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;
     }

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,6 +445,11 @@ namespace storm {
                         }
                     }
                 }
+                STORM_LOG_ASSERT(lowerDiff >= storm::utility::zero<ValueType>(), "Expected non-negative lower diff.");
+                STORM_LOG_ASSERT(upperDiff >= storm::utility::zero<ValueType>(), "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<ValueType>(*lowerX, *upperX, precision, this->getSettings().getRelativeTerminationCriterion())) {

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) {

src/storm/storage/dd/BisimulationDecomposition.cpp

@@ -81,7 +81,7 @@ namespace storm {
                     auto durationSinceLastMessage = std::chrono::duration_cast<std::chrono::seconds>(now - timeOfLastMessage).count();
                     if (static_cast<uint64_t>(durationSinceLastMessage) >= showProgressDelay) {
                         auto durationSinceStart = std::chrono::duration_cast<std::chrono::seconds>(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();
                     }
                 }

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::milliseconds>(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;

src/storm/utility/constants.cpp

@@ -40,6 +40,10 @@ namespace storm {
             return a == zero<ValueType>();
         }
         
+        bool isAlmostZero(double const& a) {
+            return a < 1e-15;
+        }
+
         template<typename ValueType>
         bool isConstant(ValueType const&) {
             return true;

src/storm/utility/constants.h

@@ -40,6 +40,8 @@ namespace storm {
         template<typename ValueType>
         bool isZero(ValueType const& a);
 
+        bool isAlmostZero(double const& a);
+        
         template<typename ValueType>
         bool isConstant(ValueType const& a);
 

src/storm/utility/vector.h

@@ -792,10 +792,10 @@ namespace storm {
             template<class T>
             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<T>(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.