function [rte,TC] = vrpexchange(rte,C,varargin)
%VRPEXCHANGE Two-vertex exchange procedure for VRP improvement.
% [rte,TC] = vrpexchange(rte,C,cap,twin,rtefeas,h)
% Two vertices from different routes are simultaneousely exchanged; e.g.,
% vertex m from route i is exchanged with vertex n from route j:
%             rte{i} = [rte{i}(1:m-1) rte{i}(n) rte{i}(m+1:end)]
%             rte{j} = [rte{j}(1:n-1) rte{i}(m) rte{j}(n+1:end)]
%    rte = m-element cell array of m routes (can be infeasible)
%      C = n x n matrix of costs between n vertices
%    cap = {q,Q} = cell array of capacity arguments, where
%              q = n-element vector of vertex demands, with depot q(1) = 0
%              Q = maximum route load
%   twin = {ld,TW} = cell array of time window arguments
%             ld = 1, n, or (n+1)-element vector of load/unload timespans
%             TW = n or (n+1) x 2 matrix of time windows
%                = (n+1)-element cell array, if multiple windows
%rtefeas = {'rtefeasfun',P1,P2,...} = cell array specifying user-defined
%          function to test the feasibility of a single route
%      h = (optional) handle to vertex plot, e.g., h = PPLOT(XY,'.');
%          use 'Set Pause Time' on figure's Matlog menu to pause plot
%     TC = m-element vector of route costs
%
% See RTETC for information about the input parameters
%
% Example:
% vrpnc1  % Loads XY, q, Q, ld, and maxTC
% C = dists(XY,XY,2);
% h = pplot(XY,'r.');
% pplot(XY,num2cellstr(1:size(XY,1)))
% [rte,TC] = vrpinsert(C,{q,Q},{ld},{'maxTCfeas',maxTC},[]);
% [rte,TC] = vrpexchange(rte,C,{q,Q},{ld},{'maxTCfeas',maxTC},h);

% Copyright (c) 1994-2006 by Michael G. Kay
% Matlog Version 9 13-Jan-2006 (http://www.ie.ncsu.edu/kay/matlog)

% Input Error Checking ****************************************************
error(nargchk(1,6,nargin))

if length(varargin) < 4, [varargin{length(varargin)+1:4}] = deal([]); end
[cap,twin,rtefeas,h] = deal(varargin{:});

try  % Use for error checking and to store input arguments
   [rte,TC] = tsp2opt(rte,C,cap,twin,rtefeas);
catch
   errstr = lasterr;
   idx = find(double(errstr) == 10);
   error(errstr(idx(1)+1:end))
end

if isempty(h), h = NaN; end

if (~ishandle(h) & ~isnan(h)) | (ishandle(h) & (length(h) ~= 1 | ...
      ~strcmp(get(h,'Type'),'line') | ...
      ~strcmp(get(h,'LineStyle'),'none') | ...
      length(get(h,'XData')) ~= size(C,1) | ...
      length(get(h,'YData')) ~= size(C,1)))
   error('Invalid handle "h".')
end
% End (Input Error Checking) **********************************************

str = sprintf('VRP Exchange: 2-Opt Improvement, TC = %f',sum(TC));
fprintf('%s\n',str)

if ishandle(h)
   axes(get(h,'Parent'))
   title(str)
   delete(findobj(gca,'Tag','rteplot'))
   XY = [get(h,'XData')' get(h,'YData')'];
   if strcmp(get(gca,'Tag'),'proj'), XY = invproj(XY); end
   pplot(rte,XY,'m','Tag','rteplot')
   pplot(XY(1,:),'rs','Tag','rteplot')
   if isempty(pauseplot('get')), pauseplot(Inf), end
end

done = 0;
Done = logical(zeros(length(rte)));
while ~done
   done = 1;
   for i = 1:length(rte)-1
      for j = i+1:length(rte)
         if Done(i,j), continue, end  % Skip if no change
         for m = 2:length(rte{i})-1
            for n = 2:length(rte{j})-1
               rteim = [rte{i}(1:m-1) rte{j}(n) rte{i}(m+1:end)];
               [rteim,TCim] = tsp2opt(rteim);
               if isfinite(TCim)
                  rtejn = [rte{j}(1:n-1) rte{i}(m) rte{j}(n+1:end)];
                  [rtejn,TCjn] = tsp2opt(rtejn);
                  if TCim + TCjn < TC(i) + TC(j)
                     rte{i} = rteim; rte{j} = rtejn;
                     TC(i) = TCim; TC(j) = TCjn;
                     done = 0;
                     Done([i j],:) = 0; Done(:,[i j]) = 0;
                     
                     str = sprintf(['VRP Exchange: Rte(Vtx) %d(%d) ' ...
                           'and %d(%d), TC = %f'],...
                           i,rte{i}(m),j,rte{j}(n),sum(TC));
                     fprintf('%s\n',str)
                     
                     if ishandle(h)
                        title(str)
                        delete(findobj(gca,'Tag','rteplot'))
                        pplot(rte,XY,'m','Tag','rteplot')
                        pauseplot
                     end
                     
                     break
                  end
               end
            end % n
            if ~done, break, end
         end % m
         if done, Done(i,j) = 1; end  % No improvement
         if ~done, break, end
      end % j
      if ~done, break, end
   end % i
end

if ishandle(h)
   str = sprintf('VRP Exchange: Final TC = %f and %d Routes',...
      sum(TC),length(TC));
   fprintf('%s\n\n',str)
   title(str)
   delete(findobj(gca,'Tag','rteplot'))
   pplot(XY(1,:),'rs','Tag','rteplot')
   pplot(rte,XY,'m','Tag','rteplot')
   pauseplot
end