% assertions
:- discontiguous want/1, got/1,set/4.

% main drivers
set1(X) :- set(X,Y,[],L),print(X = Y),nl,println(L),print(ind),nl,independents(L,I),println(I).
eg(1)   :- set1(goal).

% querying semantics
wants     :- wants1(L), println(L).
gots      :- gots1(L), print(L).

wants1(L) :- setof(X,wants2(X),L).
gots1(L)  :- setof(X, gots2(X),L).

wants2(X/L) :- want(X),myBagof(S,X^status(X,S),L).
gots2( X/L) :- got(X), myBagof(S,X^status(X,S),L).

status(X,ungot) :-  want(X), \+ got(X).
status(X,want ) :- got(X), \+ want(X).

roots(L) :- setof(X,root(X),L).
root(X)  :- clause(set(X,_,_,_),(trig(_,_,_,_),_)).

independents(W,L) :- setof(X,W^independent(X,W),L).
independent(X=V,W) :- root(X),member(X=V0,W), V is 5+ V0*10.
 
% load-time macros
term_expansion(X = Y,Out ) :- xpand(X=Y,Out).

xpand(Head=Body0,[got(Head),Rule|Rest]) :- 
	xpand1(Body0,Arg,N0,Body,Rest,[]), N is N0,
	expand_term((set(Head,Out) --> Body,{Out is Arg/N}),Rule).

xpand1(A0 + B0, A1 + B1, A2+B2,(A,B)) --> !, xpand1(A0,A1,A2,A), xpand1(B0,B1,B2,B).
xpand1(A0 - B0, A1 - B1, A2+B2,(A,B)) --> !, xpand1(A0,A1,A2,A), xpand1(B0,B1,B2,B).
xpand1(A * N,   P, N   ,       Goal ) --> !, [want(A)], {xpand2(A=X,X*N,P,Goal)}.
xpand1(A / N,   P, 1/N ,     Goal   ) --> !, [want(A)], {xpand2(A=X,X/N,P,Goal)}.
xpand1(tr(Min,Max,Mode),X,1,{tr(Min,Max,Mode,X)}) --> [],!.
xpand1(trig(Min,Max,Mode),X,1,{trig(Min,Max,Mode,X)}) --> [],!.
xpand1(A    ,   P, 1,       Goal  ) -->    [want(A)], {xpand2(A=X,X,P,Goal)}.

xpand2(           - A = X,    P, -1*P, A=X) :- !.
xpand2(             A = X,    P,    P, A=X).

% utils

printl([H|T]) :- print(H), forall(member(X,T),format(',~p',X)).
 
println(L) :- println(L,'\t').
println(L,Tab) :- forall(member(X,L),format('~p~p\n',[Tab,X])).

myBagof(X,Y,Z) :- bagof(X,Y,Z),!.
myBagof(_,_,[]).

% DCG-based assumption space management
=(X,Y,In,     In)   :- member(X=Z,In),!,Y=Z.
=(X,Y,In,[X=Y|Out]) :- clause(set(X,_,_,_),_),!,set(X,Z0,In,Out),ones(Z0,Z),Y=Z.
=(X,Y,In,[X=Y|In])  :- Y is 1 + -2 * random(1000)/999.

% example files
tar3(F,N) :-
	tell(F),
	roots([H|T]),
	format('$~p',H),
	forall(member(X,L),format(',$~p',X)).
		printl(L),print(',>$'),print(goal),nl,
	forall(between(1,N,_),tar3a),
	told.

tar3a :-	set(goal,X,[],W0), independents(W0,W1),maplist(arg(2),W1,W), printl(W),print(','),print(X),nl.

ones(N0,N)              :- within(-1,1,N0,N).

% numerics
:- arithmetic_function(ranf/0).
ranf(X) :-  X is random(65535)/65536.

within(Bottom,Top,N,N     ) :- N >= Bottom, N=<Top.
within(Bottom,_,  N,Bottom) :- N < Bottom.
within(_,     Top,N,Top   ) :- N > Top.
	
normal(M,S,N) :- box_muller(M,S,N).
box_muller(M,S,N) :- w(W0,X), W is sqrt((-2.0 * log(W0))/W0), Y1 is X * W, N is M + Y1*S.

w(W,X) :-
	X1 is 2.0 * ranf - 1,
	X2 is 2.0 * ranf - 1,
	W0 is X1*X1 + X2*X2,
	(W0  >= 1.0 -> w(W,X) ; W0=W, X = X1).

trCDF(A,B,C,X,Y) :- X =< C -> Y is (X - A)^2/((B-A)*(C-A)) | Y is 1 -  (B - X)^2/((B-A)*(B-C)).

trig(A,B,C,Z) :- tr(A,B,C,Z0), Z is 2*(Z0 - A)/(B-A) - 1.
tr(A,B,C,_)   :- (C < A; C > B),!,print(bad(tr(min = A,max = B, mode = C))),nl,fail.
tr(A,A,_,A)   :- !.
tr(A,B,C,Z)   :- R is ranf, tr1(R,0.01, A,B,C,A,B,Z).

tr1(R,P,A,B,C,Min,Max,Z) :-
	Half is Min + (Max - Min)/2,
	trCDF(A,B,C,Half,Sample),
	myCompare(Op,P,R,Sample),
	tr2(Op,R,P,A,B,C,Min,Half,Max,Half,Z).

tr2(=,_,_,_,_,_,_  ,_   ,_  ,Z,Z).
tr2(>,R,P,A,B,C,_  ,Half,Max,_,Z) :- tr1(R,P,A,B,C,Half,Max,Z).
tr2(<,R,P,A,B,C,Min,Half,_  ,_,Z) :- tr1(R,P,A,B,C,Min,Half,Z).

myCompare(=, P,X,Y) :- abs(X - Y) =< P,!.
myCompare(Op,_,X,Y) :- compare(Op,X,Y).

testTR :-
	tell(dat),
	forall(between(1,1000,_), (tr(10,20,10.0001,X),print(X),nl)),
	told,
	shell("sort -n dat > dat.sorted").


goal = reqQual.
	
userExper 	= trig(4,7,6).
userInvolve = trig(4,6,5).

reqQualUserInput = userExper + userInvolve.

novel		= trig(6,8,7).

reqStab = reqQualUserInput + novel.
 
sysDocQual 	= trig(6,8,7).
reqCplx 	= trig(3,5,4).
reqProbSpc = reqStab+ sysDocQual + reqCplx.


reqQual = reqProbSpc + reqDevStaff + reqProcRigor.

reqDevTrnOvr = trig(5,7,6).
reqDevStfLev = trig(8,10,9).
reqDevResAvail = reqDevTrnOvr + reqDevStfLev. 

reqDevExpr = trig(5,7,6).
reqDevDomain = trig(7,10,9).
reqDevStaff = reqDevExpr + reqDevDomain + reqDevResAvail.

reqDevSched  = trig(4,6,5).
reqDevBudget = trig(6,10,8).

reqExtCstr = reqDevSched + reqDevBudget.

reqDevProcEff = trig(6,10,8).
reqDevQualOrg = trig(5,7,6).
reqProcAdhere = reqExtCstr + reqDevProcEff + reqDevQualOrg.

reqDevTool = trig(6,10,8).
devCMM     = trig(7, 9,8).
reqProcRigor = devCMM + reqProcAdhere + reqDevTool.