For review: mercury linear inequation solver/termination analysis mods.

[Thomas Charles CONWAY]

Hi

The following diff contains changes necessary to make use of a
linear inequality solver written in Mercury rather than calling
lpsolve.

It does not address the issue of whether the solver belongs in the
compiler or in a separate library. I've taken the former approach.

-- 
ZZ:wq!
^X^C
Thomas Conway               				      conway at cs.mu.oz.au
AD DEUM ET VINUM	  			      Every sword has two edges.


Modify termination analysis to use a linear inequation solver
written in Mercury rather than invoking an external 3rd party one.

compiler/term_pass1.m:
	Remove all the code for writing out the system of equations,
	invoking the solver and parsing the output.
	Add code to convert the system of equations into the format
	expected by the solver, and call the solver.

compiler/term_errors.m:
	the 'lpsolve_failed' constructor goes from arity 1 to arity 0.
	modify the output of error messages accordingly.

compiler/lp.m:
	NEW file. Contains the linear inequation solver based on
	the standard simplex method.

cvs diff: Diffing .
Index: term_errors.m
===================================================================
RCS file: /home/staff/zs/imp/mercury/compiler/term_errors.m,v
retrieving revision 1.2
diff -u -r1.2 term_errors.m
--- term_errors.m	1997/10/09 09:39:12	1.2
+++ term_errors.m	1997/10/12 23:24:32
@@ -85,7 +85,7 @@
 	;	not_subset(pred_proc_id, bag(var), bag(var))
 	;	not_dag
 	;	no_eqns
-	;	lpsolve_failed(eqn_soln)
+	;	lpsolve_failed
 	;	call_in_single_arg(pred_proc_id)
 			% single argument analysis did not find a head
 			% variable that was decreasing in size. 
@@ -473,49 +473,12 @@
 	io__write_string("  and functions in this SCC\n").
 
 term_errors__output_2(_PredId, _ProcId, _Module, _ConstErrorOutput, ForHLDSDump,
-		Context - lpsolve_failed(Reason)) -->
+		Context - lpsolve_failed) -->
 	io__write_string("the constraint solver "),
-	(
- 		{ Reason = optimal },
-		{ error("term_errors:output_2 Unexpected return value from lp_solve") }
-	;
-		{ Reason = infeasible },
-		io__write_string("found the\n"),
-		maybe_write_string(ForHLDSDump, "% "),
-		prog_out__write_context(Context),
-		io__write_string("  constraints that the analysis produced to be infeasible\n")
-	;
-		{ Reason = unbounded },
-		io__write_string("found that the\n"),
-		maybe_write_string(ForHLDSDump, "% "),
-		prog_out__write_context(Context),
-		io__write_string(
-			"  constraints were not strong enough to put a\n"),
-		maybe_write_string(ForHLDSDump, "% "),
-		prog_out__write_context(Context),
-		io__write_string(
-			"  bound on the value of the change constants\n")
-	;
-		{ Reason = failure },
-		io__write_string("was unable to\n"),
-		maybe_write_string(ForHLDSDump, "% "),
-		prog_out__write_context(Context),
-		io__write_string(
-			"  solve the constraints that were produced\n")
-	;
-		{ Reason = fatal_error },
-		io__set_exit_status(1),
-		io__write_string("was unable to create and/or\n"),
-		maybe_write_string(ForHLDSDump, "% "),
-		prog_out__write_context(Context),
-		io__write_string("  remove the temporary files it required\n")
-	;
-		{ Reason = parse_error },
-		{ error("term_errors:output_2 Unexpected return value from lp_solve") }
-	;
-		{ Reason = solved(_) },
-		{ error("term_errors:output_2 Unexpected return value from lp_solve") }
-	).
+	io__write_string("found the\n"),
+	maybe_write_string(ForHLDSDump, "% "),
+	prog_out__write_context(Context),
+	io__write_string("  constraints that the analysis produced to be infeasible\n").
 
 % call_in_single_arg will only be printed out as the second part of the
 % single_arg_failed(NormErr, SingleErr) error.  Therefore, the following
Index: term_pass1.m
===================================================================
RCS file: /home/staff/zs/imp/mercury/compiler/term_pass1.m,v
retrieving revision 1.2
diff -u -r1.2 term_pass1.m
--- term_pass1.m	1997/10/09 09:39:15	1.2
+++ term_pass1.m	1997/10/17 03:37:44
@@ -49,8 +49,8 @@
 :- implementation.
 
 :- import_module int, list, bag, require, bool, std_util, char, map, string.
-:- import_module hlds_pred, hlds_goal, hlds_data. 
-:- import_module term_errors, mode_util, type_util, term.
+:- import_module hlds_pred, hlds_goal, hlds_data, float.
+:- import_module term_errors, mode_util, type_util, term, varset, lp.
 
 %------------------------------------------------------------------------------
 
@@ -856,437 +856,94 @@
 :- mode solve_equations(in, in, in, in, in, out, di, uo) is det.
 
 solve_equations(Module0, Equations, PPIds, UsedArgs, Context, Module) -->
-	io__tmpnam(ConstraintFile),
-	solve_equations_create_constraint_file(Equations, PPIds, 
-		ConstraintFile, ConstraintSucc),
-	( { ConstraintSucc = yes } ->
-		solve_equations_solve_constraint_file(ConstraintFile, Soln)
-	;
-		% Failed to create the constraint file.
-		{ Soln = fatal_error }
-	),
-	% The equations have been solved, now put the
-	% results into the module_info.
-	( { Soln = solved(SolutionList) } ->
-		% The solver successfully solved the constraints.
-		{ module_info_preds(Module0, PredTable0) },
-		{ proc_inequalities_set_module(SolutionList, UsedArgs, 
-			PredTable0, PredTable) },
-		{ module_info_set_preds(Module0, PredTable, 
-			Module) }
-	; { Soln = optimal } ->
-		% All 'optimal' results should have been
-		% changed into a list of solutions in
-		% solve_equations.  
-		{ error("term_pass1__solve_equations: Unexpected Value\n")}
-	;
-		% The constraint solver failed to solve the
-		% set of constraints - set the termination
-		% constant to infinity.
-		{ Error = Context - lpsolve_failed(Soln) },
-		{ set_pred_proc_ids_const(PPIds, 
-			inf(Error), Module0, Module) }
-	).
-
-:- pred solve_equations_solve_constraint_file(string, eqn_soln, 
-	io__state, io__state).
-:- mode solve_equations_solve_constraint_file(in, out, di, uo) is det.
-solve_equations_solve_constraint_file(ConstraintFile, Soln) -->
-	io__tmpnam(OutputFile),
-	% run lp_solve
-	solve_equations_run_lp_solve(ConstraintFile, OutputFile, MaybeResult),
-	( 
-		{ MaybeResult = yes(Result) },
-		( { Result = optimal } ->
-			solve_equations_process_output_file(OutputFile, Soln0)
-		;
-			% lp_solve failed to solve the constraints.
-			% This could be for a number of reasons,
-			% and the value of Result will represent
-			% the reason.
-			{ Soln0 = Result }
-		)
-	;
-		{ MaybeResult = no },
-		{ Soln0 = fatal_error }
-	),
-
-	% Remove and close all temporary files.
-	solve_equations_remove_file(ConstraintFile, Success0),
-	solve_equations_remove_file(OutputFile, Success1),
-	{ bool__or(Success0, Success1, Success) },
-	{ ( 
-		Success = yes,
-		Soln = Soln0
-	;
-		Success = no,
-		Soln = fatal_error
-	) }.
-
-% This runs lp_solve, and outputs an error message and returns `no' 
-% if the system call fails.  On success, it returns the return value of
-% lp_solve after the return value has been changed from an integer into a
-% lpsolve_ret_val type.
-:- pred solve_equations_run_lp_solve(string, string, maybe(eqn_soln),
-	io__state, io__state).
-:- mode solve_equations_run_lp_solve(in, in, out, di, uo) is det.
-solve_equations_run_lp_solve(ConstraintFile, OutputFile, MaybeResult) -->
-	io__get_environment_var("MERCURY_LPSOLVE", Lpsolve0),
-	( { Lpsolve0 = yes(Lpsolve1) } ->
-		{ Lpsolve = Lpsolve1 }
-	;
-		{ Lpsolve = "lp_solve" }
-	),
-	{ string__append_list([ Lpsolve, " <",
-		ConstraintFile, " > ",
-		OutputFile], Command) },
-	io__call_system(Command, Res0),
-	% Test the return value
-	( 
-		{ Res0  = ok(RetVal) },
-		{ lpsolve_ret_val(RetVal, Result) },
-		{ MaybeResult = yes(Result) }
-	;
-		{ Res0 = error(Error0) },
-		io__progname_base("term_pass1.m", ProgName),
-		{ io__error_message(Error0, Msg0) },
-		io__write_strings([
-			ProgName,
-			": Error with system call `",
-			Command,
-			"': ",
-			Msg0,
-			"\n" ]),
-		io__set_exit_status(1),
-		{ MaybeResult = no }
-	).
-
-:- pred solve_equations_remove_file(string, bool, io__state, io__state).
-:- mode solve_equations_remove_file(in, out, di, uo) is det.
-solve_equations_remove_file(File, Success) -->
-	io__remove_file(File, Res1),
-	( { Res1 = error(Error1) } ->
-		io__progname_base("term_pass1.m", ProgName),
-		{ io__error_message(Error1, Msg1) },
-		io__write_strings([
-			ProgName,
-			": Error deleting temporary file `",
-			File,
-			"' : ",
-			Msg1,
-			"\n" ]),
-		io__set_exit_status(1),
-		{ Success = no }
-	;
-		{ Success = yes }
-	).
-
-% This predicate creates the constraint file in a format suitable for
-% lp_solve.  This really shouldn't be called with Equations=[] as lp_solve
-% exits with an error if it is called without any constraints.
-:- pred solve_equations_create_constraint_file(list(term_pass1__equation),
-	list(pred_proc_id), string, bool, io__state, io__state).
-:- mode solve_equations_create_constraint_file(in, in, in, out, di, uo) is det.
-
-solve_equations_create_constraint_file(Equations, PPIds, 
-		ConstraintFile, Success) -->
-	io__open_output(ConstraintFile, Res),
-	( 	{ Res = error(Error) },
-		% error message and quit
-		io__progname_base("termination.m", ProgName),
-		{ io__error_message(Error, Msg) },
-	
-		io__write_string(ProgName),
-		io__write_string(": cannot open temporary file `"),
-		io__write_string(ConstraintFile),
-		io__write_string("' for output: "),
-		io__write_string(Msg),
-		io__write_string("\n"),
-		io__set_exit_status(1),
-		{ Success = no }
-	;
-		{ Res = ok(Stream) },
-		( { Equations = [] } ->
-			{ Success = no }
-		;
-			io__set_output_stream(Stream, OldStream),
-			% create the constraint file
-			output_equations(Equations, PPIds, Success),
-	
-			io__set_output_stream(OldStream, _),
-			io__close_output(Stream)
-		)
-	).
-
-% Prepare to parse the output from lp_solve.
-:- pred solve_equations_process_output_file(string, eqn_soln, 
-	io__state, io__state).
-:- mode solve_equations_process_output_file(in, out, di, uo) is det.
-solve_equations_process_output_file(OutputFile, Soln) -->
-	io__open_input(OutputFile, OutputRes),
-	( 	{ OutputRes = error(Error) },
-		io__progname_base("term_pass1.m", ProgName),
-		{ io__error_message(Error, Msg) },
-
-		io__write_string(ProgName),
-		io__write_string(": cannot open temporary file `"),
-		io__write_string(OutputFile),
-		io__write_string("' for input: "),
-		io__write_string(Msg),
-		io__write_string("\n"),
-		io__set_exit_status(1),
-		{ Soln = fatal_error }
-	;
-		{ OutputRes = ok(Stream) },
-		io__set_input_stream(Stream, OldStream),
-		% need to interpret it now
-		solve_equations_parse_output_file(Soln),
-		io__set_input_stream(OldStream, _),
-		io__close_input(Stream)
-	).
-
-% Parse the output from lp_solve.
-:- pred solve_equations_parse_output_file(eqn_soln, io__state, io__state).
-:- mode solve_equations_parse_output_file(out, di, uo) is det.
-solve_equations_parse_output_file(Soln) -->
-	io__read_line(Res1),
-	( { Res1 = ok(_) } ->
-		solve_equations_parse_output_file_2(Soln0, MaybeBVal),
-		( { Soln0 = solved(Result0) } ->
-			( { MaybeBVal = yes(BVal) } ->
-				{ solve_equations_output_file_2(Result0, BVal,
-					Result) },
-				{ Soln = solved(Result) }
-			;
-				{ Soln = parse_error }
-			)
-		;
-			io__write_string(
-				"parse_output_file returned not solved\n"),
-			{ Soln = parse_error }
-		)
-	;
-		{ Soln = parse_error }
-	).
-
-:- pred solve_equations_output_file_2(list(pair(pred_proc_id, int)), int,
-	list(pair(pred_proc_id, int))).
-:- mode solve_equations_output_file_2(in, in, out) is det.
-solve_equations_output_file_2([], _, []). 
-solve_equations_output_file_2([X | Xs], BVal, [Y | Ys]) :-
-	X = PPId - XVal, 	% pair deconstruction
-	YVal is XVal - BVal,	% subtraction
-	Y = PPId - YVal,	% pair construction
-	solve_equations_output_file_2(Xs, BVal, Ys).
-		
-:- pred solve_equations_parse_output_file_2(eqn_soln, maybe(int), io__state,
-	io__state).
-:- mode solve_equations_parse_output_file_2(out, out, di, uo) is det.
-solve_equations_parse_output_file_2(Soln, BVal) -->
-	io__read_line(Res1),
-	( { Res1 = ok([ X | Xs ]) } ->
-		( 
-			{ X = 'a' },
-			{ char_list_remove_int(Xs, PredInt, Xs1) },
-			{ Xs1 = [ '_' | Xs2 ] },
-			{ char_list_remove_int(Xs2, ProcInt, Xs3) },
-			{ char_list_remove_whitespace(Xs3, Xs4) },
-			{ char_list_remove_int(Xs4, Value, _Xs5) }
-		->
-			% Have found a solution.
-			{ pred_id_to_int(PredId, PredInt) },
-			{ proc_id_to_int(ProcId, ProcInt) },
-			solve_equations_parse_output_file_2(Soln0, BVal),
-			( { Soln0 = solved(SolnList) } ->
-				{ NewSoln = proc(PredId, ProcId) - Value },
-				{ Soln = solved([NewSoln | SolnList ]) }
-			;
-				{ Soln = Soln0 }
-			)
-		; % else if
-			{ X = 'b' },
-			{ char_list_remove_whitespace(Xs, Xs1) },
-			{ char_list_remove_int(Xs1, Value, _Xs2) }
-		->
-			solve_equations_parse_output_file_2(Soln, _Bval),
-			{ BVal = yes(Value) }
-		;
-			{ Soln = parse_error },
-			{ BVal = no }
+	(
+		{ convert_equations(Equations, Varset, LPEquations,
+			Objective, PPVars) }
+	->
+		{ map__values(PPVars, AllVars0) },
+		{ list__sort_and_remove_dups(AllVars0, AllVars) },
+		lp_solve(LPEquations, min, Objective, Varset, AllVars, Soln),
+		(
+			{ Soln = unsatisfiable },
+			{ Error = Context - lpsolve_failed },
+			{ set_pred_proc_ids_const(PPIds, 
+				inf(Error), Module0, Module) }
+		;
+			{ Soln = satisfiable(_ObjVal, SolnVal) },
+			{ list__map(lookup_coeff(PPVars, SolnVal), PPIds,
+				SolutionList) },
+			{ module_info_preds(Module0, PredTable0) },
+			{ proc_inequalities_set_module(SolutionList, UsedArgs, 
+				PredTable0, PredTable) },
+			{ module_info_set_preds(Module0, PredTable, 
+				Module) }
 		)
-	; { Res1 = eof } ->
-		{ Soln = solved([]) },
-		{ BVal = no }
-	;
-		{ Soln = parse_error },
-		{ BVal = no }
-	).
-
-:- pred char_list_remove_int(list(char), int, list(char)).
-:- mode char_list_remove_int(in, out, out) is semidet.
-char_list_remove_int([X | Xs], Int, ListOut) :-
-	char__is_digit(X),
-	char__to_int(X, Int0),
-	char__to_int('0', IntValueofZero),
-	Int1 is Int0 - IntValueofZero,   
-	char_list_remove_int_2(Xs, Int1, Int, ListOut).
-
-:- pred char_list_remove_int_2(list(char), int, int, list(char)).
-:- mode char_list_remove_int_2(in, in, out, out) is semidet.
-
-char_list_remove_int_2([], Int, Int, []).
-char_list_remove_int_2([X | Xs], Int0, Int, ListOut) :-
-	( char__is_digit(X) ->
-		char__to_int('0', IntValueofZero),
-		char__to_int(X, Int1),
-		Int2 is Int0 * 10 + Int1 - IntValueofZero,
-		char_list_remove_int_2(Xs, Int2, Int, ListOut)
-	;
-		ListOut = [ X | Xs ],
-		Int = Int0
-	).
-		
-:- pred char_list_remove_whitespace(list(char), list(char)).
-:- mode char_list_remove_whitespace(in, out) is det.
-char_list_remove_whitespace([], []).
-char_list_remove_whitespace([ X | Xs ], Out) :-
-	( char__is_whitespace(X) ->
-		char_list_remove_whitespace(Xs, Out)
-	;
-		Out = [ X | Xs ]
-	).
-
-:- pred lpsolve_ret_val(int, eqn_soln).
-:- mode lpsolve_ret_val(in, out) is det.
-lpsolve_ret_val(Int, Result) :-
-	( Int = 0	->	Result = optimal
-	; Int = 2	->	Result = infeasible
-	; Int = 3	->	Result = unbounded
-	; 			Result = failure
-	).
-
-%-----------------------------------------------------------------------------%
-% These predicates are used to output a list of equations in a form
-% suitable for lp_solve.  
-:- pred output_equations(list(term_pass1__equation), list(pred_proc_id),
-	bool, io__state , io__state).
-:- mode output_equations(in, in, out, di, uo) is det.
-
-output_equations(Xs, PPIds, Success) --> 
-	% output: 'max: # b - PPIds'
-	io__write_string("max: "),
-	{ list__length(PPIds, Length) },
-	io__write_int(Length),
-	io__write_string(" b "),
-	output_eqn_2(PPIds),
-	io__write_string(";\n"),
-
-	output_equations_2(Xs, 1, Success).
-
-:- pred output_equations_2(list(term_pass1__equation), int,
-	bool, io__state , io__state).
-:- mode output_equations_2(in, in, out, di, uo) is det.
-
-output_equations_2([], _Count, yes) --> [].
-output_equations_2([ X | Xs ], Count, Succ) --> 
-	output_eqn(X, Count, Succ0),
-	( { Succ0 = yes } ->
-		{ Count1 is Count + 1 },
-		output_equations_2(Xs, Count1, Succ)
 	;
-		{ Succ = Succ0 }
+		{ Error = Context - lpsolve_failed },
+		{ set_pred_proc_ids_const(PPIds, inf(Error), Module0, Module) }
 	).
 
-:- pred output_eqn(term_pass1__equation, int, bool, io__state,
-	io__state).
-:- mode output_eqn(in, in, out, di, uo) is det.
-
-% each constraint is of the form:
-% c#: # b + PPId - (PPIdList) > Const;
-% each PPId is printed as `aPredId_ProcId'
-% As each PPId can be negative, and lp_solve allows only positive
-% variables, we introduce a dummy variable 'b' such that 
-% Actual PPId value = returned PPId value - b, where the returned value is 
-% always non-negative
-output_eqn(Eqn, Count, Succ) -->
-	{ Eqn = eqn(Const, PPId, PPIdList0) },
-	{ list__length(PPIdList0, Length) },
-	% As there are `Length' negative PPIds, and 1 positive PPId, and
-	% each PPId is replaced with `PPId - b', the multiplier of b is
-	% Length - 1.
-	{ BMultiplier is Length - 1 },
-	( { list__delete_first(PPIdList0, PPId, PPIdList1) } ->
-		{ Deleted = yes },
-		{ PPIdList = PPIdList1 }
+:- pred convert_equations(list(term_pass1__equation), varset,
+		lp__equations, objective, map(pred_proc_id, var)).
+:- mode convert_equations(in, out, out, out, out) is semidet.
+
+convert_equations(TermEquations, Varset, LPEquations, Objective, PPVars) :-
+	varset__init(Varset0),
+	map__init(PredProcVars0),
+	convert_equations2(TermEquations, PredProcVars0, PPVars,
+		Varset0, Varset, [], LPEquations),
+	map__values(PPVars, Vars),
+	Convert = lambda([Var::in, Coeff::out] is det,
+	(
+		Coeff = Var - 1.0
+	)),
+	list__map(Convert, Vars, Objective).
+
+:- pred convert_equations2(list(term_pass1__equation),
+		map(pred_proc_id, var), map(pred_proc_id, var), varset, varset,
+		lp__equations, lp__equations).
+:- mode convert_equations2(in, in, out, in, out, in, out) is semidet.
+
+convert_equations2([], PPVars, PPVars, Varset, Varset, LPEqns, LPEqns).
+convert_equations2([Eqn|Eqns], PPVars0, PPVars, Varset0, Varset,
+		LPEqns0, LPEqns) :-
+	Eqn = eqn(set(IntConst), ThisPPId, PPIds),
+	int__to_float(IntConst, FloatConst),
+	LPEqn = eqn(Coeffs, (>=), FloatConst),
+	pred_proc_var(ThisPPId, ThisVar, Varset0, Varset2, PPVars0, PPVars1),
+	Coeffs = [ThisVar - 1.0|RestCoeffs],
+	Convert = lambda([PPId::in, Coeff::out, Pair0::in, Pair::out] is det,
+	(
+		Pair0 = VS0 - PPV0,
+		pred_proc_var(PPId, Var, VS0, VS, PPV0, PPV),
+		Coeff = Var - (-1.0),
+		Pair = VS - PPV
+	)),
+	list__map_foldl(Convert, PPIds, RestCoeffs, Varset2 - PPVars1,
+		Varset3 - PPVars2),
+	convert_equations2(Eqns, PPVars2, PPVars, Varset3, Varset,
+						[LPEqn|LPEqns0], LPEqns).
+
+:- pred lookup_coeff(map(pred_proc_id, var), map(var, float), pred_proc_id,
+		pair(pred_proc_id, int)).
+:- mode lookup_coeff(in, in, in, out) is det.
+
+lookup_coeff(PPIds, Soln, PPId, PPId - ICoeff) :-
+	map__lookup(PPIds, PPId, Var),
+	map__lookup(Soln, Var, Coeff),
+	ceiling_to_int(Coeff, ICoeff).
+
+:- pred pred_proc_var(pred_proc_id, var, varset, varset,
+		map(pred_proc_id, var), map(pred_proc_id, var)).
+:- mode pred_proc_var(in, out, in, out, in, out) is det.
+
+pred_proc_var(PPId, Var, Varset0, Varset, PPVars0, PPVars) :-
+	( map__search(PPVars0, PPId, Var0) ->
+		Var = Var0,
+		Varset = Varset0,
+		PPVars = PPVars0
 	;
-		{ Deleted = no },
-		{ PPIdList = PPIdList0 }
-	),
-
-	( 
-		{ Length = 1 },
-		{ Deleted = yes }
-	->
-		% There is nothing on the left hand side  of the
-		% constraint, as there was only one element in the list,
-		% and it was deleted.  Therefore the left hand side of the
-		% equation is PPId - PPId which lpsolve does not process.
-		% Constraints of this type should all be removed by 
-		% proc_inequalities_scc_remove_useless_offsets.
-		{ Succ = no }
-	;
-		% output 'c#: '
-		io__write_char('c'),
-		io__write_int(Count),
-		io__write_string(": "),
-	
-		% maybe output '# b '
-		( { BMultiplier = 0 } ->
-			[]
-		;	
-			io__write_int(BMultiplier),
-			io__write_string(" b ")
-		),
-		
-		% maybe output ' + PPId'
-		( { Deleted = yes } ->
-			[]
-		;
-			io__write_string(" + a"),
-			output_eqn_ppid(PPId)
-		),
-
-		% output `PPIdList' 
-		output_eqn_2(PPIdList),
-
-		% output ` > Const;'
-		io__write_string(" > "),
-		( { Const = set(Int) } ->
-			io__write_int(Int), 
-			{ Succ = yes }
-		;
-			{ Succ = no }
-		),
-		io__write_string(";\n")
+		varset__new_var(Varset0, Var, Varset),
+		map__set(PPVars0, PPId, Var, PPVars)
 	).
-	
-% Outputs each of the pred proc id's in the form: ' - aPredId_ProcId'
-:- pred output_eqn_2(list(pred_proc_id), io__state, io__state).
-:- mode output_eqn_2(in, di, uo) is det.
-
-output_eqn_2([]) --> [].
-output_eqn_2([ X | Xs ]) --> 
-	io__write_string(" - a"),
-	output_eqn_ppid(X),
-	output_eqn_2(Xs).
-
-% Outputs a pred proc id in the form "PredId_ProcId"
-:- pred output_eqn_ppid(pred_proc_id, io__state, io__state).
-:- mode output_eqn_ppid(in, di, uo) is det.
-output_eqn_ppid(proc(PredId, ProcId)) -->
-	{ pred_id_to_int(PredId, PredInt) },
-	{ proc_id_to_int(ProcId, ProcInt) },
-	io__write_int(PredInt),
-	io__write_char('_'),
-	io__write_int(ProcInt).
 
cvs diff: Diffing notes
%-----------------------------------------------------------------------------%
% Copyright (C) 1994-1995, 1997 The University of Melbourne.
% This file may only be copied under the terms of the GNU General
% Public License - see the file COPYING in the Mercury distribution.
%-----------------------------------------------------------------------------%
%
% file: lp.m
% main author: conway,	Oct 1997
%
% This module implements a linear constraint solver that finds an
% optimal solution to a set of linear [in]equalities with respect
% to some objective function. It does this using the simplex method.
%
% The form of an [in]equation is
% 	a1.x1 + a2.x2 + ... + an.xn {=<,=,>=} b
% where alll the numbers are floats.
% By defaut, there is an additional constraint on each of the `xn's:
%	xn >= 0
% If you want xn to take on any value, you can include it in the list
% of URS (UnRestricted in Sign) variables.
%
% The objective function is simply a weighted sum of the variables.
%
% The `x's are represented by `term__var's. The varset from which
% they are allocated is passed to the solver because it needs to
% introduce new variables as part of the solving algorithm.
%
%------------------------------------------------------------------------------%
:- module lp.

:- interface.

%------------------------------------------------------------------------------%

:- import_module float, io, list, map, std_util, term, varset.

:- type coeff	==	pair(var, float).

:- type equation
	--->	eqn(list(coeff), operator, float).

:- type operator
	--->	(=<) ; (=) ; (>=) .

:- type equations	==	list(equation).

:- type objective	==	list(coeff).

:- type direction
	--->	max ; min .

:- type lp__result
	--->	unsatisfiable
	;	satisfiable(float, map(var, float))
	.

%------------------------------------------------------------------------------%

	% lp_solve(Inequations, MaxOrMin, Objective, Varset, URSVars,
	%		Result, IO0, IO)
	% maximize (or minimize - depending on `MaxOrMin') `Objective'
	% subject to the constraints `Inequations'. The variables in
	% the objective and inequations are from `Varset' which is passed
	% so that the solver can allocate fresh variables as required.
	% URSVars is the list of variable that are unrestricted in sign.
	% lp_solve binds Result either to `unsatisfiable' if the there
	% was no optimum value of the objective function (ie the
	% constraints were inconsistent, or the objective function
	% is unbounded by the constraints), or `satisfiable(ObjVal,
	% MapFromObjVarsToVals)'.

:- pred lp_solve(equations, direction, objective, varset, list(var),
		lp__result, io__state, io__state).
:- mode lp_solve(in, in, in, in, in, out, di, uo) is det.

%------------------------------------------------------------------------------%
%------------------------------------------------------------------------------%

:- implementation.

:- import_module bool, int, require, set.

:- type lp_info
	---> lp(
		varset,
		map(var, pair(var)),	% variables with urs.
		list(var),		% slack variables
		list(var)		% artificial variables
	).

lp_solve(Eqns, Dir, Obj, Varset0, URSVars, Result, IO0, IO) :-
	lp_info_init(Varset0, URSVars, Info0),
	lp_solve2(Eqns, Dir, Obj, Result, IO0, IO, Info0, _).

:- pred lp_solve2(equations, direction, objective, lp__result,
		io__state, io__state, lp_info, lp_info).
:- mode lp_solve2(in, in, in, out, di, uo, in, out) is det.

lp_solve2(Eqns0, Dir, Obj0, Result, IO0, IO, Info0, Info) :-
	standardize_equations(Eqns0, Eqns, Info0, Info1),
	(
		Dir = max,
		negate_equation(eqn(Obj0, (=), 0.0), eqn(Obj1, _, _))
	;
		Dir = min,
		Obj1 = Obj0
	),
	simplify_coeffs(Obj1, Obj2),
	get_urs_vars(URS, Info1, _),
	expand_urs_vars(Obj2, URS, Obj),
	list__length(Eqns, Rows),
	collect_vars(Eqns, Obj, Vars),
	set__to_sorted_list(Vars, VarList),
	list__length(VarList, Cols),
	map__init(VarNums0),
	number_vars(VarList, 0, VarNums0, VarNums),
	get_art_vars(ArtVars, Info1, Info2),
	get_urs_vars(URSVars, Info2, Info),
	init_tableau(Rows, Cols, VarNums, URSVars, Tableau0),
	insert_equations(Eqns, 1, Cols, VarNums,
		Tableau0, Tableau1),
	(
		ArtVars = [_|_],
		% There are one or more artificial variables, so we use
		% the two-phase method for solving the system.
		two_phase(Obj0, Obj, Dir, ArtVars, VarNums, Tableau1, Result,
			IO0, IO)
	;
		ArtVars = [],
		one_phase(Obj0, Obj, Dir, VarNums, Tableau1, Result, IO0, IO)
	).

%------------------------------------------------------------------------------%

:- pred one_phase(list(coeff), list(coeff), direction, map(var, int), tableau,
		lp__result, io__state, io__state).
:- mode one_phase(in, in, in, in, in, out, di, uo) is det.

one_phase(Obj0, Obj, Dir, VarNums, Tableau0, Result, IO0, IO) :-
	insert_coeffs(Obj, 0, VarNums, Tableau0, Tableau1),
	GetObjVar = lambda([V::out] is nondet, (
		list__member(X, Obj0),
		X = V - _Cof
	)),
	solutions(GetObjVar, ObjVars),
	optimize(ObjVars, Tableau1, _, Result0, IO0, IO),
	(
		Dir = max,
		Result = Result0
	;
		Dir = min,
		(
			Result0 = unsatisfiable,
			Result = Result0
		;
			Result0 = satisfiable(NOptVal, OptCoffs),
			OptVal is -NOptVal,
			Result = satisfiable(OptVal, OptCoffs)
		)
	).

%------------------------------------------------------------------------------%

:- pred two_phase(list(coeff), list(coeff), direction, list(var), map(var, int),
		tableau, lp__result, io__state, io__state).
:- mode two_phase(in, in, in, in, in, in, out, di, uo) is det.

two_phase(Obj0, Obj, Dir, ArtVars, VarNums, Tableau0, Result, IO0, IO) :-
		% Do phase 1:
		%	minimize the sum of the artificial variables
	construct_art_objective(ArtVars, ArtObj),
	insert_coeffs(ArtObj, 0, VarNums, Tableau0, Tableau1a),
	ensure_zero_obj_coeffs(ArtVars, Tableau1a, Tableau1b),
	optimize(ArtVars, Tableau1b, Tableau1c,
		Res0, IO0, IO1),
	(
		Res0 = unsatisfiable,
		Result = unsatisfiable,
		IO = IO1
	;
		Res0 = satisfiable(Val, _ArtRes),
		( Val \= 0.0 ->
			Result = unsatisfiable,
			IO = IO1
		;
			fix_basis_and_rem_cols(ArtVars, Tableau1c, Tableau2),
				% Do phase 2:
				%	insert the real objective,
				%	zero the objective coefficients of the
				%	basis variables,
				%	optimize the objective.
			insert_coeffs(Obj, 0, VarNums, Tableau2, Tableau3),
			get_basis_vars(Tableau3, BasisVars),
			ensure_zero_obj_coeffs(BasisVars,
					Tableau3, Tableau4),
			GetObjVar = lambda([V::out] is nondet, (
				list__member(X, Obj0),
				X = V - _Cof
			)),
			solutions(GetObjVar, ObjVars),
			optimize(ObjVars, Tableau4, _, Res1, IO1, IO),
			(
				Dir = max,
				Result = Res1
			;
				Dir = min,
				(
					Res1 = unsatisfiable,
					Result = Res1
				;
					Res1 = satisfiable(NOptVal, OptCoffs),
					OptVal is -NOptVal,
					Result = satisfiable(OptVal, OptCoffs)
				)
			)
		)
	).

%------------------------------------------------------------------------------%

:- pred construct_art_objective(list(var), list(coeff)).
:- mode construct_art_objective(in, out) is det.

construct_art_objective([], []).
construct_art_objective([V|Vs], [V - (1.0)|Rest]) :-
	construct_art_objective(Vs, Rest).

%------------------------------------------------------------------------------%

:- pred standardize_equations(equations, equations, lp_info, lp_info).
:- mode standardize_equations(in, out, in, out) is det.

standardize_equations(Eqns0, Eqns) -->
	list__map_foldl(standardize_equation, Eqns0, Eqns).

:- pred standardize_equation(equation, equation, lp_info, lp_info).
:- mode standardize_equation(in, out, in, out) is det.

standardize_equation(Eqn0, Eqn) -->
	{ Eqn0 = eqn(Coeffs0, (=<), Const0) },
	( { Const0 < 0.0 } ->
		{ negate_equation(Eqn0, Eqn1) },
		standardize_equation(Eqn1, Eqn) 
	;
		new_slack_var(Var),
		{ Coeffs = [Var - 1.0|Coeffs0] },
		{ simplify(eqn(Coeffs, (=<), Const0), Eqn1) },
		get_urs_vars(URS),
		{ expand_urs_vars_e(Eqn1, URS, Eqn) }
	).

standardize_equation(Eqn0, Eqn) -->
	{ Eqn0 = eqn(Coeffs0, (=), Const0) },
	( { Const0 < 0.0 } ->
		{ negate_equation(Eqn0, Eqn1) },
		standardize_equation(Eqn1, Eqn)
	;
		new_art_var(Var),
		{ Coeffs = [Var - 1.0|Coeffs0] },
		{ simplify(eqn(Coeffs, (=<), Const0), Eqn1) },
		get_urs_vars(URS),
		{ expand_urs_vars_e(Eqn1, URS, Eqn) }
	).

standardize_equation(Eqn0, Eqn) -->
	{ Eqn0 = eqn(Coeffs0, (>=), Const0) },
	( { Const0 < 0.0 } ->
		{ negate_equation(Eqn0, Eqn1) },
		standardize_equation(Eqn1, Eqn)
	;
		new_slack_var(SVar),
		new_art_var(AVar),
		{ Coeffs = [SVar - (-1.0), AVar - (1.0)|Coeffs0] },
		{ simplify(eqn(Coeffs, (>=), Const0), Eqn1) },
		get_urs_vars(URS),
		{ expand_urs_vars_e(Eqn1, URS, Eqn) }
	).

:- pred negate_equation(equation, equation).
:- mode negate_equation(in, out) is det.

negate_equation(eqn(Coeffs0, Op0, Const0), eqn(Coeffs, Op, Const)) :-
	(
		Op0 = (=<), Op = (>=)
	;
		Op0 = (=), Op = (=)
	;
		Op0 = (>=), Op = (=<)
	),
	Neg = lambda([Pair0::in, Pair::out] is det, (
		Pair0 = V - X0,
		X is -X0,
		Pair = V - X
	)),
	list__map(Neg, Coeffs0, Coeffs),
	Const is -Const0.

:- pred simplify(equation, equation).
:- mode simplify(in, out) is det.

simplify(eqn(Coeffs0, Op, Const), eqn(Coeffs, Op, Const)) :-
	simplify_coeffs(Coeffs0, Coeffs).

:- pred simplify_coeffs(list(coeff), list(coeff)).
:- mode simplify_coeffs(in, out) is det.

simplify_coeffs(Coeffs0, Coeffs) :-
	map__init(CoeffMap0),
	AddCoeff = lambda([Pair::in, Map0::in, Map::out] is det, (
		Pair = Var - Coeff,
		add_var(Map0, Var, Coeff, Map)
	)),
	list__foldl(AddCoeff, Coeffs0, CoeffMap0, CoeffMap),
	map__to_assoc_list(CoeffMap, Coeffs).

:- pred add_var(map(var, float), var, float, map(var, float)).
:- mode add_var(in, in, in, out) is det.

add_var(Map0, Var, Coeff, Map) :-
	( map__search(Map0, Var, Acc0) ->
		Acc1 = Acc0
	;
		Acc1 = 0.0
	),
	Acc is Acc1 + Coeff,
	map__set(Map0, Var, Acc, Map).

:- pred expand_urs_vars_e(equation, map(var, pair(var)), equation).
:- mode expand_urs_vars_e(in, in, out) is det.

expand_urs_vars_e(eqn(Coeffs0, Op, Const), Vars, eqn(Coeffs, Op, Const)) :-
	expand_urs_vars(Coeffs0, Vars, Coeffs).

:- pred expand_urs_vars(list(coeff), map(var, pair(var)), list(coeff)).
:- mode expand_urs_vars(in, in, out) is det.

expand_urs_vars(Coeffs0, Vars, Coeffs) :-
	expand_urs_vars(Coeffs0, Vars, [], Coeffs1),
	list__reverse(Coeffs1, Coeffs).

:- pred expand_urs_vars(list(coeff), map(var, pair(var)),
		list(coeff), list(coeff)).
:- mode expand_urs_vars(in, in, in, out) is det.

expand_urs_vars([], _Vars, Coeffs, Coeffs).
expand_urs_vars([Var - Coeff|Rest], Vars, Coeffs0, Coeffs) :-
	( map__search(Vars, Var, PVar - NVar) ->
		NCoeff is -Coeff,
		Coeffs1 = [NVar - NCoeff, PVar - Coeff|Coeffs0]
	;
		Coeffs1 = [Var - Coeff|Coeffs0]
	),
	expand_urs_vars(Rest, Vars, Coeffs1, Coeffs).

%------------------------------------------------------------------------------%

:- pred collect_vars(equations, objective, set(var)).
:- mode collect_vars(in, in, out) is det.

collect_vars(Eqns, Obj, Vars) :-
	GetVar = lambda([Var::out] is nondet, (
		(
			list__member(Eqn, Eqns),
			Eqn = eqn(Coeffs, _, _),
			list__member(Pair, Coeffs),
			Pair = Var - _
		;
			list__member(Pair, Obj),
			Pair = Var - _
		)
	)),
	solutions(GetVar, VarList),
	set__list_to_set(VarList, Vars).

:- pred number_vars(list(var), int, map(var, int), map(var, int)).
:- mode number_vars(in, in, in, out) is det.

number_vars([], _, VarNums, VarNums).
number_vars([Var|Vars], N, VarNums0, VarNums) :-
	map__det_insert(VarNums0, Var, N, VarNums1),
	N1 is N + 1,
	number_vars(Vars, N1, VarNums1, VarNums).

:- pred insert_equations(equations, int, int, map(var, int), tableau, tableau).
:- mode insert_equations(in, in, in, in, in, out) is det.

insert_equations([], _, _, _, Tableau, Tableau).
insert_equations([Eqn|Eqns], Row, ConstCol, VarNums, Tableau0, Tableau) :-
	Eqn = eqn(Coeffs, _Op, Const),
	insert_coeffs(Coeffs, Row, VarNums, Tableau0, Tableau1),
	set_index(Tableau1, Row, ConstCol, Const, Tableau2),
	Row1 is Row + 1,
	insert_equations(Eqns, Row1, ConstCol, VarNums, Tableau2, Tableau).

:- pred insert_coeffs(list(coeff), int, map(var, int), tableau, tableau).
:- mode insert_coeffs(in, in, in, in, out) is det.

insert_coeffs([], _Row, _VarNums, Tableau, Tableau).
insert_coeffs([Coeff|Coeffs], Row, VarNums, Tableau0, Tableau) :-
	Coeff = Var - Const,
	map__lookup(VarNums, Var, Col),
	set_index(Tableau0, Row, Col, Const, Tableau1),
	insert_coeffs(Coeffs, Row, VarNums, Tableau1, Tableau).

%------------------------------------------------------------------------------%

:- pred optimize(list(var), tableau, tableau, lp__result,
		io__state, io__state).
:- mode optimize(in, in, out, out, di, uo) is det.

optimize(ObjVars, A0, A, Result) -->
	simplex(A0, A, Res0),
	(
		{ Res0 = no },
		{ Result = unsatisfiable }
	;
		{ Res0 = yes },
		{ rhs_col(A, M) },
		{ index(A, 0, M, ObjVal) },
		{ extract_objective(ObjVars, A, ObjMap) },
		{ Result = satisfiable(ObjVal, ObjMap) }
	).

:- pred extract_objective(list(var), tableau, map(var, float)).
:- mode extract_objective(in, in, out) is det.

extract_objective(ObjVars, Tab, Res) :-
	map__init(Res0),
	list__foldl(extract_obj_var(Tab), ObjVars, Res0, Res).

:- pred extract_obj_var(tableau, var, map(var, float), map(var, float)).
:- mode extract_obj_var(in, in, in, out) is det.

extract_obj_var(Tab, Var, Map0, Map) :-
	urs_vars(Tab, Vars),
	( map__search(Vars, Var, Pos - Neg) ->
		extract_obj_var2(Tab, Pos, PosVal),
		extract_obj_var2(Tab, Neg, NegVal),
		Val is PosVal - NegVal
	;
		extract_obj_var2(Tab, Var, Val)
	),
	map__set(Map0, Var, Val, Map).

:- pred extract_obj_var2(tableau, var, float).
:- mode extract_obj_var2(in, in, out) is det.

extract_obj_var2(Tab, Var, Val) :-
	var_col(Tab, Var, Col),
	GetCell = lambda([Val0::out] is nondet, (
		all_rows(Tab, Row),
		index(Tab, Row, Col, 1.0),
		rhs_col(Tab, RHS),
		index(Tab, Row, RHS, Val0)
	)),
	solutions(GetCell, Solns),
	( Solns = [Val1] ->
		Val = Val1
	;
		Val = 0.0
	).

:- pred simplex(tableau, tableau, bool, io__state, io__state).
:- mode simplex(in, out, out, di, uo) is det.

simplex(A0, A, Result, IO0, IO) :-
	AllColumns = all_cols(A0),
	MinAgg = lambda([Col::in, Min0::in, Min::out] is det, (
		(
			Min0 = no,
			index(A0, 0, Col, MinVal),
			( MinVal < 0.0 ->
				Min = yes(Col - MinVal)
			;
				Min = no
			)
		;
			Min0 = yes(_ - MinVal0),
			index(A0, 0, Col, CellVal),
			( CellVal < MinVal0 ->
				Min = yes(Col - CellVal)
			;
				Min = Min0
			)
		)
	)),
	aggregate(AllColumns, MinAgg, no, MinResult),
	(
		MinResult = no,
		A = A0,
		IO = IO0,
		Result = yes
	;
		MinResult = yes(Q - _Val),
		AllRows = all_rows(A0),
		MaxAgg = lambda([Row::in, Max0::in, Max::out] is det, (
			(
				Max0 = no,
				index(A0, Row, Q, MaxVal),
				( MaxVal > 0.0 ->
					rhs_col(A0, RHSC),
					index(A0, Row, RHSC, MVal),
					CVal is MVal/MaxVal,
					Max = yes(Row - CVal)
				;
					Max = no
				)
			;
				Max0 = yes(_ - MaxVal0),
				index(A0, Row, Q, CellVal),
				rhs_col(A0, RHSC),
				index(A0, Row, RHSC, MVal),
				MaxVal1 is MVal/CellVal,
				( CellVal > 0.0, MaxVal1 =< MaxVal0 ->
					Max = yes(Row - MaxVal1)
				;
					Max = Max0
				)
			)
		)),
		aggregate(AllRows, MaxAgg, no, MaxResult),
		(
			MaxResult = no,
			A = A0,
			IO = IO0,
			Result = no
		;
			MaxResult = yes(P - _),
			pivot(P, Q, A0, A1),
			simplex(A1, A, Result, IO0, IO)
		)
	).

%------------------------------------------------------------------------------%

:- pred ensure_zero_obj_coeffs(list(var), tableau, tableau).
:- mode ensure_zero_obj_coeffs(in, in, out) is det.

ensure_zero_obj_coeffs([], Tableau, Tableau).
ensure_zero_obj_coeffs([V|Vs], Tableau0, Tableau) :-
	var_col(Tableau0, V, Col),
	index(Tableau0, 0, Col, Val),
	( Val = 0.0 ->
		ensure_zero_obj_coeffs(Vs, Tableau0, Tableau)
	;
		FindOne = lambda([P::out] is nondet, (
			all_rows(Tableau0, R),
			index(Tableau0, R, Col, ValF0),
			ValF0 \= 0.0,
			P = R - ValF0
		)),
		solutions(FindOne, Ones),
		(
			Ones = [Row - Fac0|_],
			Fac is -Val/Fac0,
			row_op(Fac, Row, 0, Tableau0, Tableau1),
			ensure_zero_obj_coeffs(Vs, Tableau1, Tableau)
		;
			Ones = [],
			error("problem with artificial variable")
		)
	).

:- pred fix_basis_and_rem_cols(list(var), tableau, tableau).
:- mode fix_basis_and_rem_cols(in, in, out) is det.

fix_basis_and_rem_cols([], Tab, Tab).
fix_basis_and_rem_cols([V|Vs], Tab0, Tab) :-
	var_col(Tab0, V, Col),
	BasisAgg = lambda([R::in, Ones0::in, Ones::out] is det, (
		index(Tab0, R, Col, Val),
		( Val = 0.0 ->
			Ones = Ones0
		;
			Ones = [Val - R|Ones0]
		)
	)),
	aggregate(all_rows(Tab0), BasisAgg, [], Res),
	(
		Res = [1.0 - Row]
	->
		PivGoal = lambda([Col1::out] is nondet, (
			all_cols(Tab0, Col1),
			Col \= Col1,
			index(Tab0, Row, Col1, Zz),
			Zz \= 0.0
		)),
		solutions(PivGoal, PivSolns),
		(
			PivSolns = [],
			remove_col(Col, Tab0, Tab0a),
			remove_row(Row, Tab0a, Tab1)
		;
			PivSolns = [Col2|_],
			pivot(Row, Col2, Tab0, Tab0a),
			remove_col(Col, Tab0a, Tab1)
		)
	;
		Tab1 = Tab0
	),
	remove_col(Col, Tab1, Tab2),
	fix_basis_and_rem_cols(Vs, Tab2, Tab).

%------------------------------------------------------------------------------%

:- type cell	--->	cell(int, int).

:- pred pivot(int, int, tableau, tableau).
:- mode pivot(in, in, in, out) is det.

pivot(P, Q, A0, A) :-
	index(A0, P, Q, Apq),
	MostCells = lambda([Cell::out] is nondet, (
		all_rows0(A0, J),
		J \= P,
		all_cols0(A0, K),
		K \= Q,
		Cell = cell(J, K)
	)),
	ScaleCell = lambda([Cell::in, T0::in, T::out] is det, (
		Cell = cell(J, K),
		index(T0, J, K, Ajk),
		index(T0, J, Q, Ajq),
		index(T0, P, K, Apk),
		NewAjk is Ajk - Apk * Ajq / Apq,
		set_index(T0, J, K, NewAjk, T)
	)),
	aggregate(MostCells, ScaleCell, A0, A1),
	QColumn = lambda([Cell::out] is nondet, (
		all_rows0(A1, J),
		Cell = cell(J, Q)
	)),
	Zero = lambda([Cell::in, T0::in, T::out] is det, (
		Cell = cell(J, K),
		set_index(T0, J, K, 0.0, T)
	)),
	aggregate(QColumn, Zero, A1, A2),
	PRow = all_cols0(A2),
	ScaleRow = lambda([K::in, T0::in, T::out] is det, (
		index(T0, P, K, Apk),
		NewApk is Apk / Apq,
		set_index(T0, P, K, NewApk, T)
	)),
	aggregate(PRow, ScaleRow, A2, A3),
	set_index(A3, P, Q, 1.0, A).

:- pred row_op(float, int, int, tableau, tableau).
:- mode row_op(in, in, in, in, out) is det.

row_op(Scale, From, To, A0, A) :-
	AllCols = all_cols0(A0),
	AddRow = lambda([Col::in, T0::in, T::out] is det, (
		index(T0, From, Col, X),
		index(T0, To, Col, Y),
		Z is Y + (Scale * X),
		set_index(T0, To, Col, Z, T)
	)),
	aggregate(AllCols, AddRow, A0, A).

%------------------------------------------------------------------------------%

:- type tableau
	---> tableau(
		int,
		int,
		map(var, int),
		map(var, pair(var)),
		list(int),	% shunned rows
		list(int),	% shunned cols
		map(pair(int), float)
	).

:- pred init_tableau(int::in, int::in, map(var, int)::in, 
		map(var, pair(var))::in, tableau::out) is det.

init_tableau(Rows, Cols, VarNums, URSVars, Tableau) :-
	map__init(Cells),
	Tableau = tableau(Rows, Cols, VarNums, URSVars, [], [], Cells).

:- pred index(tableau, int, int, float).
:- mode index(in, in, in, out) is det.

index(Tableau, J, K, R) :-
	Tableau = tableau(_, _, _, _, SR, SC, Cells),
	(
		( list__member(J, SR)
		; list__member(K, SC)
		)
	->
		error("attempt to address shunned row/col")
	;
		true
	),
	(
		map__search(Cells, J - K, R0)
	->
		R = R0
	;
		R = 0.0
	).

:- pred set_index(tableau, int, int, float, tableau).
:- mode set_index(in, in, in, in, out) is det.

set_index(Tableau0, J, K, R, Tableau) :-
	Tableau0 = tableau(Rows, Cols, VarNums, URS, SR, SC, Cells0),
	(
		( list__member(J, SR)
		; list__member(K, SC)
		)
	->
		error("attempt to write shunned row/col")
	;
		true
	),
	( R = 0.0 ->
		map__delete(Cells0, J - K, Cells)
	;
		map__set(Cells0, J - K, R, Cells)
	),
	Tableau = tableau(Rows, Cols, VarNums, URS, SR, SC, Cells).

:- pred rhs_col(tableau, int).
:- mode rhs_col(in, out) is det.

rhs_col(tableau(_, RHS, _, _, _, _, _), RHS).

:- pred all_rows0(tableau, int).
:- mode all_rows0(in, out) is nondet.

all_rows0(Tableau, Row) :-
	Tableau = tableau(Rows, _Cols, _, _, SR, _, _),
	between(0, Rows, Row),
	\+ list__member(Row, SR).

:- pred all_rows(tableau, int).
:- mode all_rows(in, out) is nondet.

all_rows(Tableau, Row) :-
	Tableau = tableau(Rows, _Cols, _, _, SR, _, _),
	between(1, Rows, Row),
	\+ list__member(Row, SR).

:- pred all_cols0(tableau, int).
:- mode all_cols0(in, out) is nondet.

all_cols0(Tableau, Col) :-
	Tableau = tableau(_Rows, Cols, _, _, _, SC, _),
	between(0, Cols, Col),
	\+ list__member(Col, SC).

:- pred all_cols(tableau, int).
:- mode all_cols(in, out) is nondet.

all_cols(Tableau, Col) :-
	Tableau = tableau(_Rows, Cols, _, _, _, SC, _),
	Cols1 is Cols - 1,
	between(0, Cols1, Col),
	\+ list__member(Col, SC).

:- pred var_col(tableau, var, int).
:- mode var_col(in, in, out) is det.

var_col(Tableau, Var, Col) :-
	Tableau = tableau(_, _, VarCols, _, _, _, _),
	map__lookup(VarCols, Var, Col).

:- pred urs_vars(tableau, map(var, pair(var))).
:- mode urs_vars(in, out) is det.

urs_vars(Tableau, URS) :-
	Tableau = tableau(_, _, _, URS, _, _, _).

:- pred remove_row(int, tableau, tableau).
:- mode remove_row(in, in, out) is det.

remove_row(R, Tableau0, Tableau) :-
	Tableau0 = tableau(Rows, Cols, VarNums, URS, SR, SC, Cells),
	Tableau = tableau(Rows, Cols, VarNums, URS, [R|SR], SC, Cells).

:- pred remove_col(int, tableau, tableau).
:- mode remove_col(in, in, out) is det.

remove_col(C, Tableau0, Tableau) :-
	Tableau0 = tableau(Rows, Cols, VarNums, URS, SR, SC, Cells),
	Tableau = tableau(Rows, Cols, VarNums, URS, SR, [C|SC], Cells).

:- pred get_basis_vars(tableau, list(var)).
:- mode get_basis_vars(in, out) is det.

get_basis_vars(Tab, Vars) :-
	BasisCol = lambda([C::out] is nondet, (
		all_cols(Tab, C),
		NonZeroGoal = lambda([P::out] is nondet, (
			all_rows(Tab, R),
			index(Tab, R, C, Z),
			Z \= 0.0,
			P = R - Z
		)),
		solutions(NonZeroGoal, Solns),
		Solns = [_ - 1.0]
	)),
	solutions(BasisCol, Cols),
	BasisVars = lambda([V::out] is nondet, (
		list__member(Col, Cols),
		Tab = tableau(_, _, VarCols, _, _, _, _),
		map__member(VarCols, V, Col)
	)),
	solutions(BasisVars, Vars).

%------------------------------------------------------------------------------%

:- import_module string.

:- pred show_tableau(tableau, io__state, io__state).
:- mode show_tableau(in, di, uo) is det.

/*
show_tableau(_N, _M, _Tableau) --> [].
*/

show_tableau(Tableau) -->
	{ Tableau = tableau(N, M, _, _, _, _, _) },
	{ string__format("Tableau (%d, %d):\n", [i(N), i(M)], Str) },
	io__write_string(Str),
	aggregate(all_rows0(Tableau), show_row(Tableau)).

:- pred show_row(tableau, int, io__state, io__state).
:- mode show_row(in, in, di, uo) is det.

show_row(Tableau, Row) -->
	aggregate(all_cols0(Tableau), show_cell(Tableau, Row)),
	io__write_string("\n").

:- pred show_cell(tableau, int, int, io__state, io__state).
:- mode show_cell(in, in, in, di, uo) is det.

show_cell(Tableau, Row, Col) -->
	{ index(Tableau, Row, Col, Val) },
	{ string__format("%2.2f\t", [f(Val)], Str) },
	io__write_string(Str).

%------------------------------------------------------------------------------%

:- pred lp_info_init(varset, list(var), lp_info).
:- mode lp_info_init(in, in, out) is det.

lp_info_init(Varset0, URSVars, LPInfo) :-
	Introduce = lambda([Orig::in, VP0::in, VP::out] is det, (
		VP0 = VS0 - VM0,
		varset__new_var(VS0, V1, VS1),
		varset__new_var(VS1, V2, VS),
		map__set(VM0, Orig, V1 - V2, VM),
		VP = VS - VM
	)),
	map__init(URSMap0),
	list__foldl(Introduce, URSVars, Varset0 - URSMap0, Varset - URSMap),
	LPInfo = lp(Varset, URSMap, [], []).

:- pred new_slack_var(var::out, lp_info::in, lp_info::out) is det.

new_slack_var(Var) -->
	get_varset(Varset0),
	{ varset__new_var(Varset0, Var, Varset) },
	set_varset(Varset),
	get_slack_vars(Vars),
	set_slack_vars([Var|Vars]).

:- pred new_art_var(var::out, lp_info::in, lp_info::out) is det.

new_art_var(Var) -->
	get_varset(Varset0),
	{ varset__new_var(Varset0, Var, Varset) },
	set_varset(Varset),
	get_art_vars(Vars),
	set_art_vars([Var|Vars]).

:- pred get_varset(varset::out, lp_info::in, lp_info::out) is det.

get_varset(Varset, Info, Info) :-
	Info = lp(Varset, _URSVars, _Slack, _Art).

:- pred set_varset(varset::in, lp_info::in, lp_info::out) is det.

set_varset(Varset, Info0, Info) :-
	Info0 = lp(_Varset, URSVars, Slack, Art),
	Info  = lp(Varset, URSVars, Slack, Art).

:- pred get_urs_vars(map(var, pair(var))::out, lp_info::in, lp_info::out) is det.

get_urs_vars(URSVars, Info, Info) :-
	Info = lp(_Varset, URSVars, _Slack, _Art).

:- pred set_urs_vars(map(var, pair(var))::in, lp_info::in, lp_info::out) is det.

set_urs_vars(URSVars, Info0, Info) :-
	Info0 = lp(Varset, _URSVars, Slack, Art),
	Info  = lp(Varset, URSVars, Slack, Art).

:- pred get_slack_vars(list(var)::out, lp_info::in, lp_info::out) is det.

get_slack_vars(Slack, Info, Info) :-
	Info = lp(_Varset, _URSVars, Slack, _Art).

:- pred set_slack_vars(list(var)::in, lp_info::in, lp_info::out) is det.

set_slack_vars(Slack, Info0, Info) :-
	Info0 = lp(Varset, URSVars, _Slack, Art),
	Info  = lp(Varset, URSVars, Slack, Art).

:- pred get_art_vars(list(var)::out, lp_info::in, lp_info::out) is det.

get_art_vars(Art, Info, Info) :-
	Info = lp(_Varset, _URSVars, _Slack, Art).

:- pred set_art_vars(list(var)::in, lp_info::in, lp_info::out) is det.

set_art_vars(Art, Info0, Info) :-
	Info0 = lp(Varset, URSVars, Slack, _Art),
	Info  = lp(Varset, URSVars, Slack, Art).

%------------------------------------------------------------------------------%

:- pred between(int, int, int).
:- mode between(in, in, out) is nondet.

between(Min, Max, I) :-
	Min =< Max,
	(
		I = Min
	;
		Min1 is Min + 1,
		between(Min1, Max, I)
	).

%------------------------------------------------------------------------------%
%------------------------------------------------------------------------------%