Bug fix (CTMC cumulative rewards with rewards on self-loops) (sparse engine).

git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@829 bbc10eb1-c90d-0410-af57-cb519fbb1720

prism/src/sparse/PS_NondetReachReward.cc

@@ -82,6 +82,10 @@ jboolean min		// min or max probabilities (true = min, false = max)
 	// timing stuff
 	long start1, start2, start3, stop;
 	double time_taken, time_for_setup, time_for_iters;
+	// adversary stuff
+	bool adv = true, adv_loop = false;
+	FILE *fp_adv = NULL;
+	int adv_l, adv_h;
 	// misc
 	int i, j, k, k_r, l1, h1, l2, h2, l2_r, h2_r, iters;
 	double d1, d2, kb, kbt;
@@ -165,7 +169,13 @@ jboolean min		// min or max probabilities (true = min, false = max)
 	done = false;
 	PS_PrintToMainLog(env, "\nStarting iterations...\n");
 
-	while (!done && iters < max_iters) {
+	// open file to store adversary (if required)
+	if (adv) {
+		fp_adv = fopen("adv.tra", "w");
+		fprintf(fp_adv, "%d ?\n", n);
+	}
+	
+	while ((!done && iters < max_iters) || adv_loop) {
 	
 		iters++;
 
@@ -220,16 +230,18 @@ jboolean min		// min or max probabilities (true = min, false = max)
 					d2 += non_zeros[k] * soln[cols[k]];
 				}
 				// see if this value is the min/max so far
-				if (min) {
-					if (first | d2 < d1) d1 = d2;
-				} else {
-					if (first | d2 > d1) d1 = d2;
+				if (first || min&&(d2<d1) || !min&&(d2>d1)) {
+					d1 = d2;
+					if (adv_loop) { adv_l = l2; adv_h = h2; }
 				}
 				first = false;
 			}
 			// set vector element
 			// (if there were no choices from this state, reward is zero)
 			soln2[i] = (h1 > l1) ? d1 : 0;
+			// store adversary info (if required)
+			if (adv_loop) if (h1 > l1)
+				for (k = adv_l; k < adv_h; k++) fprintf(fp_adv, "%d %d %g\n", i, cols[k], non_zeros[k]);
 		}
 		
 		// check convergence
@@ -262,6 +274,9 @@ jboolean min		// min or max probabilities (true = min, false = max)
 		soln = soln2;
 		soln2 = tmpsoln;
 		
+		// if we're done, but adversary generation is required, go round once more
+		if (done && adv) adv_loop = !adv_loop;
+		
 //		PS_PrintToMainLog(env, "%.2f %.2f sec\n", ((double)(util_cpu_time() - start3)/1000), ((double)(util_cpu_time() - start2)/1000)/iters);
 	}
 	
@@ -281,6 +296,11 @@ jboolean min		// min or max probabilities (true = min, false = max)
 		for (i = 0; i < n; i++) if (inf_vec[i] > 0) soln[i] = HUGE_VAL;
 	}
 	
+	// close file to store adversary (if required)
+	if (adv) {
+		fclose(fp_adv);
+	}
+	
 	// free memory
 	Cudd_RecursiveDeref(ddman, a);
 	Cudd_RecursiveDeref(ddman, state_rewards);

prism/src/sparse/PS_NondetUntil.cc

@@ -76,10 +76,14 @@ jboolean min		// min or max probabilities (true = min, false = max)
 	// timing stuff
 	long start1, start2, start3, stop;
 	double time_taken, time_for_setup, time_for_iters;
+	// adversary stuff
+	bool adv = true, adv_loop = false;
+	FILE *fp_adv = NULL;
+	int adv_l, adv_h;
 	// misc
 	int i, j, k, l1, h1, l2, h2, iters;
 	double d1, d2, kb, kbt;
-	bool done;
+	bool done, first;
 	
 	// start clocks	
 	start1 = start2 = util_cpu_time();
@@ -138,7 +142,13 @@ jboolean min		// min or max probabilities (true = min, false = max)
 	done = false;
 	PS_PrintToMainLog(env, "\nStarting iterations...\n");
 	
-	while (!done && iters < max_iters) {
+	// open file to store adversary (if required)
+	if (adv) {
+		fp_adv = fopen("adv.tra", "w");
+		fprintf(fp_adv, "%d ?\n", n);
+	}
+	
+	while ((!done && iters < max_iters) || adv_loop) {
 		
 		iters++;
 		
@@ -158,6 +168,7 @@ jboolean min		// min or max probabilities (true = min, false = max)
 		h1 = h2 = 0;
 		for (i = 0; i < n; i++) {
 			d1 = min ? 2 : -1;
+			first = true;
 			if (!use_counts) { l1 = row_starts[i]; h1 = row_starts[i+1]; }
 			else { l1 = h1; h1 += row_counts[i]; }
 			for (j = l1; j < h1; j++) {
@@ -167,15 +178,18 @@ jboolean min		// min or max probabilities (true = min, false = max)
 				for (k = l2; k < h2; k++) {
 					d2 += non_zeros[k] * soln[cols[k]];
 				}
-				if (min) {
-					if (d2 < d1) d1 = d2;
-				} else {
-					if (d2 > d1) d1 = d2;
+				if (first || min&&(d2<d1) || !min&&(d2>d1)) {
+					d1 = d2;
+					if (adv_loop) { adv_l = l2; adv_h = h2; }
 				}
+				first = false;
 			}
 			// set vector element
 			// (if no choices, use value of yes)
 			soln2[i] = (h1 > l1) ? d1 : yes_vec[i];
+			// store adversary info (if required)
+			if (adv_loop) if (h1 > l1)
+				for (k = adv_l; k < adv_h; k++) fprintf(fp_adv, "%d %d %g\n", i, cols[k], non_zeros[k]);
 		}
 		
 		// check convergence
@@ -208,6 +222,9 @@ jboolean min		// min or max probabilities (true = min, false = max)
 		soln = soln2;
 		soln2 = tmpsoln;
 		
+		// if we're done, but adversary generation is required, go round once more
+		if (done && adv) adv_loop = !adv_loop;
+		
 //		PS_PrintToMainLog(env, "%.2f %.2f sec\n", ((double)(util_cpu_time() - start3)/1000), ((double)(util_cpu_time() - start2)/1000)/iters);
 	}
 	
@@ -219,6 +236,11 @@ jboolean min		// min or max probabilities (true = min, false = max)
 	// print iterations/timing info
 	PS_PrintToMainLog(env, "\nIterative method: %d iterations in %.2f seconds (average %.6f, setup %.2f)\n", iters, time_taken, time_for_iters/iters, time_for_setup);
 	
+	// close file to store adversary (if required)
+	if (adv) {
+		fclose(fp_adv);
+	}
+	
 	// free memory
 	Cudd_RecursiveDeref(ddman, a);
 	free_nd_sparse_matrix(ndsm);

prism/src/sparse/PS_StochCumulReward.cc

@@ -157,13 +157,7 @@ jdouble time		// time bound
 	}
 	
 	// combine state/transition rewards into a single vector
-	// for efficiency, we keep diagonal stuff separately as an array
-	// maths is as follows, where: c=state costs, q=unif const, P=unif matrix, C=trans costs, 1=vect of 1sm, Q/R = gen/rate matrix, D=row sums of R
-	// new state rewards = c + q(P.C)1
-	// = c + q((I+Q/q).C)1
-	// = c + q((I+(R-D)/q).C)1
-	// = c + q((I-D/q).C+(R/q.C))1
-	// = c + (q(I-D/q).C)1 + (R.C)1
+	// new state rewards = c + (R.C)1
 	// first, multiply transition rates by transition rewards and sum rows
 	// = (R.C)1
 	Cudd_Ref(trans);
@@ -174,25 +168,14 @@ jdouble time		// time bound
 	// = c + (R.C)1
 	Cudd_Ref(state_rewards);
 	tmp = DD_Apply(ddman, APPLY_PLUS, tmp, state_rewards);
-	soln2 = mtbdd_to_double_vector(ddman, tmp, rvars, num_rvars, odd);
-	Cudd_RecursiveDeref(ddman, tmp);
-	// now get diagonal transition costs
-	Cudd_Ref(trans_rewards);
-	tmp = DD_Apply(ddman, APPLY_TIMES, trans_rewards, DD_Identity(ddman, rvars, cvars, num_rvars));
-	tmp = DD_SumAbstract(ddman, tmp, cvars, num_cvars);
 	soln = mtbdd_to_double_vector(ddman, tmp, rvars, num_rvars, odd);
 	Cudd_RecursiveDeref(ddman, tmp);
-	// and then multiply by q and modified diagonals
-	// = (q(I-D/q).C)1
-	for (i = 0; i < n; i++) soln[i] *= (unif * ((!compact_d)?diags[i]:diags_dist->dist[diags_dist->ptrs[i]]));
-	// finally sum two arrays
-	// = c + (q(I-D/q).C)1 + (R.C)1
-	for (i = 0; i < n; i++) soln[i] += soln2[i];
 	
 	// create solution/iteration vectors
 	PS_PrintToMainLog(env, "Allocating iteration vectors... ");
 	// soln has already been created and initialised to rewards vector as required
-	// soln2 has also already been created; initial values unimportant
+	// need to create soln2 and sum
+	soln2 = new double[n];
 	sum = new double[n];
 	kb = n*8.0/1024.0;
 	kbt += 3*kb;