Triconnected_Planar.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6e2aa13f",
   "metadata": {},
   "source": [
    "Computing the resistance of a triconnected planar graph.  We use the equivalence of resistance networks to\n",
    "tilings by rectangles observed by @Lisanne.  If we add an edge from the top vertex to the bottom vertex (a \"battery\" edge) the resulting graph is a biconnected planar graph.  Conversely, corresponding to every biconnected\n",
    "planar graph is a rectangle dissection.  So it suffices to run through the planar biconnected graphs.\n",
    "\n",
    "We will compute the rational functions giving the aspect ratio of the rectangle tiled by the smaller rectangles of fixed aspect ratio.\n",
    "\n",
    "If the graph is just biconnected, then removing two vertices will separate a piece that has a top vertex and bottom vertex and no other connection to the larger graph.  We can replace this by a single resistor having a specified rational function resistance.  If we compute the possible aspect ratios for rectangles tiled by fewer than \\(n\\) vertices, then we can plug in the precomputed resistances.  So it suffices to consider only 3-connected planar graphs.\n",
    "\n",
    "Furthermore, as if we take the dual of a graph and invert the resistance of each edge, this is equivalent to rotating the rectangle by \\(90^{\\circ}\\), so it suffices to only consider one of a pair of dual 3-connected planar graphs.\n",
    "\n",
    "Once we have computed the rational functions representing the aspect ratio of the tiled rectangle, we then compute the real positive roots of numerator-denominator, the values where the tiled rectangle is a square.  Finally, we check that none of the edges end up having zero current."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c12d8e28",
   "metadata": {},
   "outputs": [],
   "source": [
    "def equivalent_resistance(graph):\n",
    "    \"\"\"\n",
    "    Compute the rational function describing the resistance of the graph in terms of the resistances of \n",
    "    the edges of the graph.\n",
    "    \n",
    "    INPUT:\n",
    "    \n",
    "    - `graph` The underlying graph.  This graph is assumed to be at least 3-connected and a planar graph with an \n",
    "         embedding.  The graph will have one edge labelled \"b\" (battery) indicating the top and bottom vertex, the\n",
    "         rest of the edges labelled \"r\" (resistor)\n",
    "    \"\"\"\n",
    "    n = graph.num_edges() - 1\n",
    "    \n",
    "    # set up the underlying field of rational functions\n",
    "    R = FractionField(PolynomialRing(QQ, 'x', n))\n",
    "    \n",
    "    #set up a dictionary of edges and source and sink\n",
    "    edge_lookup = dict()\n",
    "    pos = 0\n",
    "    source = 0\n",
    "    sink = 0\n",
    "    for a, b, label in graph.edges(sort=True):\n",
    "        if label != 'b':\n",
    "            edge_lookup[(a, b)] = pos\n",
    "            edge_lookup[(b, a)] = pos\n",
    "            pos += 1\n",
    "        else:\n",
    "            if a < b:\n",
    "                source = a\n",
    "                sink = b\n",
    "            else:\n",
    "                source = b\n",
    "                sink = a\n",
    "    \n",
    "    vects = []\n",
    "    \n",
    "    # flows out of the source\n",
    "    source_current = [0 for _ in range(n)]\n",
    "    for v in graph.neighbor_iterator(source):\n",
    "        if v != sink:\n",
    "            e = edge_lookup[(source, v)]\n",
    "            source_current[e] = 1 if v > source else -1\n",
    "    vects.append(source_current)\n",
    "\n",
    "    # flows through non source sink vertices\n",
    "    for v in graph.vertices(sort=True):\n",
    "        if v != source and v != sink:\n",
    "            temp_current = [0 for _ in range(n)]\n",
    "            for v1 in graph.neighbor_iterator(v):\n",
    "                e = edge_lookup[(v, v1)]\n",
    "                temp_current[e] = 1 if v1 > v else -1\n",
    "            vects.append(temp_current)\n",
    "    \n",
    "    source_sink_paths = []\n",
    "    # flows around the faces\n",
    "    for circuit in graph.faces():\n",
    "        temp_voltage = [0 for _ in range(n)]\n",
    "        if (source, sink) in circuit:\n",
    "            for (v1,v2) in circuit:\n",
    "                if (v1, v2) != (source, sink):\n",
    "                    e = edge_lookup[(v1, v2)]\n",
    "                    temp_voltage[e] = -R.gen(e) if v1 < v2 else R.gen(e)\n",
    "            source_sink_paths.append(temp_voltage)\n",
    "        elif (sink, source) in circuit:\n",
    "            for (v1,v2) in circuit:\n",
    "                if (v1, v2) != (sink, source):\n",
    "                    e = edge_lookup[(v1, v2)]\n",
    "                    temp_voltage[e] = R.gen(e) if v1 < v2 else -R.gen(e)\n",
    "            source_sink_paths.append(temp_voltage)\n",
    "        else:\n",
    "            for (v1, v2) in circuit:\n",
    "                e = edge_lookup[(v1, v2)]\n",
    "                temp_voltage[e] = R.gen(e) if v1 < v2 else -R.gen(e)\n",
    "            vects.append(temp_voltage)\n",
    "            \n",
    "    # find the current flows through each resistor if we apply unit current from source to sink.\n",
    "    m = matrix(vects)\n",
    "    external_current = vector([1 if x == 0 else 0 for x in range(n)]) \n",
    "    internal_current = m.inverse()*external_current\n",
    "\n",
    "    # the voltage drop is the sum of I*R over a path from sourc to sink\n",
    "    external_voltage = internal_current.dot_product(vector(source_sink_paths[0]))\n",
    "    ev_sign = external_voltage.numerator().coefficients()[0].sign()\n",
    "    ev = ev_sign * external_voltage.numerator() / (ev_sign * external_voltage.denominator())\n",
    "    return (ev, internal_current)       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ed7655e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n",
      "number of base graphs =  1\n",
      "5\n",
      "number of base graphs =  3\n",
      "6\n",
      "number of base graphs =  9\n",
      "7\n",
      "number of base graphs =  16\n",
      "8\n",
      "number of base graphs =  16\n",
      "9\n",
      "number of base graphs =  16\n",
      "Counter({12: 52, 11: 11, 10: 8, 8: 2, 9: 2, 6: 1})\n"
     ]
    }
   ],
   "source": [
    "# compute the triconnected planar graphs with 12 or fewer edges\n",
    "# we only keep one of a graph and its dual\n",
    "\n",
    "from collections import Counter\n",
    "from collections import defaultdict\n",
    "\n",
    "base_graphs=[]\n",
    "labelled_graphs=[]\n",
    "graph_determinants=defaultdict(list)\n",
    "\n",
    "for n in range(4,10):\n",
    "    print(n)\n",
    "    # generate graphs with n vertices, up to 12 edges that are biconnected with minimal degree 3\n",
    "    gen = graphs.nauty_geng(\"{0} 0:12 -C -d3\".format(n))\n",
    "    try:\n",
    "        for g in gen:\n",
    "            # check for triconnected and planar\n",
    "            if g.vertex_connectivity() > 2 and g.is_planar(set_embedding=True):\n",
    "                gd = g.planar_dual()\n",
    "                # if we have already stored the dual, ignore this\n",
    "                if not any(gd.is_isomorphic(g1) for g1 in base_graphs):\n",
    "                    base_graphs.append(g)\n",
    "    except ValueError:\n",
    "        pass\n",
    "        \n",
    "    print(\"number of base graphs = \",len(base_graphs))\n",
    "    \n",
    "# Compute the polar nets that can be obtained from the above graphs up to isomorphism.\n",
    "# Label one edge to be the battery, and check if it is not isomorphic to a previous arrangement\n",
    "for g in base_graphs:\n",
    "    h = g.copy()\n",
    "    for u,v,l in g.edges(sort=False):\n",
    "        h.set_edge_label(u, v, \"r\")\n",
    "    g_edges = g.edges(sort=False)\n",
    "    u,v,l = g_edges[0]\n",
    "    g2 = h.copy()\n",
    "    g2.set_edge_label(u, v, \"b\")\n",
    "    temp_list = [(u,v,g2)]\n",
    "    for u,v,l in g_edges[1:]:\n",
    "        g2 = h.copy()\n",
    "        g2.set_edge_label(u,v,\"b\")\n",
    "        if not any(g2.is_isomorphic(g3,edge_labels=True) for _,_,g3 in temp_list):\n",
    "            temp_list.append((u,v,g2))\n",
    "    labelled_graphs.extend(temp_list)\n",
    "counts = Counter(g.size() for _,_,g in labelled_graphs)\n",
    "print(counts)\n",
    "for _,_,g in labelled_graphs:\n",
    "    if g.size() <= 13:\n",
    "        graph_determinants[g.num_edges()-1].append(equivalent_resistance(g))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "951748c0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1   []   1\n",
      "2   3\n",
      "2   []   1\n",
      "3   10\n",
      "3   []   3\n",
      "4   38\n",
      "4   []   11\n",
      "5   164\n",
      "5   [1]   51\n",
      "6   787\n",
      "6   []   245\n",
      "7   4258\n",
      "7   []   1372\n",
      "8   25767\n",
      "8   []   8522\n",
      "9   172504\n",
      "9   []   58347\n",
      "10   1258949\n",
      "10   []   433926\n",
      "[1, 3, 10, 38, 164, 787, 4258, 25767, 172504, 1258949]\n",
      "[1, 1, 3, 11, 51, 245, 1372, 8522, 58347, 433926]\n"
     ]
    }
   ],
   "source": [
    "from itertools import combinations_with_replacement as cwr, product as itprod\n",
    "\n",
    "R.<w>=PolynomialRing(ZZ)\n",
    "K=FractionField(R)\n",
    "\n",
    "# A list where resistances[n] gives a rational function describing the aspect ratio of the enclosing rectangle\n",
    "# and a description of the topology of the subrectangles (see construct_biconnected for a description)\n",
    "resistances = [[],[(K(w),(1,0)), (K(1/w),(-1,0))]]\n",
    "\n",
    "def construct_biconnected(n):\n",
    "    \"\"\"\n",
    "    Compute all the rational functions that can occur as the aspect ratio of a rectangle made up of\n",
    "    n rectangles of aspect ratio w.  This assumes that the `resistances` array has been filled out to `n-1`.\n",
    "    \n",
    "    This is accomplished by either joining two polar nets in series or parallel. Or else building\n",
    "    a larger graph by replacing each edge of one of the triconnected graphs with a polar net\n",
    "    \"\"\"\n",
    "    # `s` is the set of possible rational functions\n",
    "    s = set()\n",
    "    # The results. Each element of `l` a tuple consisting of \n",
    "    #  - The rational function\n",
    "    #  - The coarse structure\n",
    "    #     + (1,0) for a single vertical rectangle (-1,0) for a horizontal rectangle.\n",
    "    #     + (2,0) for 2 subrectangles in series (-2,0) for 2 rectangles in parallel.\n",
    "    #     + (k,i) for an arrangement scaffolded by a triconnected graph with k edges\n",
    "    #                 the graph `graph_determinants[k][i]` in particular.\n",
    "    #     + (-k,i) the triconnected scaffold dual to `graph_deteminants[k][i]`\n",
    "    # - The fine structure:  A list of `k` edges in the form (j,i) representing `resistances[j][i]`\n",
    "    l = []\n",
    "    # compute the series parallel compositions\n",
    "    for i in range(1,floor(n/2)+1):\n",
    "        j = n - i\n",
    "        for i1,r1 in enumerate(resistances[i]):\n",
    "            for i2,r2 in enumerate(resistances[j]):\n",
    "                r3 = r1[0] + r2[0]\n",
    "                if r3 not in s:\n",
    "                    s.add(r3)\n",
    "                    l.append((r3,(2,0),[(i,i1),(j,i2)]))\n",
    "                    s.add(K(1)/r3)\n",
    "                    l.append((K(1)/r3,(-2,0),[(i,i1),(j,i2)]))\n",
    "    # Compute the constructions with more complicated triconnected scaffolds\n",
    "    if n >= 5:\n",
    "        for comb in range(5,n+1):\n",
    "            if comb in graph_determinants:\n",
    "                # run through the different ways of covering the `comb` edges with subrectangles\n",
    "                for p in iter(Partitions(n,min_length=comb, max_length=comb)):\n",
    "                    parms = [[[(size,x) for x in pos] for pos in cwr(range(len(resistances[size])), count)]\n",
    "                            for size, count in enumerate(p.to_exp(), start=1) if count != 0]\n",
    "                    for xparms in itprod(*parms):\n",
    "                        # flat parms is a list of subrectangles that can be assigned to the edges\n",
    "                        flat_parms = tuple(x for y in xparms for x in y)\n",
    "                        # twist parms is a specific permutation of `flat_parms`\n",
    "                        #   (`Arrangements` is a list of permutations of a multiset)\n",
    "                        for twist_parms in Arrangements(flat_parms,comb).list():\n",
    "                            # `poly_parms` is the list of rational functions for the chosen rectangles\n",
    "                            poly_parms = [resistances[size][pos][0] for size,pos in twist_parms]\n",
    "                            for det_index, det in enumerate(graph_determinants[comb]):\n",
    "                                # func is the rational function representing the aspect ratio.\n",
    "                                func = det[0](*list(poly_parms))\n",
    "                                func.reduce()\n",
    "                                if func not in s:\n",
    "                                    # check that none of the edges is identically zero\n",
    "                                    if all([f(*poly_parms) for f in det[1]]):\n",
    "                                        s.add(func)\n",
    "                                        l.append((func,(comb,det_index),twist_parms))\n",
    "                                        s.add(K(1)/func)\n",
    "                                        l.append((K(1)/func, (-comb,det_index), twist_parms))\n",
    "    return l\n",
    "   \n",
    "\n",
    "def validate_resistances(ratio, graph, edges):\n",
    "    \"\"\"\n",
    "    INPUT:\n",
    "    - `ratio`  An the aspect ratio of the tiling rectangles as an algebraic number.\n",
    "    - `graph`  A tuple containing the number of edges and index of the scaffolding graph.\n",
    "    - `edges`  A list of tuples containing the size and index of the edges of the scaffold.\n",
    "    \n",
    "    OUTPUT:\n",
    "    `True` if none of the small rectangles has size identically zero\n",
    "    `False` otherwise.\n",
    "    \"\"\"\n",
    "    # series parallel graphs have current flow through all resistors.\n",
    "    if abs(graph[0]) >= 5:\n",
    "        # The compute the proportional flow on each edge\n",
    "        poly_parms = [resistances[size][pos][0](ratio) for size,pos in edges]\n",
    "        g = graph_determinants[abs(graph[0])][graph[1]]\n",
    "        # Check that none of the edges have zero flow.\n",
    "        if not all(f(*poly_parms) for f in g[1]):\n",
    "            return False\n",
    "    for size, index in edges:\n",
    "        # recursively check the smaller rectangles.\n",
    "        if size >= 5:\n",
    "            r = resistances[size][index]\n",
    "            if not validate_resistances(ratio, r[1], r[2]):\n",
    "                return False\n",
    "    return True\n",
    "\n",
    "def count_roots(n):\n",
    "    \"\"\"\n",
    "    INPUT:\n",
    "    - `n` The number of edges.\n",
    "    \n",
    "    OUTPUT:\n",
    "      A dictionary mapping real positive algebraic numbers such that a square can be tiled by `n` rectangles\n",
    "      having the algebric aspect ratio to an index in the resistances array describing the arrangement of the\n",
    "      rectangles.\n",
    "    \"\"\"\n",
    "    # The algebraic numbers rejected because they cause degenrate rectangles.\n",
    "    rej = set()\n",
    "    # The algebraic number giving valid soulutions.\n",
    "    s = dict()\n",
    "    foundone = False\n",
    "    for index,r in enumerate(resistances[n]):\n",
    "        f = r[0].numerator() - r[0].denominator()\n",
    "        # The real roots where the aspect ratio is one (a square)\n",
    "        for rt,_ in f.roots(AA):\n",
    "            if rt not in s and rt >= AA(1):\n",
    "                #if it is new and greater than or equal to one ...\n",
    "                if abs(r[1][0]) >= 5:\n",
    "                    # ... and does not cause a degenerate rectangle\n",
    "                    if not validate_resistances(rt, r[1], r[2]):\n",
    "                        if rt >= AA(1):\n",
    "                            rej.add(rt)\n",
    "                        continue\n",
    "                s[rt] = index\n",
    "    print(n, \" \", [x for x in rej if x not in s], \" \", len(s))\n",
    "    return s\n",
    "\n",
    "max_rectangles = 11\n",
    "root_collection = [dict(), count_roots(1)]\n",
    "for i in range(2, max_rectangles):\n",
    "    l = construct_biconnected(i)\n",
    "    print(i, \" \", len(l)/2)\n",
    "    resistances.append(l)\n",
    "    s = count_roots(i)\n",
    "    root_collection.append(s)\n",
    "    \n",
    "print([len(resistances[i])/2 for i in range(1,max_rectangles)]) \n",
    "\n",
    "print([len(root_collection[i]) for i in range(1, max_rectangles)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c08da4f5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.0.beta3",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}