#########################################################################
## SCRIPT FILE:   Maya the Traveler:  Kyiv, Paris, Rome, Vaasa. 
#########################################################################

## The distance matrix C, and its vector form c augmented with u:
M = 10000;                              
C = [   M   2023   1674   1497;
     2023      M   1105   1967;
     1674   1105      M   2429;
     1497   1967   2429      M];
c = [C'(:); 0; 0; 0];

## The first 4*4 columns are for the tour matrix X in the vector form
## and the last 3 columns are for the auxiliary vector u:
A = [ 1  0  0  0  1  0  0  0  1  0  0  0  1  0  0  0    0  0  0;
      0  1  0  0  0  1  0  0  0  1  0  0  0  1  0  0    0  0  0;
      0  0  1  0  0  0  1  0  0  0  1  0  0  0  1  0    0  0  0;
      0  0  0  1  0  0  0  1  0  0  0  1  0  0  0  1    0  0  0;

      1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0    0  0  0;
      0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0    0  0  0;
      0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0    0  0  0;
      0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1    0  0  0;

      0  0  0  0  0  0  4  0  0  0  0  0  0  0  0  0    1 -1  0;
      0  0  0  0  0  0  0  4  0  0  0  0  0  0  0  0    1  0 -1;
      0  0  0  0  0  0  0  0  0  4  0  0  0  0  0  0   -1  1  0;
      0  0  0  0  0  0  0  0  0  0  0  4  0  0  0  0    0  1 -1;
      0  0  0  0  0  0  0  0  0  0  0  0  0  4  0  0   -1  0  1;
      0  0  0  0  0  0  0  0  0  0  0  0  0  0  4  0    0 -1  1];
##   11 12 13 14 21 22 23 24 31 32 33 34 41 42 43 44     

## The bounds and other glpk parameters, and finally call glpk:
L = 50;                                         # For u upper bounds.
b = [1 1 1 1 1 1 1 1 3 3 3 3 3 3]';
lb = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
ub = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 L L L]';
ctype = "SSSSSSSSUUUUUU";
vtype = "IIIIIIIIIIIIIIIICCC";
param.msglev=0;
[x_min,tot_cost,stat] = glpk(c,A,b,lb,ub,ctype,vtype,1,param);

## Get tour.  A tour always starts from the city 1 (why not?).
X = reshape(x_min(1:16),4,4);                   # x_min to matrix form.
tour = [1];                                     # Start from city 1.
for k=1:3                                       # Find the next 3 cities.
  tour = [tour; find(X(tour(k),:))];            # Append city from row k.
endfor

## Print the outputs
tot_cost
tour
