Last active
March 29, 2019 18:02
-
-
Save EnricoMiccoli/d7e5ead0825523aa4daa5bc75e5ed62f to your computer and use it in GitHub Desktop.
Very Large Resistor Lattice, as in xkcd 356
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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))) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment