/* Edmonds-Random.cpp -- An example program for finding maximum matchings on graphs.
 * James S. Plank

Revision 1.1
December 5, 2017

James S. Plank
Department of Electrical Engineering and Computer Science
University of Tennessee
Knoxville, TN 37996
jplank@utk.edu

Copyright (c) 2017, James S. Plank
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

 - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.

 - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in
   the documentation and/or other materials provided with the
   distribution.

 - Neither the name of the University of Tennessee nor the names of its
   contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY
WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/


/* The first 488 lines here are from Edmonds.cpp.  The rest is a main() that generates
   a graph of a random size, and then prints the graph and its matching. */

#include <string>
#include <vector>
#include <list>
#include <algorithm>
#include <map>
#include <set>
#include <iostream>
#include <cstdio>
#include <cstdlib>
using namespace std;

#define VIT(i, v) for (i = 0; i < v.size(); i++) 
#define IT(it, ds) for (it = ds.begin(); it != ds.end(); it++)
#define FUP(i, n) for (i = 0; i < n; i++)

class Vertex {
  public:
    string name;
    list <class Edge *> edges;
    class Blossom *blossom;
    class Blossom *in_blossom;
    class Blossom *tmp;
    int exposed;
    int marked;
    int d_to_root;
    Vertex *root;
    class Edge *edge_to_root;
    class Edge *matching_edge;
    int in_forest;
    list <Vertex *>::iterator vertices_iterator;
    multimap <double, Vertex *>::iterator unmarked_iterator;

    Vertex(string n);
};

class Edge {
   public:
     string name;
     Vertex *v1;
     Vertex *v2;
     Edge *backedge;
     int marked;
     class Blossom *blossom;
     list <Edge *>::iterator v1_iterator;
     list <Edge *>::iterator matching_iterator;
     list <Edge *>::iterator edge_list_iterator;
};

class Blossom {
   public:
     list <Vertex *> vertices;   /* The order is in the order of the cycle */
     list <Edge *> cycle;
};

typedef list <Vertex *>::iterator LVIT;

class Graph {
  public:
    map <string, Vertex *> v_by_name;
    list <Vertex *> vertices;
    list <Edge *> edges;
    list <Edge *> matching;
    map <string, Vertex *> forest;

    void Add_Edge(string v1, string v2);
    int Find_Augmenting_Path();
    void Add_Edge_To_Matching(Edge *e);  
    Edge *Create_Edge_And_Backedge(Vertex *v1, Vertex *v2);
    void Lift(Vertex *vb);
    void Find_Matching();
};

void Graph::Add_Edge(string vv1, string vv2)
{
  Vertex *v1, *v2;
  map <string, Vertex *>::iterator vit;

  vit = v_by_name.find(vv1);
  if (vit == v_by_name.end()) {
    v1 = new Vertex(vv1);
    vertices.push_front(v1);
    v1->vertices_iterator = vertices.begin();
    v_by_name[vv1] = v1;
  } else {
    v1 = vit->second;
  }

  vit = v_by_name.find(vv2);
  if (vit == v_by_name.end()) {
    v2 = new Vertex(vv2);
    vertices.push_front(v2);
    v2->vertices_iterator = vertices.begin();
    v_by_name[vv2] = v2;
  } else {
    v2 = vit->second;
  }
  Create_Edge_And_Backedge(v1, v2);
}

Edge *Create_Edge(Vertex *v1, Vertex *v2)
{
  Edge *e;

  e = new Edge;
  e->v1 = v1;
  e->v2 = v2;
  e->blossom = NULL;
  e->name = "[" + v1->name + ":" + v2->name + "]";
  v1->edges.push_front(e);
  e->v1_iterator = v1->edges.begin();
  return e;
}

Edge *Graph::Create_Edge_And_Backedge(Vertex *v1, Vertex *v2)
{
  Edge *e1, *e2;

  e1 = Create_Edge(v1, v2);
  e2 = Create_Edge(v2, v1);
  e1->backedge = e2;
  e2->backedge = e1;
  edges.push_front(e1);
  e1->edge_list_iterator = edges.begin();
  edges.push_front(e2);
  e2->edge_list_iterator = edges.begin();
  return e1;
}

Vertex::Vertex(string n)
{
  name = n;
  exposed = 1;
  blossom = NULL;
  in_blossom = NULL;
  tmp = NULL;
}

Vertex *Find_Head_node(Vertex *v1, Vertex*v2)
{
  while(v1->d_to_root > v2->d_to_root) v1 = v1->edge_to_root->v2;
  while(v2->d_to_root > v1->d_to_root) v2 = v2->edge_to_root->v2;
  while (v1 != v2) {
    v1 = v1->edge_to_root->v2;
    v2 = v2->edge_to_root->v2;
  }
  return v1;
}

void Graph::Add_Edge_To_Matching(Edge *e)
{
  Edge *be;

  // cout << "Add_Edge_To_Matching(" << e->name << ")\n";
  matching.push_front(e);
  e->matching_iterator = matching.begin();
  be = e->backedge;
  matching.push_front(be);
  be->matching_iterator = matching.begin();
  e->v1->matching_edge = e;
  be->v1->matching_edge = be;

  if (e->v1->root == e->v1) {
    e->v1->exposed = 0;
  } else {
    matching.erase(e->v1->edge_to_root->matching_iterator);
    matching.erase(e->v1->edge_to_root->backedge->matching_iterator);
    Add_Edge_To_Matching(e->v1->edge_to_root->v2->edge_to_root->backedge);
  }

  if (e->v1->root != e->v2->root) { /* This code only gets executed on the initial call, not the recursive ones.  */
    if (e->v2->root == e->v2) {
      e->v2->exposed = 0;
    } else {
      matching.erase(e->v2->edge_to_root->matching_iterator);
      matching.erase(e->v2->edge_to_root->backedge->matching_iterator);
      Add_Edge_To_Matching(e->v2->edge_to_root->v2->edge_to_root->backedge);
    }
  }

  /* When we return, the augmenting path is done, so don't worry about 
     roots or edges to roots right now.  Fix that in a later iteration. */
}

int Graph::Find_Augmenting_Path()
{
  int retval;
  Vertex *v1, *v2, *vhead, *vb, *v;
  Blossom *b;
  Edge *e1, *e, *be1;
  list <Edge *> tmplist;
  list <Vertex *>::iterator vit;
  list <Edge *>::iterator mit;
  list <Edge *>::iterator eit;
  list <Edge *>::iterator v1eit;
  multimap <double, Vertex *> unmarked_v;
  multimap <double, Vertex *>::iterator uvit;

  // cout << "Finding Augmenting Path\n";
  // Print_Matching();

  forest.clear();

  for (vit = vertices.begin(); vit != vertices.end(); vit++) {
    v1 = *vit;
    v1->tmp = NULL;
    v1->in_blossom = NULL;
    if (v1->exposed) {
      v1->in_forest = 1;
      v1->d_to_root = 0;
      v1->root = v1;
      v1->edge_to_root = NULL;
      v1->marked = 0;
      v1->unmarked_iterator = unmarked_v.insert(make_pair(drand48(), v1));
      forest[v1->name] = v1;
    } else {
      v1->in_forest = 0;
    }
  }


  for (mit = edges.begin(); mit != edges.end(); mit++) {
    e1 = *mit;
    e1->marked = 0;
    e1->blossom = NULL;
  }

  for (mit = matching.begin(); mit != matching.end(); mit++) {
    e1 = *mit;
    e1->marked = 1;
  }

  /* 1. Find an unmarked vertex v in F s.t. distance v->root is even */

  while (1) {
    if (unmarked_v.empty()) return 0;
    uvit = unmarked_v.begin();
    v1 = uvit->second;
 
    // cout << "Considering unmarked vertex " << v1->name << " d_to_root: " << v1->d_to_root << endl;
    if (v1->d_to_root%2 == 0) {

    /* Traverse v1's edges to find unmarked edges */
    
      for (v1eit = v1->edges.begin(); v1eit != v1->edges.end(); v1eit++) {
        e1 = *v1eit;
        if (e1->marked == 0) {
          // cout << "Considering unmarked edge " << e1->name << endl;

          /* Case 3: Has an important by-product -- this will add a vertex to the unmarked vertices pool.  */
          if (!e1->v2->in_forest) {
            // cout << "Case 3 -- Adding " << e1->v2->name << "(" << e1->v2->exposed << ") & " << e1->v2->matching_edge->v2->name << " to " << e1->v1->name << "'s tree.\n";
   

            e1->v2->root = v1->root;
            e1->v2->d_to_root = v1->d_to_root + 1;
            e1->v2->edge_to_root = e1->backedge;
            e1->v2->in_forest = 1;
            
            v2 = e1->v2->matching_edge->v2;
            v2->root = v1->root;
            v2->d_to_root = v1->d_to_root + 2;
            v2->edge_to_root = v2->matching_edge;
            v2->in_forest = 1;
            v2->unmarked_iterator = unmarked_v.insert(make_pair(drand48(), v2));

            e1->marked = 1;
            e1->backedge->marked = 1;
      
          /* Case 4 */
          } else if (e1->v2->d_to_root % 2 == 1) {
            // cout << "Case 4 -- Mark edge " << e1->name << endl;
            e1->marked = 1;
            e1->backedge->marked = 1;

          /* Case 5: We have found an augmenting path from v1->root to v2->root.
             We need to remove edges along that path that are currently matching, and 
             remove the others.  This will be done with a recursive procedure called
             add_edge_to_matching. */
          } else if (v1->root != e1->v2->root) { 
            // cout << "Case 5 -- Adding Edge " << e1->name << " to Matching\n";
            Add_Edge_To_Matching(e1);
            return 1;
          
          } else {
            // cout << "Case 6 -- Discovered a blossom\n";
            v2 = e1->v2;
            vb = new Vertex("*-"+v1->name);
            vertices.push_front(vb);
            vb->vertices_iterator = vertices.begin();
            b = new Blossom;
            vb->blossom = b;
            vhead = Find_Head_node(v1, v2); 

            /* Create the cycle on b's edge list */
            
            b->cycle.push_back(e1->backedge);
            for (v = v1; v != vhead; v = v->edge_to_root->v2) {
              b->cycle.push_back(v->edge_to_root);
            } 
            for (v = v2; v != vhead; v = v->edge_to_root->v2) {
              tmplist.push_front(v->edge_to_root->backedge);
            } 
            for (eit = tmplist.begin(); eit != tmplist.end(); eit++) {
              b->cycle.push_back(*eit);
            }

            /* Create the vertex list in the same order as the cycle */
            for (eit = b->cycle.begin(); eit != b->cycle.end(); eit++) {
              e = *eit;
              b->vertices.push_back(e->v1);
            }
            tmplist.clear();

            /* Get rid of all matching edges for the blossom */

            for (vit = b->vertices.begin(); vit != b->vertices.end(); vit++) {
              v = *vit;
              v->in_blossom = b;
              if (!v->exposed) {
                matching.erase(v->matching_edge->matching_iterator);
                v->exposed = 1;
              }
            } 

            /* Now, traverse the vertices, and delete each edge from g.edges.
               If an edge goes from inside the blossom to outside the blossom,
               then if you should, create the new edges to/from vb.  Manage
               matching edges properly. */

            for (vit = b->vertices.begin(); vit != b->vertices.end(); vit++) {
              v = *vit;
              vertices.erase(v->vertices_iterator);
              for (eit = v->edges.begin(); eit != v->edges.end(); eit++) {
                e = *eit;
                edges.erase(e->edge_list_iterator);
                e->blossom = b;
                if (e->v2->in_blossom != b) {
                  if (e->v2->tmp != b) {
                    be1 = Create_Edge_And_Backedge(vb, e->v2);
                    e->v2->tmp = b;
                  }
                  if (!e->v2->exposed && e->v2->matching_edge == e->backedge) {
                    be1 = *(e->v2->edges.begin());
                    be1 = be1->backedge;
                    matching.erase(e->backedge->matching_iterator);
                    matching.push_front(be1);
                    be1->matching_iterator = matching.begin();
                    matching.push_front(be1->backedge);
                    be1->backedge->matching_iterator = matching.begin();
                    vb->matching_edge = be1;
                    e->v2->matching_edge = be1->backedge;
                    vb->exposed = 0;
                  }
                  e->v2->edges.erase(e->backedge->v1_iterator);
                  edges.erase(e->backedge->edge_list_iterator);
                }
              }
            } 
            retval = Find_Augmenting_Path();
            Lift(vb);
            return retval;
          }
        }
      }
    }

    unmarked_v.erase(uvit);
    v1->marked = 1;
  }
}

void Graph::Lift(Vertex *vb)
{
  list <Vertex *>::iterator vit;
  list <Edge *>::iterator eit;
  Blossom *b;
  Vertex *connecting_vertex;
  Vertex *v;
  Edge *e;
  int do_match;

  // cout << "Lifting " << vb->name << endl;
  // Print_Graph();

  connecting_vertex = NULL;

  b = vb->blossom;
  if (b == NULL) return;

  /* Restore all vertices and edges to their rightful places. */

  for (vit = b->vertices.begin(); vit != b->vertices.end(); vit++) {
    v = *vit;
    vertices.push_front(v);
    v->vertices_iterator = vertices.begin();
  }

  for (vit = b->vertices.begin(); vit != b->vertices.end(); vit++) {
    v = *vit;
    for (eit = v->edges.begin(); eit != v->edges.end(); eit++) {
      e = *eit;
      edges.push_front(e);
      e->edge_list_iterator = edges.begin();
      e->blossom = NULL;
      if (e->v2->in_blossom != b) {
        e->v2->edges.push_front(e->backedge);
        e->backedge->v1_iterator = e->v2->edges.begin();
        edges.push_front(e->backedge);
        e->backedge->edge_list_iterator = edges.begin();
        if (!vb->exposed && e->v2 == vb->matching_edge->v2) {
          matching.erase(vb->matching_edge->matching_iterator);
          matching.erase(e->v2->matching_edge->matching_iterator);
          vb->exposed = 1;
          matching.push_front(e);
          e->matching_iterator = matching.begin();
          e->v1->matching_edge = e;
          e->v1->exposed = 0;
          connecting_vertex = e->v1;
          matching.push_front(e->backedge);
          e->backedge->matching_iterator = matching.begin();
          e->v2->matching_edge = e->backedge;
        } 
      }
    }
  }

  for (vit = b->vertices.begin(); vit != b->vertices.end(); vit++) {
    v = *vit;
    v->in_blossom = NULL;
    v->tmp = NULL;
  }

  vertices.erase(vb->vertices_iterator);
  for (eit = vb->edges.begin(); eit != vb->edges.end(); eit++) {
    e = *eit;
    e->v2->edges.erase(e->backedge->v1_iterator);
    edges.erase(e->edge_list_iterator);
    edges.erase(e->backedge->edge_list_iterator);
  }

  if (connecting_vertex == NULL) {
    eit = b->cycle.begin();
    connecting_vertex = (*eit)->v1;
  } else {
    for (eit = b->cycle.begin(); (*eit)->v1 != connecting_vertex; eit++) ;
  }

  e = *eit;
  do_match = 0;
  while (e->v2 != connecting_vertex) {
    if (do_match) {
      matching.push_front(e);
      e->matching_iterator = matching.begin();
      e->v1->exposed = 0;
      e->v1->matching_edge = e;

      matching.push_front(e->backedge);
      e->backedge->matching_iterator = matching.begin();
      e->v2->exposed = 0;
      e->v2->matching_edge = e->backedge;
    }
    do_match = 1-do_match;
    eit++;
    if (eit == b->cycle.end()) eit = b->cycle.begin();
    e = *eit;
  }

  delete b;

  for (eit = vb->edges.begin(); eit != vb->edges.end(); eit++) {
    e = *eit;
    delete e->backedge;
    delete e;
  }
  delete vb;
  
}

void Graph::Find_Matching()
{
  while (Find_Augmenting_Path()) {
  }
}

void usage(const string s)
{
  fprintf(stderr, "usage: edmonds nodes edges seed(-1 to read) print(y|n|g)\n");
  if (s != "") fprintf(stderr, "%s\n", s.c_str());
  exit(1);
}

int main(int argc, char **argv)
{
  int nodes, edges, print, i, j;
  long seed;
  Graph *g;
  string p, s;
  vector <double> X, Y;
  vector <string> names;
  char buf[100];
  multimap <double, int> D;
  multimap <double, int>::iterator dit;
  map <int, int> E;                     // Edges.  second is whether it's in the matching.
  map <int, int>::iterator eit;
  double d, x, y;
  list <Edge *>::iterator mit;
  Edge *e;

  /* Parse and error check the command line. */

  if (argc != 5) usage("");
  nodes = atoi(argv[1]);
  edges = atoi(argv[2]);
  if (sscanf(argv[3], "%ld", &seed) == 0) usage("Seed must be a number");
  if (seed != -1 && nodes <= 0) usage("Nodes must be > 0");
  if (seed != -1 && edges <= 0) usage("Edges must be > 0");

  p = argv[4];
  if (p != "y" && p != "n" && p != "g") usage("print must be y, n or g");
  print = p[0];
  
  if (seed == -1) {
    nodes = 0;
    while (cin >> x >> y) {
      X.push_back(x);
      Y.push_back(y);
      sprintf(buf, "%d", nodes);
      names.push_back(buf);
      nodes++;
    }
    cin.clear();
    cin >> s;
    while (cin >> i) E[i] = 0;
    edges = E.size();
  } else {

    /* Create nodes with random X and Y values from 0 to 10. */
    srand48(seed);
  
    for (i = 0; i < nodes; i++) {
      X.push_back(drand48()*10.0);
      Y.push_back(drand48()*10.0);
      sprintf(buf, "%d", i);
      names.push_back(buf);
    }
  
    /* Now, create D to be a multimap of potential edges, ordered by distance (squared). */
  
    for (i = 0; i < nodes; i++) {
      for (j = 0; j < i; j++) {
        d = (X[i]-X[j])*(X[i]-X[j])+(Y[i]-Y[j])*(Y[i]-Y[j]);
        D.insert(make_pair(-d, i*nodes+j));
        if (D.size() > edges) D.erase(D.begin());
      }
    }
    for (dit = D.begin(); dit != D.end(); dit++) E[dit->second] = 0;
  }

  if (print == 'g') {
    for (i = 0; i < X.size(); i++) printf("%.4lf %.4lf\n", X[i], Y[i]);
    printf("Edges\n");
    for (eit = E.begin(); eit != E.end(); eit++) printf("  %d\n", eit->first);
    exit(0);
  }


  g = new Graph;

  for (eit = E.begin(); eit != E.end(); eit++) {
    i = eit->first/nodes;
    j = eit->first%nodes;
    g->Add_Edge(names[i], names[j]);
  }

  g->Find_Matching();
  for (mit = g->matching.begin(); mit != g->matching.end(); mit++) {
    e = *mit;
    i = atoi(e->v1->name.c_str());
    j = atoi(e->v2->name.c_str());
    if (j < i) {
      if (E.find(i*nodes+j) == E.end()) {
      } else {
        E[i*nodes+j] = 1;
      }
    }
  }
  
  if (print == 'y') {
    printf("newgraph\n");
    printf("xaxis nodraw min 0 max 10 size 5\n");
    printf("yaxis nodraw min 0 max 10 size 5\n");
    for (eit = E.begin(); eit != E.end(); eit++) {
      i = eit->first/nodes;
      j = eit->first%nodes;
      printf("newline linethickness 2 color %d 0 0 pts %.4lf %.4lf  %.4lf %.4lf\n", 
              eit->second, X[i], Y[i], X[j], Y[j]);
    }
    printf("newcurve marktype circle pts\n");
    for (i = 0; i < X.size(); i++) printf("  %.4lf %.4lf\n", X[i], Y[i]);
  }
  return 0;
}