/*******************************************************************************
+
+  LEDA 3.5
+
+  _g_generate.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
#include <LEDA/graph.h>
#include <LEDA/ugraph.h>
#include <LEDA/vector.h>
#include <LEDA/matrix.h>
#include <LEDA/array2.h>

#include <ctype.h>
#include <math.h>

//------------------------------------------------------------------------------
// some graph generators
//
// S. Naeher  (1995-1996)
//------------------------------------------------------------------------------

void complete_graph(graph& G, int n, bool directed)
{ 
  G.clear();

  node* V = new node[n];

  for(int i=0;i<n;i++) V[i] = G.new_node();

  if (directed)
    { //memory_allocate_block(sizeof(node_struct),n);
      for(int i=0;i<n;i++) 
        for(int j=0;j<n;j++) G.new_edge(V[i],V[j]);
     }
  else
    { //memory_allocate_block(sizeof(edge_struct),n*n/2);
      for(int i=0; i<n; i++) 
        for(int j=i+1; j<n; j++) G.new_edge(V[i],V[j]);
     }

  delete[] V;
}



void grid_graph(graph& G, int n) 
{ node_array<double> xcoord; 
  node_array<double> ycoord;
  grid_graph(G,xcoord,ycoord,n);
 }

void grid_graph(graph& G, node_array<double>& xcoord, 
                          node_array<double>& ycoord, int n) 
{
  array2<node>  A(n,n);
  node v;
  int N = n*n;
  int x;
  int y;

  double d = 1.0/(n+1);

  G.clear();

  xcoord.init(G,N,0);
  ycoord.init(G,N,0);

  for(y=0; y<n; y++)
    for(x=0; x<n; x++)
      { A(x,y) = v = G.new_node();
        xcoord[v] = (x+1)*d;
        ycoord[v] = (y+1)*d;
       }

  for(x=0; x<n; x++)
    for(y=0; y<n; y++)
       { if (x < n-1) G.new_edge(A(x,y),A(x+1,y));
         if (y < n-1) G.new_edge(A(x,y),A(x,y+1));
        }
}


void complete_bigraph(graph& G, int n1, int n2, list<node>& A, list<node>& B)
{ 
  G.clear();

  while (n1--)  A.append(G.new_node());
  while (n2--)  B.append(G.new_node());

  list_item a,b;
  forall_items(a,A) 
    forall_items(b,B) 
      G.new_edge(A[a],B[b]);
    
 }

void user_graph(graph& G)
{ int  n = read_int("|V| = ");
  int  i,j;

  node* V = new node[n];
  for(j=0; j<n; j++) V[j] = G.new_node();

  for(j=0; j<n; j++) 
  { list<int> il;
    int ok = false;
    while (!ok)
    { ok = true;
      cout << "edges from [" << j << "] to: ";
      il.read();
      forall(i,il) 
        if (i < 0 || i >= n) 
        { ok=false;
          cout << "illegal node " << i << "\n";
         }
     }
    forall(i,il) G.new_edge(V[j],V[i]);
  }
  G.print();
  if (Yes("save graph ? ")) G.write(read_string("file: "));

  delete[] V;
}



void test_graph(graph& G)
{ 
  G.clear();
  char c;

  do c = read_char("graph: f(ile) r(andom) c(omplete) p(lanar) u(ser): ");
  while (c!='f' && c!='r' && c!='c' && c!='p'&& c!='u');   

  switch (c) {

   case 'f' : { G.read(read_string("file: "));
                break;
               }

   case 'u' : { user_graph(G);
                break;
               }

   case 'c' : { complete_graph(G,read_int("|V| = "));
                break;
               }

   case 'r' : { int n = read_int("|V| = ");
                int m = read_int("|E| = ");
                random_graph(G,n,m);
                break;
               }

   case 'p' : { random_planar_graph(G,read_int("|V| = "));
                break;
               }
   }//switch

}

void test_ugraph(ugraph& G)
{ 
  G.clear();
  char c;

  do c = read_char("graph: f(ile) r(andom) c(omplete) p(lanar) u(ser): ");
  while (c!='f' && c!='r' && c!='c' && c!='p'&& c!='u');   

  int  i;
  node v;
  
  switch (c) {

  case 'f' : { G.read(read_string("file: "));
               break;
              }

   case 'u' : { int  n = read_int("|V| = ");
                int  j = 0;
                node* V = new node[n];
                for(i=0; i<n; i++) V[i] = G.new_node();
                forall_nodes(v,G)
                  { list<int> il;
                    cout << "edges from " << j++ << " to: ";
                    il.read();
                    forall(i,il) 
                      if (i >= 0 && i < n) G.new_edge(v,V[i]);
                      else cerr << "illegal node " << i << " (ignored)\n";
                   }
                G.print();
                if (Yes("save graph ? ")) G.write(read_string("file: "));
                delete[] V;
                break;
               }

   case 'c' : { int n = read_int("|V| = ");
                complete_graph(G,n);
                break;
               }

   case 'r' : { int n = read_int("|V| = ");
                int m = read_int("|E| = ");
                random_graph(G,n,m);
                break;
               }

   }//switch

}





void test_bigraph(graph& G, list<node>& A, list<node>& B)
{
  int a,b;
  int n1 = 0;
  int n2 = 0;
  char c;

  do c = read_char("bipartite graph: f(ile) r(andom) c(omplete) u(ser): ");
  while (c!='f' && c!='r' && c!='c' && c!='u');   

  A.clear();
  B.clear();
  G.clear();

  if (c!='f') 
   { n1 = read_int("|A| = ");
     n2 = read_int("|B| = ");
    }
  
  
  switch (c) {

  case 'f' : { G.read(read_string("file: "));
               node v;
               forall_nodes(v,G) 
               if (G.outdeg(v) > 0) A.append(v);
               else B.append(v);

               break;
              }

   case 'u' : { node* AV = new node[n1+1];
                node* BV = new node[n2+1];

                for(a=1; a<=n1; a++)  A.append(AV[a] = G.new_node());
                for(b=1; b<=n2; b++)  B.append(BV[b] = G.new_node());

                for(a=1; a<=n1; a++)
                { list<int> il;
                  cout << "edges from " << a << " to: ";
                  il.read();
                  forall(b,il) 
                    if (b<=n2) G.new_edge(AV[a],BV[b]);
                    else break;
                  if (b>n2) break;
                 }
                delete[] AV;
                delete[] BV;
                break;
               }

   case 'c' : complete_bigraph(G,n1,n2,A,B);
              break;

   case 'r' : { int m = read_int("|E| = ");
                random_bigraph(G,n1,n2,m,A,B);
                break;
               }

       } // switch

}



void cmdline_graph(graph& G, int argc, char** argv)
{ 
  // construct graph from cmdline arguments

  if (argc == 1)           // no arguments 
     { test_graph(G);
       return;
      }
  else 
     if (argc == 2)       // one argument
        { if (isdigit(argv[1][0])) 
             { cout << "complete graph |V| = " << argv[1];
               newline;
               newline;
               complete_graph(G,atoi(argv[1]));
              }
          else 
             { cout << "reading graph from file " << argv[1];
               newline;
               newline;
               G.read(argv[1]);
              }
          return;
         }
     else
        if (argc == 3 && isdigit(argv[1][0]) && isdigit(argv[1][0])) 
           { cout << "random graph |V| = " << argv[1] << "  |E| = " << argv[2];
             newline;
             newline;
             random_graph(G,atoi(argv[1]),atoi(argv[2]));
             return;
            }

  error_handler(1,"cmdline_graph: illegal arguments");
}




//------------------------------------------------------------------------------
// triangulated planar graph
//------------------------------------------------------------------------------


struct triang_point {

double x;
double y;
node   v;

LEDA_MEMORY(triang_point)

triang_point(double a=0, double b = 0) { x = a; y = b; v = nil; }
triang_point(const triang_point& p)    { x = p.x; y = p.y; v = p.v; }
~triang_point() {};

friend bool right_turn(const triang_point& a, const triang_point& b, const triang_point& c)
{ return (a.y-b.y)*(a.x-c.x)+(b.x-a.x)*(a.y-c.y) > 0; }

friend bool left_turn(const triang_point& a, const triang_point& b, const triang_point& c)
{ return (a.y-b.y)*(a.x-c.x)+(b.x-a.x)*(a.y-c.y) < 0; }


friend bool operator==(const triang_point& a, const triang_point& b)
{ return a.x == b.x && a.y == b.y; }

friend ostream& operator<<(ostream& out, const triang_point& p)
{ return out << p.x << " " << p.y; }

friend istream& operator>>(istream& in, triang_point& p)
{ return in >> p.x >> p.y; }

friend int compare(const triang_point& p, const triang_point& q)
{ int c = compare(p.x,q.x);
  if (c==0) c = compare(p.y,q.y);
  return c;
 }

};



void triangulated_planar_graph(graph& G, list<node>& outer_face,
                                         node_array<double>& xcoord,
                                         node_array<double>& ycoord, int n)
{ 
  G.clear();

  list<triang_point> L;

  while(n--)
  { double x = rand_int(0,1000000)/1000000.0; 
    double y = rand_int(0,1000000)/1000000.0; 
    L.append(triang_point(x,y));
   }

  L.sort();  // sort triang_points lexicographically

  list<triang_point> CH;
  list_item last;
  triang_point p,q;


  // eliminate multiple triang_points

  list_item it;
  forall_items(it,L)
  { list_item it1 = L.succ(it);
    while (it1 != nil && L[it1] == L[it])
    { L.del(it1);
      it1 = L.succ(it);
     }
   }

  n = L.length();

  xcoord.init(G,n,0);
  ycoord.init(G,n,0);

  forall_items(it,L)
  { node v = G.new_node();
    xcoord[v] = L[it].x;
    ycoord[v] = L[it].y;
    L[it].v = v;
   }

  // initialize convex hull with first two points

  p = L.pop();
  CH.append(p);

  while (L.head() == p) L.pop();

  q = L.pop();
  last = CH.append(q);

  G.new_edge(p.v,q.v);


  // scan remaining points

  forall(p,L)
  {


    node v = p.v;

    G.new_edge(v,CH[last].v);

    // compute upper tangent (p,up)

    list_item up = last;
    list_item it = CH.cyclic_succ(up);

    while (left_turn(CH[it],CH[up],p))
    { up = it;
      it = CH.cyclic_succ(up);
      G.new_edge(v,CH[up].v);
     }


    // compute lower tangent (p,low)

    list_item low = last;
    it = CH.cyclic_pred(low);

    while (right_turn(CH[it],CH[low],p))
    { low = it;
      it = CH.cyclic_pred(low);
      G.new_edge(v,CH[low].v);
     }


    // remove all points between up and low

    if (up != low)
    { it = CH.cyclic_succ(low);

      while (it != up)
      { CH.del(it);
        it = CH.cyclic_succ(low);
       }
     }

    // insert new point

    last = CH.insert(p,low);

   }

  outer_face.clear();
  forall(p,CH) outer_face.append(p.v);

}


void triangulated_planar_graph(graph& G, int m)
{ node_array<double> xcoord;
  node_array<double> ycoord;
  list<node> L;
  triangulated_planar_graph(G,L,xcoord,ycoord,m);
 }


static bool tutte_embed(const graph& G, const node_array<bool>& fixed,
                        node_array<double>& xpos, node_array<double>& ypos)
{ node v,w;
  edge e;

  list<node> other_nodes;
  forall_nodes(v,G) 
      if(!fixed[v]) other_nodes.append(v);

  node_array<int> ind(G);          // position of v in other_nodes and A

  int i = 0;
  forall(v,other_nodes) ind[v] = i++;

  int n = other_nodes.size();   // #other nodes
  vector coord(n);              // coordinates (first x then y)
  vector rhs(n);                // right hand side
  matrix A(n,n);                // equations

  // initialize non-zero entries in matrix A
  forall(v,other_nodes)
  {
    double one_over_d = 1.0/double(G.degree(v));
    forall_inout_edges(e,v)
    {
      // get second node of e
      w = (v == source(e)) ? target(e) : source(e);
      if(!fixed[w]) A(ind[v],ind[w]) = one_over_d;
    }
    A(ind[v],ind[v]) = -1;
  }

  if(!A.det()) return false;

  // compute right hand side for x coordinates
  forall(v,other_nodes)
  { rhs[ind[v]] = 0;
    double one_over_d = 1.0/double(G.degree(v));
    forall_inout_edges(e,v)
    { // get second node of e
      w = (v == source(e)) ? target(e) : source(e);
      if(fixed[w]) rhs[ind[v]] -= (one_over_d*xpos[w]);
    }
  }

  // compute x coordinates
  coord = A.solve(rhs);
  forall(v,other_nodes) xpos[v] = coord[ind[v]];

  // compute right hand side for y coordinates
  forall(v,other_nodes)
  { rhs[ind[v]] = 0;
    double one_over_d = 1.0/double(G.degree(v));
    forall_inout_edges(e,v)
    { // get second node of e
      w = (v == source(e)) ? target(e) : source(e);
      if(fixed[w]) rhs[ind[v]] -= (one_over_d*ypos[w]);
    }
  }

  // compute y coordinates
  coord = A.solve(rhs);
  forall(v,other_nodes) ypos[v] = coord[ind[v]];

  return true;
}


void triangulated_planar_graph(graph& G, node_array<double>& xcoord,
                                         node_array<double>& ycoord, int n)
{ list<node> L;
  triangulated_planar_graph(G,L,xcoord,ycoord,n);

  if (n > 128) return;

  node_array<bool> fixed(G,false);

  double step  = 6.2832/L.length();
  double alpha = 0;
  node v;
  forall(v,L)
  { xcoord[v] = cos(alpha);
    ycoord[v] = sin(alpha);
    alpha+=step;
    fixed[v] = true;
  }
 
  tutte_embed(G,fixed,xcoord,ycoord);
 }


