%%% TESTED IN B Prolog, November 8, 2012


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FOR SICSTUS PROLOG
%:- use_module(library(clpfd)).
%:- use_module(library(timeout)).
%:- use_module(library(random)).
%:- use_module(library(lists)).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% PREDICATO PRINCIPALE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Tempo in millisecondi
%%% Tecnica b_in, b_up, loc

minimizza(Soluzioni,Tempo,Tecnica):-
     vincola(Vars,Costo),
     retractall(initialtime(_)),
     statistics(runtime,[T,_]),
     assert(initialtime(T)), 
     ricerca(Vars,Costo,Soluzioni,Tempo,Tecnica).

%%% Per prova con il protein folding

minimizza(Primary,Soluzioni,Tempo,Tecnica):-
     constrain(Primary,Tertiary,Energy),
     retractall(initialtime(_)),
     statistics(runtime,[T,_]),
     assert(initialtime(T)), 
     ricerca(Tertiary,Energy,Soluzioni,Tempo,Tecnica).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% ESEMPIO DI PROBLEMA VINCOLATO con COSTO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vincola(Vars,Costo) :-
     Vars=[A,B,C,D,E,F],
     domain(Vars,-100,100),
     A-2*B+3*C-D*D+4*E-5*F #> 20,
     Costo #= 2*A - 3*B + 4*C - D - E + F. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Tre opzioni di ricerca: builtin, vincolata da noi,
%%% local search
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% In b_in NON CI INTERESSANO il tempo massimo ed il 
%%% numero di soluzioni trovate.

ricerca(Vars,Costo,_,_, b_in ):-
     labeling([minimize(Costo)],Vars),!,
     %%% output
     write('Ottimo: '),write(Costo),nl,
     write('Tempo totale di computazione: '),  
     statistics(runtime,[_,F]),write(F),write('ms'),nl.
     
ricerca(Vars,Costo,Sol,Tempo, b_up):-
    (var(Sol),Sol=1;nonvar(Sol)),        %%% default number of solutions = 1
    (var(Tempo),Tempo=1000;nonvar(Tempo)), %%% default time 1 second
     retractall(contate(_)),
     retractall(best(_,_)),
     assert(contate(0)),
     assert(best(1000000,_)),
     time_out(bottomup(Sol,Vars,Costo),Tempo,Out),
    (Out = time_out, nl, write('Warning: terminato per timeout'),nl;
     Out = success),
     write('*********************************************************'),nl,
     fail.

%%% In loc NON CI INTERESSA il 
%%% numero di soluzioni trovate.

ricerca(Vars,Costo,_Sol,Tempo, loc) :-
    (var(Tempo),Tempo=1000;nonvar(Tempo),true), %%% default time 1 second
     retractall(contate(_)),  assert(contate(0)),
     retractall(best(_,_)),   
     length(Vars,N),      
     length(Faked,N), 
     assert(best(1000000,Faked)),
     write('*** LNS search started '),nl,
     time_out(local(Vars,Costo),Tempo,_),
     write('*********************************************************'),nl,
     fail.

ricerca(_,_,_,_,_) :-
     best(1000000,_),!,
     write('Tempo insufficiente per la prima soluzione'),nl.

ricerca(_,_,_,_,_) :-
     retract(best(Val,_Sol)),
     write('Ottimo Trovato: '),write(Val),nl,
     write('Tempo totale di computazione: '),  
     statistics(runtime,[Tf,_]),
     initialtime(Ti),DT is Tf - Ti,write(DT),write('ms'),nl.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% BOTTOMUP (myconstraint)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
bottomup(0,_,_) :-
     write('Trovate tutte le soluzioni chieste'),nl.
bottomup(N,Vars,Costo) :-
     N > 0,
    (unasoluzione_bu(Vars,Costo), fail;
      N1 is N-1,
      bottomup(N1,Vars,Costo)).  
      
unasoluzione_bu(Vars,Costo) :-
     best(Lim,_),
     Costo #< Lim,  %%% <---------------
     write('Sto cercando sol. con Costo < '), write(Lim),nl,
     labeling([ffc,down],Vars), !, %%% Il cut ferma la ricerca di
     retract(best(_,_)),           %%% altre sol con lo stesso limite
     retract(contate(N)),M is N + 1,
     assert(contate(M)),
     %%%% STAMPE
     write('Trovata soluzione #'), write(M), write(' di Costo = '),write(Costo),
     write(' ---  In tempo: '), 
     statistics(runtime,[Tf,_]),
     initialtime(Ti),DT is Tf - Ti,write(DT),write('ms'),nl,
     assert(best(Costo,Vars)).  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% LOCAL SEARCH (local)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   
local(Vars,Costo) :-
     \+ altra_sol(Vars,Costo),
     local(Vars,Costo).

altra_sol(Vars,Costo) :-
     best(Lim,BestSol),
     Costo #< Lim, %%% Hill-climbing 
     random_partial_assignment(Vars,BestSol),
     %write('Sto cercando sol. con Costo < '), write(Lim),nl,
     %%% SE NON SIAMO NEL PRIMO CASO TROVIAMO L'OTTIMO LOCALE
     (ground(BestSol),labeling([ffc,down,minimize(Costo)],Vars);
      \+ ground(BestSol),
      write('Sto cercando la prima'),write(Lim),nl,
      labeling([ffc,down],Vars)),
     !,                            %%% Il cut ferma la ricerca di
     retract(best(_,_)),           %%% altre sol con lo stesso Lim
     retract(contate(N)),M is N + 1,
     assert(contate(M)),
     %%%% STAMPE
     write('Trovata soluzione #'), write(M), write(' di Costo = '),write(Costo),
     write(' ---  In tempo: '), 
     statistics(runtime,[Tf,_]),
     initialtime(Ti),DT is Tf - Ti,write(DT),write('ms'),nl,
     assert(best(Costo,Vars)),
     fail.  

%%%% Immagino di lasciare fissare ai valori precedenti
%%%  alcune variabili (le altre restano libere)

random_partial_assignment(Vars,BestSol) :-
     length(Vars,N),
     mask(N,Random),
     binds(Random,Vars,BestSol).

%%% Creo una maschera random di 0 e 1
%%% Il primo argomento di random bit sposta la percentuale di 1.
%%% Se 0 -> sempre 0, se 100 sempre 1, se 50 circa la meta' 

mask(0,[]) :- !.
mask(N,[A|R]) :- 
     randombit(50,A),
     N1 is N-1,
     mask(N1,R).

%%% La uso poi per un assegnamento parziale (random)    
     
binds([],[],[]).
binds([0|R],[_|Vars],[_|OldSol]) :-
     binds(R, Vars,OldSol).
binds([1|R],[X|Vars],[X|OldSol]) :-
     binds(R, Vars,OldSol).
                         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


constrain(Primary,Tertiary,Energy) :- 
       length(Primary,N),
       M is 2*N,
       length(Tertiary,M),
       Min is 0, %  N - integer(sqrt(N)), % 0
       Max is M, %  N + integer(sqrt(N)),   % 2*N
       domain(Tertiary,Min,Max),
       starting_point(Tertiary,N),
       avoid_self_loops(Tertiary),
       next_constraints(Tertiary),
       energy_constraint(Primary,Tertiary,Energy).
  
starting_point([N,N,N,N1|T],N) :-
      N1 is N + 1.
%      all_up([N,N,N,N1|T], Up, N, N),
%      lex_chain([[N,N,N,N1|T],Up]).   % Mats - break symmetry

avoid_self_loops(Tertiary) :-
       positions_to_integers(Tertiary, ListaInteri),
        all_distinct(ListaInteri).

positions_to_integers([X,Y|R], [I|S]):-
       I  #= X*100+Y,
       positions_to_integers(R,S).
positions_to_integers([],[]).

next_constraints([_,_]) .
next_constraints([X1,Y1,X2,Y2|C]) :-
       next(X1,Y1,X2,Y2),
       next_constraints([X2,Y2|C]).

next(X1,Y1,X2,Y2):-
%       domain([Dx,Dy],0,1),
%       Dx #= abs(X1-X2),
%       Dy #= abs(Y1-Y2),
%       Dx + Dy #= 1.
abs(X1-X2)+ abs(Y1-Y2) #= 1.


energy_constraint([_],_,0).
energy_constraint([A,B|Primary],[XA,YA,XB,YB|Tertiary],E) :-
       energy_contribution_of_A(0,A,XA,YA,Primary,Tertiary,E1),
       energy_constraint([B|Primary],[XB,YB|Tertiary],E2),
       E #= E1 + E2.

energy_contribution_of_A(_,_,_,_,[],[],0).
energy_contribution_of_A(0,A,XA,YA,[_|Primary],[_,_|Tertiary],E):-
       energy_contribution_of_A(1,A,XA,YA,Primary,Tertiary,E).
energy_contribution_of_A(1,A,XA,YA,[B|Primary],[XB,YB|Tertiary],E):-
       energy(A,XA,YA,B,XB,YB,C),
       energy_contribution_of_A(0,A,XA,YA,Primary,Tertiary,E1),
       E #= E1 + C.

energy(h,XA,YA,h,XB,YB,C) :-
       C in -1..0,
%       DX #= abs(XA - XB),
%       DY #= abs(YA - YB),
%       1 #= DX + DY   #<=> C #= -1.
 abs(XA - XB) + abs(YA - YB) #= 1  #<=> C #= -1.


energy(h,_,_,p,_,_,0).
energy(p,_,_,h,_,_,0).
energy(p,_,_,p,_,_,0).

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

all_up([], [], _, _).
all_up([_,_|Tert], [X,Y|Up], X, Y) :-
    Y1 is Y+1,
    all_up(Tert, Up, X, Y1).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% RANDOM PRIMITIVES (reimplemented)

randombit(Perc, R) :-
  (random > (100-Perc)*16383//100, !, R=0;
   R=1).

random_subseq([],[]).
random_subseq([A|R],T) :-
   random_subseq(R,S),
   randombit(F),
   (F==1 -> T = [A|S];
    F==0 -> T = S).

%%% AUXILIARY: 
%%%    used to find the maximum value returned 
%%%    by random (not found in the manual)     
%%% Apparently 32767
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

metaran(L,T) :-
  metaran(1,L,T).

metaran(N,L,T) :-
   (random > L, T=N;
    N1 is N+1, metaran(N1,L,T)). 

%%% To check "randombit":

manyrandom(N) :-
   manyrandom(N,L),
   %write(L),nl,
   sum(L,#=,Ones),
   format("Number ~d of 1 in the list of ~d elements ~n",[Ones,L]).
   
manyrandom(0,[]) :-!.
manyrandom(N,[A|L]) :-
    randombit(50,A),
    N1 is N-1,
    manyrandom(N1,L).

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