From b0fd3c17308dc447637c587ea629082e852f3a13 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Thu, 30 Nov 2017 16:34:56 +0100
Subject: [PATCH 01/19] started to rework Fox-Glynn

---
 .../csl/helper/SparseCtmcCslHelper.cpp        |   4 +-
 src/storm/utility/numerical.h                 | 148 ++++++++++++++++++
 2 files changed, 151 insertions(+), 1 deletion(-)

diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
index 41b768058..2d270c720 100644
--- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
@@ -632,7 +632,9 @@ namespace storm {
                 }
                 
                 // Use Fox-Glynn to get the truncation points and the weights.
-                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
+//                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
+                
+                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
                 STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult));
                 
                 // Scale the weights so they add up to one.
diff --git a/src/storm/utility/numerical.h b/src/storm/utility/numerical.h
index bf9a533b4..78be4f6d6 100644
--- a/src/storm/utility/numerical.h
+++ b/src/storm/utility/numerical.h
@@ -2,10 +2,13 @@
 #include <tuple>
 #include <cmath>
 
+#include <boost/math/constants/constants.hpp>
+
 #include "storm/utility/macros.h"
 #include "storm/utility/constants.h"
 #include "storm/exceptions/InvalidArgumentException.h"
 #include "storm/exceptions/OutOfRangeException.h"
+#include "storm/exceptions/PrecisionExceededException.h"
 
 namespace storm {
     namespace utility {
@@ -125,7 +128,152 @@ namespace storm {
                     
                     return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights);
                 }
+ 
+            }
+
+            template<typename ValueType>
+            std::tuple<uint64_t, uint64_t, ValueType> foxGlynnFinder(ValueType lambda, ValueType underflow, ValueType overflow, ValueType accuracy) {
+                uint64_t m = static_cast<uint64_t>(lambda);
+                int64_t L = 0;
+                uint64_t R = 0;
+                uint64_t k = 4;
+                
+                ValueType pi = boost::math::constants::pi<ValueType>();
+
+                if (m < 25) {
+                    STORM_LOG_THROW(std::exp(-static_cast<double>(m)) >= underflow, storm::exceptions::PrecisionExceededException, "Underflow in Fox-Glynn.");
+                    L = 0;
+                } else {
+                    ValueType rhsb = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1/lambda) * std::exp(1 / (8 * lambda)));
+                    k = 4;
+                    
+                    do {
+                        L = m - std::ceil(k * std::sqrt(lambda) + 0.5);
+                        if (L <= 0) {
+                            break;
+                        }
+                        if (std::exp(-(k*k) / 2) / k <= rhsb) {
+                            break;
+                        }
+                        ++k;
+                    } while (true);
+                }
+                
+                uint64_t mmax;
+                uint64_t maxr;
+                ValueType r2l;
+                ValueType rhsa;
+
+                if (m < 400) {
+                    mmax = 400;
+                    r2l = std::sqrt(2 * 400);
+                    rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/400) * std::sqrt(2.0) * std::exp(1.0 / 16));
+                    maxr = 400 + static_cast<uint64_t>(std::ceil((400 + 1) / 2.0));
+                } else {
+                    mmax = m;
+                    r2l = std::sqrt(2 * lambda);
+                    rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/lambda) * std::sqrt(2.0) * std::exp(1.0 / 16));
+                    maxr = m + static_cast<uint64_t>(std::ceil((lambda + 1) / 2.0));
+                }
+                
+                k = 4;
+                do {
+                    R = maxr + static_cast<uint64_t>(std::ceil(k * r2l + 0.5));
+                    STORM_LOG_THROW(R <= maxr, storm::exceptions::PrecisionExceededException, "Fox-Glynn: cannot bound right tail.");
+                    ValueType dkl = storm::utility::one<ValueType>();
+                    if (dkl * std::exp(-(k * k / 2.0)) / k <= rhsa) {
+                        break;
+                    }
+                    ++k;
+                } while (true);
+                
+                ValueType wm = std::pow<ValueType>(10, -10) * overflow * (R - L + 1);
+
+                if (m >= 25) {
+                    ValueType lcm = -1 - (1.0/12*25) - std::log(2 * pi) - 0.5 * std::log(static_cast<ValueType>(m));
+                    uint64_t i = m - L;
+                    ValueType llb;
+                    if (i <= static_cast<uint64_t>(L)) {
+                        llb = -static_cast<ValueType>(i)*(i + 1)*(((2 * i + 1) / (6 * lambda)) + 0.5) * (1 / lambda) + lcm;
+                    } else {
+                        llb = std::max(i * std::log(1 - (i / (m + 1.0))), -lambda);
+                    }
+                    STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at lower bound.");
+                    if (m >= 400) {
+                        llb = lcm - std::pow(R - m + 1, 2) / (2 * lambda);
+                        STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at upper bound.");
+                    }
+                }
+                
+                return std::make_tuple(static_cast<uint64_t>(L), R, wm);
+            }
+            
+            template<typename ValueType>
+            std::tuple<uint64_t, uint64_t, ValueType, std::vector<ValueType>> foxGlynnWeighter(ValueType lambda, ValueType accuracy) {
+                
+                STORM_LOG_THROW(lambda > storm::utility::zero<ValueType>(), storm::exceptions::InvalidArgumentException, "Fox-Glynn: Lambda must be positive.");
+                
+                ValueType underflow = std::numeric_limits<ValueType>::min();
+                ValueType overflow = std::numeric_limits<ValueType>::max();
+
+                // Initialize.
+                uint64_t m = static_cast<uint64_t>(lambda);
+                std::tuple<uint64_t, uint64_t, ValueType> finderResult = foxGlynnFinder(lambda, underflow, overflow, accuracy);
+                uint64_t L = std::get<0>(finderResult);
+                uint64_t R = std::get<1>(finderResult);
+                ValueType wm = std::get<2>(finderResult);
+                std::vector<ValueType> w(R - L + 1);
+                w.back() = wm;
+
+                // Down.
+                for (uint64_t j = m; j > L; --j) {
+                    w[j - L - 1] = (j / lambda) * w[j - L];
+                }
+
+                // Up.
+                if (lambda < 400) {
+                    // Special.
+                    STORM_LOG_THROW(R <= 600, storm::exceptions::PrecisionExceededException, "Potential underflow in Fox-Glynn.");
+                    for (uint64_t j = m; j < R; ++j) {
+                        auto q = lambda / (j + 1);
+                        if (w[j - L] > underflow / q) {
+                            w[j - L + 1] = q * w[j - L];
+                        } else {
+                            R = j;
+                            break;
+                        }
+                    }
+                } else {
+                    for (uint64_t j = m; j < R; ++j) {
+                        w[j - L + 1] = (lambda / (j + 1)) * w[j - L];
+                    }
+                }
+
+                // Compute W.
+                ValueType W = storm::utility::zero<ValueType>();
+                uint64_t s = L;
+                uint64_t t = R;
+                
+                while (s < t) {
+                    if (w[s - L] <= w[t - L]) {
+                        W += w[s - L];
+                        ++s;
+                    } else {
+                        W += w[t - L];
+                        --t;
+                    }
+                }
+                W += w[s - L];
+                
+                // Return.
+                return std::make_tuple(L, R, W, w);
             }
+            
+            template<typename ValueType>
+            std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynn(ValueType lambda, ValueType accuracy) {
+                return foxGlynnWeighter(lambda, accuracy);
+            }
+                
         }
     }
 }

From 48b0a40d8a6a19f067547065fc74f52cd131a378 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Thu, 30 Nov 2017 16:44:58 +0100
Subject: [PATCH 02/19] fix typo

---
 src/storm/utility/numerical.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/storm/utility/numerical.h b/src/storm/utility/numerical.h
index 78be4f6d6..946bdd4a5 100644
--- a/src/storm/utility/numerical.h
+++ b/src/storm/utility/numerical.h
@@ -178,7 +178,7 @@ namespace storm {
                 
                 k = 4;
                 do {
-                    R = maxr + static_cast<uint64_t>(std::ceil(k * r2l + 0.5));
+                    R = mmax + static_cast<uint64_t>(std::ceil(k * r2l + 0.5));
                     STORM_LOG_THROW(R <= maxr, storm::exceptions::PrecisionExceededException, "Fox-Glynn: cannot bound right tail.");
                     ValueType dkl = storm::utility::one<ValueType>();
                     if (dkl * std::exp(-(k * k / 2.0)) / k <= rhsa) {

From 27558e2140045edbec3216bd78da02e46e8ecbe3 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Thu, 30 Nov 2017 22:23:11 +0100
Subject: [PATCH 03/19] started c++ifying David Jansen's implementation of
 Fox-Glynn

---
 src/storm/utility/numerical.cpp | 290 ++++++++++++++++++++++++++++++++
 src/storm/utility/numerical.h   | 272 +-----------------------------
 2 files changed, 298 insertions(+), 264 deletions(-)
 create mode 100644 src/storm/utility/numerical.cpp

diff --git a/src/storm/utility/numerical.cpp b/src/storm/utility/numerical.cpp
new file mode 100644
index 000000000..1143d8bfd
--- /dev/null
+++ b/src/storm/utility/numerical.cpp
@@ -0,0 +1,290 @@
+#include <vector>
+#include <tuple>
+#include <cmath>
+
+#include <boost/math/constants/constants.hpp>
+
+#include "storm/utility/macros.h"
+#include "storm/utility/constants.h"
+#include "storm/exceptions/InvalidArgumentException.h"
+#include "storm/exceptions/OutOfRangeException.h"
+#include "storm/exceptions/PrecisionExceededException.h"
+
+namespace storm {
+    namespace utility {
+        namespace numerical {
+
+            template<typename ValueType>
+            FoxGlynnResult<ValueType>::FoxGlynnResult() : left(0), right(0), totalWeight(0) {
+                // Intentionally left empty.
+            }
+            
+            template<typename ValueType>
+            FoxGlynnResult<ValueType> foxGlynnFinder(ValueType lambda, ValueType epsilon) {
+                int left, right;
+                FoxGlynn * pFG;
+                
+                /* tau is only used in underflow checks, which we are going to do in the
+                 logarithm domain. */
+                tau = log(tau);
+                /* In error bound comparisons, we always compare with epsilon*sqrt_2_pi.*/
+                epsilon *= sqrt_2_pi;
+                /****Compute pFG->left truncation point****/
+                if ( m < 25 )
+                {
+                    /* for lambda below 25 the exponential can be smaller than tau */
+                    /* if that is the case we expect underflows and warn the user */
+                    if ( -lambda <= tau )
+                    {
+                        printf("ERROR: Fox-Glynn: 0 < lambda < 25, underflow near Poi(%g,"
+                               "0) = %.2e. The results are UNRELIABLE.\n",
+                               lambda, exp(-lambda));
+                    }
+                    /* zero is used as left truncation point for lambda <= 25 */
+                    left = 0;
+                }
+                /* compute the left truncation point for lambda >= 25 */
+                /* for lambda < 25 we use zero as left truncation point */
+                else
+                {
+                    const double bl = (1 + 1 / lambda) * exp((1/lambda) * 0.125);
+                    const double sqrt_lambda = sqrt(lambda);
+                    int k;
+                    /*Start looking for the left truncation point*/
+                    /* start search at k=4 (taken from original Fox-Glynn paper) */
+                    /* increase the left truncation point until we fulfil the error
+                     condition */
+                    
+                    for ( k = 4 ; TRUE ; k++ ) {
+                        double max_err;
+                        
+                        left = m - (int) ceil(k*sqrt_lambda + 0.5);
+                        /* for small lambda the above calculation can yield negative truncation points, crop them here */
+                        if ( left <= 0 ) {
+                            left = 0;
+                            break;
+                        }
+                        /* Note that Propositions 2-4 in Fox--Glynn mix up notation: they
+                         write Phi where they mean 1 - Phi. (In Corollaries 1 and 2, phi is
+                         used correctly again.) */
+                        max_err = bl * exp(-0.5 * (k*k)) / k;
+                        if ( max_err * 2 <= epsilon ) {
+                            /* If the error on the left hand side is smaller, we can be
+                             more lenient on the right hand side. To this end, we now set
+                             epsilon to the part of the error that has not yet been eaten
+                             up by the left-hand truncation. */
+                            epsilon -= max_err;
+                            break;
+                        }
+                    }
+                    /*Finally the left truncation point is found*/
+                }
+                
+                /****Compute pFG->right truncation point****/
+                {
+                    double lambda_max;
+                    int m_max, k;
+                    /*According to Fox-Glynn, if lambda < 400 we should take lambda = 400,
+                     otherwise use the original value. This is for computing the right truncation point*/
+                    if ( m < 400 ) {
+                        lambda_max = lambda_400;
+                        m_max = 400;
+                        epsilon *= 0.662608824988162441697980;
+                        /* i.e. al = (1+1/400) * exp(1/16) * sqrt_2; epsilon /= al; */
+                    } else {
+                        lambda_max = lambda;
+                        m_max = m;
+                        epsilon *= (1 - 1 / (lambda + 1)) * 0.664265347050632847802225;
+                        /* i.e. al = (1+1/lambda) * exp(1/16) * sqrt_2; epsilon /= al; */
+                    }
+                    /* find right truncation point */
+                    
+                    /* This loop is a modification to the original Fox-Glynn paper.
+                     The search for the right truncation point is only terminated by
+                     the error condition and not by the stop index from the FG paper.
+                     This can yield more accurate results if necessary.*/
+                    for ( k = 4 ; TRUE ; k++ )
+                    {
+                        /* dkl_inv is between 1 - 1e-33 and 1 if lambda_max >= 400 and
+                         k >= 4; this will always be rounded to 1.0. We therefore leave the
+                         factor out.
+                         double dkl_inv=1 - exp(-266/401.0 * (k*sqrt(2*lambda_max) + 1.5));
+                         */
+                        if ( k * epsilon /* actually: "k * (dkl_inv*epsilon/al)", which
+                                          has been calculated above */  >= exp(-0.5*(k*k)) )
+                            break;
+                    }
+                    right = m_max + (int) ceil(k * sqrt(2 * lambda_max) + 0.5);
+                    if ( right > m_max + (int) ceil((lambda_max + 1) * 0.5) ) {
+                        printf("ERROR: Fox-Glynn: right = %d >> lambda = %g, cannot "
+                               "bound the right tail. The results are "
+                               "UNRELIABLE.\n", right, lambda_max);
+                    }
+                }
+                
+                /*Time to set the initial value for weights*/
+                pFG = calloc((size_t) 1, sizeof(FoxGlynn) +
+                             (right - left) * sizeof(pFG->weights[0]));
+                if ( NULL == pFG ) {
+                    err_msg_3(err_MEMORY, "finder(%d,%g,_,%g,_)", m, lambda, omega, NULL);
+                }
+                pFG->right = right;
+                pFG->left = left;
+                pFG->weights[m - left] = omega / ( 1.0e+10 * (right - left) );
+                
+                if ( m >= 25 )
+                {
+                    /* perform underflow check */
+                    double result, log_c_m_inf;
+                    int i;
+                    
+                    /* we are going to compare with tau - log(w[m]) */
+                    tau -= log(pFG->weights[m - left]);
+                    /*We take the c_m_inf = 0.14627 / sqrt( m ), as for lambda >= 25
+                     c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf*/
+                    /* Note that m-lambda is in the interval (-1,0],
+                     and -1/(12*m) is in [-1/(12*25),0).
+                     So, exp(m-lambda - 1/(12*m)) is in (exp(-1-1/(12*25)),exp(0)).
+                     Therefore, we can improve the lower bound on c_m to
+                     exp(-1-1/(12*25)) / sqrt(2*pi) = ~0.14627. Its logarithm is
+                     -1 - 1/(12*25) - log(2*pi) * 0.5 = ~ -1.922272 (rounded towards
+                     -infinity). */
+                    log_c_m_inf = -1.922272 - log((double) m) * 0.5;
+                    
+                    /* We use FG's Proposition 6 directly (and not Corollary 4 i and ii),
+                     as k_prime may be too large if pFG->left == 0. */
+                    i = m - left;
+                    if ( i <= left /* equivalent to 2*i <= m,
+                                    equivalent to i <= lambda/2 */ )
+                    {
+                        /* Use Proposition 6 (i). Note that Fox--Glynn are off by one in
+                         the proof of this proposition; they sum up to i-1, but should have
+                         summed up to i. */
+                        result = log_c_m_inf
+                        - i * (i+1) * (0.5 + (2*i+1)/(6*lambda)) / lambda;
+                    }
+                    else
+                    {
+                        /* Use Corollary 4 (iii). Note that k_prime <= sqrt(m+1)/m is a
+                         misprint for k_prime <= m/sqrt(m+1), which is equivalent to
+                         left >= 0, which holds trivially. */
+                        result = -lambda;
+                        if ( 0 != left ) {
+                            /* also use Proposition 6 (ii) */
+                            double result_1 = log_c_m_inf + i * log(1 - i/(double) (m+1));
+                            /*Take the maximum*/
+                            if ( result_1 > result )
+                                result = result_1;
+                        }
+                    }
+                    if ( result <= tau )
+                    {
+                        const int log10_result = (int) floor(result * log10_e);
+                        printf("ERROR: Fox-Glynn: lambda >= 25, underflow near Poi(%g,%d)"
+                               " <= %.2fe%+d. The results are UNRELIABLE.\n",
+                               lambda, left, exp(result - log10_result/log10_e),
+                               log10_result);
+                    }
+                    
+                    /*We still have to perform an underflow check for the right truncation point when lambda >= 400*/
+                    if ( m >= 400 )
+                    {
+                        /* Use Proposition 5 of Fox--Glynn */
+                        i = right - m;
+                        result = log_c_m_inf - i * (i + 1) / (2 * lambda);
+                        if ( result <= tau )
+                        {
+                            const int log10_result = (int) floor(result * log10_e);
+                            printf("ERROR: Fox-Glynn: lambda >= 400, underflow near "
+                                   "Poi(%g,%d) <= %.2fe%+d. The results are "
+                                   "UNRELIABLE.\n", lambda, right,
+                                   exp(result - log10_result / log10_e),
+                                   log10_result);
+                        }
+                    }
+                }
+                
+                return pFG;
+            }
+            
+            template<typename ValueType>
+            FoxGlynnResult<ValueType> foxGlynnWeighter(ValueType lambda, ValueType epsilon) {
+                /*The magic m point*/
+                const uint64_t m = (int) floor(lambda);
+                int j, t;
+                FoxGlynn * pFG;
+                
+                pFG = finder(m, lambda, tau, omega, epsilon);
+                if ( NULL == pFG ) {
+                    err_msg_4(err_CALLBY, "weighter(%g,%g,%g,%g)", lambda,
+                              tau, omega, epsilon, NULL);
+                }
+                
+                /*Fill the left side of the array*/
+                for ( j = m - pFG->left ; j > 0 ; j-- )
+                    pFG->weights[j-1] = (j+pFG->left) / lambda * pFG->weights[j];
+                
+                t = pFG->right - pFG->left;
+                /*Fill the right side of the array, have two cases lambda < 400 & lambda >= 400*/
+                if ( m < 400 )
+                {
+                    /*Perform the underflow check, according to Fox-Glynn*/
+                    if( pFG->right > 600 )
+                    {
+                        printf("ERROR: Fox-Glynn: pFG->right > 600, underflow is possible\n");
+                        freeFG(pFG);
+                        return NULL;
+                    }
+                    /*Compute weights*/
+                    for ( j = m - pFG->left ; j < t ; j++ )
+                    {
+                        double q = lambda / (j + 1 + pFG->left);
+                        if ( pFG->weights[j] > tau / q )
+                        {
+                            pFG->weights[j + 1] = q * pFG->weights[j];
+                        }else{
+                            t = j;
+                            pFG->right = j + pFG->left;
+                            break; /*It's time to compute W*/
+                        }
+                    }
+                }else{
+                    /*Compute weights*/
+                    for ( j = m - pFG->left ; j < t ; j++ )
+                        pFG->weights[j+1] = lambda / (j+1 + pFG->left) * pFG->weights[j];
+                }
+                
+                /*It is time to compute the normalization weight W*/
+                pFG->total_weight = 0.0;
+                j = 0;
+                /* t was set above */
+                
+                while( j < t )
+                {
+                    if ( pFG->weights[j] <= pFG->weights[t] )
+                    {
+                        pFG->total_weight += pFG->weights[j];
+                        j++;
+                    }else{
+                        pFG->total_weight += pFG->weights[t];
+                        t--;
+                    }
+                }
+                pFG->total_weight += pFG->weights[j];
+                
+                /* printf("Fox-Glynn: ltp = %d, rtp = %d, w = %10.15le \n", pFG->left, pFG->right, pFG->total_weight); */
+                
+                return pFG;
+            }
+            
+            template<typename ValueType>
+            FoxGlynnResult<ValueType> foxGlynn(ValueType lambda, ValueType epsilon) {
+                STORM_LOG_THROW(lambda > 0, storm::exceptions::InvalidArgumentException, "Fox-Glynn requires positive lambda.");
+                return foxGlynnWeighter(lambda, epsilon);
+            }
+
+            template FoxGlynnResult<double> foxGlynn(double lambda, double epsilon);
+            
+        }
+    }
+}
diff --git a/src/storm/utility/numerical.h b/src/storm/utility/numerical.h
index 946bdd4a5..233bb8a74 100644
--- a/src/storm/utility/numerical.h
+++ b/src/storm/utility/numerical.h
@@ -1,278 +1,22 @@
 #include <vector>
-#include <tuple>
-#include <cmath>
-
-#include <boost/math/constants/constants.hpp>
-
-#include "storm/utility/macros.h"
-#include "storm/utility/constants.h"
-#include "storm/exceptions/InvalidArgumentException.h"
-#include "storm/exceptions/OutOfRangeException.h"
-#include "storm/exceptions/PrecisionExceededException.h"
+#include <cstdint>
 
 namespace storm {
     namespace utility {
         namespace numerical {
-            template<typename ValueType>
-            std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> getFoxGlynnCutoff(ValueType lambda, ValueType overflow, ValueType accuracy) {
-                STORM_LOG_THROW(lambda != storm::utility::zero<ValueType>(), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero.");
-                
-                // This code is a modified version of the one in PRISM. According to their implementation, for lambda
-                // smaller than 400, we compute the result using the naive method.
-                if (lambda < 400) {
-                    ValueType eToPowerMinusLambda = std::exp(-lambda);
-                    ValueType targetValue = (1 - accuracy) / eToPowerMinusLambda;
-                    std::vector<ValueType> weights;
-                    
-                    ValueType exactlyKEvents = 1;
-                    ValueType atMostKEvents = exactlyKEvents;
-                    weights.push_back(exactlyKEvents * eToPowerMinusLambda);
-                    
-                    uint_fast64_t k = 1;
-                    do {
-                        exactlyKEvents *= lambda / k;
-                        atMostKEvents += exactlyKEvents;
-                        weights.push_back(exactlyKEvents * eToPowerMinusLambda);
-                        ++k;
-                    } while (atMostKEvents < targetValue);
-                    
-                    return std::make_tuple(0, k - 1, 1.0, weights);
-                } else {
-                    STORM_LOG_THROW(accuracy >= 1e-10, storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: the accuracy must not be below 1e-10.");
-                    
-                    // Factor from Fox&Glynn's paper. The paper does not explain where it comes from.
-                    ValueType factor = 1e+10;
-                    
-                    // Now start the Finder algorithm to find the truncation points.
-                    ValueType m = std::floor(lambda);
-                    uint_fast64_t leftTruncationPoint = 0, rightTruncationPoint = 0;
-                    {
-                        // Factors used by the corollaries explained in Fox & Glynns paper.
-                        // Square root of pi.
-                        ValueType sqrtpi = 1.77245385090551602729;
-                        
-                        // Square root of 2.
-                        ValueType sqrt2 = 1.41421356237309504880;
-                        
-                        // Set up a_\lambda, b_\lambda, and the square root of lambda.
-                        ValueType aLambda = 0, bLambda = 0, sqrtLambda = 0;
-                        if (m < 400) {
-                            sqrtLambda = std::sqrt(400.0);
-                            aLambda = (1.0 + 1.0 / 400.0) * std::exp(0.0625) * sqrt2;
-                            bLambda = (1.0 + 1.0 / 400.0) * std::exp(0.125 / 400.0);
-                        } else {
-                            sqrtLambda = std::sqrt(lambda);
-                            aLambda = (1.0 + 1.0 / lambda) * std::exp(0.0625) * sqrt2;
-                            bLambda = (1.0 + 1.0 / lambda) * std::exp(0.125 / lambda);
-                        }
-                        
-                        // Use Corollary 1 from the paper to find the right truncation point.
-                        uint_fast64_t k = 4;
-                        
-                        ValueType dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0) * (k * sqrt2 * sqrtLambda + 1.5)));
-                        
-                        // According to David Jansen the upper bound can be ignored to achieve more accurate results.
-                        // Right hand side of the equation in Corollary 1.
-                        while ((accuracy / 2.0) < (aLambda * dkl * std::exp(-(k*k / 2.0)) / (k * sqrt2 * sqrtpi))) {
-                            ++k;
-                            
-                            // d(k,Lambda) from the paper.
-                            dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0)*(k * sqrt2 * sqrtLambda + 1.5)));
-                        }
-                        
-                        // Left hand side of the equation in Corollary 1.
-                        rightTruncationPoint = static_cast<uint_fast64_t>(std::ceil((m + k * sqrt2 * sqrtLambda + 1.5)));
-                        
-                        // Use Corollary 2 to find left truncation point.
-                        k = 4;
-                        
-                        // Right hand side of the equation in Corollary 2.
-                        while ((accuracy / 2.0) < ((bLambda * std::exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) {
-                            ++k;
-                        }
-                        
-                        // Left hand side of the equation in Corollary 2.
-                        leftTruncationPoint = static_cast<uint_fast64_t>(std::max(std::trunc(m - k * sqrtLambda - 1.5), storm::utility::zero<ValueType>()));
-                        
-                        // Check for underflow.
-                        STORM_LOG_THROW(std::trunc((m - k * sqrtLambda - 1.5)) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point.");
-                    }
-                    
-                    std::vector<ValueType> weights(rightTruncationPoint - leftTruncationPoint + 1);
-                    weights[m - leftTruncationPoint] = overflow / (factor * (rightTruncationPoint - leftTruncationPoint));
-                    for (uint_fast64_t j = m; j > leftTruncationPoint; --j) {
-                        weights[j - 1 - leftTruncationPoint] = (j / lambda) * weights[j - leftTruncationPoint];
-                    }
-                    
-                    for (uint_fast64_t j = m; j < rightTruncationPoint; ++j) {
-                        weights[j + 1 - leftTruncationPoint] = (lambda / (j + 1)) * weights[j - leftTruncationPoint];
-                    }
-                    
-                    
-                    // Compute the total weight and start with smallest to avoid roundoff errors.
-                    ValueType totalWeight = storm::utility::zero<ValueType>();
-                    
-                    uint_fast64_t s = leftTruncationPoint;
-                    uint_fast64_t t = rightTruncationPoint;
-                    while (s < t) {
-                        if (weights[s - leftTruncationPoint] <= weights[t - leftTruncationPoint]) {
-                            totalWeight += weights[s - leftTruncationPoint];
-                            ++s;
-                        } else {
-                            totalWeight += weights[t - leftTruncationPoint];
-                            --t;
-                        }
-                    }
-                    
-                    totalWeight += weights[s - leftTruncationPoint];
-                    
-                    return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights);
-                }
- 
-            }
-
-            template<typename ValueType>
-            std::tuple<uint64_t, uint64_t, ValueType> foxGlynnFinder(ValueType lambda, ValueType underflow, ValueType overflow, ValueType accuracy) {
-                uint64_t m = static_cast<uint64_t>(lambda);
-                int64_t L = 0;
-                uint64_t R = 0;
-                uint64_t k = 4;
-                
-                ValueType pi = boost::math::constants::pi<ValueType>();
-
-                if (m < 25) {
-                    STORM_LOG_THROW(std::exp(-static_cast<double>(m)) >= underflow, storm::exceptions::PrecisionExceededException, "Underflow in Fox-Glynn.");
-                    L = 0;
-                } else {
-                    ValueType rhsb = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1/lambda) * std::exp(1 / (8 * lambda)));
-                    k = 4;
-                    
-                    do {
-                        L = m - std::ceil(k * std::sqrt(lambda) + 0.5);
-                        if (L <= 0) {
-                            break;
-                        }
-                        if (std::exp(-(k*k) / 2) / k <= rhsb) {
-                            break;
-                        }
-                        ++k;
-                    } while (true);
-                }
-                
-                uint64_t mmax;
-                uint64_t maxr;
-                ValueType r2l;
-                ValueType rhsa;
-
-                if (m < 400) {
-                    mmax = 400;
-                    r2l = std::sqrt(2 * 400);
-                    rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/400) * std::sqrt(2.0) * std::exp(1.0 / 16));
-                    maxr = 400 + static_cast<uint64_t>(std::ceil((400 + 1) / 2.0));
-                } else {
-                    mmax = m;
-                    r2l = std::sqrt(2 * lambda);
-                    rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/lambda) * std::sqrt(2.0) * std::exp(1.0 / 16));
-                    maxr = m + static_cast<uint64_t>(std::ceil((lambda + 1) / 2.0));
-                }
-                
-                k = 4;
-                do {
-                    R = mmax + static_cast<uint64_t>(std::ceil(k * r2l + 0.5));
-                    STORM_LOG_THROW(R <= maxr, storm::exceptions::PrecisionExceededException, "Fox-Glynn: cannot bound right tail.");
-                    ValueType dkl = storm::utility::one<ValueType>();
-                    if (dkl * std::exp(-(k * k / 2.0)) / k <= rhsa) {
-                        break;
-                    }
-                    ++k;
-                } while (true);
-                
-                ValueType wm = std::pow<ValueType>(10, -10) * overflow * (R - L + 1);
 
-                if (m >= 25) {
-                    ValueType lcm = -1 - (1.0/12*25) - std::log(2 * pi) - 0.5 * std::log(static_cast<ValueType>(m));
-                    uint64_t i = m - L;
-                    ValueType llb;
-                    if (i <= static_cast<uint64_t>(L)) {
-                        llb = -static_cast<ValueType>(i)*(i + 1)*(((2 * i + 1) / (6 * lambda)) + 0.5) * (1 / lambda) + lcm;
-                    } else {
-                        llb = std::max(i * std::log(1 - (i / (m + 1.0))), -lambda);
-                    }
-                    STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at lower bound.");
-                    if (m >= 400) {
-                        llb = lcm - std::pow(R - m + 1, 2) / (2 * lambda);
-                        STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at upper bound.");
-                    }
-                }
-                
-                return std::make_tuple(static_cast<uint64_t>(L), R, wm);
-            }
-            
             template<typename ValueType>
-            std::tuple<uint64_t, uint64_t, ValueType, std::vector<ValueType>> foxGlynnWeighter(ValueType lambda, ValueType accuracy) {
-                
-                STORM_LOG_THROW(lambda > storm::utility::zero<ValueType>(), storm::exceptions::InvalidArgumentException, "Fox-Glynn: Lambda must be positive.");
-                
-                ValueType underflow = std::numeric_limits<ValueType>::min();
-                ValueType overflow = std::numeric_limits<ValueType>::max();
+            struct FoxGlynnResult {
+                FoxGlynnResult();
 
-                // Initialize.
-                uint64_t m = static_cast<uint64_t>(lambda);
-                std::tuple<uint64_t, uint64_t, ValueType> finderResult = foxGlynnFinder(lambda, underflow, overflow, accuracy);
-                uint64_t L = std::get<0>(finderResult);
-                uint64_t R = std::get<1>(finderResult);
-                ValueType wm = std::get<2>(finderResult);
-                std::vector<ValueType> w(R - L + 1);
-                w.back() = wm;
-
-                // Down.
-                for (uint64_t j = m; j > L; --j) {
-                    w[j - L - 1] = (j / lambda) * w[j - L];
-                }
-
-                // Up.
-                if (lambda < 400) {
-                    // Special.
-                    STORM_LOG_THROW(R <= 600, storm::exceptions::PrecisionExceededException, "Potential underflow in Fox-Glynn.");
-                    for (uint64_t j = m; j < R; ++j) {
-                        auto q = lambda / (j + 1);
-                        if (w[j - L] > underflow / q) {
-                            w[j - L + 1] = q * w[j - L];
-                        } else {
-                            R = j;
-                            break;
-                        }
-                    }
-                } else {
-                    for (uint64_t j = m; j < R; ++j) {
-                        w[j - L + 1] = (lambda / (j + 1)) * w[j - L];
-                    }
-                }
-
-                // Compute W.
-                ValueType W = storm::utility::zero<ValueType>();
-                uint64_t s = L;
-                uint64_t t = R;
-                
-                while (s < t) {
-                    if (w[s - L] <= w[t - L]) {
-                        W += w[s - L];
-                        ++s;
-                    } else {
-                        W += w[t - L];
-                        --t;
-                    }
-                }
-                W += w[s - L];
-                
-                // Return.
-                return std::make_tuple(L, R, W, w);
+                uint64_t left;
+                uint64_t right;
+                ValueType totalWeight;
+                std::vector<ValueType> weights;
             }
             
             template<typename ValueType>
-            std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynn(ValueType lambda, ValueType accuracy) {
-                return foxGlynnWeighter(lambda, accuracy);
-            }
+            FoxGlynnResult<ValueType> foxGlynn(ValueType lambda, ValueType epsilon);
                 
         }
     }

From 70818dd9ddbb817c5b222ffa1eadc00aa6c7cf4e Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Fri, 1 Dec 2017 11:20:04 +0100
Subject: [PATCH 04/19] finished c++ifying David Jansen's implementation of
 Fox-Glynn

---
 .../csl/helper/SparseCtmcCslHelper.cpp        |  22 +-
 src/storm/utility/numerical.cpp               | 355 +++++++++---------
 src/storm/utility/numerical.h                 |   2 +-
 3 files changed, 182 insertions(+), 197 deletions(-)

diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
index 2d270c720..90d556133 100644
--- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
@@ -634,19 +634,19 @@ namespace storm {
                 // Use Fox-Glynn to get the truncation points and the weights.
 //                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
                 
-                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
-                STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult));
+                storm::utility::numerical::FoxGlynnResult<ValueType> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
+                STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << foxGlynnResult.left << ", right=" << foxGlynnResult.right);
                 
                 // Scale the weights so they add up to one.
-                for (auto& element : std::get<3>(foxGlynnResult)) {
-                    element /= std::get<2>(foxGlynnResult);
+                for (auto& element : foxGlynnResult.weights) {
+                    element /= foxGlynnResult.totalWeight;
                 }
                 
                 // If the cumulative reward is to be computed, we need to adjust the weights.
                 if (useMixedPoissonProbabilities) {
                     ValueType sum = storm::utility::zero<ValueType>();
                     
-                    for (auto& element : std::get<3>(foxGlynnResult)) {
+                    for (auto& element : foxGlynnResult.weights) {
                         sum += element;
                         element = (1 - sum) / uniformizationRate;
                     }
@@ -656,10 +656,10 @@ namespace storm {
                 
                 // Initialize result.
                 std::vector<ValueType> result;
-                uint_fast64_t startingIteration = std::get<0>(foxGlynnResult);
+                uint_fast64_t startingIteration = foxGlynnResult.left;
                 if (startingIteration == 0) {
                     result = values;
-                    storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]);
+                    storm::utility::vector::scaleVectorInPlace(result, foxGlynnResult.weights.front());
                     ++startingIteration;
                 } else {
                     if (useMixedPoissonProbabilities) {
@@ -674,9 +674,9 @@ namespace storm {
                 std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(env, std::move(uniformizedMatrix), storm::solver::LinearEquationSolverTask::Multiply);
                 solver->setCachingEnabled(true);
                 
-                if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) {
+                if (!useMixedPoissonProbabilities && foxGlynnResult.left > 1) {
                     // Perform the matrix-vector multiplications (without adding).
-                    solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1);
+                    solver->repeatedMultiply(values, addVector, foxGlynnResult.left - 1);
                 } else if (useMixedPoissonProbabilities) {
                     std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; };
                     
@@ -691,10 +691,10 @@ namespace storm {
                 // multiplication, scale and add the result.
                 ValueType weight = 0;
                 std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; };
-                for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) {
+                for (uint_fast64_t index = startingIteration; index <= foxGlynnResult.right; ++index) {
                     solver->repeatedMultiply(values, addVector, 1);
                     
-                    weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
+                    weight = foxGlynnResult.weights[index - foxGlynnResult.left];
                     storm::utility::vector::applyPointwise(result, values, result, addAndScale);
                 }
                 
diff --git a/src/storm/utility/numerical.cpp b/src/storm/utility/numerical.cpp
index 1143d8bfd..cec474ff7 100644
--- a/src/storm/utility/numerical.cpp
+++ b/src/storm/utility/numerical.cpp
@@ -1,13 +1,11 @@
-#include <vector>
-#include <tuple>
-#include <cmath>
+#include "storm/utility/numerical.h"
 
+#include <cmath>
 #include <boost/math/constants/constants.hpp>
 
 #include "storm/utility/macros.h"
 #include "storm/utility/constants.h"
 #include "storm/exceptions/InvalidArgumentException.h"
-#include "storm/exceptions/OutOfRangeException.h"
 #include "storm/exceptions/PrecisionExceededException.h"
 
 namespace storm {
@@ -15,79 +13,93 @@ namespace storm {
         namespace numerical {
 
             template<typename ValueType>
-            FoxGlynnResult<ValueType>::FoxGlynnResult() : left(0), right(0), totalWeight(0) {
+            FoxGlynnResult<ValueType>::FoxGlynnResult() : left(0), right(0), totalWeight(storm::utility::zero<ValueType>()) {
                 // Intentionally left empty.
             }
             
+            /*!
+             * The following implementation of Fox and Glynn's algorithm is taken from David Jansen's patched version
+             * in MRMC, which is based on his paper:
+             *
+             * https://pms.cs.ru.nl/iris-diglib/src/getContent.php?id=2011-Jansen-UnderstandingFoxGlynn
+             *
+             * We have only adapted the code to match more of C++'s and our coding guidelines.
+             */
+            
             template<typename ValueType>
             FoxGlynnResult<ValueType> foxGlynnFinder(ValueType lambda, ValueType epsilon) {
-                int left, right;
-                FoxGlynn * pFG;
+                ValueType tau = std::numeric_limits<ValueType>::min();
+                ValueType omega = std::numeric_limits<ValueType>::max();
+                ValueType const sqrt_2_pi = boost::math::constants::root_two_pi<ValueType>();
+                ValueType const log10_e = std::log10(boost::math::constants::e<ValueType>());
+                
+                uint64_t m = static_cast<uint64_t>(lambda);
+                
+                int64_t left = 0;
+                int64_t right = 0;
                 
-                /* tau is only used in underflow checks, which we are going to do in the
-                 logarithm domain. */
+                // tau is only used in underflow checks, which we are going to do in the logarithm domain.
                 tau = log(tau);
-                /* In error bound comparisons, we always compare with epsilon*sqrt_2_pi.*/
+                
+                // In error bound comparisons, we always compare with epsilon*sqrt_2_pi.
                 epsilon *= sqrt_2_pi;
-                /****Compute pFG->left truncation point****/
-                if ( m < 25 )
-                {
-                    /* for lambda below 25 the exponential can be smaller than tau */
-                    /* if that is the case we expect underflows and warn the user */
-                    if ( -lambda <= tau )
-                    {
-                        printf("ERROR: Fox-Glynn: 0 < lambda < 25, underflow near Poi(%g,"
-                               "0) = %.2e. The results are UNRELIABLE.\n",
-                               lambda, exp(-lambda));
+                
+                // Compute left truncation point.
+                if (m < 25) {
+                    // For lambda below 25 the exponential can be smaller than tau. If that is the case we expect
+                    // underflows and warn the user.
+                    if (-lambda <= tau) {
+                        STORM_LOG_WARN("Fox-Glynn: 0 < lambda < 25, underflow near Poi(" << lambda << ", 0) = " << std::exp(-lambda) << ". The results are unreliable.");
                     }
-                    /* zero is used as left truncation point for lambda <= 25 */
+                    
+                    // Zero is used as left truncation point for lambda <= 25.
                     left = 0;
-                }
-                /* compute the left truncation point for lambda >= 25 */
-                /* for lambda < 25 we use zero as left truncation point */
-                else
-                {
-                    const double bl = (1 + 1 / lambda) * exp((1/lambda) * 0.125);
-                    const double sqrt_lambda = sqrt(lambda);
-                    int k;
-                    /*Start looking for the left truncation point*/
-                    /* start search at k=4 (taken from original Fox-Glynn paper) */
-                    /* increase the left truncation point until we fulfil the error
-                     condition */
+                } else {
+                    // Compute the left truncation point for lambda >= 25 (for lambda < 25 we use zero as left truncation point).
+
+                    ValueType const bl = (1 + 1 / lambda) * std::exp((1/lambda) * 0.125);
+                    ValueType const sqrt_lambda = std::sqrt(lambda);
+                    int64_t k;
                     
-                    for ( k = 4 ; TRUE ; k++ ) {
-                        double max_err;
+                    // Start looking for the left truncation point:
+                    // * start search at k=4 (taken from original Fox-Glynn paper)
+                    // * increase the left truncation point until we fulfil the error condition
+                    
+                    for (k = 4;; ++k) {
+                        ValueType max_err;
+                        
+                        left = m - static_cast<int64_t>(std::ceil(k*sqrt_lambda + 0.5));
                         
-                        left = m - (int) ceil(k*sqrt_lambda + 0.5);
-                        /* for small lambda the above calculation can yield negative truncation points, crop them here */
-                        if ( left <= 0 ) {
+                        // For small lambda the above calculation can yield negative truncation points, crop them here.
+                        if (left <= 0) {
                             left = 0;
                             break;
                         }
-                        /* Note that Propositions 2-4 in Fox--Glynn mix up notation: they
-                         write Phi where they mean 1 - Phi. (In Corollaries 1 and 2, phi is
-                         used correctly again.) */
+                        
+                        // Note that Propositions 2-4 in Fox--Glynn mix up notation: they write Phi where they mean
+                        // 1 - Phi. (In Corollaries 1 and 2, phi is used correctly again.)
                         max_err = bl * exp(-0.5 * (k*k)) / k;
-                        if ( max_err * 2 <= epsilon ) {
-                            /* If the error on the left hand side is smaller, we can be
-                             more lenient on the right hand side. To this end, we now set
-                             epsilon to the part of the error that has not yet been eaten
-                             up by the left-hand truncation. */
+                        if (max_err * 2 <= epsilon) {
+                            // If the error on the left hand side is smaller, we can be more lenient on the right hand
+                            // side. To this end, we now set epsilon to the part of the error that has not yet been eaten
+                            // up by the left-hand truncation.
                             epsilon -= max_err;
                             break;
                         }
                     }
-                    /*Finally the left truncation point is found*/
+                    
+                    // Finally the left truncation point is found.
                 }
                 
-                /****Compute pFG->right truncation point****/
+                // Compute right truncation point.
                 {
-                    double lambda_max;
-                    int m_max, k;
-                    /*According to Fox-Glynn, if lambda < 400 we should take lambda = 400,
-                     otherwise use the original value. This is for computing the right truncation point*/
-                    if ( m < 400 ) {
-                        lambda_max = lambda_400;
+                    ValueType lambda_max;
+                    int64_t m_max, k;
+                    
+                    // According to Fox-Glynn, if lambda < 400 we should take lambda = 400, otherwise use the original
+                    // value. This is for computing the right truncation point.
+                    if (m < 400) {
+                        lambda_max = 400;
                         m_max = 400;
                         epsilon *= 0.662608824988162441697980;
                         /* i.e. al = (1+1/400) * exp(1/16) * sqrt_2; epsilon /= al; */
@@ -97,184 +109,157 @@ namespace storm {
                         epsilon *= (1 - 1 / (lambda + 1)) * 0.664265347050632847802225;
                         /* i.e. al = (1+1/lambda) * exp(1/16) * sqrt_2; epsilon /= al; */
                     }
-                    /* find right truncation point */
                     
-                    /* This loop is a modification to the original Fox-Glynn paper.
-                     The search for the right truncation point is only terminated by
-                     the error condition and not by the stop index from the FG paper.
-                     This can yield more accurate results if necessary.*/
-                    for ( k = 4 ; TRUE ; k++ )
-                    {
-                        /* dkl_inv is between 1 - 1e-33 and 1 if lambda_max >= 400 and
-                         k >= 4; this will always be rounded to 1.0. We therefore leave the
-                         factor out.
-                         double dkl_inv=1 - exp(-266/401.0 * (k*sqrt(2*lambda_max) + 1.5));
-                         */
-                        if ( k * epsilon /* actually: "k * (dkl_inv*epsilon/al)", which
-                                          has been calculated above */  >= exp(-0.5*(k*k)) )
+                    // Find right truncation point.
+                    
+                    // This loop is a modification to the original Fox-Glynn paper.
+                    // The search for the right truncation point is only terminated by  the error condition and not by
+                    // the stop index from the FG paper. This can yield more accurate results if necessary.
+                    for (k = 4;; ++k) {
+                        // dkl_inv is between 1 - 1e-33 and 1 if lambda_max >= 400 and k >= 4; this will always be
+                        // rounded to 1.0. We therefore leave the factor out.
+                        // double dkl_inv=1 - exp(-266/401.0 * (k*sqrt(2*lambda_max) + 1.5));
+                        
+                        // actually: "k * (dkl_inv*epsilon/al) >= exp(-0.5 * k^2)", but epsilon has been changed appropriately.
+                        if (k * epsilon >= exp(-0.5*(k*k))) {
                             break;
+                        }
                     }
-                    right = m_max + (int) ceil(k * sqrt(2 * lambda_max) + 0.5);
-                    if ( right > m_max + (int) ceil((lambda_max + 1) * 0.5) ) {
-                        printf("ERROR: Fox-Glynn: right = %d >> lambda = %g, cannot "
-                               "bound the right tail. The results are "
-                               "UNRELIABLE.\n", right, lambda_max);
+                    right = m_max + static_cast<int64_t>(std::ceil(k * std::sqrt(2 * lambda_max) + 0.5));
+                    if (right > m_max + static_cast<int64_t>(std::ceil((lambda_max + 1) * 0.5))) {
+                        STORM_LOG_WARN("Fox-Glynn: right = " << right << " >> lambda = " << lambda_max << ", cannot bound the right tail. The results are unreliable.");
                     }
                 }
                 
-                /*Time to set the initial value for weights*/
-                pFG = calloc((size_t) 1, sizeof(FoxGlynn) +
-                             (right - left) * sizeof(pFG->weights[0]));
-                if ( NULL == pFG ) {
-                    err_msg_3(err_MEMORY, "finder(%d,%g,_,%g,_)", m, lambda, omega, NULL);
-                }
-                pFG->right = right;
-                pFG->left = left;
-                pFG->weights[m - left] = omega / ( 1.0e+10 * (right - left) );
+                // Time to set the initial value for weights.
+                FoxGlynnResult<ValueType> fgresult;
+                fgresult.left = static_cast<uint64_t>(left);
+                fgresult.right = static_cast<uint64_t>(right);
+                fgresult.weights.resize(fgresult.right - fgresult.left + 1);
+
+                fgresult.weights[m - left] = omega / (1.0e+10 * (right - left));
                 
-                if ( m >= 25 )
-                {
-                    /* perform underflow check */
-                    double result, log_c_m_inf;
-                    int i;
+                if (m >= 25) {
+                    // Perform underflow check.
+                    ValueType result, log_c_m_inf;
+                    int64_t i;
                     
-                    /* we are going to compare with tau - log(w[m]) */
-                    tau -= log(pFG->weights[m - left]);
-                    /*We take the c_m_inf = 0.14627 / sqrt( m ), as for lambda >= 25
-                     c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf*/
-                    /* Note that m-lambda is in the interval (-1,0],
-                     and -1/(12*m) is in [-1/(12*25),0).
-                     So, exp(m-lambda - 1/(12*m)) is in (exp(-1-1/(12*25)),exp(0)).
-                     Therefore, we can improve the lower bound on c_m to
-                     exp(-1-1/(12*25)) / sqrt(2*pi) = ~0.14627. Its logarithm is
-                     -1 - 1/(12*25) - log(2*pi) * 0.5 = ~ -1.922272 (rounded towards
-                     -infinity). */
+                    // we are going to compare with tau - log(w[m]).
+                    tau -= std::log(fgresult.weights[m - left]);
+                    
+                    // We take the c_m_inf = 0.14627 / sqrt( m ), as for lambda >= 25
+                    // c_m = 1 / ( sqrt( 2.0 * pi * m ) ) * exp( m - lambda - 1 / ( 12.0 * m ) ) => c_m_inf.
+                    // Note that m-lambda is in the interval (-1,0], and -1/(12*m) is in [-1/(12*25),0).
+                    // So, exp(m-lambda - 1/(12*m)) is in (exp(-1-1/(12*25)),exp(0)).
+                    // Therefore, we can improve the lower bound on c_m to exp(-1-1/(12*25)) / sqrt(2*pi) = ~0.14627.
+                    // Its logarithm is -1 - 1/(12*25) - log(2*pi) * 0.5 = ~ -1.922272 (rounded towards -infinity).
                     log_c_m_inf = -1.922272 - log((double) m) * 0.5;
                     
-                    /* We use FG's Proposition 6 directly (and not Corollary 4 i and ii),
-                     as k_prime may be too large if pFG->left == 0. */
+                    // We use FG's Proposition 6 directly (and not Corollary 4 i and ii), as k_prime may be too large
+                    // if pFG->left == 0.
                     i = m - left;
-                    if ( i <= left /* equivalent to 2*i <= m,
-                                    equivalent to i <= lambda/2 */ )
-                    {
-                        /* Use Proposition 6 (i). Note that Fox--Glynn are off by one in
-                         the proof of this proposition; they sum up to i-1, but should have
-                         summed up to i. */
+                    
+                    // Equivalent to 2*i <= m, equivalent to i <= lambda/2.
+                    if (i <= left) {
+                        // Use Proposition 6 (i). Note that Fox--Glynn are off by one in the proof of this proposition;
+                        // they sum up to i-1, but should have summed up to i. */
                         result = log_c_m_inf
                         - i * (i+1) * (0.5 + (2*i+1)/(6*lambda)) / lambda;
-                    }
-                    else
-                    {
-                        /* Use Corollary 4 (iii). Note that k_prime <= sqrt(m+1)/m is a
-                         misprint for k_prime <= m/sqrt(m+1), which is equivalent to
-                         left >= 0, which holds trivially. */
+                    } else {
+                        // Use Corollary 4 (iii). Note that k_prime <= sqrt(m+1)/m is a misprint for k_prime <= m/sqrt(m+1),
+                        // which is equivalent to left >= 0, which holds trivially.
                         result = -lambda;
-                        if ( 0 != left ) {
-                            /* also use Proposition 6 (ii) */
+                        if (left != 0) {
+                            // Also use Proposition 6 (ii).
                             double result_1 = log_c_m_inf + i * log(1 - i/(double) (m+1));
-                            /*Take the maximum*/
-                            if ( result_1 > result )
+                            
+                            // Take the maximum.
+                            if (result_1 > result) {
                                 result = result_1;
+                            }
                         }
                     }
-                    if ( result <= tau )
-                    {
-                        const int log10_result = (int) floor(result * log10_e);
-                        printf("ERROR: Fox-Glynn: lambda >= 25, underflow near Poi(%g,%d)"
-                               " <= %.2fe%+d. The results are UNRELIABLE.\n",
-                               lambda, left, exp(result - log10_result/log10_e),
-                               log10_result);
+                    if (result <= tau) {
+                        int64_t const log10_result = static_cast<int64_t>(std::floor(result * log10_e));
+                        STORM_LOG_WARN("Fox-Glynn: lambda >= 25, underflow near Poi(" << lambda << "," << left << ") <= " << std::exp(result - log10_result/log10_e) << log10_result << ". The results are unreliable.");
                     }
                     
-                    /*We still have to perform an underflow check for the right truncation point when lambda >= 400*/
-                    if ( m >= 400 )
-                    {
-                        /* Use Proposition 5 of Fox--Glynn */
+                    // We still have to perform an underflow check for the right truncation point when lambda >= 400.
+                    if (m >= 400) {
+                        // Use Proposition 5 of Fox--Glynn.
                         i = right - m;
                         result = log_c_m_inf - i * (i + 1) / (2 * lambda);
-                        if ( result <= tau )
-                        {
-                            const int log10_result = (int) floor(result * log10_e);
-                            printf("ERROR: Fox-Glynn: lambda >= 400, underflow near "
-                                   "Poi(%g,%d) <= %.2fe%+d. The results are "
-                                   "UNRELIABLE.\n", lambda, right,
-                                   exp(result - log10_result / log10_e),
-                                   log10_result);
+                        if (result <= tau) {
+                            int64_t const log10_result = static_cast<int64_t>(std::floor(result * log10_e));
+                            STORM_LOG_WARN("Fox-Glynn: lambda >= 25, underflow near Poi(" << lambda << "," << right << ") <= " << std::exp(result - log10_result/log10_e) << log10_result << ". The results are unreliable.");
                         }
                     }
                 }
                 
-                return pFG;
+                return fgresult;
             }
             
             template<typename ValueType>
             FoxGlynnResult<ValueType> foxGlynnWeighter(ValueType lambda, ValueType epsilon) {
-                /*The magic m point*/
-                const uint64_t m = (int) floor(lambda);
-                int j, t;
-                FoxGlynn * pFG;
+                ValueType tau = std::numeric_limits<ValueType>::min();
+
+                // The magic m point.
+                uint64_t m = static_cast<uint64_t>(lambda);
+                int64_t j, t;
+
+                FoxGlynnResult<ValueType> result = foxGlynnFinder(lambda, epsilon);
                 
-                pFG = finder(m, lambda, tau, omega, epsilon);
-                if ( NULL == pFG ) {
-                    err_msg_4(err_CALLBY, "weighter(%g,%g,%g,%g)", lambda,
-                              tau, omega, epsilon, NULL);
+                // Fill the left side of the array.
+                for (j = m - result.left; j > 0; --j) {
+                    result.weights[j - 1] = (j + result.left) / lambda * result.weights[j];
                 }
                 
-                /*Fill the left side of the array*/
-                for ( j = m - pFG->left ; j > 0 ; j-- )
-                    pFG->weights[j-1] = (j+pFG->left) / lambda * pFG->weights[j];
+                t = result.right - result.left;
                 
-                t = pFG->right - pFG->left;
-                /*Fill the right side of the array, have two cases lambda < 400 & lambda >= 400*/
-                if ( m < 400 )
-                {
-                    /*Perform the underflow check, according to Fox-Glynn*/
-                    if( pFG->right > 600 )
-                    {
-                        printf("ERROR: Fox-Glynn: pFG->right > 600, underflow is possible\n");
-                        freeFG(pFG);
-                        return NULL;
-                    }
-                    /*Compute weights*/
-                    for ( j = m - pFG->left ; j < t ; j++ )
-                    {
-                        double q = lambda / (j + 1 + pFG->left);
-                        if ( pFG->weights[j] > tau / q )
-                        {
-                            pFG->weights[j + 1] = q * pFG->weights[j];
-                        }else{
+                // Fill the right side of the array, have two cases lambda < 400 & lambda >= 400.
+                if (m < 400) {
+                    // Perform the underflow check, according to Fox-Glynn.
+                    STORM_LOG_THROW(result.right <= 600, storm::exceptions::PrecisionExceededException, "Fox-Glynn: " << result.right << " > 600, underflow is possible.");
+
+                    // Compute weights.
+                    for (j = m - result.left; j < t; ++j) {
+                        ValueType q = lambda / (j + 1 + result.left);
+                        if (result.weights[j] > tau / q) {
+                            result.weights[j + 1] = q * result.weights[j];
+                        } else {
                             t = j;
-                            pFG->right = j + pFG->left;
-                            break; /*It's time to compute W*/
+                            result.right = j + result.left;
+                            
+                            // It's time to compute W.
+                            break;
                         }
                     }
-                }else{
-                    /*Compute weights*/
-                    for ( j = m - pFG->left ; j < t ; j++ )
-                        pFG->weights[j+1] = lambda / (j+1 + pFG->left) * pFG->weights[j];
+                } else {
+                    // Compute weights.
+                    for (j = m - result.left; j < t; ++j) {
+                        result.weights[j + 1] = lambda / (j + 1 + result.left) * result.weights[j];
+                    }
                 }
                 
-                /*It is time to compute the normalization weight W*/
-                pFG->total_weight = 0.0;
+                // It is time to compute the normalization weight W.
+                result.totalWeight = storm::utility::zero<ValueType>();
                 j = 0;
-                /* t was set above */
                 
-                while( j < t )
-                {
-                    if ( pFG->weights[j] <= pFG->weights[t] )
-                    {
-                        pFG->total_weight += pFG->weights[j];
+                // t was set above.
+                while(j < t) {
+                    if (result.weights[j] <= result.weights[t]) {
+                        result.totalWeight += result.weights[j];
                         j++;
-                    }else{
-                        pFG->total_weight += pFG->weights[t];
+                    } else {
+                        result.totalWeight += result.weights[t];
                         t--;
                     }
                 }
-                pFG->total_weight += pFG->weights[j];
+                result.totalWeight += result.weights[j];
                 
-                /* printf("Fox-Glynn: ltp = %d, rtp = %d, w = %10.15le \n", pFG->left, pFG->right, pFG->total_weight); */
+                STORM_LOG_TRACE("Fox-Glynn: ltp = " << result.left << ", rtp = " << result.right << ", w = " << result.totalWeight << ".");
                 
-                return pFG;
+                return result;
             }
             
             template<typename ValueType>
diff --git a/src/storm/utility/numerical.h b/src/storm/utility/numerical.h
index 233bb8a74..7bd18cebe 100644
--- a/src/storm/utility/numerical.h
+++ b/src/storm/utility/numerical.h
@@ -13,7 +13,7 @@ namespace storm {
                 uint64_t right;
                 ValueType totalWeight;
                 std::vector<ValueType> weights;
-            }
+            };
             
             template<typename ValueType>
             FoxGlynnResult<ValueType> foxGlynn(ValueType lambda, ValueType epsilon);

From cd34e3d67eb7a03fccd9adb02200cecb4222103a Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Sat, 2 Dec 2017 13:08:14 +0100
Subject: [PATCH 05/19] fixed issue in rational search preventing convergence
 in many cases

---
 src/storm/solver/NativeLinearEquationSolver.cpp        |  4 ++--
 .../solver/SymbolicMinMaxLinearEquationSolver.cpp      |  8 +++++---
 .../solver/SymbolicNativeLinearEquationSolver.cpp      | 10 ++++++----
 src/storm/storage/dd/Add.cpp                           |  6 +++++-
 src/storm/storage/dd/cudd/InternalCuddAdd.cpp          |  7 +++++++
 src/storm/storage/dd/cudd/InternalCuddAdd.h            |  5 +++++
 src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp      |  7 +++++++
 src/storm/storage/dd/sylvan/InternalSylvanAdd.h        |  2 ++
 8 files changed, 39 insertions(+), 10 deletions(-)

diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp
index c9974a700..d6de61266 100644
--- a/src/storm/solver/NativeLinearEquationSolver.cpp
+++ b/src/storm/solver/NativeLinearEquationSolver.cpp
@@ -776,9 +776,9 @@ namespace storm {
         template<typename ValueType>
         template<typename RationalType, typename ImpreciseType>
         bool NativeLinearEquationSolver<ValueType>::sharpen(uint64_t precision, storm::storage::SparseMatrix<RationalType> const& A, std::vector<ImpreciseType> const& x, std::vector<RationalType> const& b, std::vector<RationalType>& tmp) {
-            for (uint64_t p = 0; p <= precision; ++p) {
+            for (uint64_t p = 1; p <= precision; ++p) {
                 storm::utility::kwek_mehlhorn::sharpen(p, x, tmp);
-                
+
                 if (NativeLinearEquationSolver<RationalType>::isSolution(A, tmp, b)) {
                     return true;
                 }
diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
index b92ed6d9d..b1501fd32 100644
--- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
@@ -142,9 +142,10 @@ namespace storm {
         template<typename RationalType, typename ImpreciseType>
         storm::dd::Add<DdType, RationalType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
             
-            // Storage for the rational sharpened vector.
+            // Storage for the rational sharpened vector and the power iteration intermediate vector.
             storm::dd::Add<DdType, RationalType> sharpenedX;
-            
+            storm::dd::Add<DdType, ImpreciseType> currentX = x;
+
             // The actual rational search.
             uint64_t overallIterations = 0;
             uint64_t valueIterationInvocations = 0;
@@ -153,7 +154,7 @@ namespace storm {
             bool relative = env.solver().minMax().getRelativeTerminationCriterion();
             SolverStatus status = SolverStatus::InProgress;
             while (status == SolverStatus::InProgress && overallIterations < maxIter) {
-                typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter);
+                typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, currentX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter);
                 
                 ++valueIterationInvocations;
                 STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << viResult.iterations << " iterations.");
@@ -170,6 +171,7 @@ namespace storm {
                 if (isSolution) {
                     status = SolverStatus::Converged;
                 } else {
+                    currentX = viResult.values;
                     precision /= storm::utility::convertNumber<ValueType, uint64_t>(10);
                 }
             }
diff --git a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp b/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
index 760683bf5..4252b1483 100644
--- a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
@@ -174,10 +174,10 @@ namespace storm {
             
             storm::dd::Add<DdType, RationalType> sharpenedX;
             
-            for (uint64_t p = 1; p < precision; ++p) {
+            for (uint64_t p = 1; p <= precision; ++p) {
                 sharpenedX = x.sharpenKwekMehlhorn(p);
+
                 isSolution = rationalSolver.isSolutionFixedPoint(sharpenedX, rationalB);
-                
                 if (isSolution) {
                     break;
                 }
@@ -190,8 +190,9 @@ namespace storm {
         template<typename RationalType, typename ImpreciseType>
         storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
             
-            // Storage for the rational sharpened vector.
+            // Storage for the rational sharpened vector and the power iteration intermediate vector.
             storm::dd::Add<DdType, RationalType> sharpenedX;
+            storm::dd::Add<DdType, ImpreciseType> currentX = x;
             
             // The actual rational search.
             uint64_t overallIterations = 0;
@@ -201,7 +202,7 @@ namespace storm {
             bool relative = env.solver().native().getRelativeTerminationCriterion();
             SolverStatus status = SolverStatus::InProgress;
             while (status == SolverStatus::InProgress && overallIterations < maxIter) {
-                typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter - overallIterations);
+                typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(currentX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter - overallIterations);
                 
                 ++powerIterationInvocations;
                 STORM_LOG_TRACE("Completed " << powerIterationInvocations << " power iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations.");
@@ -218,6 +219,7 @@ namespace storm {
                 if (isSolution) {
                     status = SolverStatus::Converged;
                 } else {
+                    currentX = result.values;
                     precision /= storm::utility::convertNumber<ValueType, uint64_t>(10);
                 }
             }
diff --git a/src/storm/storage/dd/Add.cpp b/src/storm/storage/dd/Add.cpp
index a62600562..635d16aed 100644
--- a/src/storm/storage/dd/Add.cpp
+++ b/src/storm/storage/dd/Add.cpp
@@ -941,7 +941,7 @@ namespace storm {
         
         template<DdType LibraryType, typename ValueType>
         std::ostream& operator<<(std::ostream& out, Add<LibraryType, ValueType> const& add) {
-            out << "ADD with " << add.getNonZeroCount() << " nnz, " << add.getNodeCount() << " nodes, " << add.getLeafCount() << " leaves" << std::endl;
+            out << "ADD [" << add.getInternalAdd().getStringId() << "] with " << add.getNonZeroCount() << " nnz, " << add.getNodeCount() << " nodes, " << add.getLeafCount() << " leaves" << std::endl;
             std::vector<std::string> variableNames;
             for (auto const& variable : add.getContainedMetaVariables()) {
                 variableNames.push_back(variable.getName());
@@ -993,15 +993,19 @@ namespace storm {
         }
 		
         template class Add<storm::dd::DdType::CUDD, double>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::CUDD, double> const& add);
         template class Add<storm::dd::DdType::CUDD, uint_fast64_t>;
 
         template class Add<storm::dd::DdType::Sylvan, double>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::Sylvan, double> const& add);
         template class Add<storm::dd::DdType::Sylvan, uint_fast64_t>;
 
 #ifdef STORM_HAVE_CARL
         template class Add<storm::dd::DdType::CUDD, storm::RationalNumber>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::CUDD, storm::RationalNumber> const& add);
 
         template class Add<storm::dd::DdType::Sylvan, storm::RationalNumber>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::Sylvan, storm::RationalNumber> const& add);
 		template class Add<storm::dd::DdType::Sylvan, storm::RationalFunction>;
 
         template Add<storm::dd::DdType::CUDD, storm::RationalNumber> Add<storm::dd::DdType::CUDD, storm::RationalNumber>::toValueType<storm::RationalNumber>() const;
diff --git a/src/storm/storage/dd/cudd/InternalCuddAdd.cpp b/src/storm/storage/dd/cudd/InternalCuddAdd.cpp
index 83ace48ae..9645ebdf3 100644
--- a/src/storm/storage/dd/cudd/InternalCuddAdd.cpp
+++ b/src/storm/storage/dd/cudd/InternalCuddAdd.cpp
@@ -420,6 +420,13 @@ namespace storm {
         DdNode* InternalAdd<DdType::CUDD, ValueType>::getCuddDdNode() const {
             return this->getCuddAdd().getNode();
         }
+        
+        template<typename ValueType>
+        std::string InternalAdd<DdType::CUDD, ValueType>::getStringId() const {
+            std::stringstream ss;
+            ss << this->getCuddDdNode();
+            return ss.str();
+        }
 
         template<typename ValueType>
         Odd InternalAdd<DdType::CUDD, ValueType>::createOdd(std::vector<uint_fast64_t> const& ddVariableIndices) const {
diff --git a/src/storm/storage/dd/cudd/InternalCuddAdd.h b/src/storm/storage/dd/cudd/InternalCuddAdd.h
index 5871e20dc..2262d5da0 100644
--- a/src/storm/storage/dd/cudd/InternalCuddAdd.h
+++ b/src/storm/storage/dd/cudd/InternalCuddAdd.h
@@ -629,6 +629,11 @@ namespace storm {
              */
             DdNode* getCuddDdNode() const;
             
+            /*!
+             * Retrieves a string representation of an ID for thid ADD.
+             */
+            std::string getStringId() const;
+            
         private:
             /*!
              * Performs a recursive step to perform the given function between the given DD-based vector and the given
diff --git a/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp b/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
index b3723a4c4..0173f6aac 100644
--- a/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
+++ b/src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp
@@ -1151,6 +1151,13 @@ namespace storm {
             return sylvanMtbdd;
         }
         
+        template<typename ValueType>
+        std::string InternalAdd<DdType::Sylvan, ValueType>::getStringId() const {
+            std::stringstream ss;
+            ss << this->getSylvanMtbdd().GetMTBDD();
+            return ss.str();
+        }
+        
         template class InternalAdd<DdType::Sylvan, double>;
         template class InternalAdd<DdType::Sylvan, uint_fast64_t>;
 
diff --git a/src/storm/storage/dd/sylvan/InternalSylvanAdd.h b/src/storm/storage/dd/sylvan/InternalSylvanAdd.h
index 79a29f187..dc65e240a 100644
--- a/src/storm/storage/dd/sylvan/InternalSylvanAdd.h
+++ b/src/storm/storage/dd/sylvan/InternalSylvanAdd.h
@@ -620,6 +620,8 @@ namespace storm {
              */
             static ValueType getValue(MTBDD const& node);
             
+            std::string getStringId() const;
+            
         private:
             /*!
              * Recursively builds the ODD from an ADD.

From e67c04d2d6b7cbdfa213248b59e76147e7e79712 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Sat, 2 Dec 2017 13:09:43 +0100
Subject: [PATCH 06/19] update changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index cd45d0c7b..064dea378 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,6 +20,7 @@ The releases of major and minor versions contain an overview of changes since th
 - Performance improvements for sparse model building
 - Performance improvements for conditional properties on MDPs
 - Automatically convert MA without probabilistic states into CTMC
+- Fixed implemention of Fox and Glynn' algorithm
 - storm-pars: support for welldefinedness constraints in mdps.
 - storm-dft: split DFT settings into IO settings and fault tree settings
 - storm-dft: removed obsolete explicit model builder for DFTs

From a8caaf83c2d77fd4fc3f93883b3e76a01995445e Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Tue, 5 Dec 2017 16:28:14 +0100
Subject: [PATCH 07/19] made passing carl to sylvan more robust

---
 resources/3rdparty/CMakeLists.txt        | 2 +-
 resources/3rdparty/sylvan/CMakeLists.txt | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/resources/3rdparty/CMakeLists.txt b/resources/3rdparty/CMakeLists.txt
index ddef97000..1fde47b26 100644
--- a/resources/3rdparty/CMakeLists.txt
+++ b/resources/3rdparty/CMakeLists.txt
@@ -418,7 +418,7 @@ ExternalProject_Add(
         DOWNLOAD_COMMAND ""
         PREFIX "sylvan"
         SOURCE_DIR ${STORM_3RDPARTY_SOURCE_DIR}/sylvan
-        CMAKE_ARGS -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DGMP_LOCATION=${GMP_LIB_LOCATION}  -DGMP_INCLUDE=${GMP_INCLUDE_DIR}  -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DSYLVAN_BUILD_DOCS=OFF -DSYLVAN_BUILD_EXAMPLES=OFF -DCMAKE_BUILD_TYPE=${SYLVAN_BUILD_TYPE} -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DUSE_CARL=ON -Dcarl_INCLUDE_DIR=${carl_INCLUDE_DIR} -DSYLVAN_PORTABLE=${STORM_PORTABLE} -Dcarl_LIBRARIES=${carl_LIBRARIES} -DBoost_INCLUDE_DIRS=${Boost_INCLUDE_DIRS} -DBUILD_SHARED_LIBS=OFF -DSYLVAN_BUILD_TESTS=OFF
+        CMAKE_ARGS -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DGMP_LOCATION=${GMP_LIB_LOCATION}  -DGMP_INCLUDE=${GMP_INCLUDE_DIR}  -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DSYLVAN_BUILD_DOCS=OFF -DSYLVAN_BUILD_EXAMPLES=OFF -DCMAKE_BUILD_TYPE=${SYLVAN_BUILD_TYPE} -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DUSE_CARL=ON -Dcarl_DIR=${carl_DIR} -DSYLVAN_PORTABLE=${STORM_PORTABLE} -Dcarl_LIBRARIES=${carl_LIBRARIES} -DBoost_INCLUDE_DIRS=${Boost_INCLUDE_DIRS} -DBUILD_SHARED_LIBS=OFF -DSYLVAN_BUILD_TESTS=OFF
         BINARY_DIR ${STORM_3RDPARTY_BINARY_DIR}/sylvan
         BUILD_IN_SOURCE 0
         INSTALL_COMMAND ""
diff --git a/resources/3rdparty/sylvan/CMakeLists.txt b/resources/3rdparty/sylvan/CMakeLists.txt
index 4104f9219..60d0186c7 100644
--- a/resources/3rdparty/sylvan/CMakeLists.txt
+++ b/resources/3rdparty/sylvan/CMakeLists.txt
@@ -61,8 +61,8 @@ include_directories("${Boost_INCLUDE_DIRS}")
 
 if(USE_CARL)
     add_definitions(-DSYLVAN_HAVE_CARL)
-    include_directories("${carl_INCLUDE_DIR}")
     message(STATUS "Sylvan - Using CArL.")
+    find_package(carl REQUIRED PATHS ${carl_DIR} NO_DEFAULT_PATH)
 else()
     message(STATUS "Sylvan - Not using CArL.")
 endif()

From f6eadc14ca28fecf8e52fd8e2ac044083486391f Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 09:43:12 +0100
Subject: [PATCH 08/19] adding proper carl_DIR when using shipped carl

---
 resources/3rdparty/CMakeLists.txt | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/resources/3rdparty/CMakeLists.txt b/resources/3rdparty/CMakeLists.txt
index 1fde47b26..b69832778 100644
--- a/resources/3rdparty/CMakeLists.txt
+++ b/resources/3rdparty/CMakeLists.txt
@@ -280,7 +280,8 @@ else()
 
 	add_dependencies(resources carl)
     set(carl_INCLUDE_DIR "${STORM_3RDPARTY_BINARY_DIR}/carl/include/")
-	set(carl_LIBRARIES ${STORM_3RDPARTY_BINARY_DIR}/carl/lib/libcarl${DYNAMIC_EXT})
+    set(carl_DIR "${STORM_3RDPARTY_BINARY_DIR}/carl/")
+    set(carl_LIBRARIES ${STORM_3RDPARTY_BINARY_DIR}/carl/lib/libcarl${DYNAMIC_EXT})
     set(STORM_HAVE_CARL ON)
 
     message(STATUS "Storm - Linking with shipped carl ${carl_VERSION} (include: ${carl_INCLUDE_DIR}, library ${carl_LIBRARIES}, CARL_USE_CLN_NUMBERS: ${CARL_USE_CLN_NUMBERS}, CARL_USE_GINAC: ${CARL_USE_GINAC}).")
@@ -418,7 +419,7 @@ ExternalProject_Add(
         DOWNLOAD_COMMAND ""
         PREFIX "sylvan"
         SOURCE_DIR ${STORM_3RDPARTY_SOURCE_DIR}/sylvan
-        CMAKE_ARGS -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DGMP_LOCATION=${GMP_LIB_LOCATION}  -DGMP_INCLUDE=${GMP_INCLUDE_DIR}  -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DSYLVAN_BUILD_DOCS=OFF -DSYLVAN_BUILD_EXAMPLES=OFF -DCMAKE_BUILD_TYPE=${SYLVAN_BUILD_TYPE} -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DUSE_CARL=ON -Dcarl_DIR=${carl_DIR} -DSYLVAN_PORTABLE=${STORM_PORTABLE} -Dcarl_LIBRARIES=${carl_LIBRARIES} -DBoost_INCLUDE_DIRS=${Boost_INCLUDE_DIRS} -DBUILD_SHARED_LIBS=OFF -DSYLVAN_BUILD_TESTS=OFF
+        CMAKE_ARGS -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DGMP_LOCATION=${GMP_LIB_LOCATION}  -DGMP_INCLUDE=${GMP_INCLUDE_DIR}  -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DSYLVAN_BUILD_DOCS=OFF -DSYLVAN_BUILD_EXAMPLES=OFF -DCMAKE_BUILD_TYPE=${SYLVAN_BUILD_TYPE} -DCMAKE_POSITION_INDEPENDENT_CODE=ON -DUSE_CARL=ON -Dcarl_DIR=${carl_DIR} -DSYLVAN_PORTABLE=${STORM_PORTABLE} -DBoost_INCLUDE_DIRS=${Boost_INCLUDE_DIRS} -DBUILD_SHARED_LIBS=OFF -DSYLVAN_BUILD_TESTS=OFF
         BINARY_DIR ${STORM_3RDPARTY_BINARY_DIR}/sylvan
         BUILD_IN_SOURCE 0
         INSTALL_COMMAND ""

From c1102209e83fdcbc9d571e0cb021a26caa883512 Mon Sep 17 00:00:00 2001
From: Sebastian Junges <sebastian.junges@rwth-aachen.de>
Date: Wed, 6 Dec 2017 11:32:57 +0100
Subject: [PATCH 09/19] update changelog

---
 CHANGELOG.md | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 064dea378..9976e61cb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,9 +4,10 @@ Changelog
 This changelog lists only the most important changes. Smaller (bug)fixes as well as non-mature features are not part of the changelog.
 The releases of major and minor versions contain an overview of changes since the last major/minor update.
 
+Version 1.2.x
+-------------
 
-
-### Version 1.2
+### Version 1.2.0
 - C++ api changes: Building model takes BuilderOptions instead of extended list of Booleans, does not depend on settings anymore.
 - storm-cli-utilities now contains cli related stuff, instead of storm-lib
 - Symbolic (MT/BDD) bisimulation

From 945e360bcfdde271cbcf25b619dcc11e4f7f8937 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 14:39:22 +0100
Subject: [PATCH 10/19] Increased minimal carl version to 17.10

---
 resources/3rdparty/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/resources/3rdparty/CMakeLists.txt b/resources/3rdparty/CMakeLists.txt
index b69832778..3cfefe1cc 100644
--- a/resources/3rdparty/CMakeLists.txt
+++ b/resources/3rdparty/CMakeLists.txt
@@ -208,7 +208,7 @@ include(${STORM_3RDPARTY_SOURCE_DIR}/include_cudd.cmake)
 #############################################################
 
 set(STORM_HAVE_CARL OFF)
-set(CARL_MINVERSION "17.08")
+set(CARL_MINVERSION "17.10")
 if (NOT STORM_FORCE_SHIPPED_CARL)
     if (NOT "${STORM_CARL_DIR_HINT}" STREQUAL "")
 		find_package(carl QUIET PATHS ${STORM_CARL_DIR_HINT} NO_DEFAULT_PATH)

From dbcc49f48f3e57750484945c4431c5cc8a4517aa Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 14:39:40 +0100
Subject: [PATCH 11/19] Use new carl release 17.12

---
 resources/3rdparty/carl/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/resources/3rdparty/carl/CMakeLists.txt b/resources/3rdparty/carl/CMakeLists.txt
index dcab4ec30..4b819682d 100644
--- a/resources/3rdparty/carl/CMakeLists.txt
+++ b/resources/3rdparty/carl/CMakeLists.txt
@@ -8,7 +8,7 @@ message(STORM_3RDPARTY_BINARY_DIR: ${STORM_3RDPARTY_BINARY_DIR})
 
 ExternalProject_Add(carl-config
 	GIT_REPOSITORY https://github.com/smtrat/carl
-	GIT_TAG 17.10
+	GIT_TAG 17.12
 	PREFIX here
 	SOURCE_DIR source_dir
 	BINARY_DIR ${STORM_3RDPARTY_BINARY_DIR}/carl

From 0481ca38551b26981b02908ec031a589099405da Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 15:36:08 +0100
Subject: [PATCH 12/19] Fixed deprecated getType()

---
 src/storm/adapters/Smt2ExpressionAdapter.h | 2 +-
 src/storm/solver/SmtlibSmtSolver.cpp       | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/storm/adapters/Smt2ExpressionAdapter.h b/src/storm/adapters/Smt2ExpressionAdapter.h
index f1b96b71e..5abaf290d 100644
--- a/src/storm/adapters/Smt2ExpressionAdapter.h
+++ b/src/storm/adapters/Smt2ExpressionAdapter.h
@@ -133,7 +133,7 @@ namespace storm {
                         STORM_LOG_DEBUG("Declaring the variable " + variableString);
                         declaredVariables.back().insert(variableString);
                         std::string varDeclaration = "( declare-fun " + variableString + " () ";
-                        switch (variableToCheck.getType()){
+                        switch (variableToCheck.type()){
                             case carl::VariableType::VT_BOOL:
                                 varDeclaration += "Bool";
                                 break;
diff --git a/src/storm/solver/SmtlibSmtSolver.cpp b/src/storm/solver/SmtlibSmtSolver.cpp
index 69394ec8a..aad46d81e 100644
--- a/src/storm/solver/SmtlibSmtSolver.cpp
+++ b/src/storm/solver/SmtlibSmtSolver.cpp
@@ -103,7 +103,7 @@ namespace storm {
         }
         
         void SmtlibSmtSolver::add(const storm::RationalFunctionVariable& variable, bool value){
-            STORM_LOG_THROW((variable.getType()==carl::VariableType::VT_BOOL), storm::exceptions::IllegalArgumentException, "Tried to add a constraint that consists of a non-boolean variable.");
+            STORM_LOG_THROW((variable.type()==carl::VariableType::VT_BOOL), storm::exceptions::IllegalArgumentException, "Tried to add a constraint that consists of a non-boolean variable.");
             std::set<storm::RationalFunctionVariable> variableSet;
             variableSet.insert(variable);
             std::vector<std::string> const varDeclarations = expressionAdapter->checkForUndeclaredVariables(variableSet);

From 24c0e23db2608f35d29056087545af8842b0afbe Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 16:25:25 +0100
Subject: [PATCH 13/19] Added checklist for new release

---
 doc/checklist_new_release.md | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)
 create mode 100644 doc/checklist_new_release.md

diff --git a/doc/checklist_new_release.md b/doc/checklist_new_release.md
new file mode 100644
index 000000000..c279ffa95
--- /dev/null
+++ b/doc/checklist_new_release.md
@@ -0,0 +1,30 @@
+The following steps should be performed before releasing a new storm version.
+Note that in most case a simultaneous release of [carl](https://github.com/smtrat/carl), [storm](https://github.com/moves-rwth/storm), [pycarl](https://github.com/moves-rwth/pycarl/) and [stormpy](https://github.com/moves-rwth/stormpy/) is preferred.
+
+1. Update `CHANGELOG.md`
+   To get all the commits from an author since the last tag execute:
+   ```console
+   git log last_tag..HEAD --author "author_name"
+   ```
+
+2. Update used carl version:
+  * Update `GIT_TAG` in `resources/3rdparty/carl/CMakeLists.txt`
+  * Maybe update `CARL_MINVERSION` in `resources/3rdparty/CMakeLists.txt`
+
+3. Check that storm builds without errors and all tests are successful
+   * [Travis](https://travis-ci.org/moves-rwth/storm) should run successfully
+
+4. Set new storm version:
+   * Set new storm version in `version.cmake`
+
+5. Set new tag in git
+   ```console
+   git tag new_version
+   git push origin new_version
+   ```
+
+6. [Add new release](https://github.com/moves-rwth/storm/releases/new) in Github
+
+7. Update [Homebrew formula](https://github.com/moves-rwth/homebrew-storm)
+
+8. Announce new storm version on [website](http://www.stormchecker.org/news.html)

From 90a42d775fea3f55dda6dd45ef59145806f40f88 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 16:32:55 +0100
Subject: [PATCH 14/19] Use markdown in CHANGELOG

---
 CHANGELOG.md | 20 +++++++++-----------
 1 file changed, 9 insertions(+), 11 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9976e61cb..8218b4c83 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,9 +7,9 @@ The releases of major and minor versions contain an overview of changes since th
 Version 1.2.x
 -------------
 
-### Version 1.2.0
-- C++ api changes: Building model takes BuilderOptions instead of extended list of Booleans, does not depend on settings anymore.
-- storm-cli-utilities now contains cli related stuff, instead of storm-lib
+### Version 1.2.0 (2017/12)
+- C++ api changes: Building model takes `BuilderOptions` instead of extended list of Booleans, does not depend on settings anymore.
+- `storm-cli-utilities` now contains cli related stuff, instead of `storm-lib`
 - Symbolic (MT/BDD) bisimulation
 - Fixed issue related to variable names that can not be used in Exprtk.
 - DRN parser improved
@@ -22,9 +22,9 @@ Version 1.2.x
 - Performance improvements for conditional properties on MDPs
 - Automatically convert MA without probabilistic states into CTMC
 - Fixed implemention of Fox and Glynn' algorithm
-- storm-pars: support for welldefinedness constraints in mdps.
-- storm-dft: split DFT settings into IO settings and fault tree settings
-- storm-dft: removed obsolete explicit model builder for DFTs
+- `storm-pars`: support for welldefinedness constraints in mdps.
+- `storm-dft`: split DFT settings into IO settings and fault tree settings
+- `storm-dft`: removed obsolete explicit model builder for DFTs
 - Features for developers:
 	* Solvers can now expose requirements
 	* unbounded reachability and reachability rewards now correctly respect solver requirements
@@ -36,10 +36,9 @@ Version 1.1.x
 -------------
 
 ### Version 1.1.0 (2017/8)
-
 - Support for long-run average rewards on MDPs and Markov automata using a value-iteration based approach.
 - Storm can now check MDPs and Markov Automata (i.e. MinMax equation systems) via Linear Programming.
-- Parametric model checking is now handled in a separated library/executable called storm-pars.
+- Parametric model checking is now handled in a separated library/executable called `storm-pars`.
 - Wellformedness constraints on PMCs:
     * include constraints from rewards
     * are in smtlib2
@@ -50,15 +49,14 @@ Version 1.1.x
 - Support for parsing/building models given in the explicit input format of IMCA.
 - Storm now overwrites files if asked to write files to a specific location.
 - Changes in build process to accommodate for changes in carl. Also, more robust against issues with carl.
-- USE_POPCNT removed in favor of FORCE_POPCNT. The popcnt instruction is used if available due to march=native, unless portable is set.
-    Then, using FORCE_POPCNT enables the use of the SSE 4.2 instruction
+- `USE_POPCNT` removed in favor of `FORCE_POPCNT`. The popcnt instruction is used if available due to `march=native`, unless portable is set.
+  Then, using `FORCE_POPCNT` enables the use of the SSE 4.2 instruction
 
 
 Version 1.0.x
 -------------
 
 ### Version 1.0.1 (2017/4)
-
 - Multi-objective model checking support now fully included
 - Several improvements in parameter lifting
 - Several improvements in JANI parsing

From 240899d32dc6125766d8f5f10a83b735c582b998 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 16:49:05 +0100
Subject: [PATCH 15/19] Storm release 1.2.0

---
 version.cmake | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/version.cmake b/version.cmake
index 21ee2f621..d76de7858 100644
--- a/version.cmake
+++ b/version.cmake
@@ -1,4 +1,4 @@
 set(STORM_VERSION_MAJOR 1)
-set(STORM_VERSION_MINOR 1)
+set(STORM_VERSION_MINOR 2)
 set(STORM_VERSION_PATCH 0)
 

From 4e34754ac1ac137813ccf4ef000cda1f44a9d1d8 Mon Sep 17 00:00:00 2001
From: Matthias Volk <matthias.volk@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 17:29:13 +0100
Subject: [PATCH 16/19] Updated release steps for setting tag on GitHub

---
 doc/checklist_new_release.md | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/doc/checklist_new_release.md b/doc/checklist_new_release.md
index c279ffa95..a0462b5d8 100644
--- a/doc/checklist_new_release.md
+++ b/doc/checklist_new_release.md
@@ -19,11 +19,17 @@ Note that in most case a simultaneous release of [carl](https://github.com/smtra
 
 5. Set new tag in git
    ```console
-   git tag new_version
+   git tag -a new_version
    git push origin new_version
    ```
+   Next we push the tag to GitHub. This step requires the GitHub repo to to be configured as a remote.
+   ```console
+   git remote add github https://github.com/moves-rwth/storm.git
+   git push github new_version
+   ```
+   The new tag should now be visible on [GitHub](https://github.com/moves-rwth/storm/tags)
 
-6. [Add new release](https://github.com/moves-rwth/storm/releases/new) in Github
+6. [Add new release](https://github.com/moves-rwth/storm/releases/new) in GitHub
 
 7. Update [Homebrew formula](https://github.com/moves-rwth/homebrew-storm)
 

From df86b6c8154f9287775f01d3c4744b690c29bc25 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Wed, 6 Dec 2017 19:57:38 +0100
Subject: [PATCH 17/19] fixing issue related to relevant value restriction in
 conditional properties

---
 .../prctl/helper/SparseDtmcPrctlHelper.cpp    | 30 +++++++++++++++++++
 .../prctl/helper/SparseDtmcPrctlHelper.h      |  5 +++-
 2 files changed, 34 insertions(+), 1 deletion(-)

diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
index 0d0ec5928..781024096 100644
--- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
@@ -453,6 +453,22 @@ namespace storm {
                 
                 return result;
             }
+
+            template<typename ValueType, typename RewardModelType>
+            storm::storage::BitVector SparseDtmcPrctlHelper<ValueType, RewardModelType>::BaierTransformedModel::getNewRelevantStates() const {
+                storm::storage::BitVector newRelevantStates(transitionMatrix.get().getRowCount());
+                for (uint64_t i = 0; i < this->beforeStates.getNumberOfSetBits(); ++i) {
+                    newRelevantStates.set(i);
+                }
+                return newRelevantStates;
+            }
+            
+            template<typename ValueType, typename RewardModelType>
+            storm::storage::BitVector SparseDtmcPrctlHelper<ValueType, RewardModelType>::BaierTransformedModel::getNewRelevantStates(storm::storage::BitVector const& oldRelevantStates) const {
+                storm::storage::BitVector result = oldRelevantStates % this->beforeStates;
+                result.resize(transitionMatrix.get().getRowCount());
+                return result;
+            }
             
             template<typename ValueType, typename RewardModelType>
             std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeConditionalProbabilities(Environment const& env, storm::solver::SolveGoal<ValueType>&& goal, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
@@ -472,6 +488,13 @@ namespace storm {
                         
                         // Now compute reachability probabilities in the transformed model.
                         storm::storage::SparseMatrix<ValueType> const& newTransitionMatrix = transformedModel.transitionMatrix.get();
+                        storm::storage::BitVector newRelevantValues;
+                        if (goal.hasRelevantValues()) {
+                            newRelevantValues = transformedModel.getNewRelevantStates(goal.relevantValues());
+                        } else {
+                            newRelevantValues = transformedModel.getNewRelevantStates();
+                        }
+                        goal.setRelevantValues(std::move(newRelevantValues));
                         std::vector<ValueType> conditionalProbabilities = computeUntilProbabilities(env, std::move(goal), newTransitionMatrix, newTransitionMatrix.transpose(), storm::storage::BitVector(newTransitionMatrix.getRowCount(), true), transformedModel.targetStates.get(), qualitative, linearEquationSolverFactory);
                         
                         storm::utility::vector::setVectorValues(result, transformedModel.beforeStates, conditionalProbabilities);
@@ -498,6 +521,13 @@ namespace storm {
                         
                         // Now compute reachability probabilities in the transformed model.
                         storm::storage::SparseMatrix<ValueType> const& newTransitionMatrix = transformedModel.transitionMatrix.get();
+                        storm::storage::BitVector newRelevantValues;
+                        if (goal.hasRelevantValues()) {
+                            newRelevantValues = transformedModel.getNewRelevantStates(goal.relevantValues());
+                        } else {
+                            newRelevantValues = transformedModel.getNewRelevantStates();
+                        }
+                        goal.setRelevantValues(std::move(newRelevantValues));
                         std::vector<ValueType> conditionalRewards = computeReachabilityRewards(env, std::move(goal), newTransitionMatrix, newTransitionMatrix.transpose(), transformedModel.stateRewards.get(), transformedModel.targetStates.get(), qualitative, linearEquationSolverFactory);
                         storm::utility::vector::setVectorValues(result, transformedModel.beforeStates, conditionalRewards);
                     }
diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
index 6564c4f5d..0c0712723 100644
--- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
+++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
@@ -59,7 +59,10 @@ namespace storm {
                     BaierTransformedModel() : noTargetStates(false) {
                         // Intentionally left empty.
                     }
-                    
+
+                    storm::storage::BitVector getNewRelevantStates(storm::storage::BitVector const& oldRelevantStates) const;
+                    storm::storage::BitVector getNewRelevantStates() const;
+
                     storm::storage::BitVector beforeStates;
                     boost::optional<storm::storage::SparseMatrix<ValueType>> transitionMatrix;
                     boost::optional<storm::storage::BitVector> targetStates;

From 905ae821f36b7eda508188bbc9f926c297eed114 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Fri, 8 Dec 2017 11:47:16 +0100
Subject: [PATCH 18/19] extended SMT-based minimal label set generator so that
 it can deal with lower-bounded properties (however loosing the minimality
 property in some sense)

---
 .../SMTMinimalLabelSetGenerator.h             | 159 +++++++++++++-----
 src/storm/logic/ComparisonType.h              |  57 ++++---
 .../storage/sparse/PrismChoiceOrigins.cpp     |   2 +-
 src/storm/storage/sparse/PrismChoiceOrigins.h |   2 +-
 4 files changed, 155 insertions(+), 65 deletions(-)

diff --git a/src/storm/counterexamples/SMTMinimalLabelSetGenerator.h b/src/storm/counterexamples/SMTMinimalLabelSetGenerator.h
index 008a81e0f..e8f2464a9 100644
--- a/src/storm/counterexamples/SMTMinimalLabelSetGenerator.h
+++ b/src/storm/counterexamples/SMTMinimalLabelSetGenerator.h
@@ -17,6 +17,9 @@
 #include "storm/settings/SettingsManager.h"
 #include "storm/settings/modules/CoreSettings.h"
 
+#include "storm/utility/macros.h"
+#include "storm/exceptions/NotSupportedException.h"
+
 #include "storm/utility/counterexamples.h"
 #include "storm/utility/cli.h"
 
@@ -596,28 +599,6 @@ namespace storm {
                 // cuts.
                 std::unique_ptr<storm::solver::SmtSolver> localSolver(new storm::solver::Z3SmtSolver(program.getManager()));
                 storm::expressions::ExpressionManager const& localManager = program.getManager();
-//
-//                // Create a context and register all variables of the program with their correct type.
-//                z3::context localContext;
-//                z3::solver localSolver(localContext);
-//                std::map<std::string, z3::expr> solverVariables;
-//                for (auto const& booleanVariable : program.getGlobalBooleanVariables()) {
-//                    solverVariables.emplace(booleanVariable.getName(), localContext.bool_const(booleanVariable.getName().c_str()));
-//                }
-//                for (auto const& integerVariable : program.getGlobalIntegerVariables()) {
-//                    solverVariables.emplace(integerVariable.getName(), localContext.int_const(integerVariable.getName().c_str()));
-//                }
-//
-//                for (auto const& module : program.getModules()) {
-//                    for (auto const& booleanVariable : module.getBooleanVariables()) {
-//                        solverVariables.emplace(booleanVariable.getName(), localContext.bool_const(booleanVariable.getName().c_str()));
-//                    }
-//                    for (auto const& integerVariable : module.getIntegerVariables()) {
-//                        solverVariables.emplace(integerVariable.getName(), localContext.int_const(integerVariable.getName().c_str()));
-//                    }
-//                }
-//
-//                storm::adapters::Z3ExpressionAdapter expressionAdapter(localContext, false, solverVariables);
 
                 // Then add the constraints for bounds of the integer variables..
                 for (auto const& integerVariable : program.getGlobalIntegerVariables()) {
@@ -1592,9 +1573,7 @@ namespace storm {
              * Returns the submdp obtained from removing all choices that do not originate from the specified filterLabelSet.
              * Also returns the Labelsets of the submdp
              */
-            static std::pair<storm::models::sparse::Mdp<T>, std::vector<boost::container::flat_set<uint_fast64_t>>> restrictMdpToLabelSet(storm::models::sparse::Mdp<T> const& mdp,  std::vector<boost::container::flat_set<uint_fast64_t>> const& labelSets, boost::container::flat_set<uint_fast64_t> const& filterLabelSet) {
-
-                STORM_LOG_THROW(mdp.getNumberOfChoices() == labelSets.size(), storm::exceptions::InvalidArgumentException, "The given number of labels does not match the number of choices.");
+            static std::pair<storm::models::sparse::Mdp<T>, std::vector<boost::container::flat_set<uint_fast64_t>>> restrictMdpToLabelSet(storm::models::sparse::Mdp<T> const& mdp,  boost::container::flat_set<uint_fast64_t> const& filterLabelSet) {
 
                 std::vector<boost::container::flat_set<uint_fast64_t>> resultLabelSet;
                 storm::storage::SparseMatrixBuilder<T> transitionMatrixBuilder(0, mdp.getTransitionMatrix().getColumnCount(), 0, true, true, mdp.getTransitionMatrix().getRowGroupCount());
@@ -1604,7 +1583,8 @@ namespace storm {
                 for(uint_fast64_t state = 0; state < mdp.getNumberOfStates(); ++state) {
                     bool stateHasValidChoice = false;
                     for (uint_fast64_t choice = mdp.getTransitionMatrix().getRowGroupIndices()[state]; choice < mdp.getTransitionMatrix().getRowGroupIndices()[state + 1]; ++choice) {
-                        bool choiceValid = std::includes(filterLabelSet.begin(), filterLabelSet.end(), labelSets[choice].begin(), labelSets[choice].end());
+                        auto const& choiceLabelSet = mdp.getChoiceOrigins()->asPrismChoiceOrigins().getCommandSet(choice);
+                        bool choiceValid = std::includes(filterLabelSet.begin(), filterLabelSet.end(), choiceLabelSet.begin(), choiceLabelSet.end());
 
                         // If the choice is valid, copy over all its elements.
                         if (choiceValid) {
@@ -1615,7 +1595,7 @@ namespace storm {
                             for (auto const& entry : mdp.getTransitionMatrix().getRow(choice)) {
                                 transitionMatrixBuilder.addNextValue(currentRow, entry.getColumn(), entry.getValue());
                             }
-                            resultLabelSet.push_back(labelSets[choice]);
+                            resultLabelSet.push_back(choiceLabelSet);
                             ++currentRow;
                         }
                     }
@@ -1636,8 +1616,6 @@ namespace storm {
             }
 
         public:
-         
-            
             /*!
              * Computes the minimal command set that is needed in the given MDP to exceed the given probability threshold for satisfying phi until psi.
              *
@@ -1645,9 +1623,8 @@ namespace storm {
              * @param mdp The MDP in which to find the minimal command set.
              * @param phiStates A bit vector characterizing all phi states in the model.
              * @param psiStates A bit vector characterizing all psi states in the model.
-             * @param probabilityThreshold The probability value that must be achieved or exceeded.
-             * @param strictBound A flag indicating whether the probability must be achieved (in which case the flag must be set) or strictly exceeded
-             * (if the flag is set to false).
+             * @param probabilityThreshold The threshold that is to be achieved or exceeded.
+             * @param strictBound Indicates whether the threshold needs to be achieved (true) or exceeded (false).
              * @param checkThresholdFeasible If set, it is verified that the model can actually achieve/exceed the given probability value. If this check
              * is made and fails, an exception is thrown.
              */
@@ -1672,7 +1649,7 @@ namespace storm {
 
                 // (0) Obtain the label sets for each choice.
                 // The label set of a choice corresponds to the set of prism commands that induce the choice.
-                STORM_LOG_THROW(mdp.hasChoiceOrigins(), storm::exceptions::InvalidArgumentException, "Restriction to minimal command set is impossible for model without choice origns.");
+                STORM_LOG_THROW(mdp.hasChoiceOrigins(), storm::exceptions::InvalidArgumentException, "Restriction to minimal command set is impossible for model without choice origins.");
                 STORM_LOG_THROW(mdp.getChoiceOrigins()->isPrismChoiceOrigins(), storm::exceptions::InvalidArgumentException, "Restriction to command set is impossible for model without prism choice origins.");
                 storm::storage::sparse::PrismChoiceOrigins const& choiceOrigins = mdp.getChoiceOrigins()->asPrismChoiceOrigins();
                 std::vector<boost::container::flat_set<uint_fast64_t>> labelSets;
@@ -1698,8 +1675,8 @@ namespace storm {
                 RelevancyInformation relevancyInformation = determineRelevantStatesAndLabels(mdp, labelSets, phiStates, psiStates);
                 
                 // (3) Create a solver.
-                std::shared_ptr<storm::expressions::ExpressionManager> manager(new storm::expressions::ExpressionManager());
-                std::unique_ptr<storm::solver::SmtSolver> solver(new storm::solver::Z3SmtSolver(*manager));
+                std::shared_ptr<storm::expressions::ExpressionManager> manager = std::make_shared<storm::expressions::ExpressionManager>();
+                std::unique_ptr<storm::solver::SmtSolver> solver = std::make_unique<storm::solver::Z3SmtSolver>(*manager);
                 
                 // (4) Create the variables for the relevant commands.
                 VariableInformation variableInformation = createVariables(manager, mdp, psiStates, relevancyInformation, includeReachabilityEncoding);
@@ -1747,7 +1724,7 @@ namespace storm {
                     // Restrict the given MDP to the current set of labels and compute the reachability probability.
                     modelCheckingClock = std::chrono::high_resolution_clock::now();
                     commandSet.insert(relevancyInformation.knownLabels.begin(), relevancyInformation.knownLabels.end());
-                    auto subMdpChoiceOrigins = restrictMdpToLabelSet(mdp, labelSets, commandSet);
+                    auto subMdpChoiceOrigins = restrictMdpToLabelSet(mdp, commandSet);
                     storm::models::sparse::Mdp<T> const& subMdp = subMdpChoiceOrigins.first;
                     std::vector<boost::container::flat_set<uint_fast64_t>> const& subLabelSets = subMdpChoiceOrigins.second;
   
@@ -1811,25 +1788,94 @@ namespace storm {
 #endif
             }
             
-            static std::shared_ptr<PrismHighLevelCounterexample> computeCounterexample(Environment const& env,storm::prism::Program program, storm::models::sparse::Mdp<T> const& mdp, std::shared_ptr<storm::logic::Formula const> const& formula) {
+            static void extendCommandSetLowerBound(storm::models::sparse::Mdp<T> const& mdp, boost::container::flat_set<uint_fast64_t>& commandSet, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                auto startTime = std::chrono::high_resolution_clock::now();
+                
+                // Create sub-MDP that only contains the choices allowed by the given command set.
+                storm::models::sparse::Mdp<T> subMdp = restrictMdpToLabelSet(mdp, commandSet).first;
+                
+                // Then determine all prob0E(psi) states that are reachable in the sub-MDP.
+                storm::storage::BitVector reachableProb0EStates = storm::utility::graph::getReachableStates(subMdp.getTransitionMatrix(), subMdp.getInitialStates(), phiStates, psiStates);
+                
+                // Create a queue of reachable prob0E(psi) states so we can check which commands need to be added
+                // to give them a strategy that avoids psi states.
+                std::queue<uint_fast64_t> prob0EWorklist;
+                for (auto const& e : reachableProb0EStates) {
+                    prob0EWorklist.push(e);
+                }
+                
+                // As long as there are reachable prob0E(psi) states, we add commands so they can stay within
+                // prob0E(states).
+                while (!prob0EWorklist.empty()) {
+                    uint_fast64_t state = prob0EWorklist.front();
+                    prob0EWorklist.pop();
+                    
+                    // Now iterate over the original choices of the prob0E(psi) state and add at least one.
+                    bool hasLabeledChoice = false;
+                    uint64_t smallestCommandSetSize = 0;
+                    uint64_t smallestCommandChoice = mdp.getTransitionMatrix().getRowGroupIndices()[state];
+                    
+                    // Determine the choice with the least amount of commands (bad heuristic).
+                    for (uint64_t choice = smallestCommandChoice; choice < mdp.getTransitionMatrix().getRowGroupIndices()[state + 1]; ++choice) {
+                        bool onlyProb0ESuccessors = true;
+                        for (auto const& successorEntry : mdp.getTransitionMatrix().getRow(choice)) {
+                            if (!psiStates.get(successorEntry.getColumn())) {
+                                onlyProb0ESuccessors = false;
+                                break;
+                            }
+                        }
+                        
+                        if (onlyProb0ESuccessors) {
+                            auto const& labelSet = mdp.getChoiceOrigins()->asPrismChoiceOrigins().getCommandSet(choice);
+                            hasLabeledChoice |= !labelSet.empty();
+                            
+                            if (smallestCommandChoice == 0 || labelSet.size() < smallestCommandSetSize) {
+                                smallestCommandSetSize = labelSet.size();
+                                smallestCommandChoice = choice;
+                            }
+                        }
+                    }
+                    
+                    if (hasLabeledChoice) {
+                        // Take all labels of the selected choice.
+                        auto const& labelSet = mdp.getChoiceOrigins()->asPrismChoiceOrigins().getCommandSet(smallestCommandChoice);
+                        commandSet.insert(labelSet.begin(), labelSet.end());
+                        
+                        // Check for which successor states choices need to be added
+                        for (auto const& successorEntry : mdp.getTransitionMatrix().getRow(smallestCommandChoice)) {
+                            if (!storm::utility::isZero(successorEntry.getValue())) {
+                                if (!reachableProb0EStates.get(successorEntry.getColumn())) {
+                                    reachableProb0EStates.set(successorEntry.getColumn());
+                                    prob0EWorklist.push(successorEntry.getColumn());
+                                }
+                            }
+                        }
+                    }
+                }
+                
+                auto endTime = std::chrono::high_resolution_clock::now();
+                std::cout << std::endl << "Extended command for lower bounded property to size " << commandSet.size() << " in " << std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count() << "ms." << std::endl;
+            }
+            
+            static std::shared_ptr<PrismHighLevelCounterexample> computeCounterexample(Environment const& env, storm::prism::Program program, storm::models::sparse::Mdp<T> const& mdp, std::shared_ptr<storm::logic::Formula const> const& formula) {
 #ifdef STORM_HAVE_Z3
                 std::cout << std::endl << "Generating minimal label counterexample for formula " << *formula << std::endl;
                 
                 STORM_LOG_THROW(formula->isProbabilityOperatorFormula(), storm::exceptions::InvalidPropertyException, "Counterexample generation does not support this kind of formula. Expecting a probability operator as the outermost formula element.");
                 storm::logic::ProbabilityOperatorFormula const& probabilityOperator = formula->asProbabilityOperatorFormula();
                 STORM_LOG_THROW(probabilityOperator.hasBound(), storm::exceptions::InvalidPropertyException, "Counterexample generation only supports bounded formulas.");
-                storm::logic::ComparisonType comparisonType = probabilityOperator.getComparisonType();
-                STORM_LOG_THROW(comparisonType == storm::logic::ComparisonType::Less || comparisonType == storm::logic::ComparisonType::LessEqual, storm::exceptions::InvalidPropertyException, "Counterexample generation only supports formulas with an upper probability bound.");
                 STORM_LOG_THROW(probabilityOperator.getSubformula().isUntilFormula() || probabilityOperator.getSubformula().isEventuallyFormula(), storm::exceptions::InvalidPropertyException, "Path formula is required to be of the form 'phi U psi' for counterexample generation.");
-                
+
+                storm::logic::ComparisonType comparisonType = probabilityOperator.getComparisonType();
                 bool strictBound = comparisonType == storm::logic::ComparisonType::Less;
-                double threshold = probabilityOperator.getThresholdAs<double>();
+                double threshold = probabilityOperator.getThresholdAs<T>();
                 
                 storm::storage::BitVector phiStates;
                 storm::storage::BitVector psiStates;
                 storm::modelchecker::SparseMdpPrctlModelChecker<storm::models::sparse::Mdp<T>> modelchecker(mdp);
                 
                 if (probabilityOperator.getSubformula().isUntilFormula()) {
+                    STORM_LOG_THROW(!storm::logic::isLowerBound(comparisonType), storm::exceptions::NotSupportedException, "Lower bounds in counterexamples are only supported for eventually formulas.");
                     storm::logic::UntilFormula const& untilFormula = probabilityOperator.getSubformula().asUntilFormula();
                     
                     std::unique_ptr<storm::modelchecker::CheckResult> leftResult = modelchecker.check(env, untilFormula.getLeftSubformula());
@@ -1851,12 +1897,43 @@ namespace storm {
                     psiStates = subQualitativeResult.getTruthValuesVector();
                 }
                 
+                bool lowerBoundedFormula = false;
+                if (storm::logic::isLowerBound(comparisonType)) {
+                    // If the formula specifies a lower bound, we need to modify the phi and psi states.
+                    // More concretely, we convert P(min)>lambda(F psi) to P(max)<(1-lambda)(G !psi) = P(max)<(1-lambda)(!psi U prob0E(psi))
+                    // where prob0E(psi) is the set of states for which there exists a strategy \sigma_0 that avoids
+                    // reaching psi states completely.
+
+                    // This means that from all states in prob0E(psi) we need to include labels such that \sigma_0
+                    // is actually included in the resulting MDP. This prevents us from guaranteeing the minimality of
+                    // the returned counterexample, so we warn about that.
+                    STORM_LOG_WARN("Generating counterexample for lower-bounded property. The resulting command set need not be minimal.");
+
+                    // Modify bound appropriately.
+                    comparisonType = storm::logic::invertPreserveStrictness(comparisonType);
+                    threshold = storm::utility::one<T>() - threshold;
+                    
+                    // Modify the phi and psi states appropriately.
+                    storm::storage::BitVector statesWithProbability0E = storm::utility::graph::performProb0E(mdp.getTransitionMatrix(), mdp.getTransitionMatrix().getRowGroupIndices(), mdp.getBackwardTransitions(), phiStates, psiStates);
+                    phiStates = ~psiStates;
+                    psiStates = std::move(statesWithProbability0E);
+
+                    // Remember our transformation so we can add commands to guarantee that the prob0E(a) states actually
+                    // have a strategy that voids a states.
+                    lowerBoundedFormula = true;
+                }
+                
                 // Delegate the actual computation work to the function of equal name.
                 auto startTime = std::chrono::high_resolution_clock::now();
                 auto commandSet = getMinimalCommandSet(env, program, mdp, phiStates, psiStates, threshold, strictBound, true, storm::settings::getModule<storm::settings::modules::CounterexampleGeneratorSettings>().isEncodeReachabilitySet());
                 auto endTime = std::chrono::high_resolution_clock::now();
                 std::cout << std::endl << "Computed minimal command set of size " << commandSet.size() << " in " << std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count() << "ms." << std::endl;
                 
+                // Extend the command set properly.
+                if (lowerBoundedFormula) {
+                    extendCommandSetLowerBound(mdp, commandSet, phiStates, psiStates);
+                }
+                
                 return std::make_shared<PrismHighLevelCounterexample>(program.restrictCommands(commandSet));
 
 #else
diff --git a/src/storm/logic/ComparisonType.h b/src/storm/logic/ComparisonType.h
index bf932fc2b..19475960f 100644
--- a/src/storm/logic/ComparisonType.h
+++ b/src/storm/logic/ComparisonType.h
@@ -4,32 +4,45 @@
 #include <iostream>
 
 namespace storm {
-	namespace logic {
-            enum class ComparisonType { Less, LessEqual, Greater, GreaterEqual };
-
-            inline bool isStrict(ComparisonType t) {
-                return (t == ComparisonType::Less || t == ComparisonType::Greater);
-            }
-
-            inline bool isLowerBound(ComparisonType t) {
-                return (t == ComparisonType::Greater || t == ComparisonType::GreaterEqual);
+    namespace logic {
+        enum class ComparisonType { Less, LessEqual, Greater, GreaterEqual };
+        
+        inline bool isStrict(ComparisonType t) {
+            return (t == ComparisonType::Less || t == ComparisonType::Greater);
+        }
+        
+        inline bool isLowerBound(ComparisonType t) {
+            return (t == ComparisonType::Greater || t == ComparisonType::GreaterEqual);
+        }
+        
+        inline ComparisonType invert(ComparisonType t) {
+            switch(t) {
+                case ComparisonType::Less:
+                    return ComparisonType::GreaterEqual;
+                case ComparisonType::LessEqual:
+                    return ComparisonType::Greater;
+                case ComparisonType::Greater:
+                    return ComparisonType::LessEqual;
+                case ComparisonType::GreaterEqual:
+                    return ComparisonType::Less;
             }
+        }
         
-            inline ComparisonType invert(ComparisonType t) {
-                switch(t) {
-                    case ComparisonType::Less:
-                        return ComparisonType::GreaterEqual;
-                    case ComparisonType::LessEqual:
-                        return ComparisonType::Greater;
-                    case ComparisonType::Greater:
-                        return ComparisonType::LessEqual;
-                    case ComparisonType::GreaterEqual:
-                        return ComparisonType::Less;
-                }
+        inline ComparisonType invertPreserveStrictness(ComparisonType t) {
+            switch(t) {
+                case ComparisonType::Less:
+                    return ComparisonType::Greater;
+                case ComparisonType::LessEqual:
+                    return ComparisonType::GreaterEqual;
+                case ComparisonType::Greater:
+                    return ComparisonType::Less;
+                case ComparisonType::GreaterEqual:
+                    return ComparisonType::LessEqual;
             }
+        }
         
-            std::ostream& operator<<(std::ostream& out, ComparisonType const& comparisonType);
-	}
+        std::ostream& operator<<(std::ostream& out, ComparisonType const& comparisonType);
+    }
 }
 
 #endif /* STORM_LOGIC_COMPARISONTYPE_H_ */
diff --git a/src/storm/storage/sparse/PrismChoiceOrigins.cpp b/src/storm/storage/sparse/PrismChoiceOrigins.cpp
index 5e19dc654..366b9524b 100644
--- a/src/storm/storage/sparse/PrismChoiceOrigins.cpp
+++ b/src/storm/storage/sparse/PrismChoiceOrigins.cpp
@@ -111,4 +111,4 @@ namespace storm {
             }
         }
     }
-}
\ No newline at end of file
+}
diff --git a/src/storm/storage/sparse/PrismChoiceOrigins.h b/src/storm/storage/sparse/PrismChoiceOrigins.h
index ad5e50444..2a004d633 100644
--- a/src/storm/storage/sparse/PrismChoiceOrigins.h
+++ b/src/storm/storage/sparse/PrismChoiceOrigins.h
@@ -67,4 +67,4 @@ namespace storm {
             };
         }
     }
-}
\ No newline at end of file
+}

From a982af3348d0f1489d2c99b7efa7c7603742d4be Mon Sep 17 00:00:00 2001
From: TimQu <tim.quatmann@cs.rwth-aachen.de>
Date: Sat, 9 Dec 2017 15:56:20 +0100
Subject: [PATCH 19/19] setting lower bounds  for equation solvers via
 move-reference

---
 src/storm/solver/AbstractEquationSolver.cpp | 5 +++++
 src/storm/solver/AbstractEquationSolver.h   | 5 +++++
 2 files changed, 10 insertions(+)

diff --git a/src/storm/solver/AbstractEquationSolver.cpp b/src/storm/solver/AbstractEquationSolver.cpp
index 4d0a34abc..d17a9041d 100644
--- a/src/storm/solver/AbstractEquationSolver.cpp
+++ b/src/storm/solver/AbstractEquationSolver.cpp
@@ -139,6 +139,11 @@ namespace storm {
             lowerBounds = values;
         }
         
+        template<typename ValueType>
+        void AbstractEquationSolver<ValueType>::setLowerBounds(std::vector<ValueType>&& values) {
+            lowerBounds = std::move(values);
+        }
+        
         template<typename ValueType>
         void AbstractEquationSolver<ValueType>::setUpperBounds(std::vector<ValueType> const& values) {
             upperBounds = values;
diff --git a/src/storm/solver/AbstractEquationSolver.h b/src/storm/solver/AbstractEquationSolver.h
index 64eab4d1d..d5000f8c7 100644
--- a/src/storm/solver/AbstractEquationSolver.h
+++ b/src/storm/solver/AbstractEquationSolver.h
@@ -117,6 +117,11 @@ namespace storm {
              */
             void setLowerBounds(std::vector<ValueType> const& values);
             
+            /*!
+             * Sets lower bounds for the solution that can potentially be used by the solver.
+             */
+            void setLowerBounds(std::vector<ValueType>&& values);
+            
             /*!
              * Sets upper bounds for the solution that can potentially be used by the solver.
              */