From 1f16008b75d426dc3acff10a697ec90c0f39af98 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Thu, 12 Oct 2017 19:59:34 +0200
Subject: [PATCH] added proper exception handling to sylvan-based sharpening
---
.../3rdparty/sylvan/src/storm_wrapper.cpp | 13 +++--
resources/3rdparty/sylvan/src/sylvan_mtbdd.c | 54 +++++++++++++++++++
resources/3rdparty/sylvan/src/sylvan_mtbdd.h | 8 +++
.../3rdparty/sylvan/src/sylvan_mtbdd_storm.c | 10 +++-
.../3rdparty/sylvan/src/sylvan_obj_storm.cpp | 8 ++-
.../solver/EigenLinearEquationSolver.cpp | 2 -
.../IterativeMinMaxLinearEquationSolver.cpp | 2 +-
.../solver/NativeLinearEquationSolver.cpp | 25 +--------
.../SymbolicMinMaxLinearEquationSolver.cpp | 2 +-
.../SymbolicNativeLinearEquationSolver.cpp | 2 +-
src/storm/storage/SparseMatrix.cpp | 32 +++++------
11 files changed, 109 insertions(+), 49 deletions(-)
diff --git a/resources/3rdparty/sylvan/src/storm_wrapper.cpp b/resources/3rdparty/sylvan/src/storm_wrapper.cpp
index 287b56f3e..33316bd0b 100644
--- a/resources/3rdparty/sylvan/src/storm_wrapper.cpp
+++ b/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;
}
diff --git a/resources/3rdparty/sylvan/src/sylvan_mtbdd.c b/resources/3rdparty/sylvan/src/sylvan_mtbdd.c
index 1eef5fa19..964b538be 100755
--- a/resources/3rdparty/sylvan/src/sylvan_mtbdd.c
+++ b/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;
diff --git a/resources/3rdparty/sylvan/src/sylvan_mtbdd.h b/resources/3rdparty/sylvan/src/sylvan_mtbdd.h
index 41cb2528e..facb9618b 100755
--- a/resources/3rdparty/sylvan/src/sylvan_mtbdd.h
+++ b/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).
diff --git a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c
index 55c78e4f5..011bb4868 100644
--- a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.c
+++ b/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)
diff --git a/resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp b/resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp
index 974e60a18..efd50edb5 100644
--- a/resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp
+++ b/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
diff --git a/src/storm/solver/EigenLinearEquationSolver.cpp b/src/storm/solver/EigenLinearEquationSolver.cpp
index 91b0d7bda..c285423e1 100644
--- a/src/storm/solver/EigenLinearEquationSolver.cpp
+++ b/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();
}
diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
index 56a7b0189..696425cbc 100644
--- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
+++ b/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;
}
}
diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp
index d17a2eb4d..e77d23d94 100644
--- a/src/storm/solver/NativeLinearEquationSolver.cpp
+++ b/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;
}
}
diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
index e2d8f931b..c6c6ec261 100644
--- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
@@ -214,7 +214,7 @@ namespace storm {
if (isSolution) {
status = SolverStatus::Converged;
} else {
- precision = precision / 100;
+ precision /= 10;
}
}
diff --git a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp b/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
index 209623ef2..21d06a678 100644
--- a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
@@ -253,7 +253,7 @@ namespace storm {
if (isSolution) {
status = SolverStatus::Converged;
} else {
- precision = precision / 100;
+ precision /= 10;
}
}
diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp
index 643802b56..c8f53a342 100644
--- a/src/storm/storage/SparseMatrix.cpp
+++ b/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;
}
}