From 7c612640daa4e3a751ad87bc750154af7558071b Mon Sep 17 00:00:00 2001 From: Joachim Klein Date: Fri, 21 Jul 2017 12:51:01 +0000 Subject: [PATCH] PS iteration methods: use MeasureSupNorm. MultiObjective remains todo. Uses common code that was refactored to Measures.h git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@12079 bbc10eb1-c90d-0410-af57-cb519fbb1720 --- prism/include/Measures.h | 30 ++++++++++-------------- prism/src/sparse/PS_JOR.cc | 25 +++++++++----------- prism/src/sparse/PS_NondetReachReward.cc | 23 ++++++++---------- prism/src/sparse/PS_NondetUntil.cc | 23 ++++++++---------- prism/src/sparse/PS_Power.cc | 25 +++++++++----------- prism/src/sparse/PS_ProbTransient.cc | 19 +++++++-------- prism/src/sparse/PS_SOR.cc | 19 +++++++-------- prism/src/sparse/PS_StochBoundedUntil.cc | 19 +++++++-------- prism/src/sparse/PS_StochCumulReward.cc | 19 +++++++-------- prism/src/sparse/PS_StochTransient.cc | 19 +++++++-------- 10 files changed, 95 insertions(+), 126 deletions(-) diff --git a/prism/include/Measures.h b/prism/include/Measures.h index 8bd20617..8bb092eb 100644 --- a/prism/include/Measures.h +++ b/prism/include/Measures.h @@ -66,24 +66,18 @@ public: inline void measure(double v1, double v2) { double x; - // special case: one of the values is infinite (this can happen e.g. for non-converging - // iterations when the values grow extremely large and overflow to infinity). - if (isinf(v1) || isinf(v2)) { - x = std::numeric_limits::infinity(); - } else { - // compute absolute of difference - x = fabs(v2 - v1); - if (relative) { - // for relative mode: divide by second value - // We take the absolute value of v2 to ensure that - // x remains non-negative. v2 can become negative e.g. - // during iterations with over-relaxation with large - // omega values. - // Note: if v2 is 0, then x will become +inf for x>0 and NaN for x=0, i.e., v1=v2=0 - // In the later case, the max computation below will ignore the NaN, - // as NaN > y is false for all y - x /= fabs(v2); - } + // compute absolute of difference + x = fabs(v2 - v1); + if (relative) { + // for relative mode: divide by second value + // We take the absolute value of v2 to ensure that + // x remains non-negative. v2 can become negative e.g. + // during iterations with over-relaxation with large + // omega values. + // Note: if v2 is 0, then x will become +inf for x>0 and NaN for x=0, i.e., v1=v2=0 + // In the later case, the max computation below will ignore the NaN, + // as NaN > y is false for all y + x /= fabs(v2); } // sup_norm = max { x, sup_norm } diff --git a/prism/src/sparse/PS_JOR.cc b/prism/src/sparse/PS_JOR.cc index a02b8b5e..412d7329 100644 --- a/prism/src/sparse/PS_JOR.cc +++ b/prism/src/sparse/PS_JOR.cc @@ -36,6 +36,7 @@ #include "PrismSparseGlob.h" #include "jnipointer.h" #include "prism.h" +#include "Measures.h" #include "ExportIterations.h" #include #include @@ -87,9 +88,11 @@ jdouble omega // omega (over-relaxation parameter) double time_taken, time_for_setup, time_for_iters; // misc int i, j, l, h, iters; - double d, x, sup_norm, kb, kbt; + double d, kb, kbt; bool done; - + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); + // exception handling around whole function try { @@ -214,7 +217,7 @@ jdouble omega // omega (over-relaxation parameter) time_for_setup = (double)(stop - start2)/1000; start2 = stop; start3 = stop; - + // start iterations iters = 0; done = false; @@ -281,21 +284,15 @@ jdouble omega // omega (over-relaxation parameter) iterationExport->exportVector(soln2, n, 0); // check convergence - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param) { done = true; } - + // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { - PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, ", %.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_NondetReachReward.cc b/prism/src/sparse/PS_NondetReachReward.cc index db2cf467..2779c63c 100644 --- a/prism/src/sparse/PS_NondetReachReward.cc +++ b/prism/src/sparse/PS_NondetReachReward.cc @@ -39,6 +39,7 @@ #include "jnipointer.h" #include "ExportIterations.h" #include +#include "Measures.h" #include //------------------------------------------------------------------------------ @@ -101,9 +102,11 @@ jboolean min // min or max probabilities (true = min, false = max) int num_actions; // misc int i, j, k, k_r, l1, h1, l2, h2, l2_r, h2_r, iters; - double d1, d2, x, sup_norm, kb, kbt; + double d1, d2, kb, kbt; bool done, first; - + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); + // exception handling around whole function try { @@ -323,21 +326,15 @@ jboolean min // min or max probabilities (true = min, false = max) iterationExport->exportVector(soln2, n, 0); // check convergence - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param) { done = true; } - + // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { - PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, ", %.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_NondetUntil.cc b/prism/src/sparse/PS_NondetUntil.cc index 957a8258..c7dad9eb 100644 --- a/prism/src/sparse/PS_NondetUntil.cc +++ b/prism/src/sparse/PS_NondetUntil.cc @@ -37,6 +37,7 @@ #include "PrismNativeGlob.h" #include "PrismSparseGlob.h" #include "jnipointer.h" +#include "Measures.h" #include "ExportIterations.h" #include #include @@ -97,9 +98,11 @@ jlong _strat // strategy storage int num_actions; // misc int i, j, k, l1, h1, l2, h2, iters; - double d1, d2, x, sup_norm, kb, kbt; + double d1, d2, kb, kbt; bool done, first; - + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); + // exception handling around whole function try { @@ -212,7 +215,7 @@ jlong _strat // strategy storage time_for_setup = (double)(stop - start2)/1000; start2 = stop; start3 = stop; - + // start iterations iters = 0; done = false; @@ -283,21 +286,15 @@ jlong _strat // strategy storage iterationExport->exportVector(soln2, n, 0); // check convergence - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param) { done = true; } // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { - PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, ", %.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_Power.cc b/prism/src/sparse/PS_Power.cc index 352d2203..94d6eb1e 100644 --- a/prism/src/sparse/PS_Power.cc +++ b/prism/src/sparse/PS_Power.cc @@ -36,6 +36,7 @@ #include "PrismSparseGlob.h" #include "jnipointer.h" #include "prism.h" +#include "Measures.h" #include "ExportIterations.h" #include #include @@ -85,9 +86,11 @@ jboolean transpose // transpose A? (i.e. solve xA=x not Ax=x?) double time_taken, time_for_setup, time_for_iters; // misc int i, j, l, h, iters; - double d, x, sup_norm, kb, kbt; + double d, kb, kbt; bool done; - + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); + // exception handling around whole function try { @@ -162,7 +165,7 @@ jboolean transpose // transpose A? (i.e. solve xA=x not Ax=x?) time_for_setup = (double)(stop - start2)/1000; start2 = stop; start3 = stop; - + // start iterations iters = 0; done = false; @@ -223,21 +226,15 @@ jboolean transpose // transpose A? (i.e. solve xA=x not Ax=x?) iterationExport->exportVector(soln2, n, 0); // check convergence - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param) { done = true; } - + // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { - PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, (measure.isRelative()?"relative ":""), measure.value()); PS_PrintToMainLog(env, ", %.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_ProbTransient.cc b/prism/src/sparse/PS_ProbTransient.cc index 032f16af..f3b778d4 100644 --- a/prism/src/sparse/PS_ProbTransient.cc +++ b/prism/src/sparse/PS_ProbTransient.cc @@ -36,6 +36,7 @@ #include "sparse.h" #include "PrismSparseGlob.h" #include "jnipointer.h" +#include "Measures.h" #include //------------------------------------------------------------------------------ @@ -77,7 +78,9 @@ jint time // time // misc bool done; int i, j, l, h, iters; - double d, x, sup_norm, kb, kbt; + double d, kb, kbt; + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); // exception handling around whole function try { @@ -185,15 +188,9 @@ jint time // time // check for steady state convergence if (do_ss_detect) { - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param) { done = true; } } @@ -201,7 +198,7 @@ jint time // time // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { PS_PrintToMainLog(env, "Iteration %d (of %d): ", iters, time); - if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, "%.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_SOR.cc b/prism/src/sparse/PS_SOR.cc index a0d7fe76..ace9666b 100644 --- a/prism/src/sparse/PS_SOR.cc +++ b/prism/src/sparse/PS_SOR.cc @@ -36,6 +36,7 @@ #include "PrismSparseGlob.h" #include "jnipointer.h" #include "prism.h" +#include "Measures.h" #include "ExportIterations.h" #include #include @@ -88,8 +89,10 @@ jboolean forwards // forwards or backwards? double time_taken, time_for_setup, time_for_iters; // misc int i, j, fb, l, h, iters; - double d, x, sup_norm, kb, kbt; + double d, kb, kbt; bool done; + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); // exception handling around whole function try { @@ -214,7 +217,7 @@ jboolean forwards // forwards or backwards? time_for_setup = (double)(stop - start2)/1000; start2 = stop; start3 = stop; - + // start iterations iters = 0; done = false; @@ -224,7 +227,7 @@ jboolean forwards // forwards or backwards? iters++; - sup_norm = 0.0; + measure.reset(); // store local copies of stuff double *non_zeros; @@ -282,11 +285,7 @@ jboolean forwards // forwards or backwards? } // compute norm for convergence // (note we must do this inside the loop because we only store one vector for sor/gauss-seidel) - x = fabs(d - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= d; - } - if (x > sup_norm) sup_norm = x; + measure.measure(soln[i], d); // set vector element soln[i] = d; } @@ -295,13 +294,13 @@ jboolean forwards // forwards or backwards? iterationExport->exportVector(soln, n, 0); // check convergence - if (sup_norm < term_crit_param) { + if (measure.value() < term_crit_param) { done = true; } // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { - PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + PS_PrintToMainLog(env, "Iteration %d: max %sdiff=%f", iters, measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, ", %.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_StochBoundedUntil.cc b/prism/src/sparse/PS_StochBoundedUntil.cc index 6b11d46e..b7c78501 100644 --- a/prism/src/sparse/PS_StochBoundedUntil.cc +++ b/prism/src/sparse/PS_StochBoundedUntil.cc @@ -37,6 +37,7 @@ #include "PrismSparseGlob.h" #include "jnipointer.h" #include "prism.h" +#include "Measures.h" #include //------------------------------------------------------------------------------ @@ -88,7 +89,9 @@ jlong __jlongpointer mu // probs for multiplying bool done; int j, l, h; long i, iters, num_iters; - double d, x, sup_norm, max_diag, weight, kb, kbt, unif, term_crit_param_unif; + double d, x, max_diag, weight, kb, kbt, unif, term_crit_param_unif; + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); // exception handling around whole function try { @@ -276,15 +279,9 @@ jlong __jlongpointer mu // probs for multiplying // check for steady state convergence if (do_ss_detect) { - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param_unif) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param_unif) { done = true; } } @@ -310,7 +307,7 @@ jlong __jlongpointer mu // probs for multiplying // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { PS_PrintToMainLog(env, "Iteration %d (of %d): ", iters, fgw.right); - if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, "%.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_StochCumulReward.cc b/prism/src/sparse/PS_StochCumulReward.cc index 52bd606e..1c142aa0 100644 --- a/prism/src/sparse/PS_StochCumulReward.cc +++ b/prism/src/sparse/PS_StochCumulReward.cc @@ -36,6 +36,7 @@ #include "sparse.h" #include "PrismSparseGlob.h" #include "jnipointer.h" +#include "Measures.h" #include //------------------------------------------------------------------------------ @@ -85,7 +86,9 @@ jdouble time // time bound bool done; int j, l, h; long i, iters, num_iters; - double d, x, sup_norm, max_diag, weight, kb, kbt, unif, term_crit_param_unif; + double d, max_diag, weight, kb, kbt, unif, term_crit_param_unif; + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); // exception handling around whole function try { @@ -285,15 +288,9 @@ jdouble time // time bound // check for steady state convergence if (do_ss_detect) { - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param_unif) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param_unif) { done = true; } } @@ -320,7 +317,7 @@ jdouble time // time bound // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { PS_PrintToMainLog(env, "Iteration %d (of %d): ", iters, fgw.right); - if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, "%.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); } diff --git a/prism/src/sparse/PS_StochTransient.cc b/prism/src/sparse/PS_StochTransient.cc index abdfa89c..938350d8 100644 --- a/prism/src/sparse/PS_StochTransient.cc +++ b/prism/src/sparse/PS_StochTransient.cc @@ -36,6 +36,7 @@ #include "sparse.h" #include "PrismSparseGlob.h" #include "jnipointer.h" +#include "Measures.h" #include //------------------------------------------------------------------------------ @@ -81,7 +82,9 @@ jdouble time // time bound bool done; int j, l, h; long i, iters, num_iters; - double d, x, sup_norm, max_diag, weight, kb, kbt, unif, term_crit_param_unif; + double d, max_diag, weight, kb, kbt, unif, term_crit_param_unif; + // measure for convergence termination check + MeasureSupNorm measure(term_crit == TERM_CRIT_RELATIVE); // exception handling around whole function try { @@ -255,15 +258,9 @@ jdouble time // time bound // check for steady state convergence if (do_ss_detect) { - sup_norm = 0.0; - for (i = 0; i < n; i++) { - x = fabs(soln2[i] - soln[i]); - if (term_crit == TERM_CRIT_RELATIVE) { - x /= soln2[i]; - } - if (x > sup_norm) sup_norm = x; - } - if (sup_norm < term_crit_param_unif) { + measure.reset(); + measure.measure(soln, soln2, n); + if (measure.value() < term_crit_param_unif) { done = true; } } @@ -289,7 +286,7 @@ jdouble time // time bound // print occasional status update if ((util_cpu_time() - start3) > UPDATE_DELAY) { PS_PrintToMainLog(env, "Iteration %d (of %d): ", iters, fgw.right); - if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", (term_crit == TERM_CRIT_RELATIVE)?"relative ":"", sup_norm); + if (do_ss_detect) PS_PrintToMainLog(env, "max %sdiff=%f, ", measure.isRelative()?"relative ":"", measure.value()); PS_PrintToMainLog(env, "%.2f sec so far\n", ((double)(util_cpu_time() - start2)/1000)); start3 = util_cpu_time(); }