//==============================================================================
//	
//	Copyright (c) 2002-
//	Authors:
//	* Dave Parker <david.parker@comlab.ox.ac.uk> (University of Oxford)
//	* Hongyang Qu <hongyang.qu@comlab.ox.ac.uk> (University of Oxford)
//	
//------------------------------------------------------------------------------
//	
//	This file is part of PRISM.
//	
//	PRISM is free software; you can redistribute it and/or modify
//	it under the terms of the GNU General Public License as published by
//	the Free Software Foundation; either version 2 of the License, or
//	(at your option) any later version.
//	
//	PRISM is distributed in the hope that it will be useful,
//	but WITHOUT ANY WARRANTY; without even the implied warranty of
//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	See the
//	GNU General Public License for more details.
//	
//	You should have received a copy of the GNU General Public License
//	along with PRISM; if not, write to the Free Software Foundation,
//	Inc., 59 Temple Place, Suite 330, Boston, MA	02111-1307	USA
//	
//==============================================================================

// includes
#include "PrismSparse.h"
#include <math.h>
#include <util.h>
#include <cudd.h>
#include <dd.h>
#include <odd.h>
#include <dv.h>
#include "sparse.h"
#include "PrismSparseGlob.h"
#include "jnipointer.h"
#include <new>
#include "eclipse.h"

typedef double REAL;

//------------------------------------------------------------------------------

JNIEXPORT jdouble __jlongpointer JNICALL Java_sparse_PrismSparse_PS_1NondetMultiReach
(
		JNIEnv *env,
		jclass cls,
		jlong __jlongpointer t, // trans matrix
		jlong __jlongpointer od, // odd
		jlong __jlongpointer rv, // row vars
		jint num_rvars,
		jlong __jlongpointer cv, // col vars
		jint num_cvars,
		jlong __jlongpointer ndv, // nondet vars
		jint num_ndvars,
		jlongArray _targets, // target state sets
		jintArray _relops, // target relops (0:=?, 1:>, 2:>=)
		jdoubleArray _bounds, // target probability bounds
		jlong __jlongpointer m, // 'maybe' states
		jlong __jlongpointer _start, // initial state(s)
		jint checkReach, // reachability analysis or LTL model checking

		jlong __jlongpointer m_r, // 'maybe' states for the reward formula
		jlong __jlongpointer trr // transition rewards
		// Reward paramters for more complicated cases
		//jlongArray _rewards,	// sets of rewards for reward formulas
		//jlongArray _targets_r,	// target state sets for reward formulas
		//jdoubleArray _bounds_r,	// target probability bounds for reward formulas
		//jlongArray _maybe_r,	// maybe state sets for reward formulas
)
{
	// cast function parameters
	DdNode *trans = jlong_to_DdNode(t); // trans matrix
	ODDNode *odd = jlong_to_ODDNode(od); // reachable states
	DdNode **rvars = jlong_to_DdNode_array(rv); // row vars
	DdNode **cvars = jlong_to_DdNode_array(cv); // col vars
	DdNode **ndvars = jlong_to_DdNode_array(ndv); // nondet vars
	DdNode *maybe = jlong_to_DdNode(m); // 'maybe' states
	//DdNode *maybe_r = jlong_to_DdNode(m_r); 	// 'maybe' states for the reward formula
	DdNode *start = jlong_to_DdNode(_start); // initial state(s)
	DdNode *trans_rewards = NULL; // transition rewards

	// target info	
	jlong *target_ptrs = NULL;
	DdNode **targets = NULL;
	jint *relops = NULL;
	double *bounds = NULL;
	// mtbdds
	DdNode *a = NULL;
	// model stats
	int num_targets, n, nc, nc_r;
	long nnz, nnz_r;
	// sparse matrix
	NDSparseMatrix *ndsm = NULL, *ndsm_r = NULL;
	// vectors
	double *yes1_vec = NULL, *yes2_vec = NULL, *maybe_vec = NULL/*, *maybe_vec_r = NULL*/;
	// timing stuff
	long start1, start2, start3, stop, stop2;
	double time_taken, time_for_setup, time_for_lp;
	// lp stuff
	REAL *arr_reals = NULL;// *lp_soln = NULL;
	int *arr_ints = NULL;
	bool selfloop;
	int arr_size, count, res;
	// misc
	int i, j, k, k_r, l1, h1, l2, h2, l2_r, h2_r, start_index;
	double d1, d2, kb, kbt;
	bool done, first;
	jclass vn_cls;
	jmethodID vn_mid;
	int n_maybe_r = 0; // the number of maybe states for rewards

	double lp_result = 0.0;

	/********************** Multiple rewards ****************************
	 // Variables for rewards
	 jlong *rewards_ptrs = env->GetLongArrayElements(_rewards, 0);;
	 DdNode **rewards;	// transition rewards
	 int num_rewards = (int)env->GetArrayLength(_rewards);
	 rewards = new DdNode*[num_rewards];
	 for (i = 0; i < num_rewards; i++) 
	 rewards[i] = jlong_to_DdNode(rewards_ptrs[i]);

	 // Variables for reward formulas
	 jlong *targets_r_ptrs = env->GetLongArrayElements(_targets_r, 0);;
	 DdNode **targets_r;	// target sets for reward formulas
	 int num_targets_r = (int)env->GetArrayLength(_targets_r);
	 targets_r = new DdNode*[num_targets_r];
	 double *bounds_r = env->GetDoubleArrayElements(_bounds_r, 0);
	 for (i = 0; i < num_targets_r; i++) 
	 targets_r[i] = jlong_to_DdNode(targets_r_ptrs[i]);
	 DdNode* yes_r[num_targets_r];
	 double* yes_vecs_r[num_targets_r];
	 for(i=0; i<num_targets_r; i++) 
	 yes_r[i] = targets_r[i];
	 ********************************************************************/

	// extract arrays of target info from function parameters
	num_targets = (int)env->GetArrayLength(_targets);
	target_ptrs = env->GetLongArrayElements(_targets, 0);
	targets = new DdNode*[num_targets];
	for (i = 0; i < num_targets; i++) targets[i] = jlong_to_DdNode(target_ptrs[i]);
	relops = env->GetIntArrayElements(_relops, 0);
	bounds = env->GetDoubleArrayElements(_bounds, 0);

	// Display some target info (just a test)
	PS_PrintToMainLog(env, "\n%d Targets:\n", num_targets);
	for (i = 0; i < num_targets; i++) {
		PS_PrintToMainLog(env, "#%d: ", i);
		switch (relops[i]) {
			case 0: PS_PrintToMainLog(env, "Pmax=?"); break;
			case 1: PS_PrintToMainLog(env, ">%g", bounds[i]); break;
			case 2: PS_PrintToMainLog(env, ">=%g", bounds[i]); break;
			case 3: PS_PrintToMainLog(env, "Rmax=?"); break;
			case 4: PS_PrintToMainLog(env, ">=%g", bounds[i]); break;
		}
		PS_PrintToMainLog(env, " (%.0f states)\n", DD_GetNumMinterms(ddman, targets[i], num_rvars));
	}

	//DdNode *yes1 = targets[0];
	//DdNode *yes2 = targets[1];
	DdNode* yes[num_targets];
	double* yes_vecs[num_targets];
	for(i=0; i<num_targets; i++)
	yes[i] = targets[i];

	// exception handling around whole function
	try {

		// start clocks	
		start1 = start2 = util_cpu_time();

		// get a - filter out rows
		Cudd_Ref(trans);
		Cudd_Ref(maybe);
		if(checkReach) // For reachability, we only use maybe
		a = DD_Apply(ddman, APPLY_TIMES, trans, maybe);
		else { // For LTL formulas, we need both maybe and yes
			DdNode *maybe_yes = maybe;
			for(i=0; i<num_targets; i++) {
				Cudd_Ref(yes[i]);
				maybe_yes = DD_Or(ddman, maybe_yes, yes[i]);
			}
			a = DD_Apply(ddman, APPLY_TIMES, trans, maybe_yes);
		}

		// get number of states
		n = odd->eoff + odd->toff;

		// build sparse matrix
		PS_PrintToMainLog(env, "\nBuilding sparse matrix... ");
		ndsm = build_nd_sparse_matrix(ddman, a, rvars, cvars, num_rvars, ndvars, num_ndvars, odd);
		// get number of transitions/choices
		nnz = ndsm->nnz;
		nc = ndsm->nc;
		kb = ndsm->mem;
		kbt = kb;
		// print out info
		PS_PrintToMainLog(env, "[n=%d, nc=%d, nnz=%d, k=%d] ", n, nc, nnz, ndsm->k);
		PS_PrintMemoryToMainLog(env, "[", kb, "]\n");

		// get vectors for yes/maybe
		PS_PrintToMainLog(env, "Creating vectors for yes... ");
		for(i=0; i<num_targets; i++) {
			yes_vecs[i] = mtbdd_to_double_vector(ddman, yes[i], rvars, num_rvars, odd);
			kb = n*8.0/1024.0;
			kbt += kb;
			PS_PrintMemoryToMainLog(env, "[", kb, "]\n");
		}
		/*PS_PrintToMainLog(env, "Creating vectors for yes2... ");
		 yes2_vec = mtbdd_to_double_vector(ddman, yes2, rvars, num_rvars, odd);
		 kb = n*8.0/1024.0;
		 kbt += kb;
		 PS_PrintMemoryToMainLog(env, "[", kb, "]\n");*/
		PS_PrintToMainLog(env, "Creating vector for maybe... ");
		maybe_vec = mtbdd_to_double_vector(ddman, maybe, rvars, num_rvars, odd);
		kb = n*8.0/1024.0;
		kbt += kb;
		PS_PrintMemoryToMainLog(env, "[", kb, "]\n");

		// get index of single (first) initial state
		start_index = get_index_of_first_from_bdd(ddman, start, rvars, num_rvars, odd);
		PS_PrintToMainLog(env, "Initial state index: %1d\n", start_index);

		// print total memory usage
		PS_PrintMemoryToMainLog(env, "TOTAL: [", kbt, "]\n");

		// store local copies of sparse matrix stuff
		double *non_zeros = ndsm->non_zeros;
		unsigned char *row_counts = ndsm->row_counts;
		int *row_starts = (int *)ndsm->row_counts;
		unsigned char *choice_counts = ndsm->choice_counts;
		int *choice_starts = (int *)ndsm->choice_counts;
		bool use_counts = ndsm->use_counts;
		unsigned int *cols = ndsm->cols;

		/***************** set up rewards *******************/
		/*double *non_zeros_r;
		 unsigned char *choice_counts_r;
		 int *choice_starts_r;
		 bool use_counts_r;
		 unsigned int *cols_r;
		 
		 if(relops[0] > 2) {
		 trans_rewards	= jlong_to_DdNode(trr);	// transition rewards		
		 Cudd_Ref(trans_rewards);
		 trans_rewards = DD_Apply(ddman, APPLY_TIMES, trans_rewards, maybe);
		 PS_PrintToMainLog(env, "Building sparse matrix (transition rewards)... ");
		 ndsm_r = build_sub_nd_sparse_matrix(ddman, a, trans_rewards, rvars, cvars, num_rvars, ndvars, num_ndvars, odd);
		 nnz_r = ndsm_r->nnz;
		 nc_r = ndsm_r->nc;
		 // print out info
		 PS_PrintToMainLog(env, "[n=%d, nc=%d, nnz=%d, k=%d] ", n, nc_r, nnz_r, ndsm_r->k);
		 kb = (nnz_r*12.0+nc_r*4.0+n*4.0)/1024.0;
		 kbt += kb;
		 PS_PrintMemoryToMainLog(env, "[", kb, "]\n");
		 non_zeros_r = ndsm_r->non_zeros;
		 choice_counts_r = ndsm_r->choice_counts;
		 choice_starts_r = (int *)ndsm_r->choice_counts;
		 use_counts_r = ndsm_r->use_counts;
		 cols_r = ndsm_r->cols;
		 //for(i=0; i<nnz_r; i++)
		 //	printf("non_zeros_r[%1d] = %g\n", i, non_zeros_r[i]);
		 / *PS_PrintToMainLog(env, "Creating vector for reward maybe... ");
		 maybe_vec_r = mtbdd_to_double_vector(ddman, maybe_r, rvars, num_rvars, odd);
		 for(i=0; i<n; i++)
		 if(maybe_vec_r[i] > 0)
		 n_maybe_r ++;
		 kb = n*8.0/1024.0;
		 kbt += kb;
		 PS_PrintMemoryToMainLog(env, "[", kb, "]\n");	* /
		 //maybe_vec_r = maybe_vec;
		 }*/
		/***************** end of set up *******************/

		// Set up LP problem
		PS_PrintToMainLog(env, "\nBuilding LP problem...\n");
		// Compute number of vars/constraints
		int maybe_nc = 0;
		for (i = 0; i < n; i++) {
			if (maybe_vec[i]> 0) {
				h1 = use_counts ? row_counts[i] : (row_starts[i+1] - row_starts[i]);
				maybe_nc += h1;
			}
		}
		int num_yes12 = 0;
		int num_maybe = 0;
		for (i = 0; i < n; i++) {
			//if (yes1_vec[i] > 0 || yes2_vec[i] > 0) num_yes12++;
			for(j=0; j<num_targets; j++)
			if(yes_vecs[j][i]> 0) {
				num_yes12++;
				break;
			}
			if (maybe_vec[i]> 0 ) num_maybe++;
		}
		int num_lp_vars = maybe_nc + num_yes12;
		//printf("num_lp_vars = %d + %d = %d\n", maybe_nc, num_yes12, num_lp_vars);
		//int num_constr = num_maybe + num_yes12;
		//printf("num_constr = %d + %d = %d\n", num_maybe, num_yes12, num_constr);

		/*yes1_vec = yes_vecs[0];
		 yes2_vec = yes_vecs[1];
		 
		 start3 = util_cpu_time();

		 if ((lp=make_lp(n, 0)) == NULL) throw "Could not create LP problem";
		 // Lower verbosity: Warnings and errors only
		 //set_verbose(lp, IMPORTANT);

		 // Compute array size: max num elements in a row + 1
		 // or number of elements in maybe, which ever is greater
		 arr_size = 0;
		 h1 = h2 = 0;
		 for (i = 0; i < n; i++) {
		 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++) {
		 h2 = use_counts ? choice_counts[j] : (choice_starts[j+1] - choice_starts[j]);
		 if (h2 > arr_size) arr_size = h2;
		 }
		 }
		 arr_size++;
		 count = 0;
		 for (i = 0; i < n; i++) if (maybe_vec[i] > 0) count++;
		 if (count > arr_size) arr_size = count;
		 // Allocate vectors for use in building LP problem
		 arr_reals = new REAL[arr_size];
		 arr_ints = new int[arr_size];
		 
		 // Set objective function for LP problem
		 set_maxim(lp); 
		 count = 0;
		 / *	for (i = 0; i < n; i++) if (maybe_vec[i] > 0) {
		 arr_reals[count] = 1.0;
		 arr_ints[count] = i+1;
		 count++;
		 }
		 set_obj_fnex(lp, count, arr_reals, arr_ints);* /
		 
		 // Add columns to LP problem
		 h1 = h2 = 0;
		 for (i = 0; i < n; i++) {
		 printf("%d: %s%s%s\n", i, maybe_vec[i]>0?"1":"0", yes1_vec[i]>0?"1":"0", yes2_vec[i]>0?"1":"0");
		 if (yes1_vec[i] > 0 || yes2_vec[i] > 0) {
		 printf("yes %d\n", i);
		 arr_ints[0] = i+1;
		 arr_reals[0] = 1;
		 count = 1;
		 if (yes2_vec[i] > 0) {
		 
		 count++;
		 }
		 add_columnex(lp, 1, arr_reals, arr_ints);
		 }
		 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++) {
		 if (!use_counts) { l2 = choice_starts[j]; h2 = choice_starts[j+1]; }
		 else { l2 = h2; h2 += choice_counts[j]; }
		 if (maybe_vec[i] > 0) {
		 arr_ints[0] = i+1;
		 arr_reals[0] = 1.0;
		 count = 1;
		 d2 = 0.0;
		 for (k = l2; k < h2; k++) {
		 printf("#%d (%d.%d): %d,%f\n", count, i, j, cols[k], non_zeros[k]);
		 // Ignore no cols
		 printf("%s%s%s\n", maybe_vec[cols[k]]>0?"1":"0", yes1_vec[cols[k]]>0?"1":"0", yes2_vec[cols[k]]>0?"1":"0");
		 if (maybe_vec[cols[k]] == 0 && yes1_vec[cols[k]] == 0 && yes2_vec[cols[k]] == 0) continue;
		 // Normal case
		 if (1) {
		 arr_ints[count] = cols[k]+1;
		 arr_reals[count] = -non_zeros[k];
		 count++;
		 }
		 }
		 for (k = 0; k < count; k++) printf(" %f", arr_reals[k]); printf("\n");
		 for (k = 0; k < count; k++) printf(" %d", arr_ints[k]); printf("\n");
		 add_columnex(lp, count, arr_reals, arr_ints);
		 }
		 }
		 }
		 
		 print_lp(lp);
		 
		 // get setup time
		 stop = util_cpu_time();
		 time_for_setup = (double)(stop - start3)/1000;
		 start2 = stop;

		 // Solve LP problem
		 PS_PrintToMainLog(env, "Solving LP problem...\n");
		 res = solve(lp);
		 
		 //Get LP solving time
		 stop2 = util_cpu_time();  
		 time_for_lp = (double)(stop2 - start2)/1000;
		 
		 if (res != 0) throw "Could not solve LP problem";
		 get_ptr_variables(lp, &lp_soln);
		 //print_solution(lp, 2);
		 
		 // Reuse yes vector to store/return soln
		 for (i = 0; i < n; i++) {
		 //if (yes2_vec[i] == 0.0) yes2_vec[i] = lp_soln[i];
		 }
		 
		 // stop clocks
		 stop = util_cpu_time();
		 time_for_lp = (double)(stop - start2)/1000;
		 time_taken = (double)(stop - start3)/1000;
		 
		 // print timing info
		 PS_PrintToMainLog(env, "\nLP problem solved in %.2f seconds (setup %.2f, lpsolve %.2f)\n", time_taken, time_for_setup, time_for_lp);*/

		/*
		 * Hongyang: Use Eclipse with COIN CBC solver 
		 */
		FILE *f; /* FILE pointer */

		// Generate LP for the objectives
		start3 = util_cpu_time();

		ec_refs Vars;
		ec_ref Cost;
		pword varlist;
		//int yes_vec[n];
		//int map_var[n+1]; 
		int *yes_vec;
		int *map_var;

		yes_vec = new int[n];
		map_var = new int[n+1];

		arr_ints = new int[n];
		int sp;

		int yes_nc = 0; // the total number of choices of all target states
		/*int *map_var_r;

		 if(relops[0] > 2) {
		 map_var_r = new int[n];
		 k = 0;
		 for(i=0; i<n; i++) 
		 if(maybe_vec_r[i] > 0)
		 map_var_r[i] = k++;
		 else
		 map_var_r[i] = -1;
		 }*/

		if(checkReach) {
			PS_PrintToMainLog(env, "Check reachability...\n");
			// Maps states to ECLiPSe variables. 
			// Each yes state has a coresponding variable.
			// Each maybe state has one or more variables.

			// merge all yes vectors into one
			int flag = 0;
			for(i=0; i<n; i++) {
				flag = 0;
				for(j=0; j<num_targets; j++)
				if(yes_vecs[j][i]> 0) {
					yes_vec[i] = 1;
					flag = 1;
					break;
				}
				if(!flag)
				yes_vec[i] = 0;
			}

			//if (arr_reals) delete[] arr_reals;
			PS_PrintToMainLog(env, "Number of LP variables = %1d\n", num_lp_vars + n_maybe_r);
			arr_reals = new REAL[num_lp_vars + n_maybe_r];

			// create map vector
			count = 0;
			for(i=0; i<n; i++) {
				if(yes_vec[i]> 0)
				map_var[i] = count++;
				else {
					if(maybe_vec[i]> 0) {
						map_var[i] = count;
						count += use_counts ? row_counts[i] : (row_starts[i+1] - row_starts[i]);
					} else
					map_var[i] = count;
				}
			}
			map_var[n] = num_lp_vars;

			/*printf("map_var = ");
			 for(i=0; i<=n; i++)
			 printf(" %1d", map_var[i]);
			 printf("\n");
			 fflush(stdout);*/

			// Initialize ECLiPSe
			ec_init();

			// Load Eplex library
			ec_post_string("lib(eplex)");

			// Create ECLiPSe variables
			Vars = ec_refs_create_newvars(num_lp_vars + n_maybe_r);
			Cost = ec_ref_create_newvar();
			varlist = ec_listofrefs(Vars);

			// Define the domain of the variables
			ec_post_goal(ec_term(
							ec_did("::",2),
							varlist,
							ec_term(ec_did("..",2), ec_double(0.0), ec_double(1e20))));

			stop = util_cpu_time();
			time_for_setup = (double)(stop - start3)/1000;
			start2 = stop;

			// Add constraints to LP problem
			int x;
			l1 = h1 = 0;
			for (i = 0; i < n; i++) {
				// Add constraints for maybe states
				if (maybe_vec[i]> 0) {
					for(k=0; k<num_lp_vars; k++)
					arr_reals[k] = 0.0;

					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++) // \Sigma_{v\in V/F} y_{(v, \gamma)}
					arr_reals[map_var[i]+j-l1] = 1.0;

					count = h1-l1;
					int l3, h3;
					h3 = h2 = 0;
					for (x = 0; x < n; x++) { // h2: row index
						if(maybe_vec[x]> 0) {
							if (!use_counts) {
								l3 = row_starts[x];
								h3 = row_starts[x+1];
							} else { // Attention: there could exist a bug here!!!!!!!
								l3 = h3;
								h3 += row_counts[x];
							}
							for (j = l3; j < h3; j++) { // j: choice index; x+j: variable index 
								if (!use_counts) {
									l2 = choice_starts[j];
									h2 = choice_starts[j+1];
								} else {
									l2 = h2;
									h2 += choice_counts[j];
								}
								for(k=l2; k<h2; k++) { // k: column index
									if(cols[k] == i) {
										// get the index of the corresponding variable
										double old_value = arr_reals[map_var[x]+j-l3];
										arr_reals[map_var[x]+j-l3] -= non_zeros[k];
										//printf("map_var[%1d]+%1d-%1d = %1d\n", x, j, l3, map_var[x]+j-l3);
										if(old_value == 0.0 && arr_reals[map_var[x]+j-l3] != 0.0)
										count++;
										else if(old_value != 0.0 && arr_reals[map_var[x]+j-l3] == 0.0)
										count--;
									}
								}
							}
						} else if(use_counts) {
							l3 = h3;
							h3 += row_counts[x];
							for (j = l3; j < h3; j++)
							h2 += choice_counts[j];
						}
					}

					/*printf("		arr_reals:");
					 for (k = 0; k < num_lp_vars; k++) 
					 printf(" %1.1f", arr_reals[k]); 
					 printf(" = %1.1f", (start_index==i ? 1.0 : 0.0));
					 printf("		count = %1d\n", count);
					 fflush(stdout);*/

					pword tail = ec_nil();
					x = 0;

					for(k=num_lp_vars-1; k>=0; k--)
					if(arr_reals[k] != 0.0)
					tail = ec_list(ec_refs_get(Vars, k), tail);
					for(k=0; k<num_lp_vars; k++)
					if(arr_reals[k] != 0.0)
					arr_reals[x++] = arr_reals[k];
					ec_post_goal(ec_term(ec_did(("$="),2),
									ec_term(ec_did("*",2),
											tail,
											ec_listofdouble(count, arr_reals)),
									ec_double((start_index==i ? 1.0 : 0.0))));
					/*ec_post_goal(ec_term(ec_did(("$="),2), 
					 ec_term(ec_did("*",2), 
					 varlist, 
					 ec_listofdouble(num_lp_vars, arr_reals)), 
					 ec_double((start_index==i ? 1.0 : 0.0))));*/
				} else if(use_counts) {
					l1 = h1;
					h1 += row_counts[i];
				}
			}
			for (i = 0; i < n; i++) {
				if(yes_vec[i]> 0) {
					for(k=0; k<num_lp_vars; k++)
					arr_reals[k] = 0.0;

					count = 1;
					h1 = h2 = 0;
					for (x = 0; x < n; x++) { // x: row index
						if(maybe_vec[x]> 0) {
							if (!use_counts) {
								l1 = row_starts[x];
								h1 = row_starts[x+1];
							} else {
								l1 = h1;
								h1 += row_counts[x];
							}
							if(x == i) {
								if(use_counts)
								for (j = l1; j < h1; j++)
								h2 += choice_counts[j];
								continue;
							}
							for (j = l1; j < h1; j++) { // j: choice index; x+j: variable index 
								if (!use_counts) {
									l2 = choice_starts[j];
									h2 = choice_starts[j+1];
								} else {
									l2 = h2;
									h2 += choice_counts[j];
								}
								for(k=l2; k<h2; k++) { // k: column index
									if(cols[k] == i) {
										// get the index of the corresponding variable
										arr_reals[map_var[x]+j-l1] = -non_zeros[k];
										//printf("map_var[%1d]+%1d-%1d = %1d\n", x, j, l1, map_var[x]+j-l1);
										//printf("arr_reals[%1d] = %1.1f\n", map_var[x]+j-l1, -non_zeros[k]);
										count++;
									}
								}
							}
						} else if(use_counts) {
							l1 = h1;
							h1 += row_counts[x];
							for (j = l1; j < h1; j++)
							h2 += choice_counts[j];
						}
					}

					arr_reals[map_var[i]] = 1.0;
					/*printf("				arr_reals:");
					 for (k = 0; k < num_lp_vars; k++) 
					 printf(" %1.1f", arr_reals[k]);
					 printf(" = 0.0");
					 printf("		count = %1d\n", count);
					 fflush(stdout);*/

					pword tail = ec_nil();
					x = 0;

					for(k=num_lp_vars-1; k>=0; k--)
					if(arr_reals[k] != 0.0)
					tail = ec_list(ec_refs_get(Vars, k), tail);
					for(k=0; k<num_lp_vars; k++)
					if(arr_reals[k] != 0.0)
					arr_reals[x++] = arr_reals[k];
					ec_post_goal(ec_term(ec_did(("$="),2),
									ec_term(ec_did("*",2),
											tail,
											ec_listofdouble(count, arr_reals)),
									ec_double(0.0)));

					/*ec_post_goal(ec_term(ec_did(("$="),2), 
					 ec_term(ec_did("*",2), 
					 varlist, 
					 ec_listofdouble(num_lp_vars, arr_reals)), 
					 ec_double(0.0)));*/
				}
			}

			PS_PrintToMainLog(env, "Adding extra constraints for probabilisitic bounds\n");
			for(i=0; i<num_targets; i++) {
				pword tail = ec_nil();
				for(k=n-1; k>=0; k--)
				if(yes_vecs[i][k]> 0)
				tail = ec_list(ec_refs_get(Vars, map_var[k]), tail);
				x=0;
				for(k=0; k<n; k++)
				if(yes_vecs[i][k]> 0)
				arr_reals[x++] = 1.0;
				ec_post_goal(ec_term(ec_did("$>=",2),
								ec_term(ec_did("*",2),
										tail,
										ec_listofdouble(x, arr_reals)),
								ec_double(bounds[i])));
			}

			PS_PrintToMainLog(env, "Computing optimisation...\n");
			fflush(stdout);
			pword tail = ec_nil();
			x = 0;

			if(relops[0]> 0 && relops[0] <= 2) {
				for(k=n-1; k>=0; k--)
				if(yes_vec[k]> 0)
				tail = ec_list(ec_refs_get(Vars, map_var[k]), tail);
				for(k=0; k<n; k++)
				if(yes_vec[k]> 0)
				arr_reals[x++] = 1.0;
			} else if(relops[0] == 0) {
				for(k=n-1; k>=0; k--)
				if(yes_vecs[0][k]> 0)
				tail = ec_list(ec_refs_get(Vars, map_var[k]), tail);
				for(k=0; k<n; k++)
				if(yes_vecs[0][k]> 0)
				arr_reals[x++] = 1.0;
			} else {
				/*x = 1;
				 arr_reals[0] = 1.0;
				 tail = ec_list(ec_refs_get(Vars, map_var_r[start_index]+num_lp_vars), tail);*/
				PS_PrintToMainLog(env, "computing max function for rewards ...\n");
				double *non_zeros_r;
				unsigned char *choice_counts_r;
				int *choice_starts_r;
				bool use_counts_r;
				unsigned int *cols_r;
				trans_rewards = jlong_to_DdNode(trr); // transition rewards		
				Cudd_Ref(trans_rewards);
				trans_rewards = DD_Apply(ddman, APPLY_TIMES, trans_rewards, maybe);
				PS_PrintToMainLog(env, "Building sparse matrix (transition rewards)... ");
				ndsm_r = build_sub_nd_sparse_matrix(ddman, a, trans_rewards, rvars, cvars, num_rvars, ndvars, num_ndvars, odd);
				nnz_r = ndsm_r->nnz;
				nc_r = ndsm_r->nc;
				// print out info
				PS_PrintToMainLog(env, "[n=%d, nc=%d, nnz=%d, k=%d] ", n, nc_r, nnz_r, ndsm_r->k);
				kb = (nnz_r*12.0+nc_r*4.0+n*4.0)/1024.0;
				kbt += kb;
				PS_PrintMemoryToMainLog(env, "[", kb, "]\n");
				non_zeros_r = ndsm_r->non_zeros;
				choice_counts_r = ndsm_r->choice_counts;
				choice_starts_r = (int *)ndsm_r->choice_counts;
				use_counts_r = ndsm_r->use_counts;
				cols_r = ndsm_r->cols;

				for(i=0; i<num_lp_vars; i++)
				arr_reals[i] = 0;

				h1 = h2 = 0;
				for(i=0; i<n; i++)
				if(maybe_vec[i]> 0) {
					if (!use_counts) { // assume we always do not use counts
						l1 = row_starts[i];
						h1 = row_starts[i+1];
					} else { // Problematic
						l1 = h1;
						h1 += row_counts[i];
					}
					for (j = l1; j < h1; j++) {// \Sigma_{v\in V/F} y_{(v, \gamma)}*c_{(v, \gamma)}
						if (!use_counts) {
							l2 = choice_starts[j];
							h2 = choice_starts[j+1];
						} else {
							l2 = h2;
							h2 += choice_counts[j];
						}
						if (!use_counts_r) {
							l2_r = choice_starts_r[j];
							h2_r = choice_starts_r[j+1];
						} else {
							l2_r = h2_r;
							h2_r += choice_counts_r[j];
						}
						k_r = l2_r;
						while (k_r < h2_r && cols_r[k_r] != cols[l2])
						k_r++;
						// if there is one, add reward * prob to reward value
						if (k_r < h2_r) {
							arr_reals[map_var[i]+j-l1] = non_zeros_r[k_r];
						}
					}
				} else if(use_counts) { // count the number of choices for non maybe nodes. Is it necessary?
					l1 = h1;
					h1 += row_counts[i];
					for (j = l1; j < h1; j++) {
						h2 += choice_counts[j];
						h2_r += choice_counts_r[j];
					}
				}
				/*PS_PrintToMainLog(env, "		arr_reals:");
				 for (k = 0; k < num_lp_vars; k++) 
				 PS_PrintToMainLog(env, " %1.1f", arr_reals[k]); 
				 fflush(stdout);*/

				for(k=num_lp_vars-1; k>=0; k--)
				if(arr_reals[k] != 0.0)
				tail = ec_list(ec_refs_get(Vars, k), tail);
				for(k=0; k<num_lp_vars; k++)
				if(arr_reals[k] != 0.0)
				arr_reals[x++] = arr_reals[k];

			}

			ec_post_goal( ec_term(ec_did("optimize", 2),
							ec_term(ec_did("max", 1),
									ec_term(ec_did("*",2), tail, ec_listofdouble(x, arr_reals))),
							ec_ref_get(Cost)));
			/*ec_post_goal( ec_term(ec_did("optimize", 2),
			 ec_term(ec_did("max", 1), 
			 ec_term(ec_did("*",2), varlist, ec_listofdouble(num_lp_vars, pc))), 
			 ec_ref_get(Cost)));*/
		} else { // ------------------------------ LTL model checking ----------------------------
			PS_PrintToMainLog(env, "Check LTL formulas...\n");

			int num_map;
			int yes_count = 0; // the number of target states

			// Compute the number of lp variables for LTL model checking
			maybe_nc = 0; // number of variables for maybe states
			h1 = h2 = 0;
			count = 0;
			for (i = 0; i < n; i++) {
				if (!use_counts) {
					l1 = row_starts[i];
					h1 = row_starts[i+1];
				} else {
					l1 = h1;
					h1 += row_counts[i];
				}
				k=0;
				for (j = l1; j < h1; j++) { // j: choice index; x+j: variable index 
					if (!use_counts) {
						l2 = choice_starts[j];
						h2 = choice_starts[j+1];
					} else {
						l2 = h2;
						h2 += choice_counts[j];
					}
					if(h2-l2>1 || cols[l2] != i || non_zeros[l2]!=1.0)
					k++;
				}
				if (maybe_vec[i]> 0) {
					maybe_nc += k;
					map_var[i] = count;
					count += k;
				}
				else {
					map_var[i] = count;
					for(j=0; j<num_targets; j++)
					if(yes_vecs[j][i]> 0) {
						yes_nc += k;
						map_var[i] = count;
						yes_count ++; // each target state has one extra action
						count += k+1;
						break;
					}
				}
			}
			num_lp_vars = maybe_nc + yes_nc + yes_count;
			map_var[n] = num_lp_vars; // maybe need to be modified.		
			PS_PrintToMainLog(env, "Number of LP variables = %1d\n", num_lp_vars);

			for(i=0; i<n; i++)
			yes_vec[i] = 0;
			count = 0;
			for(i=0; i<n; i++)
			for(j=0; j<num_targets; j++)
			if(yes_vecs[j][i]> 0) { // For each state, count the number of targets it belongs to.
				yes_vec[i]++;
				count++;
			}

			/*printf("map_var = ");
			 for(i=0; i<=n; i++)
			 printf(" %1d", map_var[i]);
			 printf("\n");
			 fflush(stdout);*/

			h1 = h2 = 0;
			for(i=0; i<n; i++)
			arr_ints[i] = map_var[i+1]-map_var[i];
			for (i = 0; i < n; i++) {
				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++) { // j: choice index; x+j: variable index 
					if (!use_counts) {
						l2 = choice_starts[j];
						h2 = choice_starts[j+1];
					} else {
						l2 = h2;
						h2 += choice_counts[j];
					}
					if(h2-l2>1 || cols[l2] != i || non_zeros[l2]!=1.0) {
						if (maybe_vec[i]> 0 || yes_vec[i]> 0)
						for(k=l2; k<h2; k++) {
							if(maybe_vec[cols[k]]> 0 || yes_vec[cols[k]]> 0)
							arr_ints[cols[k]] ++;
						}
					}
				}
			}
			//int *constraints_ints[n];
			//double *constraints_reals[n];
			int *constraints_sp;
			//int constraints_sp[n];
			int *constraints_ints;
			double *constraints_reals;
			unsigned long *constraints_int_ptr;
			unsigned long *constraints_real_ptr;
			//PS_PrintToMainLog(env, "Allocating memory space		\n");
			constraints_sp = new int[n];
			constraints_int_ptr = new unsigned long[n];
			constraints_real_ptr = new unsigned long[n];
			for(i=0; i<n; i++)
			if(maybe_vec[i]> 0 || yes_vec[i]> 0) {
				//printf("i = %1d, ", i);
				constraints_ints = new int[arr_ints[i]];
				constraints_reals = new REAL[arr_ints[i]];
				constraints_sp[i] = map_var[i+1]-map_var[i];
				constraints_int_ptr[i] = (unsigned long)constraints_ints;
				constraints_real_ptr[i] = (unsigned long)constraints_reals;
				for(j=0; j<map_var[i+1]-map_var[i]; j++) {
					constraints_ints[j] = map_var[i] + j;
					constraints_reals[j] = 1.0;
				}
				for(j=map_var[i+1]-map_var[i]; j<arr_ints[i]; j++) {
					constraints_ints[j] = -1;
					constraints_reals[j] = 0;
				}

				//printf("arr_ints[%1d] = %1d\n", i, arr_ints[i]);
			}

			//if (arr_reals) delete[] arr_reals;
			arr_reals = new REAL[num_lp_vars];

			// create map vector
			/*count = 0;
			 for(i=0; i<n; i++) {
			 if(yes_vec[i] > 0 || maybe_vec[i] > 0) {
			 map_var[i] = count;	
			 count += use_counts ? row_counts[i] : (row_starts[i+1] - row_starts[i]);
			 if(yes_vec[i] > 0)
			 count ++; // Add extra transition
			 } else 
			 map_var[i] = count;
			 }
			 map_var[n] = num_lp_vars; // maybe need to be modified.*/

			/* write the LP to file LP.ecl */
			/*f = fopen("LP.ecl", "w"); /* create a file for writing * /
			 
			 if(f==NULL) 
			 PS_PrintToMainLog(env, "Cannot create LP.ecl!\n");
			 
			 fprintf(f, ":- lib(eplex).\n\n");
			 fprintf(f, "main1(Cost, Vars) :-\n");
			 fprintf(f, "Vars = [");
			 for(i=0; i<num_lp_vars; i++) {
			 if(i>0)
			 fprintf(f, ", ");
			 fprintf(f, "X%1d", i+1);
			 }
			 
			 fprintf(f, "],\n");
			 fprintf(f, "Vars :: 0.0..inf,\n\n");*/

			/*printf("map_var = ");
			 for(i=0; i<=n; i++)
			 printf(" %1d", map_var[i]);
			 printf("\n");
			 fflush(stdout);*/

			// Initialize ECLiPSe
			ec_init();

			// Load Eplex library
			ec_post_string("lib(eplex)");

			// Create ECLiPSe variables
			Vars = ec_refs_create_newvars(num_lp_vars);
			Cost = ec_ref_create_newvar();
			varlist = ec_listofrefs(Vars);

			// Define the domain of the variables
			ec_post_goal(ec_term(
							ec_did("::",2),
							varlist,
							ec_term(ec_did("..",2), ec_double(0.0), ec_double(1e20))));

			stop = util_cpu_time();
			time_for_setup = (double)(stop - start3)/1000;
			start2 = stop;

			// Add constraints to LP problem
			int x;
			//h1 = h2 = 0;
			/*for (i = 0; i < n; i++) {
			 // Add constraints for maybe and target states
			 if (maybe_vec[i] > 0 || yes_vec[i] > 0) 
			 for (j = 0; j < map_var[i+1]-map_var[i]; j++) {// \Sigma_{v\in V/F} y_{(v, \gamma)}
			 printf("constraints_ints[%1d][%1d] = %1d, ", i, j, constraints_ints[i][j]);
			 printf("constraints_reals[%1d][%1d] = %g\n", i, j, constraints_reals[i][j]);
			 }
			 printf("constraints_sp[%1d] = %1d\n\n", i, constraints_sp[i]);
			 }*/

			//count = map_var[i+1] - map_var[i];
			h1 = h2 = 0;
			for (x = 0; x < n; x++) { // h2: row index
				if(maybe_vec[x]> 0 || yes_vec[x]> 0) {
					if (!use_counts) {
						l1 = row_starts[x];
						h1 = row_starts[x+1];
					} else {
						l1 = h1;
						h1 += row_counts[x];
					}
					count = 0;
					for (j = l1; j < h1; j++) { // j: choice index; x+j: variable index 
						if (!use_counts) {
							l2 = choice_starts[j];
							h2 = choice_starts[j+1];
						} else {
							l2 = h2;
							h2 += choice_counts[j];
						}
						if(h2-l2==1 && cols[l2] == x)
						count++;
						else
						for(k=l2; k<h2; k++) { // k: column index
							// get the index of the corresponding variable
							i = cols[k];
							if(maybe_vec[i]> 0 || yes_vec[i]> 0) {
								//printf("	x = %1d, i = %1d, nonzeros = %g, j = %1d, l1 = %1d, count = %1d, map_var[x]+j-l1-count = %1d\n", 
								//		 x, i, non_zeros[k], j, l1, count, map_var[x]+j-l1-count);
								constraints_ints = (int *)constraints_int_ptr[i];
								constraints_reals = (double *)constraints_real_ptr[i];
								for(k_r=0; k_r<constraints_sp[i]; k_r++)
								if(constraints_ints[k_r]==map_var[x]+j-l1-count)
								break;

								if(k_r == constraints_sp[i]) {
									//count++;
									constraints_ints[k_r]=map_var[x]+j-l1-count;
									constraints_reals[k_r] = -non_zeros[k];
									constraints_sp[i]++;
								} else {
									constraints_reals[k_r] -= non_zeros[k];
									if(constraints_reals[k_r] == 0.0) { // In fact, this cannot happen.
										//count--;
										for(l2_r=k_r; l2_r<constraints_sp[i]-1; l2_r++) {
											constraints_reals[l2_r] = constraints_reals[l2_r+1];
											constraints_ints[l2_r] = constraints_ints[l2_r+1];
										}
										constraints_sp[i]--;
									}
								}
							}
						}
					}
				} else if(use_counts) {
					l1 = h1;
					h1 += row_counts[x];
					for (j = l1; j < h1; j++)
					h2 += choice_counts[j];
				}
			}

			/*fprintf(f, "	");
			 l2 = 0;
			 for (k = 0; k < num_lp_vars; k++) 
			 if(arr_reals[k] != 0) {
			 if(arr_reals[k] > 0 && l2 > 0)
			 fprintf(f, "+");
			 fprintf(f, "%g*X%1d", arr_reals[k], k+1); 
			 l2++;
			 }
			 fprintf(f, " $= %g,\n", (start_index==i ? 1.0 : 0.0));
			 //printf("		count = %1d\n", count);
			 fflush(f);*/

			//pword tail = ec_nil();
			//x = 0;

			/*for(k=num_lp_vars-1; k>=0; k--)
			 if(arr_reals[k] != 0.0) 
			 tail = ec_list(ec_refs_get(Vars, k), tail);
			 for(k=0; k<num_lp_vars; k++)
			 if(arr_reals[k] != 0.0) 
			 arr_reals[x++] = arr_reals[k];*/
			/*printf("arr_reals & arr_ints = ");
			 for(k=0; k<sp; k++) {
			 printf("(%1d, %g) ", arr_ints[k], arr_reals[k]);		
			 }
			 printf("\n");*/

			/*for(k=sp-1; k>=0; k--)
			 tail = ec_list(ec_refs_get(Vars, arr_ints[k]), tail);
			 ec_post_goal(ec_term(ec_did(("$="),2), 
			 ec_term(ec_did("*",2), 
			 tail, 
			 ec_listofdouble(sp, arr_reals)), 
			 ec_double((start_index==i ? 1.0 : 0.0))));*/
			/*ec_post_goal(ec_term(ec_did(("$="),2), 
			 ec_term(ec_did("*",2), 
			 varlist, 
			 ec_listofdouble(num_lp_vars, arr_reals)), 
			 ec_double((start_index==i ? 1.0 : 0.0))));*/
			//printf("output LP constraints to ECLiSPe\n");
			for (i = 0; i < n; i++)
			if (maybe_vec[i]> 0 || yes_vec[i]> 0) {
				constraints_ints = (int *)constraints_int_ptr[i];
				constraints_reals = (double *)constraints_real_ptr[i];
				pword tail = ec_nil();
				for(k=constraints_sp[i]-1; k>=0; k--)
				tail = ec_list(ec_refs_get(Vars, constraints_ints[k]), tail);
				ec_post_goal(ec_term(ec_did(("$="),2),
								ec_term(ec_did("*",2),
										tail,
										ec_listofdouble(constraints_sp[i], constraints_reals)),
								ec_double((start_index==i ? 1.0 : 0.0))));
				/*printf("i =%1d, constrains_sp = %1d, ", i, constraints_sp[i]);
				 for(k=0; k<constraints_sp[i]; k++)
				 printf("(%1d, %g) ", constraints_ints[i][k], constraints_reals[i][k]);
				 printf("\n");*/
				delete[] constraints_ints;
				delete[] constraints_reals;
			}
			delete[] constraints_sp;
			delete[] constraints_int_ptr;
			delete[] constraints_real_ptr;
			PS_PrintToMainLog(env, "Adding extra constraints for probabilisitic bounds\n");
			for(i=0; i<num_targets; i++) {
				if(relops[i]==0 || relops[i]>2)
				continue;
				//PS_PrintToMainLog(env, "num_map = %1d, i = %1d, bound = %1.1f\n", num_map, i, bounds[i]);
				for(k=0; k<num_lp_vars; k++)
				arr_reals[k] = 0.0;

				count = 0; // count the number of extra transitions for the i_th target 
				for(k=0; k<n; k++)
				if(yes_vecs[i][k]> 0) {
					arr_reals[map_var[k+1]-1] = 1.0;
					count++;
				}

				/*fprintf(f, "	");
				 l2 = 0;
				 for (k = 0; k < num_lp_vars; k++) 
				 if(arr_reals[k] != 0) {
				 if(arr_reals[k] > 0 && l2 > 0)
				 fprintf(f, "+");
				 fprintf(f, "%g*X%1d", arr_reals[k], k+1); 
				 l2++;
				 }
				 fprintf(f, " $>= %g,\n", bounds[i]);
				 fflush(f);*/

				pword tail = ec_nil();
				x = 0;

				for(k=num_lp_vars-1; k>=0; k--)
				if(arr_reals[k] != 0.0)
				tail = ec_list(ec_refs_get(Vars, k), tail);
				for(k=0; k<num_lp_vars; k++)
				if(arr_reals[k] != 0.0)
				arr_reals[x++] = arr_reals[k];
				ec_post_goal(ec_term(ec_did(("$>="),2),
								ec_term(ec_did("*",2),
										tail,
										ec_listofdouble(count, arr_reals)),
								ec_double(bounds[i])));

				/*ec_post_goal(ec_term(ec_did(("$="),2), 
				 ec_term(ec_did("*",2), 
				 varlist, 
				 ec_listofdouble(num_lp_vars, arr_reals)), 
				 ec_double(0.0)));*/
			}

			// Adding constraints for rewards
			/*if(relops[0] > 2) {
			 for(i=0; i<num_lp_vars; i++)
			 arr_reals[i] = 0;

			 h1 = h2 = h2_r = 0;
			 for(i=0; i<n; i++) 
			 if(maybe_vec_r[i] > 0) {
			 for(j=num_lp_vars; j<num_lp_vars + n_maybe_r; j++)
			 arr_reals[j] = 0;
			 arr_reals[map_var_r[i]+num_lp_vars] = 1.0;

			 if (!use_counts) { // assume we always do not use counts
			 l1 = row_starts[i]; 
			 h1 = row_starts[i+1]; 
			 } else { // Problematic
			 l1 = h1; 
			 h1 += row_counts[i]; 
			 }
			 for (j = l1; j < h1; j++) {// \Sigma_{v\in V/F} y_{(v, \gamma)}
			 if (!use_counts) { 
			 l2 = choice_starts[j]; 
			 h2 = choice_starts[j+1]; 
			 } else { 
			 l2 = h2; 
			 h2 += choice_counts[j]; 
			 }
			 if (!use_counts_r) { 
			 l2_r = choice_starts_r[j]; 
			 h2_r = choice_starts_r[j+1]; 
			 } else { 
			 l2_r = h2_r; 
			 h2_r += choice_counts_r[j]; 
			 }
			 for (k = l2; k < h2; k++) {
			 // find corresponding transition reward if any
			 // Hongyang: Can it be simplified?
			 k_r = l2_r; 
			 while (k_r < h2_r && cols_r[k_r] != cols[k]) 
			 k_r++;
			 // if there is one, add reward * prob to reward value
			 if (k_r < h2_r && (maybe_vec[cols[k]] > 0 || yes_vecs[0][cols[k]] > 0)) {
			 arr_reals[map_var[i]+j-l1] -= non_zeros_r[k_r] * non_zeros[k]; 
			 if(maybe_vec[cols[k]] > 0)
			 arr_reals[map_var_r[cols[k]]+num_lp_vars] -= non_zeros[k_r];
			 }
			 }
			 
			 }
			 printf("		arr_reals:");
			 for (k = 0; k < num_lp_vars+n_maybe_r; k++) 
			 printf(" %1.1f", arr_reals[k]); 
			 printf(" = 0\n");
			 fflush(stdout);
			 
			 pword tail = ec_nil();
			 x = 0;
			 
			 for(k=num_lp_vars+n_maybe_r-1; k>=0; k--)
			 if(arr_reals[k] != 0.0) 
			 tail = ec_list(ec_refs_get(Vars, k), tail);
			 for(k=0; k<num_lp_vars+n_maybe_r; k++)
			 if(arr_reals[k] != 0.0) 
			 arr_reals[x++] = arr_reals[k];
			 ec_post_goal(ec_term(ec_did(("$="),2), 
			 ec_term(ec_did("*",2), 
			 tail, 
			 ec_listofdouble(x, arr_reals)), 
			 ec_double(0.0)));
			 
			 }
			 }*/

			PS_PrintToMainLog(env, "Computing optimisation...\n");
			fflush(stdout);
			pword tail = ec_nil();
			x = 0;
			if(relops[0]> 0 && relops[0]<=2) {
				//fprintf(f, "	optimize(max(");
				//l2 = 0;
				for(i=0; i<n; i++)
				if(yes_vec[i]> 0) {
					tail = ec_list(ec_refs_get(Vars, map_var[i+1]-1), tail);
					arr_reals[x++] = 1.0;
					//if(l2 > 0)
					//	fprintf(f, "+");				
					//fprintf(f, "X%1d", map_var[i+1]); 
					//l2++;
				}
				//fprintf(f, "), Cost).\n");
				//fflush(f);
			} else if(relops[0] == 0) {
				//fprintf(f, "	optimize(max(");
				//l2 = 0;
				for(i=0; i<n; i++)
				if(yes_vecs[0][i]> 0) {
					tail = ec_list(ec_refs_get(Vars, map_var[i+1]-1), tail);
					arr_reals[x++] = 1.0;
					//if(l2 > 0)
					//	fprintf(f, "+");				
					//fprintf(f, "X%1d", map_var[i+1]); 
					//l2++;
				}
				//fprintf(f, "), Cost).\n");
				//fflush(f);
			} else {
				/*x = 1;
				 arr_reals[0] = 1.0;
				 tail = ec_list(ec_refs_get(Vars, map_var_r[start_index]+num_lp_vars), tail);*/
				PS_PrintToMainLog(env, "computing max function for rewards ...\n");
				double *non_zeros_r;
				unsigned char *choice_counts_r;
				int *choice_starts_r;
				bool use_counts_r;
				unsigned int *cols_r;
				trans_rewards = jlong_to_DdNode(trr); // transition rewards		
				Cudd_Ref(trans_rewards);
				trans_rewards = DD_Apply(ddman, APPLY_TIMES, trans_rewards, maybe);
				PS_PrintToMainLog(env, "Building sparse matrix (transition rewards)... ");
				ndsm_r = build_sub_nd_sparse_matrix(ddman, a, trans_rewards, rvars, cvars, num_rvars, ndvars, num_ndvars, odd);
				nnz_r = ndsm_r->nnz;
				nc_r = ndsm_r->nc;
				// print out info
				PS_PrintToMainLog(env, "[n=%d, nc=%d, nnz=%d, k=%d] ", n, nc_r, nnz_r, ndsm_r->k);
				kb = (nnz_r*12.0+nc_r*4.0+n*4.0)/1024.0;
				kbt += kb;
				PS_PrintMemoryToMainLog(env, "[", kb, "]\n");
				non_zeros_r = ndsm_r->non_zeros;
				choice_counts_r = ndsm_r->choice_counts;
				choice_starts_r = (int *)ndsm_r->choice_counts;
				use_counts_r = ndsm_r->use_counts;
				cols_r = ndsm_r->cols;

				for(i=0; i<num_lp_vars; i++)
				arr_reals[i] = 0;

				h1 = h2 = l2_r = h2_r = 0;
				for(i=0; i<n; i++)
				if(maybe_vec[i]> 0) {
					if (!use_counts) { // assume we always do not use counts
						l1 = row_starts[i];
						h1 = row_starts[i+1];
					} else { // Problematic
						l1 = h1;
						h1 += row_counts[i];
					}
					for (j = l1; j < h1; j++) {// \Sigma_{v\in V/F} y_{(v, \gamma)}*c_{(v, \gamma)}
						if (!use_counts) {
							l2 = choice_starts[j];
							h2 = choice_starts[j+1];
						} else {
							l2 = h2;
							h2 += choice_counts[j];
						}
						if (!use_counts_r) {
							l2_r = choice_starts_r[j];
							h2_r = choice_starts_r[j+1];
						} else {
							l2_r = h2_r;
							h2_r += choice_counts_r[j];
						}
						k_r = l2_r;
						while (k_r < h2_r && cols_r[k_r] != cols[l2])
						k_r++;
						// if there is one, add reward * prob to reward value
						if (k_r < h2_r) {
							arr_reals[map_var[i]+j-l1] = non_zeros_r[k_r];
						}
					}
				} else if(use_counts) { // count the number of choices for non maybe nodes. Is it necessary?
					l1 = h1;
					h1 += row_counts[i];
					for (j = l1; j < h1; j++) {
						h2 += choice_counts[j];
						h2_r += choice_counts_r[j];
					}
				}
				/*PS_PrintToMainLog(env, "		arr_reals:");
				 for (k = 0; k < num_lp_vars; k++) 
				 PS_PrintToMainLog(env, " %1.1f", arr_reals[k]); 
				 PS_PrintToMainLog(env, "\n");
				 fflush(stdout);*/

				for(k=num_lp_vars-1; k>=0; k--)
				if(arr_reals[k] != 0.0)
				tail = ec_list(ec_refs_get(Vars, k), tail);
				for(k=0; k<num_lp_vars; k++)
				if(arr_reals[k] != 0.0)
				arr_reals[x++] = arr_reals[k];
			}
			ec_post_goal( ec_term(ec_did("optimize", 2),
							ec_term(ec_did("max", 1),
									ec_term(ec_did("*",2), tail, ec_listofdouble(x, arr_reals))),
							ec_ref_get(Cost)));
			/*ec_post_goal( ec_term(ec_did("optimize", 2),
			 ec_term(ec_did("max", 1), 
			 ec_term(ec_did("*",2), varlist, ec_listofdouble(num_lp_vars, pc))), 
			 ec_ref_get(Cost)));*/
			//fclose(f);
		}

		/* write the produce mdp to file model.dot */
		printf("Writing the model to model.dot\n"); fflush(stdout);

		f = fopen("model.dot", "w"); /* create a file for writing */
		f = NULL;
		if(f==NULL) {
			PS_PrintToMainLog(env, "\nWarning: Output of graph cancelled (could not open file \"%s\").\n", "model.dot");
		} else {
			fprintf(f, "digraph model {\n");
			for(i=0; i<n; i++)
			if(i == start_index)
			fprintf(f, "	%1d [label=\"%1d\", shape=ellipse]\n", i, i);
			else if(yes_vec[i]> 0)
			fprintf(f, "	%1d [label=\"%1d\", shape=doublecircle]\n", i, i);
			else
			fprintf(f, "	%1d [label=\"%1d\", shape=circle]\n", i, i);
			count = n; l1 = h1 = l2 = h2 = 0;
			for(i=0; i<n; i++) {
				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++) {
					if (!use_counts) {
						l2 = choice_starts[j];
						h2 = choice_starts[j+1];
					} else {
						l2 = h2;
						h2 += choice_counts[j];
					}
					if(h2-l2>1) {
						fprintf(f, "	%1d [label=\"\", shape=point]\n", count);
						fprintf(f, "		%1d -> %1d\n", i, count);
						for(k=l2; k<h2; k++)
						fprintf(f, "		%1d -> %1d [label=\"%g\"]\n", count,
								cols[k], non_zeros[k]);
						count ++;
					} else
					fprintf(f, "		%1d -> %1d [label=\"%g\"]\n", i,
							cols[l2], non_zeros[l2]);
				}
			}
			fprintf(f, "}\n");
			fclose(f);
		}

		// Get the result from the solver
		if (ec_resume1(0) == PSUCCEED) { /* solve */
			stop = util_cpu_time();
			if (ec_get_double(ec_ref_get(Cost), &lp_result) == PSUCCEED)
			PS_PrintToMainLog(env, "Cost is %g\n", lp_result);
			else
			PS_PrintToMainLog(env, "Cost is ?\n");

			double *var_values;
			var_values = new double[num_lp_vars+n_maybe_r];
			count = 0;
			for (i=0; i<num_lp_vars+n_maybe_r; i++) {
				var_values[i] = 0.0;
				if (ec_get_double(ec_refs_get(Vars,i), &var_values[i]) == PSUCCEED) {
					if(var_values[i] != 0) {
						PS_PrintToMainLog(env, "X%d = %g		", i, var_values[i]);
						count++;
					}
				}
				//else
				//	PS_PrintToMainLog(env, "X%d = ?		");
				if(count == 8) {
					PS_PrintToMainLog(env, "\n");
					count = 0;
				}
			}
			if(count)
			PS_PrintToMainLog(env, "\n");

			/***********************************************************************
			 *
			 *
			 * generate adversary according the solution
			 *
			 *
			 ***********************************************************************/
			char c = '\0';
			/*if(n < 100) {
				printf("Would you like to see the adversary (Y/N)\n"); fflush(stdout);
				while (1) {
					scanf("%c", &c);
					if(c == 'Y' || c == 'y' || c == 'N' || c == 'n')
					break;
				}
			}*/
			c = 'n';

			if(c == 'Y' || c == 'y') {
				long sp_nodes = 0; // pointer of dot nodes
				long extra_node = n; // index of intermediate nodes for actions
				long sp_edges = 0; // pointer of dot edges
				//int queue[n]; // search queue
				int *queue;
				long head = 0; // pointer to the head of the queue
				long tail = 0; // pointer to the tail of the queue
				//int nodes[n]; // indicate whether the node has been visited
				int *nodes;

				queue = new int[n];
				nodes = new int[n];
				for(i=0; i<n; i++)
				nodes[i] = 0;

				if(checkReach) {// Adversary for reachability
					int (*dot_nodes)[2];
					int (*dot_edges)[2];
					double *edge_probabilities;
					//printf("allocated memory for local variables:\n"); 
					//printf("n=%1d, num_lp_vars=%1d, nnz=%1d\n", n, num_lp_vars, nnz);
					//fflush(stdout);

					//int dot_nodes[n+num_lp_vars][2]; // first entry: node number, second: shape 
					//int dot_edges[nnz+num_lp_vars][2]; // first entry: source node, second: target
					//double edge_probabilities[nnz+num_lp_vars]; // edge's probability

					dot_nodes = new int[n+num_lp_vars][2];
					dot_edges = new int[nnz+num_lp_vars][2];
					edge_probabilities = new double[nnz+num_lp_vars];
					//int dot_nodes[n+maybe_nc][2]; // first entry: node number, second: shape 
					//int dot_edges[nnz+maybe_nc][2]; // first entry: source node, second: target
					//double edge_probabilities[nnz+maybe_nc]; // edge's probability

					// put the initial state in the queue and dot_nodes
					queue[tail++] = start_index;
					dot_nodes[sp_nodes][0] = start_index;
					dot_nodes[sp_nodes++][1] = 0; // shape = 0: ellipse, 1: circle, 2: point, 3: doublecircle

					// recursive procedure
					int head_node = -1;
					while(head < tail) {
						head_node = queue[head++];
						nodes[head_node] = 1;
						// If the head node is a target state, do nothing
						if(yes_vec[head_node] == 0 && maybe_vec[head_node]> 0.0) {
							// the number of branches in the current state
							h1 = map_var[head_node+1] - map_var[head_node];
							//PS_PrintToMainLog(env, "h1 = %1d\n", h1);

							double sum = 0;
							for(i=0; i<h1; i++)
							sum += var_values[map_var[head_node] + i];

							for(i=0; i<h1; i++)
							// test if the action has non-zero probability
							if(var_values[map_var[head_node] + i]> 0) {
								// get all successor states of this action
								// First: locate the action
								if (!use_counts) {
									l1 = row_starts[head_node];
								} else {
									l1 = 0;
									for(j=0; j<head_node; j++)
									l1 += row_counts[j];
								}
								l1 += i;

								// Second: find all columns for this choice
								if (!use_counts) {
									l2 = choice_starts[l1];
									h2 = choice_starts[l1+1];
								} else {
									l2 = 0;
									for(j=0; j<l1; j++)
									l2 += choice_counts[j];
									h2 = l2 + choice_counts[j];
								}

								if(h2-l2>1) {
									// add an intermediate node to dot_nodes
									dot_nodes[sp_nodes][0] = extra_node;
									dot_nodes[sp_nodes++][1] = 2;
									// add an edge to dot_edges
									dot_edges[sp_edges][0] = head_node;
									dot_edges[sp_edges][1] = extra_node;
									// add the probability to the edge
									edge_probabilities[sp_edges++] = sum == 1 ? var_values[map_var[head_node] + i] :
									var_values[map_var[head_node] + i]/sum;

									for(j=l2; j<h2; j++) {
										// add the successor state to dot_nodes
										dot_nodes[sp_nodes][0] = cols[j];
										if(yes_vec[cols[j]]> 0)
										dot_nodes[sp_nodes++][1] = 3;
										else
										dot_nodes[sp_nodes++][1] = 1;

										// add an edge to dot_edges
										dot_edges[sp_edges][0] = extra_node;
										dot_edges[sp_edges][1] = cols[j];
										// add the probability to the edge
										edge_probabilities[sp_edges++] = non_zeros[j];

										if(maybe_vec[cols[j]]> 0 && nodes[cols[j]] == 0)
										queue[tail++] = cols[j];
									}
									extra_node++;
								} else {
									// add the successor state to dot_nodes
									dot_nodes[sp_nodes][0] = cols[l2];
									if(yes_vec[cols[l2]]> 0)
									dot_nodes[sp_nodes++][1] = 3;
									else
									dot_nodes[sp_nodes++][1] = 1;
									// add an edge to dot_edges
									dot_edges[sp_edges][0] = head_node;
									dot_edges[sp_edges][1] = cols[l2];
									// add the probability to the edge
									//edge_probabilities[sp_edges++] = var_values[map_var[head_node] + i];
									edge_probabilities[sp_edges++] = sum == 1 ? var_values[map_var[head_node] + i] :
									var_values[map_var[head_node] + i]/sum;

									if(maybe_vec[cols[l2]]> 0 && nodes[cols[l2]] == 0)
									queue[tail++] = cols[l2];
								}
							}
						}
					}

					// write to file
					f = fopen("adversary.dot", "w"); /* create a file for writing */
					f = NULL;
					if(f==NULL) {
						PS_PrintToMainLog(env, "\nWarning: Output of adversary cancelled (could not open file \"%s\").\n", "adversary.dot");
					} else {
						fprintf(f, "digraph adversary {\n");
						for(i=0; i<sp_nodes; i++)
						if(dot_nodes[i][0] >= n)
						fprintf(f, "	%1d [label=\"\", shape=point]\n", dot_nodes[i][0]);
						else
						fprintf(f, "	%1d [label=\"%1d\", shape=%s]\n", dot_nodes[i][0],
								dot_nodes[i][0],
								((dot_nodes[i][1]==0)? "ellipse" :
										((dot_nodes[i][1]==1)? "circle" : "doublecircle")));
						for(i=0; i<sp_edges; i++)
						fprintf(f, "		%1d -> %1d [label=\"%g\"]\n", dot_edges[i][0],
								dot_edges[i][1], edge_probabilities[i]);
						fprintf(f, "}\n");
						fclose(f);
					}
				} else {// Adversary for LTL formulas
					int (*dot_nodes)[2];
					int (*dot_edges)[2];
					double *edge_probabilities;
					int *terminate_nodes;
					//printf("allocated memory for local variables:\n"); 
					//printf("n=%1d, num_lp_vars=%1d, nnz=%1d\n", n, num_lp_vars, nnz);
					//fflush(stdout);

					//int dot_nodes[n+num_lp_vars][2]; // first entry: node number, second: shape 
					//int dot_edges[nnz+num_lp_vars][2]; // first entry: source node, second: target
					//double edge_probabilities[nnz+num_lp_vars]; // edge's probability

					dot_nodes = new int[n+num_lp_vars][2];
					dot_edges = new int[nnz+num_lp_vars][2];
					edge_probabilities = new double[nnz+num_lp_vars];
					terminate_nodes = new int[n]; // the entry stores the node number for newly added terminate nodes

					//printf("Done\n");
					//fflush(stdout);

					for(i=0; i<n; i++)
					terminate_nodes[i] = -1;

					// put the initial state in the queue and dot_nodes
					queue[tail++] = start_index;
					dot_nodes[sp_nodes][0] = start_index;
					dot_nodes[sp_nodes++][1] = 0; // shape = 0: ellipse, 1: circle, 2: point, 3: doublecircle, 4: box

					// recursive procedure
					int head_node = -1;
					while(head < tail) {
						head_node = queue[head++];

						if(nodes[head_node] == 1)
						continue;
						nodes[head_node] = 1;
						// If the head node is a target state, do nothing
						if(yes_vec[head_node]> 0 || maybe_vec[head_node]> 0) {
							// the number of branches in the current state
							h1 = map_var[head_node+1] - map_var[head_node];

							double sum = 0;
							for(i=0; i<h1; i++)
							sum += var_values[map_var[head_node] + i];

							for(i=0; i<h1; i++)
							// test if the action has non-zero probability
							if(var_values[map_var[head_node] + i]> 0) {
								if(yes_vec[head_node]> 0 && i==h1-1) { // extra transition
									// add an intermediate node to dot_nodes
									if(terminate_nodes[head_node] < 0) {
										dot_nodes[sp_nodes][0] = extra_node;
										dot_nodes[sp_nodes++][1] = 4;
										dot_edges[sp_edges][1] = extra_node;
										terminate_nodes[head_node] = extra_node;
									} else
									dot_edges[sp_edges][1] = terminate_nodes[head_node];
									// add an edge to dot_edges
									dot_edges[sp_edges][0] = head_node;
									// add the probability to the edge
									//edge_probabilities[sp_edges++] = var_values[map_var[head_node] + i];
									edge_probabilities[sp_edges++] = sum == 1 ? var_values[map_var[head_node] + i] :
									var_values[map_var[head_node] + i]/sum;
									extra_node++;
								} else {
									// get all successor states of this action
									// First: locate the action
									if (!use_counts) {
										l1 = row_starts[head_node];
									} else {
										l1 = 0;
										for(j=0; j<head_node; j++)
										l1 += row_counts[j];
									}
									l1 += i;

									// Second: find all columns for this choice
									if (!use_counts) {
										l2 = choice_starts[l1];
										h2 = choice_starts[l1+1];
									} else {
										l2 = 0;
										for(j=0; j<l1; j++)
										l2 += choice_counts[j];
										h2 = l2 + choice_counts[j];
									}
									if(h2-l2>1) {
										// add an intermediate node to dot_nodes
										dot_nodes[sp_nodes][0] = extra_node;
										dot_nodes[sp_nodes++][1] = 2;
										// add an edge to dot_edges
										dot_edges[sp_edges][0] = head_node;
										dot_edges[sp_edges][1] = extra_node;
										// add the probability to the edge
										//edge_probabilities[sp_edges++] = var_values[map_var[head_node] + i];
										edge_probabilities[sp_edges++] = sum == 1 ? var_values[map_var[head_node] + i] :
										var_values[map_var[head_node] + i]/sum;

										for(j=l2; j<h2; j++) {
											// add the successor state to dot_nodes
											dot_nodes[sp_nodes][0] = cols[j];
											if(yes_vec[cols[j]]> 0)
											dot_nodes[sp_nodes++][1] = 3;
											else
											dot_nodes[sp_nodes++][1] = 1;

											// add an edge to dot_edges
											dot_edges[sp_edges][0] = extra_node;
											dot_edges[sp_edges][1] = cols[j];
											// add the probability to the edge
											edge_probabilities[sp_edges++] = non_zeros[j];

											if((maybe_vec[cols[j]]> 0 || yes_vec[cols[j]]> 0) &&
													nodes[cols[j]] == 0)
											queue[tail++] = cols[j];
										}
										extra_node++;
									} else {
										// add the successor state to dot_nodes
										dot_nodes[sp_nodes][0] = cols[l2];
										if(yes_vec[cols[l2]]> 0)
										dot_nodes[sp_nodes++][1] = 3;
										else
										dot_nodes[sp_nodes++][1] = 1;
										// add an edge to dot_edges
										dot_edges[sp_edges][0] = head_node;
										dot_edges[sp_edges][1] = cols[l2];
										// add the probability to the edge
										//edge_probabilities[sp_edges++] = var_values[map_var[head_node] + i];
										edge_probabilities[sp_edges++] = sum == 1 ? var_values[map_var[head_node] + i] :
										var_values[map_var[head_node] + i]/sum;

										if((maybe_vec[cols[l2]]> 0 || yes_vec[cols[l2]]> 0) &&
												nodes[cols[l2]] == 0) {
											queue[tail++] = cols[l2];
										}
									}
								}
							}
						}
					}

					printf("generating adversary file\n"); fflush(stdout);

					// write to file
					f = fopen("adversary.dot", "w"); /* create a file for writing */
					if(f==NULL) {
						PS_PrintToMainLog(env, "\nWarning: Output of adversary cancelled (could not open file \"%s\").\n", "adversary.dot");
					} else {
						fprintf(f, "digraph adversary {\n");
						for(i=0; i<sp_nodes; i++)
						if(dot_nodes[i][0] >= n) {
							if(dot_nodes[i][1] == 2)
							fprintf(f, "	%1d [label=\"\", shape=point]\n", dot_nodes[i][0]);
							else
							fprintf(f, "	%1d [label=\"\", shape=box, fillcolor=black]\n", dot_nodes[i][0]);
						} else
						fprintf(f, "	%1d [label=\"%1d\", shape=%s]\n", dot_nodes[i][0],
								dot_nodes[i][0],
								((dot_nodes[i][1]==0)? "ellipse" :
										((dot_nodes[i][1]==1)? "circle" : "doublecircle")));
						for(i=0; i<sp_edges; i++)
						fprintf(f, "		%1d -> %1d [label=\"%g\"]\n", dot_edges[i][0],
								dot_edges[i][1], edge_probabilities[i]);
						fprintf(f, "}\n");
						fclose(f);
					}

					delete[] dot_nodes;
					delete[] dot_edges;
					delete[] edge_probabilities;
					delete[] terminate_nodes;
				}

				delete[] queue;
				delete[] nodes;
			}
			delete[] var_values;
		}
		else PS_PrintToMainLog(env, "No solution\n");

		// Destroy variables and clean up
		ec_refs_destroy(Vars);
		ec_ref_destroy(Cost);
		ec_cleanup();

		//stop = util_cpu_time();
		time_for_lp = (double)(stop - start2)/1000;
		time_taken = (double)(stop - start3)/1000;
		PS_PrintToMainLog(env, "\nCBC: LP problem solved in %.2f seconds (setup %.2f, lpsolve %.2f)\n", time_taken, time_for_setup, time_for_lp);

		delete yes_vec;
		delete map_var;
		// catch exceptions: register error, free memory
	} catch (std::bad_alloc e) {
		PS_SetErrorMessage("Out of memory");
		if (yes2_vec) delete[] yes2_vec;
		yes2_vec = 0;
		if (target_ptrs) env->ReleaseLongArrayElements(_targets, target_ptrs, 0);
		if (relops) env->ReleaseIntArrayElements(_relops, relops, 0);
		if (bounds) env->ReleaseDoubleArrayElements(_bounds, bounds, 0);
	} catch (const char *err) {
		PS_SetErrorMessage(err);
		if (yes2_vec) delete yes2_vec;
		yes2_vec = 0;
		if (target_ptrs) env->ReleaseLongArrayElements(_targets, target_ptrs, 0);
		if (relops) env->ReleaseIntArrayElements(_relops, relops, 0);
		if (bounds) env->ReleaseDoubleArrayElements(_bounds, bounds, 0);
	}

	// free memory
	if (a) Cudd_RecursiveDeref(ddman, a);
	if (trans_rewards) Cudd_RecursiveDeref(ddman, trans_rewards);
	if (ndsm) delete ndsm;
	if (ndsm_r) delete ndsm_r;
	//if (yes1_vec) delete[] yes1_vec;
	delete[] maybe_vec;
	for(i=0; i<num_targets; i++)
	delete[] yes_vecs[i];
	delete[] arr_reals;
	delete[] arr_ints;
	if (target_ptrs) env->ReleaseLongArrayElements(_targets, target_ptrs, 0);
	if (relops) env->ReleaseIntArrayElements(_relops, relops, 0);
	if (bounds) env->ReleaseDoubleArrayElements(_bounds, bounds, 0);

	return lp_result;
}

//------------------------------------------------------------------------------