attempted formalisation and optimisation of Christopher Alexander-ish "arrange the related design elements into a bunch of (overlapping?) subsystems" problem

example_output.txt

<<< optimiser debug gibberish >>>

0000044880	ACCEPT 105 --> 104 (('grow', frozenset({1, 2, 14, 18, 26, 27}), 13))
{frozenset({33, 3, 6, 7, 8, 9, 10, 11, 19, 29}), frozenset({1, 2, 18, 26, 27, 13, 14}), frozenset({1, 2, 3, 12, 13, 20, 22, 23}), frozenset({16, 17, 21, 24, 25, 31}), frozenset({5, 15, 28, 29, 30}), frozenset({32, 17, 18, 4, 15})}

transation: on iteration 44880	i accepted the proposed move from a state with loss 105 to a new state with loss 104 via the move (('grow', frozenset({1, 2, 14, 18, 26, 27}), 13)) resulting in a new state containing the groups 

{frozenset({33, 3, 6, 7, 8, 9, 10, 11, 19, 29}), frozenset({1, 2, 18, 26, 27, 13, 14}), frozenset({1, 2, 3, 12, 13, 20, 22, 23}), frozenset({16, 17, 21, 24, 25, 31}), frozenset({5, 15, 28, 29, 30}), frozenset({32, 17, 18, 4, 15})}

<<< results >>>


*** an alleged solution ***

groups:
	'group_0'	frozenset({33, 3, 6, 7, 8, 9, 10, 11, 19, 29})
	'group_1'	frozenset({1, 2, 18, 26, 27, 13, 14})
	'group_2'	frozenset({1, 2, 3, 12, 13, 20, 22, 23})
	'group_3'	frozenset({16, 17, 21, 24, 25, 31})
	'group_4'	frozenset({5, 15, 28, 29, 30})
	'group_5'	frozenset({32, 17, 18, 4, 15})

groups containing each element:
	1	['group_2', 'group_1']
	2	['group_2', 'group_1']
	3	['group_0', 'group_2']
	4	['group_5']
	5	['group_4']
	6	['group_0']
	7	['group_0']
	8	['group_0']
	9	['group_0']
	10	['group_0']
	11	['group_0']
	12	['group_2']
	13	['group_2', 'group_1']
	14	['group_1']
	15	['group_4', 'group_5']
	16	['group_3']
	17	['group_3', 'group_5']
	18	['group_5', 'group_1']
	19	['group_0']
	20	['group_2']
	21	['group_3']
	22	['group_2']
	23	['group_2']
	24	['group_3']
	25	['group_3']
	26	['group_1']
	27	['group_1']
	28	['group_4']
	29	['group_4', 'group_0']
	30	['group_4']
	31	['group_3']
	32	['group_5']
	33	['group_0']

mcmc_optimise_alexander_design_elements.py

"""
Attempt to optimise some attempted formalisation
of Christopher Alexander "arrange the related design elements into
a bunch of (overlapping?) subsystems" problem.

featuring basic simulated-annealing-ish randomised local search
"""

import random
from itertools import combinations
import math


INFEASIBLE_PENALTY = 1000.0


def make_groups_containing(elements, groups):
    return {x:set([g for g in groups if x in g]) for x in elements}

class State:
    def __init__(self, elements, related_elements, unrelated_elements, groups):
        self.elements = elements
        self.related_elements = related_elements
        self.unrelated_elements = unrelated_elements
        self.groups = groups

    def clone(self):
        return State(
            elements=self.elements,
            related_elements=self.related_elements,
            unrelated_elements=self.unrelated_elements,
            groups=set(self.groups),
        )

    def mutate_grow(self, group, element):
        assert group in self.groups
        assert element in self.elements
        group_prime = frozenset(group | set([element]))
        state_prime = self.clone()
        state_prime.groups.remove(group)
        state_prime.groups.add(group_prime)
        return state_prime

    def mutate_prune(self, group, element):
        assert group in self.groups
        assert element in self.elements
        assert element in group
        group_prime = frozenset(group - set([element]))
        state_prime = self.clone()
        state_prime.groups.remove(group)
        state_prime.groups.add(group_prime)
        return state_prime

    def mutate_merge(self, left_group, right_group):
        assert left_group in self.groups
        assert right_group in self.groups
        merged_group = frozenset(left_group | right_group)
        state_prime = self.clone()
        state_prime.groups.remove(left_group)
        state_prime.groups.remove(right_group)
        state_prime.groups.add(merged_group)
        return state_prime

    def mutate_delete(self, group):
        assert group in self.groups
        state_prime = self.clone()
        state_prime.groups.remove(group)
        return state_prime

    def gen_moves(self):
        for g in self.groups:
            for u in self.elements:
                if u in g:
                    continue
                yield ('grow', g, u)
        for g in self.groups:
            for u in g:
                yield ('prune', g, u)
        for left_group, right_group in combinations(self.groups, 2):
            yield ('merge', left_group, right_group)
        for g in self.groups:
            yield ('delete', g)

    def apply_move(self, move):
        move_name, params = move[0], move[1:]
        return getattr(self, 'mutate_' + move_name)(*params)

    def loss(self):
        # gs[u] := set of all groups containing element u
        gs = make_groups_containing(self.elements, self.groups)
        acc = 0.0
        # penalise if u is in too few or too many groups
        for u in self.elements:
            n = len(gs[u])
            if n == 0:
                acc += INFEASIBLE_PENALTY
            elif n == 2:
                # small penalty to discourage element appearing in two groups
                acc += 1
            elif n > 2:
                acc += INFEASIBLE_PENALTY
        # penalise if we end up with empty groups.
        for g in self.groups:
            if not g:
                acc += INFEASIBLE_PENALTY
        # penalise if related elements arent in a common group
        for u, v in self.related_elements:
            if not gs[u] & gs[v]:
                acc += 1 
        # penalise each time unrelated elements share a group
        for u, v in self.unrelated_elements:
            n = len(gs[u] & gs[v])
            acc += n
        return acc


def optimise(initial_state, beta, max_iters):

    current_state = initial_state

    # Sample from distribution according to density
    # exp(-beta * loss) , i.e. weight states with
    # higher loss as having lower density.
    # c.f. https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm
    #
    # unlike physics we dont actually care about
    # correctly sampling from a distribution, we
    # want to roughly minimise loss rather than
    # sample from low-ish-loss states, but this gives
    # us a simulated annealing like optimisation
    # algorithm with a parameter beta that we can
    # tweak -- higher values of beta reduce our
    # willingness to accept moves that increase
    # the loss in the system (less exploration,
    # more optimisation of local minima)

    iters = 0

    while True:
        loss_i = current_state.loss()
        moves = list(current_state.gen_moves())
        while True:
            if iters > max_iters:
                return current_state
            move = random.choice(moves)
            proposed_state = current_state.apply_move(move)
            loss_f = proposed_state.loss()
            z = random.uniform(0.0, 1.0)
            a = min(1.0, math.exp(-beta * (loss_f - loss_i)))
            iters += 1
            if z <= a:
                print('%010d\tACCEPT %g --> %g (%r)' % (iters, loss_i, loss_f, move))
                current_state = proposed_state
                print(repr(current_state.groups))
                break


def normalise_relation(relation):
    u, v = relation
    return (u, v) if u < v else (v, u)

def make_initial_state(elements, related_elements):
    related_elements = set(normalise_relation((u, v)) for (u, v) in related_elements)
    unrelated_elements = [normalise_relation((u, v)) for (u, v) in combinations(elements, 2) if normalise_relation((u, v)) not in related_elements]

    singleton_groups = set([frozenset([u]) for u in elements])

    initial_state = State(
        elements = set(elements),
        related_elements = related_elements,
        unrelated_elements = unrelated_elements,
        groups = singleton_groups,
    )

    return initial_state


def make_dummy_example_of_system_containing_two_disjoint_subsystems():
    return make_initial_state(
        elements = set(['a', 'b', 'c', 'd']),
        related_elements=[('a', 'b'), ('c', 'd')],
    )


def make_alexander_example():
    relation_matrix = {
        1: [ 2, 3, 6, 12, 13, 14, 16, 20, 23, 25, 26, 27, 28, 33 ],
        2: [ 3, 4, 6, 10, 12, 13, 14, 19, 20, 21, 23, 25, 26, 27, 31, 32 ],
        3: [ 33, 29, 26, 23, 20, 19, 18, 17, 15, 13, 12, 11, 10, 9, 8, 7, 6 ],
        4: [ 32, 21, 20, 18, 17, 16, 11 ],
        5: [ 30, 29, 28, 25, 23, 20, 19, 15, 12, 8, 7 ],
        6: [ 33, 30, 29, 27, 25, 24, 19, 17, 15, 10, 9, 8, 7 ],
        7: [ 33, 29, 24, 23, 21, 19, 11, 10, 8 ],
        8: [ 31, 10, 9 ],
        9: [ 31, 29, 11 ],
        10: [ 29, 19, 11 ],
        11: [ 33, 24, 23, 21, 20, 16 ],
        12: [ 30, 23, 22, 20, 18, 13 ],
        13: [ 33, 28, 27, 26, 23, 22, 20, 18 ],
        14: [ 33, 29, 26, 19, 18, 15 ],
        15: [ 33, 32, 30, 29, 18, 17 ],
        16: [ 31, 27, 25, 24, 21, 20, 18, 17 ],
        17: [ 32, 31, 30, 29, 27, 25, 24, 21, 20, 19 ],
        18: [ 32, 31, 27, 26, 23, 22 ],
        19: [ 33, 29 ],
        20: [ 31, 23, 22 ],
        21: [ 31, 24, 23 ],
        22: [ 23 ],
        23: [ 33, 27 ],
        24: [ 25 ],
        25: [ ],
        26: [ 27 ],
        27: [ 30, 29 ],
        28: [ 30, 29 ],
        29: [ 30 ],
        30: [ ],
        31: [ ],
        32: [ ],
        33: [ ],
    }
    elements = set(relation_matrix.keys())
    related_elements = [normalise_relation((u, v)) for u, vvs in relation_matrix.items() for v in vvs]
    return make_initial_state(
        elements = elements,
        related_elements = related_elements,
    )

def main():
    initial_state = make_alexander_example()
    final_state = optimise(initial_state, beta=4.5, max_iters=50000)

    group_names = {g:'group_%d' % (i, ) for (i, g) in enumerate(final_state.groups)}

    print()
    print('*** an alleged solution ***')
    print()
    print('groups:')
    for g in final_state.groups:
        print('\t%r\t%r' % (group_names[g], g))

    print()
    print('groups containing each element:')
    gs = make_groups_containing(final_state.elements, final_state.groups)
    for u in final_state.elements:
        print('\t%r\t%r' % (u, [group_names[g] for g in gs[u]]))


if __name__ == '__main__':
    main()