added proper exception handling to sylvan-based sharpening

resources/3rdparty/sylvan/src/storm_wrapper.cpp

@@ -11,6 +11,7 @@
 #include "storm/utility/constants.h"
 #include "storm/utility/KwekMehlhorn.h"
 #include "storm/exceptions/InvalidOperationException.h"
+#include "storm/exceptions/PrecisionExceededException.h"
 
 #include <sylvan_config.h>
 #include <sylvan.h>
@@ -315,8 +316,13 @@ storm_rational_number_ptr storm_double_sharpen(double value, size_t precision) {
     std::lock_guard<std::mutex> lock(rationalNumberMutex);
 #endif
 
-    storm::RationalNumber* result_srn = new storm::RationalNumber(storm::utility::kwek_mehlhorn::sharpen<storm::RationalNumber, double>(precision, value));
-    return (storm_rational_number_ptr)result_srn;
+    try {
+        storm::RationalNumber tmp = storm::utility::kwek_mehlhorn::sharpen<storm::RationalNumber, double>(precision, value);
+        storm::RationalNumber* result_srn = new storm::RationalNumber(tmp);
+        return (storm_rational_number_ptr)result_srn;
+    } catch (storm::exceptions::PrecisionExceededException const& e) {
+        return (storm_rational_number_ptr)0;
+    }
 }
 
 storm_rational_number_ptr storm_rational_number_sharpen(storm_rational_number_ptr a, size_t precision) {
@@ -325,7 +331,8 @@ storm_rational_number_ptr storm_rational_number_sharpen(storm_rational_number_pt
 #endif
 
     storm::RationalNumber const& srn_a = *(storm::RationalNumber const*)a;
-    storm::RationalNumber* result_srn = new storm::RationalNumber(storm::utility::kwek_mehlhorn::sharpen<storm::RationalNumber, storm::RationalNumber>(precision, srn_a));
+    storm::RationalNumber tmp = storm::utility::kwek_mehlhorn::sharpen<storm::RationalNumber, storm::RationalNumber>(precision, srn_a);
+    storm::RationalNumber* result_srn = new storm::RationalNumber(tmp);
     return (storm_rational_number_ptr)result_srn;
 }
 

resources/3rdparty/sylvan/src/sylvan_mtbdd.c

@@ -900,6 +900,60 @@ TASK_IMPL_3(MTBDD, mtbdd_uapply, MTBDD, dd, mtbdd_uapply_op, op, size_t, param)
     return result;
 }
 
+/**
+ * Apply a unary operation <op> to <dd>.
+ */
+TASK_IMPL_3(MTBDD, mtbdd_uapply_fail_false, MTBDD, dd, mtbdd_uapply_op, op, size_t, param)
+{
+    /* Maybe perform garbage collection */
+    sylvan_gc_test();
+    
+    /* Count operation */
+    sylvan_stats_count(MTBDD_UAPPLY);
+    
+    /* Check cache */
+    MTBDD result;
+    if (cache_get3(CACHE_MTBDD_UAPPLY, dd, (size_t)op, param, &result)) {
+        sylvan_stats_count(MTBDD_UAPPLY_CACHED);
+        return result;
+    }
+    
+    /* Check terminal case */
+    result = WRAP(op, dd, param);
+    if (result != mtbdd_invalid) {
+        /* Store in cache */
+        if (cache_put3(CACHE_MTBDD_UAPPLY, dd, (size_t)op, param, result)) {
+            sylvan_stats_count(MTBDD_UAPPLY_CACHEDPUT);
+        }
+        
+        return result;
+    }
+    
+    /* Get cofactors */
+    mtbddnode_t ndd = MTBDD_GETNODE(dd);
+    MTBDD ddlow = node_getlow(dd, ndd);
+    MTBDD ddhigh = node_gethigh(dd, ndd);
+    
+    /* Recursive */
+    mtbdd_refs_spawn(SPAWN(mtbdd_uapply, ddhigh, op, param));
+    MTBDD low = mtbdd_refs_push(CALL(mtbdd_uapply, ddlow, op, param));
+    MTBDD high = mtbdd_refs_sync(SYNC(mtbdd_uapply));
+    mtbdd_refs_pop(1);
+    
+    if (low == mtbdd_false || high == mtbdd_false) {
+        result = mtbdd_false;
+    } else {
+        result = mtbdd_makenode(mtbddnode_getvariable(ndd), low, high);
+    }
+    
+    /* Store in cache */
+    if (cache_put3(CACHE_MTBDD_UAPPLY, dd, (size_t)op, param, result)) {
+        sylvan_stats_count(MTBDD_UAPPLY_CACHEDPUT);
+    }
+    
+    return result;
+}
+
 TASK_2(MTBDD, mtbdd_uop_times_uint, MTBDD, a, size_t, k)
 {
     if (a == mtbdd_false) return mtbdd_false;

resources/3rdparty/sylvan/src/sylvan_mtbdd.h

@@ -420,6 +420,14 @@ TASK_DECL_5(MTBDD, mtbdd_applyp, MTBDD, MTBDD, size_t, mtbdd_applyp_op, uint64_t
 TASK_DECL_3(MTBDD, mtbdd_uapply, MTBDD, mtbdd_uapply_op, size_t);
 #define mtbdd_uapply(dd, op, param) CALL(mtbdd_uapply, dd, op, param)
 
+/**
+ * Apply a unary operation <op> to <dd> and fail if one of the subresults is mtbdd_false, where failing means returning
+ * mtbdd_false.
+ * Callback <op> is consulted after the cache, thus the application to a terminal is cached.
+ */
+TASK_DECL_3(MTBDD, mtbdd_uapply_fail_false, MTBDD, mtbdd_uapply_op, size_t);
+#define mtbdd_uapply_fail_false(dd, op, param) CALL(mtbdd_uapply_fail_false, dd, op, param)
+
 /**
  * Callback function types for abstraction.
  * MTBDD mtbdd_abstract_op(MTBDD a, MTBDD b, int k).

resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c

@@ -603,7 +603,13 @@ TASK_IMPL_2(MTBDD, mtbdd_op_sharpen, MTBDD, a, size_t, p)
     
     if (mtbddnode_isleaf(na)) {
         if (mtbddnode_gettype(na) == 1) {
-            MTBDD result = mtbdd_storm_rational_number(storm_double_sharpen(mtbdd_getdouble(a), p));
+            storm_rational_number_ptr rnp = storm_double_sharpen(mtbdd_getdouble(a), p);
+            
+            // If the sharpening failed, we return mtbdd_false so this can be detected at the top level.
+            if (rnp == (storm_rational_number_ptr)0) {
+                return mtbdd_false;
+            }
+            MTBDD result = mtbdd_storm_rational_number(rnp);
             return result;
         } else if (mtbddnode_gettype(na) == srn_type) {
             return mtbdd_storm_rational_number(storm_rational_number_sharpen((storm_rational_number_ptr)mtbdd_getstorm_rational_number_ptr(a), p));
@@ -618,7 +624,7 @@ TASK_IMPL_2(MTBDD, mtbdd_op_sharpen, MTBDD, a, size_t, p)
 
 TASK_IMPL_2(MTBDD, mtbdd_sharpen, MTBDD, dd, size_t, p)
 {
-    return mtbdd_uapply(dd, TASK(mtbdd_op_sharpen), p);
+    return mtbdd_uapply_fail_false(dd, TASK(mtbdd_op_sharpen), p);
 }
 
 TASK_IMPL_2(MTBDD, mtbdd_op_to_rational_number, MTBDD, a, size_t, p)

resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp

@@ -3,6 +3,8 @@
 #include "sylvan_storm_rational_number.h"
 #include "sylvan_storm_rational_function.h"
 
+#include "storm/exceptions/PrecisionExceededException.h"
+
 /*********************************************
  Functions added to sylvan's Sylvan class.
  *********************************************/
@@ -200,7 +202,11 @@ Mtbdd::EqualNormRel(const Mtbdd& other, double epsilon) const {
 Mtbdd
 Mtbdd::SharpenKwekMehlhorn(size_t precision) const {
     LACE_ME;
-    return mtbdd_sharpen(mtbdd, precision);
+    MTBDD result = mtbdd_sharpen(mtbdd, precision);
+    if (result == mtbdd_false) {
+        throw storm::exceptions::PrecisionExceededException() << "Exceeded precision of double, consider switching to rational numbers.";
+    }
+    return result;
 }
 
 Mtbdd

src/storm/solver/EigenLinearEquationSolver.cpp

@@ -143,8 +143,6 @@ namespace storm {
         
         template<typename ValueType>
         void EigenLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
-            std::cout << A << std::endl;
-            
             eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A);
             this->clearCache();
         }

src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp

@@ -902,7 +902,7 @@ namespace storm {
                     TemporaryHelper<RationalType, ImpreciseType>::swapSolutions(rationalX, temporaryRational, x, currentX, newX);
                 } else {
                     // Increase the precision.
-                    precision = precision / 100;
+                    precision /= 10;
                 }
             }
             

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -146,15 +146,12 @@ namespace storm {
             localA.reset();
             this->A = &A;
             clearCache();
-            
-            std::cout << *this->A << std::endl;
         }
 
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
             localA = std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A));
             this->A = localA.get();
-            std::cout << *this->A << std::endl;
             clearCache();
         }
 
@@ -690,7 +687,7 @@ namespace storm {
             
             std::vector<ImpreciseType>* currentX = &x;
             std::vector<ImpreciseType>* newX = &tmpX;
-            
+
             SolverStatus status = SolverStatus::InProgress;
             uint64_t overallIterations = 0;
             uint64_t valueIterationInvocations = 0;
@@ -701,24 +698,6 @@ namespace storm {
                 typename NativeLinearEquationSolver<ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(currentX, newX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations);
                 
                 // At this point, the result of the imprecise value iteration is stored in the (imprecise) current x.
-                std::cout << "solution" << std::endl;
-                uint64_t pos = 0;
-                for (auto const& e : *currentX) {
-                    std::cout << "[" << pos << "] " << e << std::endl;
-                    ++pos;
-                }
-                std::cout << "rational b" << std::endl;
-                pos = 0;
-                for (auto const& e : rationalB) {
-                    std::cout << "[" << pos << "] " << e << std::endl;
-                    ++pos;
-                }
-                std::cout << "b" << std::endl;
-                pos = 0;
-                for (auto const& e : b) {
-                    std::cout << "[" << pos << "] " << e << std::endl;
-                    ++pos;
-                }
 
                 ++valueIterationInvocations;
                 STORM_LOG_TRACE("Completed " << valueIterationInvocations << " power iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations.");
@@ -743,7 +722,7 @@ namespace storm {
                     TemporaryHelper<RationalType, ImpreciseType>::swapSolutions(rationalX, temporaryRational, x, currentX, newX);
                 } else {
                     // Increase the precision.
-                    precision = precision / 100;
+                    precision = precision / 10;
                 }
             }
             

src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp

@@ -214,7 +214,7 @@ namespace storm {
                 if (isSolution) {
                     status = SolverStatus::Converged;
                 } else {
-                    precision = precision / 100;
+                    precision /= 10;
                 }
             }
             

src/storm/solver/SymbolicNativeLinearEquationSolver.cpp

@@ -253,7 +253,7 @@ namespace storm {
                 if (isSolution) {
                     status = SolverStatus::Converged;
                 } else {
-                    precision = precision / 100;
+                    precision /= 10;
                 }
             }
             

src/storm/storage/SparseMatrix.cpp

@@ -1377,17 +1377,19 @@ namespace storm {
                 summandIterator = summand->begin();
             }
             
-            for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) {
+            for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++summandIterator) {
+                ValueType newValue;
                 if (summand) {
-                    *resultIterator = *summandIterator;
-                    ++summandIterator;
+                    newValue = *summandIterator;
                 } else {
-                    *resultIterator = storm::utility::zero<ValueType>();
+                    newValue = storm::utility::zero<ValueType>();
                 }
                 
                 for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
-                    *resultIterator += it->getValue() * vector[it->getColumn()];
+                    newValue += it->getValue() * vector[it->getColumn()];
                 }
+                
+                *resultIterator = newValue;
             }
         }
         
@@ -1404,20 +1406,18 @@ namespace storm {
             }
             
             for (; resultIterator != resultIteratorEnd; --rowIterator, --resultIterator, --summandIterator) {
+                ValueType newValue;
                 if (summand) {
-                    *resultIterator = *summandIterator;
-                    std::cout << "row[" << std::distance(rowIndications.begin(), rowIterator) << "]: " << *resultIterator << " because of summand" << std::endl;
+                    newValue = *summandIterator;
                 } else {
-                    *resultIterator = storm::utility::zero<ValueType>();
+                    newValue = storm::utility::zero<ValueType>();
                 }
                 
                 for (ite = this->begin() + *rowIterator - 1; it != ite; --it) {
-                    std::cout << "row[" << std::distance(rowIndications.begin(), rowIterator) << "]: " << *resultIterator << std::endl;
-                    std::cout << "row[" << std::distance(rowIndications.begin(), rowIterator) << "]: op " << *resultIterator << "  +  " << (it->getValue() * vector[it->getColumn()]) << " (= " << it->getValue() << " * " << vector[it->getColumn()] << ") = " << (*resultIterator + (it->getValue() * vector[it->getColumn()])) << std::endl;
-                    *resultIterator = *resultIterator + (it->getValue() * vector[it->getColumn()]);
-//                    std::cout << "row[" << std::distance(rowIndications.begin(), rowIterator) << "]: " << *resultIterator << " because of " << it->getValue() << " * " << vector[it->getColumn()] << " = " << (it->getValue() * vector[it->getColumn()]) << std::endl;
+                    newValue += (it->getValue() * vector[it->getColumn()]);
                 }
-                std::cout << "row[" << std::distance(rowIndications.begin(), rowIterator) << "] final value " << *resultIterator << std::endl;
+                
+                *resultIterator = newValue;
             }
         }
         
@@ -1447,11 +1447,13 @@ namespace storm {
                 }
                 
                 for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++summandIterator) {
-                    *resultIterator = summand ? *summandIterator : storm::utility::zero<ValueType>();
+                    ValueType newValue = summand ? *summandIterator : storm::utility::zero<ValueType>();
                     
                     for (ite = columnsAndEntries.begin() + *(rowIterator + 1); it != ite; ++it) {
-                        *resultIterator += it->getValue() * x[it->getColumn()];
+                        newValue += it->getValue() * x[it->getColumn()];
                     }
+                    
+                    *resultIterator = newValue;
                 }
             }