% (c) Florence Benoy, Andy King and Fred Mesnard
%
% See the TPLP paper "Calculating Convex Hulls with
% a Linear Solver" for explanation and discussion
%

:- module(hull, [project/3, convex_hull/6]).
:- use_module(library(clpq)).

project(Xs, Cxs, ProjectCxs) :-
    call_residue(copy_term(Xs-Cxs, CpyXs-CpyCxs), _),
    tell_cs(CpyCxs), 
    prepare_dump(CpyXs, Xs, Zs, DumpCxs, ProjectCxs),
    dump(Zs, Vs, DumpCxs), Xs = Vs.

tell_cs([]).
tell_cs([C|Cs]) :-  {C}, tell_cs(Cs).

prepare_dump([], [], [], Cs, Cs).
prepare_dump([X|Xs], YsIn, ZsOut, CsIn, CsOut) :-
    (ground(X) ->
        YsIn  = [Y|Ys],
        ZsOut = [_|Zs],
        CsOut = [Y=X|Cs]
    ;
        YsIn  = [_|Ys],
        ZsOut = [X|Zs],
        CsOut = Cs
    ),
    prepare_dump(Xs, Ys, Zs, CsIn, Cs).

convex_hull(Xs, Cxs, Ys, Cys, Zs, Czs) :-
    scale(Cxs, Sig1, [], C1s),
    scale(Cys, Sig2, C1s, C2s),
    add_vect(Xs, Ys, Zs, C2s, C3s),
    project(Zs, [Sig1 >= 0, Sig2 >= 0, Sig1+Sig2 = 1|C3s], Czs).

scale([], _, Cs, Cs).
scale([C1|C1s], Sig, C2s, C3s) :-
    C1 =.. [RelOp, A1, B1],
    C2 =.. [RelOp, A2, B2],
    mul_exp(A1, Sig, A2),
    mul_exp(B1, Sig, B2),
    scale(C1s, Sig, [C2|C2s], C3s).

mul_exp(E1, Sigma, E2) :- once(mulexp(E1, Sigma, E2)).

mulexp(  X,   _,     X) :- var(X).
mulexp(N*X,   _,   N*X) :- ground(N), var(X).
mulexp( -X, Sig,    -Y) :- mulexp(X, Sig, Y).
mulexp(A+B, Sig,   C+D) :- mulexp(A, Sig, C), mulexp(B, Sig, D).
mulexp(A-B, Sig,   C-D) :- mulexp(A, Sig, C), mulexp(B, Sig, D).
mulexp(  N, Sig, N*Sig) :- ground(N).

add_vect([], [], [], Cs, Cs).
add_vect([U|Us], [V|Vs], [W|Ws], C1s, C2s) :-
    add_vect(Us, Vs, Ws, [W = U+V|C1s], C2s).
