Initial commit

.gitignore

@@ -0,0 +1,149 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# Emacs
+.dir-locals.el
+.dir-locals.el~
+
+# PyCharm
+.idea/
+
+# Mac OS
+.DS_Store
+
+# Stored Trees, Dominions, and Homomorphisms
+Tree0/
+Tree1/
+Tree2/
+Tree3/
+test2.py
+
+# Stored trees and dominions
+output/

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2025 Charlotte Aten
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,80 @@
+# Discrete neural nets
+This repository contains code and examples of implementing the notion of a
+discrete neural net and polymorphic learning in Python. For more information on
+these notions, see
+[the corresponding preprint](https://arxiv.org/abs/2308.00677) on the arXiv.
+
+### Project structure
+The scripts that define basic components of the system are in the `src` folder.
+These are:
+
+* `arithmetic_operations.py`: Definitions of arithmetic operations modulo some
+    positive integer. These are used to test the basic functionality of the
+    `NeuralNet` class.
+* `dominion.py`: Tools for creating dominions, a combinatorial object used in
+    the definition of the dominion polymorphisms in `polymorphisms.py`.
+* `dominion_setup.py`: Utilities for creating files of trees, dominions with
+    those trees as constraint graphs, and the data for the corresponding
+    polymorphisms.
+* `graphs.py`: Utilities for creating and storing simple graphs, including
+    randomly-generated trees.
+* `mnist_training_binary.py`: Describes how to manufacture binary relations
+    from the MNIST dataset which can be passed as arguments into the
+    polymorphisms in `polymorphisms.py`.
+* `neural_net.py`: Definition of the `NeuralNet` class, including feeding
+    forward and learning.
+* `operations.py`: Definitions pertaining to the `Operation` class, whose
+    objects are to be thought of as operations in the sense of universal
+    algebra/model theory.
+* `polymorphisms.py`: Definitions of polymorphisms of the Hamming graph, as
+    well as a neighbor function for the learning algorithm implemented in
+    `neural_net.py`.
+* `random_neural_net.py`: Tools for making `NeuralNet` objects with
+    randomly-chosen architectures and activation functions.
+* `relations.py`: Definitions pertaining to the `Relation` class, whose objects
+    are relations in the sense of model theory.
+
+The scripts that run various tests and example applications of the system are
+in the `tests` folder. These are:
+
+* `src.py`: This script allows horizontal imports from the sibling `src`
+    folder. (That is, it adds it to the system `PATH` variable.)
+* `test_binary_relation_polymorphisms`: Examples of the basic functionality for
+    the polymorphisms defined in `polymorphisms.py` when applied to binary
+    relations.
+* `test_dominion.py`: Examples of constructing and displaying dominions as
+    defined in `dominion.py`.
+* `test_dominion_setup.py`: Create trees and dominions for use with dominion
+    polymorphisms.
+* `test_graphs.py`: Examples of creating graphs (including random trees) as
+    defined in `graphs.py`.
+* `test_mnist_training_binary.py`: Verification that MNIST training data is
+    being loaded correctly from the training dataset.
+* `test_neural_net.py`: Examples of creating `NeuralNet`s using activation
+    functions from `arithmetic_operations.py` and the `RandomOperation` from
+    `random_neural_net.py`.
+* `test_relations.py`: Examples of the basic functionality for the `Relation`s
+    defined in `relations.py`.
+
+### Environment
+This project should run on any Python3 environment without configuration. It
+assumes that there is a project folder which contains these subdirectories:
+`src` (for source code), `tests` (for tests of basic functionality and
+examples), and `output` (for output json, image files, etc.). The `output`
+folder is in the `.gitignore`, so it should not be seen on cloning. It will be
+created when a script that needs to use it is run.
+
+### TODO
+* Reincorporate the polymorphisms for the higher-arity analogues of the
+    Hamming graph which Lillian coded.
+
+### Thanks
+Thanks to all the contributors to the original incarnation of this repository:
+* Rachel Dennis
+* Hussein Khalil
+* Lillian Stolberg
+* Kevin Xue
+* Andrey Yao
+
+Thanks also to the University of Rochester and the University of Colorado
+Boulder for supporting this project.

src/arithmetic_operations.py

@@ -0,0 +1,63 @@
+"""
+Arithmetic operations for use as neural net activation functions
+"""
+from operations import Operation
+
+
+class ModularAddition(Operation):
+    """
+    Addition modulo a positive integer.
+    """
+
+    def __init__(self, order, cache_values=False):
+        """
+        Create the addition operation modulo a given positive integer.
+
+        Arguments:
+            order (int): The modulus for performing addition.
+            cache_values (bool): Whether to memoize the operation.
+        """
+
+        # Complain if the order is nonpositive.
+        assert order > 0 
+        Operation.__init__(self, 2, lambda *x: (x[0] + x[1]) % order, 
+        cache_values)
+
+
+class ModularMultiplication(Operation):
+    """
+    Multiplication modulo a positive integer.
+    """
+
+    def __init__(self, order, cache_values=False):
+        """
+        Create the multiplication operation modulo a given positive integer.
+
+        Arguments:
+            order (int): The modulus for performing multiplication.
+            cache_values (bool): Whether to memoize the operation.
+        """
+
+        # Complain if the order is nonpositive.
+        assert order > 0 
+        Operation.__init__(self, 2, lambda *x: (x[0] * x[1]) % order, 
+        cache_values)
+
+
+class ModularNegation(Operation):
+    """
+    Negation modulo a positive integer.
+    """
+
+    def __init__(self, order, cache_values=False):
+        """
+        Create the negation operation modulo a given positive integer.
+
+        Arguments:
+            order (int): The modulus for performing negation.
+            cache_values (bool): Whether to memoize the operation.
+        """
+
+        # Complain if the order is nonpositive.
+        assert order > 0
+        Operation.__init__(self, 1, lambda *x: (-x) % order, cache_values)

src/dominion.py

@@ -0,0 +1,150 @@
+"""
+Dominion
+
+Tools for creating 2-dimensional dominions
+"""
+import random, pathlib
+from matplotlib import pyplot as plt
+import output
+
+
+class Dominion:
+    """
+    A dominion, which is a square array of entries with the property that every
+    2 by 2 subarray has at most two distinct entries. Higher-dimensional
+    analogues may be implemented in the future.
+
+    Attributes:
+        labels (frozenset): The labels which may appear as entries in the
+            dominion.
+        array (tuple of tuple): The array of entries belonging to the dominion.
+    """
+
+    def __init__(self, labels, array):
+        """
+        Create a dominion with a given array of labels.
+
+        Argument:
+            labels (iterable): The labels which may appear as entries in the
+                dominion.
+            array (iterable of iterable): The array of entries belonging to the
+                dominion.
+        """
+
+        self.labels = frozenset(labels)
+        self.array = tuple(tuple(row) for row in array)
+
+    def show(self):
+        """
+        Display a textual representation of the dominion in question.
+        """
+
+        for row in self.array:
+            print(row)
+
+    def __repr__(self):
+        return "A Dominion of size {} with {} possible labels.".format(
+        len(self.array), len(self.labels))
+
+    def __str__(self):
+        labels = '{' + ', '.join(map(str, self.labels)) + '}'
+        return "A Dominion of size {} with labels from {}.".format(
+        len(self.array), labels)
+
+    def draw(self, color_map, filename):
+        """
+        Render an image from a given dominion and color map.
+
+        Arguments:
+            color_map (string): The name of a color map.
+            filename (string): The name of the resulting file.
+        """
+
+        plt.imsave(output.path + '//{}.png'.format(filename), \
+        self.array, cmap=color_map)
+
+
+def new_row(row, labels, constraint_graph=None):
+    """
+    Construct a new row for a dominion with a given collection of labels and a
+    graph constraining which labels can appear together.
+
+    Arguments:
+        row (tuple): A tuple of labels representing a row of a dominion.
+        labels (iterable): The pixel labels used in the dominion. The entries
+            of `row` should come from this.
+        constraint_graph (Graph): The graph determining which labels can appear
+            next to each other. The vertices of `constraint_graph` should be
+            the entries of `labels`. The default value `None` behaves as though
+            the graph is the complete graph on the vertex set whose members are
+            the entries of `labels'.
+    Returns:
+        tuple: A new row which is permitted to follow `row` in a dominion with
+            the given labels and constraints.
+    """
+
+    partial_row = []
+    n = len(row)
+    for i in range(n):
+        if i == 0:
+            left_candidates = frozenset((row[0],))
+            right_candidates = frozenset((row[0], row[1]))
+        elif i == n - 1:
+            left_candidates = frozenset(
+            (row[n - 2], row[n - 1], partial_row[n - 2]))
+            right_candidates = frozenset((row[n - 1],))
+        else:
+            left_candidates = frozenset(
+            (row[i - 1], row[i], partial_row[i - 1]))
+            right_candidates = frozenset((row[i], row[i + 1]))
+        # If either side already has two candidates, we must choose from the
+        # intersection of the two sides.
+        candidates = left_candidates.intersection(right_candidates)
+        # Otherwise, it must be that both the left and right sides have only a
+        # single member. In this case, we may also choose an adjacent vertex on
+        # the constraint graph.
+        if len(left_candidates) == 1 and len(right_candidates) == 1:
+            if constraint_graph is None:
+                candidates = labels
+            else:
+                candidates = candidates.union(constraint_graph.neighbors(
+                tuple(candidates)[0]))
+        # Add a random candidate.
+        random_candidate = random.sample(list(candidates), 1)
+        partial_row += random_candidate
+    return tuple(partial_row)
+
+
+def random_dominion(size, labels, constraint_graph=None):
+    """
+    Create a random dominion given a size, collection of labels, and constraint
+    graph.
+
+    Arguments:
+        size (int): The number of rows (and columns) of the dominion.
+        labels (iterable): The pixel labels used in the dominion. The entries
+            of `row` should come from this.
+        constraint_graph (Graph): The graph determining which labels can appear
+            next to each other. The vertices of `constraint_graph` should be
+            the entries of `labels`. The default value `None` behaves as though
+            the graph is the complete graph on the vertex set whose members are
+            the entries of `labels'.
+    Returns:
+        Dominion: The randomly-generated dominion.
+    """
+
+    partial_dominion = [[random.choice(labels)]]
+    for _ in range(size - 1):
+        if constraint_graph is None:
+            new_label = random.choice(labels)
+        else:
+            new_label = random.choice(
+                tuple(constraint_graph.neighbors(
+                partial_dominion[0][-1])) + (partial_dominion[0][-1],))
+        partial_dominion[0].append(new_label)
+
+    for _ in range(size - 1):
+        next_row = new_row(partial_dominion[-1], labels, constraint_graph)
+        partial_dominion.append(next_row)
+
+    return Dominion(labels, partial_dominion)

src/dominion_setup.py

@@ -0,0 +1,102 @@
+"""
+Dominion setup
+
+Create files describing trees, dominions, and corresponding polymorphisms
+"""
+import random, json
+import output
+from graphs import random_tree, load_graph_from_file
+from dominion import random_dominion
+from relations import random_relation, random_adjacent_relation, Relation
+
+
+def grow_forest(filename, num_of_trees, num_of_vertices):
+    """
+    Add a specified number of trees to a given file.
+
+    Arguments:
+        filename (str): The name of the output file.
+        num_of_trees (int): The number of trees to be created.
+        num_of_vertices (int): How many vertices each of these trees should
+            have.
+    """
+
+    for _ in range(num_of_trees):
+        T = random_tree(range(num_of_vertices))
+        T.write_to_file(filename)
+
+
+def build_dominions(tree_filename, dominion_filename, num_of_dominions,
+    dominion_size):
+    """
+    Use the trees stored in a given file as constraint graphs for creating
+    dominions. These dominions are then stored in their own file, along with a
+    note about which tree was used to create them.
+
+    Arguments:
+        tree_filename (str): The name of the file where trees are stored.
+        dominion_filename (str): The name of the output file.
+        num_of_dominions (int): The number of dominions to be created.
+        dominion_size (int): The number of rows (and columns) of the dominions.
+    """
+
+    with open(output.path + '//{}.json'.format(tree_filename), 'r') \
+    as read_file:
+        num_of_trees = sum(1 for _ in read_file)
+    for _ in range(num_of_dominions):
+        tree_number = random.randrange(num_of_trees)
+        T = load_graph_from_file(tree_filename, tree_number)
+        D = random_dominion(dominion_size, tuple(T.vertices), T)
+        with open(output.path + '//{}.json'.format(dominion_filename), 'a') \
+        as write_file:
+            json.dump((tree_number, D.array), write_file)
+            write_file.write('\n')
+
+
+def find_homomorphisms(tree_filename, homomorphism_filename, universe_size):
+    """
+    Produce a file detailing homomorphisms from a given family of trees to a
+    given Hamming graph.
+
+    Arguments:
+        tree_filename (str): The name of the file where trees are stored.
+        homomorphism_filename (str): The name of the output file.
+        universe_size (int): The number of elements in the universe of the
+            relations to be produced.
+    """
+
+    with open(output.path + '//{}.json'.format(tree_filename), 'r') \
+    as read_file:
+        num_of_trees = sum(1 for _ in read_file)
+    for tree_number in range(num_of_trees):
+        T = load_graph_from_file(tree_filename, tree_number)
+        # Choose a root of the tree and build a list of (parent, child) node
+        # pairs.
+        unexplored_vertices = list(T.vertices)
+        next_vertices_to_check = [unexplored_vertices.pop()]
+        explored_vertices = set()
+        pairs = []
+        while unexplored_vertices:
+            next_vertex = next_vertices_to_check.pop()
+            new_neighbors = frozenset(
+            T.neighbors(next_vertex)).difference(explored_vertices)
+            for neighbor in new_neighbors:
+                pairs.append((next_vertex, neighbor))
+                unexplored_vertices.remove(neighbor)
+                next_vertices_to_check.append(neighbor)
+            explored_vertices.add(next_vertex)
+        # Create a list whose entries will become the images of each label
+        # under the homomorphism. Initialize every spot to 0.
+        homomorphism_values = len(T.vertices)*[0]
+        homomorphism_values[pairs[0][0]] = random_relation(universe_size)
+        # Starting all homomorphisms at empty relation for an experiment.
+        # homomorphism_values[pairs[0][0]] = Relation([], 28, 2)
+        for (parent, child) in pairs:
+            homomorphism_values[child] = \
+            random_adjacent_relation(homomorphism_values[parent])
+        homomorphism_values = tuple(tuple(rel.tuples)
+        for rel in homomorphism_values)
+        with open(output.path + '//{}.json'.format(homomorphism_filename),
+        'a') as write_file:
+            json.dump((tree_number, homomorphism_values), write_file)
+            write_file.write('\n')

src/graphs.py

@@ -0,0 +1,130 @@
+"""
+Graphs and trees
+"""
+import itertools, json, random
+import output
+
+
+def take_other_element(p, e):
+    """
+    Given a pair {a,e} and an element e, return a.
+
+    Arguments:
+        p (frozenset): The pair of elements, one of which is meant to be `e`.
+        e (Object): The element in question.
+    """
+
+    for x in p:
+        if x != e:
+            return x
+
+
+class Graph:
+    """
+    A simple graph. That is, a set of vertices with unordered pairs of vertices
+    as edges.
+
+    Attributes:
+        vertices (frozenset): The vertices of the graph.
+        edges (frozenset of frozenset): The unordered pairs of vertices
+            constituting the edges of the graph.
+    """
+
+    def __init__(self, vertices=frozenset(), edges=frozenset()):
+        """
+        Create a graph with given vertices and edges.
+
+        Arguments:
+            vertices (iterable): The vertices of the graph.
+            edges (iterable of iterable): The unordered pairs of vertices
+                constituting the edges of the graph.
+        """
+
+        self.vertices = frozenset(vertices)
+        self.edges = frozenset(frozenset(edge) for edge in edges)
+
+    def neighbors(self, vertex):
+        """
+        Construct an iterator through the neighbors of a vertex in the graph.
+
+        Argument:
+            vertex (Object): The vertex for which we find neighbors.
+
+        Returns:
+            iterator: The neighbors of the vertex in question.
+        """
+
+        return (take_other_element(edge, vertex) \
+        for edge in self.edges if vertex in edge)
+
+    def __repr__(self):
+
+        return "A Graph with {} vertices and {} edges.".format( \
+        len(self.vertices), len(self.edges))
+
+    def __str__(self):
+
+        vertices = '{' + ', '.join(map(str, self.vertices)) + '}'
+        edges = '{' + ', '.join('{' + ', '.join(map(str, edge)) + '}' \
+        for edge in self.edges) + '}'
+        return "A Graph with vertex set {} and edge set {}.".format( \
+        vertices, edges)
+
+    def write_to_file(self, filename):
+        """
+        Write a Graph to a json file. A file with the appropriate name will be
+        created if it doesn't already exist. Note that the target directory
+        does need to exist before this method is called. The Graph will be
+        appended to the next line of the file if it already exists.
+
+        Argument:
+            filename (str): The name of the output file.
+        """
+
+        with open(output.path + '//{}.json'.format(filename), 'a') \
+        as write_file:
+            # The Graph is rendered as a pair of lists, since frozensets are
+            # not serializable in json.
+            json.dump((tuple(self.vertices), tuple(map(tuple, self.edges))), \
+            write_file)
+            write_file.write('\n')
+
+
+def load_graph_from_file(filename, graph_number):
+    """
+    Create a Graph by reading from a json file.
+
+    Attributes:
+        filename (str): The name of the json file containing the Graph.
+        graph_number (int): The line number in the file describing the
+            desired Graph.
+    """
+
+    with open(output.path + '//{}.json'.format(filename), 'r') as read_file:
+        unprocessed_graph = \
+        itertools.islice(read_file, graph_number, graph_number+1).__next__()
+    return Graph(*json.loads(unprocessed_graph))
+
+
+def random_tree(vertices):
+    """
+    Create a random tree as a Graph object.
+
+    Argument:
+        vertices (iterable): The collection of vertices in the tree.
+            Should be nonempty.
+
+    Returns:
+        Graph: The randomly-created tree.
+    """
+
+    unplaced_vertices = set(vertices)
+    root_vertex = unplaced_vertices.pop()
+    placed_vertices = [root_vertex]
+    edges = set()
+    while unplaced_vertices:
+        new_vertex = unplaced_vertices.pop()
+        old_vertex = random.choice(placed_vertices)
+        edges.add((old_vertex, new_vertex))
+        placed_vertices.append(new_vertex)
+    return Graph(vertices, edges)

src/mnist_training_binary.py

@@ -0,0 +1,193 @@
+"""
+Modified MNIST training set for binary image classification
+"""
+import json
+import pathlib
+from relations import Relation, random_relation
+from itertools import product
+
+
+def import_mnist_data(data_type):
+    """
+    Create an iterator for MNIST data. The resulting JSON files have each line
+    representing a greyscale image of a handwritten digit. Each line is a
+    dictionary whose keys are integers between 1 and 255, or the string
+    'label'. The values associated to the integer keys are lists of the
+    coordinates at which the greyscale value for the MNIST image is equal to
+    the key. For example, if the greyscale value 25 is found at coordinate
+    [2,14], then the key 25 would be associated to a list of pairs, one of
+    which is [2,14]. Any pairs not belonging to a value in the dictionary are
+    assumed to be assigned greyscale value 0. The 'label' key has an integer
+    value between 0 and 9, indicated the intended handwritten digit for the
+    corresponding image.
+
+    Argument:
+        data_type (str): Either 'train' or 'test', depending on which data one
+            would like to convert.
+
+    Yields:
+        dict: The dictionary of data specifying a greyscale image and its
+            intended handwritten digit.
+    """
+
+    with open(str(pathlib.Path.cwd().parent) + \
+    '//..//JSONforMNIST//{}_data.json'.format(data_type), 'r') as read_file:
+        for line in read_file:
+            data = json.loads(line)
+            # By default, all the integer keys in the dictionary returned from
+            # the JSON file will be converted to strings. Let's undo this.
+            cleaned_data = \
+            {int(key): data[key] for key in data if key != 'label'}
+            cleaned_data['label'] = data['label']
+            yield cleaned_data
+
+
+def greyscale_to_binary(image, cutoff=127):
+    """
+    Convert a greyscale image from the MNIST training set to a binary relation.
+
+    Arguments:
+        image (dict): A dictionary representing a greyscale image as described
+            in `import_mnist_data`.
+        cutoff (int): Any pixel coordinates in `image` which are over this
+            value will be taken to be in the relation.
+
+    Returns:
+        Relation: A binary relation on a universe of size 28 whose pairs are
+            those coordinates from `image` which are at least as large as
+            `cutoff`.
+    """
+
+    pairs = []
+    for val in range(cutoff, 256):
+        if val in image:
+            pairs += image[val]
+    return Relation(pairs, 28)
+
+
+def mnist_binary_relations(data_type, cutoff=127):
+    """
+    Create an iterator for binary relations coming from MNIST data.
+
+    Arguments:
+        data_type (str): Either 'train' or 'test', depending on which data one
+            would like to examine.
+        cutoff: Any pixel coordinates in a greyscale image which are over this
+            value will be taken to be in the corresponding relation.
+
+    Yields:
+        tuple: A binary relation corresponding to a greyscale image from an
+            MNIST dataset and its corresponding integer label.
+    """
+
+    data = import_mnist_data(data_type)
+    for dic in data:
+        yield greyscale_to_binary(dic, cutoff), dic['label']
+
+
+def build_training_data(pairs, data_type, cutoff=127):
+    """
+    Create an iterable of pairs for training or testing a discrete neural net
+    using the MNIST datasets. Either the train data or the test data from MNIST
+    may be used.
+
+    The following values provided in `pairs` will be substituted for binary
+    relations:
+        0, 1, 2, 3, 4, 5, 6, 7, 8, or 9: These `int`s will be replaced with a
+            corresponding handwritten digit from MNIST dataset.
+        'Empty', 'Full': These strings will be replaced with the empty binary
+            relation and the full binary relation, respectively.
+
+    Arguments:
+        pairs (iterable of tuple): A sequence of pairs, the first entry being a
+            tuple of inputs and the second entry being a tuple of outputs. It
+            is assumed that all the first-entry tuples have the same length,
+            which is the number of input nodes in a neural net to be
+            trained/tested on such data. Similarly, the second-entry tuples
+            are assumed to have the same length, which is the number of output
+            nodes in a neural net to be trained/tested on such data. See the
+            description above for possible values that these tuples may
+            contain.
+        data_type (str): Either 'train' or 'test', depending on which data one
+            would like to examine.
+        cutoff (int): Any pixel coordinates in a greyscale image which are over
+            this value will be taken to be in the corresponding relation.
+
+    Yields:
+        tuple: A pair whose first entry is a dictionary indicating that a tuple
+            of binary relations is to be fed into a discrete neural net as the
+            inputs `x0`, `x1`, `x2`, etc. and whose second entry is a tuple of
+            binary relations which should appear as the corresponding outputs.
+    """
+
+    # Create a dictionary for the substitutions described above. The images
+    # corresponding to the digits will be updated dynamically from the MNIST
+    # training data.
+    substitution_dic = {i: None for i in range(10)}
+    substitution_dic['Empty'] = Relation(tuple(), 28, 2)
+    substitution_dic['Full'] = Relation(product(range(28), repeat=2), 28)
+    substitution_dic['One pixel'] = Relation(((0, 0),), 28)
+    # Load the MNIST data
+    data = mnist_binary_relations(data_type, cutoff)
+    # Initialize the images corresponding to the digits.
+    for i in range(10):
+        # For each digit, we try to find a candidate image.
+        while not substitution_dic[i]:
+            # We pull the next image from MNIST.
+            new_image = next(data)
+            # If an image for that digit hasn't been found yet, regardless of
+            # whether it was the one we intended to look for, that image will
+            # be added as the one representing its digit in `substitution_dic`.
+            if not substitution_dic[new_image[1]]:
+                substitution_dic[new_image[1]] = new_image[0]
+    for pair in pairs:
+        # Update one of the digits using the next values from MNIST.
+        new_image = next(data)
+        substitution_dic[new_image[1]] = new_image[0]
+        # Choose a new random image.
+        substitution_dic['Random'] = random_relation(28)
+        yield {'x{}'.format(i): substitution_dic[pair[0][i]]
+            for i in range(len(pair[0]))}, \
+            tuple(substitution_dic[pair[1][i]] for i in range(len(pair[1])))
+
+
+def binary_mnist_zero_one(quantity_of_zeroes, data_type, \
+quantity_of_ones=None, cutoff=127):
+    """
+    Create a data set for training a discrete neural net to recognize
+    handwritten zeroes and ones. Zeroes are labeled with the empty relation and
+    ones are labeled with the full relation.
+
+    Arguments:
+        quantity_of_zeroes (int): The number of examples of handwritten zeroes
+            to show.
+        data_type (str): Either 'train' or 'test', depending on which data one
+            would like to examine.
+        quantity_of_ones (int): The number of examples of handwritten ones to
+            show.
+        cutoff (int): Any pixel coordinates in a greyscale image which are over
+            this value will be taken to be in the corresponding relation.
+
+    Returns:
+        iterable: An iterable of training data where handwritten zeroes and
+            ones are mapped to full and empty relations.
+    """
+
+    # If the number of ones to use is not specified, it is assumed to be the
+    # same as the number of zeroes.
+    if quantity_of_ones is None:
+        quantity_of_ones = quantity_of_zeroes
+    pairs = [((0,), ('Empty',)) for _ in range(quantity_of_zeroes)]
+    pairs += [((1,), ('Full',)) for _ in range(quantity_of_ones)]
+    return build_training_data(pairs, data_type, cutoff)
+
+
+def experiment_mnist_zero_one(quantity_of_zeroes, data_type, \
+quantity_of_ones=None, cutoff=127):
+    # If the number of ones to use is not specified, it is assumed to be the
+    # same as the number of zeroes.
+    if quantity_of_ones is None:
+        quantity_of_ones = quantity_of_zeroes
+    pairs = [(('Random',), (0,)) for _ in range(quantity_of_zeroes)]
+    pairs += [(('Random',), (1,)) for _ in range(quantity_of_ones)]
+    return build_training_data(pairs, data_type, cutoff)

src/neural_net.py

@@ -0,0 +1,207 @@
+"""
+Discrete neural net
+"""
+import random
+from copy import copy
+import numpy
+
+
+class Neuron:
+    """
+    A neuron in a neural net.
+
+    Attributes:
+        activation_func (Operation): The activation function of the neuron.
+        inputs (list of Neuron): The neurons which act as inputs to the neuron
+            in question.
+    """
+
+    def __init__(self, activation_func, inputs):
+        """
+        Construct a neuron for use in a neural net.
+
+        Arguments:
+            activation_func (Operation): The activation function of the neuron.
+            inputs ((tuple of str) or (list of Neuron)): The neurons which act
+                as inputs to the neuron in question.
+        """
+
+        self.activation_func = activation_func
+        self.inputs = inputs
+
+
+class Layer:
+    """
+    A layer in a neural net.
+
+    Attribute:
+        neurons ((tuple of str) or (list of Neuron)): If `neurons` is a tuple
+            of str then we take the corresponding Layer object to be an input
+            layer for a neural net, with the entries of `neurons` being
+            distinct variable names for the arguments to the neural net.
+    """
+
+    def __init__(self, neurons):
+        """
+        Construct a layer with a given collection of neurons.
+
+        Argument:
+            neurons ((tuple of str) or (list of Neuron)): If `neurons` is a
+                tuple of str then we take the corresponding Layer object to be
+                an input layer for a neural net, with the entries of `neurons`
+                being distinct variable names for the arguments to the neural
+                net.
+        """
+
+        self.neurons = neurons
+
+
+def zero_one_loss(x, y):
+    """
+    Compute the 0-1 loss for a given pair of tuples.
+    The input tuples should have the same length.
+
+    Arguments:
+        x (tuple): A tuple of outputs from feeding forward through a neural
+            net.
+        y (tuple): A tuple of target outputs from a training set.
+
+    Returns:
+        int: Either 0 (the tuples agree) or 1 (the tuples do not agree).
+    """
+
+    return 1 - (x == y)
+
+
+class NeuralNet:
+    """
+    A (discrete) neural net.
+
+    Attribute:
+        architecture (list of Layer): The layers of the neural net, starting
+            with the input layer, whose neurons should be a list of distinct
+            variable names. Later layers should consist of Neurons carrying
+            activation functions.
+    """
+
+    def __init__(self, architecture):
+        """
+        Construct a neural net with a given architecture.
+
+        Argument:
+            architecture (list of Layer): The layers of the neural net,
+                starting with the input layer, whose neurons should be a list
+                of distinct variable names. Later layers should consist of
+                Neurons carrying activation functions.
+        """
+
+        self.architecture = architecture
+
+    def feed_forward(self, x):
+        """
+        Feed the values `x` forward through the neural net.
+
+        Argument:
+            x (dict of str: object): An assignment of variable names to values.
+
+        Returns:
+            tuple: The current values of each of the output layer neurons after
+                feeding forward.
+        """
+
+        # A copy is made so as to not modify the training data.
+        current_vals = copy(x)
+        for layer in self.architecture[1:]:
+            for neuron in layer.neurons:
+                tup = tuple(current_vals[input_neuron] for
+                input_neuron in neuron.inputs)
+                current_vals[neuron] = neuron.activation_func(*tup)
+        return tuple(current_vals[neuron] for
+        neuron in self.architecture[-1].neurons)
+
+    def empirical_loss(self, training_pairs, loss_func=zero_one_loss):
+        """
+        Calculate the current empirical loss of the neural net with respect to
+        the training pairs and loss function.
+
+        Argument:
+            training_pairs (iterable): Training pairs (x,y) where x is a
+                dictionary of inputs and y is a tuple of outputs.
+            loss_func (function): The loss function to use for training. The
+                default is the 0-1 loss.
+
+        Returns:
+            numpy.float64: The empirical loss. This is a float between 0 and 1,
+                with 0 meaning our model is perfect on the training set and 1
+                being complete failure.
+        """
+
+        # Create a tuple of loss function values for each pair in our training
+        # set, then average them.
+        return numpy.average(tuple(loss_func(self.feed_forward(x), y) for
+        (x, y) in training_pairs))
+
+    def training_step(self, training_pairs, neighbor_func,
+    loss_func=zero_one_loss):
+        """
+        Perform one step of training the neural net using the given training
+        pairs, neighbor function, and loss function. At each step a random
+        non-input neuron is explored. The neighbor function tells us which
+        other activation functions we should try in place of the one already
+        present at that neuron. We use the loss function and the training pairs
+        to determine which of these alternative activation functions we should
+        use at the given neuron instead.
+
+        Arguments:
+            training_pairs (iterable): Training pairs (x,y) where x is a
+                dictionary of inputs and y is a tuple of outputs.
+            neighbor_func (function): A function which takes an Operation as
+                input and returns an iterable of Operations as output.
+            loss_func (function): The loss function to use for training. The
+                default is the 0-1 loss.
+        """
+
+        # Select a random non-input layer from the neural net.
+        layer = random.choice(self.architecture[1:])
+        # Choose a random neuron from that layer.
+        neuron = random.choice(layer.neurons)
+        # Store a list of all the adjacent operations given by the supplied
+        # neighbor function.
+        ops = list(neighbor_func(neuron.activation_func))
+        # Also keep a list of the empirical loss associated with each of the
+        # operations in `ops`.
+        emp_loss = []
+        # Try each of the operations in `ops`.
+        for neighbor_op in ops:
+            # Change the activation function of `neuron` to the current
+            # candidate under consideration.
+            neuron.activation_func = neighbor_op
+            # Add the corresponding empirical loss (the average of the loss
+            # values) to the list of empirical losses.
+            emp_loss.append(self.empirical_loss(training_pairs, loss_func))
+        # Conclude the training step by changing the activation function of
+        # `neuron` to the candidate activation function which results in the
+        # lowest empirical loss.
+        neuron.activation_func = ops[emp_loss.index(min(emp_loss))]
+
+    def train(self, training_pairs, neighbor_func, iterations,
+    loss_func=zero_one_loss, report_loss=False):
+        """
+        Train the neural net by performing the training step repeatedly.
+
+        Arguments:
+            training_pairs (iterable): Training pairs (x,y) where x is a
+                dictionary of inputs and y is a tuple of outputs.
+            neighbor_func (function): A function which takes an Operation as
+                input and returns an iterable of Operations as output.
+            loss_func (function): The loss function to use for training. The
+                default is the 0-1 loss.
+            iterations (int): The number of training steps to perform.
+            report_loss (bool): Whether to print the final empirical loss after
+                the training has concluded.
+        """
+
+        for _ in range(iterations):
+            self.training_step(training_pairs, neighbor_func, loss_func)
+        if report_loss:
+            print(self.empirical_loss(training_pairs, loss_func))

src/operations.py

@@ -0,0 +1,146 @@
+"""
+Operations for use as neural net activation functions
+"""
+
+
+class Operation:
+    """
+    A finitary operation.
+
+    Unlike `Relation`s, the objects of the `Operation` class do not have an 
+    explicit reference to their universes. This is because in applications the 
+    universe is often more structured than an initial section of the natural 
+    numbers, so storing or type-checking this is expensive in general.
+
+    Attributes:
+        arity (int): The number of arguments the operation takes. This 
+            quantity should be at least 0. A 0-ary Operation takes empty 
+            tuples as arguments. See the method __getitem__ below for more 
+            information on this.
+        func (function or constant): The function which is used to compute the 
+            output value of the Operation when applied to some inputs.
+        cache_values (bool): Whether to store already-computed values of the 
+            Operation in memory.
+        values (dict): If `cache_values` is True then this attribute will keep 
+            track of which input-output pairs have already been computed for 
+            this Operation so that they may be reused. This can be replaced by 
+            another object that can be indexed.
+    """
+
+    def __init__(self, arity, func, cache_values=True):
+        """
+        Create a finitary operation.
+
+        Arguments:
+            arity (int): The number of arguments the operation takes. This 
+                quantity should be at least 0. A 0-ary Operation takes empty
+                tuples as arguments. See the method __getitem__ below for more
+                information on this.
+            func (function): The function which is used to compute the output
+                value of the Operation when applied to some inputs. If the
+                arity is 0, pass a constant, not a function, here.
+            cache_values (bool): Whether to store already-computed values of
+                the Operation in memory.
+        """
+
+        self.arity = arity
+        self.func = func
+        self.cache_values = cache_values
+        if self.cache_values:
+            self.values = {}
+
+    def __call__(self, *tup):
+        """
+        Compute the value of the Operation on given inputs.
+
+        Argument:
+            tup (tuple of int): The tuple of inputs to plug in to the
+                Operation.
+        """
+
+        if self.arity == 0:
+            return self.func
+        if self.cache_values:
+            if tup not in self.values:
+                self.values[tup] = self.func(*tup)
+            return self.values[tup]
+        return self.func(*tup)
+
+    def __getitem__(self, ops):
+        """
+        Form the generalized composite with a collection of operations.
+        The generalized composite of an operation f of arity k with k-many
+        operations g_i of arity n is an n-ary operation f[g_1,...,g_k]
+        where we evaluate as
+        (f[g_1,...,g_k])(x_1,...,x_n)=f(g_1(x_1,...,x_n),...,g_k(x_1,...,x_n)).
+
+        Composite operations are not memoized, but if their constituent
+        operations are memoized then the composite will perform the appropriate
+        lookups when called rather than recomputing those values from scratch.
+
+        Currently, this will not work when applied to a 0-ary operation.
+
+        Arguments:
+            ops (Operation | iterable of Operation): The operations with which
+                to form the generalized composite. This should have length
+                `self.arity` and all of its entries should have the same
+                arities.
+
+        Returns:
+            Operation: The result of composing the operations in question.
+        """
+
+        assert self.arity > 0
+        # When a single operation is being passed we turn it into a list.
+        if isinstance(ops, Operation):
+            ops = [ops]
+        assert len(ops) == self.arity
+        arities = frozenset(op.arity for op in ops)
+        assert len(arities) == 1
+        new_arity = tuple(arities)[0]
+
+        def composite(*tup):
+            """
+            Evaluate the composite operation.
+
+            Arguments:
+                *tup: A tuple of arguments to the composite operation. The
+                    length of this should be the arity of the operations g_i.
+
+            Returns:
+                object: The result of applying the generalized composite
+                    operation to the arguments.
+            """
+
+            return self(*(op(*tup) for op in ops))
+
+        return Operation(new_arity, composite, cache_values=False)
+
+
+class Identity(Operation):
+    """
+    An identity operation.
+    """
+
+    def __init__(self):
+        Operation.__init__(self, 1, lambda *x: x[0], cache_values=False)
+
+
+class Projection(Operation):
+    """
+    A projection operation.
+    """
+
+    def __init__(self, arity, coordinate):
+        Operation.__init__(self, arity, lambda *x: x[coordinate],
+        cache_values=False)
+
+
+class Constant(Operation):
+    """
+    An operation whose value is `constant` for all inputs. The default arity
+    is 0, in which case the correct way to evaluate is as f[()], not f[].
+    """
+
+    def __init__(self, constant, arity=0, cache_values=False):
+        Operation.__init__(self, arity, lambda *x: constant, cache_values)

src/output.py

@@ -0,0 +1,11 @@
+"""
+Use project output folder
+"""
+import os, pathlib
+
+# The output folder will be a sibling to `src`.
+path = str(pathlib.Path.cwd().parent / 'output')
+
+# Create the output folder if it does not exist.
+if not os.path.exists(path):
+    os.makedirs(path)

src/polymorphisms.py

@@ -0,0 +1,303 @@
+"""
+Polymorphisms
+"""
+from relations import Relation
+from operations import Operation, Projection
+import random
+import numpy
+from pathlib import Path
+import json
+import itertools
+
+
+def quarter_turn(rel):
+    """
+    Rotate a binary relation by a quarter turn counterclockwise.
+
+    Args:
+        rel (Relation): The binary relation to be rotated.
+
+    Returns:
+        Relation: The same relation rotated by a quarter turn counterclockwise.
+    """
+
+    return Relation(((rel.universe_size - tup[1], tup[0]) for tup in rel),
+    rel.universe_size, rel.arity)
+
+
+class RotationAutomorphism(Operation):
+    """
+    An automorphism of the Hamming graph obtained by rotating an image.
+    """
+
+    def __init__(self, k=1):
+        """
+        Create a rotation automorphism.
+
+        Argument:
+            k (int): The number of quarter turns by which to rotate the image
+                counterclockwise.
+        """
+
+        def func(x):
+            for _ in range(k % 4):
+                x = quarter_turn(x)
+            return x
+
+        Operation.__init__(self, 1, func=func)
+
+
+class ReflectionAutomorphism(Operation):
+    """
+    An automorphism of the Hamming graph obtained by reflecting an image across
+    its vertical axis of symmetry.
+    """
+
+    def __init__(self):
+        """
+        Create a reflection automorphism.
+        """
+
+        Operation.__init__(self, 1,
+        lambda rel: Relation(((rel.universe_size - tup[0], tup[1])
+        for tup in rel), rel.universe_size, rel.arity))
+
+
+class SwappingAutomorphism(Operation):
+    """
+    An automorphism of the Hamming graph obtained by taking the componentwise
+    xor with a fixed relation.
+    """
+
+    def __init__(self, b):
+        """
+        Create a swapping automorphism for a given relation.
+
+        Argument:
+            b (Relation): The fixed relation used for swapping. This must have
+                the same universe and arity as the argument passed to the
+                automorphism.
+        """
+
+        Operation.__init__(self, 1, lambda a: a ^ b)
+
+
+class BlankingEndomorphism(Operation):
+    """
+    An endomorphism of the Hamming graph obtained by taking the intersection
+    with a fixed relation.
+    """
+
+    def __init__(self, b):
+        """
+        Create a blanking endomorphism for a relation.
+
+        Argument:
+            b (Relation): The fixed relation used for blanking pixels.
+        """
+
+        Operation.__init__(self, 1, lambda a: a & b)
+
+
+def indicator_polymorphism(tup, a, b):
+    """
+    Perform an indicator polymorphism where the output is either an empty
+    relation or a relation containing a single tuple.
+
+    Args:
+        tup (tuple of int): The single tuple in question.
+        a (iterable of Relation): A sequence of relations, thought of as inputs
+            to the polymorphism.
+        b (iterable of Relation): A sequence of Relations with which dot
+            products are to be taken, thought of as constants. This should be
+            the same length as `a`.
+
+    Returns:
+        Relation: The relation obtained by applying the indicator polymorphism.
+    """
+
+    a = tuple(a)
+    universe_size = a[0].universe_size
+    if all(rel[0].dot(rel[1]) for rel in zip(a, b)):
+        return Relation((tup,), universe_size)
+    else:
+        return Relation(tuple(), universe_size, len(tup))
+
+
+class IndicatorPolymorphism(Operation):
+    """
+    Create a polymorphism of the Hamming graph by taking dot products with
+    fixed relations.
+    """
+
+    def __init__(self, tup, b):
+        """
+        Create an indicator polymorphism where the output is either an empty
+        relation or a relation containing a single tuple.
+
+        Args:
+            tup (tuple of int): The single tuple in question.
+            b (iterable of Relation): A sequence of Relations with which dot
+                products are to be taken, thought of as constants. Should
+                contain at least one entry.
+        """
+
+        Operation.__init__(self, len(b),
+        lambda *a: indicator_polymorphism(tup, a, b))
+
+
+def dominion_polymorphism(dominion_filename, dominion_num,
+    homomorphism_filename, a, b):
+    """
+    Perform a dominion polymorphism by using the given files.
+
+    Args:
+        dominion_filename (str): The name of the file where dominions are
+            stored.
+        dominion_num (int): The number of the dominion to use.
+        homomorphism_filename (str): The name of the file where homomorphisms
+            are stored.
+        a (Relation): The first argument to the polymorphism.
+        b (Relation): The second argument to the polymorphism.
+
+    Returns:
+        Relation: The result of performing the dominion polymorphism.
+    """
+
+    with open(str(Path(__file__).parent.resolve()) +
+    '//..//output//{}.json'.format(dominion_filename), 'r') as read_file:
+        unprocessed_dominion = itertools.islice(read_file, dominion_num,
+        dominion_num + 1).__next__()
+        tree_num, dominion = json.loads(unprocessed_dominion)
+    label = dominion[len(a)][len(b)]
+    with open(str(Path(__file__).parent.resolve()) +
+    '//..//output//{}.json'.format(homomorphism_filename), 'r') as read_file:
+        for line in read_file:
+            line_tree_num, values = json.loads(line)
+            # Note that this always takes the first homomorphism in a given
+            # file with a given tree number, even if more than one is present.
+            if line_tree_num == tree_num:
+                return Relation(values[label], a.universe_size, a.arity)
+
+
+class DominionPolymorphism(Operation):
+    """
+    A dominion polymorphism.
+    """
+
+    def __init__(self, dominion_filename, dominion_num, homomorphism_filename):
+        """
+        Create a dominion polymorphism which uses a specified dominion in
+        memory.
+
+        Arguments:
+            dominion_filename (str): The name of the file where dominions are
+                stored.
+            dominion_num (int): The number of the dominion to use.
+            homomorphism_filename (str): The name of the file where
+                homomorphisms are stored.
+        """
+
+        Operation.__init__(self, 2,
+        lambda a, b: dominion_polymorphism(dominion_filename, dominion_num,
+        homomorphism_filename, a, b))
+
+
+def polymorphism_neighbor_func(op, num_of_neighbors, constant_relations,
+    dominion_filename=None, homomorphism_filename=None):
+    """
+    Find the neighbors of a given polymorphism of the Hamming graph. Currently,
+    this assumes relations are all binary. There is also an implicit assumption
+    here that dominion polymorphisms should be binary operations. This could be
+    changed as well, but likely is not necessary.
+
+    Arguments:
+        homomorphism_filename:
+        op (Operation): A Hamming graph polymorphism operation.
+        num_of_neighbors (int): The amount of possible neighbors to generate.
+        constant_relations (iterable of Relation): An iterable of relations to
+            use as constants. This is assumed to have nonzero length.
+        dominion_filename (None | str): The name of the file where dominions
+            are stored, or None. If None, dominion polymorphisms are not used.
+        homomorphism_filename (None | str): The name of the file where
+            homomorphisms are stored, or None. If None, dominion polymorphisms
+            are not used.
+
+    Yields:
+        Operation: A neighboring operation to the given one.
+    """
+
+    endomorphisms = []
+    endomorphisms += [RotationAutomorphism(k) for k in range(4)]
+    endomorphisms.append(ReflectionAutomorphism())
+    endomorphisms.append('Swapping')
+    endomorphisms.append('Blanking')
+    constant_relations = tuple(constant_relations)
+    universe_size = constant_relations[0].universe_size
+    arity = constant_relations[0].arity
+    yield op
+    for _ in range(num_of_neighbors):
+        twist = random.randint(0, 1)
+        if twist:
+            endomorphisms_to_use = (op.arity + 1)*[RotationAutomorphism(0)]
+            twist_spot = random.randint(0, op.arity - 1)
+            endomorphisms_to_use[twist_spot] = random.choice(endomorphisms)
+            for i in range(len(endomorphisms_to_use)):
+                if endomorphisms_to_use[i] == 'Blanking':
+                    endomorphisms_to_use[i] = \
+                    BlankingEndomorphism(random.choice(constant_relations))
+                if endomorphisms_to_use[i] == 'Swapping':
+                    endomorphisms_to_use[i] = \
+                    SwappingAutomorphism(random.choice(constant_relations))
+            for i in range(len(endomorphisms_to_use) - 1):
+                endomorphisms_to_use[i] = \
+                endomorphisms_to_use[i][Projection(op.arity, i)]
+            yield endomorphisms_to_use[-1][op[endomorphisms_to_use[:-1]]]
+        else:
+            if op.arity == 1:
+                random_endomorphism = random.choice(endomorphisms)
+                if random_endomorphism == 'Blanking':
+                    random_endomorphism = \
+                    BlankingEndomorphism(random.choice(constant_relations))
+                if random_endomorphism == 'Swapping':
+                    random_endomorphism = \
+                    SwappingAutomorphism(random.choice(constant_relations))
+                yield random_endomorphism
+            if op.arity >= 2:
+                # Choose a dominion polymorphism 50% of the time, if they are
+                # available.
+                if random.randint(0, 1) and dominion_filename \
+                and homomorphism_filename:
+                    with open(str(Path(__file__).parent.resolve()) + \
+                    '//..//output//{}.json'.format(dominion_filename), 'r') \
+                    as read_file:
+                        num_of_dominions = sum(1 for _ in read_file)
+                    dominion_num = random.randint(0, num_of_dominions - 1)
+                    yield DominionPolymorphism(dominion_filename, dominion_num,
+                    homomorphism_filename)
+                else:
+                    # The universe size and relation arity for the indicator
+                    # polymorphisms is read off from the `constant_relations`.
+                    yield IndicatorPolymorphism(tuple(
+                    random.randrange(universe_size) for _ in range(arity)),
+                    random.choices(constant_relations, k=op.arity))
+
+
+def hamming_loss(x, y):
+    """
+    Compute the average Hamming loss between two iterables of relations. It is
+    assumed that `x` and `y` have the same length and that corresponding pairs
+    of relations in `x` and `y` have the same arity so that their symmetric
+    difference may be taken. In practice, this is always used when all the
+    relations belonging to `x` and `y` have the same arity.
+
+    Args:
+        x (iterable of Relation): A sequence of relations.
+        y (iterable of Relation): Another sequence of relations.
+
+    Returns:
+        numpy.float64: The average size of the symmetric difference of
+            corresponding pairs of relations in `x` and `y`.
+    """
+
+    return numpy.average(tuple(len(rel0 ^ rel1) for (rel0, rel1) in zip(x, y)))

src/random_neural_net.py

@@ -0,0 +1,112 @@
+"""
+Tools for creating random neural nets
+"""
+import random
+from neural_net import Neuron, Layer, NeuralNet
+from operations import Operation
+
+
+class RandomOperation(Operation):
+    """
+    A random operation. The values of the operation on its arguments are chosen
+    randomly and lazily, but they are memoized so that the operation is
+    well-defined.
+    """
+
+    def __init__(self, order, arity):
+        """
+        Create a random operation of a given arity and order.
+
+        Arguments:
+            order (int): The size of the universe.
+            arity (int): The arity of the operation.
+        """
+
+        if arity == 0:
+            # For a nullary operation, we choose a random member of the
+            # universe to be the corresponding constant.
+            random_constant = random.randint(0, order - 1)
+            Operation.__init__(self, 0, random_constant)
+        else:
+            Operation.__init__(self, arity,
+            lambda *x: random.randint(0, order - 1))
+
+
+class RandomNeuron(Neuron):
+    """
+    A random neuron. The activation function of the neuron will be chosen from
+    a provided list of possibilities and the inputs of the neuron will be
+    chosen from a provided previous layer.
+    """
+
+    def __init__(self, basic_ops, previous_layer):
+        """
+        Create a random neuron.
+
+        Arguments:
+            basic_ops (dict of int: iterable): The keys of this dictionary are
+                arities and the values are iterables of `Operation`s of that
+                arity.
+            previous_layer (Layer): The preceding layer from which inputs are
+                taken.
+        """
+
+        activation_func = random.choice(basic_ops[random.choice(
+        tuple(basic_ops.keys()))])
+        Neuron.__init__(self, activation_func, [random.choice(
+        previous_layer.neurons) for _ in range(activation_func.arity)])
+
+
+class RandomLayer(Layer):
+    """
+    A random layer consisting of random neurons.
+    """
+
+    def __init__(self, basic_ops, previous_layer, size):
+        """
+        Create a random layer. This takes the same dictionary of basic
+        operations as the `RandomNeuron` constructor, as well as a previous
+        layer and a desired number of nodes.
+
+        Arguments:
+            basic_ops (dict of int: iterable): The keys of this dictionary are
+                arities and the values are iterables of `Operation`s of that
+                arity.
+            previous_layer (Layer): The preceding layer from which inputs are
+                taken.
+            size (int): The number of nodes to include in the random layer.
+        """
+
+        Layer.__init__(self, [RandomNeuron(basic_ops, previous_layer)
+        for _ in range(size)])
+
+
+class RandomNeuralNet(NeuralNet):
+    """
+    A neural net whose architecture and activation functions are chosen
+    randomly.
+    """
+
+    def __init__(self, basic_ops, inputs, outputs, depth, breadth):
+        """
+        Create a random neurol net with a given collection of basic activation
+        functions. The breadth and depth of the net should be specified, as
+        well as the number of inputs/outputs, but otherwise the architecture is
+        chosen randomly.
+
+        Arguments:
+            basic_ops (dict of int: iterable): The keys of this dictionary are
+                arities and the values are iterables of `Operation`s of that
+                arity.
+            inputs (iterable of str): The names of the input neurons.
+            outputs (int): The number of output neurons.
+            depth (int): The number of layers in the neural net.
+            breadth (int): The maximum number of neurons in a layer.
+        """
+
+        architecture = [Layer(inputs)]
+        for _ in range(depth - 2):
+            architecture.append(RandomLayer(basic_ops, architecture[-1],
+            random.randint(1, breadth)))
+        architecture.append(RandomLayer(basic_ops, architecture[-1], outputs))
+        NeuralNet.__init__(self, architecture)

src/relations.py

@@ -0,0 +1,511 @@
+"""
+Relations
+"""
+from itertools import product
+from functools import wraps
+import random
+
+
+def comparison(method):
+    """
+    Check a method for an appropriate comparison by using
+    `self.comparison_check` first.
+
+    Arguments:
+        method (function): The method to which to apply this check.
+
+    Returns:
+        function: The given `method` with `self.comparison_check` being
+            called first.
+    """
+
+    @wraps(method)
+    def checked_method(self, other):
+        assert self.comparison_check(other)
+        return method(self, other)
+
+    return checked_method
+
+
+class Relation:
+    """
+    A finitary relation on a finite set.
+
+    Attributes:
+        tuples (frozenset of tuple of int): The tuples belonging to the
+            relation.
+        universe_size (int): The number of elements in the universe, which is
+            assumed to consist of an initial section of the nonnegative
+            integers.
+        arity (int): The length of each tuple in the relation. Can be inferred
+            from `tuples` unless that iterable is empty.
+    """
+
+    def __init__(self, tuples, universe_size, arity=0):
+        """
+        Construct a relation from a collection of tuples.
+
+        Arguments:
+            tuples (iterable of iterable of int): The tuples belonging to the
+                relation.
+            universe_size (int): The number of elements in the universe, which
+                is assumed to consist of an initial section of the nonnegative
+                integers.
+            arity (int): The length of each tuple in the relation. Can be
+                inferred from `tuples` unless that iterable is empty.
+        """
+
+        # Create a tuple of tuples of integers from the given iterable
+        # `tuples`.
+        tuples = tuple(tuple(entry) for entry in tuples)
+        # If `tuples` is empty then we have an empty relation and cannot infer
+        # its arity from its members. If no value is provided for the arity,
+        # it defaults to 0.
+        if tuples:
+            # We assume that all entries in `tuples` have the same length.
+            # Take one of them to get the arity of the relation.
+            self._arity = len(tuples[0])
+        else:
+            self._arity = arity
+        # Cast `tuples` to a frozenset and store it as the `tuples` attribute
+        # of the relation.
+        self._tuples = frozenset(tuples)
+        # Store the size of the universe.
+        self._universe_size = universe_size
+
+    @property
+    def tuples(self):
+        return self._tuples
+
+    @property
+    def universe_size(self):
+        return self._universe_size
+
+    @property
+    def arity(self):
+        return self._arity
+
+    def __len__(self):
+        """
+        Give the number of tuples in the relation.
+
+        Returns:
+            int: The number of tuples in `self.tuples`.
+        """
+
+        return len(self.tuples)
+
+    def __str__(self):
+        """
+        Display basic information about the relation.
+
+        Returns:
+            str: Information about the universe, arity, and size of the
+                relation.
+        """
+
+        # When the universe size is large we use ellipsis rather than write
+        # out the whole universe.
+        if self.universe_size > 10:
+            universe = '{0,...,' + str(self.universe_size - 1) + '}'
+        else:
+            universe = '{' + \
+            ','.join(map(str, range(self.universe_size))) + '}'
+        # Check whether 'tuple' needs to be pluralized.
+        if len(self) == 1:
+            plural = ''
+        else:
+            plural = 's'
+        return 'A relation on {} of arity {} containing {} tuple{}'.format(
+        universe, self.arity, len(self), plural)
+
+    def __contains__(self, tup):
+        """
+        Check whether a tuple belongs to a relation.
+
+        Argument:
+            tup (tuple of int): The tuple we are checking.
+
+        Returns:
+            bool: True when `tup` belongs to `self.tuples`, False otherwise.
+        """
+
+        return tup in self.tuples
+
+    def __iter__(self):
+        """
+        Produce an iterator for the tuples in the relation.
+
+        Returns:
+            frozenset: The set of tuples in the relation.
+        """
+
+        return iter(self.tuples)
+
+    def __bool__(self):
+        """
+        Cast a relation to a boolean value.
+
+        Returns:
+            bool: True when self.tuples is nonempty, False otherwise.
+        """
+
+        return bool(self.tuples)
+
+    def show(self, special_binary_display=None):
+        """
+        Display the members of `self.tuples`.
+
+        Argument:
+            special_binary_display (str): Show a binary relation through some
+                other method than just printing pairs. The default is None,
+                which means pairs will be printed as usual. This can be set to
+                'binary_pixels' in order to print out the binary image
+                corresponding to the relation or 'sparse' to display the
+                presence of a pair as an X and the absence of a pair as a
+                space.
+        """
+
+        if special_binary_display:
+            if special_binary_display == 'binary_pixels':
+                for row in range(self.universe_size):
+                    line = ''
+                    for column in range(self.universe_size):
+                        if (row, column) in self:
+                            line += '1'
+                        else:
+                            line += '0'
+                    print(line)
+            if special_binary_display == 'sparse':
+                for row in range(self.universe_size):
+                    line = ''
+                    for column in range(self.universe_size):
+                        if (row, column) in self:
+                            line += 'X'
+                        else:
+                            line += ' '
+                    print(line)
+            if special_binary_display == 'latex_matrix':
+                output = '\\begin{bmatrix}'
+                for row in range(self.universe_size):
+                    output += '&'.join(map(str, (int((row, column) in self)
+                    for column in range(self.universe_size))))
+                    if row != self.universe_size - 1:
+                        output += '\\\\'
+                output += '\\end{bmatrix}'
+                print(output)
+        else:
+            for tup in self:
+                print(tup)
+
+    def comparison_check(self, other):
+        """
+        Determine whether another `Relation` object is of the correct type to
+        be comparable with the relation in question.
+
+        Argument:
+            other (Relation): The other relation to which to compare this
+                relation.
+
+        Returns:
+            bool: True when the two relations have the same universe and arity,
+                False otherwise.
+        """
+
+        return self.universe_size == other.universe_size and \
+        self.arity == other.arity
+
+    def __hash__(self):
+        """
+        Find the hash value for the `Relation` object.
+
+        Returns:
+            int: The hash value of the `Relation` object.
+        """
+
+        return hash((self.tuples, self.universe_size, self.arity))
+
+    @comparison
+    def __eq__(self, other):
+        """
+        Check whether the relation is equal to another relation.
+
+        Argument:
+            other (Relation): The other relation to which to compare this
+                relation.
+
+        Returns:
+            bool: True when self.tuples is equal to other.tuples and False
+                otherwise.
+        """
+
+        return self.tuples == other.tuples
+
+    @comparison
+    def __lt__(self, other):
+        """
+        Check whether the relation is properly contained in another relation.
+
+        Argument:
+            other (Relation): The other relation to which to compare this
+                relation.
+
+        Returns:
+            bool: True when self.tuples is a proper subset of other.tuples and
+                False otherwise.
+        """
+
+        return self.tuples < other.tuples
+
+    @comparison
+    def __le__(self, other):
+        """
+        Check whether the relation is contained in another relation.
+
+        Argument:
+            other (Relation): The other relation to which to compare this
+                relation.
+
+        Returns:
+            bool: True when self.tuples is a subset of other.tuples and False
+                otherwise.
+        """
+
+        return self.tuples <= other.tuples
+
+    @comparison
+    def __gt__(self, other):
+        """
+        Check whether the relation properly contains in another relation.
+
+        Argument:
+            other (Relation): The other relation to which to compare this
+                relation.
+
+        Returns:
+            bool: True when self.tuples is a proper superset of other.tuples
+                and False otherwise.
+        """
+
+        return self.tuples > other.tuples
+
+    @comparison
+    def __ge__(self, other):
+        """
+        Check whether the relation contains in another relation.
+
+        Argument:
+            other (Relation): The other relation to which to compare this
+                relation.
+
+        Returns:
+            bool: True when self.tuples is a superset of other.tuples and False
+                otherwise.
+        """
+
+        return self.tuples >= other.tuples
+
+    def __invert__(self):
+        """
+        Create the complement of a relation. That is, a tuple in the
+        appropriate Cartesian power of the universe will belong to the
+        complement if and only if it does not belong to the given relation.
+
+        Returns:
+            Relation: The relation which is dual to the given relation in the
+                above sense.
+        """
+
+        return Relation((tup for tup in product(range(self.universe_size),
+        repeat=self.arity) if tup not in self), self.universe_size, self.arity)
+
+    @comparison
+    def __sub__(self, other):
+        """
+        Take the difference of two relations. This is the same as the set
+        difference of their sets of tuples.
+
+        Argument:
+            other (Relation): The other relation to remove from this relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+                inputs which is their set difference.
+        """
+
+        return Relation(self.tuples.difference(other.tuples),
+        self.universe_size, self.arity)
+
+    @comparison
+    def __and__(self, other):
+        """
+        Take the intersection of two relations. This is the same as bitwise
+        multiplication.
+
+        Argument:
+            other (Relation): The other relation to intersect with this
+                relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+                inputs which is their intersection.
+        """
+
+        return Relation(self.tuples.intersection(other.tuples),
+        self.universe_size, self.arity)
+
+    @comparison
+    def __or__(self, other):
+        """
+        Take the union of two relations. This is the same as bitwise
+        disjunction.
+
+        Argument:
+            other (Relation): The other relation to union with this relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+                inputs which is their union.
+        """
+
+        return Relation(self.tuples.union(other.tuples),
+        self.universe_size, self.arity)
+
+    @comparison
+    def __xor__(self, other):
+        """
+        Take the symmetric difference of two relations. This is the same as
+        bitwise addition.
+
+        Argument:
+            other (Relation): The other relation to which to add this relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+            inputs which is their symmetric difference.
+        """
+
+        return Relation(self.tuples.symmetric_difference(other.tuples),
+        self.universe_size, self.arity)
+
+    @comparison
+    def __isub__(self, other):
+        """
+        Take the set difference of two relations with augmented assignment.
+
+        Argument:
+            other (Relation): The other relation to remove from this relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+                inputs which is their set difference.
+        """
+
+        return self - other
+
+    @comparison
+    def __iand__(self, other):
+        """
+        Take the set intersection of two relations with augmented assignment.
+
+        Argument:
+            other (Relation): The other relation to remove from this relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+                inputs which is their intersection.
+        """
+
+        return self & other
+
+    @comparison
+    def __ior__(self, other):
+        """
+        Take the set union of two relations with augmented assignment.
+
+        Argument:
+            other (Relation): The other relation to union with this relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+                inputs which is their union.
+        """
+
+        return self | other
+
+    @comparison
+    def __ixor__(self, other):
+        """
+        Take the symmetric difference of two relations with augmented
+            assignment.
+
+        Argument:
+            other (Relation): The other relation to which to add this relation.
+
+        Returns:
+            Relation: The relation with the same universe and arity as the
+                inputs which is their symmetric difference.
+        """
+
+        return self ^ other
+
+    @comparison
+    def dot(self, other):
+        """
+        Take the dot product of two relations modulo 2. This is the same as
+        computing the size of the intersection of the two relations modulo 2.
+
+        Argument:
+            other (Relation): The other relation with to take the dot product.
+
+        Returns:
+            int: Either 0 or 1, depending on the parity of the number of tuples
+                in `self` and `other`.
+        """
+
+        return len(self & other) % 2
+
+
+def random_relation(universe_size):
+    """
+    Create a random binary operation whose universe is an initial section of
+    the naturals.
+
+    Argument:
+        universe_size (int): The number of elements of the universe.
+
+    Returns:
+        Relation: The randomly-created binary relation.
+    """
+
+    # Choose a random number of pairs to add to our binary relation. We might
+    # get some duplicates, so this will give an upper bound for the number of
+    # pairs in the relation.
+    rel_size = random.randint(0, universe_size ** 2 - 1)
+    rel = set()
+    for _ in range(rel_size):
+        rel.add((random.randint(0, universe_size - 1),
+        random.randint(0, universe_size - 1)))
+    return Relation(rel, universe_size, arity=2)
+
+
+def random_adjacent_relation(rel):
+    """
+    Returns a random neighbor of a given binary relation in the Hamming graph
+    by "switching" at most one tuple.
+
+    Argument:
+        rel (Relation): A binary relation.
+
+    Returns:
+        Relation: An adjacent binary relation.
+    """
+
+    size = rel.universe_size
+    # Return the original relation with probability 1/(size * size).
+    if random.randint(0, size * size - 1) == 0:
+        return rel
+    # Otherwise, flip a random tuple.
+    rand_tuple = (random.randint(0, size - 1), random.randint(0, size - 1))
+    rand_rel = Relation([rand_tuple], size, 2)
+    return rel ^ rand_rel

tests/src.py

@@ -0,0 +1,9 @@
+"""
+Import project src folder
+"""
+# Note that this depends on the `src` folder being a sibling directory to the
+# current directory.
+import sys, pathlib
+
+
+sys.path.append(str(pathlib.Path.cwd().parent / 'src'))

tests/test_binary_relation_polymorphisms.py

@@ -0,0 +1,81 @@
+"""
+Tests for binary relation polymorphisms
+"""
+import src
+from polymorphisms import RotationAutomorphism, ReflectionAutomorphism, \
+SwappingAutomorphism, BlankingEndomorphism, IndicatorPolymorphism
+from mnist_training_binary import binary_mnist_zero_one
+
+# Load some binary images from the modified MNIST training set.
+training_pairs = tuple(binary_mnist_zero_one(100, 'train'))
+
+# Create a rotation automorphism.
+rot = RotationAutomorphism()
+# Choose an image to rotate.
+print('Original image')
+img = training_pairs[24][0]['x0']
+print(type(img))
+# Display the original.
+img.show('sparse')
+# Display the rotated image.
+print('Rotated image')
+rot(img, ).show('sparse')
+# We can rotate by any number of quarter turns.
+print('Rotated half a turn')
+rot2 = RotationAutomorphism(2)
+rot2(img, ).show('sparse')
+print('Rotated three quarter turns')
+rot3 = RotationAutomorphism(3)
+rot3(img, ).show('sparse')
+
+# Create a reflection automorphism.
+refl = ReflectionAutomorphism()
+# Reflect our test image.
+print('Reflected image')
+refl(img, ).show('sparse')
+# We can compose rotations and reflections.
+print('Rotated and reflected image')
+rot(refl(img, ), ).show('sparse')
+
+# Create a swapping automorphism.
+swap = SwappingAutomorphism(training_pairs[37][0]['x0'])
+# Display the image used for swapping.
+print('Image to use for swap')
+training_pairs[37][0]['x0'].show('sparse')
+# Swap an image.
+print('Swapped image')
+swap(img, ).show('sparse')
+
+# Create a blanking endomorphism.
+blank = BlankingEndomorphism(training_pairs[37][0]['x0'])
+# Display the image used for blanking.
+print('Image to use for blanking')
+training_pairs[37][0]['x0'].show('sparse')
+# Swap an image.
+print('Blanked image')
+blank(img, ).show('sparse')
+
+# Create a binary indicator polymorphism.
+ind_pol = IndicatorPolymorphism((0, 0), \
+(training_pairs[2][0]['x0'], training_pairs[51][0]['x0']))
+# Display the images used for dot products.
+print('First image for dot product')
+training_pairs[2][0]['x0'].show('sparse')
+print('Second image used for dot product')
+training_pairs[51][0]['x0'].show('sparse')
+# Display a pair of images to which to apply the polymorphism.
+img1 = training_pairs[3][0]['x0']
+img2 = training_pairs[5][0]['x0']
+print('First input image')
+img1.show('sparse')
+print('Second input image')
+img2.show('sparse')
+# Apply the polymorphism.
+print('Image obtained from polymorphism')
+ind_pol(img1, img2).show('sparse')
+# Change one of the inputs and check the new output.
+print('New first input')
+img3 = training_pairs[34][0]['x0']
+img3.show('sparse')
+print('New image obtained from polymorphism')
+ind_pol(img3, img2).show('sparse')

tests/test_dominion.py

@@ -0,0 +1,38 @@
+"""
+Tests for creating and storing dominions
+"""
+import src
+from dominion import random_dominion
+from graphs import random_tree
+
+
+print('Create a random dominion of size 4 with labels {0,1,2}.')
+D = random_dominion(4, range(3))
+print('Show some information about this dominion.')
+print(D)
+print('Display the dominion we created.')
+D.show()
+print()
+
+print('We can also demand that only certain labels can appear next to each \
+other.')
+print('We will use a tree to illustrate this, although we could have used any \
+kind of graph we like.')
+print('The vertices of our tree will be {0,1,2,3,4,5}.')
+T = random_tree(range(6))
+print(T)
+print('Create a new random dominion of size 10 with labels {0,1,2,3,4,5} and \
+display it.')
+D = random_dominion(10, range(6), T)
+D.show()
+print()
+
+print('We can also save an image of this dominion as a file.')
+print('We\'ll use `magma` for our color map.')
+D.draw('magma', 'dominion_draw_test')
+print()
+
+print('Let\' see how long it takes to make a bigger dominion and draw it.')
+D = random_dominion(1000, range(6), T)
+D.draw('magma', 'dominion_draw_test_big')
+print()

tests/test_dominion_setup.py

@@ -0,0 +1,22 @@
+"""
+Tests for creating the files used by dominion polymorphisms
+"""
+import src
+from dominion_setup import grow_forest, build_dominions, find_homomorphisms
+
+print('Grow a forest file consisting of trees that we will use as constraint graphs.')
+print('We\'ll make 10 tress, each with 100 nodes.')
+grow_forest('forest', 10, 100)
+print()
+
+print('Create dominions with random trees from this collection.')
+# If we're going to use these dominions to create polymorphisms, then we need to have the size be one more than the
+# maximum possible Hamming weight for a given universe size for binary relations. Say our universe has 28 members.
+# The number of possible values the Hamming weight can take on is 28*28+1=785.
+print('We\'ll create 20 dominions, each of size 785. That is, universe size 28.')
+build_dominions('forest', 'dominions', 20, 785)
+print()
+
+print('Find homomorphisms from the trees we created to the Hamming graph.')
+print('We\'ll find one homomorphism for each tree we created, for the universe size 28.')
+find_homomorphisms('forest', 'homomorphisms', 28)

tests/test_graphs.py

@@ -0,0 +1,25 @@
+"""
+Graphs test
+"""
+import src
+from graphs import Graph, load_graph_from_file, random_tree
+
+print('Create a graph on the vertex set {a,b,c} with edges {a,b} and {c,b}.')
+G = Graph(('a', 'b', 'c'), (('a', 'b'), ('c', 'b')))
+print('Display the neighbors of the vertex b.')
+print(tuple(G.neighbors('b')))
+print('Display information about our graph.')
+print(G)
+print()
+
+print('Write our graph to a file.')
+print('The file will be appended to if it already exists.')
+G.write_to_file('graphs')
+print('We can read the data of this graph back from the file.')
+print(load_graph_from_file('graphs', 0))
+print()
+
+print('Create a random tree with 10 vertices.')
+T = random_tree(range(10))
+print('Display information about our tree.')
+print(T)

tests/test_mnist_training_binary.py

@@ -0,0 +1,36 @@
+"""
+Check that MNIST training/test data is functioning
+"""
+import src
+import mnist_training_binary
+
+# Create a list of 1000 training pairs.
+mnist_relations_train = mnist_training_binary.mnist_binary_relations('train')
+training_pairs = tuple(next(mnist_relations_train) for _ in range(1000))
+# Display the 59th image.
+training_pairs[59][0].show('sparse')
+# Display the corresponding label. Can you see the digit in the above array?
+print(training_pairs[59][1])
+print()
+
+# Create a list of 1000 test pairs.
+mnist_relations_test = mnist_training_binary.mnist_binary_relations('test')
+test_pairs = tuple(next(mnist_relations_test) for _ in range(1000))
+# Display the 519th image.
+test_pairs[519][0].show('sparse')
+# Display the corresponding label. Can you see the digit in the above array?
+print(test_pairs[519][1])
+print()
+
+# Create a list of 100 training pairs for use with a discrete neural net.
+zero_training_pairs = \
+tuple(mnist_training_binary.binary_mnist_zero_one(100, 'train'))
+# This digit 0 is labeled with an all-black image (all ones) to indicate it is
+# a handwritten 0.
+zero_training_pairs[0][0]['x0'].show('sparse')
+zero_training_pairs[0][1][0].show('binary_pixels')
+print()
+# This digit 1 is labeled with an all-white image (all zeroes) to indicate it
+# is not a handwritten 0.
+zero_training_pairs[100][0]['x0'].show('sparse')
+zero_training_pairs[100][1][0].show('binary_pixels')

tests/test_neural_net.py

@@ -0,0 +1,70 @@
+"""
+Discrete neural net test
+"""
+import src
+from neural_net import Neuron, Layer, NeuralNet
+import arithmetic_operations
+from random_neural_net import RandomOperation
+from itertools import product
+
+# Our neural net will have three inputs.
+layer0 = Layer(('x0', 'x1', 'x2'))
+
+# The neural net will have input and output values which are integers modulo
+# `order`.
+order = 100
+
+# The first layer has two neurons, which are initialized to carry modular
+# addition and a random operation as activation functions.
+neuron0 = Neuron(arithmetic_operations.ModularAddition(order), ('x0', 'x1'))
+neuron1 = Neuron(RandomOperation(order, 2), ('x1', 'x2'))
+layer1 = Layer([neuron0, neuron1])
+
+# The third layer has a single neuron, which is initialized to carry modular
+# multiplication.
+neuron2 = Neuron(arithmetic_operations.ModularMultiplication(5),
+[neuron0, neuron1])
+layer2 = Layer([neuron2])
+
+net = NeuralNet([layer0, layer1, layer2])
+
+# We can feed values forward and display the result.
+print(net.feed_forward({'x0': 0, 'x1': 1, 'x2': 2}))
+print()
+
+# We create a training set in an effort to teach our net how to compute
+# (x0+x1)*(x1+x2).
+# We'll do this modulo `order`.
+training_pairs = [({'x0': x[0], 'x1': x[1], 'x2': x[2]},
+(((x[0] + x[1]) * (x[1] + x[2])) % order,))
+for x in product(range(order // 2 + 1), repeat=3)]
+
+# We can check out empirical loss with respect to this training set.
+# Our loss function will just be the 0-1 loss.
+print(net.empirical_loss(training_pairs))
+print()
+
+
+def neighbor_func(op):
+    """
+    Report all the neighbors of any operation as being addition,
+    multiplication, or a random binary operation.
+    Our example only has binary operations for activation functions so we don't
+    need to be any more detailed than this.
+
+    Argument:
+        op (operation): The Operation whose neighbors we'd like to find.
+
+    Returns:
+        list of Operations: The neighboring Operations.
+    """
+
+    return [arithmetic_operations.ModularAddition(order),
+            arithmetic_operations.ModularMultiplication(order),
+            RandomOperation(order, 2)]
+
+
+# We can now begin training.
+# Usually it will only take a few training steps to learn to replace the random
+# operation with addition.
+net.train(training_pairs, neighbor_func, 5, report_loss=True)

tests/test_relations.py

@@ -0,0 +1,201 @@
+"""
+Relations test
+"""
+import src
+from relations import Relation, random_relation, random_adjacent_relation
+from itertools import product
+
+print('Create a binary relation on the set {0,1,2} whose members are the \
+pairs (0,0), (0,1), and (2,0).\n\
+Note that the pair (0,0) is repeated at the end of our list of pairs. \
+Such duplicates are ignored by the constructor for the `Relation` class.')
+print()
+R = Relation([[0, 0], [0, 1], [2, 0], [0, 0]], 3)
+
+print('We can display some basic information about the relation.')
+print(R)
+print()
+
+print('The relation has a frozenset of tuples.')
+print(R.tuples)
+print()
+
+print('In many ways it acts like a frozenset. It has a length, which is the \
+number of tuples in the relation.')
+print(len(R))
+print()
+
+print('There is a convenience function for printing the members of \
+`R.tuples`.')
+R.show()
+print()
+
+print('We can create another relation which has the same tuples and \
+universe.\n\
+Note that we can pass in any iterable of iterables as long as the innermost \
+iterables are pairs of integers.')
+S = Relation([(2, 0), (0, 1), [0, 0]], 3)
+print()
+
+print('Show the basic information about both relations we created.')
+print(R)
+print(S)
+print()
+
+print('We can check the tuples in each of these relations.')
+print('The tuples in `R`:')
+R.show()
+print('The tuples in `S`:')
+S.show()
+print()
+
+print('Observe that `R` and `S` are different objects as far as Python is \
+concerned.')
+print(R is S)
+print()
+
+print('You can also see this by printing their object ids.')
+print(id(R))
+print(id(S))
+print()
+
+print('However, `R` and `S` are equal to each other.')
+print(R == S)
+print()
+
+print('We can also do comparisons. We have that `R` is contained in `S`, but \
+this containment isn\'t proper.')
+print(R <= S)
+print(R < S)
+print()
+
+print('If we create a relation whose set of tuples is contained in those for \
+`R` but whose universe is a different size, we will see that these \
+comparisons must be between relations with the same universe and arity.')
+T = Relation({(0, 0), (0, 1)}, 2)
+print(T.tuples < R.tuples)
+try:
+    print(T < R)
+except AssertionError:
+    print('This will be printed because an AssertionError is thrown when \
+comparing two relations on different universes.')
+print()
+
+print('Naturally, comparisons in the reverse direction work, as well.')
+print(R >= S)
+print(R > S)
+print()
+
+print('Since `Relation` objects are hashable, we can use them as entries in \
+tuples, members of sets, or keys of dictionaries.')
+tup = (R, S, T)
+set_of_relations = {R, S, T}
+D = {R: 1, S: 'a', T: R}
+print()
+
+print('Arithmetic operations are also possible for relations.')
+print('We can create the bitwise complement of a relation.')
+U = ~R
+U.show()
+print()
+
+print('We can take the symmetric difference of two relations if they have the \
+same universe and arity.')
+X = Relation({(0, 0, 1), (0, 1, 1)}, 2)
+Y = Relation({(0, 0, 1), (1, 0, 1)}, 2)
+(X ^ Y).show()
+print()
+
+print('Similarly to the order comparisons, we will get an AssertionError if \
+we try to add two relations with different universes or arities.')
+Z = Relation({(0, 0, 1), (0, 1, 1)}, 3)
+try:
+    X ^ Z
+except AssertionError:
+    print('This will print since we have raised an AssertionError by trying \
+to take the symmetric difference of two relations with different universes.')
+print()
+
+print('Taking the set difference of two relations can be done as follows.')
+(X - Y).show()
+print()
+
+print('Taking the set intersection is done using the & operator. It is \
+bitwise multiplication.')
+(X & Y).show()
+print()
+
+print('Taking the set union is done using the | operator.')
+(X | Y).show()
+print()
+
+print('Any of these binary operations can be done with augmented assignment \
+as well.')
+print(X)
+X -= Y
+print(X)
+print()
+
+print('We can take the dot product of two relations modulo 2, which is the \
+same as the size of the intersection modulo 2.')
+val1 = X.dot(Y)
+val2 = len(X & Y) % 2
+print(val1, val2)
+print()
+
+print('We can check whether a given tuple belongs to a relation.')
+print((0, 0, 1) in Z)
+print((0, 1, 0) in Z)
+print()
+
+print('For binary relations, there are a few options for displaying the \
+relation.')
+W = Relation(((0, 0), (1, 1), (1, 2), (2, 0)), 3)
+W.show()
+print()
+W.show('binary_pixels')
+print()
+W.show('sparse')
+print()
+print('We can even produce LaTeX for a matrix.')
+W.show('latex_matrix')
+print()
+
+print('Let\'s show off a little bit.')
+m = 29
+# Create the circle of radius 0 in over Z/mZ.
+A = Relation(((i, j) for (i, j) in product(range(m), repeat=2)
+if (i ** 2 + j ** 2) % m == 0), m)
+# Create two translates of it.
+B = Relation(((i, j) for (i, j) in product(range(m), repeat=2)
+if ((i+2) ** 2 + j ** 2) % m == 0), m)
+C = Relation(((i, j) for (i, j) in product(range(m), repeat=2)
+if ((i+3) ** 2 + (j-1) ** 2) % m == 0), m)
+# Find all points which lie on the complement of the first circle and either \
+# of the two translates.
+D = ~A & (B | C)
+D.show('sparse')
+print()
+
+print('Note that relations are iterable.')
+for tup in Z:
+    print(tup)
+print(list(Z))
+print()
+
+print('Relations can also be used as boolean values.')
+if Z:
+    print('There are members of `Z.tuples`.')
+if Z ^ Z:
+    print('This won\'t be printed because `Z ^ Z` is empty.')
+print()
+
+print('We can create a random binary relation.')
+random_rel = random_relation(28)
+random_rel.show('sparse')
+print()
+
+print('We can also find a random relation that differs from the previous \
+one by one pixel.')
+new_random_rel = random_adjacent_relation(random_rel)
+new_random_rel.show('sparse')