Three resistors

MCD286561.py

from collections import Counter
from math import sqrt, log

import plotly.graph_objects as go

# Rational polynomials
# x || y => 1 / (1/x + 1/y) = xy / (x+y)
fns = {
	"A": (Counter({"A": 1}), Counter({"": 1})),
	"B": (Counter({"B": 1}), Counter({"": 1})),
	"C": (Counter({"": 100, "A": -1, "B": -1}), Counter({"": 1})), # 100-A-B
	"A+B": (Counter({"A": 1, "B": 1}), Counter({"": 1})),
	"A+C": (Counter({"": 100, "B": -1}), Counter({"": 1})), # A+100-A-B
	"B+C": (Counter({"": 100, "A": -1}), Counter({"": 1})), # B+100-A-B
	"A+B+C": (Counter({"": 100}), Counter({"": 1})),
	"A||B": (Counter({"AB": 1}), Counter({"A": 1, "B": 1})), # AB/(A+B)
	"A||C": (Counter({"A": 100, "AA": -1, "AB": -1}), Counter({"": 100, "B": -1})), # A(100-A-B)/(A+100-A-B)
	"B||C": (Counter({"B": 100, "BB": -1, "AB": -1}), Counter({"": 100, "A": -1})), # B(100-A-B)/(B+100-A-B)
	"A+(B||C)": (Counter({"A": 100, "B": 100, "AA": -1, "AB": -1, "BB": -1}), Counter({"": 100, "A": -1})), # A + BC/(B+C) = (AB+AC+BC)/(B+C) = (100A+100B-AA-AB-BB))/(100-A)
	"B+(A||C)": (Counter({"A": 100, "B": 100, "AA": -1, "AB": -1, "BB": -1}), Counter({"": 100, "B": -1})), # B + AC/(A+C) = (AB+AC+BC) / (100-B)
	"C+(A||B)": (Counter({"A": 100, "B": 100, "AA": -1, "AB": -1, "BB": -1}), Counter({"A": 1, "B": 1})), # C + AB/(A+B)
	"A||(B+C)": (Counter({"A": 100, "AA": -1}), Counter({"": 100})), # A(B+C) / (A+B+C) = A(100-A) / 100
	"B||(A+C)": (Counter({"B": 100, "BB": -1}), Counter({"": 100})),
	"C||(A+B)": (Counter({"A": 100, "B": 100, "AA": -1, "AB": -2, "BB": -1}), Counter({"": 100})), # (100-A-B)(A+B) / 100 = (100A+100B-AA-2AB-BB) / 100
	"A||B||C": (Counter({"AAB": 1, "ABB": 1, "AB": -100}), Counter({"AA": 1, "AB": 1, "BB": 1, "A": -100, "B": -100})), # ABC / (AB+AC+BC) = (AB(100-A-B)) / (AB + (A+B)(100-(A+B))) = (AAB + ABB - 100AB)/(AA + AB + BB - 100A - 100B)
}

A0, B0, C0 = 200/9, 300/9, 400/9


def sym_prod(a, b):
	rv = Counter()
	for ak, av in a.items():
		for bk, bv in b.items():
			rv["".join(sorted(ak + bk))] += av * bv

	return rv


def eval_poly(p, A, B):
	rv = 0
	for k, v in p.items():
		term = v
		for ch in k:
			if ch == 'A':
				term *= A
			elif ch == 'B':
				term *= B
			else:
				assert False

		rv += term

	return rv


def eval_rat(r, A, B):
	p, q = r
	n, d = eval_poly(p, A, B), eval_poly(q, A, B)
	return n * float('+inf') if d == 0 else n / d


def build_figures(key1, key2):
	n1, d1 = fns[key1]
	n2, d2 = fns[key2]

	# Counter __sub__ Counter does unhelpful cropping at zero.
	poly = sym_prod(n1, d2)
	for k, v in sym_prod(n2, d1).items():
		poly[k] -= v

	if poly["BBB"]:
		print("Can't compare", key1, "vs", key2, "until B^3 is supported")

		fn1 = fns[key1]
		fn2 = fns[key2]

		xs = [A10 / 10 for A10 in range(350)]
		ys = [B10 / 10 for B10 in range(500)]
		zs = [[eval_rat(fn1, A, B) - eval_rat(fn2, A, B) for A in xs] for B in ys]
		return (go.Contour(z=zs, x=xs, y=ys, ncontours=40, contours_coloring='lines', name = key1 + " - " + key2), )

	cxlo = []
	cylo = []
	cxhi = []
	cyhi = []
	for A10 in range(350):
		A = A10 / 10
		coeffs = [0, 0, 0]
		for k, v in poly.items():
			term = v
			idx = 0
			for ch in k:
				if ch == 'A':
					term *= A
				elif ch == 'B':
					idx += 1
				else:
					assert False

			coeffs[idx] += term

		c, b, a = coeffs
		if a == 0:
			# bB + c = 0
			if b != 0 and 0 <= -c/b <= 50:
				cxlo.append(A)
				cylo.append(-c / b)
				cxhi.append(A)
				cyhi.append(-c / b)
		else:
			d = b**2 - 4*a*c
			if d >= 0:
				if 0 <= (-b - sqrt(d)) / (2*a) <= 50:
					cxlo.append(A)
					cylo.append((-b - sqrt(d)) / (2*a))
				if 0 <= (-b + sqrt(d)) / (2*a) <= 50:
					cxhi.append(A)
					cyhi.append((-b + sqrt(d)) / (2*a))

	return (go.Scatter(x = cxlo, y = cylo, line = {"color": "green"}, name = key1 + " == " + key2 + " (lo)"),
			go.Scatter(x = cxhi, y = cyhi, line = {"color": "blue"}, name = key1 + " == " + key2 + " (hi)"))


valid_region = go.Scatter(x = [0, 100/3,  0, 0], y = [0, 100/3, 50, 0], line = {"color": "red"}, name = "Valid region")

phase_diagram = [
	build_figures("A||(B+C)", "B||C")[0],
	build_figures("A", "B||C")[0],
	build_figures("A", "B||(A+C)")[0],
	build_figures("A", "C||(A+B)")[1],
	build_figures("B", "A+(B||C)")[1],
	build_figures("B", "C||(A+B)")[1],
	build_figures("C", "B+(A||C)")[0],
	build_figures("C", "A+(B||C)")[0],
	build_figures("C", "A+B")[1],
	build_figures("A+B", "C+(A||B)")[1]
]


step = 10**(1/17)
target = [100 / step**i for i in range(17)]
target.reverse()

def sorted_values(A, B):
	return sorted(eval_rat(r, A, B) for r in fns.values())


def goal_contour_plot(name, goal_fn, minA = 10, maxA = 30, minB = 20, maxB = 40, resolution = 200):
	As = [minA + i * (maxA - minA) / resolution for i in range(resolution)]
	Bs = [minB + i * (maxB - minB) / resolution for i in range(resolution)]
	zs = [[goal_fn([log(val) - log(tgt) for val, tgt in zip(sorted_values(A, B), target)]) for A in As] for B in Bs]

	best = min(min(row) for row in zs)
	bestAB = [(As[j], Bs[i]) for i in range(len(zs)) for j in range(len(zs[i])) if zs[i][j] == best]
	print("Goal", name, "has min score", best)
	print("Best A,B", bestAB, "=>", best)
	for A, B in bestAB[:1]:
		print()
		print("A", A)
		print("B", B)
		for i, val in enumerate(sorted_values(A, B)):
			print(i, val)

	return go.Contour(z=zs, x=As, y=Bs, ncontours=80, contours_coloring='lines', name = name)


def least_squares(minA = 10, maxA = 30, minB = 20, maxB = 40, resolution = 200):
	return goal_contour_plot("Log least squares", lambda errs: sum(err**2 for err in errs), minA, maxA, minB, maxB, resolution)


def max_error(minA = 10, maxA = 30, minB = 20, maxB = 40, resolution = 200):
	return goal_contour_plot("Max error", lambda errs: max(abs(err) for err in errs), minA, maxA, minB, maxB, resolution)


def total_error(minA = 10, maxA = 30, minB = 20, maxB = 40, resolution = 200):
	return goal_contour_plot("Total error", lambda errs: sum(abs(err) for err in errs), minA, maxA, minB, maxB, resolution)


go.Figure(phase_diagram + [least_squares(26, 28, 31.5, 34.5), valid_region]).show()
#go.Figure(phase_diagram + [max_error(24, 26.5, 27, 36), valid_region]).show()
go.Figure(phase_diagram + [total_error(26, 28.5, 32, 36), valid_region]).show()