This repository has been archived on 2023-08-20. You can view files and clone it, but cannot push or open issues or pull requests.
yap-6.3/packages/cplint/approx/montecarlo.pl
RIGUZZI FABRIZIO - Dipartimento di Ingegneria db2eefd0c9 added approximated cplint
2010-03-18 16:11:21 +01:00

276 lines
8.9 KiB
Prolog

/*==============================================================================
* 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)).