% pairwise(I, (L1,L2), (V1,V2)) is true iff L1[I]=V1 AND L2[I]=V2
% Prerequisites : L1,L2 are both lists of FD_var with the same size, I,V1,V2 are FD_var
%                 --- All the FD_vars MUST have a bounded domains ---
% Warning : all the FD_var are supposed to have numeric bounds (inf,sup not managed)
% Maintain bound-consistency on V1,V2 and arc-consistency on I  (To Be Verified)

% Designed for sicstus 4.1.1
% Date : June 2010
% Authors : Arnaud Gotlieb -- INRIA Rennes Bretagne Atlantique
% usage:

% ?- domain([X1,X2,X3,V1,V2],-1,0), pairwise(I, ([X1,X2,X3],[-1,0,-1]), (V1,V2)), I = 3.
%     V1 = X3, V2 = -1.
% ?- pairwise(I, ([X1,X2,X3],[0,-1,-1]), (V1,V2)), V1 = 0, X1 #< 0.
%     I in 2..3, %%  could deduce also V2 = -1 but at the cost of an awakening - suspended ctr 
% ?- pairwise(I, ([X1,X2,X3], [Y1,Y2,Y3]), (V1,V2)), X1=0, X2=0,X3=0.
%     V1 = 0.
% ?- pairwise(I, ([0,-1,-1], [Y1,Y2,Y3]), (0,-1)).
%     I = 1, Y1 = -1.
% ?- pairwise(I, ([0,-1,-1], [Y1,Y2,Y3]), (0,-1)), Y1=0.
%    no
% ?- domain([X1,X2,X3,V1,V2],0,2), pairwise(I, ([X1,X2,X3],[0,0,0]), (V1,V2)), X1 #< 2, X2 #< 2, V1=2.
%    X1 in 0..1,
%    X2 in 0..1,
%    X3 = 2,
%    V2 = 0


:- ensure_loaded(library(clpfd)).
:- ensure_loaded(library(lists)).
:- multifile clpfd:dispatch_global/4.

pairwise(I, (L1,L2), (V1,V2)) :-
        lists:is_list(L1), lists:is_list(L2),
        minmax_suspensions(L1, Susp1, 0, Len1),
        minmax_suspensions(L2, Susp2, 0, Len2),
        Len1 == Len2,
        I in 1 .. Len1,
        !,
        lists:append(Susp1, Susp2, Susp3),
        fd_global(pw(I,L1,L2,V1,V2), state, [dom(I),minmax(V1),minmax(V2)|Susp3]).
%pairwise(_, _, _) :-
%        write('Pairwise constraint problem').


clpfd:dispatch_global(pw(I,L1,L2,V1,V2), state, state, Actions) :-
        pw_solver(I,L1,L2,V1,V2,Actions).

pw_solver(I, L1, L2, V1, V2, Actions) :-
        number(I),
        !,
        lists:nth1(I,L1,X1),
        lists:nth1(I,L2,X2),
        Actions = [exit, call(V1=X1), call(V2=X2)].
pw_solver(I, L1, L2, V1, V2, Actions) :-
%        nl,nl,
        fd_set(I, SetI), fd_set(V1, SetV1), fd_set(V2, SetV2),
        fdset_to_list(SetI, ListI),
        pw_rec(ListI, SetI, L1, L2, SetV1, SetV2, NSetI, [], L1Union, [], L2Union),
        fdset_intersection(SetV1, L1Union, NSetV1),
%        write(L1Union),nl,
        fdset_intersection(SetV2, L2Union, NSetV2),
%        write(L2Union),nl,
        build_actions(NSetI, NSetV1, NSetV2, I,L1, L2, V1,V2,Actions).

pw_rec([], SetI, L1, L2, SetV1, SetV2, SetI, L1U, L1U, L2U, L2U):- !.
pw_rec([Val|S], SetI, L1, L2, SetV1, SetV2, NSetI, L1U, NL1U, L2U, NL2U) :-
        lists:nth1(Val, L1, X1),
        lists:nth1(Val, L2, X2),
        fd_set(X1, SetX1),
        fd_set(X2, SetX2),
        test_disjoint(SetX1, SetV1, Val, SetI, SetI1, L1U, L1Ua),
        test_disjoint(SetX2, SetV2, Val, SetI1, SetI2, L2U, L2Ua),
        pw_rec(S, SetI2, L1, L2, SetV1, SetV2, NSetI, L1Ua, NL1U, L2Ua, NL2U).

test_disjoint(SetX, SetV, Val, SetI, NSetI, LU, LU) :-
        fdset_disjoint(SetX, SetV),
        !,
        fdset_del_element(SetI, Val, NSetI).
test_disjoint(SetX, SetV, Val, SetI, SetI, LU, NLU) :-
        fdset_union(LU, SetX, NLU).

build_actions(NSetI, NSetV1, NSetV2, I,L1,L2,V1,V2,Actions):-
        ( empty_fdset(NSetI) ; empty_fdset(NSetV1) ; empty_fdset(NSetV2)), 
        !,
        Actions = [fail].
build_actions(NSetI, NSetV1, NSetV2, I,L1, L2, V1,V2,Actions):-
        fdset_singleton(NSetI, Ival),
        !,
        lists:nth1(Ival, L1, X1),
        lists:nth1(Ival, L2, X2),
        Actions = [exit, call(I=Ival), call(V1=X1), call(V2=X2)].
%build_actions(NSetI, NSetV1, NSetV2, I,L1, L2, V1,V2,Actions):-
%        ( fdset_singleton(NSetV1, V1val) -> (L1 = [V1=V1val],X=true) ; L1 = [] ),
%        ( fdset_singleton(NSetV2, V2val) -> (L2 = [V2=V2val],X=true) ; L2 = [] ),
%        ( X == true -> lists:append(L1,L2,Actions) ; fail ),
%        !.
build_actions(NSetI, NSetV1, NSetV2, I,L1, L2, V1,V2,Actions):-
        Actions = [I in_set NSetI, V1 in_set NSetV1, V2 in_set NSetV2].
        
%=====================================UTILS============================================
minmax_suspensions([], [], Len, Len).
minmax_suspensions([X|L], [minmax(X)|SL], Len, Len1) :-
        LenI is Len+1,
        minmax_suspensions(L,SL, LenI,Len1).



%=================================VERSION AH===========================================

                                %solver : première partie : cas ou I devient un number : on impose les valeurs puis exit
                                %   deuxième partie : quand il y a une information qui arrive sur les domaines de L1 et L2

pairwise_solver(L1,L2,I,V1,V2,Actions):-
        (number(I) ->
            elem(L1,I,A),elem(L2,I,B),Actions = [exit|Ps1],
            fdset_singleton(Set1, V1),fdset_singleton(Set2, V2),fdset_singleton(Set3, I),
            Ps1 = [A in_set Set1, B in_set Set2, I in_set Set3]

        ;
            rang2(L1,0,V1,LR1),rang2(L2,0,V2,LR2),
            list_to_fdset(LR1,FS1),list_to_fdset(LR2,FS2),
            fdset_intersection(FS1,FS2,SI),
            Actions = [I in_set SI]
        ).                            

% rang2(Colonne,entier(0), Val , liste des rang empty) : rend la liste des position ou la valeur Val peut être prise.
% ?- rang2([0,0,2,A,B,0],0,0,LR1).
%LR1 = [1,2,4,5,6] ?

rang2([A|R],N,Val,[N1|R1]) :-
        A == Val, N1 is N+ 1, rang2(R,N1,Val,R1).
rang2([A|R],N,Val,[N1|R1]) :-
        \+(number(A)), N1 is N+ 1, rang2(R,N1,Val,R1).
rang2([A|R],N,Val,R1) :-
        A \= Val, N1 is N+ 1,rang2(R,N1,Val,R1).
rang2([],_,_,[]).