From fe6f9bff899abb35a6244dbf17da7ebca513487a Mon Sep 17 00:00:00 2001 From: Nick Downing Date: Wed, 16 Jul 2025 13:45:38 +1000 Subject: [PATCH] Initial commit, uses dense method to solve planar resistance graph problem --- resistor_solver.py | 115 +++++++++++++++++++++++++++++++ schur_complement.py | 36 ++++++++++ stillOKfor32.png | Bin 0 -> 1900 bytes stillOKfor32_100x100.png | Bin 0 -> 5163 bytes stillOKfor32_100x100_trimmed.png | Bin 0 -> 5862 bytes stillOKfor32_50x50.png | Bin 0 -> 5059 bytes 6 files changed, 151 insertions(+) create mode 100755 resistor_solver.py create mode 100755 schur_complement.py create mode 100644 stillOKfor32.png create mode 100644 stillOKfor32_100x100.png create mode 100644 stillOKfor32_100x100_trimmed.png create mode 100644 stillOKfor32_50x50.png diff --git a/resistor_solver.py b/resistor_solver.py new file mode 100755 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 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 index 0000000000000000000000000000000000000000..4b51a32cdad44fcfb8d21a4faffa2cc2952d9b09 GIT binary patch literal 1900 zcmeAS@N?(olHy`uVBq!ia0vp^CqS5k8A#4*i(3PvSkfJR9T^xl_H+M9WMyDrW(e>J zaRrJqi2nc2!tnq9e-L{VjE2BS3IUmL$E!g9au#?*7BevL9Ry*<9TT(P0tF>XTq84D?J>Ta8kIndqLUi(^Pd+}rbS%N|+ruwK~q{qOys zOm{Do3DqYk9bdEbZ~8K~iRXTLXfB-+
    PKHMK~vsan--JX5&KkF;%%`g3bxO}_P zFZ(|~+*s~^`22f9b$g?^;uk&T_m0n)oZmY<5(>_Dc%YR0dh+%|4_+-m`sENyuJL{yL+c+2p+Vn199`7%=J9OR&KoD)wI^+ueg2wqy#0M$_dV)A z3QRBA9$lV2@vY^)*scu!PanNhs%M?Qu6X@(&3>MZZZ|e?EUc>e`Tby*>4_`J!6_@% zl)l~&i(J{|VEe6W&1IK)a=TrpWh(i#uV4RHtvf_;ne)aIQoboK1N1hzhC(C2 zGGUu{|I*nfE?(0psuI<_HOtj6Pst^E_lcdpDH*Rdin2sGXXV+3Tegm>nHvnuAck%vPHDw z-~76j#R+2J-7f#9$G@I4)Wb=zM0MpKDss{ueC^PR(7Kd*Yg;e@<%(xkIj^f;dyjI z#P3)^GM~A#g<Fnm6bqd|fi~m`C zo9EE5_(l~QCtuP_%M^>2sp^7{ILzj`b}>g($XUjv-R)$loSHkirta5y`9*t4Ki~Be;}S2+993cu*_hsy-M>NYnopqar!3*)#R^FBq1TBluD9P?zU zTA0s*W3BDFGut8z8^y%swxviVx!S1a+2o!uInpJNWG=C9Qf7y!q03sMD&@-p+DRv( zQopz@Qw(Df+;-@JIZ&IR%Pp?v`yg#1h21@hKsAd_M5pk8WZWJwIhPypZdB?y!rJx@ zq(o+7xn$)8pn-)c_fnrYg3J-BYHySWDVlJ_n&%Hlu5rn|#0Tv@&TBf9foYW=B-eLl zcaIIoMFvXZ61ye?6(y?Z8I-7lY~zT5*yegqDH)_h;Zc~xA2py%lSg#ohuI)=6mzCt zlm@xT`I^>ZbC7{8D|p-1gZ;WOAL3WeDLX!FPz7pJyYT4H;b}l^GaP%1`gMv@kBY37 zi+pmb%TmgA(}}5~ldn{#Ow`uCrf^xpI!xiBf;W)Rh(4!&-lF=P+F9)!tLZxHH@^0n rwCBvO%^nk>c@vRKHJ2*a*8Jh0^?1fj{`8v_pd!)J)z4*}Q$iB}xqAi; literal 0 HcmV?d00001 diff --git a/stillOKfor32_100x100.png b/stillOKfor32_100x100.png new file mode 100644 index 0000000000000000000000000000000000000000..bac61b80a35aad35ecb0e3abb56970149b5d7d9a GIT binary patch literal 5163 zcmeHKc~BeY6&GNNEEB-7!4RiqEtmvaue92ggtSQJkc1@{VQ~aZe3aGh3f4jgs|88K z*s&qRS7MxcI;lhI(Au~$@d3E4UD7dN4hOq+L&tWVaWYJgfI}QNPHNmbF8x-*9MjBn zm}&khJG0vT-tl|y{oeb&(eAP2t$jlDlt>^DJYh2G3&4}&Z($htU8|vO;4xTRw1F!i zypWr9*_cWi;{0wJqJ0cd`aW=2Hbz~9A|H+}PY?CSjHtLfM?WkZJhJob*0$$+8&s^Z z_x{7}y>0cT`qoYF95!UO6+4CNnk(v4KN$%ryfJn-`sHrpV&{*qTDM#ny^XwHqC1qA z^Vd7QZ0w1&5f;0bGrXfuSe$mknXCzC>#uJ+^!z=41HG^0kU#tQwLhKtzJKe)uWXer z^H|cAH=6WkcKmAm7$Eq3Ppq_?~)&z*BJL*H?q4~x_!bfk6_G__h%1b557G^^j5j2qkX_gedB zfBrJ$lBMMBl_MiD@o>TY!&keyK6v5bL~iGW#DQP3myYxgzYj*w$AA&rV9v!!mqS8O zE-NkZIox3O1OiQl&rOhanuDx#Ipfrd?{-}fLky)Amnh7L*{!217-KC<7uMz#k+pVG zO^Gwsh%`PN061ulfP9Whrw8|G#Q|I#yz^qI7z&thcCC1W*#hZYEDgyeatQ((d`t~0 zUL%4uEM>zB^f^-yV5JpTaGV>LO1)mM#EVH>Y`IjXR;#55Dn(HkSil~?lOuev)04zQ z1Tplqhh!Ny$GDsjk4acv)tpu=272gmeh#R3<^B4u^EMhlexN z0FbGKzShH21a_*lfcChmS&}x?&`vICHUvdZ`@5^z%0M|3DWxlE2Qc-3sIoaB*P6_h zX%C))a>n5fcmcBKSaOVQhO9Zc@tQz6vl9XC)3|f2AJ-l*23BS>u6L2ue0nCmR?Nr8 zDHq95c;Hn@TV5_`g#wn7G7MIe6bh$dh*~a35xJ72XF-{q9*%I5G!F&9B@Dnx zLkKmgR%0+m(>7QxM`+lp28dRb)kY|-GDMDIvmn;9446uya&}fc6a}DEHd<-3VhY%% zRHsX$Rwt(4Cg71OK-Od^jdLnK&G!?y(o z#QsXrLTx^lcRjC#8DM7sFwAu`|6Gd=@p5SNzil7P{NyBBzA{mBb zava6JCqr;#mc5H&Z2tdd&2Jt^6KuJW@qqaLfGD`53hAogD!8m<0$T}!0-FLS$Y2Q` zqK2jdc>=7UimV`<!GQck=1BVwl3#HSC;S^^e-*57o>cux10)iHT%!}Uzx)$h~7X$M$UYM>0 zy5_~eyo?v7>wiXs7S}w>hyMLj4*&}ATosleC<&EoJ44*FOSbe#50NeR`o7@SUM+|}}8cUelj zNEq1|8d~2{mI5L1qQ>FG^Zn?Gn2({8{o6b3LQ!wz(U&Bh+p)P; YK(2XqP1oxDRebABhCKbL?6O_|23b2fF8}}l literal 0 HcmV?d00001 diff --git a/stillOKfor32_100x100_trimmed.png b/stillOKfor32_100x100_trimmed.png new file mode 100644 index 0000000000000000000000000000000000000000..368a3d0b1f05ffd32ad52dae3320e283bf825041 GIT binary patch literal 5862 zcmeHLd2kcg8J8~{93T!gV2EQj!VU(y(%#Z;CGa61qLvZ1!8S+m9&hEv*5TF47A71M z2*spag(Q&BW+)IyNYbPYQ%cHE6R0T$LQ}#}P6ucbU?6mo;Rpoxtz>*mGt*%P@=v{* z)$V)md%yR6zwh_$)2i8BQaCVaM3TW^80ai=l)|f|pO3}C??M}p;ia{yyizOW8cbM2bvQ{Mg-6g+hv5ok>M!fAG z_R99^@e8jUK9J}+B)R+|L|KdY4-_=!>cQkB-Yt%}y4V9ymmNkufIW+`mIB4YR&JZ)@u~{l=zAho=vCZ0Wu~9KyWk z2cEy=I@SK?NogRnDs%kkbjzg|_N|;zVVi$phE~K$^{p<0nl~ayP->_ri@*UUyHY@eWm2qt|wz9v>82Ukce`CrqcR>Pk{=>_gD~tXy ze{buAi{E9cV=j(ev|5arQnPC2gj*x;UOPUvxFqcSJpSnS=dB}0W!DU=IsWvQcTUv5 zHah$I@SUHYx;Hy+-R+FsXO9bo*&R%v{><^U_I5Qf^?5=nYuJLgHo6UGKl$o)HT&Yw zuV#MJRD1Qpw=oCY6>mM&-kzL_mfqO)$wza!xp!_BZ{IUg%~IcNF0n(oB|oD;g;lzW zInnPm@seKvrm!~vrDrhMCWHgLSPL{n03O9>H-7!m9wVYic4GzY!d!toP@@zzsi3T> zq+D#O60;EZ8jni6Rp;6taQ1NJiiJM%0!tk zG(W5~5XSLIh)tDbuGCS`2?1MnV~wT-IJ3F2vC-5>n*6HAjI%6j#t1V(pwI#h&hcq{ z81)55=@3y22MCI)63`UC579As!5`A>Mk5?Y?)&ErxLjTEzF;Q{kPma1514TiX7+l` zJ!b^9{00cp>Cg|(2$sV_HJ5^*KctEvzXABPQ9U6fv1@)Hq}D~kkwi17172tvgsbAc zmMnC-++8zt3OtH85SazZ?q#VdayMDMeA7oF;q-I_PVd6)Wqm*Ph%vNsxj2Vk4C(GU z9d@I>J}3D_MdBj23=nY9N)aeyp=p#7aS~-ki9jtR#!?i4QH%(BKskLujrWN_hl1cH z1>#sRNfdZmL@^*rC?%6DYGq{tWi1xWN?BM*!pR;GlT-z&lCSIOl@29AD60$@Ss-at zW~?kqNrD8uVQ`eBC4!PfMr2t$0wsxDfnW9VFrA8*_W*Oi=ZOsHgmXD=r`<@Hu zK#KSYv7$p_4e#>+`1I(M>-~P^hok@~LKJl&2{M!dDG8_`T6t6^ahf4mOyrqpR9)yn zzpOR#D#-Cb9wApyo)NB)tVp89^!BB(2IyIUgrS%P?Ieui2)dIn^P>+DV+9(gFaZ@= z8cGyUJjyb>j1nR#3xH&;6vp;;`XPkKa0Kz=g-DVNM$toPctb}= zhkoC8^fB=N#6Vg_&Sx;*7~$uI4#U!fr;A{tf1s!;XfOf46f8~xxgk5Kxj-2wEvsX{F4BCqgh&%ezdmRf3mg4=WKc1Lkx%1&qCN)={%-ma_ z*7~Q(hW@dU&)$aX&GDb(x0>TOV(zSk3uEF3q$H&rtUU1CmVy&Iq$I8&`?t9}ClV8z zWBIY`walB?))QB5KCy&YeQq9?e&^*4pq+A_rJiX`ods(Z5cc=>W^6)PQjN6KN z=CwGE_%~OF_UFfb(=c!Q{5{owW!SV<-D6k!4*H4M&6kD8VW0+Qeu?9~-0Bzp14WSc A!vFvP literal 0 HcmV?d00001 diff --git a/stillOKfor32_50x50.png b/stillOKfor32_50x50.png new file mode 100644 index 0000000000000000000000000000000000000000..7f8b9dd435569b018892410edc2276ba946f768b GIT binary patch literal 5059 zcmeHKeNYtV8DA=Ta0cgSF;!|Oo7GX7c5nB-Z&waI?|{3E2M3&jL`V1b6j23TT`rBH3VBTPOM2QEs|s=CPC4Rq>|B=zPkqq)68U; zY5poZv$yX)@9+6N&+~hq_nq4v-jd=hb*>tMpe(n`Sq`3h`OcUDe)B$%wt`1@y{}3v zrxI`^DsX&|g~hrE3$qf>L6CGcjpdhoKME5s*?(akT;mgm_GHFC&itVNlfG@wY^?2heyH}Mu5swooZBxw z{;YOoJj;LLJ08cWg60n&oxihn2iknKzItDJ%Zj@K6@Fvgsn>2u#a}da9Dk}pJNNp& z1bX}=`nP3&?3$ff{cxEDUifX-Z(poh++%i`FQL*uAGdb4-0!`dHk?*kvi7#8h@S^R zX+ctckhiCk)72oh%e)fm9NZM{M;7% zGp&oxpD&A`*5S^ocx2OiSI(@~q%U9b_nAuyznk}Yqx+YIPpP{9B_(ccx;sacJNMw) z^~ahH_HyJN&9M`!bwh2MSsG_f*`Bp;B^Mp&U-pN=?EZoE$-tNU@5P3W)U~Ijbx6DC zAM|wD^H1{MOWxhF@Yjo8Id}csM~l}#J=`(S)-5F4wX^dLe(v2n=d%5se?9Y<>e#L$ z*`?XtKF^^i9v$esP}hEJ&5_RcY+uekTi3wu888+R$+QF6`RoZ0-< z(&6E4>rXkgE$di3`okC0S+5VMMz&u2M|De2CYG0k-rTU@I-Qeo^II3+{OE~CmX+K| z&nTEtS>L-M1R3_f{C;osjf_vO&C(w2{rT&vb%T(jv^99)XwJ)dIqK8*<~LN#nb$S5 z{{61L>U-7(NdK4BpUtvb)8A_Q*~!5>gZtj;`?u}#XSY-dP4uwi=wn+z=Se*1uqw|g zk`_W*iV^&*RtiNx>p@Wdaw$U5wX6vH*#IB5A-8^i3W0gXhE$q7xF_OZYj{_El&z>Q z@zM3QG{GRt3)J}%2?#>0NWoGl7><#W4N-VWFqX|20xJ-))`nDhys$%vvanHW)Z%EN z#BbIk1!_1y%5Y@4vuKn8yxEW%QH+onmPjPD34>OM1~45#5E!n<^m-H^&{$nqq$D&P zTPRbcIGk*Zj`9(a7s9a2N%@7iXhRS%50B*+ig-NZ^x@d33P2A`q9T}1i({b>Hqj#{ z7H$S4qX~VbN6ZH{DptKHy6tdXu z@s4}Q5(M~AMDYS*PeO`3_du*kxydt1ITI5B?&G|Z&||eLV1V*?NT)!@<@DT68zRRi z8G+^*QW+Wv1I^)P93?1=2{p1N9ZL9VJ!;|zhA=T^z0sS}B@A zjg;Ps5_-;nTJ=_%VHlmCV{nCvp~)g48lu2*@*yg~Vv%q_nUDo1?OwMH(Q8vs#o`T8 zA_p97$SOV@mmWZUe2A?ODOpXO#iTbI%oYP~!Oc3W{)?id-0mOa={qggv|LkCU`pWW?wXcs zN(xK~Jl$RYGr80cKBw3)_!pD_A4@YE@4O8`Ga&COU!nZPHFsnr9U2)Kk;l{5v;zN6 z3M{C!Ujg-IMqF!S5H#nEe5d`cqxL2+szkS^Q1z)wGXqxFbk1u6rij~V_qn=nUMW*I z!5XM6KG=AD>Gp?ylRRJ0_VJG&X