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/mcintyre.pl
Fabrizio Riguzzi be449b3aef updated cplint
2011-10-22 16:33:04 +02:00

393 lines
9.9 KiB
Prolog

/*
LPAD and CP-Logic reasoning suite
File mcintyre.pl
Monte Carlo inference for LPADs
Copyright (c) 2011, Fabrizio Riguzzi
*/
:-dynamic rule_n/1,setting/2.
:-use_module(library(random)).
:-use_module(library(lists)).
:- set_prolog_flag(unknown,fail).
:-source.
rule_n(1).
setting(epsilon_parsing, 0.0001).
setting(compiling,true).
/* values: true, failse */
/* k
*
* numnber of samples
*
* Default value: 1000
*/
setting(k, 1000).
/* min_error
* ---------
* This parameter shows the threshold for the probability interval.
*
* Default value: 0.01
*/
setting(min_error, 0.01).
s(Goals, Samples, CPUTime, WallTime, Lower, Prob, Upper):-
solve(Goals, Samples, CPUTime, WallTime, Lower, Prob, Upper).
solve(Goals, Samples, CPUTime, WallTime, Lower, Prob, Upper) :-
% Retrieving functional parameters...
setting(k, K),
setting(min_error, MinError),
% Resetting the clocks...
statistics(cputime,[_,_]),
statistics(walltime,[_,_]),
% Performing resolution...
montecarlo_cycle(0, 0, Goals, K, MinError, Samples, Lower, Prob, Upper),
% Taking elapsed times...
statistics(cputime,[_,CT1]),
CPUTime is CT1/1000,
statistics(walltime,[_,WT1]),
WallTime is WT1/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_cycle(N0, S0, Goals, K, MinError, Samples, Lower, Prob, Upper):-!,
montecarlo(K,N0, S0, Goals, N, S),
P is S / N,
D is N - S,
Semi is 1.95996 * 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
*/
%format("Batch: samples ~d positive ~d interval ~f~n",[N,S,Int]),
%flush_output,
%((S > 5, D > 5,
((Int < MinError; Int =:= 0) ->
Samples is N,
Lower is P - Semi,
Prob is P,
Upper is P + Semi
;
montecarlo_cycle(N, S, Goals, K, MinError, Samples, Lower, Prob, Upper)
).
montecarlo(0,N,S , _Goals,N,S):-!.
montecarlo(K1,Count, Success, Goals,N1,S1):-
abolish_all_tables,
eraseall(exp),
(Goals->
Valid=1
;
Valid=0
),
N is Count + 1,
S is Success + Valid,
K2 is K1-1,
montecarlo(K2,N, S, Goals, N1,S1).
member_eq(Item, [Head|_Tail]) :-
Item==Head, !.
member_eq(Item, [_Head|Tail]) :-
member_eq(Item, Tail).
list2and([X], X) :-
X\= (_, _) , !.
list2and([H|T], (H, Ta)) :- !,
list2and(T, Ta).
list2or([X], X) :-
X \= ( _ ; _ ) , !.
list2or([H|T], (H ; Ta)) :- !,
list2or(T, Ta).
get_var_n(R,S,Probs,V):-
(v(R,S,V)->
true
;
length(Probs,L),
add_var(L,Probs,V),
assert(v(R,S,V))
).
sample_head(_HeadList,R,VC,NH):-
recorded(exp,(R,VC,NH),_),!.
sample_head(HeadList,R,VC,NH):-
sample(HeadList,NH),
recorda(exp,(R,VC,NH),_).
generate_rules_fact([],_HeadList,_VC,_R,_Probs,_N,[]).
generate_rules_fact([(Head:_P1),(null:_P2)],_HeadList,VC,R,Probs,N,[Clause]):-!,
(Clause=(Head:-(sample_head(Probs,R,VC,NH),NH=N))).
generate_rules_fact([(Head:_P)|T],HeadList,VC,R,Probs,N,[Clause|Clauses]):-
(Clause=(Head:-(sample_head(Probs,R,VC,NH),NH=N))),
N1 is N+1,
generate_rules_fact(T,HeadList,VC,R,Probs,N1,Clauses).
generate_clause(Head,Body,_HeadList,VC,R,Probs,_BDDAnd,N,_Builtin,Clause):-
(ground(VC)->
Clause=(Head:-(sample_head(Probs,R,VC,NH),NH=N,Body))
;
Clause=(Head:-(Body,sample_head(Probs,R,VC,NH),NH=N))
).
generate_rules([],_Body,_HeadList,_VC,_R,_Probs,_BDDAnd,_N,_Builtin,[]).
generate_rules([(Head:_P1),(null:_P2)],Body,HeadList,VC,R,Probs,BDDAnd,N,Builtin,[Clause]):-!,
generate_clause(Head,Body,HeadList,VC,R,Probs,BDDAnd,N,Builtin,Clause).
generate_rules([(Head:_P)|T],Body,HeadList,VC,R,Probs,BDDAnd,N,Builtin,[Clause|Clauses]):-
generate_clause(Head,Body,HeadList,VC,R,Probs,BDDAnd,N,Builtin,Clause),
N1 is N+1,
generate_rules(T,Body,HeadList,VC,R,Probs,BDDAnd,N1,Builtin,Clauses).
/* set(Par,Value) can be used to set the value of a parameter */
set(Parameter,Value):-
retract(setting(Parameter,_)),
assert(setting(Parameter,Value)).
extract_vars(Variable, Var0, Var1) :-
var(Variable), !,
(member_eq(Variable, Var0) ->
Var1 = Var0;
append(Var0, [Variable], Var1)).
extract_vars(Term, Var0, Var1) :-
Term=..[_F|Args],
extract_vars_list(Args, Var0, Var1).
extract_vars_list([], Var, Var).
extract_vars_list([Term|Tail], Var0, Var1) :-
extract_vars(Term, Var0, Var),
extract_vars_list(Tail, Var, Var1).
difference([],_,[]).
difference([H|T],L2,L3):-
member_eq(H,L2),!,
difference(T,L2,L3).
difference([H|T],L2,[H|L3]):-
difference(T,L2,L3).
process_head(HeadList, GroundHeadList) :-
ground_prob(HeadList), !,
process_head_ground(HeadList, 0, GroundHeadList).
process_head(HeadList, HeadList).
/* process_head_ground([Head:ProbHead], Prob, [Head:ProbHead|Null])
* ----------------------------------------------------------------
*/
process_head_ground([Head:ProbHead], Prob, [Head:ProbHead|Null]) :-!,
ProbLast is 1 - Prob - ProbHead,
setting(epsilon_parsing, Eps),
EpsNeg is - Eps,
ProbLast > EpsNeg,
(ProbLast > Eps ->
Null = [null:ProbLast];
Null = []).
process_head_ground([Head:ProbHead|Tail], Prob, [Head:ProbHead|Next]) :-
ProbNext is Prob + ProbHead,
process_head_ground(Tail, ProbNext, Next).
ground_prob([]).
ground_prob([_Head:ProbHead|Tail]) :-
ground(ProbHead), % Succeeds if there are no free variables in the term ProbHead.
ground_prob(Tail).
get_probs([], []).
get_probs([_H:P|T], [P1|T1]) :-
P1 is P,
get_probs(T, T1).
builtin(_A is _B).
builtin(_A > _B).
builtin(_A < _B).
builtin(_A >= _B).
builtin(_A =< _B).
builtin(_A =:= _B).
builtin(_A =\= _B).
builtin(true).
builtin(false).
builtin(_A = _B).
builtin(_A==_B).
builtin(_A\=_B).
builtin(_A\==_B).
builtin(length(_L, _N)).
builtin(member(_El, _L)).
builtin(average(_L, _Av)).
builtin(max_list(_L, _Max)).
builtin(min_list(_L, _Max)).
builtin(nth0(_, _, _)).
builtin(nth(_, _, _)).
builtin(eraseall(_Id)).
builtin(recordzifnot(_Id, _Item, _)).
parse(FileIn,FileOut):-
open(FileIn,read,SI),
read_clauses(SI,C),
close(SI),
process_clauses(C,[],C1),
open(FileOut,write,SO),
write_clauses(C1,SO),
close(SO).
process_clauses([end_of_file],C,C).
process_clauses([H|T],C0,C1):-
(expand_term(H,H1)->
true
;
H1=H
),
(H1=[_|_]->
append(C0,H1,C2)
;
append(C0,[H1],C2)
),
process_clauses(T,C2,C1).
read_clauses(S,[Cl|Out]):-
read_term(S,Cl,[]),
(Cl=end_of_file->
Out=[]
;
read_clauses(S,Out)
).
write_clauses([],_).
write_clauses([H|T],S):-
write(S,H),
write(S,'.'),
nl(S),
write_clauses(T,S).
sample(HeadList, HeadId) :-
random(Prob),
sample(HeadList, 0, 0, Prob, HeadId), !.
sample([HeadProb|Tail], Index, Prev, Prob, HeadId) :-
Succ is Index + 1,
Next is Prev + HeadProb,
(Prob =< Next ->
HeadId = Index;
sample(Tail, Succ, Next, Prob, HeadId)).
get_next_rule_number(R):-
retract(rule_n(R)),
R1 is R+1,
assert(rule_n(R1)).
user:term_expansion((Head :- Body), Clauses):-
% disjunctive clause with more than one head atom
setting(compiling,true),
Head = (_;_), !,
list2or(HeadListOr, Head),
process_head(HeadListOr, HeadList),
get_next_rule_number(R),
get_probs(HeadList,Probs),
extract_vars((Head:-Body),[],VC),
generate_rules(HeadList,Body,HeadList,VC,R,Probs,_BDDAnd,0,_Builtin,Clauses).
user:term_expansion((Head :- Body), Clauses) :-
% disjunctive clause with a single head atom
setting(compiling,true),
((Head:-Body) \= ((user:term_expansion(_,_) ):- _ )),
Head = (H:_), !,
list2or(HeadListOr, Head),
process_head(HeadListOr, HeadList),
get_next_rule_number(R),
get_probs(HeadList,Probs),
extract_vars((Head:-Body),[],VC),
generate_clause(H,Body,HeadList,VC,R,Probs,_BDDAnd,0,_Builtin,Clauses).
user:term_expansion(Head,Clauses) :-
% disjunctive fact with more than one head atom
setting(compiling,true),
Head = (_;_),!,
list2or(HeadListOr, Head),
process_head(HeadListOr, HeadList),
extract_vars((Head),[],VC),
get_next_rule_number(R),
get_probs(HeadList,Probs),
generate_rules_fact(HeadList,HeadList,VC,R,Probs,0,Clauses).
user:term_expansion(Head,Clause) :-
% disjunctive fact with a single head atom
setting(compiling,true),
(Head \= ((term_expansion(_,_)) :- _ )),
Head = (Head1:_), !,
list2or(HeadListOr, Head),
process_head(HeadListOr, HeadList),
get_probs(HeadList,Probs),
extract_vars((Head),[],VC),
get_next_rule_number(R),
(Clause=(Head1:-(sample_head(Probs,R,VC,NH),NH=0))).