(*****************************************************************************
 *           Procedures for generating and counting Laman graphs             *
 *                  (by Christoph Koutschan, June 20 2016)                   *
 *****************************************************************************)

(* We represent a graph by the integer obtained by flattening the upper triangle
 * of its adjacency matrix and interpreting this binary sequence as an integer.
 * Conversely, G2Mat retrieves the adjacency matrix from this integer representation.
 *)
Mat2G[mat_] := FromDigits[Flatten[MapIndexed[Drop[#1, #2[[1]]]&, mat]], 2];
G2Mat[g_Integer, n_Integer] := (# + Transpose[#])&[PadLeft[
  Table[Take[#, {(2n-i)*(i-1)/2+1, (2n-i-1)*i/2}], {i, n}]&[PadLeft[IntegerDigits[g, 2], n*(n-1)/2]], {n, n}]];
G2Mat[g_Integer] := G2Mat[g, Floor[(3 + Sqrt[8*Floor[Log[2, g]] + 1]) / 2]];


(* Input is an undirected graph represented by its symmetric n*n adjacency matrix.
 * Among all graphs that differ from g only by relabelling the vertices, compute a
 * unique representative; we chose this representative to be the lexicographically
 * largest among all equivalent graphs with decreasing vertex neighborhoods.
 * "Lexicographically" in the sense that we interpret the adjacency matrix as a
 * {0,1}-word of length n^2 (not by flattening but by taking blocks of increasing size).
 *)
GraphNormalForm[gr1_List] :=
Module[{gr = gr1, rn = Range[Length[gr1]], nbs, perm, b, p1},

  (* Apply some permutation to the graph s.t. the neighborhoods of its vertices are
     decreasing (by length and lexicographically). By a "neighborhood" of a vertex
     we mean the sorted list of valencies of the neighboring vertices. *)
  nbs = Total /@ gr; nbs = Sort[Pick[nbs, #, 1]]& /@ gr;
  {nbs, perm} = Transpose[Reverse[Sort[Transpose[{nbs, rn}]]]];
  gr = gr[[perm,perm]];

  (* Among all graphs whose vertices are ordered as above, pick the one with
     "lexicographically largest" matrix. *)
  b = Map[Last, SplitBy[Transpose[{nbs, rn}], First], {2}];
  perm = {{}};
  Do[
    perm = Join @@ Outer[Join, perm, Permutations[b[[k]]], 1];
    p1 = FromDigits[Flatten[gr[[#,#]]], 2]& /@ perm;
    perm = Pick[perm, p1, Max[p1]];
  , {k, Length[b]}];
  Return[gr[[#,#]]& @ perm[[1]]];
];


(* Given g, an encoded graph with n-1 vertices, compute all graphs (with n vertices)
 * that can be obtained from g by a single Henneberg move (type I or type II).
 *)
Henneberg[g_Integer] := Henneberg[Floor[(3 + Sqrt[8*Floor[Log[2, g]] + 1]) / 2], g];
Henneberg[n_Integer, g_Integer] :=
Module[{gr, ed},
  gr = PadRight[G2Mat[g, n - 1], {n, n}];
  ed = Select[Position[gr, 1], Less @@ # &];
  Union[Flatten[{

    (* Perform all possible type I Henneberg moves. *)
    Table[Mat2G[GraphNormalForm[ReplacePart[gr, {{v1, n} -> 1, {n, v1} -> 1, {v2, n} -> 1, {n, v2} -> 1}]]],
      {v1, n - 1}, {v2, v1 + 1, n - 1}],

    (* Perform all possible type II Henneberg moves. *)
    Table[
      {v1, v2} = ed[[k]];
      Mat2G[GraphNormalForm[ReplacePart[gr, {{v1, v2} -> 0, {v2, v1} -> 0, {v1, n} -> 1, {n, v1} -> 1,
        {v2, n} -> 1, {n, v2} -> 1, {#, n} -> 1, {n, #} -> 1}]]]& /@ Complement[Range[n - 1], {v1, v2}],
      {k, Length[ed]}]
  }]]];


(* Similar as before, but only with Henneberg type I moves. *)
Henneberg1[g_Integer] := Henneberg1[Floor[(3 + Sqrt[8*Floor[Log[2, g]] + 1]) / 2], g];
Henneberg1[n_Integer, g_] :=
Module[{gr},
  gr = PadRight[G2Mat[g, n - 1], {n, n}];
  Union[Flatten[
    (* Perform all possible type I Henneberg moves. *)
    Table[Mat2G[GraphNormalForm[ReplacePart[gr, {{v1, n} -> 1, {n, v1} -> 1, {v2, n} -> 1, {n, v2} -> 1}]]],
      {v1, n - 1}, {v2, v1 + 1, n - 1}]
  ]]];


(* Compute a list of Laman graphs with n vertices. *)
LamanGraphs[3] = {7};
LamanGraphs[n_ /; n > 3] := LamanGraphs[n] = Union @@ (Henneberg[n, #]& /@ LamanGraphs[n - 1]);

(* Compute a list of H1-Laman graphs (i.e., obtained by only H1 moves) with n vertices. *)
H1LamanGraphs[3] = {7};
H1LamanGraphs[n_ /; n > 3] := H1LamanGraphs[n] = Union @@ (Henneberg1[n, #]& /@ H1LamanGraphs[n - 1]);