From d776fd0ceac356756f568ef4a838932b30db13f8 Mon Sep 17 00:00:00 2001 From: Nick Downing Date: Fri, 18 Jul 2025 22:32:37 +1000 Subject: [PATCH] In /scripts/net_channels.py add resistor solver to determine channel resistance --- scripts/net_channels.py | 206 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 192 insertions(+), 14 deletions(-) diff --git a/scripts/net_channels.py b/scripts/net_channels.py index 1336d66..3a00b44 100755 --- a/scripts/net_channels.py +++ b/scripts/net_channels.py @@ -1,8 +1,22 @@ #!/usr/bin/env python3 +#import PIL.Image import numpy +import scipy.sparse +import scipy.sparse.linalg import sys +# palette for resistor solver diagnostics (usually commented) +#PALETTE = numpy.array( +# [ +# [0x00, 0x00, 0x00], # colour 0 = black (insulator) +# [0x00, 0x00, 0xff], # colour 1 = blue (sink) +# [0xff, 0x00, 0x00], # colour 2 = red (source) +# [0xff, 0xff, 0xff], # colour 3 = white (conductor) +# ], +# numpy.uint8 +#) + if len(sys.argv) < 3: print( f'usage: {sys.argv[0]:s} nets.txt poly_colour,diff_colour,buried_colour,channel_colour' @@ -16,20 +30,23 @@ diff_colours = {diff_colour, buried_colour} with open(nets_txt) as fin: line = fin.readline() assert len(line) + _, ys, xs = [int(i) for i in line.split()] line = fin.readline() assert len(line) n_items = int(line.rstrip()) items = [] - channel_nets = {} + channel_items = {} for i in range(n_items): line = fin.readline() assert len(line) net, colour, y, x0, x1 = [int(i) for i in line.split()] + if colour == channel_colour: + if net not in channel_items: + channel_items[net] = ([], []) # items, reports (report may be swapped) + channel_items[net][0].append(len(items)) items.append((net, colour, y, x0, x1)) - if colour == channel_colour and net not in channel_nets: - channel_nets[net] = [] line = fin.readline() assert len(line) @@ -41,21 +58,182 @@ with open(nets_txt) as fin: item0, item1 = [int(i) for i in line.split()] net0 = items[item0][0] net1 = items[item1][0] - if net0 in channel_nets: - channel_nets[net0].append(item1) - if net1 in channel_nets: - channel_nets[net1].append(item0) - -for channel_net, adjacent_items in sorted(channel_nets.items()): - #adjacent_nets = set([items[i][0] for i in adjacent_items]) - #adjacent_colours = set([items[i][1] for i in adjacent_items]) + if net0 in channel_items: + channel_items[net0][1].append((item0, item1)) + if net1 in channel_items: + channel_items[net1][1].append((item1, item0)) + +def resistor_solver(image): + 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 + # (when node is electrode, positive I means current is entering the array) + # resistor R connected n1 -> n2: + # causes a current of (V(n1) - V(n2)) / R, positive when V(n1) > V(n2) + # => when current is positive, it should increase I(n1) and decrease I(n2) + # => I(n1) += 1/R V(n1) - 1/R V(n2), I(n2) -= 1/R V(n1) - 1/R V(n2) + v = [] + i = [] + j = [] + for n1 in numpy.nonzero(image_horiz.reshape((n,)))[0]: + n2 = n1 + 1 + v.extend([1, -1, -1, 1]) + i.extend([n1, n1, n2, n2]) + j.extend([n1, n2, n1, n2]) + for n1 in numpy.nonzero(image_vert.reshape((n,)))[0]: + n2 = n1 + xs + v.extend([1, -1, -1, 1]) + i.extend([n1, n1, n2, n2]) + j.extend([n1, n2, n1, n2]) + v = numpy.array(v, numpy.int32) + i = numpy.array(i, numpy.int32) + j = numpy.array(j, numpy.int32) + A = scipy.sparse.coo_array( + (v, (i, j)), + shape = (n, n), + dtype = numpy.int32 + ).tocsr() + + # partition nodes into groups: + # 0: black (insulator) + # 1: blue (sink) + # 2: red (source) + # 3: white (conductor) + g0 = (image == 0).reshape((n,)) + n0 = numpy.sum(g0) + g1 = (image == 1).reshape((n,)) + n1 = numpy.sum(g1) + g2 = (image == 2).reshape((n,)) + n2 = numpy.sum(g2) + g3 = (image == 3).reshape((n,)) + n3 = numpy.sum(g3) + #print('n0', n0, 'n1', n1, 'n2', n2, 'n3', n3) + + # partition A into 4 x 4 submatrices by groups, V and I into subvectors + #A00 = A[g0, :][:, g0] + #A01 = A[g0, :][:, g1] + #A02 = A[g0, :][:, g2] + #A03 = A[g0, :][:, g3] + #A10 = A[g1, :][:, g0] + #A11 = A[g1, :][:, g1] + #A12 = A[g1, :][:, g2] + #A13 = A[g1, :][:, g3] + #A20 = A[g2, :][:, g0] + #A21 = A[g2, :][:, g1] + A22 = A[g2, :][:, g2] + A23 = A[g2, :][:, g3] + #A30 = A[g3, :][:, g0] + #A31 = A[g3, :][:, g1] + A32 = A[g3, :][:, g2] + A33 = A[g3, :][:, g3] + + # solve: + # A00 V0 + A01 V1 + A02 V2 + A03 V3 = I0 + # A10 V0 + A11 V1 + A12 V2 + A13 V3 = I1 + # A20 V0 + A21 V1 + A22 V2 + A23 V3 = I2 + # A30 V0 + A31 V1 + A32 V2 + A33 V3 = I3 + # where V0 = all 0, V1 = all 0, V2 = all 1, I0 = all 0, I2 = all 0 + # => unknowns are I1, I2, V3 + #A02_V2 = numpy.sum(A02.astype(numpy.int32), 1) + #A12_V2 = numpy.sum(A12.astype(numpy.int32), 1) + A22_V2 = numpy.sum(A22.astype(numpy.int32), 1) + A32_V2 = numpy.sum(A32.astype(numpy.int32), 1) + #print('solve') + V3 = scipy.sparse.linalg.spsolve( + A33.astype(numpy.double), + -A32_V2.astype(numpy.double) + ) + #print('solve done') + + # we only care about sum(I1) and sum(I2), and we know sum(I1) = -sum(I2) + sum_I2 = numpy.sum(A22_V2, 0) + numpy.sum(A23, 0) @ V3 + #print('R', 1. / sum_I2, 'squares') + return 1. / sum_I2 + +image = numpy.zeros((ys, xs), numpy.uint8) +for ( + channel_net, + (channel_items, channel_reports) +) in sorted(channel_items.items()): + #adjacent_nets = set([items[item1][0] for _, item1 in channel_reports]) + #adjacent_colours = set([items[item1][1] for _, item1 in channel_reports]) #print('channel_net', channel_net, 'adjacent_nets', adjacent_nets, 'adjacent_colours', adjacent_colours) poly_nets = set( - [items[i][0] for i in adjacent_items if items[i][1] == poly_colour] + [ + items[item1][0] + for _, item1 in channel_reports + if items[item1][1] == poly_colour + ] ) diff_nets = set( - [items[i][0] for i in adjacent_items if items[i][1] in diff_colours] + [ + items[item1][0] + for _, item1 in channel_reports + if items[item1][1] in diff_colours + ] ) - if len(poly_nets) != 1 or len(diff_nets) != 2: + if len(poly_nets) == 1 and len(diff_nets) == 2: + # find channel resistance by creating image for resistor solver + # conductor is filled in colour 3, then electrodes in colours 1 and 2 + # bounding box is tracked in order to pass the smallest possible image + # into the solver, and to clear this region from the large image after + min_x = xs + max_x = 0 + min_y = ys + max_y = 0 + for item in channel_items: + _, _, y, x0, x1 = items[item] + if x0 < min_x: + min_x = x0 + if x1 > max_x: + max_x = x1 + if y < min_y: + min_y = y + if y + 1 > max_y: + max_y = y + 1 + image[y, x0:x1] = 3 # white (conductor) + sink_net, source_net = tuple(diff_nets) + electrodes = {sink_net: 1, source_net: 2} + for item0, item1 in channel_reports: + net, colour, y1, x2, x3 = items[item1] + if colour in diff_colours: + # trim neighbouring item to only include direct neighbouring pixels + _, _, y0, x0, x1 = items[item0] + if x2 < x0 - 1: + x2 = x0 - 1 + if x3 > x1 + 1: + x3 = x1 + 1 + + if x2 < min_x: + min_x = x2 + if x3 > max_x: + max_x = x3 + if y1 < min_y: + min_y = y1 + if y1 + 1 > max_y: + max_y = y1 + 1 + image[y1, x2:x3] = electrodes[net] + + # diagnostics (usually commented), need to make directory first + #image_pil = PIL.Image.new('P', (max_x - min_x, max_y - min_y), None) + #image_pil.frombytes(image[min_y:max_y, min_x:max_x].tobytes()) + #image_pil.putpalette(list(PALETTE.reshape((12,)))) + #image_pil.save(f'channels/{channel_net:d}.png') + + R = resistor_solver(image[min_y:max_y, min_x:max_x]) + image[min_y:max_y, min_x:max_x] = 0 + + print('channel_net', channel_net, 'poly_nets', poly_nets, 'diff_nets', diff_nets, 'R', R) + else: print('channel_net', channel_net, 'poly_nets', poly_nets, 'diff_nets', diff_nets) -- 2.34.1