sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

started to rework Fox-Glynn

tempestpy_adaptions
dehnert 7 years ago
parent
78fb458154
commit
b0fd3c1730
2 changed files with 151 additions and 1 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 4
      src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
  2. 148
      src/storm/utility/numerical.h

4
src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
View File

@ -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.

148
src/storm/utility/numerical.h
View File

@ -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);
}
}
}
}
Write Preview
Loading…
Cancel
Save