MDPNonlinearEqn.run

# Call solver and give it options

# TODO run simulation with estimated parameters

# Load model and data
model "MarkovActv.mod";
data "MarkovActv.dat";

# Is Cauchy distribution used ? 
fix IS_CAUCHY := 1;

# logging options
let debug_log := 1;

# include the code that define the state and choice set
include MDPStateAction.run

# Define the component MDP problem for each individual
problem componentMDP {n in PERS}:

# Choose the objective function
likelihood0,

# List the variables
{t in 0..H, j in AUW[n]} EV[n,t,j], 
{t in 0..H, j in AUW[n]} actvUtil[n,t,j],
{(t,j) in X[n], (k,h) in DA[n,t,j]} sumActvUtil[n,t,j,k,h],
{(t,j) in X[n], (k,h) in D[n,t,j]} choiceUtil[n,t,j,k,h],
{(t,j) in X[n], (k,h) in D[n,t,j]} choiceProb[n,t,j,k,h],
{j in ACTV} Um[n,j],
{j in ACTV} b[n,j],
{j in ACTV} c[n,j],
# EV, actvUtil, choiceUtil, choiceProb, VoT, theta, Uw, xi, gamma,

# List the constraints
{(t,j) in X[n]} Bellman_Eqn[n,t,j],
Bellman_EqnH[n];

# Solve the component MDPs separately
for {n in PERS} {

	# Set the current problem
	problem componentMDP[n];

	fix {j in ACTV} Um[n,j] := Um0[n,j];
	fix {j in ACTV}  b[n,j] := b0[n,j];
	fix {j in ACTV}  c[n,j] := c0[n,j];

	# Set at a trivial initial value
	# let {t in 0..H, j in ACTVWORK[n]} EV[n,t,j] := initEV;

	# Specify KNITRO solver options:
	option solver "C:\Ziena\KNITRO900\knitroampl\knitroampl.exe";
	option knitro_options "alg=2 hessopt=1 outlev=3 maxit=300 xtol=0.0000000001 wantsol=1";
	# option presolve 0;

	# Solve command
	solve;

	# Solution status
	printf "Component MDP %d\n", n > DATA/MDP.sol;
	if match (solve_message, "Locally optimal solution") > 0 then { 
		printf "%s\n", "Optimal solution found" > DATA/MDP.sol;
		printf "%5.0f \n", 0 > DATA/MDP.sol ; 
	}
	else if match (solve_message, "Iteration limit reached") > 0 then { 
		printf "%s\n", "Iteration limit reached" > DATA/MDP.sol;
		printf "%5.0f \n", 400 > DATA/MDP.sol ; 
	}
	else if match (solve_message, "Evaluation error") > 0 then { 
		printf "%s\n", "Evaluation error" > DATA/MDP.sol;
		printf "%5.0f \n", 502 > DATA/MDP.sol ; 
	}
	else {
		printf "%s\n", "No solution" > DATA/MDP.sol;
		printf "%5.0f \n", 1000 > DATA/MDP.sol ; 
	}

	# display _solve_time > DATA/MDP.sol;

	# write the structural parameters
	display beta, VoT, theta > DATA/MDP.sol;
	display Um, b, c > DATA/MDP.sol;

	# also display the structural parameters
	printf "\n\nComponent MDP %d", n;
	printf "\n\nstructural parameters (all fixed):\n\n";
	display beta, VoT, theta;
	display Um, b, c;
	# display Uw, xi, gamma, lambda;
}


# Output commands
option display_round 6, display_width 120;

# write the activity utility and travel cost
# display actvUtil > DATA/actvUtil.dat;
# display sumActvUtil > DATA/sumActvUtil.dat;
# display sumTravelCost > DATA/sumTravelCost.dat;

# write the value function
# display EV > DATA/MDP.sol;
# write the choice-specific utility
# display choiceUtil > DATA/MDP.sol;
# write the choice probability
# display choiceProb > DATA/MDP.sol;

# write EV to a AMPL data file EUtil.dat
# display EV > DATA/EUtil.dat;
# write choiceProb to a AMPL data file Prob.dat
# display choiceProb > DATA/Prob.dat;


# export flow variables
include ExportNetworkFlows.run;

# export EV to file EUtil.m
printf "EV = zeros(%d, %d, %d);\n", N, H, M > DATA/EUtil.m;
for {t in TIME} {
	printf "EV(:, %3d, :) = [", t+1 > DATA/EUtil.m;
	for {n in PERS} {
		for {j in AUW[n]} {
			if (t,j) in X[n] then
				printf " %8.2f", EV[n,t,j] > DATA/EUtil.m;
			else
				printf " %8s", "nan" > DATA/EUtil.m;
		}
		printf ";" > DATA/EUtil.m;
	}
	printf "];\n" > DATA/EUtil.m;
}

# debug output
if debug_log == 1 then {

# export choiceProb to file choiceProb.m
printf "Pr = zeros(%d, %d, %d, %d, %d);\n", N, H, M, M, DH > DATA/choiceProb.m;
for {n in PERS} {
	for {(t,j) in X[n]} {
		printf "Pr(%d, %d, %d, :, :) = [\n", n, t+1, j > DATA/choiceProb.m;
		for {h in 1..DH} {
			for {k in AUW[n]} {
				if (k,h) in D[n,t,j] then
					if choiceProb[n,t,j,k,h] >= -0.01 and
					   choiceProb[n,t,j,k,h] <=  1.01 then
						printf " %8.5f", choiceProb[n,t,j,k,h] > DATA/choiceProb.m;
					else
						printf " %8s", "nan" > DATA/choiceProb.m;
				else
					printf " %8.5f", 0.0 > DATA/choiceProb.m;
			}
			printf "\n" > DATA/choiceProb.m;
		}
		printf "]';\n" > DATA/choiceProb.m;
	}
}

# export exp(choiceProb) to file expChoiceProb.txt
printf "ePr = zeros(%d, %d, %d, %d, %d);\n", N, H, M, M, DH > DATA/expChoiceProb.txt;
for {n in PERS} {
	for {(t,j) in X[n]} {
		printf "ePr(%d, %d, %d, :, :) = [\n", n, t+1, j > DATA/expChoiceProb.txt;
		for {h in 1..DH} {
			for {k in AUW[n]} {
				if (k,h) in D[n,t,j] then
					printf " %8.2f", ( theta * ( choiceUtil[n,t,j,k,h] + 
												 beta**h * EV[n,(t+h),k]) - 
									   theta * EV[n,t,j]) > DATA/expChoiceProb.txt;
				else
					printf " %8s", "nan" > DATA/expChoiceProb.txt;
			}
			printf "\n" > DATA/expChoiceProb.txt;
		}
		printf "]';\n" > DATA/expChoiceProb.txt;
	}
}

# end of log
}

# export fixed parameters to MDP.m
printf "beta = %f;\n", beta > DATA/MDP.m;
printf "theta = %f;\n", theta > DATA/MDP.m;

# export the structural parameters to MDP.m
printf "VoT0 = %f;\n", VoT > DATA/MDP.m;
printf "Um0 = [ " > DATA/MDP.m;
for {n in PERS} {
	for {j in ALLACTV} {
		printf "%f ", Um[n,j] > DATA/MDP.m;
	}
	printf "\n" > DATA/MDP.m;
}
printf "];\n" > DATA/MDP.m;
printf "b0 = [ " > DATA/MDP.m;
for {n in PERS} {
	for {j in ALLACTV} {
		printf "%f ", b[n,j] > DATA/MDP.m;
	}
	printf "\n" > DATA/MDP.m;
}
printf "];\n" > DATA/MDP.m;
printf "c0 = [ " > DATA/MDP.m;
for {n in PERS} {
	for {j in ALLACTV} {
		printf "%f ", c[n,j] > DATA/MDP.m;
	}
	printf "\n" > DATA/MDP.m;
}
printf "];\n" > DATA/MDP.m;