/*

  Benchmarking TSP problem in Picat.


  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import cp,util,benchmark_labeling.

main => go.

go=>        
  tsp_test(glpk, Matrix),
  % tsp_test(chip, Matrix),
  % selection(VariableSelect),
  % choice(ValueSelect),
  % testing a smaller grid
  VariableSelect = [ff,ffc,ffd,forward],
  ValueSelect = [split,up],
  Timeout = 3000, % millis
  Goal = $tsp(Matrix, Cities, Costs,Cost),
  benchmark(VariableSelect,ValueSelect,Timeout,Goal),
  nl.


%%
%% Ah, the found solution in Cities,Costs, and Cost "caches" in furher tries!!!!
%%

% TSP using a matrix, element/1 and circuit/1 constraints.
tsp(Matrix,Cities, Costs,Cost,VarSel,ValSel) =>


   % println(matrix=Matrix),

   Len = Matrix.length,
   Cities = new_list(Len),
   Cities :: 1..Len,

   % calculate upper and lower bounds of the Costs list
   Dists = [Matrix[I,J] : I in 1..Len, J in 1..Len, Matrix[I,J] > 0],
   MinDist = min(Dists),
   MaxDist = max(Dists),

   % println([mindist=MinDist,maxDist=MaxDist]),
   Costs = new_list(Len),
   Costs :: MinDist..MaxDist,

   Cost #= sum(Costs),
   Cost #> 0,

   all_different(Cities), % might give a speedup
   circuit(Cities),

   % connect cities and costs 
   foreach({Row,City,C} in zip(Matrix,Cities,Costs))
     element(City,Row,C)
   end,

   % println(cost=Cost), % show the inferred domain

   % Vars = Cities ++ Costs, % slower
   Vars = Costs ++ Cities,

   % solve($[min(Cost),report(println(cost=Cost))],Vars), % 3.7s/175883 backtracks for the GLPK problem
   % solve($[ffd,split,min(Cost),report(println(cost=Cost))],Vars), % 3.7s/0 backtracks for the GLPK problem
   % solve($[forward,min(Cost),report(println(cost=Cost))],Vars), % 2.4s/0 backtracks for the GLPK problme

   % solve($[forward,up,min(Cost),report(println(cost=Cost))],Vars), % 3.7s/0 backtracks for the GLPK problme

   % solve($[forward,split,min(Cost),report(println(cost=Cost))],Vars), % 2.3s/0 backtracks for the GLPK problem
   % solve($[split,min(Cost),report(println(cost=Cost))],Vars), % s/0 backtracks for the GLPK problem
   solve($[ValSel,ValSel,min(Cost)],Vars).




% This problem is from 
% GLPK:s example tsp.mod
% (via http://www.hakank.org/minizinc/tsp.mzn)
% """
% These data correspond to the symmetric instance ulysses16 from:
% Reinelt, G.: TSPLIB - A travelling salesman problem library.
% ORSA-Journal of the Computing 3 (1991) 376-84;
% http://elib.zib.de/pub/Packages/mp-testdata/tsp/tsplib 
% 
% The optimal solution is 6859
% """
tsp_test(glpk, Matrix) =>
  Matrix = 
  [[0,509,501,312,1019,736,656,60,1039,726,2314,479,448,479,619,150],
   [509,0,126,474,1526,1226,1133,532,1449,1122,2789,958,941,978,1127,542],
   [501,126,0,541,1516,1184,1084,536,1371,1045,2728,913,904,946,1115,499],
   [312,474,541,0,1157,980,919,271,1333,1029,2553,751,704,720,783,455],
   [1019,1526,1516,1157,0,478,583,996,858,855,1504,677,651,600,401,1033],
   [736,1226,1184,980,478,0,115,740,470,379,1581,271,289,261,308,687],
   [656,1133,1084,919,583,115,0,667,455,288,1661,177,216,207,343,592],
   [60,532,536,271,996,740,667,0,1066,759,2320,493,454,479,598,206],
   [1039,1449,1371,1333,858,470,455,1066,0,328,1387,591,650,656,776,933],
   [726,1122,1045,1029,855,379,288,759,328,0,1697,333,400,427,622,610],
   [2314,2789,2728,2553,1504,1581,1661,2320,1387,1697,0,1838,1868,1841,1789,2248],
   [479,958,913,751,677,271,177,493,591,333,1838,0,68,105,336,417],
   [448,941,904,704,651,289,216,454,650,400,1868,68,0,52,287,406],
   [479,978,946,720,600,261,207,479,656,427,1841,105,52,0,237,449],
   [619,1127,1115,783,401,308,343,598,776,622,1789,336,287,237,0,636],
   [150,542,499,455,1033,687,592,206,933,610,2248,417,406,449,636,0]
  ].

%
% This problem is from the SICStus example 
% ./library/clpfd/examples/tsp.pl
% The "chip" examples 
%
tsp_test(chip,Matrix) =>
 Matrix = 
       [[0,205,677,581,461,878,345],
        [205,0,882,427,390,1105,540],
        [677,882,0,619,316,201,470],
        [581,427,619,0,412,592,570],
        [461,390,316,412,0,517,190],
        [878,1105,201,592,517,0,691],
        [345,540,470,570,190,691,0]].


generate_random(N,MaxVal) = Matrix =>
   _ = random2(),
   Matrix1 = { {1 + random() mod MaxVal : _ in 1..N} : _ in 1..N},
   foreach(I in 1..N)
     Matrix1[I,I] := 0
   end,
   Matrix = Matrix1.array_matrix_to_list_matrix().