Very Large Resistor Lattice, as in xkcd 356

vlrl-356.py

import sys
import numpy as np
import csv
import logging
from math import ceil
logging.basicConfig(level=logging.WARNING)

OUTFILE = "vll.csv"
RESISTANCE = 1
CURRENT = 1

# Netlist generator
def formattuple(coords):
    return "i{}j{}".format(coords[0],coords[1])

def component(ctype, value, acoord, bcoord):
    anode = formattuple(acoord)
    bnode = formattuple(bcoord)
    assert ctype in ['A','R']
    if ctype == 'R':
        typename = 'r'
        typecode = 'R'
    if ctype == 'A':
        typename = 'a'
        typecode = 'A'
    name = "{}{}.{}".format(typename, anode, bnode)
    return "{},{},{},{},{}\n".format(name, typecode, value, anode, bnode)

def current(acoord, bcoord):
    return component('A', RESISTANCE, acoord, bcoord)

def resistor(acoord, bcoord):
    return component('R', CURRENT, acoord, bcoord)

# CSV parsing
NCOL = 0 # component name
TCOL = 1 # type of component
VCOL = 2 # value of component, eg resistance
ACOL = 3 # node connected to first lead, currents enter here
BCOL = 4 # node connected to second lead

NODE_TYPES = ["A", "R"]

def find_ground_node(degrees):
    ground = max(degrees.keys(), key=(lambda x: degrees[x]))
    logging.debug("ground node-> {}".format(ground))
    return ground

def read_netlist(netlist_path):
    components = {}
    # We will need to iterate over components twice
    # in the same order, so we save keys
    component_keys = []
    with open(netlist_path, 'r') as infile:
        netlist = csv.reader(infile)
        nums = {}
        nums["components"] = 0
        nums["anomalies"] = 0
        nums["be"] = 0  # number of branch equations
        nums["kcl"] = 0 # number of non-ground nodes
        degrees = {}
        anomnum = {}
        for component in netlist:
            key = component[NCOL]
            assert type(key) == str
            component_keys.append(key)
            components[key] = [None] * 8

            newcomp = components[key]
            newcomp[NCOL] = key
            assert component[TCOL] in NODE_TYPES
            newcomp[TCOL] = component[TCOL]

            try:
                newcomp[VCOL] = float(component[VCOL])
            except ValueError:
                print((
                    "Bad input: expected a number"
                    "for component value of {},"
                    "got {} instead.").format(
                        component[NCOL],
                        component[VCOL])
                    )
                exit(1)

            newcomp[ACOL] = component[ACOL]
            newcomp[BCOL] = component[BCOL]
            assert len(component) == 5

            nums["components"] += 1
            curnodes = component[ACOL:BCOL+1]
            newnodes = [key for key in curnodes
                    if key not in degrees]
            for node in newnodes:
                degrees[node] = 0
            for node in curnodes:
                degrees[node] += 1

    ground = find_ground_node(degrees)

    i = 0
    nodenum = {}
    for node in [k for k in degrees.keys() if k != ground]:
        nodenum[node] = i
        i += 1
    assert len(nodenum) == len(degrees) - 1
    logging.debug("nodenum={}".format(nodenum))
    nums["kcl"] = len(nodenum)
    nums["be"] = nums["anomalies"]
    logging.debug("nums={}".format(nums))
    logging.debug("anomnum={}".format(anomnum))
    # From now on nums shall become immutable

    state = [nums, degrees, anomnum, components, component_keys, ground, nodenum] 
    return state

def build_coefficients(state):
    [nums, degrees, anomnum, components, component_keys, ground, nodenum] = state
    n = nums["kcl"] + nums["be"] # number of unknowns
    G = np.zeros(shape=(n, n))
    A = np.zeros(shape=(n, 1))
    currents = []
    for key in component_keys: # preserve order of iteration
        component = components[key]
        anode = component[ACOL]
        bnode = component[BCOL]
        if anode != ground:
            i = nodenum[anode]
            assert i >= 0 and i <= nums["kcl"]
        if bnode != ground:
            j = nodenum[bnode]
            assert j >= 0 and j <= nums["kcl"]
        if component[TCOL] == "R":
            try:
                conductance = 1 / component[VCOL]
            except ZeroDivisionError:
                print("Model error: resistors can't have null resistance")
                exit(1)
            if anode != ground:
                G[i,i] += conductance
            if bnode != ground:
                G[j,j] += conductance
            if bnode != ground and anode != ground:
                G[i,j] -= conductance
                G[j,i] -= conductance
        elif component[TCOL] == "A":
            current = component[VCOL]
            if anode != ground:
                A[i] += current
            if bnode != ground:
                A[j] -= current
        else:
            exit(1)

    logging.debug("currents={}".format(currents))
    logging.debug("G=\n{}".format(G))
    logging.debug("A=\n{}".format(A))
    return [G, A, currents]

def solve_system(G, A):
    try:
        e = np.linalg.solve(G, A)
    except np.linalg.linalg.LinAlgError:
        print("Model error: matrix is singular")
        logging.error(G)
        exit(1)
    return e

def print_solution(e, nodenum, nums, currents):
    nodelabel = dict(reversed(x) for x in nodenum.items())
    print("Ground node: {}".format(ground))
    i = 0
    for potential in e[0:nums["kcl"]]:
        print("e({}) = {}".format(nodelabel[i], potential))
        i += 1
    i = 0
    for current in e[nums["kcl"]:nums["kcl"]+nums["be"]]:
        print("i({}) = {}".format(currents[i], current))
        i += 1

def experiment(n):
    assert n >= 3
    assert n % 2 == 1
    I = n
    J = I + 1
    with open(OUTFILE, 'w') as outfile:
        for i in range(I):
            for j in range(J):
                here = (i,j)
                right = (i,j+1)
                left = (i,j-1)
                up = (i-1,j)
                down = (i+1, j)
                if j < J-1:
                    outfile.write(resistor(here, right))
                if i < I-1:
                    outfile.write(resistor(here, down))

        c = ceil(I/2)
        anode = (c, c-1)
        bnode = (c-1,c+1)
        outfile.write(current(anode, bnode))

    state = read_netlist(OUTFILE)
    [nums, degrees, anomnum, components, component_keys, ground, nodenum] = state

    [G, A, currents] = build_coefficients(state)
    e = solve_system(G, A)

    aname = formattuple(anode)
    bname = formattuple(bnode)
    if aname == ground:
        ea = 0
    else:
        anum = nodenum[aname]
        ea = e[anum]
    if bname == ground:
        eb = 0
    else:
        bnum = nodenum[bname]
        eb = e[bnum]

    deltav = ea - eb
    equivalent_resistance = abs(deltav) / CURRENT

    return equivalent_resistance[0]

if __name__ == '__main__':
    i = int(sys.argv[1])
    print("{},{}".format(i, experiment(i)))