function noder = graph_dist_pairing_greedy ( maxit, nodeb, noder, nnode, ...
  tol, xb, xr, yb, yr )

%*****************************************************************************80
%
%% graph_dist_pairing_greedy() pairs two sets of nodes using the least total distance.
%
%  Discussion:
%
%    The method is iterative, and is not guaranteed to find the best
%    possible arrangement.  This is particulary true because it is a
%    "local" method, which only considers pairwise switches of the
%    red nodes that reduce the total distance.  This means that a
%    "locally minimizing" pairing might be found which is not the
%    global minimizer.
%
%    On the other hand, in the absence of a theoretical plan for how
%    to reach the global minimizer, the brute force search would
%    require that EVERY possible pairing be considered, and its total
%    distance computed.  This means that a total of factorial(NNODE)
%    graphs would have to be generated.
%
%    The approach used here, on each iterative step, looks at a
%    maximum of NNODE * (NNODE-1) graphs, which represents a
%    significantly more efficient method.
%
%    It would not be hard to extend this approach to a method which
%    considers switches of THREE red nodes at a time, though the
%    work there involve looking at NNODE * (NNODE-1) * (NNODE-2)
%    graphs, and as we increase the number of graphs we examine,
%    we begin to approach the NNODE% rate for the brute force
%    algorithm.
%
%    It also would not be hard to extend this method to a case where
%    there are three sets of nodes, arranged in triples, and again
%    the total distance is to be minimized.
%
%
%    If it is suspected that the pairing returned is only
%    a local minimizer, then the user is advised to restart the
%    calculation after randomly permuting the entries of NODER, so that
%    the routine starts from a different point in the space of graphs.
%
%    The routine is given:
%
%      an initial ordering of the black and red nodes, so that
%      ( NODEB(I), NODER(I) ) represents the I-th pair,
%
%      the X and Y coordinates of the black and red nodes,
%
%      a maximum number of iterations, and a relative distance
%      decrease requirement,
%
%    and computes:
%
%      a new ordering of the red nodes, contained in NODER, which should
%      reduce the total distance between corresponding red and black
%      nodes.
%
%    The approach can be applied to a variety of problems including:
%
%    1) We are given two sets of NNODE points, which we will call the
%       "red" and "black " groups, and the (X,Y) coordinates of each
%       point.  We may imagine these points as forming the two sets of
%       nodes of a bipartite graph lying in the (X,Y) plane.  We wish
%       to choose a pairing of red and black nodes which results in
%       the shortest total arc length.
%
%    2) We are given two sets of NNODE complex quantities, which we
%       believe are approximations to the same (unknown) set of
%       quantities.  We wish to arrange this data into NNODE pairs,
%       each containing a unique element from each set of data, which
%       minimizes the sum of squares of the discrepancies between the
%       pairs.  In particular, the two sets of data might be two
%       separate estimates of the complex eigenvalues of a matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    05 November 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer MAXIT, the maximum number of iterations
%    allowed.  Each iteration considers, one at a time, each black node, and
%    seeks to switch its red neighbor for another red neighbor that
%    reduces the total distance.
%
%    integer NODEB(NNODE), the "labels" of the black
%    nodes.  You probably want to just set NODEB(I) = I, for i = 1 to NNODE.
%    The entries in NODEB will not be changed.
%
%    integer NODER(NNODE), the initial "labels" of the red nodes.
%    You probably want to just set the input value of NODER(I) = I,
%    for i = 1 to NNODE.  The values of ( NODEB(I), NODER(I) ) contain the
%    labels of the I-th pair of black and red nodes.
%
%    integer NNODE, the number of nodes in the black and red sets.
%
%    real TOL.
%    TOL is the relative decrease that the user demands in the
%    total distance, after each iterative step.  If we denote
%    the distance before the iterative step as OLDTOT, and the
%    distance after the iterative step as TOTAL, then the
%    routine will try another iterative step as long as "enough"
%    progress was made on this step.  Enough progress was made
%    whenever OLDTOT - TOTAL < TOL * TOTAL
%
%    real XB(NNODE), the X coordinates of the black nodes.
%
%    real XR(NNODE), the X coordinates of the red nodes.
%
%    real YB(NNODE), the Y coordinates of the black nodes.
%
%    real YR(NNODE), the Y coordinates of the red nodes.
%
%  Output:
%
%    integer NODER(NNODE), the updated "labels" of the red nodes.
%    If a better pairing of the nodes has been found,
%    this will be reflected in the newly permuted values of NODER.
%

%
%  Compute the total distance of the starting pairing.
%
  total = 0.0;

  for indx = 1 : nnode

    nodeb1 = nodeb(indx);
    noder1 = noder(indx);

    total = total + sqrt ( ...
      ( xb(nodeb1) - xr(noder1) )^2 + ( yb(nodeb1) - yr(noder1) )^2 );

  end

  fprintf ( 1, '\n' );
%
%  Begin the iterations.
%
  for it = 1 : maxit

    if ( total == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'graph_dist_pairing_greedy(): Early termination.\n' );
      fprintf ( 1, '  Total discrepancy is low enough.\n' );
      return
    end
%
%  Save the current total distance for comparison at the end of the
%  iteration.
%
    oldtot = total;
    nswap = 0;
%
%  Consider each black node, by running through indices INDX1 = 1
%  through NNODE of the NODEB array.
%
    for indx1 = 1 : nnode
%
%  Get the labels of the current INDX1-th pair of black and red nodes.
%
      nodeb1 = nodeb(indx1);
      noder1 = noder(indx1);
%
%  Look at the black node with INDX2 = 1 through NNODE, but ignore
%  the case where INDX1 = INDX2.
%
      for indx2 = 1 : nnode
%
%  Get the labels of the current INDX2-th pair of black and red nodes.
%
        nodeb2 = nodeb(indx2);
        noder2 = noder(indx2);

        if ( indx2 ~= indx1 )
%
%  Compute the total distance between (NODEB1,NODER1) and
%  (NODEB2,NODER2), and compare it to the total where we switch the
%  red nodes.
%
          dist1 = sqrt ( ( xb(nodeb1) - xr(noder1) )^2 ...
                       + ( yb(nodeb1) - yr(noder1) )^2 ) ...
                + sqrt ( ( xb(nodeb2) - xr(noder2) )^2 ...
                       + ( yb(nodeb2) - yr(noder2) )^2 );

          dist2 = sqrt ( ( xb(nodeb1) - xr(noder2) )^2 ...
                       + ( yb(nodeb1) - yr(noder2) )^2 ) ...
                + sqrt ( ( xb(nodeb2) - xr(noder1) )^2 ...
                       + ( yb(nodeb2) - yr(noder1) )^2 );
%
%  If the new arrangement is any shorter, take it, by shuffling the
%  red nodes only, and update the total distance.
%
          if ( dist2 < dist1 )
            node_temp = noder(indx1);
            noder(indx1) = noder(indx2);
            noder(indx2) = node_temp;
            nswap = nswap + 1;
          end

        end

      end

    end
%
%  We've checked all pairs of nodes,
%  print the new total distance, and see if we may
%  continue, or should give up.
%
    total = 0.0;
    for indx1 = 1 : nnode

      nodeb1 = nodeb(indx1);
      noder1 = noder(indx1);

      total = total + sqrt ( ( xb(nodeb1) - xr(noder1) )^2 ...
                           + ( yb(nodeb1) - yr(noder1) )^2 );

    end

    fprintf ( 1, '  On step %d\n', it );
    fprintf ( 1, '  discrepancy = %g\n', total );
    fprintf ( 1, '  Swaps made was %d\n', nswap );

    if ( oldtot - total <= tol * oldtot )

      temp = ( oldtot - total ) / oldtot;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'graph_dist_pairing_greedy(): Warning:\n' );
      fprintf ( 1, '  The relative change in the discrepancy \n' );
      fprintf ( 1, '  was only %g\n', temp );
      fprintf ( 1, '  which is less than the tolerance TOL = %g\n', tol );
      fprintf ( 1, '  Bailing out of the iteration.\n' );
      fprintf ( 1, '\n' );
      return

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'graph_dist_pairing_greedy(): Note:\n' );
  fprintf ( 1, '  The discrepancy has decreased by at least the\n' );
  fprintf ( 1, '  tolerance on every step.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Increasing the number of iterations might\n' );
  fprintf ( 1, '  provide further improvement at this rate.\n' );

  return
end
