Mondrian floorplan enumerator

mondrian_plans.sage

#!/usr/bin/sage

def partial_spans(idx, mask, blocked_rowspans):
	if not mask:
		yield tuple(), blocked_rowspans
		return

	while (mask & 1) == 0:
		idx += 1
		mask >>= 1

	r0 = idx
	rs = 1
	while mask & 1:
		if (r0,rs) not in blocked_rowspans:
			for tail, brs in partial_spans(r0 + rs, mask >> 1, blocked_rowspans | set([(r0,rs)])):
				yield ( (r0,rs), ) + tail, brs

		rs += 1
		mask >>= 1


def generate_candidates(w, h, max_rooms = None):
	# w: number of columns;
	# h: number of rows;
	# returns tuples (row0, rowspan, col0, colspan) with canonicalisation to quotient by the symmetries of the square
	if max_rooms is None:
		max_rooms = w * h

	def rot90(rooms):
		return sorted((w-c0-cs,cs,r0,rs) for (r0,rs,c0,cs) in rooms)


	def fliph(rooms):
		return sorted((r0,rs,w-c0-cs,cs) for (r0,rs,c0,cs) in rooms)


	def flipv(rooms):
		return sorted((h-r0-rs,rs,c0,cs) for (r0,rs,c0,cs) in rooms)


	def canonical(rooms):
		rooms = sorted(rooms)

		if w == h:
			equiv = [rooms]
			for _ in range(3):
				equiv.append(rot90(equiv[-1]))
			equiv += [fliph(plan) for plan in equiv]
		else:
			equiv = [rooms, fliph(rooms), flipv(rooms), fliph(flipv(rooms))]

		return tuple(min(equiv))


	def inner(open_rooms, closed_rooms, blocked_rowspans, blocked_colspans, x):
		# We must open at least one room per column for it not to be redundant.
		if len(open_rooms) + len(closed_rooms) + (w - x) > max_rooms:
			return

		if x == w:
			# Final sanity checks
			if any(room for room in open_rooms if room[3] == w):
				# Found a room which spans all columns
				return

			closed_rooms += open_rooms

			# Consecutive rows must have different column spans to avoid redundancy
			cols_by_row = [set() for _ in range(h)]
			for (r0,rs,c0,cs) in closed_rooms:
				for r in range(r0, r0 + rs):
					cols_by_row[r].add( (c0, cs) )

			for r in range(1, h):
				if cols_by_row[r-1] == cols_by_row[r]:
					return

			yield closed_rooms
			return

		# Each room is extended or not. If not, we need to open new rooms.
		# We want to keep at least one to avoid a slice, and close at least one to avoid a redundant column.
		for keeper_mask in range(1, (1 << len(open_rooms)) - 1):
			available_row_mask = int(0)
			next_open_rooms = []
			next_closed_rooms = list(closed_rooms)
			next_blocked_colspans = set(blocked_colspans)
			is_blocked = False
			for room_idx, room in enumerate(open_rooms):
				keep = (keeper_mask >> room_idx) & 1
				(r0, rs, c0, cs) = room
				if keep:
					next_open_rooms.append( (r0, rs, c0, cs+1) )
				else:
					next_closed_rooms.append(room)
					if room[2:] in next_blocked_colspans:
						is_blocked = True
					next_blocked_colspans.add(room[2:])

					available_row_mask |= ((int(1) << rs) - int(1)) << r0

			if is_blocked:
				continue

			# Generate new open rooms to cover the available row mask
			for new_rooms, next_blocked_rowspans in partial_spans(0, available_row_mask, blocked_rowspans):
				if (0,h) in new_rooms:
					continue

				yield from inner(tuple(next_open_rooms) + tuple((r0,rs,x,1) for r0, rs in new_rooms), tuple(next_closed_rooms), next_blocked_rowspans, next_blocked_colspans, x+1)


	seen_canonical = set()
	for row_spans, blocked_rowspans in partial_spans(0, (1 << h) - 1, set()):
		for plan in inner(tuple((r0,rs,0,1) for r0, rs in row_spans), tuple(), blocked_rowspans, set(), 1):
			plan = canonical(plan)
			# Canonicalise more for efficiency of manual postprocessing than efficiency of search.
			if plan not in seen_canonical:
				yield plan
				seen_canonical.add(plan)


def find_solutions(base_ring, poly_ring, I):
	# I: ideal
	def find_solutions_inner(basis, subs = None):
		# base_ring: e.q. QQ, AA, QQbar
		# basis: the basis of an ideal
		if subs is None:
			subs = dict()

		if not basis:
			yield tuple(), subs
			return

		pending_multivariate = []
		for i, eq in enumerate(basis):
			eq = poly_ring(eq).subs(subs)
			if eq == base_ring(0):
				continue

			# Is eq a non-zero constant?
			# Previously I tried `eq in base_ring`, but that can require evaluating an expensive `eq.minpoly()`.
			# It seems that when base_ring is QQ we can get eq as a constant polynomial over poly_ring, but when
			# base_ring is AA we get it as an element of base_ring, so we need two tests.
			if eq.parent() != poly_ring or not eq.variables():
				# If we get here it's a non-zero constant, and we're finding solutions from an ideal, so this is an obstruction.
				return

			if eq.is_univariate():
				# Consider solutions between 0 and 1 exclusive.
				v = eq.variable(0)
				Rz.<z> = PolynomialRing(base_ring)
				for root, _ in eq.univariate_polynomial(Rz).roots(base_ring):
					if root <= base_ring(0) or root >= base_ring(1):
						continue
					yield from find_solutions_inner(basis[i+1:] + tuple(pending_multivariate), subs | { v: root })

				return

			pending_multivariate.append(eq)

		# No univariate polynomial found.
		yield pending_multivariate, subs


	solutions = list(find_solutions_inner(tuple(I.groebner_basis())))
	if all(len(unresolved) == 0 for unresolved, resolved in solutions):
		yield from solutions
		return

	# This is generally more expensive, so we only do it if the straightforward approach fails.
	for pi in I.minimal_associated_primes():
		yield from find_solutions_inner(tuple(pi.basis))


def validate_noncongruence(x, y, plan):
	rv = dict()
	sizes = set()
	for room in plan:
		r0,rs,c0,cs = room
		x0 = sum(x[c] for c in range(c0))
		y0 = sum(y[r] for r in range(r0))
		w = sum(x[c] for c in range(c0,c0+cs))
		h = sum(y[r] for r in range(r0,r0+rs))
		rv[room] = (x0, y0, w, h)
		sizes.add((min(w, h), max(w, h)))

	return rv if len(sizes) == len(plan) else None


def analyse(plan, include_algebraic_solutions = False):
	w = max(c0+cs for (r0,rs,c0,cs) in plan)
	h = max(r0+rs for (r0,rs,c0,cs) in plan)

	R = PolynomialRing(QQ, ",".join(f"x{i}" for i in range(w)) + "," + ",".join(f"y{j}" for j in range(h)), order = "invlex")
	RA = PolynomialRing(QQbar, ",".join(f"x{i}" for i in range(w)) + "," + ",".join(f"y{j}" for j in range(h)))
	x = R.gens()[:w]
	y = R.gens()[w:]

	rect_areas = [sum(y[r] for r in range(r0,r0+rs)) * sum(x[c] for c in range(c0,c0+cs)) for (r0,rs,c0,cs) in plan]
	I = ideal([sum(x)-1, sum(y)-1] + [area - 1 / len(rect_areas) for area in rect_areas])
	for soln, subs in find_solutions(QQ, R, I):
		if len(subs) == len(R.gens()):
			layout = validate_noncongruence([subs[xi] for xi in x], [subs[yj] for yj in y], plan)
			if layout is not None:
				for room, (x0, y0, w, h) in layout.items():
					print(f"\t{room} at ({x0}, {y0}) of size {w} x {h}")
				print("=================")
				exit(0)
		else:
			print("\t", soln)
			print("", flush = True)

	if include_algebraic_solutions:
		for soln, subs in find_solutions(AA, RA, I):
			# Triggering calculations here has a big performance benefit in calculating the minpolys of the sums
			for v in subs.values():
				v.minpoly()

			if len(subs) == len(R.gens()):
				layout = validate_noncongruence([subs[RA(xi)] for xi in x], [subs[RA(yj)] for yj in y], plan)
				if layout is not None:
					plan_degree = max(max(w.minpoly().degree(), h.minpoly().degree()) for (x0,y0,w,h) in layout.values())
					print(f"{len(plan)} rooms, degree {plan_degree}")
					for room, (x0, y0, w, h) in layout.items():
						print(f"\t{room} at ({x0}, {y0}) of size {w} x {h} with minpolys {w.minpoly()} and {h.minpoly()}")
					print("", flush = True)
			else:
				print("\t", soln)
				for root in subs.values():
					print("\t\t", root, "is a root of", root.minpoly())
				print("", flush = True)


max_rooms = 12

for L in range(6, 2*(max_rooms-2) + 1):
	for w in range(3, L // 2 + 1):
		h = L - w

		if max_rooms is not None and max(w, h) > max_rooms - 2:
			continue

		print("Progress", w, h, flush = True)
		if w <= h:
			for plan in generate_candidates(w, h, max_rooms):
				analyse(plan, True)