/*==============================================================================
 *	LPAD and CP-Logic reasoning suite
 *	File: montecarlo.pl
 *	Solves LPADs with Monte Carlo (main predicate: solve(Goals, Prob, Samples, ResTime, BddTime)
 *	Copyright (c) 2009, Stefano Bragaglia
 *============================================================================*/
 
/* EXTERNAL FILE
 * -------------
 * The following libraries are required by the program to work fine.
 */

:- dynamic rule/4, def_rule/2, randx/1, randy/1, randz/1.

:- use_module(library(lists)).
:- use_module(library(random)).
:- use_module(library(ugraphs)).
:- use_module(params).

% :- source.
% :- yap_flag(single_var_warnings, on).





/* SOLVING PREDICATES
 * ------------------
 * The predicates in this section solve any given problem with several class of 
 * algorithms.
 *
 * Note: the original predicates (no more need and eligible to be deleted) have 
 *       been moved to the end of the file.
 */

/* Starting seed. */
randx(1).
randy(1).
randz(1).

/* newsample
 * ---------
 * This predicate programmatically generates and sets a new seed for the 
 * algorithm.
 */
newsample :-
	retract(randx(X)),
	randy(Y),
	randz(Z),
	(X =< 30269 ->
		SuccX is X + 1,
		assert(randx(SuccX));
		assert(randx(1)),
		retract(randy(_)),
		(Y =< 30307 ->
			SuccY is Y + 1,
			assert(randy(SuccY));
			assert(randy(1)),
			retract(randz(_)),
			(Z =< 30323 ->
				SuccZ is Z + 1,
				assert(randz(SuccZ));
				assert(randz(1))))),
	setrand(rand(X, Y, Z)).


/* solve(Goals, Samples, Time, Lower, Prob, Upper)
 * -----------------------------------------------
 * This predicate calls the Monte Carlo solving method for the current problem.
 * It requires the Goals to fulfil and returns the number of Samples considered,
 * the Time required, the extimated Probability and the Lower and Upper bounds.
 *
 * INPUT
 *  - Goals: given list of goals to fulfil.
 *
 * OUTPUT
 *  - Samples: number of samples considered to solve the problem.
 *  - Time: time required to solve the problem.
 *  - Lower: lower end of the confidence interval.
 *  - Prob: extimated probability.
 *  - Upper: upper end of the confidence interval.
 */
solve(Goals, Samples, Time, Lower, Prob, Upper) :-
	% Retrieving functional parameters...
	setting(k, K),
	setting(min_error, MinError),
	% Resetting the clocks...
	statistics(walltime, [_, _]),	
	% Performing resolution...
	montecarlo(0, 0, Goals, K, MinError, Samples, Lower, Prob, Upper),	
	% Taking elapsed times...
	statistics(walltime, [_, ElapTime]),
	% Setting values...
	Time is ElapTime/1000.



/* montecarlo(Count, Success, Goals, K, MinError, Samples, Lower, Prob, Upper)
 * ---------------------------------------------------------------------------
 * This tail recursive predicate solves the problem currently in memory with a
 * Monte Carlo approach. 
 * It requires the number of samples and successes so far (Count and Success),
 * the desired list of Goals to fulfil, the number K of samples to consider at
 * once and the threshold MinError for the binomial proportion confidence 
 * interval.
 * It returns the total number of Samples considered, the Lower and Upper ends 
 * of the the binomial proportion confidence interval and the extimated Prob.
 * 
 * INPUT
 *  - Count: number of samples considered so far.
 *  - Success: number of successfull samples considered so far.
 *  - Goals: list of goals to fulfil.
 *  - K: number of samples to consider at once.
 *  - MinError: threshold for the binomial proportion confidence interval.
 * 
 * OUTPUT
 *  - Samples: total number of samples considered.
 *  - Lower: lower end of the the binomial proportion confidence interval.
 *  - Prob: extimated probability.
 *  - Upper: upper end of the the binomial proportion confidence interval.
 *
 * NB: This method is based on the binomial proportion confidence interval and
 *     the Brown's rule of thumb to avoid the case the sample proportion is 
 *     exactly 0.0 or 1.0 and doesn't make use of BDDs.
 */
montecarlo(Count, Success, Goals, K, MinError, Samples, Lower, Prob, Upper) :-
	/* Decomment the following line if you want to test the algorithm with an
	   incremental seed for each sample.
	newsample,
	*/
	main(Goals, [], _Explan, Valid),
	N is Count + 1,
	S is Success + Valid,
	(N mod K =:= 0 ->
%		format("Advancement: ~t~d/~d~30+~n", [S, N]),
		P is S / N,
		D is N - S,
		Semi is 2 * sqrt(P * (1 - P) / N),
		Int is 2 * Semi,
	/*   N * P > 5;   N * S / N > 5;   S > 5
	 *   N (1 - P) > 5;   N (1 - S / N) > 5;   N (N - S) / N > 5;   N - S > 5
	 */
		((S > 5, D > 5, (Int < MinError; Int =:= 0)) ->
			Samples is N,
			Lower is P - Semi,
			Prob is P,
			Upper is P + Semi;
			montecarlo(N, S, Goals, K, MinError, Samples, Lower, Prob, Upper));
		montecarlo(N, S, Goals, K, MinError, Samples, Lower, Prob, Upper)).



/* null
 * ----
 * This is dummy predicate to use sparingly when needed. 
 * Typical uses are as spying predicate during tracing or as dead branch in 
 * ( -> ; ) predicate.
 */	
null.



/* main(Goals, Explan0, Explan1, Valid)
 * ------------------------------------
 * This tail recursive predicate looks for a solution to the given Goals
 * starting from the given Explan0 and returns the final Explan and 1 (0 otherwise) if it is a 
 * Valid sample for Montecarlo.
 */
main([], Explan, Explan, 1). 

main([\+ Goal|Tail], Explan0, Explan1, Valid) :-
	builtin(Goal), !,
	(call((\+ Goal)) ->
		main(Tail, Explan0, Explan1, Valid);
		Explan1 = Explan0,
		Valid = 0).
	
main([Goal|Tail], Explan0, Explan1, Valid) :-
	builtin(Goal), !,
	(call(Goal) ->
		main(Tail, Explan0, Explan1, Valid);
		Explan1 = Explan0,
		Valid = 0).

main([Goal|Tail], Explan0, Explan1, Valid) :-
	findall((IsSample, Goals, Step), explore([Goal|Tail], Explan0, IsSample, Goals, Step), List), 
	cycle(List, Explan0, Explan1, Valid).

	

/* explore([Goal|Tail], Explan, Valid, Goals, Step)
 * ------------------------------------------------
 * This predicate looks for a Body and the Step to reach it from the given Goal 
 * and Explan and returns 1 (0 otherwise) if they are a Valid sample for
 * Montecarlo. 
 * Please note that Body and Step are meaningfull only when Valid is 1.
 *
 * This comment has to be fixed.
 */
explore([Goal|Tail], _Explan, 1, Goals, []) :-
	def_rule(Goal, Body),
	append(Body, Tail, Goals).

explore([Goal|Tail], Explan, Valid, Goals, Step) :- 
	findrule(Goal, Explan, Valid, Body, (HeadID, RuleID, Subst)),
	append(Body, Tail, Goals),
	(member_eq((HeadID, RuleID, Subst), Explan) ->
		Step = [];
		Step = [(HeadID, RuleID, Subst)]).



/* findrule(Goal, Explan, Valid, Body, (HeadID, RuleID, Subst))
 * ---------------------------------------------------------------
 * This predicate finds a rule that matches with the given Goal and Explan and
 * returns 1 (0 otherwise) if it is a Valid sample for Montecarlo.
 * If the sample is Valid, the other return parameters are also meaningfull and
 * are the Body and (RuleID, Subst, HeadIS) of the rule that matches with the 
 * given Goal and Explan.
 *
 * This comment has to be fixed.
 */
findrule(Goal, Explan, Valid, Body, (HeadId, RuleId, Subst)) :-
	rule(Goal, _Prob, Required, RuleId, Subst, _Heads, HeadsList, Body),
	sample(HeadsList, HeadId),
	not_already_present_with_a_different_head(HeadId, RuleId, Subst, Explan),
	(HeadId =:= Required ->
		Valid = 1;
		Valid = 0).
	


/* sample(Heads, RuleId, HeadId, Subst)
 * ------------------------------------
 * This tail recursive predicate samples a random head among the given Heads of
 * the given RuleId and returns its HeadId and Subst.
 */
sample(HeadList, HeadId) :-
	random(Prob), 
	sample(HeadList, 0, 0, Prob, HeadId), !.

sample([_HeadTerm:HeadProb|Tail], Index, Prev, Prob, HeadId) :-
	Succ is Index + 1,
	Next is Prev + HeadProb,
	(Prob =< Next ->
		HeadId = Index;
		sample(Tail, Succ, Next, Prob, HeadId)).



/* cycle([(IsSample, Body, [Step])|Tail], Explan0, Explan1, Found)
 * -----------------------------------------------------------------
 * This tail recursive predicate analyzes the given Body and Step to reach it
 * and returns 0 as it's not a Valid sample for Montecarlo. 
 * If it is Valid, it looks for a solution to the Body and the given Goals 
 * starting from the Step and the given Explan and returns 1 if it finds a
 * Valid one.
 * If it does not find it, it considers the next Body and Step and returns their
 * Valid value.
 *
 * NB: This comment needs to be updated.
 */
cycle([], Explan, Explan, 0).

cycle([(0, _Goals, Step)|Tail], Explan0, Explan1, IsSample) :- !,
	append(Step, Explan0, Explan2),
	cycle(Tail, Explan2, Explan1, IsSample).

cycle([(1, Goals, Step)|Tail], Explan0, Explan1, IsSample) :- 
	append(Step, Explan0, Explan),
	main(Goals, Explan, Explan2, Valid),
	(Valid == 1 ->
		Explan1 = Explan2,
		IsSample = 1;
		cycle(Tail, Explan2, Explan1, IsSample)).