Playing around with Miquel's theorem

README.md

Experiments involving Miquel's theorem. Joint work with Finn Klein.

miquel.py

#!/usr/bin/env python3

from typing import Tuple, FrozenSet, Set, Dict
import random
import argparse
import concurrent.futures
import queue
import os
import sys

# For testing purposes
basicConfiguration = {
    frozenset(('A', 'F', 'B', 'D')),
    frozenset(('C', 'E', 'B', 'F')),
    frozenset(('E', 'H', 'D', 'B')),
    frozenset(('H', 'G', 'A', 'D')),
    frozenset(('G', 'C', 'F', 'A')),
}

basicNumerical = {
    frozenset((1, 6, 2, 4)),
    frozenset((3, 5, 2, 6)),
    frozenset((5, 8, 4, 2)),
    frozenset((8, 7, 1, 4)),
    frozenset((7, 3, 6, 1)),
}

def choose(n : int, k : int) -> int:
    (p, q) = (1, 1)
    for i in range(k):
        p *= n-i
        q *= k-i
    return p // q

def randomConfiguration(n : int) -> FrozenSet:
    """Return a random configuration of 4 points chosen from n points."""
    return frozenset(random.sample(range(n), 4))

def orderedTuples(n : int, k : int):
    """All ordered k-tuples of numbers in range(n)."""
    if k == 0:
        yield []
    else:
        for last in range(k - 1, n):
            for t in orderedTuples(last, k-1):
                t.append(last)
                yield t

def randomConfigurations(n : int, k : int) -> Set[FrozenSet]:
    """Return k random configurations of 4 points chosen from n points."""
    cs = tuple(map(tuple, orderedTuples(n, 4)))
    return set(frozenset(t) for t in random.sample(cs, k))

def deduce(concyclic : FrozenSet[FrozenSet], verbose=False) -> Set[FrozenSet]:
    """Given a set of concyclic quadruples, deduce new ones using Miquel's theorem."""
    # h1: AFBD
    # h2: CEBF
    # h3: EHDB
    # h4: HGAD
    # h5: GCFA
    # => CEHG
    if verbose: print("Deducing from {0}".format(concyclic))
    newConfigs = set()
    for h1 in concyclic:
        for h2 in concyclic:
            bf = frozenset.intersection(h1, h2)
            # h1 and h2 must have two points in common
            if len(bf) != 2: continue
            (b, f) = tuple(bf) # maybe have to be swapped
            (a, d) = tuple(h1.difference(bf)) # maybe have to be swapped
            (c, e) = tuple(h2.difference(bf)) # maybe have to be swapped
            for h3 in concyclic:
                # h3 must contain precisely one, b or f
                if (b in h3 and f not in h3):
                    pass
                elif (b not in h3 and f in h3):
                    (b, f) = (f, b)
                else:
                    continue
                # we have b and f
                # b is in h3, f is not in h3
                if e in h3 and c not in h3:
                    pass
                elif e not in h3 and c in h3:
                    (c, e) = (e, c)
                else:
                    continue
                # we have c and e
                if c > e: continue
                # e is in h3, c is not in h3
                if d in h3 and a not in h3:
                    pass
                elif d not in h3 and a in h3:
                    (a, d) = (d, a)
                else: continue
                # we have a and d
                # d is in h3, f is not in h3
                (h,) = tuple(h3.difference((d, b, e)))
                if c > h: continue
                # we have h
                for h4 in concyclic:
                    if a in h4 and d in h4 and h in h4:
                        (g,) = tuple(h4.difference((a, d, h)))
                    else:
                        continue
                    if (c > g or e > g): continue
                    h5 = frozenset((c,f,a,g))
                    if len(h5) != 4: continue
                    h6 = frozenset((c, e, h, g))
                    if len(h6) != 4: continue
                    if h5 in concyclic and h6 not in concyclic:
                        if verbose: print("{0} => {1}".format((h1,h2,h3,h4,h5), h6))
                        newConfigs.add(h6)
    return newConfigs

def deduceAll(concyclic : FrozenSet[FrozenSet], verbose=False) -> Tuple[int, FrozenSet[FrozenSet]]:
    configs = set(concyclic)
    foundNew = True
    k = 0
    while foundNew:
        if verbose:
            print("============= Round {0} ({1} statements) =============n".format(k, len(configs)))
            print("\n".join(map(str,configs)))
        newConfigs = deduce(frozenset(configs), verbose=verbose)
        foundNew = any(c not in configs for c in newConfigs)
        configs |= newConfigs
        k += 1
    return (k, frozenset(configs))

def runSample(configurations, circles, verbose=False):
    axioms = frozenset(frozenset(t) for t in random.sample(configurations, circles))
    steps, deductiveClosure = deduceAll(axioms, verbose=verbose)
    return (len(deductiveClosure), steps)

if __name__ == '__main__':

    parser = argparse.ArgumentParser(description = "Generate images of complex zeroes")
    parser.add_argument('--samples', dest='samples', default=1, type=int, help='number of samples')
    parser.add_argument('--verbose', dest='verbose', action='store_true', default=False, help='print information')
    parser.add_argument('--points', dest='points', required=True, type=int, help='the number of points')
    parser.add_argument('--circles', dest='circles', required=True, type=int, help='the number of circles')
    args = parser.parse_args()

    print("Using {0} points and {1} circles.".format(
        args.points,
        args.circles))
    print("Creating {0} configurations.".format(choose(args.points, 4)))
    allConfigurations = tuple(orderedTuples(args.points, 4))

    samples = queue.Queue()

    num_workers = (os.cpu_count() or 2) - 1
    print("Using {0} workers.".format(num_workers))
    with concurrent.futures.ProcessPoolExecutor(max_workers=num_workers) as executor:
        futures = []
        try:
            print("Submitting tasks...")
            futures = [
                executor.submit(runSample, allConfigurations, args.circles, args.verbose)
                for _s in range(args.samples)
            ]

            # Collect results as threads complete
            print("Collecting tasks...")
            k = 0
            for future in concurrent.futures.as_completed(futures):
                print("Sampling {0}/{1} ...".format(k, args.samples), end='\r', flush=True)
                k += 1
                samples.put(future.result())
        except KeyboardInterrupt:
            # Handle Ctrl-C and shutdown the pool
            print("You pressed Ctrl-C, terminating all processes...")
            for future in futures: future.cancel()
            executor.shutdown(wait=True)
            sys.exit(1)

    # Now update the histogram with all collected results
    hist : Dict[int, Tuple[int, int]] = dict() # histogram of results

    while not samples.empty():
        outputSize, steps = samples.get()
        freq, stepsSum = hist.get(outputSize, (0,0))
        hist[outputSize] = (freq + 1, stepsSum + steps)

    print('Distribution of output sizes and average closure iterations:')
    for (outputSize, (freq, stepsSize )) in sorted(hist.items()):
        print("{0} (frequency = {1}, steps = {2:0.2f})".format(
            outputSize,
            freq,
            stepsSize/freq))