diff --git a/R.py b/R.py
index 426941a..fcd6d4a 100644
--- a/R.py
+++ b/R.py
@@ -2,24 +2,22 @@
 Modeling the logic R
 """
 from logic import (
-    Conjunction,
-    Disjunction,
-    Implication,
-    Logic,
-    Negation,
     PropositionalVariable,
     Rule,
+    Logic,
+    Implication,
+    Conjunction,
+    Negation,
+    Disjunction,
+    Rule,
 )
-from model import Model, ModelFunction, ModelValue, satisfiable
+from model import Model, ModelFunction, ModelValue, has_vsp, satisfiable
 from generate_model import generate_model
-from vsp import has_vsp
 
 
 # ===================================================
 
-"""
-Defining the Logic of R
-"""
+# Defining the logic of R
 
 x = PropositionalVariable("x")
 y = PropositionalVariable("y")
@@ -61,13 +59,12 @@ logic_rules = implication_rules | negation_rules | conjunction_rules | disjuncti
 
 operations = {Negation, Conjunction, Disjunction, Implication}
 
-R_logic = Logic(operations, logic_rules, "R")
+R_logic = Logic(operations, logic_rules)
 
 # ===============================
 
-"""
-Example 2-Element Model of R
-"""
+# Example Model of R
+
 
 a0 = ModelValue("a0")
 a1 = ModelValue("a1")
@@ -90,14 +87,14 @@ mconjunction = ModelFunction(2, {
     (a0, a0): a0,
     (a0, a1): a0,
     (a1, a0): a0,
-    (a1, a1): a1
+    (a1, a1): a1  
 })
 
 mdisjunction = ModelFunction(2, {
     (a0, a0): a0,
     (a0, a1): a1,
     (a1, a0): a1,
-    (a1, a1): a1
+    (a1, a1): a1  
 })
 
 
@@ -106,7 +103,7 @@ designated_values = {a1}
 logical_operations = {
     mnegation, mimplication, mconjunction, mdisjunction
 }
-R_model_2 = Model(carrier_set, logical_operations, designated_values, "R2")
+R_model_2 = Model(carrier_set, logical_operations, designated_values)
 
 interpretation = {
     Negation: mnegation,
@@ -115,36 +112,24 @@ interpretation = {
     Implication: mimplication
 }
 
-print(R_model_2)
-
 
 # =================================
 
-"""
-Generate models of R of a specified size
-"""
-
-print("*" * 30)
+# Generate models of R of a given size
 
 model_size = 2
-print("Generating models of Logic", R_logic.name, "of size", model_size)
-solutions = generate_model(R_logic, model_size, print_model=False)
+solutions = generate_model(R_logic, model_size, print_model=True)
 
-print(f"Found {len(solutions)} satisfiable models")
+print(f"There are {len(solutions)} satisfiable models of element length {model_size}")
 
 for model, interpretation in solutions:
-    print(has_vsp(model, interpretation))
+    print("Has VSP?", has_vsp(model, interpretation))
 
-print("*" * 30)
+print("-" * 5)
 
 ######
 
-"""
-Showing the smallest model for R that has the
-variable sharing property.
-
-This model has 6 elements.
-"""
+# Smallest model for R that has the variable sharing property
 
 a0 = ModelValue("a0")
 a1 = ModelValue("a1")
@@ -163,7 +148,7 @@ mnegation = ModelFunction(1, {
     a2: a2,
     a1: a4,
     a0: a5
-}, "¬")
+})
 
 mimplication = ModelFunction(2, {
     (a5, a5): a5,
@@ -207,7 +192,7 @@ mimplication = ModelFunction(2, {
     (a0, a2): a5,
     (a0, a1): a5,
     (a0, a0): a5
-}, "→")
+})
 
 
 mconjunction = ModelFunction(2, {
@@ -252,7 +237,7 @@ mconjunction = ModelFunction(2, {
     (a0, a2): a0,
     (a0, a1): a0,
     (a0, a0): a0
-}, "∧")
+})
 
 mdisjunction = ModelFunction(2, {
     (a5, a5): a5,
@@ -296,12 +281,12 @@ mdisjunction = ModelFunction(2, {
     (a0, a2): a2,
     (a0, a1): a1,
     (a0, a0): a0
-}, "∨")
+})
 
 logical_operations = {
     mnegation, mimplication, mconjunction, mdisjunction
 }
-R_model_6 = Model(carrier_set, logical_operations, designated_values, "R6")
+R_model_6 = Model(carrier_set, logical_operations, designated_values)
 
 interpretation = {
     Negation: mnegation,
@@ -311,5 +296,7 @@ interpretation = {
 }
 
 print(R_model_6)
-print(f"Model {R_model_6.name} satisfies logic {R_logic.name}?", satisfiable(R_logic, R_model_6, interpretation))
-print(has_vsp(R_model_6, interpretation))
+
+print("Satisfiable", satisfiable(R_logic, R_model_6, interpretation))
+
+print("Has VSP?", has_vsp(R_model_6, interpretation))
diff --git a/README.md b/README.md
index 0c12b56..6e0ed93 100644
--- a/README.md
+++ b/README.md
@@ -4,9 +4,3 @@ This repository is mostly an experiment to help
 me better understand matrix models.
 
 You're likely better off using [arranstewart/magic](https://github.com/arranstewart/magic).
-
-We support output from magic using the ugly data format.
-
-```
-python3 parse_magic.py < UGLY_FILE_FROM_MAGIC
-```
diff --git a/generate_model.py b/generate_model.py
index 1306adc..2048b6e 100644
--- a/generate_model.py
+++ b/generate_model.py
@@ -1,15 +1,11 @@
 """
-Generate all the models for a given logic
-with a specified number of elements.
+File which generates all the models
 """
 from common import set_to_str
 from logic import Logic, Operation, Rule, get_operations_from_term
-from model import (
-    Interpretation, ModelValue, Model,
-    ModelFunction, satisfiable
-)
+from model import ModelValue, Model, satisfiable, ModelFunction, ModelOrderConstraint
 from itertools import combinations, chain, product
-from typing import Set, List, Tuple
+from typing import Set, List, Dict, Tuple
 
 def possible_designations(iterable):
     """Powerset without the empty and complete set"""
@@ -17,11 +13,6 @@ def possible_designations(iterable):
     return chain.from_iterable(combinations(s, r) for r in range(1, len(s)))
 
 def possible_functions(operation, carrier_set):
-    """
-    Create every possible input, output pair
-    for a given model function based on an
-    operation and a carrier set.
-    """
     arity = operation.arity
 
     inputs = list(product(carrier_set, repeat=arity))
@@ -35,19 +26,12 @@ def possible_functions(operation, carrier_set):
         yield ModelFunction(arity, new_function, operation.symbol)
 
 
-def only_rules_with(rules: Set[Rule], operation: Operation) -> List[Rule]:
-    """
-    Filter the list of rules in a logic to those
-    that only contain the logical operation specified.
-    """
-    result_rules: List[Rule] = []
+def only_rules_with(rules: Set[Rule], operation: Operation) -> Set[Rule]:
+    result_rules = []
     for rule in rules:
         is_valid = True
-        # Go through every term in the premises and conclusion
         for t in (rule.premises | {rule.conclusion,}):
             t_operations = get_operations_from_term(t)
-            # Make sure there's only one operation
-            # and that it matches the operation specified
             if len(t_operations) > 1:
                 is_valid = False
                 break
@@ -64,32 +48,23 @@ def only_rules_with(rules: Set[Rule], operation: Operation) -> List[Rule]:
 
 def possible_interpretations(
         logic: Logic, carrier_set: Set[ModelValue],
-        designated_values: Set[ModelValue]):
-    """
-    Consider every possible interpretation of operations
-    within the specified logic given the carrier set of
-    model values, and the set of designated values.
-    """
-    operations: List[Operation] = []
-    model_functions: List[List[ModelFunction]] = []
+        designated_values: Set[ModelValue], ordering: Set[ModelOrderConstraint]):
+    operations = []
+    model_functions = []
 
     for operation in logic.operations:
         operations.append(operation)
         candidate_functions = list(possible_functions(operation, carrier_set))
-        passed_functions: List[ModelFunction] = []
+        passed_functions = []
         """
-        Discard candidate functions that don't pass
-        the rules that only contain the given
-        logical operation.
+        Only consider functions that at least pass
+        in the rules with the operation by itself.
         """
         restricted_rules = only_rules_with(logic.rules, operation)
         if len(restricted_rules) > 0:
             small_logic = Logic({operation,}, restricted_rules)
-            # Add candidate functions whose small model
-            # and logic are satisfied given the restricted
-            # rule set.
             for f in candidate_functions:
-                small_model = Model(carrier_set, {f,}, designated_values)
+                small_model = Model(carrier_set, {f,}, designated_values, ordering)
                 interp = {operation: f}
                 if satisfiable(small_logic, small_model, interp):
                     passed_functions.append(f)
@@ -103,42 +78,45 @@ def possible_interpretations(
             )
         model_functions.append(passed_functions)
 
-    # The model_functions variables contains
-    # the candidate functions for each operation.
-
     functions_choice = product(*model_functions)
-    # Assign a function to each operation
     for functions in functions_choice:
         assert len(operations) == len(functions)
-        interpretation: Interpretation = dict()
+        interpretation = dict()
         for operation, function in zip(operations, functions):
             interpretation[operation] = function
         yield interpretation
 
-def generate_model(
-        logic: Logic, number_elements: int, num_solutions: int = -1,
-        print_model=False) -> List[Tuple[Model, Interpretation]]:
-    """
-    Generate the specified number of models that
-    satisfy a logic of a certain size.
-    """
+def generate_model(logic: Logic, number_elements: int, num_solutions: int = -1, print_model=False) -> List[Tuple[Model, Dict[Operation, ModelFunction]]]:
     assert number_elements > 0
     carrier_set = {
         ModelValue("a" + str(i)) for i in range(number_elements)
     }
 
+    ordering = set()
+
+    # a(0) is less than all other elements
+    a0 = ModelValue("a0")
+    for v in carrier_set:
+        if v != a0:
+            ordering.add(a0 < v)
+
+    # Every other element is less than a(n - 1)
+    an = ModelValue(f"a{number_elements-1}")
+    for v in carrier_set:
+        if an != v:
+            ordering.add(v < an)
+
     possible_designated_values = possible_designations(carrier_set)
 
-    solutions: List[Tuple[Model, Interpretation]] = []
-
+    solutions: List[Tuple[Model, Dict[Operation, ModelFunction]]] = []
     for designated_values in possible_designated_values:
         designated_values = set(designated_values)
         print("Considering models for designated values", set_to_str(designated_values))
+        possible_interps = possible_interpretations(logic, carrier_set, designated_values, ordering)
 
-        possible_interps = possible_interpretations(logic, carrier_set, designated_values)
         for interpretation in possible_interps:
             is_valid = True
-            model = Model(carrier_set, set(interpretation.values()), designated_values)
+            model = Model(carrier_set, set(interpretation.values()), designated_values, ordering)
             # Iteratively test possible interpretations
             # by adding one axiom at a time
             for rule in logic.rules:
@@ -149,6 +127,7 @@ def generate_model(
 
             if is_valid:
                 solutions.append((model, interpretation))
+                # satisfied_models.append(model)
                 if print_model:
                     print(model, flush=True)
 
diff --git a/logic.py b/logic.py
index 39eb95d..b38bbdd 100644
--- a/logic.py
+++ b/logic.py
@@ -1,4 +1,4 @@
-from typing import Optional, Set, Tuple
+from typing import Any, Set, Tuple
 from functools import lru_cache
 
 class Operation:
@@ -9,10 +9,10 @@ class Operation:
         def immutable(self, name, value):
             raise Exception("Operations are immutable")
         self.__setattr__ = immutable
-
+    
     def __hash__(self):
         return self.hashed_value
-
+    
     def __call__(self, *args):
         # Ensure the arity is met
         assert len(args) == self.arity
@@ -20,7 +20,7 @@ class Operation:
         for t in args:
             assert isinstance(t, Term)
         return OpTerm(self, args)
-
+        
 
 class Term:
     def __init__(self):
@@ -35,7 +35,7 @@ class PropositionalVariable(Term):
         def immutable(self, name, value):
             raise Exception("Propositional variables are immutable")
         self.__setattr__ = immutable
-
+    
     def __hash__(self):
         return self.hashed_value
 
@@ -53,10 +53,10 @@ class OpTerm(Term):
         def immutable(self, name, value):
             raise Exception("Terms are immutable")
         self.__setattr__ = immutable
-
+    
     def __hash__(self):
         return self.hashed_value
-
+    
     def __str__(self):
         arg_strs = [str(a) for a in self.arguments]
         return self.operation.symbol + "(" + ",".join(arg_strs) + ")"
@@ -72,13 +72,13 @@ class Inequation:
         self.antecedant = antecedant
         self.consequent = consequent
     def __str__(self):
-        return str(self.antecedant) + "≤" + str(self.consequent)
+        return str(self.antecedant) + "≤" + str(self.consequent) 
 
 class InequalityRule:
     def __init__(self, premises : Set[Inequation], conclusion: Inequation):
         self.premises = premises
         self.conclusion = conclusion
-
+    
     def __str__(self):
         str_premises = [str(p) for p in self.premises]
         str_premises2 = "{" + ",".join(str_premises) + "}"
@@ -88,22 +88,16 @@ class Rule:
     def __init__(self, premises : Set[Term], conclusion: Term):
         self.premises = premises
         self.conclusion = conclusion
-
+    
     def __str__(self):
         str_premises = [str(p) for p in self.premises]
         str_premises2 = "{" + ",".join(str_premises) + "}"
         return str(str_premises2) + "=>" + str(self.conclusion)
 
 class Logic:
-    def __init__(self,
-            operations: Set[Operation], rules: Set[Rule],
-            name: Optional[str] = None):
+    def __init__(self, operations: Set[Operation], rules: Set[Rule]):
         self.operations = operations
         self.rules = rules
-        self.name = str(abs(hash((
-            frozenset(operations),
-            frozenset(rules)
-        ))))[:5] if name is None else name
 
 
 def get_prop_var_from_term(t: Term) -> Set[PropositionalVariable]:
@@ -113,7 +107,7 @@ def get_prop_var_from_term(t: Term) -> Set[PropositionalVariable]:
     result: Set[PropositionalVariable] = set()
     for arg in t.arguments:
         result |= get_prop_var_from_term(arg)
-
+    
     return result
 
 @lru_cache
@@ -124,7 +118,7 @@ def get_propostional_variables(rules: Tuple[Rule]) -> Set[PropositionalVariable]
         # Get all vars in premises
         for premise in rule.premises:
             vars |= get_prop_var_from_term(premise)
-
+    
         # Get vars in conclusion
         vars |= get_prop_var_from_term(rule.conclusion)
 
@@ -133,9 +127,9 @@ def get_propostional_variables(rules: Tuple[Rule]) -> Set[PropositionalVariable]
 def get_operations_from_term(t: Term) -> Set[Operation]:
     if isinstance(t, PropositionalVariable):
         return set()
-
+    
     result: Set[Operation] = {t.operation,}
     for arg in t.arguments:
         result |= get_operations_from_term(arg)
-
+    
     return result
diff --git a/model.py b/model.py
index 46406b5..63cd62e 100644
--- a/model.py
+++ b/model.py
@@ -1,19 +1,19 @@
 """
-Matrix model semantics and satisfiability of
-a given logic.
+Defining what it means to be a model
 """
 from common import set_to_str
 from logic import (
-    get_propostional_variables, Logic,
-    Operation, PropositionalVariable, Term
+PropositionalVariable, get_propostional_variables, Logic, Term,
+Operation, Conjunction, Disjunction, Implication
 )
-from collections import defaultdict
-from functools import cached_property, lru_cache, reduce
-from itertools import chain, combinations_with_replacement, permutations, product
-from typing import Dict, List, Optional, Set, Tuple
+from typing import Set, Dict, Tuple, Optional
+from functools import lru_cache
+from itertools import combinations, chain, product, permutations
+from copy import deepcopy
 
 
-__all__ = ['ModelValue', 'ModelFunction', 'Model', 'Interpretation']
+__all__ = ['ModelValue', 'ModelFunction', 'Model']
+
 
 
 class ModelValue:
@@ -29,17 +29,20 @@ class ModelValue:
         return self.hashed_value
     def __eq__(self, other):
         return isinstance(other, ModelValue) and self.name == other.name
-    def __deepcopy__(self, _):
+    def __lt__(self, other):
+        assert isinstance(other, ModelValue)
+        return ModelOrderConstraint(self, other)
+    def __deepcopy__(self, memo):
         return ModelValue(self.name)
 
+
 class ModelFunction:
     def __init__(self, arity: int, mapping, operation_name = ""):
         self.operation_name = operation_name
         self.arity = arity
 
-        # Transform the mapping such that the
-        # key is always a tuple of model values
-        corrected_mapping: Dict[Tuple[ModelValue], ModelValue] = {}
+        # Correct input to always be a tuple
+        corrected_mapping = dict()
         for k, v in mapping.items():
             if isinstance(k, tuple):
                 assert len(k) == arity
@@ -53,21 +56,7 @@ class ModelFunction:
 
         self.mapping = corrected_mapping
 
-    @cached_property
-    def domain(self):
-        result_set: Set[ModelValue] = set()
-        for args in self.mapping.keys():
-            for v in args:
-                result_set.add(v)
-        return result_set
-
     def __str__(self):
-        if self.arity == 1:
-            return unary_function_str(self)
-        elif self.arity == 2:
-            return binary_function_str(self)
-
-        # Default return dictionary representation
         str_dict = dict()
         for k, v in self.mapping.items():
             inputstr = "(" + ", ".join(str(ki) for ki in k) + ")"
@@ -77,34 +66,19 @@ class ModelFunction:
     def __call__(self, *args):
         return self.mapping[args]
 
+    # def __eq__(self, other):
+    #     return isinstance(other, ModelFunction) and self.name == other.name and self.arity == other.arity
 
-def unary_function_str(f: ModelFunction) -> str:
-    assert isinstance(f, ModelFunction) and f.arity == 1
-    sorted_domain = sorted(f.domain, key=lambda v : v.name)
-    header_line = f" {f.operation_name} | " + " ".join((str(v) for v in sorted_domain))
-    sep_line = "-" + ("-" * len(f.operation_name)) + "-+-" +\
-         ("-" * len(sorted_domain)) +\
-         ("-" * reduce(lambda sum, v : sum + len(v.name), sorted_domain, 0))
-    data_line = (" " * (len(f.operation_name) + 2)) + "| " + " ".join((str(f.mapping[(v,)]) for v in sorted_domain))
-    return "\n".join((header_line, sep_line, data_line)) + "\n"
-
-def binary_function_str(f: ModelFunction) -> str:
-    assert isinstance(f, ModelFunction) and f.arity == 2
-    sorted_domain = sorted(f.domain, key=lambda v : v.name)
-    max_col_width = max(chain((len(v.name) for v in sorted_domain), (len(f.operation_name),)))
-    header_line = f" {f.operation_name} " +\
-         (" " * (max_col_width - len(f.operation_name))) + "| " +\
-         " ".join((str(v) for v in sorted_domain))
-    sep_line = "-" + ("-" * max_col_width) + "-+-" +\
-         ("-" * len(sorted_domain)) +\
-         ("-" * reduce(lambda sum, v : sum + len(v.name), sorted_domain, 0))
-    data_lines = ""
-    for row_v in sorted_domain:
-        data_line = f" {row_v.name} | " + " ".join((str(f.mapping[(row_v, col_v)]) for col_v in sorted_domain))
-        data_lines += data_line + "\n"
-    return "\n".join((header_line, sep_line, data_lines))
-
-Interpretation = Dict[Operation, ModelFunction]
+class ModelOrderConstraint:
+    # a < b
+    def __init__(self, a: ModelValue, b: ModelValue):
+        self.a = a
+        self.b = b
+    def __hash__(self):
+        return hash(self.a) * hash(self.b)
+    def __eq__(self, other):
+        return isinstance(other, ModelOrderConstraint) and \
+            self.a == other.a and self.b == other.b
 
 class Model:
     def __init__(
@@ -112,46 +86,32 @@ class Model:
             carrier_set: Set[ModelValue],
             logical_operations: Set[ModelFunction],
             designated_values: Set[ModelValue],
-            name: Optional[str] = None
+            ordering: Optional[Set[ModelOrderConstraint]] = None
     ):
         assert designated_values <= carrier_set
         self.carrier_set = carrier_set
         self.logical_operations = logical_operations
         self.designated_values = designated_values
-        self.name = str(abs(hash((
-            frozenset(carrier_set),
-            frozenset(logical_operations),
-            frozenset(designated_values)
-        ))))[:5] if name is None else name
+        self.ordering = ordering  if ordering is not None else set()
+        # TODO: Make sure ordering is "valid"
+        # That is: transitive, etc.
 
     def __str__(self):
-        result = ("=" * 25) + f"""
-Model Name: {self.name}
-Carrier Set: {set_to_str(self.carrier_set)}
+        result = f"""Carrier Set: {set_to_str(self.carrier_set)}
 Designated Values: {set_to_str(self.designated_values)}
 """
         for function in self.logical_operations:
             result += f"{str(function)}\n"
 
-        return result + ("=" * 25) + "\n"
+        return result
 
 
-def evaluate_term(
-        t: Term, f: Dict[PropositionalVariable, ModelValue],
-        interpretation: Dict[Operation, ModelFunction]) -> ModelValue:
-    """
-    Given a term in a logic, mapping
-    between terms and model values,
-    as well as an interpretation
-    of operations to model functions,
-    return the evaluated model value.
-    """
-
+def evaluate_term(t: Term, f: Dict[PropositionalVariable, ModelValue], interpretation: Dict[Operation, ModelFunction]) -> ModelValue:
     if isinstance(t, PropositionalVariable):
         return f[t]
 
     model_function = interpretation[t.operation]
-    model_arguments: List[ModelValue] = []
+    model_arguments = []
     for logic_arg in t.arguments:
         model_arg = evaluate_term(logic_arg, f, interpretation)
         model_arguments.append(model_arg)
@@ -161,15 +121,11 @@ def evaluate_term(
 def all_model_valuations(
         pvars: Tuple[PropositionalVariable],
         mvalues: Tuple[ModelValue]):
-    """
-    Given propositional variables and model values,
-    produce every possible mapping between the two.
-    """
 
     all_possible_values = product(mvalues, repeat=len(pvars))
 
     for valuation in all_possible_values:
-        mapping: Dict[PropositionalVariable, ModelValue] = {}
+        mapping: Dict[PropositionalVariable, ModelValue] = dict()
         assert len(pvars) == len(valuation)
         for pvar, value in zip(pvars, valuation):
             mapping[pvar] = value
@@ -181,92 +137,98 @@ def all_model_valuations_cached(
         mvalues: Tuple[ModelValue]):
     return list(all_model_valuations(pvars, mvalues))
 
-
-def satisfiable(logic: Logic, model: Model, interpretation: Dict[Operation, ModelFunction]) -> bool:
+def rule_ordering_satisfied(model: Model, interpretation: Dict[Operation, ModelFunction]) -> bool:
     """
-    Determine whether a model satisfies a logic
-    given an interpretation.
+    Currently testing whether this function helps with runtime...
     """
-    pvars = tuple(get_propostional_variables(tuple(logic.rules)))
-    mappings = all_model_valuations_cached(pvars, tuple(model.carrier_set))
+    if Conjunction in interpretation:
+        possible_inputs = ((a, b) for (a, b) in product(model.carrier_set, model.carrier_set))
+        for a, b in possible_inputs:
+            output = interpretation[Conjunction](a, b)
+            if a < b in model.ordering and output != a:
+                print("RETURNING FALSE")
+                return False
+            if b < a in model.ordering and output != b:
+                print("RETURNING FALSE")
+                return False
 
-    for mapping in mappings:
-        # Make sure that the model satisfies each of the rules
-        for rule in logic.rules:
-            # The check only applies if the premises are designated
-            premise_met = True
-            premise_ts: Set[ModelValue] = set()
-
-            for premise in rule.premises:
-                premise_t = evaluate_term(premise, mapping, interpretation)
-                # As soon as one premise is not designated,
-                # move to the next rule.
-                if premise_t not in model.designated_values:
-                    premise_met = False
-                    break
-                # If designated, keep track of the evaluated term
-                premise_ts.add(premise_t)
-
-            if not premise_met:
-                continue
-
-            # With the premises designated, make sure the consequent is designated
-            consequent_t = evaluate_term(rule.conclusion, mapping, interpretation)
-            if consequent_t not in model.designated_values:
+    if Disjunction in interpretation:
+        possible_inputs = ((a, b) for (a, b) in product(model.carrier_set, model.carrier_set))
+        for a, b in possible_inputs:
+            output = interpretation[Disjunction](a, b)
+            if a < b in model.ordering and output != b:
+                print("RETURNING FALSE")
+                return False
+            if b < a in model.ordering and output != a:
+                print("RETURNING FALSE")
                 return False
 
     return True
 
 
 
+def satisfiable(logic: Logic, model: Model, interpretation: Dict[Operation, ModelFunction]) -> bool:
+    pvars = tuple(get_propostional_variables(tuple(logic.rules)))
+    mappings = all_model_valuations_cached(pvars, tuple(model.carrier_set))
+
+    # NOTE: Does not look like rule ordering is helping for finding
+    # models of R...
+    if not rule_ordering_satisfied(model, interpretation):
+        return False
+
+    for mapping in mappings:
+        for rule in logic.rules:
+            premise_met = True
+            premise_ts = set()
+            for premise in rule.premises:
+                premise_t = evaluate_term(premise, mapping, interpretation)
+                if premise_t not in model.designated_values:
+                    premise_met = False
+                    break
+                premise_ts.add(premise_t)
+
+            if not premise_met:
+                continue
+
+            consequent_t = evaluate_term(rule.conclusion, mapping, interpretation)
+
+            if consequent_t not in model.designated_values:
+                return False
+
+    return True
+
+from itertools import combinations_with_replacement
+from collections import defaultdict
+
 def model_closure(initial_set: Set[ModelValue], mfunctions: Set[ModelFunction]):
-    """
-    Given an initial set of model values and a set of model functions,
-    compute the complete set of model values that are closed
-    under the operations.
-    """
     closure_set: Set[ModelValue] = initial_set
-    last_new: Set[ModelValue] = initial_set
-    changed: bool = True
+    last_new = initial_set
+    changed = True
 
     while changed:
         changed = False
-        new_elements: Set[ModelValue] = set()
-        old_closure: Set[ModelValue] = closure_set - last_new
+        new_elements = set()
+        old_closure = closure_set - last_new
 
         # arity -> args
         cached_args = defaultdict(list)
 
-        # Pass elements into each model function
         for mfun in mfunctions:
 
-            # If a previous function shared the same arity,
-            # we'll use the same set of computed arguments
-            # to pass into the model functions.
+            # Use cached args if this arity was looked at before
             if mfun.arity in cached_args:
                 for args in cached_args[mfun.arity]:
-                    # Compute the new elements
-                    # given the cached arguments.
                     element = mfun(*args)
                     if element not in closure_set:
                         new_elements.add(element)
-
-                # We don't need to compute the arguments
-                # thanks to the cache, so move onto the
-                # next function.
+                # Move onto next function
                 continue
 
-            # At this point, we don't have cached arguments, so we need
-            # to compute this set.
-
-            # Each argument must have at least one new element to not repeat
-            # work. We'll range over the number of new model values within our
-            # argument.
+            # Iterate over how many new elements would be within the arguments
+            # NOTE: To not repeat work, there must be at least one new element
             for num_new in range(1, mfun.arity + 1):
                 new_args = combinations_with_replacement(last_new, r=num_new)
                 old_args = combinations_with_replacement(old_closure, r=mfun.arity - num_new)
-                # Determine every possible ordering of the concatenated
-                # new and old model values.
                 for new_arg, old_arg in product(new_args, old_args):
                     for args in permutations(new_arg + old_arg):
                         cached_args[mfun.arity].append(args)
diff --git a/parse_magic.py b/parse_magic.py
index 5c7b746..4b95ba4 100644
--- a/parse_magic.py
+++ b/parse_magic.py
@@ -49,39 +49,27 @@ def parse_matrices(infile: TextIO) -> List[Tuple[Model, Dict]]:
 
                     for mimplication in results:
                         logical_operations = {
-                            mnegation, mimplication
+                            mnegation, mconjunction, mdisjunction,
+                            mimplication
                         }
-                        model = Model(carrier_set, logical_operations, designated_values, name=str(len(solutions)))
+                        model = Model(carrier_set, logical_operations, designated_values)
                         interpretation = {
                             Negation: mnegation,
+                            Conjunction: mconjunction,
+                            Disjunction: mdisjunction,
                             Implication: mimplication
                         }
-                        if mconjunction is not None:
-                            logical_operations.add(mconjunction)
-                            interpretation[Conjunction] = mconjunction
-                        if mdisjunction is not None:
-                            logical_operations.add(mdisjunction)
-                            interpretation[Disjunction] = mdisjunction
-
                         solutions.append((model, interpretation))
-                        print(f"Parsed Matrix {model.name}")
 
     return solutions
 
 def carrier_set_from_size(size: int):
-    """
-    Construct a carrier set of model values
-    based on the desired size.
-    """
     return {
         mvalue_from_index(i) for i in range(size + 1)
     }
 
 def parse_size(infile: TextIO) -> Optional[int]:
-    """
-    Parse the line representing the matrix size.
-    NOTE: Elements are represented in hexidecimal.
-    """
+    # Elements are represented in hexidecimal
     size = int(next(infile), 16)
     if size == -1:
         return None
@@ -89,9 +77,6 @@ def parse_size(infile: TextIO) -> Optional[int]:
     return size
 
 def parse_negation(infile: TextIO, size: int) -> Optional[ModelFunction]:
-    """
-    Parse the line representing the negation table.
-    """
     line = next(infile).strip()
     if line == '-1':
         return None
@@ -105,29 +90,18 @@ def parse_negation(infile: TextIO, size: int) -> Optional[ModelFunction]:
         y = parse_mvalue(j)
         mapping[(x, )] = y
 
-    return ModelFunction(1, mapping, "¬")
+    return ModelFunction(1, mapping, "Negation")
 
 
 def mvalue_from_index(i: int):
-    """
-    Given an index, return the hexidecimal
-    representation of the model value.
-    """
     return ModelValue(f"a{hex(i)[-1]}")
 
 def parse_mvalue(x: str) -> ModelValue:
-    """
-    Parse an element and return the model value.
-    """
     return mvalue_from_index(int(x, 16))
 
 def determine_cresult(size: int, ordering: Dict[ModelValue, ModelValue], a: ModelValue, b: ModelValue) -> ModelValue:
-    """
-    Determine what a ∧ b should be given the ordering table.
-    """
     for i in range(size + 1):
         c = mvalue_from_index(i)
-        
         if not ordering[(c, a)]:
             continue
         if not ordering[(c, b)]:
@@ -145,10 +119,9 @@ def determine_cresult(size: int, ordering: Dict[ModelValue, ModelValue], a: Mode
         if not invalid:
             return c
 
+    print(a, "&", b, "is not defined")
+
 def determine_dresult(size: int, ordering: Dict[ModelValue, ModelValue], a: ModelValue, b: ModelValue) -> ModelValue:
-    """
-    Determine what a ∨ b should be given the ordering table.
-    """
     for i in range(size + 1):
         c = mvalue_from_index(i)
         if not ordering[(a, c)]:
@@ -170,9 +143,6 @@ def determine_dresult(size: int, ordering: Dict[ModelValue, ModelValue], a: Mode
             return c
 
 def parse_order(infile: TextIO, size: int) -> Optional[Tuple[ModelFunction, ModelFunction]]:
-    """
-    Parse the line representing the ordering table
-    """
     line = next(infile).strip()
     if line == '-1':
         return None
@@ -200,30 +170,16 @@ def parse_order(infile: TextIO, size: int) -> Optional[Tuple[ModelFunction, Mode
         for j in range(size + 1):
             y = mvalue_from_index(j)
 
-            cresult = determine_cresult(size, omapping, x, y)
-            if cresult is None:
-                print("[Warning] Conjunction and Disjunction are not well-defined")
-                print(f"{x} ∧ {y} = ??")
-                return None, None
-            cmapping[(x, y)] = cresult
-
-            dresult = determine_dresult(size, omapping, x, y)
-            if dresult is None:
-                print("[Warning] Conjunction and Disjunction are not well-defined")
-                print(f"{x} ∨ {y} = ??")
-                return None, None
-            dmapping[(x, y)] = dresult
+            cmapping[(x, y)] = determine_cresult(size, omapping, x, y)
+            dmapping[(x, y)] = determine_dresult(size, omapping, x, y)
 
 
-    mconjunction = ModelFunction(2, cmapping, "∧")
-    mdisjunction = ModelFunction(2, dmapping, "∨")
+    mconjunction = ModelFunction(2, cmapping, "Conjunction")
+    mdisjunction = ModelFunction(2, dmapping, "Disjunction")
 
     return mconjunction, mdisjunction
 
 def parse_designated(infile: TextIO, size: int) -> Optional[Set[ModelValue]]:
-    """
-    Parse the line representing which model values are designated.
-    """
     line = next(infile).strip()
     if line == '-1':
         return None
@@ -242,10 +198,6 @@ def parse_designated(infile: TextIO, size: int) -> Optional[Set[ModelValue]]:
 
 
 def parse_implication(infile: TextIO, size: int) -> Optional[List[ModelFunction]]:
-    """
-    Parse the line representing the list of implication
-    tables.
-    """
     line = next(infile).strip()
     if line == '-1':
         return None
@@ -271,7 +223,7 @@ def parse_implication(infile: TextIO, size: int) -> Optional[List[ModelFunction]
 
                 mapping[(x, y)] = r
 
-        mimplication = ModelFunction(2, mapping, "→")
+        mimplication = ModelFunction(2, mapping, "Implication")
         mimplications.append(mimplication)
 
     return mimplications
@@ -281,5 +233,9 @@ if __name__ == "__main__":
     solutions: List[Model] = parse_matrices(sys.stdin)
     print(f"Parsed {len(solutions)} matrices")
     for i, (model, interpretation) in enumerate(solutions):
-        print(model)
-        print(has_vsp(model, interpretation))
+        # print(model)
+        if has_vsp(model, interpretation):
+            print(model)
+            print("Has VSP")
+        else:
+            print("Model", i, "does not have VSP")
diff --git a/vsp.py b/vsp.py
index 23d6beb..ff142af 100644
--- a/vsp.py
+++ b/vsp.py
@@ -3,64 +3,40 @@ Check to see if the model has the variable
 sharing property.
 """
 from itertools import chain, combinations, product
-from typing import Dict, List, Optional, Set, Tuple
-from common import set_to_str
+from typing import Dict, Set, Tuple, List
 from model import (
-    Model, model_closure, ModelFunction, ModelValue
+    Model, ModelFunction, ModelValue, model_closure
 )
 from logic import Implication, Operation
 
-def preseed(
-        initial_set: Set[ModelValue],
-        cache:List[Tuple[Set[ModelValue], Set[ModelValue]]]):
+def preseed(initial_set: Set[ModelValue], cache:List[Tuple[Set[ModelValue], Set[ModelValue]]]):
     """
-    Given a cache of previous model_closure calls,
-    use this to compute an initial model closure
-    set based on the initial set.
+    Cache contains caches of model closure calls:
+    Ex:
+    {1, 2, 3} -> {....}
 
-    Basic Idea:
-    Let {1, 2, 3} -> X be in the cache.
     If {1,2,3} is a subset of initial set,
-    then X is the subset of the output of model_closure.
+    then {....} is the subset of the output of model_closure.
 
-    This is used to speed up subsequent calls to model_closure
+    We'll use the output to speed up the saturation procedure
     """
     candidate_preseed: Tuple[Set[ModelValue], int] = (None, None)
 
     for i, o in cache:
         if i < initial_set:
             cost = len(initial_set - i)
-            # If i is a subset with less missing elements than
-            # the previous candidate, then it's the new candidate.
             if candidate_preseed[1] is None or cost < candidate_preseed[1]:
                 candidate_preseed = o, cost
 
     same_set = candidate_preseed[1] == 0
     return candidate_preseed[0], same_set
 
-class VSP_Result:
-    def __init__(
-            self, has_vsp: bool, model_name: Optional[str] = None,
-            subalgebra1: Optional[Set[ModelValue]] = None,
-            subalgebra2: Optional[Set[ModelValue]] = None):
-        self.has_vsp = has_vsp
-        self.model_name = model_name
-        self.subalgebra1 = subalgebra1
-        self.subalgebra2 = subalgebra2
-
-    def __str__(self):
-        if not self.has_vsp:
-            return f"Model {self.model_name} does not have the variable sharing property."
-        return f"""Model {self.model_name} has the variable sharing property.
-Subalgebra 1: {set_to_str(self.subalgebra1)}
-Subalgebra 2: {set_to_str(self.subalgebra2)}
-"""
-
-def has_vsp(model: Model, interpretation: Dict[Operation, ModelFunction]) -> VSP_Result:
+def has_vsp(model: Model, interpretation: Dict[Operation, ModelFunction]) -> bool:
     """
-    Checks whether a model has the variable
-    sharing property.
+    Tells you whether a model has the
+    variable sharing property.
     """
+
     impfunction = interpretation[Implication]
 
     # Compute I the set of tuples (x, y) where
@@ -133,6 +109,6 @@ def has_vsp(model: Model, interpretation: Dict[Operation, ModelFunction]) -> VSP
                 break
 
         if falsified:
-            return VSP_Result(True, model.name, carrier_set_left, carrier_set_right)
+            return True
 
-    return VSP_Result(False, model.name)
+    return False