%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Agostino Dovier, August 29, 2012
%% TESTED with BProlog
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% MAIN: Look for a tree of T nodes (including leaves and root)
%%% that justifies ALL but (AT MOST) K characters.

phylo_decision(Characters, K, T, Tree):-
        statistics(runtime,_),
        size(Characters,M,N),
        tree_domain(T,Tree),
        ancestor(N,M,AncestorCharacters),
        constrain(T, Characters, AncestorCharacters, Tree, Compatible),
        addcompatibility(Compatible,M,K),
        format("Start searching ~n",[]),
        term_variables(Compatible,VC),
        term_variables(AncestorCharacters,VAC),
        append(Tree,VC,TEMP),append(TEMP,VAC,VARS),
        labeling([], VARS), !,
        statistics(runtime,[_,Time]),
        %%% DEBUG
        format("Sol in ~d ms ~n",[Time]),
        print_sol(Tree,Characters, AncestorCharacters,Compatible).
        
phylo_decision(_,_,_,_) :-         
        statistics(runtime,[_,Time]),
		format("Unsat, stated in ~d ms ~n",[Time]).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Variables and Domains
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Tree = [L1,...,LN | Internal Nodes ]  
%%% Tree has T elements. T1 is T+1    
%%% Assume that T > N 
%%% Tree = [X1,...,XN,...,XT] Xi=j means that Xi is "son" of Xj   
%%% Xi=T1 means isolated (useless node) or root                
                                                                 
tree_domain(T, Tree) :- 
        T1 is T+1, 
        length(Tree,T), 
        domain(Tree,1,T1).
        
%%% ancestor generates the matrix of M x N
%%% Boolean characters
%%% M is the number of (Boolean) characters
%%% N is the number of taxa (leaves)

ancestor(0,_,[]).
ancestor(I,M,[C|Characters]):-
        I > 0,  I1 is I-1,
        length(C,M),
        domain(C,0,1), %%% Boolean characters
        ancestor(I1,M,Characters).
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
constrain(T, InputCharacters, AncestorCharacters, Tree, Compatible) :-
       size(InputCharacters,M,N),
       N1 is N+1, T1 is T+1,           
       %%%%%%%%%%%%%%%%%%%%%%%% REDUCE SYMMETRIES:
       Tree = [N1|_],       %%% Set the first value 
       append(_,[T1],Tree), %%% Set the last element as the root        
       monotone(Tree),      %%% Set parents in weak increasing order
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       binary(N1,Tree,T1),  %%% The tree is binary (N1 is the first internal node)
       count(T1,Tree,#=,1), %%% Unique root (no forests)
       append(InputCharacters, AncestorCharacters, Characters),
       %%% OPTIONAL CONSTRAINT:
       %%% alldifferenttuples(Characters), %%% All taxa and internal nodes have at least one difference in one character
       %%% Descendant and compatible constraints
       son_constraints(1,N1, N1, T1,Characters,Tree),
       build_compatible_vector(1,2,N1,N1,T1,Characters,Tree,Compatible,M).           

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CONSTRAINTS (in order)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% monotone states that parents are chosen in weak increasing order
%%% Can be replaced by a global constraint
%%% It removes some symmetrical equivalent solutions

monotone([]).     
monotone([_]).
monotone([A,B|R]) :-
      A #=< B,
      monotone([B|R]).     
    
%%%% binary states that Tree is a binary constraint
%%%% (the number of nodes that have a certain parent "I" is = 2)

binary(T1,_,T1) :-!.
binary(I,Tree,T1) :-
      count(I,Tree,#=,2),
      I1 is I+1,
      binary(I1,Tree,T1). 

%%%    alldifferent (on tuples) 
%%%    Implemented viewing a tuple as the binary representation of
%%%    a number and using the built-in global constraint all_different

alldifferenttuples(Characters) :-
       build_list(Characters,Weighted),
       all_different(Weighted).

%%% AUX: Converting tuples to integers:

build_list([], []).
build_list([C|Characters], [W|Weighted]) :-
     build_list_aux(C, W),
     build_list(Characters, Weighted).
build_list_aux([], 0).
build_list_aux([C|Cs], W) :-
     W #= C + 2*Ws,
     build_list_aux(Cs, Ws).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% son_constraints states that the sons differ in exactly one position
%%% from the father. Moreover, that 0->1 is allowed, while 1->0 is not.
%%% (this is a cladistic constraint... can be eliminated for other)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Loop control
son_constraints(T1,_,_, T1, _, _ ) :- !.
son_constraints(I,T1, N1, T1,Characters,Tree) :- !,
     I1 is I+1,
     J is max(N1,I1),
     son_constraints(I1,J,  N1, T1,Characters,Tree).     

%%% Main case:  Tree[i] =  J   
son_constraints(I,J, N1, T1,Characters,Tree) :- 
     nth1(I,Tree,TreeI),
     nth1(I,Characters,CharI),
     nth1(J,Characters,CharJ),
     %%%
     Flag in 0..1,
     (TreeI #= J) #<=> Flag,
     %%%% Optional constraint for Cladistic characters
     %%%% can be removed if we don't know the ``order'' of evolution    
     increasing(CharI,CharJ,Flag),
     %%%
     difference(CharI,CharJ,DeltaIJ),
     %%% HERE it can be #= 1 in case of ``Hamming'' Trees
     Flag #=> DeltaIJ #=< 1, 
     %%% recursive call
     J1 is max(N1,J+1),
     son_constraints(I,J1, N1, T1,Characters,Tree).

%%%%% AUX

difference([],[],0).
difference([A|R],[B|S],D) :-
    %%% D #= abs(A - B) + C, 
    D #= A - B + C, %% no need of abs(..), in case of cladistic
    difference(R,S,C).

increasing([],[],_).
increasing([A|R],[B|S],F) :-
    F #=> A #>= B, 
    increasing(R,S,F).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% cross_constraint generates a "Compatible" vector.
%%% If Compatible[i]=0 then the tree is compatible 
%%% with character i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% loop control
build_compatible_vector(T1,_,_, _,T1, _,_,Compatible,M) :- !,
     length(Compatible,M),
     domain(Compatible,0,0).
build_compatible_vector(I,T1,_, N1, T1,Characters,Tree, Compatible,M) :- !,
     I1 is I+1,
     J is I1+1,
     K is max(N1,J+1),
     build_compatible_vector(I1,J,K,  N1, T1,Characters, Tree, Compatible,M).
build_compatible_vector(I,J,T1, N1, T1,Characters,Tree,Compatible,M) :- !,
     J1 is J+1,
     K is max(J1+1,N1),
     build_compatible_vector(I,J1,K, N1, T1,Characters,Tree, Compatible,M).
     
%%% Main case:  Tree[i] = Tree[j] = K   
%%% It identifies a change and filter out wrong cases

build_compatible_vector(I,J,K, N1, T1,Characters,Tree,CompatibleNew,M) :- 
     nth1(I,Tree,TreeI),      nth1(I,Characters,CharI),
     nth1(J,Tree,TreeJ),      nth1(J,Characters,CharJ),
     nth1(K,Characters,CharK),
     %%%
     compatibility_analysis(CharI,CharJ,CharK,Mask),
     Flag in 0..1,
     (TreeI #= K #/\ TreeJ #= K) #<=> Flag,
     compatibility_sum(CompatibleOld,Mask,Flag,CompatibleNew), 
     %%% recursive call
     K1 is K+1,
     build_compatible_vector(I,J,K1, N1, T1,Characters,Tree,CompatibleOld,M).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% compatibility analysis, using a table, 
%%% computes the local compatibility in a node
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
compatibility_analysis([],[],[],[]).
compatibility_analysis([CI|CharI],[CJ|CharJ],[CK|CharK],[D|Mask]) :-
     sontable(TABLE),
     [ (CI,CJ,CK,D) ] in TABLE,
     compatibility_analysis(CharI,CharJ,CharK,Mask).

%%%%% Consider (CI,CJ,CK). 8 cases are possible
% i     0 0  0 -> ok  -> no change: 0
% ii    0 0  1 -> NOT ALLOWED (decreasing)
% iii   0 1  0 -> ok -> change: 1
% iv    0 1  1 -> NOT ALLOWED (decreasing)
% v     1 0  0 -> ok -> change: 1
% vi    1 0  1 -> NOT ALLOWED (decreasing)
% vii   1 1  0 -> NOT ALLOWED (double change) 
% viii  1 1  1 -> ok -> no change
%%%%%%%%%%%%%%%%%%%   
%%% These case are stored in the following table
%%% No change:0. Allowed change:1. NOT ALLOWED: 2

sontable([
( 0,0,0,0 ),
( 0,0,1,2 ),
( 0,1,0,1 ),
( 0,1,1,2 ),
( 1,0,0,1 ),
( 1,0,1,2 ),
( 1,1,0,2 ),
( 1,1,1,0 )
]).

%%%% Here the sum of the changes per characters is made

compatibility_sum([],_,_,[]). 
compatibility_sum([CO|CompatibleOld],[D|Mask],FlagIJK,[CN|CompatibleNew]) :-
     CN #= CO + D * FlagIJK,
     compatibility_sum(CompatibleOld,Mask,FlagIJK,CompatibleNew).     
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% If Compatible[i] <= 1 then that character is compatible
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

addcompatibility(Compatible,M,K) :-
    count(0,Compatible,#=,Zeros),
    count(1,Compatible,#=,Ones),
    M - Zeros - Ones  #=<  K.

%%% AUX
% Compute input size

size(Characters,M,N) :-
     length(Characters,N),
     Characters=[C|_],
     length(C,M).

%%%%% PRINTING AND EXAMPLES
				
print_sol(Tree,C1,C2,Compatible) :-
    discover_compatible(1,Compatible), 
    append(C1,C2,Chars),
    length(Tree,T),
    print_sol_aux(1,T,Tree,Chars).

discover_compatible(_,[]) :- nl.
discover_compatible(I,[CI|R]) :- 
     (CI==0 -> format("Caracter ~d never changed~n",[I]);
      CI==1 -> format("Caracter ~d changed once~n",[I]);
      true  -> nop),
      I1 is I+1,
      discover_compatible(I1,R).

nop.     
        
print_sol_aux(A,B,_,_) :- 
    A>B, !,
    nl.
    
print_sol_aux(I,T,Tree,Chars) :-
    nth1(I,Tree,F),
    nth1(I,Chars,CI),
    write(I),   write(' :'),
    write(CI),  write(' -- '),
    (F > T  ->  write(root);
     F=< T  -> (write(F),  write(' :'),
                nth1(F,Chars,CF),
                write(CF))), 
    nl,
    I1 is I+1,
    print_sol_aux(I1,T,Tree,Chars).
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SOME INSTANCES. call as (eg)
%%% :- phylo(L,0),phylo_decision(L,1,9,Tree).    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Hand/Father 1,9 ok 0,9 no

phylo([ [1,0], [1,1], [0,1], [0,0], [0,0] ], 0).

%%% Camin/Sokal example 3,7 ok

phylo([[1,0,1,1], [0,1,1,1], [0,1,0,1], [1,1,0,0], [1,0,1,0]],1).

%%% Good example 0,5

phylo([[0,1,0,1], [0,0,1,1], [1,0,0,0]],2).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
