From 835e668ab02068ed24f568951b44ff953acd07d2 Mon Sep 17 00:00:00 2001 From: Dave Parker Date: Wed, 1 Nov 2006 20:25:16 +0000 Subject: [PATCH] Fixed uniformisation-based methods to use epsilon/8 instead of epsilon. git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@101 bbc10eb1-c90d-0410-af57-cb519fbb1720 --- prism/src/hybrid/PH_StochBoundedUntil.cc | 11 +++++++---- prism/src/hybrid/PH_StochCumulReward.cc | 11 +++++++---- prism/src/hybrid/PH_StochTransient.cc | 11 +++++++---- prism/src/mtbdd/PM_StochBoundedUntil.cc | 11 +++++++---- prism/src/mtbdd/PM_StochCumulReward.cc | 11 +++++++---- prism/src/mtbdd/PM_StochTransient.cc | 11 +++++++---- prism/src/sparse/PS_StochBoundedUntil.cc | 11 +++++++---- prism/src/sparse/PS_StochCumulReward.cc | 11 +++++++---- prism/src/sparse/PS_StochTransient.cc | 11 +++++++---- 9 files changed, 63 insertions(+), 36 deletions(-) diff --git a/prism/src/hybrid/PH_StochBoundedUntil.cc b/prism/src/hybrid/PH_StochBoundedUntil.cc index b09abeb7..addcc1f0 100644 --- a/prism/src/hybrid/PH_StochBoundedUntil.cc +++ b/prism/src/hybrid/PH_StochBoundedUntil.cc @@ -93,7 +93,7 @@ jint mu // probs for multiplying // misc bool done; int i, j, iters, num_iters; - double x, kb, kbt, max_diag, weight; + double x, kb, kbt, max_diag, weight, term_crit_param_unif; // start clocks start1 = start2 = util_cpu_time(); @@ -190,9 +190,12 @@ jint mu // probs for multiplying // print total memory usage PH_PrintToMainLog(env, "TOTAL: [%.1f KB]\n", kbt); + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PH_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -239,7 +242,7 @@ jint mu // probs for multiplying case TERM_CRIT_ABSOLUTE: done = true; for (i = 0; i < n; i++) { - if (fabs(soln2[i] - soln[i]) > term_crit_param) { + if (fabs(soln2[i] - soln[i]) > term_crit_param_unif) { done = false; break; } @@ -248,7 +251,7 @@ jint mu // probs for multiplying case TERM_CRIT_RELATIVE: done = true; for (i = 0; i < n; i++) { - if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param) { + if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param_unif) { done = false; break; } diff --git a/prism/src/hybrid/PH_StochCumulReward.cc b/prism/src/hybrid/PH_StochCumulReward.cc index 1bce34ef..d3ca5575 100644 --- a/prism/src/hybrid/PH_StochCumulReward.cc +++ b/prism/src/hybrid/PH_StochCumulReward.cc @@ -92,7 +92,7 @@ jdouble time // time bound // misc bool done; int i, j, l, h, iters, num_iters; - double d, max_diag, weight, kb, kbt; + double d, max_diag, weight, kb, kbt, term_crit_param_unif; // start clocks start1 = start2 = util_cpu_time(); @@ -205,9 +205,12 @@ jdouble time // time bound // print total memory usage PH_PrintToMainLog(env, "TOTAL: [%.1f KB]\n", kbt); + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PH_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -267,7 +270,7 @@ jdouble time // time bound case TERM_CRIT_ABSOLUTE: done = true; for (i = 0; i < n; i++) { - if (fabs(soln2[i] - soln[i]) > term_crit_param) { + if (fabs(soln2[i] - soln[i]) > term_crit_param_unif) { done = false; break; } @@ -276,7 +279,7 @@ jdouble time // time bound case TERM_CRIT_RELATIVE: done = true; for (i = 0; i < n; i++) { - if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param) { + if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param_unif) { done = false; break; } diff --git a/prism/src/hybrid/PH_StochTransient.cc b/prism/src/hybrid/PH_StochTransient.cc index eede421a..91cd24be 100644 --- a/prism/src/hybrid/PH_StochTransient.cc +++ b/prism/src/hybrid/PH_StochTransient.cc @@ -88,7 +88,7 @@ jdouble time // time bound // misc bool done; int i, j, iters, num_iters; - double kb, kbt, max_diag, weight; + double kb, kbt, max_diag, weight, term_crit_param_unif; // start clocks start1 = start2 = util_cpu_time(); @@ -169,9 +169,12 @@ jdouble time // time bound // print total memory usage PH_PrintToMainLog(env, "TOTAL: [%.1f KB]\n", kbt); + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PH_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -218,7 +221,7 @@ jdouble time // time bound case TERM_CRIT_ABSOLUTE: done = true; for (i = 0; i < n; i++) { - if (fabs(soln2[i] - soln[i]) > term_crit_param) { + if (fabs(soln2[i] - soln[i]) > term_crit_param_unif) { done = false; break; } @@ -227,7 +230,7 @@ jdouble time // time bound case TERM_CRIT_RELATIVE: done = true; for (i = 0; i < n; i++) { - if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param) { + if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param_unif) { done = false; break; } diff --git a/prism/src/mtbdd/PM_StochBoundedUntil.cc b/prism/src/mtbdd/PM_StochBoundedUntil.cc index 1fe8bf3c..967bb76a 100644 --- a/prism/src/mtbdd/PM_StochBoundedUntil.cc +++ b/prism/src/mtbdd/PM_StochBoundedUntil.cc @@ -67,7 +67,7 @@ jint mu // probs for multiplying double time_taken, time_for_setup, time_for_iters; // misc int i, iters, num_iters; - double x, max_diag, weight, unif; + double x, max_diag, weight, unif, term_crit_param_unif; bool done, combine; // METHOD 1 or METHOD 2? (combine rate matrix and diagonals or keep separate?) @@ -171,9 +171,12 @@ jint mu // probs for multiplying PM_PrintToMainLog(env, "[nodes=%d] [%.1f Kb]\n", i, i*20.0/1024.0); } + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PM_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -253,12 +256,12 @@ jint mu // probs for multiplying // check for steady state convergence if (do_ss_detect) switch (term_crit) { case TERM_CRIT_ABSOLUTE: - if (DD_EqualSupNorm(ddman, tmp, sol, term_crit_param)) { + if (DD_EqualSupNorm(ddman, tmp, sol, term_crit_param_unif)) { done = true; } break; case TERM_CRIT_RELATIVE: - if (DD_EqualSupNormRel(ddman, tmp, sol, term_crit_param)) { + if (DD_EqualSupNormRel(ddman, tmp, sol, term_crit_param_unif)) { done = true; } break; diff --git a/prism/src/mtbdd/PM_StochCumulReward.cc b/prism/src/mtbdd/PM_StochCumulReward.cc index df4a976f..176d24d0 100644 --- a/prism/src/mtbdd/PM_StochCumulReward.cc +++ b/prism/src/mtbdd/PM_StochCumulReward.cc @@ -66,7 +66,7 @@ jdouble time // time bound // misc bool done; int i, iters, num_iters; - double max_diag, weight, unif; + double max_diag, weight, unif, term_crit_param_unif; // start clocks start1 = start2 = util_cpu_time(); @@ -120,9 +120,12 @@ jdouble time // time bound // set up sum vector sum = DD_Constant(ddman, 0); + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PM_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -172,12 +175,12 @@ jdouble time // time bound // check for steady state convergence if (do_ss_detect) switch (term_crit) { case TERM_CRIT_ABSOLUTE: - if (DD_EqualSupNorm(ddman, tmp, sol, term_crit_param)) { + if (DD_EqualSupNorm(ddman, tmp, sol, term_crit_param_unif)) { done = true; } break; case TERM_CRIT_RELATIVE: - if (DD_EqualSupNormRel(ddman, tmp, sol, term_crit_param)) { + if (DD_EqualSupNormRel(ddman, tmp, sol, term_crit_param_unif)) { done = true; } break; diff --git a/prism/src/mtbdd/PM_StochTransient.cc b/prism/src/mtbdd/PM_StochTransient.cc index 78fc582d..d339fc75 100644 --- a/prism/src/mtbdd/PM_StochTransient.cc +++ b/prism/src/mtbdd/PM_StochTransient.cc @@ -63,7 +63,7 @@ jdouble time // time bound double time_taken, time_for_setup, time_for_iters; // misc int i, iters, num_iters; - double max_diag, weight, unif; + double max_diag, weight, unif, term_crit_param_unif; bool done, combine; // METHOD 1 or METHOD 2? (combine rate matrix and diagonals or keep separate?) @@ -152,9 +152,12 @@ jdouble time // time bound PM_PrintToMainLog(env, "[nodes=%d] [%.1f Kb]\n", i, i*20.0/1024.0); } + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PM_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -229,12 +232,12 @@ jdouble time // time bound // check for steady state convergence if (do_ss_detect) switch (term_crit) { case TERM_CRIT_ABSOLUTE: - if (DD_EqualSupNorm(ddman, tmp, sol, term_crit_param)) { + if (DD_EqualSupNorm(ddman, tmp, sol, term_crit_param_unif)) { done = true; } break; case TERM_CRIT_RELATIVE: - if (DD_EqualSupNormRel(ddman, tmp, sol, term_crit_param)) { + if (DD_EqualSupNormRel(ddman, tmp, sol, term_crit_param_unif)) { done = true; } break; diff --git a/prism/src/sparse/PS_StochBoundedUntil.cc b/prism/src/sparse/PS_StochBoundedUntil.cc index 7d466fd8..0a237536 100644 --- a/prism/src/sparse/PS_StochBoundedUntil.cc +++ b/prism/src/sparse/PS_StochBoundedUntil.cc @@ -79,7 +79,7 @@ jint mu // probs for multiplying // misc bool done; int i, j, l, h, iters, num_iters; - double d, x, max_diag, weight, kb, kbt, unif; + double d, x, max_diag, weight, kb, kbt, unif, term_crit_param_unif; // start clocks start1 = start2 = util_cpu_time(); @@ -180,9 +180,12 @@ jint mu // probs for multiplying // print total memory usage PS_PrintToMainLog(env, "TOTAL: [%.1f KB]\n", kbt); + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PS_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -266,7 +269,7 @@ jint mu // probs for multiplying case TERM_CRIT_ABSOLUTE: done = true; for (i = 0; i < n; i++) { - if (fabs(soln2[i] - soln[i]) > term_crit_param) { + if (fabs(soln2[i] - soln[i]) > term_crit_param_unif) { done = false; break; } @@ -275,7 +278,7 @@ jint mu // probs for multiplying case TERM_CRIT_RELATIVE: done = true; for (i = 0; i < n; i++) { - if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param) { + if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param_unif) { done = false; break; } diff --git a/prism/src/sparse/PS_StochCumulReward.cc b/prism/src/sparse/PS_StochCumulReward.cc index 77e53648..f932efdc 100644 --- a/prism/src/sparse/PS_StochCumulReward.cc +++ b/prism/src/sparse/PS_StochCumulReward.cc @@ -77,7 +77,7 @@ jdouble time // time bound // misc bool done; int i, j, l, h, iters, num_iters; - double d, max_diag, weight, kb, kbt, unif; + double d, max_diag, weight, kb, kbt, unif, term_crit_param_unif; // start clocks start1 = start2 = util_cpu_time(); @@ -195,9 +195,12 @@ jdouble time // time bound // print total memory usage PS_PrintToMainLog(env, "TOTAL: [%.1f KB]\n", kbt); + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PS_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -293,7 +296,7 @@ jdouble time // time bound case TERM_CRIT_ABSOLUTE: done = true; for (i = 0; i < n; i++) { - if (fabs(soln2[i] - soln[i]) > term_crit_param) { + if (fabs(soln2[i] - soln[i]) > term_crit_param_unif) { done = false; break; } @@ -302,7 +305,7 @@ jdouble time // time bound case TERM_CRIT_RELATIVE: done = true; for (i = 0; i < n; i++) { - if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param) { + if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param_unif) { done = false; break; } diff --git a/prism/src/sparse/PS_StochTransient.cc b/prism/src/sparse/PS_StochTransient.cc index ae81612d..5f83e838 100644 --- a/prism/src/sparse/PS_StochTransient.cc +++ b/prism/src/sparse/PS_StochTransient.cc @@ -73,7 +73,7 @@ jdouble time // time bound // misc bool done; int i, j, l, h, iters, num_iters; - double d, max_diag, weight, kb, kbt, unif; + double d, max_diag, weight, kb, kbt, unif, term_crit_param_unif; // start clocks start1 = start2 = util_cpu_time(); @@ -158,9 +158,12 @@ jdouble time // time bound // print total memory usage PS_PrintToMainLog(env, "TOTAL: [%.1f KB]\n", kbt); + // compute new termination criterion parameter (epsilon/8) + term_crit_param_unif = term_crit_param / 8.0; + // compute poisson probabilities (fox/glynn) PS_PrintToMainLog(env, "\nUniformisation: q.t = %f x %f = %f\n", unif, time, unif * time); - fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param); + fgw = fox_glynn(unif * time, 1.0e-300, 1.0e+300, term_crit_param_unif); for (i = fgw.left; i <= fgw.right; i++) { fgw.weights[i-fgw.left] /= fgw.total_weight; } @@ -244,7 +247,7 @@ jdouble time // time bound case TERM_CRIT_ABSOLUTE: done = true; for (i = 0; i < n; i++) { - if (fabs(soln2[i] - soln[i]) > term_crit_param) { + if (fabs(soln2[i] - soln[i]) > term_crit_param_unif) { done = false; break; } @@ -253,7 +256,7 @@ jdouble time // time bound case TERM_CRIT_RELATIVE: done = true; for (i = 0; i < n; i++) { - if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param) { + if (fabs((soln2[i] - soln[i])/soln2[i]) > term_crit_param_unif) { done = false; break; }