% Card generation

% Classic setup
setup([[red, green, blue],
       [1, 2, 3],
       [circle, triangle, square],
       [empty, dashed, filled]]).

% Alternative setup
% setup([[red, green, blue, black],
%        [1, 2, 3, 4],
%        [circle, triangle, square, pentagon]]).

card(C) :- setup(L), maplist(member, C, L).

cards(L) :- findall(C, card(C), L).

% Set property

set([[]|_]).
set(L) :-
    heads_tails(L, H, T),
    ( all_same(H) ; all_diff(H) ),
    set(T).

heads_tails([], [], []).
heads_tails([[X|Xs]|L], [X|H], [Xs|T]) :- heads_tails(L, H, T).

all_same([_]) :- !.
all_same([X,X|Xs]) :- all_same([X|Xs]).

all_diff([_]) :- !.
all_diff([X|Xs]) :- forall(member(Y, Xs), X \= Y), all_diff(Xs).

% Probability experiments

% What is the probability P of having a set in K cards?
%   (using M Monte-Carlo iterations)
probability(K, M, P) :-
    cards(C), length(C, N),
    monte_carlo(M, C, N, K, S),
    P is S / M.

monte_carlo(0, _, _, _, 0) :- !.
monte_carlo(M, Cards, NumCards, K, SetCount) :-
    M > 0, M1 is M - 1,
    monte_carlo(M1, Cards, NumCards, K, S),
    randset(K, NumCards, Indices),
    findall(X, (member(I, Indices), nth1(I, Cards, X)), L),
    ( contains_set(L) -> SetCount is S + 1 ; SetCount = S ).

contains_set(L) :-
    setup([X|_]), length(X, N),
    select_n(N, L, S),
    set(S).

select_n(0, _, []).
select_n(N, [X|Xs], [X|L]) :- N > 0, N1 is N - 1, select_n(N1, Xs, L).
select_n(N, [_|Xs], L) :- N > 0, select_n(N, Xs, L).

% Test

% Classic setup
% ?- probability(12, 10000, P).
% P = 0.9665.
%
% ?- probability(15, 10000, P).
% P = 0.9998.
%
% Alternative setup
% ?- probability(12, 10000, P).
% P = 0.5199.
%
% ?- probability(16, 10000, P).
% P = 0.936.
%
% ?- probability(18, 10000, P).
% P = 0.9882.
%
% ?- probability(20, 10000, P).
% P = 0.9987.

% So we should use 16 of the alternative cards,
% adding another 4 when there is no set.