%
% TSP - Traveling Salesman Problem with time windows
%
% Given N cities, a matrix of distances between any two of them,
% and some requirements for when each city can be visited,
% what is a shortest round trip through all N cities?
%
% ?- solve.
% Found a solution with cost 3952
% Found a solution with cost 3753
% Found a solution with cost 3613
% Found a solution with cost 3537
% Found a solution with cost 3473
% Found a solution with cost 3449
% Found a solution with cost 3431
% Found a solution with cost 3359
% Found a solution with cost 3346
% Found a solution with cost 3336
% Found a solution with cost 3323
% Found no solution with cost 2022.0 .. 3322.0
% 3323 : [](2, 14, 4, 5, 6, 12, 13, 11, 10, 1, 9, 7, 8, 3)
% times : [](3323, 153, 740, 1029, 1520, 1920, 2102, 2499, 2675, 2951, 2632, 1939, 2226, 529)
% Yes (0.80s cpu)
% 
% Author: Joachim Schimpf, 2021
% Creative Commons Attribution 4.0 International License.
%

% Use the interval constraint solver and generic B&B
:- lib(ic).
:- lib(branch_and_bound).


% The basic TSP model takes as input a distance matrix Dist, and
% maintains a successor array Next, and a distance array Legs which
% contains the distance travelled from each node to its successor.
% To enable time window handling, we add an array Time giving the arrival
% time at each node, which can then be constrained according to requirements.
% Arrival time at the successor of node I can be computed as arrival time
% at I plus travel time from I to its successor (for simplicity, assume
% that distance = time, and that we start at node 1 at time 0). 
% We then minimise Cost, which is the total distance travelled.
tsp_tw(Dist, Next, Time, Cost) :-
        dim(Dist, [N,N]),               % get distance matrix size

        dim(Next, [N]),                 % successor variables
        Next #:: 1..N,
        circuit(Next),                  % they must form a circuit

        dim(Legs, [N]),                 % Legs[I] is distance I->Next[I]
        dim(Time, [N]),                 % Time[I] is arrival time at I
        ( for(I,1,N), param(Dist,Next,Legs,Time) do
            element(I, Next, J),
            element(J, Dist[I], DistIJ),
            element(I, Legs, DistIJ),

            ( I==1 -> TimeAtI = 0 ; element(I, Time, TimeAtI) ),
            element(J, Time, TimeAtJ),
            TimeAtJ #= TimeAtI + DistIJ
        ),

        % ... here we could further constrain Time[] ...

        Cost #= sum(Legs),              % total distance travelled
                                        % search and minimize
        minimize(search(Legs,0,smallest,indomain_min,complete,[]), Cost).


% Solve an example
solve :-
        data(burma14, Dist),
        tsp_tw(Dist, Next, Time, Cost),
        writeln(Cost:Next),
        writeln(times:Time).

data(burma14, [](
        [](0, 153, 510, 706, 966, 581, 455, 70, 160, 372, 157, 567, 342, 398),
        [](153, 0, 422, 664, 997, 598, 507, 197, 311, 479, 310, 581, 417, 376),
        [](510, 422, 0, 289, 744, 390, 437, 491, 645, 880, 618, 374, 455, 211),
        [](706, 664, 289, 0, 491, 265, 410, 664, 804, 1070, 768, 259, 499, 310),
        [](966, 997, 744, 491, 0, 400, 514, 902, 990, 1261, 947, 418, 635, 636),
        [](581, 598, 390, 265, 400, 0, 168, 522, 634, 910, 593, 19, 284, 239),
        [](455, 507, 437, 410, 514, 168, 0, 389, 482, 757, 439, 163, 124, 232),
        [](70, 197, 491, 664, 902, 522, 389, 0, 154, 406, 133, 508, 273, 355),
        [](160, 311, 645, 804, 990, 634, 482, 154, 0, 276, 43, 623, 358, 498),
        [](372, 479, 880, 1070, 1261, 910, 757, 406, 276, 0, 318, 898, 633, 761),
        [](157, 310, 618, 768, 947, 593, 439, 133, 43, 318, 0, 582, 315, 464),
        [](567, 581, 374, 259, 418, 19, 163, 508, 623, 898, 582, 0, 275, 221),
        [](342, 417, 455, 499, 635, 284, 124, 273, 358, 633, 315, 275, 0, 247),
        [](398, 376, 211, 310, 636, 239, 232, 355, 498, 761, 464, 221, 247, 0)
    )
).