%
% Travelling Salesman
%
/* 
 needs data of the form and distance(Id1, Id2, D).
Using tsUSA.pl :

Cities    Length    Tour
 5         120       [1, 2, 5, 4, 3, 1]
 6         141       [1, 2, 5, 6, 3, 4, 1]
 7         142       [1, 2, 5, 6, 7, 3, 4, 1]
 8         151       [1, 2, 5, 6, 7, 8, 3, 4, 1]
 9         161       [1, 2, 5, 6, 7, 8, 9, 3, 4, 1]
10         200       [1, 2, 5, 6, 7, 8, 10, 9, 3, 4, 1]
11         244       [1, 2, 5, 6, 7, 8, 10, 11, 9, 3, 4, 1]
12         258       [1, 2, 5, 6, 7, 8, 10, 12, 11, 9, 3, 4, 1]
13         335       [1, 2, 5, 6, 7, 8, 10, 13, 12, 11, 9, 3, 4, 1]
14         406       [1, 2, 5, 6, 7, 8, 10, 14, 13, 12, 11, 9, 3, 4, 1]

*/
%
% "Naive" Prolog version
%
salesman_P([P1|Ps], [P1|Tour], Length) :-
	min_search(perm_lengthP(Ps, Tour, P1, P1, 0, Length), Length).
	
perm_lengthP([], [Last], Current, Last, Dist, Length) :-
	distance(Current, Last, D),
	Length is Dist + D.
perm_lengthP(List, [X|Xs], Current, Last, Dist, Length) :-
	select(X, List, Newlist),  % from library(lists)
	distance(Current, X, D),
	Newdist is Dist + D,
	perm_lengthP(Newlist, Xs, X, Last, Newdist, Length).
	

min_search(Goal,Objective) :-
	current_prolog_flag(max_tagged_integer,MaxInt),
	nb_setval('$min_search',[MaxInt,Goal]),  % initialize
	Goal,                                    % generate solution(s)
	nb_getval('$min_search',[Current,_]),
	Objective < Current,                     % if better than current best, save it
	nb_setval('$min_search',[Objective,Goal]),
	fail.
min_search(Goal,Objective) :-
	nb_getval('$min_search',[Objective,Goal]).	% unify goal with best solution

%
% K contraint version
%
salesman_K([P1|Ps], [P1|Tour], Length) :-
	Length::integer(0,_),
	min_ratchet(perm_lengthK(Ps, Tour, P1, P1, 0, Length), Length).

perm_lengthK([], [Last], Current, Last, Dist, Length) :-
	distance(Current, Last, D),  % Last = C1
	Length is Dist + D.
perm_lengthK(List, [X|Xs], Current, Last, Dist, Length) :-
	select(X, List, Newlist),  % from library(lists)
	distance(Current, X, D),
	Newdist is Dist + D, {Newdist =< Length},  % fail when Newdist exceeds Length
	perm_lengthK(Newlist, Xs, X, Last, Newdist, Length).

%
% Find the lowest value of Objective generated by Goal subject to Constraint
%
min_ratchet(Goal,Objective) :- 
	{Objective <= Bound},              % define "bound" constraint (must be inclusion to avoid back flow)
	Goal,                              % generate a solution
	nb_setval('$min_ratchet', Goal),   % save solution
	range(Objective,[_,UB]),
	range(Bound,[LB,_]),
	NUB is UB-1,                       % assume integer values
	nb_setbounds(Bound,[LB,NUB]),      % set new bound
	fail.                              % and backtrack
min_ratchet(Goal,Objective) :- 
	catch(nb_getval('$min_ratchet', Goal),_,fail),  % retrieve solution
	nb_delete('$min_ratchet').
	
%
% KT contraint version
%
salesman_KT(Pts, [P1|Tour], Length) :-
	length(Pts,Count),                 % number of points
	length(Distances,Count),           % list of distances
	symbolic_sum(Distances,SumD),      % SumD = symbolic sum of distances
	[Length,Distances]::integer(0,_),  % all are positive integers
	{Length==SumD},                    % constrain Length to be  sum of distances
	Pts = [P1|Ps],
	min_ratchet(perm_lengthKT(Ps, Tour, P1, P1, Distances), Length).

symbolic_sum([X],X) :- !.
symbolic_sum([X|Xs], X+Sum) :-	symbolic_sum(Xs,Sum).

perm_lengthKT([], [Last], Current, Last, [D]) :-
	distance(Current, Last, D).
perm_lengthKT(List, [X|Xs], Current, Last, [D|Ds]) :-
	select(X, List, Newlist),  % from library(lists)
	distance(Current, X, D),
	perm_lengthKT(Newlist, Xs, X, Last, Ds).
	
%
% Leg based version
%
salesman_L(Pts, Tour, Length) :-
	length(Pts,Count),
	maplist(make_pt, Pts, Pointlist),       % make `pt` terms from points list
	leg_list(Pointlist, Leglist), 
	sort(Leglist, Ascending),               % on leg length (distance)
	setup_leg_constraints(Pts, Leglist),    % 2 legs to any point
	Length::integer(0,_),                   % Tour length (distance)
	min_ratchet(select_legs(Count, Ascending, Selected, 0, Length),Length),
	Pts = [P1|_],                           % sequence the tour from first point
	sequence_legs(Selected, Tour, P1).

% add connected var C to endpoint
make_pt(N, pt(N,Cn)).

% form the list of N*(N-1)/2 legs
leg_list([], []).
leg_list([C|Cs], DL) :- 
	leg_from(Cs, C, DL, EL),
	leg_list(Cs, EL).
	
% legs from list of points to a point
leg_from([], P, L, L).
leg_from([X|Xs], Y, [leg(D,P,Acc,X,Y)|Ls], E) :-
	dist(X,Y,D),              % distance between X  and Y
	P::boolean,               % boolean variable, true if leg part of tour
	leg_from(Xs, Y, Ls, E).

dist(pt(P1,_), pt(P2,_), D) :- 	distance(P1,P2,D).

% constrain number of enabled legs to any one point to be 2

setup_leg_constraints([], Legs).
setup_leg_constraints([P|Ps], Legs) :-
	 % only two legs including N are allowed (total degree in tour is exactly 2)
	incident(Legs, P, PNs),  % list of booleans of legs terminating at point P
	symbolic_sum(PNs, S), {S == 2},
	setup_leg_constraints(Ps, Legs).

incident([], _, []).
incident([X|Xs], N, [P|Ys]) :- 
	incid(X,N,P), !,  % point N is an endpoint of leg X, P is leg enabled 
	incident(Xs,N,Ys).
incident([X|Xs], N, Ys) :- 
	incident(Xs, N, Ys).

incid(leg(D,P,Acc, pt(N,_),_), N, P).
incid(leg(D,P,Acc, _,pt(N,_)), N, P).

% select legs for the tour
select_legs(1, [leg(D,1,Acc,X,Y)|_ ], [leg(D,X,Y)], Length, Total) :- !,  % last leg
	Total is Length+D.
select_legs(N, [leg(D,S,Acc,X,Y)|Ls], [leg(D,X,Y)|Rs], Length, Total) :- 
	connect_(X,Y),                  % connect X & Y 
	S=1,                            % mark leg as selected
	NewLength is Length+D,          % accumulate length
	estimate(N,Ls,NewLength,Est),   % best estimate for remaining
	{Est =< Total},                 % must be =< Total
	N1 is N - 1,                    % number of remaining legs
	select_legs(N1, Ls, Rs, NewLength, Total).
select_legs(N, [leg(_,0,_,_,_)|Ls], R, Length, Total) :-  % else skip leg, force boolean false
	select_legs(N, Ls, R, Length, Total).

connect_(pt(_,C1),pt(_,C2)) :- 
	C1 \== C2,  % not already connected
	C1=C2.      % mark as connected (same)
	
estimate(1,_, D,D) :- !.             % N=1, then done
estimate(N, [leg(D,S,_,_,_)|Ls], D1, Est) :- S \== 0, !, % if not unselectable	
	DD is D1 + D,                                        % accumulate
	N1 is N - 1,
	estimate(N1, Ls, DD, Est).                           % estimate rest
estimate(N, [_|Ls], D, Est) :-                           % else skip leg
	estimate(N, Ls, D, Est).

% sequence selected legs into a list of point values	
sequence_legs([], [P], P) :- !.     % last leg returns to origin ??
sequence_legs(List, [P|Ps], P) :-
	delete_leg(List, leg(D, pt(P,_),pt(Next,_)), Rest), !,
	sequence_legs(Rest, Ps, Next).
	
delete_leg([X|Xs], Y, Xs) :- 
	match_leg(X,Y), !.
delete_leg([X|Xs], Y, [X|Ys]) :- 
	delete_leg(Xs,Y,Ys).
	
match_leg(leg(D,A,B), leg(D,A,B)).
match_leg(leg(D,A,B), leg(D,B,A)).

%
% per leg accumulator
%
salesman_LA(Pts, Tour, Length) :-
	length(Pts,Count),
	maplist(make_pt, Pts, Pointlist),
	leg_list(Pointlist, Leglist),           % create Leglist
	sort(Leglist, Ascending),               % and sort by distance
	setup_leg_constraints(Pts, Leglist),
	acc_constraints(Ascending, Length),
	crosses_constraints(Ascending),
	min_ratchet(select_legsA(Count, Ascending, Selected),Length),
	Pts = [P1|_],                           % start from first point
	sequence_legs(Selected, Tour, P1).

acc_constraints([], 0).
acc_constraints([leg(D,P,Acc,C1,C2)|Ls], Acc) :-
	Acc::integer(0,_),
	{Acc == P*D+Acc1},
	acc_constraints(Ls, Acc1).

crosses_constraints([]).
crosses_constraints([L0|Ls]) :-
	crosses_constraints(Ls,L0),
	crosses_constraints(Ls).
	
crosses_constraints([],_).
crosses_constraints([L|Ls],L0) :-
	cross_constrain(L,L0),
	crosses_constraints(Ls,L0).

cross_constrain(leg(_,S0,_,pt(P0,_),pt(P1,_)),leg(_,S1,_,pt(P2,_),pt(P3,_))) :-
	\+ adjacent(P0,P1,P2,P3),              % legs don't share an endpoint
	point(P0,Cs0), point_coordinates(Cs0,X0,Y0), 
	point(P1,Cs1), point_coordinates(Cs1,X1,Y1),
	point(P2,Cs2), point_coordinates(Cs2,X2,Y2),
	point(P3,Cs3), point_coordinates(Cs3,X3,Y3),
	\+ disjoint_axis(X0,X1,X2,X3),         % optimization: X overlap?
	\+ disjoint_axis(Y0,Y1,Y2,Y3),         % optimization: Y overlap?
	[A,B]::real(0,1),
	{ X0 + A*(X1-X0) == X2 + B*(X3-X2) },  % X intersection 
	{ Y0 + A*(Y1-Y0) == Y2 + B*(Y3-Y2) },  % Y intersection
	!,                                     % success, apply constraint
	{S0 + S1 =< 1}.                        % if legs cross, both can't be selected
cross_constrain(_,_).                      % no constraints if legs don't cross

point_coordinates(geo(X,Y),X,Y).

adjacent(P,_,P,_).
adjacent(P,_,_,P).
adjacent(_,P,P,_).
adjacent(_,P,_,P).	

% basic overlap detection
disjoint_axis(V0,V1,V2,V3) :- max(V0,V1) < min(V2,V3).
disjoint_axis(V0,V1,V2,V3) :- min(V0,V1) > max(V2,V3).

select_legsA(1, [leg(D,1,Acc,P,Q)|_], [leg(D,P,Q)] ) :- !.
select_legsA(N, [leg(D,S,Acc,P,Q)|Ls], R):-
	%% constrain Acc if leg selectable but not yet determined
	(var(S) -> (estimate(N, Ls, D, TD), {Acc >= TD}) ; true),  % connect estimate to net
	sel(S, leg(D,P,Q), N, N1, R,R1),  % ¿leg select?
	select_legsA(N1, Ls, R1).

sel(S, leg(D,P,Q), N, N1,[leg(D,P,Q)|Ls],Ls) :-
	connect_(P,Q), S=1,  % quick connect test before setting boolean
	N1 is N - 1.
sel(0, Leg, N,N, L,L).   % Leg not selected

%
% support for random tour tests
%
tour_list(ListId,Length,Points) :- 
	length(Points,Length), 
	tour_list(ListId,Ps), 
	append(Points,_,Ps).

tour_list(1,[35, 37, 19, 29, 17, 32, 14, 38, 6, 36, 22, 8, 16, 5, 3, 7, 34, 40, 15, 31]).
tour_list(2,[13, 3, 10, 7, 24, 4, 1, 8, 25, 11, 20, 14, 29, 5, 16, 30, 40, 31, 38, 15]).
tour_list(3,[11, 19, 8, 38, 40, 17, 1, 39, 23, 3, 13, 21, 41, 10, 30, 14, 28, 26, 37, 22]).
