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

Made Fox-Glynn (hopefully) work. 

Former-commit-id: b07c88d122
main
dehnert 11 years ago
parent
3434405cf4
commit
cde9786dfa
2 changed files with 79 additions and 93 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. 13
      src/utility/cli.h
  2. 159
      src/utility/numerical.h

13
src/utility/cli.h
View File

@ -75,6 +75,9 @@ log4cplus::Logger printer;
#include "src/exceptions/InvalidSettingsException.h"
#include "src/exceptions/InvalidTypeException.h"
// FIXME: remove this
#include "src/utility/numerical.h"
namespace storm {
namespace utility {
namespace cli {
@ -507,6 +510,16 @@ namespace storm {
initializeFileLogging();
}
auto fgresult = storm::utility::numerical::getFoxGlynnCutoff<double>(500, 1e-300, 1e+300, 1e-6);
std::cout << "left: " << std::get<0>(fgresult) << " and right: " << std::get<1>(fgresult) << std::endl;
std::cout << "weight: " << std::get<2>(fgresult) << std::endl;
int pos = 0;
for (auto const& element : std::get<3>(fgresult)) {
std::cout << "elem[" << pos << "]: " << element << std::endl;
++pos;
}
storm::settings::modules::GeneralSettings const& settings = storm::settings::generalSettings();
// If we have to build the model from a symbolic representation, we need to parse the representation first.

159
src/utility/numerical.h
View File

@ -3,32 +3,29 @@
#include <cmath>
#include "src/utility/macros.h"
#include "src/utility/ConstantsComparator.h"
#include "src/utility/constants.h"
#include "src/exceptions/InvalidArgumentException.h"
namespace storm {
namespace utility {
namespace numerical {
template<typename ValueType>
std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> getFoxGlynnCutoff(ValueType lambda, ) {
std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> getFoxGlynnCutoff(ValueType lambda, ValueType underflow, ValueType overflow, ValueType accuracy) {
storm::utility::ConstantsComparator<ValueType> comparator;
ValueType underflow, overflow, accuracy = 0;
STORM_LOG_THROW(!comparator.isZero(lambda), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero.");
// This code is more or less taken from PRISM. According to their implementation, for lambda smaller than
// 400, we compute the result using the naive method.
if (lambda < 400) {
// 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 < 25) {
ValueType eToPowerMinusLambda = std::exp(-lambda);
ValueType targetValue = (1 - (this->accuracy / 2.0)) / eToPowerMinusLambda;
ValueType targetValue = (1 - accuracy) / eToPowerMinusLambda;
std::vector<ValueType> weights;
ValueType exactlyKEvents = 1;
ValueType atMostKEvents = exactlyKEvents;
weights.push_back(exactlyKEvents * eToPowerMinusLambda);
uint_fast64_t currentK = 1;
uint_fast64_t k = 1;
do {
exactlyKEvents *= lambda / k;
atMostKEvents += exactlyKEvents;
@ -40,118 +37,94 @@ namespace storm {
} 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;
ValueType m = std::floor(lambda);
// 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.
Type sqrtpi = 1.77245385090551602729;
ValueType sqrtpi = 1.77245385090551602729;
// Square root of 2.
Type sqrt2 = 1.41421356237309504880;
ValueType sqrt2 = 1.41421356237309504880;
// Square root of lambda.
Type sqrtLambda = std::sqrt(lambda);
Type a_Lambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2; // a_\lambda.
Type b_Lambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda); // b_\lambda.
}
}
// Use Fox Glynn algorithm for lambda>400.
if (accuracy < 1e-10) {
LOG4CPLUS_ERROR(logger, "Given Value accuracy must at least be 1e-10.");
throw storm::exceptions::InvalidArgumentException("Error while computing FoxGlynn values. accuracy < 1e-10.");
}
// Factor from Fox&Glynns paper. The paper does not explain where it comes from.
Type factor = 1e+10;
// Run the FINDER algorithm.
Type m = floor(lambda);
{
// Factores used by the corollaries explained in Fox&Glynns paper.
Type sqrtpi = 1.77245385090551602729; // Squareroot of PI.
Type sqrt2 = 1.41421356237309504880; // Squareroot of 2.
Type sqrtLambda = sqrt(lambda); // Sqareroot of Lambda.
Type a_Lambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2; // a_\lambda.
Type b_Lambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda); // b_\lambda.
// 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) * exp(0.0625) * sqrt2;
bLambda = (1.0 + 1.0 / 400.0) * exp(0.125 / 400.0);
} else {
sqrtLambda = std::sqrt(lambda);
aLambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2;
bLambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda);
}
// Use Corollary 1 from the paper to find the right truncation point.
Type k = 4;
res = a_Lambda*dkl*exp(-k*k / 2.0) / (k*sqrt2*sqrtpi);
/* Normally: Iterate between the bounds stated above to find the right truncationpoint.
* This is a modification to the Fox-Glynn paper. The search for the right truncationpoint is only
* terminated by the error condition and not by the upper bound. Thus the calculation can be more accurate.
*/
while (res > accuracy / 2.0)
{
k++;
Type dkl = 1.0 / (1 - exp(-(2.0 / 9.0)*(k*sqrt2*sqrtLambda + 1.5))); // d(k,Lambda) from the paper.
Type res = a_Lambda*dkl*exp(-k*k / 2.0) / (k*sqrt2*sqrtpi); // Right hand side of the equation in Corollary 1.
uint_fast64_t k = 3;
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.
this->rightTrunk = (int)ceil(m + k*sqrt2*sqrtLambda + 1.5);
rightTruncationPoint = static_cast<uint_fast64_t>(std::ceil((m + k * sqrt2 * sqrtLambda + 1.5)));
// Use Corollary 2 to find left truncation point.
Type res;
k = 4;
k = 3;
do
{
res = b_Lambda*exp(-k*-k / 2.0) / (k*sqrt2*sqrtpi); // Right hand side of the equation in Corollary 2
k++;
} while (res > accuracy / 2.0);
// Right hand side of the equation in Corollary 2.
while ((accuracy / 2.0) < ((bLambda * exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) {
++k;
}
this->leftTrunk = (int)(m - k*sqrtLambda - 1.5); // Left hand side of the equation in Corollary 2.
// Left hand side of the equation in Corollary 2.
leftTruncationPoint = std::max(static_cast<uint_fast64_t>(std::trunc(m - k * sqrtLambda - 1.5)), uint_fast64_t(0));
// Check for underflow.
if ((int)(m - k*sqrtLambda - 1.5) < 0.0) {
LOG4CPLUS_ERROR(logger, "Underflow while computing left truncation point.");
throw storm::exceptions::OutOfRangeException("Error while computing FoxGlynn values. Underflow of left Truncation point.");
}
STORM_LOG_THROW(static_cast<uint_fast64_t>(std::trunc((m - k * sqrtLambda - 1.5))) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point.");
}
//End of FINDER algorithm.
// Use WEIGHTER algorithm to determine weights.
// Down from m
for (Type j = m; j > this->leftTrunk; j--)
this->weights.at(j - 1 - this->leftTrunk) = (j / lambda)*this->weights.at(j - this->leftTrunk);
// Up from m
for (Type j = m; j < this->rightTrunk; j++)
this->weights.at(j + 1 - this->leftTrunk) = (lambda / (j + 1))*this->weights.at(j - this->leftTrunk);
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];
}
// Compute total weight.
// Add up weights from smallest to largest to avoid roundoff errors.
this->totalWeight = 0.0;
Type s = this->leftTrunk;
Type t = this->rightTrunk;
while (s < t)
{
if (this->weights.at(s - this->leftTrunk) <= this->weights.at(t - this->leftTrunk))
{
this->totalWeight += this->weights.at(s - this->leftTrunk);
s++;
}
else
{
this->totalWeight += this->weights.at(t - this->leftTrunk);
t--;
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;
}
}
this->totalWeight += this->weights.at(s - this->leftTrunk);
totalWeight += weights[s - leftTruncationPoint];
LOG4CPLUS_INFO(logger, "Left truncationpoint: " << this->leftTrunk << ".");
LOG4CPLUS_INFO(logger, "Right truncationpoint: " << this->rightTrunk << ".");
LOG4CPLUS_INFO(logger, "Total Weight:" << this->totalWeight << ".");
LOG4CPLUS_INFO(logger, "10. Weight: " << this->weights.at(9) << ".");
return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights);
}
}
}
}
Write Preview
Loading…
Cancel
Save