Mondrian floorplan enumerator

mondrian_floorplans.sage

#!/usr/bin/sage

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

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

	if idx in banned:
		return

	r0 = idx
	rs = 1
	while mask & 1:
		if (r0,rs) not in blocked_rowspans:
			for tail, brs in partial_spans(r0 + rs, mask >> 1, banned, 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)
			bad_corners = set()
			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:])

					bad_corners.add(r0 + rs)
					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, bad_corners, 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(), 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, 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(k - v for k, v in subs.items()), subs
		return

	eq = basis[0].subs(subs)
	# We're finding solutions from an ideal, so constant 0 is irrelevant and any other constant is an obstruction.
	if eq.is_constant():
		if eq == base_ring(0):
			yield from find_solutions(base_ring, basis[1:], subs)

		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(base_ring, basis[1:], subs | { v: root })
	else:
		# Not seen yet in practice, but we might have a more interesting variety, in which case we propagate the ideal.
		for tail, ss in find_solutions(base_ring, basis[1:], subs):
			yield (basis[0],) + tail, ss


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)))
	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])
	printed_plan = False
	for pi in I.minimal_associated_primes():
		if pi.is_trivial():
			continue

		for soln, subs in find_solutions(QQ, pi.basis):
			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:
					if not printed_plan:
						if verbose: print(plan)
						printed_plan = True

					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(poly) for poly in pi.basis]):
				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)
		if w <= h:
			for plan in generate_candidates(w, h, max_rooms):
				analyse(plan, True)