Merge branch 'foxglynn'

src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp

@@ -632,19 +632,21 @@ 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);
-                STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult));
+//                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);
+                
+                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;
                     }
@@ -654,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) {
@@ -672,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; };
                     
@@ -689,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);
                 }
                 

src/storm/utility/numerical.cpp

@@ -0,0 +1,275 @@
+#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/PrecisionExceededException.h"
+
+namespace storm {
+    namespace utility {
+        namespace numerical {
+
+            template<typename ValueType>
+            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) {
+                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 = log(tau);
+                
+                // In error bound comparisons, we always compare with epsilon*sqrt_2_pi.
+                epsilon *= sqrt_2_pi;
+                
+                // 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.
+                    left = 0;
+                } 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;
+                    
+                    // 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));
+                        
+                        // For small lambda the above calculation can yield negative truncation points, crop them here.
+                        if (left <= 0) {
+                            left = 0;
+                            break;
+                        }
+                        
/* /////////////////////////////////////////////////////////////////////////
 * File:        rangelib/error/exceptions.hpp
 *
 * Purpose:     Range exceptions.
 *
 * Created:     30th December 2005
 * Updated:     10th August 2009
 *
 * Home:        http://stlsoft.org/
 *
 * Copyright (c) 2005-2009, Matthew Wilson and Synesis Software
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * - Redistributions of source code must retain the above copyright notice, this
 *   list of conditions and the following disclaimer.
 * - Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * - Neither the name(s) of Matthew Wilson and Synesis Software nor the names of
 *   any contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 * ////////////////////////////////////////////////////////////////////// */
/** \file rangelib/error/exceptions.hpp
 *
 * \brief [C++ only] Definition of the rangelib::range_exception and
 *   rangelib::empty_range_exception exception classes.
 *
 *   (\ref group__library__error "Error" Library).
 */
#ifndef RANGELIB_INCL_RANGELIB_ERROR_HPP_EXCEPTIONS
#define RANGELIB_INCL_RANGELIB_ERROR_HPP_EXCEPTIONS
#ifndef STLSOFT_DOCUMENTATION_SKIP_SECTION
# define RANGELIB_VER_RANGELIB_ERROR_HPP_EXCEPTIONS_MAJOR     2
# define RANGELIB_VER_RANGELIB_ERROR_HPP_EXCEPTIONS_MINOR     0
# define RANGELIB_VER_RANGELIB_ERROR_HPP_EXCEPTIONS_REVISION  2
# define RANGELIB_VER_RANGELIB_ERROR_HPP_EXCEPTIONS_EDIT      17
#endif /* !STLSOFT_DOCUMENTATION_SKIP_SECTION */
/* /////////////////////////////////////////////////////////////////////////
 * Auto-generation and compatibility
 */
/*
[Incompatibilies-start]
[Incompatibilies-end]
[DocumentationStatus:Ready]
 */
/* /////////////////////////////////////////////////////////////////////////
 * Includes
 */
#ifndef RANGELIB_INCL_RANGELIB_HPP_RANGELIB
# include <rangelib/rangelib.hpp>
#endif /* !RANGELIB_INCL_RANGELIB_HPP_RANGELIB */
#ifndef STLSOFT_INCL_STLSOFT_UTIL_HPP_EXCEPTION_STRING
# include <stlsoft/util/exception_string.hpp>
#endif /* !STLSOFT_INCL_STLSOFT_UTIL_HPP_EXCEPTION_STRING */
#ifndef STLSOFT_INCL_STDEXCEPT
# define STLSOFT_INCL_STDEXCEPT
# include <stdexcept>
#endif /* !STLSOFT_INCL_STDEXCEPT */
#ifdef STLSOFT_UNITEST
# include <string.h>
#endif /* STLSOFT_UNITEST */
/* /////////////////////////////////////////////////////////////////////////
 * Namespace
 */
#ifndef RANGELIB_NO_NAMESPACE
# if defined(_STLSOFT_NO_NAMESPACE) || \
     defined(STLSOFT_DOCUMENTATION_SKIP_SECTION)
/* There is no stlsoft namespace, so must define ::rangelib */
namespace rangelib
{
# else
/* Define stlsoft::rangelib_project */
namespace stlsoft
{
namespace rangelib_project
{
# endif /* _STLSOFT_NO_NAMESPACE */
#endif /* !RANGELIB_NO_NAMESPACE */
/* /////////////////////////////////////////////////////////////////////////
 * Classes
 */
/** \brief General exception class for RangeLib failures.
 *
 * \ingroup group__library__error
 *
 */
class range_exception
#if defined(STLSOFT_COMPILER_IS_DMC)
    : public std::exception
#else /* ? compiler */
    : public stlsoft_ns_qual_std(exception)
#endif /* compiler */
{
/// \name Types
/// @{
private:
    typedef stlsoft_ns_qual(exception_string)   string_type;
public:
    typedef range_exception                     class_type;
#if defined(STLSOFT_COMPILER_IS_DMC)
    typedef std::exception                      parent_class_type;
#else /* ? compiler */
    typedef stlsoft_ns_qual_std(exception)      parent_class_type;
#endif /* compiler */
/// @}
/// \name Construction
/// @{
public:
    ss_explicit_k range_exception(char const* reason = NULL)
        : m_reason((NULL == reason) ? "" : reason)
    {}
    /// Destructor
    ///
    /// \note This does not do have any implementation, but is required to placate
    /// the Comeau and GCC compilers, which otherwise complain about mismatched
    /// exception specifications between this class and its parent
    virtual ~range_exception() stlsoft_throw_0()
    {}
/// @}
/// \name Accessors
/// @{
public:
    virtual char const* what() const stlsoft_throw_0()
    {
        return m_reason.empty() ? this->real_what_() : m_reason.c_str();
    }
/// @}
/// \name Implementation
/// @{
private:
    virtual char const* real_what_() const throw()
    {
        return "Range exception";
    }
/// @}
/// \name Members
/// @{
private:
    const string_type   m_reason;
/// @}
/// \name Not to be implemented
/// @{
private:
    class_type& operator =(class_type const&);
/// @}
};
/** \brief Indicates that an operation requiring non-empty range was invoked
 *    on an empty range.
 *
 * \ingroup group__library__error
 *
 */
class empty_range_exception
    : public range_exception
{
public:
    typedef range_exception         parent_class_type;
    typedef empty_range_exception  class_type;
/// \name Construction
/// @{
public:
    ss_explicit_k empty_range_exception(char const* reason = NULL)
        : parent_class_type(reason)
    {}
    /// Destructor
    ///
    /// \note This does not do have any implementation, but is required to placate
    /// the Comeau and GCC compilers, which otherwise complain about mismatched
    /// exception specifications between this class and its parent
    virtual ~empty_range_exception() stlsoft_throw_0()
    {}
/// @}
/// \name Implementation
/// @{
private:
    virtual char const* real_what_() const throw()
    {
        return "Range was empty";
    }
/// @}
};
////////////////////////////////////////////////////////////////////////////
// Unit-testing
#ifdef STLSOFT_UNITTEST
# include "./unittest/exceptions_unittest_.h"
#endif /* STLSOFT_UNITTEST */
/* ////////////////////////////////////////////////////////////////////// */
#ifndef RANGELIB_NO_NAMESPACE
# if defined(_STLSOFT_NO_NAMESPACE) || \
     defined(STLSOFT_DOCUMENTATION_SKIP_SECTION)
} // namespace rangelib
# else
} // namespace rangelib_project
} // namespace stlsoft
# endif /* _STLSOFT_NO_NAMESPACE */
#endif /* !RANGELIB_NO_NAMESPACE */
/* ////////////////////////////////////////////////////////////////////// */
#endif /* !RANGELIB_INCL_RANGELIB_ERROR_HPP_EXCEPTIONS */
/* ///////////////////////////// end of file //////////////////////////// */