%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% AGOSTINO DOVIER
%%% TESTED IN B PROLOG
%%% AUGUST 31 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Simplifying assumptions:
%%% Assume Gs all of the same length
%%% Assume all of them have at least a "2"
%%% Assume L > 1

haplo_decision(Gs, L):-
        statistics(runtime,_),
		size(Gs,Ng),
        generate_haplo(Hs,L,Ng),
        constraint(Gs,Hs),
        %%%    Hs = [[0,0,0,1,1,1,0]|_],
        term_variables(Hs, Vars), 		
        format("Constraints posted. Start searching~n",[]),
        labeling([ffc], Vars), !,
        statistics(runtime,[_,T]),
		format("Sol in ~d ms ~n",[T]),
        print_sol(Hs).
        
haplo_decision(_,_) :-         
        statistics(runtime,[_,T]),
		format("Unsat, stated in ~d ms ~n",[T]).

%
size([G|_],Ng) :-
     length(G,Ng).

%
generate_haplo(Hs,L,N) :-
     length(Hs,L),
     generate_haplo_aux(Hs,N).      

generate_haplo_aux([],_).        
generate_haplo_aux([H | Hs],N) :-
     length(H,N),domain(H,0,1),
     generate_haplo_aux(Hs,N).

constraint(Gs,Hs) :-
	 %lex_chain(Hs,[op(#<)]),
     strict_order(Hs),
     constraint_aux(Gs, Hs).  

%%% Not built-in in B Prolog

strict_order([]).
strict_order([_]).
strict_order([A,B|R]) :-
     strict_order_aux(A,B,1),
     strict_order([B|R]).

strict_order_aux([],[],0).
strict_order_aux([A|B],[C|D],Fo) :-
     strict_order_aux(B,D,Fi),
     ordtable(T),
     [(A,C,Fi,Fo)] in T.
     %A #= C #=> Fo #= Fi,
     %A #< C #=> Fo #= 1,
     %A #> C #=> Fo #= 0.
 
ordtable( [
% A=C
(0,0, 0,0),     
(0,0, 1,1),     
(1,1, 0,0),     
(1,1, 1,1),     
% A < C
(0,1, 0,1),     
(0,1, 1,1),
% A > C    
(1,0, 0,0),     
(1,0, 1,0)
]).
     
        
%% Apply constraints to all genotypes %%
		
constraint_aux([], _ ).
constraint_aux([G|Gs], Hs ):-
		single_gen(G, Hs, FlagsforG),
%% 1. OR of the flags
        %1 #= max(FlagsforG), %%% OR
%% 2. Using the sum        
		sum(FlagsforG,#=,1), 
%% 3. Being more relaxed
		%sum(FlagsforG,#>,0), 
%% Solo in questo caso        
        constraint_aux(Gs, Hs).
		
%% Find two haplotypes for G

single_gen(_, [], []).
single_gen(G, [H|R], F):-
		single_gen_aux(G, H, R, F1),
		single_gen(G, R, F2),
  		append(F1,F2,F).
    
single_gen_aux(_, _, [], []).
single_gen_aux(G, H1, [H2|R], [F|Flags]) :-
		solve(G, H1, H2, FpairH1H2), 
        F #= min(FpairH1H2), %%% AND BIT x BIT
 		single_gen_aux(G, H1, R, Flags).
		
%% solve a single genotype
		
solve([], [], [], []).
solve([G|Gs], [H|Hs], [K|Ks], [F|Flags]):- 
		solve(Gs, Hs, Ks, Flags),
%%%% 1 DIRECT
        (G=2 -> ( F #<=>  H #\= K);
		 G<2 -> ( F #<=> (H#=G #/\ K#=G))).
%%%% 2 Model sum
%        F #<=> (G #= H+K).


%%%%% PRINTING AND EXAMPLES
				
print_sol(Hs) :- 
        %foreach(H,Hs) do (print(H),nl).
        forall(member(H,Hs), (print(H),nl)).


%% length 5, 6

input(1,Gs) :-
Gs = [ [2,1,1,0,2,0,2],
       [2,0,0,2,1,1,0],
       [2,0,1,2,2,0,2],
       [2,2,2,2,1,2,2]].
   
input(2,Gs) :-
Gs = [ [2,2,2,2],
       [2,2,0,0],
       [0,2,2,0],
       [0,0,2,2] ].

input(3,Gs) :-
Gs = [ [2,2,2,2,0],
       [2,2,0,0,0],
       [0,2,2,0,1],
       [0,1,2,0,1],
       [0,0,2,2,1] ].

%%% 10 no, 11 yes

input(4,Gs) :-
Gs = [ [2,1,1,0,2,0,2],
       [2,0,0,2,1,1,0],
       [2,0,1,2,2,0,2],
       [2,2,2,2,1,2,2],
       [1,2,0,1,2,0,1],
       [1,2,0,2,2,0,1],
       [1,2,2,2,1,2,2],
       [1,1,2,2,1,2,2],
       [1,1,2,1,1,2,2],
       [1,0,2,2,1,0,1] ].
