/*
\subsubsection{normal/3}
*/
:- arithmetic_function(normal/2).

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).
/*
\subsubsection{beta/1}
Simple beta distributions ($0 \le x \le 1, mean=B$).
Only B = 0.1, 0.2, 0.3, ... 0.9, 1 is supported.
*/
:- arithmetic_function(beta/1).

beta(B,X) :- beta1(B,X),!.
beta(B,X) :- B1 is 1 - B, beta1(B1,Y),X is 1 - Y.

beta1(0.50,X) :- X is ranf.
beta1(0.60,X) :- X is ranf^0.67.
beta1(0.67,X) :- X is ranf^0.5.
beta1(0.75,X) :- X is ranf^0.33.
beta1(0.80,X) :- X is ranf^0.25.
beta1(0.9,X)  :- X is ranf^(1/9).
beta1(1,1).
/*
\subsubsection{gamma/2}

Gamma ranges $0 \le x \le \infty$.
with $mean = Mean$, $skew= Alpha$.
At low skew, (e.g. x=2), x always near mean
At $skew \ge 20$, gamma becomes normal

This only supports integer values for $Mean$ and $Alpha$
*/
:- arithmetic_function(gamma/2).
gamma(Mean,Alpha,Out) :-
	Beta is Mean/Alpha,
	(Alpha > 20
	-> 	Mean is Alpha * Beta,
		Sd is sqrt(Alpha*Beta*Beta),
		Out is normal(Mean,Sd,Out)
	;	gamma(Alpha,Beta,0,Out)).

gamma(0,_,X,X) :- !.
gamma(Alpha,Beta, In, Gamma) :-
	Temp is In + ( -1 * Beta * log(1-ranf)),
	Alpha1 is Alpha - 1,
	gamma(Alpha1,Beta,Temp,Gamma).

/*
\subsection{one/3} 
Pick one thing from a list.
*/
one(Order,L,X)  :- 
	one(Order,L,X,_).

one(l2r,L0,X,L) :- any( L0,X,L).
one(ror,L0,X,L) :- rone(L0,X,L).
/* 
Pick in left-to-right order.
*/
any([H|T],H,T).
any([H|T],X,[H|T1]) :- any(T,X,T1).
/*
Pick in a random order.
*/
rone(L0,X,L) :- length(L0,N), rone(L0,N,X,L,_).

rone([H],_,H,[],0) :- !.
rone([H|T],N0,X,Out,N) :-
	N   is N0 - 1,
	Pos is random(N0) + 1,
	oneLessN(Pos,[H|T],One,Rest),
	(X=One, Out=Rest
	; Out=[One|Tail], N1 is N0 - 1, rone(Rest,N1,X,Tail,_)).

oneLessN(1,[H|T],H,T) :- !.
oneLessN(N0,[H|T],X,[H|Rest]) :- N is N0 - 1, oneLessN(N,T,X,Rest).	
/* 
In the worst case, $oneLessN$,  is linear on length of list.
I spent some time trying to improve on this using
btrees that were built via a pre-sorting of the values.
problem was that, on deletion, I could not ensure that
the balance will remain, unless I got kinda complicated.
so I did some  timing tests on the rone using $oneLessN$
and it pulls items at random from lists of
length 80 in $10^{-4}$ seconds so, what the heck.
*/
