Initial commit, uses dense method to solve planar resistance graph problem

diff --git a/resistor_solver.py b/resistor_solver.py
new file mode 100755 (executable)
index 0000000..af87228
--- /dev/null
+++ b/resistor_solver.py
@@ -0,0 +1,115 @@
+#!/usr/bin/env python3
+
+# see https://github.com/Zofz/Resist2D.git
+
+import PIL.Image
+import numpy
+import sys
+
+PIL.Image.warnings.simplefilter('ignore', PIL.Image.DecompressionBombWarning)
+
+if len(sys.argv) < 3:
+  print(
+    f'usage: {sys.argv[0]:s} image_in image_out'
+  )
+  sys.exit(1)
+image_in = sys.argv[1]
+image_out = sys.argv[2]
+
+image_pil = PIL.Image.open(image_in)
+image_pil.load()
+#print('mode', image_pil.mode)
+#print('palette', image_pil.palette)
+#print('width', image_pil.width)
+#print('height', image_pil.height)
+image = numpy.array(image_pil)
+#print('shape', image.shape)
+#print('dtype', image.dtype)
+assert len(image.shape) == 2 and image.dtype == numpy.uint8
+assert numpy.all(image < 4)
+# colours are 0 = black, 1 = blue, 2 = red, 3 = white
+#print('xxx', image_pil.palette.tobytes()[:12])
+ys, xs = image.shape
+n = ys * xs
+
+# image_horiz has true pixel for left terminal of each horizontal resistor
+# image_vert has true pixel for top terminal of each vertical resistor
+image_nz = image != 0
+image_horiz = numpy.zeros_like(image_nz)
+image_horiz[:, :-1] = numpy.logical_and(image_nz[:, :-1], image_nz[:, 1:])
+image_vert = numpy.zeros_like(image_nz)
+image_vert[:-1, :] = numpy.logical_and(image_nz[:-1, :], image_nz[1:, :])
+
+# KCL is enforced as follows:
+# A V = I where V = voltage at node, I = net current out of node
+# I is current "unaccounted for", i.e. current that cancels R currents
+# resistor R connected n1 -> n2:
+#   causes a current of (V(n2) - V(n1)) / R, positive when V(n2) > V(n1)
+#   => when current is positive, it should decrease I(n1) and increase I(n2)
+#   => I(n1) -= 1/R V(n2) - 1/R V(n1), I(n2) += 1/R V(n2) - 1/R V(n1)
+#   => I(n1) += 1/R V(n1) - 1/R V(n1), I(n2) += -1/R V(n1) + 1/R V(n2)
+A = numpy.zeros((n, n), numpy.int8)
+for n1 in numpy.nonzero(image_horiz.reshape((n,)))[0]:
+  n2 = n1 + 1
+  A[n1, n1] += 1
+  A[n1, n2] -= 1
+  A[n2, n1] -= 1
+  A[n2, n2] += 1
+for n1 in numpy.nonzero(image_vert.reshape((n,)))[0]:
+  n2 = n1 + ys
+  A[n1, n1] += 1
+  A[n1, n2] -= 1
+  A[n2, n1] -= 1
+  A[n2, n2] += 1
+
+# partition nodes into 3 groups:
+# 0: black + blue (electrode), 1: red (electrode), 2: white (conductor)
+# note: black pixels are part of electrode to give them a defined voltage
+g0 = (image < 2).reshape((n,))
+n0 = numpy.sum(g0)
+g1 = (image == 2).reshape((n,))
+n1 = numpy.sum(g1)
+g2 = (image == 3).reshape((n,))
+n2 = numpy.sum(g2)
+print('n0', n0, 'n1', n1, 'n2', n2)
+
+# partition A into 3 x 3 submatrices by groups, V and I into subvectors
+#A00 = A[g0, :][:, g0]
+#A01 = A[g0, :][:, g1]
+#A02 = A[g0, :][:, g2]
+#A10 = A[g1, :][:, g0]
+A11 = A[g1, :][:, g1]
+A12 = A[g1, :][:, g2]
+#A20 = A[g2, :][:, g0]
+A21 = A[g2, :][:, g1]
+A22 = A[g2, :][:, g2]
+
+# solve:
+# A00 V0 + A01 V1 + A02 V2 = I0
+# A10 V0 + A11 V1 + A12 V2 = I1
+# A20 V0 + A21 V1 + A22 V2 = I2
+# where V0 = all 0, V1 = all 1, I2 = all 0 => unknowns are I0, I1, V2
+#A01_V1 = numpy.sum(A01.astype(numpy.int32), 1)
+A11_V1 = numpy.sum(A11.astype(numpy.int32), 1)
+A21_V1 = numpy.sum(A21.astype(numpy.int32), 1)
+print('solve')
+V2 = numpy.linalg.solve(
+  A22.astype(numpy.double),
+  -A21_V1.astype(numpy.double)
+)
+print('solve done')
+
+# we only care about sum(I0) and sum(I1), and we know sum(I0) = -sum(I1)
+sum_I1 = numpy.sum(A11_V1, 0) + numpy.sum(A12, 0) @ V2
+print('R', 1. / sum_I1, 'squares')
+
+V = numpy.zeros((n,), numpy.double)
+V[g1] = 1.
+V[g2] = V2
+
+image = V.reshape((ys, xs))
+mask = image >= .04045
+image[~mask] /= 12.92
+image[mask] = ((image[mask] + .055) / 1.055) ** 2.4
+image = numpy.round(image * 0xff).astype(numpy.uint8)
+PIL.Image.fromarray(image).save(image_out)

diff --git a/schur_complement.py b/schur_complement.py
new file mode 100755 (executable)
index 0000000..d1d02eb
--- /dev/null
+++ b/schur_complement.py
@@ -0,0 +1,36 @@
+#!/usr/bin/env python3
+
+import numpy
+import numpy.linalg
+
+A = numpy.array(
+  [
+    # N1   N2   NA   NB
+    [ 1.,  0.,  0., -1.], # R1
+    [ 0.,  1., -1.,  0.], # R2
+    [ 1., -1.,  0.,  0.], # R3
+    [ 0.,  0., -1.,  1.], # R4
+    [ 0., -1.,  1.,  0.], # R5
+    [ 0.,  1.,  0., -1.], # R6
+    [ 1.,  0., -1.,  0.], # R7
+  ],
+  numpy.double
+)
+C = numpy.diag(
+  1. / numpy.array([1., 2., 3., 4., 5., 6., 7.], numpy.double)
+)
+
+W = A.transpose() @ C @ A
+
+# set last 2 node voltages, solve for last 2 node currents
+# uses Schur complement which allows us to have an n - 2 by n - 2 inversion
+P = W[:-2, :-2]
+Q = W[-2:, :-2] # == W[:-2, -2:].transpose()
+R = W[-2:, -2:]
+excite = numpy.array([1., 0.], numpy.double)
+print((R - Q @ numpy.linalg.inv(P) @ Q.transpose()) @ excite)
+
+# set last 2 node currents, solve for all node voltages (last node = GND)
+I = numpy.array([0., 0., 1., -1], numpy.double)
+V = numpy.linalg.solve(W[:-1, :-1], I[:-1])
+print('V', V)

diff --git a/stillOKfor32.png b/stillOKfor32.png
new file mode 100644 (file)
index 0000000..4b51a32
Binary files /dev/null and b/stillOKfor32.png differ

diff --git a/stillOKfor32_100x100.png b/stillOKfor32_100x100.png
new file mode 100644 (file)
index 0000000..bac61b8
Binary files /dev/null and b/stillOKfor32_100x100.png differ

diff --git a/stillOKfor32_100x100_trimmed.png b/stillOKfor32_100x100_trimmed.png
new file mode 100644 (file)
index 0000000..368a3d0
Binary files /dev/null and b/stillOKfor32_100x100_trimmed.png differ

diff --git a/stillOKfor32_50x50.png b/stillOKfor32_50x50.png
new file mode 100644 (file)
index 0000000..7f8b9dd
Binary files /dev/null and b/stillOKfor32_50x50.png differ