%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% METABAYES WITH LINEAR SAMPLING - Version 3.0
%
% Version 3 uses rational numbers to represent intervals. This avoids
% 	erros associated with lack of precision, by using Yap's
% 	arbirary precision calculations for integers.
%
% Consider we require a sample of size n. If this is drawn in a
% regular fashion then we need to find hypotheses whose cumulative
% probability is as follows.
% 
%         p = 1/(n+1),2/(n+1),..,n/(n+1)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

:- [rational].

sample(0,_,_,_,_,[]) :- !.
sample(N,Min,Pos,Neg,BK,[Hyp1|Hyps]) :-
	Prob is_rnl Min+1/(N+1)*(1-Min), Prob lt_rnl 1, !,
	%nl, write('SAMPLE '), write(N), write(' Prob='), write(Prob),  nl, nl,
/*	once(prove(Pos,Neg,Prob-Min,[0,1],[_,Max],BK,Hyp)),
	N1 is N-1,
(Hyp==fail->
	true;
	printprog(Hyp)
),
*/

(prove(Pos,Neg,Prob-Min,[0,1],[_,Max0],BK,Hyp)->
(Max0=Min ->
    Max=Prob,Hyp1=fail;
    Max=Max0,Hyp1=Hyp,printprog(Hyp1)
);
    
    Max=Prob,Hyp1=fail
),
N1 is N-1,
	sample(N1,Max,Pos,Neg,BK,Hyps).
sample(_,_,_,_,_,[]) :- !.
	%nl, write('Population exhausted, Prob=1.0'), nl, !.
	


prove([],_,_,Int,Int,Prog,Prog). 
%prove(__,_,Int,Int,fail,fail) :- !. 
prove([Atom|As],Neg,Prob,Int1,Int2,Prog1,Prog2) :-
        once(d_prove([Atom],Prog1)), !,       % Atom already proveable
        once(prove(As,Neg,Prob,Int1,Int2,Prog1,Prog2)).
prove([Atom|As],Neg,Prob,Int1,Int2,Prog1,Prog2) :-
		% Atom is a list of Constant/Variables
		% Program is list of RuleName/MetaSubstitution pairs
	findall([metasub(RuleName,MetaSub),Body],
                (metarule(RuleName,MetaSub,(Atom :- Body),OrderTest),
                OrderTest,
		\+break_constraint(metasub(RuleName,MetaSub),Prog1),
		consistent(Neg,[metasub(RuleName,MetaSub)|Prog1])),
                Matches), 
%write('CandidateMatches: '),write(Atom-Matches),nl,
Matches=[_|_],
        chooseabduction(Matches,Body1,Prob,Int1,Int3,Prog1,Prog3), %!,
	once(prove(Body1,Neg,Prob,Int3,Int4,Prog3,Prog4)),
	once(prove(As,Neg,Prob,Int4,Int2,Prog4,Prog2)). 
%prove([_|_],_,_,Int,Int,_,fail):- fail. 

chooseabduction([],_,_,_,_,_,_) :- !, fail.
chooseabduction(Matches,Body,Prob-Min,[Min1,Max1],[Min2,Max2],Prog1,Prog2) :- 
    Min lt_rnl Max1, %-***  
	Mterm =.. [m|Matches],
	functor(Mterm,_,N),
	%write('N='), write(N), write(',[Min1,Max1]='), write([Min1,Max1]),write('Prob='), write(Prob), nl,
%-***  
    (Max1 lt_rnl Prob ->
        MaxIndex=N;
        [Nm,Dm] is_rnl (N*(Prob-Min1)/(Max1-Min1)),
        MaxIndex is ceiling(Nm/Dm)
    ),
	
    (Min lt_rnl Min1 ->
        MinIndex=1;
        [NmMin,DmMin] is_rnl (N*(Min-Min1)/(Max1-Min1)),
        MinIndex is ceiling(NmMin/DmMin)
    ),
    interval_rev(MinIndex,MaxIndex,Indexes),
    element(M,Indexes), %write(M-Indexes),nl,
%-***  
	arg(M,Mterm,[MetaSub,Body]),
	%write([Min1,Max1]), write('N='), write(N), write('M='), write(M),write(MetaSub),
	Min2 is_rnl ((M-1)/N)*(Max1-Min1)+Min1,
	Max2 is_rnl (M/N)*(Max1-Min1)+Min1,
	%write([Min2,Max2]), nl,
	abduce(MetaSub,Prog1,Prog2). %, !. %-***

% produce the integers between MinIndex and MaxIndex, but the list is in a reversed order 

interval_rev(Index1,Index2,[]):-Index1>Index2,!,fail.
interval_rev(Index,Index,[Index]):-!.
interval_rev(MinIndex,MaxIndex,[MaxIndex|Indexes]):-
    NMaxIndex is MaxIndex-1,
    interval_rev(MinIndex,NMaxIndex,Indexes). 

abduce(MetaClause,Prog,Prog) :- 
    exist(MetaClause,Prog), !.
abduce(MetaClause,Prog1,[MetaClause|Prog1]):-
    max_kBound(KBound),
    length([MetaClause|Prog1],N),N=<KBound.




element(H,[H|_]).
element(H,[_|T]) :- element(H,T).

append([],L,L).
append([H|T],L,[H|R]) :- append(T,L,R).

% Test consistency against negative examples

consistent([],Prog) :- !.
consistent([E|Neg],Prog) :-
	not(d_prove([E],Prog)),
	consistent(Neg,Prog).
%
% Printing predicates

printprogs([]) :- !.
printprogs([Prog|Progs]) :-
	printprog(Prog), nl,nl,
	printprogs(Progs).

printprog(fail) :- !. % nl,write('FAIL'),
printprog(Ms) :-
	converts(Ms,Cs), nl, sort(Cs,Cs1),
	numbervars(Cs1,0,_), printclauses(Cs1), !.
printprog(_) :- nl, write('FAIL'), !.

converts([],[]) :- !.
converts([metasub(RuleName,MetaSub)|MIs],[Clause|Cs]) :-
	metarule(RuleName,MetaSub,Clause,_),
	converts(MIs,Cs), !.

printclauses([]) :- nl, !.
printclauses([C|Cs]) :-
	printclause(C), nl,
	printclauses(Cs).

printclause((Head :- [])) :-
	printatom(Head), write('.').
printclause((Head :- Body)) :-
	printatom(Head), write(' : - '),
	printatoms(Body).

printatom(List) :- Atom =.. List, write(Atom).

printatoms([A]) :- printatom(A),  write('.'), !.
printatoms([A|As]) :- printatom(A), write(', '), printatoms(As), !.

% Filter out failed program constructions

failfree([]) :- !.
failfree([fail]) :- !, fail.
failfree([Term|Terms]) :-
	Term =.. [_|Terms1],
	failfree(Terms1),
	failfree(Terms), !.

filterfails([],[]) :- !.
filterfails([H|T1],[H|T2]) :-
	failfree(H), !, filterfails(T1,T2).
filterfails([_|T1],L2) :-
	filterfails(T1,L2).

maxhyp([X],X).
maxhyp([H1|T],MaxH) :-
        maxhyp(T,H2),
        maxhyp(H1,H2,MaxH), !.

maxhyp(Hyp1,Hyp2,Hyp1) :-
	prior(Hyp1,P1), prior(Hyp2,P2),
	P1>P2, !.
maxhyp(_,Hyp,Hyp).

bayes_predict(X,Hs,BK,Prob) :-
        bayes_predict1(X,Hs,BK,PosProb,NegProb),
        Prob is PosProb/(PosProb+NegProb), !.

bayes_predict1(X,[],BK,0,0) :- !.
bayes_predict1(X,[Hyp|Hs],BK,Prob+PosProb,NegProb) :-
	prior(Hyp,Prob),
        append(Hyp,BK,Prog),
        d_prove([X],Prog), !,
        bayes_predict1(X,Hs,BK,PosProb,NegProb), !.
bayes_predict1(X,[Hyp|Hs],BK,PosProb,Prob+NegProb) :-
	prior(Hyp,Prob),
        bayes_predict1(X,Hs,BK,PosProb,NegProb), !.

/* moved to domain specific file
prior(Hyp,1/L^2) :-
	length(Hyp,L),
	P is 1/L^2.
*/

print_predicts(TestSet,Hs,BK) :-
        predset(TestSet,Set),
        element(X,Set),
        bayes_predict(X,Hs,BK,Prob),
        write(X), write('  BayesProb= '), write(Prob), write(', '),
        entropy(Prob,Ent), write('Entropy= '), writeln(Ent),
        fail.
print_predicts(_,_,_).

entropy(0.0,0.0) :- !.
entropy(1.0,0.0) :- !.
entropy(P,Ent) :-
        Ent is -P*log(P)/log(2)-(1-P)*log(1-P)/log(2).