Repository: lab-v2/pyreason (language: python)

Issue:
Read Graph from NetworkX Graph object
Feature Request.

While using PyReason as a Python SDK, it would be useful to construct a NetworkX graph as an object programmatically, then pass that into the py engine.  

Allow both `pr.load_graph(path_to_graph)` and `pr.load_graph(G)` where `G = nx.Graph()` to avoid having to lay the XML to disk.
Expose YAML object class in SDK
Feature Request

While using the Python SDK, it would be useful to construct the `fact.yaml`, `rules.yaml` etc programmatically. If these are translated into classes, perhaps expose that class and allow `pr.load_rules(R)` where `R = pr.RulesClass()`


Discussion / hints:
Thanks for this @buddha314. I agree this is important and user friendly, we'll open a PR soon.
Thanks again @buddha314 we'll get this working soon too

Fix patch:
diff --git a/pyproject.toml b/pyproject.toml
index f1bfca4d..0e54f447 100755
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,3 +1,8 @@
 [build-system]
 requires = ['setuptools>=42']
 build-backend = 'setuptools.build_meta'
+
+[tool.pytest.ini_options]
+pythonpath = [
+  "."
+]
diff --git a/pyreason/__init__.py b/pyreason/__init__.py
index 7432abd6..59d34294 100755
--- a/pyreason/__init__.py
+++ b/pyreason/__init__.py
@@ -3,7 +3,8 @@
 package_path = os.path.abspath(os.path.dirname(__file__))
 cache_path = os.path.join(package_path, 'cache')
 cache_status_path = os.path.join(package_path, '.cache_status.yaml')
-os.environ['NUMBA_CACHE_DIR'] =  cache_path
+os.environ['NUMBA_CACHE_DIR'] = cache_path
+
 
 from pyreason.pyreason import *
 import yaml
@@ -15,17 +16,11 @@
 if not cache_status['initialized']:
     print('Imported PyReason for the first time. Initializing ... this will take a minute')
     graph_path = os.path.join(package_path, 'examples', 'hello-world', 'friends.graphml')
-    labels_path = os.path.join(package_path, 'examples', 'hello-world', 'labels.yaml')
-    facts_path = os.path.join(package_path, 'examples', 'hello-world', 'facts.yaml')
-    rules_path = os.path.join(package_path, 'examples', 'hello-world', 'rules.yaml')
-    ipl_path = os.path.join(package_path, 'examples', 'hello-world', 'ipl.yaml')
 
     settings.verbose = False
-    load_graph(graph_path)
-    load_labels(labels_path)
-    load_facts(facts_path)
-    load_rules(rules_path)
-    load_inconsistent_predicate_list(ipl_path)
+    load_graphml(graph_path)
+    add_rule('popular(x) <-1 popular(y), Friends(x,y), owns(y,z), owns(x,z)', 'popular_rule')
+    add_fact(Fact('popular-fact', 'Mary', 'popular', [1, 1], 0, 2))
     reason(timesteps=2)
 
     # Update cache status
diff --git a/pyreason/pyreason.py b/pyreason/pyreason.py
index 978424b6..c0166975 100755
--- a/pyreason/pyreason.py
+++ b/pyreason/pyreason.py
@@ -1,19 +1,22 @@
 # This is the file that will be imported when "import pyreason" is called. All content will be run automatically
+import networkx as nx
 import numba
 import time
 import sys
 import warnings
-import numpy as np
+import pandas as pd
 import memory_profiler as mp
-from typing import List, Type
+from typing import List, Type, Callable, Tuple
 
 from pyreason.scripts.utils.output import Output
 from pyreason.scripts.utils.filter import Filter
 from pyreason.scripts.program.program import Program
 from pyreason.scripts.utils.graphml_parser import GraphmlParser
 import pyreason.scripts.utils.yaml_parser as yaml_parser
+import pyreason.scripts.utils.rule_parser as rule_parser
 import pyreason.scripts.numba_wrapper.numba_types.label_type as label
 import pyreason.scripts.numba_wrapper.numba_types.rule_type as rule
+from pyreason.scripts.facts.fact import Fact
 import pyreason.scripts.numba_wrapper.numba_types.fact_node_type as fact_node
 import pyreason.scripts.numba_wrapper.numba_types.fact_edge_type as fact_edge
 import pyreason.scripts.numba_wrapper.numba_types.interval_type as interval
@@ -35,6 +38,7 @@ def __init__(self):
         self.__inconsistency_check = True
         self.__static_graph_facts = True
         self.__store_interpretation_changes = True
+        self.__parallel_computing = False
 
     @property
     def verbose(self) -> bool:
@@ -144,6 +148,15 @@ def store_interpretation_changes(self) -> bool:
         """
         return self.__store_interpretation_changes
 
+    @property
+    def parallel_computing(self) -> bool:
+        """Returns whether to use multiple CPU cores for inference. This will disable cacheing and pyreason will have
+        to be re-compiled at each run - but after compilation it will be faster. Default is False
+
+        :return: bool
+        """
+        return self.__parallel_computing
+
     @verbose.setter
     def verbose(self, value: bool) -> None:
         """Set verbose mode. Default is True
@@ -305,6 +318,19 @@ def store_interpretation_changes(self, value: bool) -> None:
         else:
             self.__store_interpretation_changes = value
 
+    @parallel_computing.setter
+    def parallel_computing(self, value: bool) -> None:
+        """Whether to use multiple CPU cores for inference. This will disable cacheing and pyreason will have
+        to be re-compiled at each run - but after compilation it will be faster. Default is False
+
+        :param value: Whether to make inference run on parallel hardware (multiple CPU cores)
+        :raises TypeError: If not bool raise error
+        """
+        if not isinstance(value, bool):
+            raise TypeError('value has to be a bool')
+        else:
+            self.__parallel_computing = value
+
 
 # VARIABLES
 __graph = None
@@ -322,6 +348,8 @@ def store_interpretation_changes(self, value: bool) -> None:
 __specific_graph_node_labels = None
 __specific_graph_edge_labels = None
 
+__annotation_functions = []
+
 __timestamp = ''
 __program = None
 
@@ -340,8 +368,16 @@ def reset():
     __edge_labels = None
 
 
+def reset_rules():
+    """
+    Resets rules to none
+    """
+    global __rules
+    __rules = None
+
+
 # FUNCTIONS
-def load_graph(path: str) -> None:
+def load_graphml(path: str) -> None:
     """Loads graph from GraphMl file path into program
 
     :param path: Path for the GraphMl file
@@ -359,7 +395,29 @@ def load_graph(path: str) -> None:
         __non_fluent_graph_facts_edge = numba.typed.List.empty_list(fact_edge.fact_type)
         __specific_graph_node_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(numba.types.string))
         __specific_graph_edge_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(numba.types.Tuple((numba.types.string, numba.types.string))))
+
+
+def load_graph(graph: nx.DiGraph) -> None:
+    """Load a networkx DiGraph into pyreason
+
+    :param graph: Networkx DiGraph object to load into pyreason
+    :type graph: nx.DiGraph
+    :return: None
+    """
+    global __graph, __graphml_parser, __non_fluent_graph_facts_node, __non_fluent_graph_facts_edge, __specific_graph_node_labels, __specific_graph_edge_labels, settings
     
+    # Load graph
+    __graph = __graphml_parser.load_graph(graph)
+
+    # Graph attribute parsing
+    if settings.graph_attribute_parsing:
+        __non_fluent_graph_facts_node, __non_fluent_graph_facts_edge, __specific_graph_node_labels, __specific_graph_edge_labels = __graphml_parser.parse_graph_attributes(settings.static_graph_facts)
+    else:
+        __non_fluent_graph_facts_node = numba.typed.List.empty_list(fact_node.fact_type)
+        __non_fluent_graph_facts_edge = numba.typed.List.empty_list(fact_edge.fact_type)
+        __specific_graph_node_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(numba.types.string))
+        __specific_graph_edge_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(numba.types.Tuple((numba.types.string, numba.types.string))))
+
 
 def load_labels(path: str) -> None:
     """Load labels from YAML file path into program
@@ -401,105 +459,86 @@ def load_inconsistent_predicate_list(path: str) -> None:
     __ipl = yaml_parser.parse_ipl(path)
 
 
-def add_rules_from_text(rule_text: str, name: str) -> None:
+def add_rule(rule_text: str, name: str, infer_edges: bool = False, set_static: bool = False, immediate_rule: bool = False) -> None:
     """Add a rule to pyreason from text format. This format is not as modular as the YAML format.
-    1. It is not possible to specify target criteria
-    2. It is not possible to specify delta_t. delta_t=0 by default.
-    3. It is not possible to specify thresholds. Threshold is greater than or equal to 1 by default
-    4. It is not possible to add edges between subsets of nodes
-    5. It is not possible to have an annotation function. We set to [1,1] by default
-    6. It is not possible to have weights for different clauses. Weights are 1 by default with bias 0
+    1. It is not possible to specify thresholds. Threshold is greater than or equal to 1 by default
+    2. It is not possible to have an annotation function.
+    3. It is not possible to have weights for different clauses. Weights are 1 by default with bias 0
+    TODO: Add threshold class where we can pass this as a parameter
+    TODO: Add weights as a parameter
+    TODO: Add annotation function and bounds as a parameter
 
     Example:
-    `'pred1(x,y) <- pred2(a, b), pred3(b, c)'`
+    `'pred1(x,y) : [0.2, 1] <- pred2(a, b) : [1,1], pred3(b, c)'`
 
     :param rule_text: The rule in text format
     :param name: The name of the rule. This will appear in the rule trace
+    :param infer_edges: Whether to infer new edges after edge rule fires
+    :param set_static: Whether to set the atom in the head as static if the rule fires. The bounds will no longer change
+    :param immediate_rule: Whether the rule is immediate. Immediate rules check for more applicable rules immediately after being applied
     """
     global __rules
 
-    # First remove all spaces from line
-    r = rule_text.replace(' ', '')
-
-    # Separate into head and body
-    head, body = r.split('<-')
-
-    # Separate clauses in body
-    body = body[:-1].replace(')', '))') + ')'
-    body = body.split('),')
-
-    # Find the target predicate
-    idx = head.find('(')
-    target = head[:idx]
-    target = label.Label(target)
-
-    # Variable(s) in the head of the rule
-    head_variables = head[idx + 1:-1].split(',')
-
-    # Target criteria is empty for this text format to pyreason rule
-    target_criteria = numba.typed.List.empty_list(numba.types.Tuple((label.label_type, interval.interval_type)))
-
-    # Get the variables in the body
-    body_predicates = []
-    body_variables = []
-    for clause in body:
-        idx = clause.find('(')
-        body_predicates.append(clause[:idx])
-        body_variables.append(clause[idx+1:-1].split(','))
-
-    # Replace the variables in the body with source/target if they match the variables in the head
-    head_var_map = {0: 'source', 1: 'target'}
-    for i in range(len(body_variables)):
-        for j in range(len(body_variables[i])):
-            # Loop through the head variables and see if there's a match
-            for k in range(len(head_variables)):
-                if body_variables[i][j] == head_variables[k]:
-                    body_variables[i][j] = head_var_map[k] if len(head_variables) == 2 else 'target'
-
-    # Start setting up neigh_criteria
-    # neigh_criteria = [c1, c2, c3, c4]
-    # thresholds = [t1, t2, t3, t4]
-
-    # Array of thresholds to keep track of for each neighbor criterion. Form [(comparison, (number/percent, total/available), thresh)]
-    thresholds = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, numba.types.UniTuple(numba.types.string, 2), numba.types.float64)))
-
-    # Array to store clauses for nodes: node/edge, [subset]/[subset1, subset2], label, interval
-    neigh_criteria = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, numba.types.UniTuple(numba.types.string, 2), label.label_type, interval.interval_type)))
-
-    # Loop though clauses
-    for predicate, variables in zip(body_predicates, body_variables):
-        # Neigh criteria
-        clause_type = 'node' if len(variables) == 1 else 'edge'
-        subset = (variables[0], variables[0]) if clause_type == 'node' else (variables[0], variables[1])
-        l = label.Label(predicate)
-        bnd = interval.closed(1, 1)
-        neigh_criteria.append((clause_type, subset, l, bnd))
-
-        # Threshold
-        quantifier = 'greater_equal'
-        quantifier_type = ('number', 'total')
-        thresh = 1
-        thresholds.append((quantifier, quantifier_type, thresh))
-
-    # Cannot add edges
-    edges = ('', '', label.Label(''))
-
-    # Bound to set atom if rule fires
-    bnd = interval.closed(1, 1)
-    ann_fn = ''
-    ann_label = label.Label('')
-
-    weights = np.ones(len(body_predicates), dtype=np.float64)
-    weights = np.append(weights, 0)
-
-    r = rule.Rule(name, target, target_criteria, numba.types.uint16(0), neigh_criteria, bnd, thresholds, ann_fn, ann_label, weights, edges)
-
+    r = rule_parser.parse_rule(rule_text, name, infer_edges, set_static, immediate_rule)
     # Add to collection of rules
     if __rules is None:
         __rules = numba.typed.List.empty_list(rule.rule_type)
     __rules.append(r)
 
 
+def add_rules_from_file(file_path: str, infer_edges: bool = False) -> None:
+    """ Add a set of rules from a text file
+
+    :param file_path: Path to the text file containing rules
+    :type file_path: str
+    :param infer_edges: Whether to infer edges on these rules if an edge doesn't exist between head variables and the body of the rule is satisfied
+    :type infer_edges: bool
+    :return: None
+    """
+    with open(file_path, 'r') as file:
+        rules = [line.rstrip() for line in file if line.rstrip() != '' and line.rstrip()[0] != '#']
+
+    rule_offset = 0 if __rules is None else len(__rules)
+    for i, r in enumerate(rules):
+        add_rule(r, f'rule_{i+rule_offset}', infer_edges)
+
+
+def add_fact(pyreason_fact: Fact) -> None:
+    """Add a PyReason fact to the program.
+
+    :param pyreason_fact: PyReason fact created using pr.Fact(...)
+    :return: None
+    """
+    global __node_facts, __edge_facts
+
+    if pyreason_fact.type == 'node':
+        f = fact_node.Fact(pyreason_fact.name, pyreason_fact.component, pyreason_fact.label, pyreason_fact.interval, pyreason_fact.t_lower, pyreason_fact.t_upper, pyreason_fact.static)
+        if __node_facts is None:
+            __node_facts = numba.typed.List.empty_list(fact_node.fact_type)
+        __node_facts.append(f)
+
+    else:
+        f = fact_edge.Fact(pyreason_fact.name, pyreason_fact.component, pyreason_fact.label, pyreason_fact.interval, pyreason_fact.t_lower, pyreason_fact.t_upper, pyreason_fact.static)
+        if __edge_facts is None:
+            __edge_facts = numba.typed.List.empty_list(fact_edge.fact_type)
+        __edge_facts.append(f)
+
+
+def add_annotation_function(function: Callable) -> None:
+    """Function to add annotation functions to PyReason. The added functions can be used in rules
+
+    :param function: Function to be added. This has to be under a numba `njit` decorator. function has signature: two parameters as input -- annotations, weights
+    :type function: Callable
+    :return: None
+    """
+    global __annotation_functions
+    # Make sure that the functions are jitted so that they can be passed around in other jitted functions
+    # TODO: Remove if necessary
+    # assert hasattr(function, 'nopython_signatures'), 'The function to be added has to be under a `numba.njit` decorator'
+
+    __annotation_functions.append(function)
+
+
 def reason(timesteps: int=-1, convergence_threshold: int=-1, convergence_bound_threshold: float=-1, again: bool=False, node_facts: List[Type[fact_node.Fact]]=None, edge_facts: List[Type[fact_edge.Fact]]=None):
     """Function to start the main reasoning process. Graph and rules must already be loaded.
 
@@ -564,8 +603,10 @@ def _reason(timesteps, convergence_threshold, convergence_bound_threshold):
 
     if __node_facts is None and __edge_facts is None:
         if settings.verbose:
-            warnings.warn('Facts yaml file has not been loaded. Use `load_facts`. Only graph attributes will be used as facts\n')
+            warnings.warn('Facts have not been loaded. Use `add_fact` or `load_facts`. Only graph attributes will be used as facts\n')
+    if __node_facts is None:
         __node_facts = numba.typed.List.empty_list(fact_node.fact_type)
+    if __edge_facts is None:
         __edge_facts = numba.typed.List.empty_list(fact_edge.fact_type)
 
     if __ipl is None:
@@ -586,8 +627,10 @@ def _reason(timesteps, convergence_threshold, convergence_bound_threshold):
         else:
             __specific_edge_labels[label_name] = edges
 
-    all_node_facts = numba.typed.List(__node_facts)
-    all_edge_facts = numba.typed.List(__edge_facts)
+    all_node_facts = numba.typed.List.empty_list(fact_node.fact_type)
+    all_edge_facts = numba.typed.List.empty_list(fact_edge.fact_type)
+    all_node_facts.extend(numba.typed.List(__node_facts))
+    all_edge_facts.extend(numba.typed.List(__edge_facts))
     all_node_facts.extend(__non_fluent_graph_facts_node)
     all_edge_facts.extend(__non_fluent_graph_facts_edge)
 
@@ -595,8 +638,11 @@ def _reason(timesteps, convergence_threshold, convergence_bound_threshold):
     if not settings.store_interpretation_changes:
         settings.atom_trace = False
 
+    # Convert list of annotation functions into tuple to be numba compatible
+    annotation_functions = tuple(__annotation_functions)
+
     # Setup logical program
-    __program = Program(__graph, all_node_facts, all_edge_facts, __rules, __ipl, settings.reverse_digraph, settings.atom_trace, settings.save_graph_attributes_to_trace, settings.canonical, settings.inconsistency_check, settings.store_interpretation_changes)
+    __program = Program(__graph, all_node_facts, all_edge_facts, __rules, __ipl, annotation_functions, settings.reverse_digraph, settings.atom_trace, settings.save_graph_attributes_to_trace, settings.canonical, settings.inconsistency_check, settings.store_interpretation_changes, settings.parallel_computing)
     __program.available_labels_node = __node_labels
     __program.available_labels_edge = __edge_labels
     __program.specific_node_labels = __specific_node_labels
@@ -630,7 +676,7 @@ def _reason_again(timesteps, convergence_threshold, convergence_bound_threshold,
 
 
 def save_rule_trace(interpretation, folder: str='./'):
-    """Saves the trace of the program. This includes every change that has occured to the interpretation. If `atom_trace` was set to true
+    """Saves the trace of the program. This includes every change that has occurred to the interpretation. If `atom_trace` was set to true
     this gives us full explainability of why interpretations changed
 
     :param interpretation: the output of `pyreason.reason()`, the final interpretation
@@ -644,6 +690,22 @@ def save_rule_trace(interpretation, folder: str='./'):
     output.save_rule_trace(interpretation, folder)
 
 
+def get_rule_trace(interpretation) -> Tuple[pd.DataFrame, pd.DataFrame]:
+    """Returns the trace of the program as 2 pandas dataframes (one for nodes, one for edges).
+    This includes every change that has occurred to the interpretation. If `atom_trace` was set to true
+    this gives us full explainability of why interpretations changed
+
+    :param interpretation: the output of `pyreason.reason()`, the final interpretation
+    :returns two pandas dataframes (nodes, edges) representing the changes that occurred during reasoning
+    """
+    global __timestamp, settings
+
+    assert settings.store_interpretation_changes, 'store interpretation changes setting is off, turn on to save rule trace'
+
+    output = Output(__timestamp)
+    return output.get_rule_trace(interpretation)
+
+
 def filter_and_sort_nodes(interpretation, labels: List[str], bound: interval.Interval=interval.closed(0,1), sort_by: str='lower', descending: bool=True):
     """Filters and sorts the node changes in the interpretation and returns as a list of Pandas dataframes that are easy to access
 
diff --git a/pyreason/scripts/facts/fact.py b/pyreason/scripts/facts/fact.py
new file mode 100644
index 00000000..1004cec3
--- /dev/null
+++ b/pyreason/scripts/facts/fact.py
@@ -0,0 +1,47 @@
+import pyreason.scripts.numba_wrapper.numba_types.interval_type as interval
+import pyreason.scripts.numba_wrapper.numba_types.label_type as label
+
+from typing import Tuple
+from typing import List
+from typing import Union
+
+
+class Fact:
+    def __init__(self, name: str, component: Union[str, Tuple[str, str]], attribute: str, bound: Union[interval.Interval, List[float]], start_time: int, end_time: int, static: bool = False):
+        """Define a PyReason fact that can be loaded into the program using `pr.add_fact()`
+
+        :param name: The name of the fact. This will appear in the trace so that you know when it was applied
+        :type name: str
+        :param component: The node or edge that whose attribute you want to change
+        :type component: str | Tuple[str, str]
+        :param attribute: The attribute you would like to change for the specified node/edge
+        :type attribute: str
+        :param bound: The bound to which you'd like to set the attribute corresponding to the specified node/edge
+        :type bound: interval.Interval | List[float]
+        :param start_time: The timestep at which this fact becomes active
+        :type start_time: int
+        :param end_time: The last timestep this fact is active
+        :type end_time: int
+        :param static: If the fact should be active for the entire program. In which case `start_time` and `end_time` will be ignored
+        :type static: bool
+        """
+        self.name = name
+        self.t_upper = end_time
+        self.t_lower = start_time
+        self.component = component
+        self.label = attribute
+        self.interval = bound
+        self.static = static
+
+        # Check if it is a node fact or edge fact
+        if isinstance(self.component, str):
+            self.type = 'node'
+        else:
+            self.type = 'edge'
+
+        # Set label to correct type
+        self.label = label.Label(attribute)
+
+        # Set bound to correct type
+        if isinstance(bound, list):
+            self.interval = interval.closed(*bound)
diff --git a/pyreason/scripts/facts/fact_edge.py b/pyreason/scripts/facts/fact_edge.py
index 09283ce7..935d40f4 100755
--- a/pyreason/scripts/facts/fact_edge.py
+++ b/pyreason/scripts/facts/fact_edge.py
@@ -28,9 +28,12 @@ def get_time_upper(self):
         return self._t_upper
 
     def __str__(self):
-        net_diff_dict = {}
-        net_diff_dict["component"] = str(self._component)
-        net_diff_dict["label"] = str(self._label)
-        net_diff_dict["confidence"] = self._interval.to_str()
-        net_diff_dict["time"] = '[' + str(self._t_lower) + ',' + str(self._t_upper) + ']'
-        return str(net_diff_dict)
\ No newline at end of file
+        fact = {
+            "type": 'pyreason edge fact',
+            "name": self._name,
+            "component": self._component,
+            "label": self._label,
+            "confidence": self._interval,
+            "time": '[' + str(self._t_lower) + ',' + str(self._t_upper) + ']'
+        }
+        return fact
diff --git a/pyreason/scripts/facts/fact_node.py b/pyreason/scripts/facts/fact_node.py
index 09283ce7..92e97c88 100755
--- a/pyreason/scripts/facts/fact_node.py
+++ b/pyreason/scripts/facts/fact_node.py
@@ -28,9 +28,12 @@ def get_time_upper(self):
         return self._t_upper
 
     def __str__(self):
-        net_diff_dict = {}
-        net_diff_dict["component"] = str(self._component)
-        net_diff_dict["label"] = str(self._label)
-        net_diff_dict["confidence"] = self._interval.to_str()
-        net_diff_dict["time"] = '[' + str(self._t_lower) + ',' + str(self._t_upper) + ']'
-        return str(net_diff_dict)
\ No newline at end of file
+        fact = {
+            "type": 'pyreason node fact',
+            "name": self._name,
+            "component": self._component,
+            "label": self._label,
+            "confidence": self._interval,
+            "time": '[' + str(self._t_lower) + ',' + str(self._t_upper) + ']'
+        }
+        return fact
diff --git a/pyreason/scripts/interpretation/interpretation.py b/pyreason/scripts/interpretation/interpretation.py
index a8dbf563..f7a271c6 100755
--- a/pyreason/scripts/interpretation/interpretation.py
+++ b/pyreason/scripts/interpretation/interpretation.py
@@ -1,21 +1,41 @@
 import pyreason.scripts.numba_wrapper.numba_types.world_type as world
 import pyreason.scripts.numba_wrapper.numba_types.label_type as label
 import pyreason.scripts.numba_wrapper.numba_types.interval_type as interval
-import pyreason.scripts.annotation_functions.annotation_functions as ann_fn
+from pyreason.scripts.interpretation.interpretation_dict import InterpretationDict
 
 import numba
-from numba import objmode
+from numba import objmode, prange
 
 
 # Types for the dictionaries
 node_type = numba.types.string
 edge_type = numba.types.UniTuple(numba.types.string, 2)
 
-# Type for storing list of qualified nodes
+# Type for storing list of qualified nodes/edges
 list_of_nodes = numba.types.ListType(node_type)
-
-# Type for storing int tuple
-int_tuple = numba.types.UniTuple(numba.types.int64, 2)
+list_of_edges = numba.types.ListType(edge_type)
+
+# Type for facts to be applied
+facts_to_be_applied_node_type = numba.types.Tuple((numba.types.uint16, node_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean))
+facts_to_be_applied_edge_type = numba.types.Tuple((numba.types.uint16, edge_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean))
+
+# Type for returning list of applicable rules for a certain rule
+# node/edge, annotations, qualified nodes, qualified edges, edges to be added
+node_applicable_rule_type = numba.types.Tuple((
+	node_type,
+	numba.types.ListType(numba.types.ListType(interval.interval_type)),
+	numba.types.ListType(numba.types.ListType(node_type)),
+	numba.types.ListType(numba.types.ListType(edge_type)),
+	numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type))
+))
+
+edge_applicable_rule_type = numba.types.Tuple((
+	edge_type,
+	numba.types.ListType(numba.types.ListType(interval.interval_type)),
+	numba.types.ListType(numba.types.ListType(node_type)),
+	numba.types.ListType(numba.types.ListType(edge_type)),
+	numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type))
+))
 
 
 class Interpretation:
@@ -24,9 +44,10 @@ class Interpretation:
 	specific_node_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(node_type))
 	specific_edge_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(edge_type))
 
-	def __init__(self, graph, ipl, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes):
+	def __init__(self, graph, ipl, annotation_functions, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes):
 		self.graph = graph
 		self.ipl = ipl
+		self.annotation_functions = annotation_functions
 		self.reverse_graph = reverse_graph
 		self.atom_trace = atom_trace
 		self.save_graph_attributes_to_rule_trace = save_graph_attributes_to_rule_trace
@@ -43,10 +64,10 @@ def __init__(self, graph, ipl, reverse_graph, atom_trace, save_graph_attributes_
 		self.rules_to_be_applied_edge_trace = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(numba.types.ListType(node_type)), numba.types.ListType(numba.types.ListType(edge_type)), numba.types.string)))
 		self.facts_to_be_applied_node_trace = numba.typed.List.empty_list(numba.types.string)
 		self.facts_to_be_applied_edge_trace = numba.typed.List.empty_list(numba.types.string)
-		self.rules_to_be_applied_node = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, node_type, label.label_type, interval.interval_type, numba.types.boolean)))
-		self.rules_to_be_applied_edge = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, edge_type, label.label_type, interval.interval_type, numba.types.boolean)))
-		self.facts_to_be_applied_node = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, node_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean)))
-		self.facts_to_be_applied_edge = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, edge_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean)))
+		self.rules_to_be_applied_node = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, node_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean)))
+		self.rules_to_be_applied_edge = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, edge_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean)))
+		self.facts_to_be_applied_node = numba.typed.List.empty_list(facts_to_be_applied_node_type)
+		self.facts_to_be_applied_edge = numba.typed.List.empty_list(facts_to_be_applied_edge_type)
 		self.edges_to_be_added_node_rule = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type)))
 		self.edges_to_be_added_edge_rule = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type)))
 
@@ -72,14 +93,32 @@ def __init__(self, graph, ipl, reverse_graph, atom_trace, save_graph_attributes_
 
 		self.interpretations_node = self._init_interpretations_node(numba.typed.List(self.graph.nodes()), self.available_labels_node, self.specific_node_labels)
 		self.interpretations_edge = self._init_interpretations_edge(numba.typed.List(self.graph.edges()), self.available_labels_edge, self.specific_edge_labels)
-		
-		# Setup graph neighbors
+
+		# Setup graph neighbors and reverse neighbors
 		self.neighbors = numba.typed.Dict.empty(key_type=node_type, value_type=numba.types.ListType(node_type))
 		for n in self.graph.nodes():
 			l = numba.typed.List.empty_list(node_type)
 			[l.append(neigh) for neigh in self.graph.neighbors(n)]
 			self.neighbors[n] = l
 
+		self.reverse_neighbors = self._init_reverse_neighbors(self.neighbors)
+
+	@staticmethod
+	@numba.njit(cache=True)
+	def _init_reverse_neighbors(neighbors):
+		reverse_neighbors = numba.typed.Dict.empty(key_type=node_type, value_type=list_of_nodes)
+		for n, neighbor_nodes in neighbors.items():
+			for neighbor_node in neighbor_nodes:
+				if neighbor_node in reverse_neighbors and n not in reverse_neighbors[neighbor_node]:
+					reverse_neighbors[neighbor_node].append(n)
+				else:
+					reverse_neighbors[neighbor_node] = numba.typed.List([n])
+			# This makes sure each node has a value
+			if n not in reverse_neighbors:
+				reverse_neighbors[n] = numba.typed.List.empty_list(node_type)
+
+		return reverse_neighbors
+
 	@staticmethod
 	@numba.njit(cache=True)
 	def _init_interpretations_node(nodes, available_labels, specific_labels):
@@ -151,7 +190,7 @@ def _init_facts(facts_node, facts_edge, facts_to_be_applied_node, facts_to_be_ap
 		return max_time
 
 	def _start_fp(self, rules, max_facts_time, verbose, again):
-		fp_cnt, t = self.reason(self.interpretations_node, self.interpretations_edge, self.tmax, self.prev_reasoning_data, rules, self.nodes, self.edges, self.neighbors, self.rules_to_be_applied_node, self.rules_to_be_applied_edge, self.edges_to_be_added_node_rule, self.edges_to_be_added_edge_rule, self.rules_to_be_applied_node_trace, self.rules_to_be_applied_edge_trace, self.facts_to_be_applied_node, self.facts_to_be_applied_edge, self.facts_to_be_applied_node_trace, self.facts_to_be_applied_edge_trace, self.available_labels_node, self.available_labels_edge, self.specific_node_labels, self.specific_edge_labels, self.ipl, self.rule_trace_node, self.rule_trace_edge, self.rule_trace_node_atoms, self.rule_trace_edge_atoms, self.reverse_graph, self.atom_trace, self.save_graph_attributes_to_rule_trace, self.canonical, self.inconsistency_check, self.store_interpretation_changes, max_facts_time, self._convergence_mode, self._convergence_delta, verbose, again)
+		fp_cnt, t = self.reason(self.interpretations_node, self.interpretations_edge, self.tmax, self.prev_reasoning_data, rules, self.nodes, self.edges, self.neighbors, self.reverse_neighbors, self.rules_to_be_applied_node, self.rules_to_be_applied_edge, self.edges_to_be_added_node_rule, self.edges_to_be_added_edge_rule, self.rules_to_be_applied_node_trace, self.rules_to_be_applied_edge_trace, self.facts_to_be_applied_node, self.facts_to_be_applied_edge, self.facts_to_be_applied_node_trace, self.facts_to_be_applied_edge_trace, self.ipl, self.rule_trace_node, self.rule_trace_edge, self.rule_trace_node_atoms, self.rule_trace_edge_atoms, self.reverse_graph, self.atom_trace, self.save_graph_attributes_to_rule_trace, self.canonical, self.inconsistency_check, self.store_interpretation_changes, max_facts_time, self.annotation_functions, self._convergence_mode, self._convergence_delta, verbose, again)
 		self.time = t - 1
 		# If we need to reason again, store the next timestep to start from
 		self.prev_reasoning_data[0] = t
@@ -161,43 +200,35 @@ def _start_fp(self, rules, max_facts_time, verbose, again):
 
 	@staticmethod
 	@numba.njit(cache=True)
-	def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data, rules, nodes, edges, neighbors, rules_to_be_applied_node, rules_to_be_applied_edge, edges_to_be_added_node_rule, edges_to_be_added_edge_rule, rules_to_be_applied_node_trace, rules_to_be_applied_edge_trace, facts_to_be_applied_node, facts_to_be_applied_edge, facts_to_be_applied_node_trace, facts_to_be_applied_edge_trace, labels_node, labels_edge, specific_labels_node, specific_labels_edge, ipl, rule_trace_node, rule_trace_edge, rule_trace_node_atoms, rule_trace_edge_atoms, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes, max_facts_time, convergence_mode, convergence_delta, verbose, again):
+	def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data, rules, nodes, edges, neighbors, reverse_neighbors, rules_to_be_applied_node, rules_to_be_applied_edge, edges_to_be_added_node_rule, edges_to_be_added_edge_rule, rules_to_be_applied_node_trace, rules_to_be_applied_edge_trace, facts_to_be_applied_node, facts_to_be_applied_edge, facts_to_be_applied_node_trace, facts_to_be_applied_edge_trace, ipl, rule_trace_node, rule_trace_edge, rule_trace_node_atoms, rule_trace_edge_atoms, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes, max_facts_time, annotation_functions, convergence_mode, convergence_delta, verbose, again):
 		t = prev_reasoning_data[0]
 		fp_cnt = prev_reasoning_data[1]
 		max_rules_time = 0
 		timestep_loop = True
-		facts_to_remove_idx = numba.typed.List.empty_list(numba.types.int64)
+		facts_to_be_applied_node_new = numba.typed.List.empty_list(facts_to_be_applied_node_type)
+		facts_to_be_applied_edge_new = numba.typed.List.empty_list(facts_to_be_applied_edge_type)
 		rules_to_remove_idx = numba.typed.List.empty_list(numba.types.int64)
 		while timestep_loop:
 			if t==tmax:
-				timestep_loop=False
+				timestep_loop = False
 			if verbose:
 				with objmode():
 					print('Timestep:', t, flush=True)
 			# Reset Interpretation at beginning of timestep if non-canonical
 			if t>0 and not canonical:
 				# Reset nodes (only if not static)
-				# General labels
 				for n in nodes:
-					for l in labels_node:
-						if not interpretations_node[n].world[l].is_static():
-							interpretations_node[n].world[l].reset()
-				# Specific labels
-				for l, ns in specific_labels_node.items():
-					for n in ns:
-						if not interpretations_node[n].world[l].is_static():
-							interpretations_node[n].world[l].reset()				
-				# Reset edges
-				# General labels
+					w = interpretations_node[n].world
+					for l in w:
+						if not w[l].is_static():
+							w[l].reset()
+
+				# Reset edges (only if not static)
 				for e in edges:
-					for l in labels_edge:
-						if not interpretations_edge[e].world[l].is_static():
-							interpretations_edge[e].world[l].reset()
-				# Specific labels
-				for l, es in specific_labels_edge.items():
-					for e in es:
-						if not interpretations_edge[e].world[l].is_static():
-							interpretations_edge[e].world[l].reset()
+					w = interpretations_edge[e].world
+					for l in w:
+						if not w[l].is_static():
+							w[l].reset()
 
 			# Convergence parameters
 			changes_cnt = 0
@@ -209,17 +240,21 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 			immediate_edge_rule_fire = False
 			immediate_rule_applied = False
 			# When delta_t = 0, we don't want to check the same rule with the same node/edge after coming back to the fp operator
-			nodes_to_skip = numba.typed.List.empty_list(int_tuple)
-			edges_to_skip = numba.typed.List.empty_list(int_tuple)
+			nodes_to_skip = numba.typed.Dict.empty(key_type=numba.types.int64, value_type=list_of_nodes)
+			edges_to_skip = numba.typed.Dict.empty(key_type=numba.types.int64, value_type=list_of_edges)
+			# Initialize the above
+			for i in range(len(rules)):
+				nodes_to_skip[i] = numba.typed.List.empty_list(node_type)
+				edges_to_skip[i] = numba.typed.List.empty_list(edge_type)
 
 			# Start by applying facts
 			# Nodes
-			facts_to_remove_idx.clear()
+			facts_to_be_applied_node_new.clear()
 			for i in range(len(facts_to_be_applied_node)):
 				if facts_to_be_applied_node[i][0]==t:
 					comp, l, bnd, static, graph_attribute = facts_to_be_applied_node[i][1], facts_to_be_applied_node[i][2], facts_to_be_applied_node[i][3], facts_to_be_applied_node[i][4], facts_to_be_applied_node[i][5]
 					# Check if bnd is static. Then no need to update, just add to rule trace, check if graph attribute and add ipl complement to rule trace as well
-					if interpretations_node[comp].world[l].is_static():
+					if l in interpretations_node[comp].world and interpretations_node[comp].world[l].is_static():
 						# Check if we should even store any of the changes to the rule trace etc.
 						# Inverse of this is: if not save_graph_attributes_to_rule_trace and graph_attribute
 						if (save_graph_attributes_to_rule_trace or not graph_attribute) and store_interpretation_changes:
@@ -264,21 +299,23 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 									changes_cnt += changes
 
 					if static:
-						facts_to_be_applied_node[i] = (numba.types.uint16(facts_to_be_applied_node[i][0]+1), comp, l, bnd, static, graph_attribute)
-					else:
-						# Add to list to be removed later
-						facts_to_remove_idx.append(i)
+						facts_to_be_applied_node_new.append((numba.types.uint16(facts_to_be_applied_node[i][0]+1), comp, l, bnd, static, graph_attribute))
 
-			# Delete facts that are not static
-			facts_to_be_applied_node[:] = numba.typed.List([facts_to_be_applied_node[i] for i in range(len(facts_to_be_applied_node)) if i not in facts_to_remove_idx])
+				# If time doesn't match, fact to be applied later
+				else:
+					facts_to_be_applied_node_new.append(facts_to_be_applied_node[i])
+
+			# Update list of facts with ones that have not been applied yet (delete applied facts)
+			facts_to_be_applied_node[:] = facts_to_be_applied_node_new.copy()
+			facts_to_be_applied_node_new.clear()
 
 			# Edges
-			facts_to_remove_idx.clear()
+			facts_to_be_applied_edge_new.clear()
 			for i in range(len(facts_to_be_applied_edge)):
 				if facts_to_be_applied_edge[i][0]==t:
 					comp, l, bnd, static, graph_attribute = facts_to_be_applied_edge[i][1], facts_to_be_applied_edge[i][2], facts_to_be_applied_edge[i][3], facts_to_be_applied_edge[i][4], facts_to_be_applied_edge[i][5]
 					# Check if bnd is static. Then no need to update, just add to rule trace, check if graph attribute, and add ipl complement to rule trace as well
-					if interpretations_edge[comp].world[l].is_static():
+					if l in interpretations_edge[comp].world and interpretations_edge[comp].world[l].is_static():
 						# Inverse of this is: if not save_graph_attributes_to_rule_trace and graph_attribute
 						if (save_graph_attributes_to_rule_trace or not graph_attribute) and store_interpretation_changes:
 							rule_trace_edge.append((numba.types.uint16(t), numba.types.uint16(fp_cnt), comp, l, interpretations_edge[comp].world[l]))
@@ -321,13 +358,15 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 									changes_cnt += changes
 
 					if static:
-						facts_to_be_applied_edge[i] = (numba.types.uint16(facts_to_be_applied_edge[i][0]+1), comp, l, bnd, static, graph_attribute)
-					else:
-						# Add to list to be removed later
-						facts_to_remove_idx.append(i)
+						facts_to_be_applied_edge_new.append((numba.types.uint16(facts_to_be_applied_edge[i][0]+1), comp, l, bnd, static, graph_attribute))
 
-			# Delete facts that are not static
-			facts_to_be_applied_edge[:] = numba.typed.List([facts_to_be_applied_edge[i] for i in range(len(facts_to_be_applied_edge)) if i not in facts_to_remove_idx])
+				# Time doesn't match, fact to be applied later
+				else:
+					facts_to_be_applied_edge_new.append(facts_to_be_applied_edge[i])
+
+			# Update list of facts with ones that have not been applied yet (delete applied facts)
+			facts_to_be_applied_edge[:] = facts_to_be_applied_edge_new.copy()
+			facts_to_be_applied_edge_new.clear()
 
 			in_loop = True
 			while in_loop:
@@ -344,16 +383,16 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 						idx = len(rules_to_be_applied_node) - 1
 
 					if i[0]==t:
-						comp, l, bnd, immediate = i[1], i[2], i[3], i[4]
+						comp, l, bnd, immediate, set_static = i[1], i[2], i[3], i[4], i[5]
 						sources, targets, edge_l = edges_to_be_added_node_rule[idx]
-						edges_added, changes = _add_edges(sources, targets, neighbors, nodes, edges, edge_l, interpretations_node, interpretations_edge)
+						edges_added, changes = _add_edges(sources, targets, neighbors, reverse_neighbors, nodes, edges, edge_l, interpretations_node, interpretations_edge)
 						changes_cnt += changes
 
 						# Update bound for newly added edges. Use bnd to update all edges if label is specified, else use bnd to update normally
 						if edge_l.value!='':
 							for e in edges_added:
 								if check_consistent_edge(interpretations_edge, e, (edge_l, bnd)):
-									u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
+									u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
 
 									update = u or update
 
@@ -367,7 +406,7 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 									if inconsistency_check:
 										resolve_inconsistency_edge(interpretations_edge, e, (edge_l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_edge, rule_trace_edge_atoms, store_interpretation_changes)
 									else:
-										u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
+										u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
 
 										update = u or update
 
@@ -379,7 +418,7 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 						else:
 							# Check for inconsistencies
 							if check_consistent_node(interpretations_node, comp, (l, bnd)):
-								u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode='rule')
+								u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode='rule')
 
 								update = u or update
 								# Update convergence params
@@ -392,7 +431,7 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 								if inconsistency_check:
 									resolve_inconsistency_node(interpretations_node, comp, (l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_node, rule_trace_node_atoms, store_interpretation_changes)
 								else:
-									u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode='rule', override=True)
+									u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode='rule', override=True)
 
 									update = u or update
 									# Update convergence params
@@ -432,16 +471,16 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 						idx = len(rules_to_be_applied_edge) - 1
 
 					if i[0]==t:
-						comp, l, bnd, immediate = i[1], i[2], i[3], i[4]
+						comp, l, bnd, immediate, set_static = i[1], i[2], i[3], i[4], i[5]
 						sources, targets, edge_l = edges_to_be_added_edge_rule[idx]
-						edges_added, changes = _add_edges(sources, targets, neighbors, nodes, edges, edge_l, interpretations_node, interpretations_edge)
+						edges_added, changes = _add_edges(sources, targets, neighbors, reverse_neighbors, nodes, edges, edge_l, interpretations_node, interpretations_edge)
 						changes_cnt += changes
 
 						# Update bound for newly added edges. Use bnd to update all edges if label is specified, else use bnd to update normally
 						if edge_l.value!='':
 							for e in edges_added:
 								if check_consistent_edge(interpretations_edge, e, (edge_l, bnd)):
-									u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
+									u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
 
 									update = u or update
 
@@ -455,7 +494,7 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 									if inconsistency_check:
 										resolve_inconsistency_edge(interpretations_edge, e, (edge_l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_edge, rule_trace_edge_atoms, store_interpretation_changes)
 									else:
-										u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
+										u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
 
 										update = u or update
 
@@ -468,7 +507,7 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 						else:
 							# Check for inconsistencies
 							if check_consistent_edge(interpretations_edge, comp, (l, bnd)):
-								u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
+								u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
 
 								update = u or update
 								# Update convergence params
@@ -481,7 +520,7 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 								if inconsistency_check:
 									resolve_inconsistency_edge(interpretations_edge, comp, (l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_edge, rule_trace_edge_atoms, store_interpretation_changes)
 								else:
-									u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, False, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
+									u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
 
 									update = u or update
 									# Update convergence params
@@ -502,7 +541,7 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 							elif not immediate_edge_rule_fire and u:
 								immediate_rule_applied = True
 								break
-								
+
 				# Remove from rules to be applied and edges to be applied lists after coming out from loop
 				rules_to_be_applied_edge[:] = numba.typed.List([rules_to_be_applied_edge[i] for i in range(len(rules_to_be_applied_edge)) if i not in rules_to_remove_idx])
 				edges_to_be_added_edge_rule[:] = numba.typed.List([edges_to_be_added_edge_rule[i] for i in range(len(edges_to_be_added_edge_rule)) if i not in rules_to_remove_idx])
@@ -521,79 +560,73 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 						immediate_rule = rule.is_immediate_rule()
 						immediate_node_rule_fire = False
 						immediate_edge_rule_fire = False
-						# Go through all nodes and check if any rules apply to them
-						# Only go through everything if the rule can be applied within the given timesteps, or we're running until convergence or it's an immediate rule.
-						# Otherwise, it's an unnecessary loop
+
+						# Only go through if the rule can be applied within the given timesteps, or we're running until convergence
 						delta_t = rule.get_delta()
-						if t+delta_t<=tmax or tmax==-1 or again:
-							for j in range(len(nodes)):
-								if (i, j) in nodes_to_skip:
-									continue
-								n = nodes[j]
-								target_criteria_satisfaction = are_satisfied_node(interpretations_node, n, rule.get_target_criteria())
-								if (target_criteria_satisfaction and is_satisfied_node(interpretations_node, n, (rule.get_target(), interval.closed(0,1)))) or (target_criteria_satisfaction and rule.get_target().value==''):
-									result, annotations, qualified_nodes, qualified_edges, edges_to_add = _is_rule_applicable(interpretations_node, interpretations_edge, neighbors, n, rule.get_neigh_criteria(), rule.get_thresholds(), reverse_graph, rule.get_annotation_function(), rule.get_annotation_label(), rule.get_edges(), atom_trace)
-									if (result and rule.get_target().value=='') or (result and not interpretations_node[n].world[rule.get_target()].is_static()):
-										bnd = influence(rule, annotations, rule.get_weights())
-										max_rules_time = max(max_rules_time, t+delta_t)
-										edges_to_be_added_node_rule.append(edges_to_add)
-										rules_to_be_applied_node.append((numba.types.uint16(t+delta_t), n, rule.get_target(), bnd, immediate_rule))
-										if atom_trace:
-											rules_to_be_applied_node_trace.append((qualified_nodes, qualified_edges, rule.get_name()))
-
-										# We apply a rule on a node/edge only once in each timestep to prevent it from being added to the to_be_added list continuously (this will improve performance
-										nodes_to_skip.append((i, j))
-
-										# Handle loop parameters for the next (maybe) fp operation
-										# If it is a t=0 rule or an immediate rule we want to go back for another fp operation to check for new rules that may fire
-										# Next fp operation we will skip this rule on this node because anyway there won't be an update
-										if delta_t==0:
-											in_loop = True
-											update = False
-										if immediate_rule and delta_t==0:
-											# immediate_rule_fire becomes True because we still need to check for more eligible rules, we're not done.
-											in_loop = True
-											update = True
-											immediate_node_rule_fire = True
-											break
+						if t + delta_t <= tmax or tmax == -1 or again:
+							applicable_node_rules = _ground_node_rule(rule, interpretations_node, interpretations_edge, nodes, neighbors, reverse_neighbors, atom_trace, reverse_graph, nodes_to_skip[i])
+							applicable_edge_rules = _ground_edge_rule(rule, interpretations_node, interpretations_edge, nodes, edges, neighbors, reverse_neighbors, atom_trace, reverse_graph, edges_to_skip[i])
+
+							# Loop through applicable rules and add them to the rules to be applied for later or next fp operation
+							for applicable_rule in applicable_node_rules:
+								n, annotations, qualified_nodes, qualified_edges, edges_to_add = applicable_rule
+								# If there is an edge to add or the predicate doesn't exist or the interpretation is not static
+								if len(edges_to_add[0]) > 0 or rule.get_target() not in interpretations_node[n].world or not interpretations_node[n].world[rule.get_target()].is_static():
+									bnd = annotate(annotation_functions, rule, annotations, rule.get_weights())
+									bnd = interval.closed(bnd[0], bnd[1])
+									max_rules_time = max(max_rules_time, t + delta_t)
+									edges_to_be_added_node_rule.append(edges_to_add)
+									rules_to_be_applied_node.append((numba.types.uint16(t + delta_t), n, rule.get_target(), bnd, immediate_rule, rule.is_static_rule()))
+									if atom_trace:
+										rules_to_be_applied_node_trace.append((qualified_nodes, qualified_edges, rule.get_name()))
+
+									# We apply a rule on a node/edge only once in each timestep to prevent it from being added to the to_be_added list continuously (this will improve performance
+									nodes_to_skip[i].append(n)
+
+									# Handle loop parameters for the next (maybe) fp operation
+									# If it is a t=0 rule or an immediate rule we want to go back for another fp operation to check for new rules that may fire
+									# Next fp operation we will skip this rule on this node because anyway there won't be an update
+									if delta_t == 0:
+										in_loop = True
+										update = False
+									if immediate_rule and delta_t == 0:
+										# immediate_rule_fire becomes True because we still need to check for more eligible rules, we're not done.
+										in_loop = True
+										update = True
+										immediate_node_rule_fire = True
+										break
 
 							# Break, apply immediate rule then come back to check for more applicable rules
 							if immediate_node_rule_fire:
 								break
 
-							# Go through all edges and check if any rules apply to them.
-							# Comment out the following lines if there are no labels or rules that deal with edges. It will be an unnecessary loop
-							for j in range(len(edges)):
-								if (i, j) in edges_to_skip:
-									continue
-								e = edges[j]
-								target_criteria_satisfaction = are_satisfied_edge(interpretations_edge, e, rule.get_target_criteria())
-								if (target_criteria_satisfaction and is_satisfied_node(interpretations_edge, e, (rule.get_target(), interval.closed(0,1)))) or (target_criteria_satisfaction and rule.get_target().value==''):
-									# Find out if rule is applicable. returns list of list, of qualified nodes and qualified edges. one for each clause
-									result, annotations, qualified_nodes, qualified_edges, edges_to_add = _is_rule_applicable_edge(interpretations_node, interpretations_edge, neighbors, e, rule.get_neigh_criteria(), rule.get_thresholds(), reverse_graph, rule.get_annotation_function(), rule.get_annotation_label(), rule.get_edges(), atom_trace)
-									if (result and rule.get_target().value=='') or (result and not interpretations_edge[e].world[rule.get_target()].is_static()):
-										bnd = influence(rule, annotations, rule.get_weights())
-										max_rules_time = max(max_rules_time, t+delta_t)
-										edges_to_be_added_edge_rule.append(edges_to_add)
-										rules_to_be_applied_edge.append((numba.types.uint16(t+delta_t), e, rule.get_target(), bnd, immediate_rule))
-										if atom_trace:
-											rules_to_be_applied_edge_trace.append((qualified_nodes, qualified_edges, rule.get_name()))
-
-										# We apply a rule on a node/edge only once in each timestep to prevent it from being added to the to_be_added list continuously (this will improve performance
-										edges_to_skip.append((i, j))
-
-										# Handle loop parameters for the next (maybe) fp operation
-										# If it is a t=0 rule or an immediate rule we want to go back for another fp operation to check for new rules that may fire
-										# Next fp operation we will skip this rule on this node because anyway there won't be an update
-										if delta_t==0:
-											in_loop = True
-											update = False
-										if immediate_rule and delta_t == 0:
-											# immediate_rule_fire becomes True because we still need to check for more eligible rules, we're not done.
-											in_loop = True
-											update = True
-											immediate_edge_rule_fire = True
-											break
+							for applicable_rule in applicable_edge_rules:
+								e, annotations, qualified_nodes, qualified_edges, edges_to_add = applicable_rule
+								# If there is an edge to add or the predicate doesn't exist or the interpretation is not static
+								if len(edges_to_add[0]) > 0 or rule.get_target() not in interpretations_edge[e].world or not interpretations_edge[e].world[rule.get_target()].is_static():
+									bnd = annotate(annotation_functions, rule, annotations, rule.get_weights())
+									bnd = interval.closed(bnd[0], bnd[1])
+									max_rules_time = max(max_rules_time, t+delta_t)
+									edges_to_be_added_edge_rule.append(edges_to_add)
+									rules_to_be_applied_edge.append((numba.types.uint16(t+delta_t), e, rule.get_target(), bnd, immediate_rule, rule.is_static_rule()))
+									if atom_trace:
+										rules_to_be_applied_edge_trace.append((qualified_nodes, qualified_edges, rule.get_name()))
+
+									# We apply a rule on a node/edge only once in each timestep to prevent it from being added to the to_be_added list continuously (this will improve performance
+									edges_to_skip[i].append(e)
+
+									# Handle loop parameters for the next (maybe) fp operation
+									# If it is a t=0 rule or an immediate rule we want to go back for another fp operation to check for new rules that may fire
+									# Next fp operation we will skip this rule on this node because anyway there won't be an update
+									if delta_t == 0:
+										in_loop = True
+										update = False
+									if immediate_rule and delta_t == 0:
+										# immediate_rule_fire becomes True because we still need to check for more eligible rules, we're not done.
+										in_loop = True
+										update = True
+										immediate_edge_rule_fire = True
+										break
 
 							# Break, apply immediate rule then come back to check for more applicable rules
 							if immediate_edge_rule_fire:
@@ -637,255 +670,666 @@ def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data
 			# Increment t
 			t += 1
 
-		return fp_cnt, t	
+		return fp_cnt, t
+
+	def add_edge(self, edge, l):
+		# This function is useful for pyreason gym, called externally
+		_add_edge(edge[0], edge[1], self.neighbors, self.reverse_neighbors, self.nodes, self.edges, l, self.interpretations_node, self.interpretations_edge)
+
+	def add_node(self, node, labels):
+		# This function is useful for pyreason gym, called externally
+		if node not in self.nodes:
+			_add_node(node, self.neighbors, self.reverse_neighbors, self.nodes, self.interpretations_node)
+			for l in labels:
+				self.interpretations_node[node].world[label.Label(l)] = interval.closed(0, 1)
+
+	def delete_edge(self, edge):
+		# This function is useful for pyreason gym, called externally
+		_delete_edge(edge, self.neighbors, self.reverse_neighbors, self.edges, self.interpretations_edge)
+
+	def get_interpretation_dict(self):
+		# This function can be called externally to retrieve a dict of the interpretation values
+		# Only values in the rule trace will be added
+
+		# Initialize interpretations for each time and node and edge
+		interpretations = {}
+		for t in range(self.tmax+1):
+			interpretations[t] = {}
+			for node in self.nodes:
+				interpretations[t][node] = InterpretationDict()
+			for edge in self.edges:
+				interpretations[t][edge] = InterpretationDict()
+
+		# Update interpretation nodes
+		for change in self.rule_trace_node:
+			time, _, node, l, bnd = change
+			interpretations[time][node][l._value] = (bnd.lower, bnd.upper)
+
+			# If canonical, update all following timesteps as well
+			if self. canonical:
+				for t in range(time+1, self.tmax+1):
+					interpretations[t][node][l._value] = (bnd.lower, bnd.upper)
+
+		# Update interpretation edges
+		for change in self.rule_trace_edge:
+			time, _, edge, l, bnd, = change
+			interpretations[time][edge][l._value] = (bnd.lower, bnd.upper)
+
+			# If canonical, update all following timesteps as well
+			if self. canonical:
+				for t in range(time+1, self.tmax+1):
+					interpretations[t][edge][l._value] = (bnd.lower, bnd.upper)
+
+		return interpretations
 
 
 @numba.njit(cache=True)
-def _is_rule_applicable(interpretations_node, interpretations_edge, neighbors, target_node, neigh_criteria, thresholds, reverse_graph, ann_fn, ann_fn_label, edges, atom_trace):
-	# Initialize dictionary where keys are strings (x1, x2 etc.) and values are lists of qualified neighbors
-	# Keep track of all the edges that are qualified
-	# If it's a node clause update (x1 or x2 etc.) qualified neighbors, if it's an edge clause update the qualified neighbors for the source and target (x1, x2)
-	# First gather all the qualified nodes for each clause
-	subsets = numba.typed.Dict.empty(key_type=numba.types.string, value_type=list_of_nodes)
-	qualified_nodes = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
-	qualified_edges = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
-	annotations = numba.typed.List.empty_list(numba.typed.List.empty_list(interval.interval_type))
-	edges_to_be_added = (numba.typed.List.empty_list(node_type), numba.typed.List.empty_list(node_type), edges[-1])
+def _ground_node_rule(rule, interpretations_node, interpretations_edge, nodes, neighbors, reverse_neighbors, atom_trace, reverse_graph, nodes_to_skip):
+	# Extract rule params
+	rule_type = rule.get_type()
+	clauses = rule.get_clauses()
+	thresholds = rule.get_thresholds()
+	ann_fn = rule.get_annotation_function()
+	rule_edges = rule.get_edges()
+
+	# We return a list of tuples which specify the target nodes/edges that have made the rule body true
+	applicable_rules = numba.typed.List.empty_list(node_applicable_rule_type)
+
+	# Return empty list if rule is not node rule and if we are not inferring edges
+	if rule_type != 'node' and rule_edges[0] == '':
+		return applicable_rules
 
 	# Steps
-	# 1. Gather all qualified nodes/edges, and check if they satisfy the thresholds. If not, break
-	# 2. Gather all annotations associated with the qualified nodes/edges
-	# 3. Assemble all constants for atom trace (this happens inside step 1 for efficiency)
-	# 4. Collect edges and nodes to be added to the graph when rule fires
-
-	# 1.
-	satisfaction = True
-	for i, clause in enumerate(neigh_criteria):
-		# Gather qualified nodes/edges
-		if clause[0]=='node':
-			if clause[1][0]=='target':
-				subset = numba.typed.List([target_node])
-			else:
-				subset = neighbors[target_node] if clause[1][0] not in subsets else subsets[clause[1][0]]
+	# 1. Loop through all nodes and evaluate each clause with that node and check the truth with the thresholds
+	# 2. Inside the clause loop it may be necessary to loop through all nodes/edges while grounding the variables
+	# 3. If the clause is true add the qualified nodes and qualified edges to the atom trace, if on. Break otherwise
+	# 4. After going through all clauses, add to the annotations list all the annotations of the specified subset. These will be passed to the annotation function
+	# 5. Finally, if there are any edges to be added, place them in the list
+
+	for piter in prange(len(nodes)):
+		target_node = nodes[piter]
+		if target_node in nodes_to_skip:
+			continue
+		# Initialize dictionary where keys are strings (x1, x2 etc.) and values are lists of qualified neighbors
+		# Keep track of qualified nodes and qualified edges
+		# If it's a node clause update (x1 or x2 etc.) qualified neighbors, if it's an edge clause update the qualified neighbors for the source and target (x1, x2)
+		subsets = numba.typed.Dict.empty(key_type=numba.types.string, value_type=list_of_nodes)
+		qualified_nodes = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+		qualified_edges = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
+		annotations = numba.typed.List.empty_list(numba.typed.List.empty_list(interval.interval_type))
+		edges_to_be_added = (numba.typed.List.empty_list(node_type), numba.typed.List.empty_list(node_type), rule_edges[-1])
+
+		satisfaction = True
+		for i, clause in enumerate(clauses):
+			# Unpack clause variables
+			clause_type = clause[0]
+			clause_label = clause[1]
+			clause_variables = clause[2]
+			clause_bnd = clause[3]
+			clause_operator = clause[4]
+
+			# Unpack thresholds
+			# This value is total/available
+			threshold_quantifier_type = thresholds[i][1][1]
+
+			# This is a node clause
+			# The groundings for node clauses are either the target node, neighbors of the target node, or an existing subset of nodes
+			if clause_type == 'node':
+				clause_var_1 = clause_variables[0]
+				subset = get_node_rule_node_clause_subset(clause_var_1, target_node, subsets, neighbors)
+
+				subsets[clause_var_1] = get_qualified_components_node_clause(interpretations_node, subset, clause_label, clause_bnd)
 
-			subsets[clause[1][0]] = get_qualified_components_node_clause(interpretations_node, subset, clause[2], clause[3])
-			if atom_trace:
-				qualified_nodes.append(numba.typed.List(subsets[clause[1][0]]))
-				qualified_edges.append(numba.typed.List.empty_list(edge_type))
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List(subsets[clause_var_1]))
+					qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qn in subsets[clause_var_1]:
+						a.append(interpretations_node[qn].world[clause_label])
+					annotations.append(a)
+
+			# This is an edge clause
+			elif clause_type == 'edge':
+				clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+				subset_source, subset_target = get_node_rule_edge_clause_subset(clause_var_1, clause_var_2, target_node, subsets, neighbors, reverse_neighbors, nodes)
+
+				# Get qualified edges
+				qe = get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, clause_bnd, reverse_graph)
+				subsets[clause_var_1] = qe[0]
+				subsets[clause_var_2] = qe[1]
 
-		elif clause[0]=='edge':
-			# Set sources for possible edges, if target use target node as source
-			if clause[1][0]=='target':
-				subset_source = numba.typed.List([target_node])
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List.empty_list(node_type))
+					qualified_edges.append(numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qe in numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])):
+						a.append(interpretations_edge[qe].world[clause_label])
+					annotations.append(a)
 			else:
-				subset_source = neighbors[target_node] if clause[1][0] not in subsets else subsets[clause[1][0]]
+				# This is a comparison clause
+				# Make sure there is at least one ground atom such that pred-num(x) : [1,1] or pred-num(x,y) : [1,1]
+				# Remember that the predicate in the clause will not contain the "-num" where num is some number.
+				# We have to remove that manually while checking
+				# Steps:
+				# 1. get qualified nodes/edges as well as number associated for first predicate
+				# 2. get qualified nodes/edges as well as number associated for second predicate
+				# 3. if there's no number in steps 1 or 2 return false clause
+				# 4. do comparison with each qualified component from step 1 with each qualified component in step 2
+
+				# It's a node comparison
+				if len(clause_variables) == 2:
+					clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+					subset_1 = get_node_rule_node_clause_subset(clause_var_1, target_node, subsets, neighbors)
+					subset_2 = get_node_rule_node_clause_subset(clause_var_2, target_node, subsets, neighbors)
+
+					# 1, 2
+					qualified_nodes_for_comparison_1, numbers_1 = get_qualified_components_node_comparison_clause(interpretations_node, subset_1, clause_label, clause_bnd)
+					qualified_nodes_for_comparison_2, numbers_2 = get_qualified_components_node_comparison_clause(interpretations_node, subset_2, clause_label, clause_bnd)
+
+				# It's an edge comparison
+				elif len(clause_variables) == 4:
+					clause_var_1_source, clause_var_1_target, clause_var_2_source, clause_var_2_target = clause_variables[0], clause_variables[1], clause_variables[2], clause_variables[3]
+					subset_1_source, subset_1_target = get_node_rule_edge_clause_subset(clause_var_1_source, clause_var_1_target, target_node, subsets, neighbors, reverse_neighbors, nodes)
+					subset_2_source, subset_2_target = get_node_rule_edge_clause_subset(clause_var_2_source, clause_var_2_target, target_node, subsets, neighbors, reverse_neighbors, nodes)
+
+					# 1, 2
+					qualified_nodes_for_comparison_1_source, qualified_nodes_for_comparison_1_target, numbers_1 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_1_source, subset_1_target, clause_label, clause_bnd, reverse_graph)
+					qualified_nodes_for_comparison_2_source, qualified_nodes_for_comparison_2_target, numbers_2 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_2_source, subset_2_target, clause_label, clause_bnd, reverse_graph)
+
+			# Check if thresholds are satisfied
+			# If it's a comparison clause we just need to check if the numbers list is not empty (no threshold support)
+			if clause_type == 'comparison':
+				if len(numbers_1) == 0 or len(numbers_2) == 0:
+					satisfaction = False
+				# Node comparison. Compare stage
+				elif len(clause_variables) == 2:
+					satisfaction, qualified_nodes_1, qualified_nodes_2 = compare_numbers_node_predicate(numbers_1, numbers_2, clause_operator, qualified_nodes_for_comparison_1, qualified_nodes_for_comparison_2)
+
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1] = qualified_nodes_1
+					subsets[clause_var_2] = qualified_nodes_2
+					qualified_comparison_nodes = numba.typed.List(qualified_nodes_1)
+					qualified_comparison_nodes.extend(qualified_nodes_2)
 
-			subset_target = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
-			if clause[1][1]=='target':
-				# Set targets for possible edges, if target use target node as target
-				for source in subset_source:
-					subset_target.append(numba.typed.List([target_node]))
-			else:
-				for source in subset_source:
-					subset_target.append(neighbors[source] if clause[1][1] not in subsets else subsets[clause[1][1]])
+					if atom_trace:
+						qualified_nodes.append(qualified_comparison_nodes)
+						qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qn in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
+				# Edge comparison. Compare stage
+				else:
+					satisfaction, qualified_nodes_1_source, qualified_nodes_1_target, qualified_nodes_2_source, qualified_nodes_2_target = compare_numbers_edge_predicate(numbers_1, numbers_2, clause_operator,
+																																										  qualified_nodes_for_comparison_1_source,
+																																										  qualified_nodes_for_comparison_1_target,
+																																										  qualified_nodes_for_comparison_2_source,
+																																										  qualified_nodes_for_comparison_2_target)
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1_source] = qualified_nodes_1_source
+					subsets[clause_var_1_target] = qualified_nodes_1_target
+					subsets[clause_var_2_source] = qualified_nodes_2_source
+					subsets[clause_var_2_target] = qualified_nodes_2_target
+
+					qualified_comparison_nodes_1 = numba.typed.List(zip(qualified_nodes_1_source, qualified_nodes_1_target))
+					qualified_comparison_nodes_2 = numba.typed.List(zip(qualified_nodes_2_source, qualified_nodes_2_target))
+					qualified_comparison_nodes = numba.typed.List(qualified_comparison_nodes_1)
+					qualified_comparison_nodes.extend(qualified_comparison_nodes_2)
 
-			qe = get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause[2], clause[3], reverse_graph)
-			subsets[clause[1][0]] = qe[0]
-			subsets[clause[1][1]] = qe[1]
+					if atom_trace:
+						qualified_nodes.append(numba.typed.List.empty_list(node_type))
+						qualified_edges.append(qualified_comparison_nodes)
 
-			if atom_trace:
-				qualified_nodes.append(numba.typed.List.empty_list(node_type))
-				qualified_edges.append(numba.typed.List(zip(subsets[clause[1][0]], subsets[clause[1][1]])))
-
-		# Check if the clause satisfies threshold
-		if thresholds[i][1][1]=='total':
-			if clause[0]=='node':
-				neigh_len = len(subset)
-			elif clause[0]=='edge':
-				neigh_len = sum([len(l) for l in subset_target])
-
-		# Available is all neighbors that have a particular label with bound inside [0,1]
-		elif thresholds[i][1][1]=='available':
-			if clause[0]=='node':
-				neigh_len = len(get_qualified_components_node_clause(interpretations_node, subset, clause[2], interval.closed(0,1)))
-			if clause[0]=='edge':
-				neigh_len = len(get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause[2], interval.closed(0,1), reverse_graph)[0])
-
-		qualified_neigh_len = len(subsets[clause[1][0]])
-		satisfaction = _satisfies_threshold(neigh_len, qualified_neigh_len, thresholds[i]) and satisfaction
-		# Exit loop if even one clause is not satisfied
-		if satisfaction==False:
-			break
-
-	# 2.
-	# Add to the subsets that will be used in annotation function
-	if satisfaction:
-		for clause in neigh_criteria:
-			if ann_fn!='':
-				a = numba.typed.List.empty_list(interval.interval_type)
-				if clause[0]=='node':
-					for qn in subsets[clause[1][0]]:
-						a.append(interpretations_node[qn].world[ann_fn_label])
-				elif clause[0]=='edge':
-					for qe in numba.typed.List(zip(subsets[clause[1][0]], subsets[clause[1][1]])):
-						a.append(interpretations_edge[qe].world[ann_fn_label])
-
-				annotations.append(a)
-
-	# 4.
-	# Collect edges to be added
-	if satisfaction:
-		source, target, _ = edges
-
-		# Edges to be added
-		if source!='' and target!='':
-			# Check if edge nodes are target
-			if source=='target':
-				edges_to_be_added[0].append(target_node)
-			elif source in subsets:
-				edges_to_be_added[0].extend(subsets[source])
-			else:
-				edges_to_be_added[0].append(source)
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qe in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
 
-			if target=='target':
-				edges_to_be_added[1].append(target_node)			
-			elif target in subsets:
-				edges_to_be_added[1].extend(subsets[target])
+			# Non comparison clause
 			else:
-				edges_to_be_added[1].append(target)
+				if threshold_quantifier_type == 'total':
+					if clause_type == 'node':
+						neigh_len = len(subset)
+					else:
+						neigh_len = sum([len(l) for l in subset_target])
+
+				# Available is all neighbors that have a particular label with bound inside [0,1]
+				elif threshold_quantifier_type == 'available':
+					if clause_type == 'node':
+						neigh_len = len(get_qualified_components_node_clause(interpretations_node, subset, clause_label, interval.closed(0,1)))
+					else:
+						neigh_len = len(get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, interval.closed(0,1), reverse_graph)[0])
+
+				qualified_neigh_len = len(subsets[clause_var_1])
+				satisfaction = _satisfies_threshold(neigh_len, qualified_neigh_len, thresholds[i]) and satisfaction
+
+			# Exit loop if even one clause is not satisfied
+			if not satisfaction:
+				break
+
+		if satisfaction:
+			# Collect edges to be added
+			source, target, _ = rule_edges
+
+			# Edges to be added
+			if source != '' and target != '':
+				# Check if edge nodes are target
+				if source == '__target':
+					edges_to_be_added[0].append(target_node)
+				elif source in subsets:
+					edges_to_be_added[0].extend(subsets[source])
+				else:
+					edges_to_be_added[0].append(source)
+
+				if target == '__target':
+					edges_to_be_added[1].append(target_node)
+				elif target in subsets:
+					edges_to_be_added[1].extend(subsets[target])
+				else:
+					edges_to_be_added[1].append(target)
 
-	return (satisfaction, annotations, qualified_nodes, qualified_edges, edges_to_be_added)
+			# node/edge, annotations, qualified nodes, qualified edges, edges to be added
+			applicable_rules.append((target_node, annotations, qualified_nodes, qualified_edges, edges_to_be_added))
+
+	return applicable_rules
 
 
 @numba.njit(cache=True)
-def _is_rule_applicable_edge(interpretations_node, interpretations_edge, neighbors, target_edge, neigh_criteria, thresholds, reverse_graph, ann_fn, ann_fn_label, edges, atom_trace):
-	# Initialize dictionary where keys are strings (x1, x2 etc.) and values are lists of qualified neighbors
-	# Keep track of all the edges that are qualified
-	# If it's a node clause update (x1 or x2 etc.) qualified neighbors, if it's an edge clause update the qualified neighbors for the source and target (x1, x2)
-	# First gather all the qualified nodes for each clause
-	subsets = numba.typed.Dict.empty(key_type=numba.types.string, value_type=list_of_nodes)
-	qualified_nodes = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
-	qualified_edges = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
-	annotations = numba.typed.List.empty_list(numba.typed.List.empty_list(interval.interval_type))
-	edges_to_be_added = (numba.typed.List.empty_list(node_type), numba.typed.List.empty_list(node_type), edges[-1])
+def _ground_edge_rule(rule, interpretations_node, interpretations_edge, nodes, edges, neighbors, reverse_neighbors, atom_trace, reverse_graph, edges_to_skip):
+	# Extract rule params
+	rule_type = rule.get_type()
+	clauses = rule.get_clauses()
+	thresholds = rule.get_thresholds()
+	ann_fn = rule.get_annotation_function()
+	rule_edges = rule.get_edges()
+
+	# We return a list of tuples which specify the target nodes/edges that have made the rule body true
+	applicable_rules = numba.typed.List.empty_list(edge_applicable_rule_type)
+
+	# Return empty list if rule is not node rule
+	if rule_type != 'edge':
+		return applicable_rules
 
 	# Steps
-	# 1. Gather all qualified nodes/edges, and check if they satisfy the thresholds. If not, break
-	# 2. Gather all annotations associated with the qualified nodes/edges
-	# 3. Assemble all constants for atom trace (this happens inside step 1 for efficiency)
-	# 4. Collect edges and nodes to be added to the graph when rule fires
-
-	# 1.
-	satisfaction = True
-	for i, clause in enumerate(neigh_criteria):
-		# Gather qualified nodes/edges
-		# For nodes gather source neighbors only, not target neighbors
-		if clause[0]=='node':
-			if clause[1][0]=='source':
-				subset = numba.typed.List([target_edge[0]])
-			elif clause[1][0]=='target':
-				subset = numba.typed.List([target_edge[1]])
-			else:
-				subset = neighbors[target_edge[0]] if clause[1][0] not in subsets else subsets[clause[1][0]]
+	# 1. Loop through all nodes and evaluate each clause with that node and check the truth with the thresholds
+	# 2. Inside the clause loop it may be necessary to loop through all nodes/edges while grounding the variables
+	# 3. If the clause is true add the qualified nodes and qualified edges to the atom trace, if on. Break otherwise
+	# 4. After going through all clauses, add to the annotations list all the annotations of the specified subset. These will be passed to the annotation function
+	# 5. Finally, if there are any edges to be added, place them in the list
+
+	for piter in prange(len(edges)):
+		target_edge = edges[piter]
+		if target_edge in edges_to_skip:
+			continue
+		# Initialize dictionary where keys are strings (x1, x2 etc.) and values are lists of qualified neighbors
+		# Keep track of qualified nodes and qualified edges
+		# If it's a node clause update (x1 or x2 etc.) qualified neighbors, if it's an edge clause update the qualified neighbors for the source and target (x1, x2)
+		subsets = numba.typed.Dict.empty(key_type=numba.types.string, value_type=list_of_nodes)
+		qualified_nodes = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+		qualified_edges = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
+		annotations = numba.typed.List.empty_list(numba.typed.List.empty_list(interval.interval_type))
+		edges_to_be_added = (numba.typed.List.empty_list(node_type), numba.typed.List.empty_list(node_type), rule_edges[-1])
+
+		satisfaction = True
+		for i, clause in enumerate(clauses):
+			# Unpack clause variables
+			clause_type = clause[0]
+			clause_label = clause[1]
+			clause_variables = clause[2]
+			clause_bnd = clause[3]
+			clause_operator = clause[4]
+
+			# Unpack thresholds
+			# This value is total/available
+			threshold_quantifier_type = thresholds[i][1][1]
+
+			# This is a node clause
+			# The groundings for node clauses are either the source, target, neighbors of the source node, or an existing subset of nodes
+			if clause_type == 'node':
+				clause_var_1 = clause_variables[0]
+				subset = get_edge_rule_node_clause_subset(clause_var_1, target_edge, subsets, neighbors)
+				
+				subsets[clause_var_1] = get_qualified_components_node_clause(interpretations_node, subset, clause_label, clause_bnd)
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List(subsets[clause_var_1]))
+					qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qn in subsets[clause_var_1]:
+						a.append(interpretations_node[qn].world[clause_label])
+					annotations.append(a)
+
+			# This is an edge clause
+			elif clause_type == 'edge':
+				clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+				subset_source, subset_target = get_edge_rule_edge_clause_subset(clause_var_1, clause_var_2, target_edge, subsets, neighbors, reverse_neighbors, nodes)
+
+				# Get qualified edges
+				qe = get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, clause_bnd, reverse_graph)
+				subsets[clause_var_1] = qe[0]
+				subsets[clause_var_2] = qe[1]
 
-			subsets[clause[1][0]] = get_qualified_components_node_clause(interpretations_node, subset, clause[2], clause[3])
-			if atom_trace:
-				qualified_nodes.append(numba.typed.List(subsets[clause[1][0]]))
-				qualified_edges.append(numba.typed.List.empty_list(edge_type))
-
-		elif clause[0]=='edge':
-			# Set sources for possible edges, if target use target node as source
-			# By default we take the subset_source to be the neighbors of the source node.
-			if clause[1][0]=='source':
-				subset_source = numba.typed.List([target_edge[0]])
-			elif clause[1][0]=='target':
-				subset_source = numba.typed.List([target_edge[1]])
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List.empty_list(node_type))
+					qualified_edges.append(numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qe in numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])):
+						a.append(interpretations_edge[qe].world[clause_label])
+					annotations.append(a)
+					
 			else:
-				subset_source = neighbors[target_edge[0]] if clause[1][0] not in subsets else subsets[clause[1][0]]
+				# This is a comparison clause
+				# Make sure there is at least one ground atom such that pred-num(x) : [1,1] or pred-num(x,y) : [1,1]
+				# Remember that the predicate in the clause will not contain the "-num" where num is some number.
+				# We have to remove that manually while checking
+				# Steps:
+				# 1. get qualified nodes/edges as well as number associated for first predicate
+				# 2. get qualified nodes/edges as well as number associated for second predicate
+				# 3. if there's no number in steps 1 or 2 return false clause
+				# 4. do comparison with each qualified component from step 1 with each qualified component in step 2
 				
-			subset_target = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
-			if clause[1][1]=='source':
-				# Set targets for possible edges, if source use source node as target
-				for source in subset_source:
-					subset_target.append(numba.typed.List([target_edge[0]]))
-			elif clause[1][1]=='target':
-				# Set targets for possible edges, if target use target node as target
-				for source in subset_source:
-					subset_target.append(numba.typed.List([target_edge[1]]))
-			else:
-				for source in subset_source:
-					subset_target.append(neighbors[source] if clause[1][1] not in subsets else subsets[clause[1][1]])
+				# It's a node comparison
+				if len(clause_variables) == 2:
+					clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+					subset_1 = get_edge_rule_node_clause_subset(clause_var_1, target_edge, subsets, neighbors)
+					subset_2 = get_edge_rule_node_clause_subset(clause_var_2, target_edge, subsets, neighbors)
+					
+					# 1, 2
+					qualified_nodes_for_comparison_1, numbers_1 = get_qualified_components_node_comparison_clause(interpretations_node, subset_1, clause_label, clause_bnd)
+					qualified_nodes_for_comparison_2, numbers_2 = get_qualified_components_node_comparison_clause(interpretations_node, subset_2, clause_label, clause_bnd)
+				
+				# It's an edge comparison
+				elif len(clause_variables) == 4:
+					clause_var_1_source, clause_var_1_target, clause_var_2_source, clause_var_2_target = clause_variables[0], clause_variables[1], clause_variables[2], clause_variables[3]
+					subset_1_source, subset_1_target = get_edge_rule_edge_clause_subset(clause_var_1_source, clause_var_1_target, target_edge, subsets, neighbors, reverse_neighbors, nodes)
+					subset_2_source, subset_2_target = get_edge_rule_edge_clause_subset(clause_var_2_source, clause_var_2_target, target_edge, subsets, neighbors, reverse_neighbors, nodes)
+
+					# 1, 2
+					qualified_nodes_for_comparison_1_source, qualified_nodes_for_comparison_1_target, numbers_1 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_1_source, subset_1_target, clause_label, clause_bnd, reverse_graph)
+					qualified_nodes_for_comparison_2_source, qualified_nodes_for_comparison_2_target, numbers_2 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_2_source, subset_2_target, clause_label, clause_bnd, reverse_graph)
+
+			# Check if thresholds are satisfied
+			# If it's a comparison clause we just need to check if the numbers list is not empty (no threshold support)
+			if clause_type == 'comparison':
+				if len(numbers_1) == 0 or len(numbers_2) == 0:
+					satisfaction = False
+				# Node comparison. Compare stage
+				elif len(clause_variables) == 2:
+					satisfaction, qualified_nodes_1, qualified_nodes_2 = compare_numbers_node_predicate(numbers_1, numbers_2, clause_operator, qualified_nodes_for_comparison_1, qualified_nodes_for_comparison_2)
+
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1] = qualified_nodes_1
+					subsets[clause_var_2] = qualified_nodes_2
+					qualified_comparison_nodes = numba.typed.List(qualified_nodes_1)
+					qualified_comparison_nodes.extend(qualified_nodes_2)
 
-			qe = get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause[2], clause[3], reverse_graph)
-			subsets[clause[1][0]] = qe[0]
-			subsets[clause[1][1]] = qe[1]
+					if atom_trace:
+						qualified_nodes.append(qualified_comparison_nodes)
+						qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qn in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
+				# Edge comparison. Compare stage
+				else:
+					satisfaction, qualified_nodes_1_source, qualified_nodes_1_target, qualified_nodes_2_source, qualified_nodes_2_target = compare_numbers_edge_predicate(numbers_1, numbers_2, clause_operator,
+																																										  qualified_nodes_for_comparison_1_source,
+																																										  qualified_nodes_for_comparison_1_target,
+																																										  qualified_nodes_for_comparison_2_source,
+																																										  qualified_nodes_for_comparison_2_target)
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1_source] = qualified_nodes_1_source
+					subsets[clause_var_1_target] = qualified_nodes_1_target
+					subsets[clause_var_2_source] = qualified_nodes_2_source
+					subsets[clause_var_2_target] = qualified_nodes_2_target
+
+					qualified_comparison_nodes_1 = numba.typed.List(zip(qualified_nodes_1_source, qualified_nodes_1_target))
+					qualified_comparison_nodes_2 = numba.typed.List(zip(qualified_nodes_2_source, qualified_nodes_2_target))
+					qualified_comparison_nodes = numba.typed.List(qualified_comparison_nodes_1)
+					qualified_comparison_nodes.extend(qualified_comparison_nodes_2)
 
-			if atom_trace:
-				qualified_nodes.append(numba.typed.List.empty_list(node_type))
-				qualified_edges.append(numba.typed.List(zip(subsets[clause[1][0]], subsets[clause[1][1]])))
-
-		# Check if the clause satisfies threshold
-		if thresholds[i][1][1]=='total':
-			if clause[0]=='node':
-				neigh_len = len(subset)
-			elif clause[0]=='edge':
-				neigh_len = sum([len(l) for l in subset_target])
-
-		# Available is all neighbors that have a particular label with bound inside [0,1]
-		elif thresholds[i][1][1]=='available':
-			if clause[0]=='node':
-				neigh_len = len(get_qualified_components_node_clause(interpretations_node, subset, clause[2], interval.closed(0,1)))
-			if clause[0]=='edge':
-				neigh_len = len(get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause[2], interval.closed(0,1), reverse_graph)[0])
-
-		qualified_neigh_len = len(subsets[clause[1][0]])
-		satisfaction = _satisfies_threshold(neigh_len, qualified_neigh_len, thresholds[i]) and satisfaction
-		# Exit loop if even one clause is not satisfied
-		if satisfaction==False:
-			break
-
-	# 2.
-	# Add to the subsets that will be used in annotation function
-	if satisfaction:
-		for clause in neigh_criteria:
-			if ann_fn!='':
-				a = numba.typed.List.empty_list(interval.interval_type)
-				if clause[0]=='node':
-					for qn in subsets[clause[1][0]]:
-						a.append(interpretations_node[qn].world[ann_fn_label])
-				elif clause[0]=='edge':
-					for qe in numba.typed.List(zip(subsets[clause[1][0]], subsets[clause[1][1]])):
-						a.append(interpretations_edge[qe].world[ann_fn_label])
-
-				annotations.append(a)
-
-	# 4.
-	# Collect edges to be added
-	if satisfaction:
-		source, target, _ = edges
-
-		# Edges to be added
-		if source!='' and target!='':
-			# Check if edge nodes are source or target
-			if source=='source':
-				edges_to_be_added[0].append(target_edge[0])
-			elif source=='target':
-				edges_to_be_added[0].append(target_edge[1])
-			elif source in subsets:
-				edges_to_be_added[0].extend(subsets[source])
-			else:
-				edges_to_be_added[0].append(source)
-
-			if target=='source':
-				edges_to_be_added[1].append(target_edge[0])
-			elif target=='target':
-				edges_to_be_added[1].append(target_edge[1])			
-			elif target in subsets:
-				edges_to_be_added[1].extend(subsets[target])
+					if atom_trace:
+						qualified_nodes.append(numba.typed.List.empty_list(node_type))
+						qualified_edges.append(qualified_comparison_nodes)
+
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qe in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
+			
+			# Non comparison clause
 			else:
-				edges_to_be_added[1].append(target)
+				if threshold_quantifier_type == 'total':
+					if clause_type == 'node':
+						neigh_len = len(subset)
+					else:
+						neigh_len = sum([len(l) for l in subset_target])
+	
+				# Available is all neighbors that have a particular label with bound inside [0,1]
+				elif threshold_quantifier_type == 'available':
+					if clause_type == 'node':
+						neigh_len = len(get_qualified_components_node_clause(interpretations_node, subset, clause_label, interval.closed(0, 1)))
+					else:
+						neigh_len = len(get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, interval.closed(0, 1), reverse_graph)[0])
+	
+				qualified_neigh_len = len(subsets[clause_var_1])
+				satisfaction = _satisfies_threshold(neigh_len, qualified_neigh_len, thresholds[i]) and satisfaction
+			
+			# Exit loop if even one clause is not satisfied
+			if not satisfaction:
+				break
+
+		# Here we are done going through each clause of the rule
+		# If all clauses we're satisfied, proceed to collect annotations and prepare edges to be added
+		if satisfaction:
+			# Collect edges to be added
+			source, target, _ = rule_edges
+
+			# Edges to be added
+			if source != '' and target != '':
+				# Check if edge nodes are source/target
+				if source == '__source':
+					edges_to_be_added[0].append(target_edge[0])
+				elif source == '__target':
+					edges_to_be_added[0].append(target_edge[1])
+				elif source in subsets:
+					edges_to_be_added[0].extend(subsets[source])
+				else:
+					edges_to_be_added[0].append(source)
+
+				if target == '__source':
+					edges_to_be_added[1].append(target_edge[0])
+				elif target == '__target':
+					edges_to_be_added[1].append(target_edge[1])
+				elif target in subsets:
+					edges_to_be_added[1].extend(subsets[target])
+				else:
+					edges_to_be_added[1].append(target)
+
+			# node/edge, annotations, qualified nodes, qualified edges, edges to be added
+			applicable_rules.append((target_edge, annotations, qualified_nodes, qualified_edges, edges_to_be_added))
+
+	return applicable_rules
 
-	return (satisfaction, annotations, qualified_nodes, qualified_edges, edges_to_be_added)
+
+@numba.njit(cache=True)
+def get_node_rule_node_clause_subset(clause_var_1, target_node, subsets, neighbors):
+	# The groundings for node clauses are either the target node, neighbors of the target node, or an existing subset of nodes
+	if clause_var_1 == '__target':
+		subset = numba.typed.List([target_node])
+	else:
+		subset = neighbors[target_node] if clause_var_1 not in subsets else subsets[clause_var_1]
+	return subset
+
+
+@numba.njit(cache=True)
+def get_node_rule_edge_clause_subset(clause_var_1, clause_var_2, target_node, subsets, neighbors, reverse_neighbors, nodes):
+	# There are 5 cases for predicate(Y,Z):
+	# 1. Either one or both of Y, Z are the target node
+	# 2. Both predicate variables Y and Z have not been encountered before
+	# 3. The source variable Y has not been encountered before but the target variable Z has
+	# 4. The target variable Z has not been encountered before but the source variable Y has
+	# 5. Both predicate variables Y and Z have been encountered before
+
+	# Case 1:
+	# Check if 1st variable or 1st and 2nd variables are the target
+	if clause_var_1 == '__target':
+		subset_source = numba.typed.List([target_node])
+
+		# If both variables are the same
+		if clause_var_2 == '__target':
+			subset_target = numba.typed.List([numba.typed.List([target_node])])
+		elif clause_var_2 in subsets:
+			subset_target = numba.typed.List([subsets[clause_var_2]])
+		else:
+			subset_target = numba.typed.List([neighbors[target_node]])
+
+	# Check if 2nd variable is the target (this means 1st variable isn't the target)
+	elif clause_var_2 == '__target':
+		subset_source = reverse_neighbors[target_node] if clause_var_1 not in subsets else subsets[clause_var_1]
+		subset_target = numba.typed.List([numba.typed.List([target_node]) for _ in subset_source])
+
+	# Case 2:
+	# We replace Y by all nodes and Z by the neighbors of each of these nodes
+	elif clause_var_1 not in subsets and clause_var_2 not in subsets:
+		subset_source = numba.typed.List(nodes)
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 3:
+	# We replace Y by the sources of Z
+	elif clause_var_1 not in subsets and clause_var_2 in subsets:
+		subset_source = numba.typed.List.empty_list(node_type)
+		subset_target = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+
+		for n in subsets[clause_var_2]:
+			sources = reverse_neighbors[n]
+			for source in sources:
+				subset_source.append(source)
+				subset_target.append(numba.typed.List([n]))
+
+	# Case 4:
+	# We replace Z by the neighbors of Y
+	elif clause_var_1 in subsets and clause_var_2 not in subsets:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 5:
+	else:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([subsets[clause_var_2] for _ in subset_source])
+
+	return subset_source, subset_target
+
+
+@numba.njit(cache=True)
+def get_edge_rule_node_clause_subset(clause_var_1, target_edge, subsets, neighbors):
+	# The groundings for node clauses are either the source, target, neighbors of the source node, or an existing subset of nodes
+	if clause_var_1 == '__source':
+		subset = numba.typed.List([target_edge[0]])
+	elif clause_var_1 == '__target':
+		subset = numba.typed.List([target_edge[1]])
+	else:
+		subset = neighbors[target_edge[0]] if clause_var_1 not in subsets else subsets[clause_var_1]
+	return subset
+
+
+@numba.njit(cache=True)
+def get_edge_rule_edge_clause_subset(clause_var_1, clause_var_2, target_edge, subsets, neighbors, reverse_neighbors, nodes):
+	# There are 5 cases for predicate(Y,Z):
+	# 1. Either one or both of Y, Z are the source or target node
+	# 2. Both predicate variables Y and Z have not been encountered before
+	# 3. The source variable Y has not been encountered before but the target variable Z has
+	# 4. The target variable Z has not been encountered before but the source variable Y has
+	# 5. Both predicate variables Y and Z have been encountered before
+	# Case 1:
+	# Check if 1st variable is the source
+	if clause_var_1 == '__source':
+		subset_source = numba.typed.List([target_edge[0]])
+
+		# If 2nd variable is source/target/something else
+		if clause_var_2 == '__source':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[0]])])
+		elif clause_var_2 == '__target':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[1]])])
+		elif clause_var_2 in subsets:
+			subset_target = numba.typed.List([subsets[clause_var_2]])
+		else:
+			subset_target = numba.typed.List([neighbors[target_edge[0]]])
+
+	# if 1st variable is the target
+	elif clause_var_1 == '__target':
+		subset_source = numba.typed.List([target_edge[1]])
+
+		# if 2nd variable is source/target/something else
+		if clause_var_2 == '__source':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[0]])])
+		elif clause_var_2 == '__target':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[1]])])
+		elif clause_var_2 in subsets:
+			subset_target = numba.typed.List([subsets[clause_var_2]])
+		else:
+			subset_target = numba.typed.List([neighbors[target_edge[1]]])
+
+	# Handle the cases where the 2nd variable is source/target but the 1st is something else (cannot be source/target)
+	elif clause_var_2 == '__source':
+		subset_source = reverse_neighbors[target_edge[0]] if clause_var_1 not in subsets else subsets[clause_var_1]
+		subset_target = numba.typed.List([numba.typed.List([target_edge[0]]) for _ in subset_source])
+
+	elif clause_var_2 == '__target':
+		subset_source = reverse_neighbors[target_edge[1]] if clause_var_1 not in subsets else subsets[clause_var_1]
+		subset_target = numba.typed.List([numba.typed.List([target_edge[1]]) for _ in subset_source])
+
+	# Case 2:
+	# We replace Y by all nodes and Z by the neighbors of each of these nodes
+	elif clause_var_1 not in subsets and clause_var_2 not in subsets:
+		subset_source = numba.typed.List(nodes)
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 3:
+	# We replace Y by the sources of Z
+	elif clause_var_1 not in subsets and clause_var_2 in subsets:
+		subset_source = numba.typed.List.empty_list(node_type)
+		subset_target = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+
+		for n in subsets[clause_var_2]:
+			sources = reverse_neighbors[n]
+			for source in sources:
+				subset_source.append(source)
+				subset_target.append(numba.typed.List([n]))
+
+	# Case 4:
+	# We replace Z by the neighbors of Y
+	elif clause_var_1 in subsets and clause_var_2 not in subsets:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 5:
+	else:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([subsets[clause_var_2] for _ in subset_source])
+
+	return subset_source, subset_target
 
 
 @numba.njit(cache=True)
@@ -899,6 +1343,20 @@ def get_qualified_components_node_clause(interpretations_node, candidates, l, bn
 	return qualified_nodes
 
 
+@numba.njit(cache=True)
+def get_qualified_components_node_comparison_clause(interpretations_node, candidates, l, bnd):
+	# Get all the qualified neighbors for a particular comparison clause and return them along with the number associated
+	qualified_nodes = numba.typed.List.empty_list(node_type)
+	qualified_nodes_numbers = numba.typed.List.empty_list(numba.types.float64)
+	for n in candidates:
+		result, number = is_satisfied_node_comparison(interpretations_node, n, (l, bnd))
+		if result:
+			qualified_nodes.append(n)
+			qualified_nodes_numbers.append(number)
+
+	return qualified_nodes, qualified_nodes_numbers
+
+
 @numba.njit(cache=True)
 def get_qualified_components_edge_clause(interpretations_edge, candidates_source, candidates_target, l, bnd, reverse_graph):
 	# Get all the qualified sources and targets for a particular clause
@@ -911,8 +1369,94 @@ def get_qualified_components_edge_clause(interpretations_edge, candidates_source
 				qualified_nodes_source.append(source)
 				qualified_nodes_target.append(target)
 
-	return (qualified_nodes_source, qualified_nodes_target)
-	
+	return qualified_nodes_source, qualified_nodes_target
+
+
+@numba.njit(cache=True)
+def get_qualified_components_edge_comparison_clause(interpretations_edge, candidates_source, candidates_target, l, bnd, reverse_graph):
+	# Get all the qualified sources and targets for a particular clause
+	qualified_nodes_source = numba.typed.List.empty_list(node_type)
+	qualified_nodes_target = numba.typed.List.empty_list(node_type)
+	qualified_edges_numbers = numba.typed.List.empty_list(numba.types.float64)
+	for i, source in enumerate(candidates_source):
+		for target in candidates_target[i]:
+			edge = (source, target) if not reverse_graph else (target, source)
+			result, number = is_satisfied_edge_comparison(interpretations_edge, edge, (l, bnd))
+			if result:
+				qualified_nodes_source.append(source)
+				qualified_nodes_target.append(target)
+				qualified_edges_numbers.append(number)
+
+	return qualified_nodes_source, qualified_nodes_target, qualified_edges_numbers
+
+
+@numba.njit(cache=True)
+def compare_numbers_node_predicate(numbers_1, numbers_2, op, qualified_nodes_1, qualified_nodes_2):
+	result = False
+	final_qualified_nodes_1 = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_2 = numba.typed.List.empty_list(node_type)
+	for i in range(len(numbers_1)):
+		for j in range(len(numbers_2)):
+			if op == '<':
+				if numbers_1[i] < numbers_2[j]:
+					result = True
+			elif op == '<=':
+				if numbers_1[i] <= numbers_2[j]:
+					result = True
+			elif op == '>':
+				if numbers_1[i] > numbers_2[j]:
+					result = True
+			elif op == '>=':
+				if numbers_1[i] >= numbers_2[j]:
+					result = True
+			elif op == '==':
+				if numbers_1[i] == numbers_2[j]:
+					result = True
+			elif op == '!=':
+				if numbers_1[i] != numbers_2[j]:
+					result = True
+
+			if result:
+				final_qualified_nodes_1.append(qualified_nodes_1[i])
+				final_qualified_nodes_2.append(qualified_nodes_2[j])
+	return result, final_qualified_nodes_1, final_qualified_nodes_2
+
+
+@numba.njit(cache=True)
+def compare_numbers_edge_predicate(numbers_1, numbers_2, op, qualified_nodes_1a, qualified_nodes_1b, qualified_nodes_2a, qualified_nodes_2b):
+	result = False
+	final_qualified_nodes_1a = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_1b = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_2a = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_2b = numba.typed.List.empty_list(node_type)
+	for i in range(len(numbers_1)):
+		for j in range(len(numbers_2)):
+			if op == '<':
+				if numbers_1[i] < numbers_2[j]:
+					result = True
+			elif op == '<=':
+				if numbers_1[i] <= numbers_2[j]:
+					result = True
+			elif op == '>':
+				if numbers_1[i] > numbers_2[j]:
+					result = True
+			elif op == '>=':
+				if numbers_1[i] >= numbers_2[j]:
+					result = True
+			elif op == '==':
+				if numbers_1[i] == numbers_2[j]:
+					result = True
+			elif op == '!=':
+				if numbers_1[i] != numbers_2[j]:
+					result = True
+
+			if result:
+				final_qualified_nodes_1a.append(qualified_nodes_1a[i])
+				final_qualified_nodes_1b.append(qualified_nodes_1b[i])
+				final_qualified_nodes_2a.append(qualified_nodes_2a[j])
+				final_qualified_nodes_2b.append(qualified_nodes_2b[j])
+	return result, final_qualified_nodes_1a, final_qualified_nodes_1b, final_qualified_nodes_2a, final_qualified_nodes_2b
+
 
 @numba.njit(cache=True)
 def _satisfies_threshold(num_neigh, num_qualified_component, threshold):
@@ -955,6 +1499,10 @@ def _update_node(interpretations, comp, na, ipl, rule_trace, fp_cnt, t_cnt, stat
 		l, bnd = na
 		updated_bnds = numba.typed.List.empty_list(interval.interval_type)
 
+		# Add label to world if it is not there
+		if l not in world.world:
+			world.world[l] = interval.closed(0, 1)
+
 		# Check if update is necessary with previous bnd
 		prev_bnd = world.world[l].copy()
 
@@ -1008,7 +1556,7 @@ def _update_node(interpretations, comp, na, ipl, rule_trace, fp_cnt, t_cnt, stat
 					updated_bnds.append(world.world[p1])
 					if store_interpretation_changes:
 						rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p1, interval.closed(lower, upper)))
-		
+
 		# Gather convergence data
 		change = 0
 		if updated:
@@ -1040,6 +1588,10 @@ def _update_edge(interpretations, comp, na, ipl, rule_trace, fp_cnt, t_cnt, stat
 		l, bnd = na
 		updated_bnds = numba.typed.List.empty_list(interval.interval_type)
 
+		# Add label to world if it is not there
+		if l not in world.world:
+			world.world[l] = interval.closed(0, 1)
+
 		# Check if update is necessary with previous bnd
 		prev_bnd = world.world[l].copy()
 
@@ -1143,6 +1695,33 @@ def is_satisfied_node(interpretations, comp, na):
 	return result
 
 
+@numba.njit(cache=True)
+def is_satisfied_node_comparison(interpretations, comp, na):
+	result = False
+	number = 0
+	l, bnd = na
+	l_str = l.value
+
+	if not (l is None or bnd is None):
+		# This is to prevent a key error in case the label is a specific label
+		try:
+			world = interpretations[comp]
+			for world_l in world.world.keys():
+				world_l_str = world_l.value
+				if l_str in world_l_str and world_l_str[len(l_str)+1:].replace('.', '').replace('-', '').isdigit():
+					# The label is contained in the world
+					result = world.is_satisfied(world_l, na[1])
+					# Find the suffix number
+					number = str_to_float(world_l_str[len(l_str)+1:])
+					break
+
+		except:
+			result = False
+	else:
+		result = True
+	return result, number
+
+
 @numba.njit(cache=True)
 def are_satisfied_edge(interpretations, comp, nas):
 	result = True
@@ -1167,24 +1746,52 @@ def is_satisfied_edge(interpretations, comp, na):
 
 
 @numba.njit(cache=True)
-def influence(rule, annotations, weights):
+def is_satisfied_edge_comparison(interpretations, comp, na):
+	result = False
+	number = 0
+	l, bnd = na
+	l_str = l.value
+
+	if not (l is None or bnd is None):
+		# This is to prevent a key error in case the label is a specific label
+		try:
+			world = interpretations[comp]
+			for world_l in world.world.keys():
+				world_l_str = world_l.value
+				if l_str in world_l_str and world_l_str[len(l_str)+1:].replace('.', '').replace('-', '').isdigit():
+					# The label is contained in the world
+					result = world.is_satisfied(world_l, na[1])
+					# Find the suffix number
+					number = str_to_float(world_l_str[len(l_str)+1:])
+					break
+
+		except:
+			result = False
+	else:
+		result = True
+	return result, number
+
+
+@numba.njit(cache=True)
+def annotate(annotation_functions, rule, annotations, weights):
 	func_name = rule.get_annotation_function()
-	if func_name=='':
-		return interval.closed(rule.get_bnd().lower, rule.get_bnd().upper)
-	elif func_name=='average':
-		return ann_fn.average(annotations, weights)
-	elif func_name=='average_lower':
-		return ann_fn.average_lower(annotations, weights)
-	elif func_name=='minimum':
-		return ann_fn.minimum(annotations, weights)
-	elif func_name=='maximum':
-		return ann_fn.maximum(annotations, weights)
+	if func_name == '':
+		return rule.get_bnd().lower, rule.get_bnd().upper
+	else:
+		with numba.objmode(annotation='Tuple((float64, float64))'):
+			for func in annotation_functions:
+				if func.__name__ == func_name:
+					annotation = func(annotations, weights)
+		return annotation
 
 
 @numba.njit(cache=True)
 def check_consistent_node(interpretations, comp, na):
 	world = interpretations[comp]
-	bnd = world.world[na[0]]
+	if na[0] in world.world:
+		bnd = world.world[na[0]]
+	else:
+		bnd = interval.closed(0, 1)
 	if (na[1].lower > bnd.upper) or (bnd.lower > na[1].upper):
 		return False
 	else:
@@ -1194,7 +1801,10 @@ def check_consistent_node(interpretations, comp, na):
 @numba.njit(cache=True)
 def check_consistent_edge(interpretations, comp, na):
 	world = interpretations[comp]
-	bnd = world.world[na[0]]
+	if na[0] in world.world:
+		bnd = world.world[na[0]]
+	else:
+		bnd = interval.closed(0, 1)
 	if (na[1].lower > bnd.upper) or (bnd.lower > na[1].upper):
 		return False
 	else:
@@ -1232,49 +1842,50 @@ def resolve_inconsistency_node(interpretations, comp, na, ipl, t_cnt, fp_cnt, at
 
 @numba.njit(cache=True)
 def resolve_inconsistency_edge(interpretations, comp, na, ipl, t_cnt, fp_cnt, atom_trace, rule_trace, rule_trace_atoms, store_interpretation_changes):
-	world = interpretations[comp]
+	w = interpretations[comp]
 	if store_interpretation_changes:
 		rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, na[0], interval.closed(0,1)))
 		if atom_trace:
-			_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[na[0]], 'Inconsistency')
+			_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), w.world[na[0]], 'Inconsistency')
 	# Resolve inconsistency and set static
-	world.world[na[0]].set_lower_upper(0, 1)
-	world.world[na[0]].set_static(True)
+	w.world[na[0]].set_lower_upper(0, 1)
+	w.world[na[0]].set_static(True)
 	for p1, p2 in ipl:
 		if p1==na[0]:
 			if atom_trace:
-				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p2], 'Inconsistency')
-			world.world[p2].set_lower_upper(0, 1)
-			world.world[p2].set_static(True)
+				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), w.world[p2], 'Inconsistency')
+			w.world[p2].set_lower_upper(0, 1)
+			w.world[p2].set_static(True)
 			if store_interpretation_changes:
 				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p2, interval.closed(0,1)))
 
 		if p2==na[0]:
 			if atom_trace:
-				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p1], 'Inconsistency')
-			world.world[p1].set_lower_upper(0, 1)
-			world.world[p1].set_static(True)
+				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), w.world[p1], 'Inconsistency')
+			w.world[p1].set_lower_upper(0, 1)
+			w.world[p1].set_static(True)
 			if store_interpretation_changes:
 				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p1, interval.closed(0,1)))
 
 
 @numba.njit(cache=True)
-def _add_node(node, neighbors, nodes, interpretations_node):
+def _add_node(node, neighbors, reverse_neighbors, nodes, interpretations_node):
 	nodes.append(node)
 	neighbors[node] = numba.typed.List.empty_list(node_type)
+	reverse_neighbors[node] = numba.typed.List.empty_list(node_type)
 	interpretations_node[node] = world.World(numba.typed.List.empty_list(label.label_type))
 
 
 @numba.njit(cache=True)
-def _add_edge(source, target, neighbors, nodes, edges, l, interpretations_node, interpretations_edge):
+def _add_edge(source, target, neighbors, reverse_neighbors, nodes, edges, l, interpretations_node, interpretations_edge):
 	# If not a node, add to list of nodes and initialize neighbors
 	if source not in nodes:
-		_add_node(source, neighbors, nodes, interpretations_node)
+		_add_node(source, neighbors, reverse_neighbors, nodes, interpretations_node)
 
 	if target not in nodes:
-		_add_node(target, neighbors, nodes, interpretations_node)
+		_add_node(target, neighbors, reverse_neighbors, nodes, interpretations_node)
 
-	# Make sure edge doesnt already exist
+	# Make sure edge doesn't already exist
 	# Make sure, if l=='', not to add the label
 	# Make sure, if edge exists, that we don't override the l label if it exists
 	edge = (source, target)
@@ -1283,6 +1894,7 @@ def _add_edge(source, target, neighbors, nodes, edges, l, interpretations_node,
 		new_edge = True
 		edges.append(edge)
 		neighbors[source].append(target)
+		reverse_neighbors[target].append(source)
 		if l.value!='':
 			interpretations_edge[edge] = world.World(numba.typed.List([l]))
 		else:
@@ -1290,21 +1902,30 @@ def _add_edge(source, target, neighbors, nodes, edges, l, interpretations_node,
 	else:
 		if l not in interpretations_edge[edge].world and l.value!='':
 			new_edge = True
-			interpretations_edge[edge].world[l] = interval.closed(0,1)
+			interpretations_edge[edge].world[l] = interval.closed(0, 1)
 
-	return (edge, new_edge)
+	return edge, new_edge
 
 
 @numba.njit(cache=True)
-def _add_edges(sources, targets, neighbors, nodes, edges, l, interpretations_node, interpretations_edge):
+def _add_edges(sources, targets, neighbors, reverse_neighbors, nodes, edges, l, interpretations_node, interpretations_edge):
 	changes = 0
 	edges_added = numba.typed.List.empty_list(edge_type)
 	for source in sources:
 		for target in targets:
-			edge, new_edge = _add_edge(source, target, neighbors, nodes, edges, l, interpretations_node, interpretations_edge)
+			edge, new_edge = _add_edge(source, target, neighbors, reverse_neighbors, nodes, edges, l, interpretations_node, interpretations_edge)
 			edges_added.append(edge)
 			changes = changes+1 if new_edge else changes
-	return (edges_added, changes)
+	return edges_added, changes
+
+
+@numba.njit(cache=True)
+def _delete_edge(edge, neighbors, reverse_neighbors, edges, interpretations_edge):
+	source, target = edge
+	edges.remove(edge)
+	del interpretations_edge[edge]
+	neighbors[source].remove(target)
+	reverse_neighbors[target].remove(source)
 
 
 @numba.njit(cache=True)
@@ -1313,3 +1934,30 @@ def float_to_str(value):
 	decimal = int(value % 1 * 1000)
 	float_str = f'{number}.{decimal}'
 	return float_str
+
+
+@numba.njit(cache=True)
+def str_to_float(value):
+	decimal_pos = value.find('.')
+	if decimal_pos != -1:
+		after_decimal_len = len(value[decimal_pos+1:])
+	else:
+		after_decimal_len = 0
+	value = value.replace('.', '')
+	value = str_to_int(value)
+	value = value / 10**after_decimal_len
+	return value
+
+
+@numba.njit(cache=True)
+def str_to_int(value):
+	if value[0] == '-':
+		negative = True
+		value = value.replace('-','')
+	else:
+		negative = False
+	final_index, result = len(value) - 1, 0
+	for i, v in enumerate(value):
+		result += (ord(v) - 48) * (10 ** (final_index - i))
+	result = -result if negative else result
+	return result
diff --git a/pyreason/scripts/interpretation/interpretation_dict.py b/pyreason/scripts/interpretation/interpretation_dict.py
new file mode 100644
index 00000000..99e1a157
--- /dev/null
+++ b/pyreason/scripts/interpretation/interpretation_dict.py
@@ -0,0 +1,58 @@
+class InterpretationDict(dict):
+    """
+    This class is specific for the interpretation for a specific timestep.
+    """
+    def __int__(self):
+        super().__init__()
+
+    def __setitem__(self, key, value):
+        assert len(value) == 2, 'Lower bound and Upper bound are required to set an Interpretation'
+        self.__dict__[key] = (value[0], value[1])
+
+    def __getitem__(self, key):
+        if key not in self.__dict__.keys():
+            return tuple((0, 1))
+        else:
+            return self.__dict__[key]
+
+    def __repr__(self):
+        return repr(self.__dict__)
+
+    def __len__(self):
+        return len(self.__dict__)
+
+    def __delitem__(self, key):
+        del self.__dict__[key]
+
+    def clear(self):
+        return self.__dict__.clear()
+
+    def copy(self):
+        return self.__dict__.copy()
+
+    def has_key(self, k):
+        return k in self.__dict__
+
+    def update(self, *args, **kwargs):
+        return self.__dict__.update(*args, **kwargs)
+
+    def keys(self):
+        return self.__dict__.keys()
+
+    def values(self):
+        return self.__dict__.values()
+
+    def items(self):
+        return self.__dict__.items()
+
+    def pop(self, *args):
+        return self.__dict__.pop(*args)
+
+    def __cmp__(self, dict_):
+        return self.__cmp__(self.__dict__, dict_)
+
+    def __contains__(self, item):
+        return item in self.__dict__
+
+    def __iter__(self):
+        return iter(self.__dict__)
diff --git a/pyreason/scripts/interpretation/interpretation_parallel.py b/pyreason/scripts/interpretation/interpretation_parallel.py
new file mode 100644
index 00000000..4c57d352
--- /dev/null
+++ b/pyreason/scripts/interpretation/interpretation_parallel.py
@@ -0,0 +1,1956 @@
+import pyreason.scripts.numba_wrapper.numba_types.world_type as world
+import pyreason.scripts.numba_wrapper.numba_types.label_type as label
+import pyreason.scripts.numba_wrapper.numba_types.interval_type as interval
+from pyreason.scripts.interpretation.interpretation_dict import InterpretationDict
+
+import numba
+from numba import objmode, prange
+
+
+# Types for the dictionaries
+node_type = numba.types.string
+edge_type = numba.types.UniTuple(numba.types.string, 2)
+
+# Type for storing list of qualified nodes/edges
+list_of_nodes = numba.types.ListType(node_type)
+list_of_edges = numba.types.ListType(edge_type)
+
+# Type for facts to be applied
+facts_to_be_applied_node_type = numba.types.Tuple((numba.types.uint16, node_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean))
+facts_to_be_applied_edge_type = numba.types.Tuple((numba.types.uint16, edge_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean))
+
+# Type for returning list of applicable rules for a certain rule
+# node/edge, annotations, qualified nodes, qualified edges, edges to be added
+node_applicable_rule_type = numba.types.Tuple((
+	node_type,
+	numba.types.ListType(numba.types.ListType(interval.interval_type)),
+	numba.types.ListType(numba.types.ListType(node_type)),
+	numba.types.ListType(numba.types.ListType(edge_type)),
+	numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type))
+))
+
+edge_applicable_rule_type = numba.types.Tuple((
+	edge_type,
+	numba.types.ListType(numba.types.ListType(interval.interval_type)),
+	numba.types.ListType(numba.types.ListType(node_type)),
+	numba.types.ListType(numba.types.ListType(edge_type)),
+	numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type))
+))
+
+
+class Interpretation:
+	available_labels_node = []
+	available_labels_edge = []
+	specific_node_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(node_type))
+	specific_edge_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(edge_type))
+
+	def __init__(self, graph, ipl, annotation_functions, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes):
+		self.graph = graph
+		self.ipl = ipl
+		self.annotation_functions = annotation_functions
+		self.reverse_graph = reverse_graph
+		self.atom_trace = atom_trace
+		self.save_graph_attributes_to_rule_trace = save_graph_attributes_to_rule_trace
+		self.canonical = canonical
+		self.inconsistency_check = inconsistency_check
+		self.store_interpretation_changes = store_interpretation_changes
+
+		# For reasoning and reasoning again (contains previous time and previous fp operation cnt)
+		self.time = 0
+		self.prev_reasoning_data = numba.typed.List([0, 0])
+
+		# Initialize list of tuples for rules/facts to be applied, along with all the ground atoms that fired the rule. One to One correspondence between rules_to_be_applied_node and rules_to_be_applied_node_trace if atom_trace is true
+		self.rules_to_be_applied_node_trace = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(numba.types.ListType(node_type)), numba.types.ListType(numba.types.ListType(edge_type)), numba.types.string)))
+		self.rules_to_be_applied_edge_trace = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(numba.types.ListType(node_type)), numba.types.ListType(numba.types.ListType(edge_type)), numba.types.string)))
+		self.facts_to_be_applied_node_trace = numba.typed.List.empty_list(numba.types.string)
+		self.facts_to_be_applied_edge_trace = numba.typed.List.empty_list(numba.types.string)
+		self.rules_to_be_applied_node = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, node_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean)))
+		self.rules_to_be_applied_edge = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, edge_type, label.label_type, interval.interval_type, numba.types.boolean, numba.types.boolean)))
+		self.facts_to_be_applied_node = numba.typed.List.empty_list(facts_to_be_applied_node_type)
+		self.facts_to_be_applied_edge = numba.typed.List.empty_list(facts_to_be_applied_edge_type)
+		self.edges_to_be_added_node_rule = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type)))
+		self.edges_to_be_added_edge_rule = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(node_type), numba.types.ListType(node_type), label.label_type)))
+
+		# Keep track of all the rules that have affected each node/edge at each timestep/fp operation, and all ground atoms that have affected the rules as well. Keep track of previous bounds and name of the rule/fact here
+		self.rule_trace_node_atoms = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(numba.types.ListType(node_type)), numba.types.ListType(numba.types.ListType(edge_type)), interval.interval_type, numba.types.string)))
+		self.rule_trace_edge_atoms = numba.typed.List.empty_list(numba.types.Tuple((numba.types.ListType(numba.types.ListType(node_type)), numba.types.ListType(numba.types.ListType(edge_type)), interval.interval_type, numba.types.string)))
+		self.rule_trace_node = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, numba.types.uint16, node_type, label.label_type, interval.interval_type)))
+		self.rule_trace_edge = numba.typed.List.empty_list(numba.types.Tuple((numba.types.uint16, numba.types.uint16, edge_type, label.label_type, interval.interval_type)))
+
+		# Nodes and edges of the graph
+		self.nodes = numba.typed.List(self.graph.nodes())
+		self.edges = numba.typed.List(self.graph.edges())
+
+		# Make sure they are correct type
+		if len(self.available_labels_node)==0:
+			self.available_labels_node = numba.typed.List.empty_list(label.label_type)
+		else:
+			self.available_labels_node = numba.typed.List(self.available_labels_node)
+		if len(self.available_labels_edge)==0:
+			self.available_labels_edge = numba.typed.List.empty_list(label.label_type)
+		else:
+			self.available_labels_edge = numba.typed.List(self.available_labels_edge)
+
+		self.interpretations_node = self._init_interpretations_node(numba.typed.List(self.graph.nodes()), self.available_labels_node, self.specific_node_labels)
+		self.interpretations_edge = self._init_interpretations_edge(numba.typed.List(self.graph.edges()), self.available_labels_edge, self.specific_edge_labels)
+
+		# Setup graph neighbors and reverse neighbors
+		self.neighbors = numba.typed.Dict.empty(key_type=node_type, value_type=numba.types.ListType(node_type))
+		for n in self.graph.nodes():
+			l = numba.typed.List.empty_list(node_type)
+			[l.append(neigh) for neigh in self.graph.neighbors(n)]
+			self.neighbors[n] = l
+
+		self.reverse_neighbors = self._init_reverse_neighbors(self.neighbors)
+
+	@staticmethod
+	@numba.njit(cache=False)
+	def _init_reverse_neighbors(neighbors):
+		reverse_neighbors = numba.typed.Dict.empty(key_type=node_type, value_type=list_of_nodes)
+		for n, neighbor_nodes in neighbors.items():
+			for neighbor_node in neighbor_nodes:
+				if neighbor_node in reverse_neighbors and n not in reverse_neighbors[neighbor_node]:
+					reverse_neighbors[neighbor_node].append(n)
+				else:
+					reverse_neighbors[neighbor_node] = numba.typed.List([n])
+			# This makes sure each node has a value
+			if n not in reverse_neighbors:
+				reverse_neighbors[n] = numba.typed.List.empty_list(node_type)
+
+		return reverse_neighbors
+
+	@staticmethod
+	@numba.njit(cache=False)
+	def _init_interpretations_node(nodes, available_labels, specific_labels):
+		interpretations = numba.typed.Dict.empty(key_type=node_type, value_type=world.world_type)
+		# General labels
+		for n in nodes:
+			interpretations[n] = world.World(available_labels)
+		# Specific labels
+		for l, ns in specific_labels.items():
+			for n in ns:
+				interpretations[n].world[l] = interval.closed(0.0, 1.0)
+
+		return interpretations
+
+	@staticmethod
+	@numba.njit(cache=False)
+	def _init_interpretations_edge(edges, available_labels, specific_labels):
+		interpretations = numba.typed.Dict.empty(key_type=edge_type, value_type=world.world_type)
+		# General labels
+		for e in edges:
+			interpretations[e] = world.World(available_labels)
+		# Specific labels
+		for l, es in specific_labels.items():
+			for e in es:
+				interpretations[e].world[l] = interval.closed(0.0, 1.0)
+
+		return interpretations
+
+	@staticmethod
+	@numba.njit(cache=False)
+	def _init_convergence(convergence_bound_threshold, convergence_threshold):
+		if convergence_bound_threshold==-1 and convergence_threshold==-1:
+			convergence_mode = 'perfect_convergence'
+			convergence_delta = 0
+		elif convergence_bound_threshold==-1:
+			convergence_mode = 'delta_interpretation'
+			convergence_delta = convergence_threshold
+		else:
+			convergence_mode = 'delta_bound'
+			convergence_delta = convergence_bound_threshold
+		return convergence_mode, convergence_delta
+
+	def start_fp(self, tmax, facts_node, facts_edge, rules, verbose, convergence_threshold, convergence_bound_threshold, again=False):
+		self.tmax = tmax
+		self._convergence_mode, self._convergence_delta = self._init_convergence(convergence_bound_threshold, convergence_threshold)
+		max_facts_time = self._init_facts(facts_node, facts_edge, self.facts_to_be_applied_node, self.facts_to_be_applied_edge, self.facts_to_be_applied_node_trace, self.facts_to_be_applied_edge_trace, self.atom_trace)
+		self._start_fp(rules, max_facts_time, verbose, again)
+
+	@staticmethod
+	@numba.njit(cache=False)
+	def _init_facts(facts_node, facts_edge, facts_to_be_applied_node, facts_to_be_applied_edge, facts_to_be_applied_node_trace, facts_to_be_applied_edge_trace, atom_trace):
+		max_time = 0
+		for fact in facts_node:
+			for t in range(fact.get_time_lower(), fact.get_time_upper() + 1):
+				max_time = max(max_time, t)
+				name = fact.get_name()
+				graph_attribute = True if name=='graph-attribute-fact' else False
+				facts_to_be_applied_node.append((numba.types.uint16(t), fact.get_component(), fact.get_label(), fact.get_bound(), fact.static, graph_attribute))
+				if atom_trace:
+					facts_to_be_applied_node_trace.append(fact.get_name())
+		for fact in facts_edge:
+			for t in range(fact.get_time_lower(), fact.get_time_upper() + 1):
+				max_time = max(max_time, t)
+				name = fact.get_name()
+				graph_attribute = True if name=='graph-attribute-fact' else False
+				facts_to_be_applied_edge.append((numba.types.uint16(t), fact.get_component(), fact.get_label(), fact.get_bound(), fact.static, graph_attribute))
+				if atom_trace:
+					facts_to_be_applied_edge_trace.append(fact.get_name())
+		return max_time
+
+	def _start_fp(self, rules, max_facts_time, verbose, again):
+		fp_cnt, t = self.reason(self.interpretations_node, self.interpretations_edge, self.tmax, self.prev_reasoning_data, rules, self.nodes, self.edges, self.neighbors, self.reverse_neighbors, self.rules_to_be_applied_node, self.rules_to_be_applied_edge, self.edges_to_be_added_node_rule, self.edges_to_be_added_edge_rule, self.rules_to_be_applied_node_trace, self.rules_to_be_applied_edge_trace, self.facts_to_be_applied_node, self.facts_to_be_applied_edge, self.facts_to_be_applied_node_trace, self.facts_to_be_applied_edge_trace, self.ipl, self.rule_trace_node, self.rule_trace_edge, self.rule_trace_node_atoms, self.rule_trace_edge_atoms, self.reverse_graph, self.atom_trace, self.save_graph_attributes_to_rule_trace, self.canonical, self.inconsistency_check, self.store_interpretation_changes, max_facts_time, self.annotation_functions, self._convergence_mode, self._convergence_delta, verbose, again)
+		self.time = t - 1
+		# If we need to reason again, store the next timestep to start from
+		self.prev_reasoning_data[0] = t
+		self.prev_reasoning_data[1] = fp_cnt
+		if verbose:
+			print('Fixed Point iterations:', fp_cnt)
+
+	@staticmethod
+	@numba.njit(cache=False)
+	def reason(interpretations_node, interpretations_edge, tmax, prev_reasoning_data, rules, nodes, edges, neighbors, reverse_neighbors, rules_to_be_applied_node, rules_to_be_applied_edge, edges_to_be_added_node_rule, edges_to_be_added_edge_rule, rules_to_be_applied_node_trace, rules_to_be_applied_edge_trace, facts_to_be_applied_node, facts_to_be_applied_edge, facts_to_be_applied_node_trace, facts_to_be_applied_edge_trace, ipl, rule_trace_node, rule_trace_edge, rule_trace_node_atoms, rule_trace_edge_atoms, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes, max_facts_time, annotation_functions, convergence_mode, convergence_delta, verbose, again):
+		t = prev_reasoning_data[0]
+		fp_cnt = prev_reasoning_data[1]
+		max_rules_time = 0
+		timestep_loop = True
+		facts_to_be_applied_node_new = numba.typed.List.empty_list(facts_to_be_applied_node_type)
+		facts_to_be_applied_edge_new = numba.typed.List.empty_list(facts_to_be_applied_edge_type)
+		rules_to_remove_idx = numba.typed.List.empty_list(numba.types.int64)
+		while timestep_loop:
+			if t==tmax:
+				timestep_loop = False
+			if verbose:
+				with objmode():
+					print('Timestep:', t, flush=True)
+			# Reset Interpretation at beginning of timestep if non-canonical
+			if t>0 and not canonical:
+				# Reset nodes (only if not static)
+				for n in nodes:
+					w = interpretations_node[n].world
+					for l in w:
+						if not w[l].is_static():
+							w[l].reset()
+
+				# Reset edges (only if not static)
+				for e in edges:
+					w = interpretations_edge[e].world
+					for l in w:
+						if not w[l].is_static():
+							w[l].reset()
+
+			# Convergence parameters
+			changes_cnt = 0
+			bound_delta = 0
+			update = False
+
+			# Parameters for immediate rules
+			immediate_node_rule_fire = False
+			immediate_edge_rule_fire = False
+			immediate_rule_applied = False
+			# When delta_t = 0, we don't want to check the same rule with the same node/edge after coming back to the fp operator
+			nodes_to_skip = numba.typed.Dict.empty(key_type=numba.types.int64, value_type=list_of_nodes)
+			edges_to_skip = numba.typed.Dict.empty(key_type=numba.types.int64, value_type=list_of_edges)
+			# Initialize the above
+			for i in range(len(rules)):
+				nodes_to_skip[i] = numba.typed.List.empty_list(node_type)
+				edges_to_skip[i] = numba.typed.List.empty_list(edge_type)
+
+			# Start by applying facts
+			# Nodes
+			facts_to_be_applied_node_new.clear()
+			for i in range(len(facts_to_be_applied_node)):
+				if facts_to_be_applied_node[i][0]==t:
+					comp, l, bnd, static, graph_attribute = facts_to_be_applied_node[i][1], facts_to_be_applied_node[i][2], facts_to_be_applied_node[i][3], facts_to_be_applied_node[i][4], facts_to_be_applied_node[i][5]
+					# Check if bnd is static. Then no need to update, just add to rule trace, check if graph attribute and add ipl complement to rule trace as well
+					if l in interpretations_node[comp].world and interpretations_node[comp].world[l].is_static():
+						# Check if we should even store any of the changes to the rule trace etc.
+						# Inverse of this is: if not save_graph_attributes_to_rule_trace and graph_attribute
+						if (save_graph_attributes_to_rule_trace or not graph_attribute) and store_interpretation_changes:
+							rule_trace_node.append((numba.types.uint16(t), numba.types.uint16(fp_cnt), comp, l, bnd))
+							if atom_trace:
+								_update_rule_trace(rule_trace_node_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), bnd, facts_to_be_applied_node_trace[i])
+							for p1, p2 in ipl:
+								if p1==l:
+									rule_trace_node.append((numba.types.uint16(t), numba.types.uint16(fp_cnt), comp, p2, interpretations_node[comp].world[p2]))
+									if atom_trace:
+										_update_rule_trace(rule_trace_node_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), interpretations_node[comp].world[p2], facts_to_be_applied_node_trace[i])
+								elif p2==l:
+									rule_trace_node.append((numba.types.uint16(t), numba.types.uint16(fp_cnt), comp, p1, interpretations_node[comp].world[p1]))
+									if atom_trace:
+										_update_rule_trace(rule_trace_node_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), interpretations_node[comp].world[p1], facts_to_be_applied_node_trace[i])
+
+					else:
+						# Check for inconsistencies (multiple facts)
+						if check_consistent_node(interpretations_node, comp, (l, bnd)):
+							mode = 'graph-attribute-fact' if graph_attribute else 'fact'
+							u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, i, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode=mode)
+
+							update = u or update
+							# Update convergence params
+							if convergence_mode=='delta_bound':
+								bound_delta = max(bound_delta, changes)
+							else:
+								changes_cnt += changes
+						# Resolve inconsistency if necessary otherwise override bounds
+						else:
+							if inconsistency_check:
+								resolve_inconsistency_node(interpretations_node, comp, (l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_node, rule_trace_node_atoms, store_interpretation_changes)
+							else:
+								mode = 'graph-attribute-fact' if graph_attribute else 'fact'
+								u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, i, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode=mode, override=True)
+
+								update = u or update
+								# Update convergence params
+								if convergence_mode=='delta_bound':
+									bound_delta = max(bound_delta, changes)
+								else:
+									changes_cnt += changes
+
+					if static:
+						facts_to_be_applied_node_new.append((numba.types.uint16(facts_to_be_applied_node[i][0]+1), comp, l, bnd, static, graph_attribute))
+
+				# If time doesn't match, fact to be applied later
+				else:
+					facts_to_be_applied_node_new.append(facts_to_be_applied_node[i])
+
+			# Update list of facts with ones that have not been applied yet (delete applied facts)
+			facts_to_be_applied_node[:] = facts_to_be_applied_node_new.copy()
+			facts_to_be_applied_node_new.clear()
+
+			# Edges
+			facts_to_be_applied_edge_new.clear()
+			for i in range(len(facts_to_be_applied_edge)):
+				if facts_to_be_applied_edge[i][0]==t:
+					comp, l, bnd, static, graph_attribute = facts_to_be_applied_edge[i][1], facts_to_be_applied_edge[i][2], facts_to_be_applied_edge[i][3], facts_to_be_applied_edge[i][4], facts_to_be_applied_edge[i][5]
+					# Check if bnd is static. Then no need to update, just add to rule trace, check if graph attribute, and add ipl complement to rule trace as well
+					if l in interpretations_edge[comp].world and interpretations_edge[comp].world[l].is_static():
+						# Inverse of this is: if not save_graph_attributes_to_rule_trace and graph_attribute
+						if (save_graph_attributes_to_rule_trace or not graph_attribute) and store_interpretation_changes:
+							rule_trace_edge.append((numba.types.uint16(t), numba.types.uint16(fp_cnt), comp, l, interpretations_edge[comp].world[l]))
+							if atom_trace:
+								_update_rule_trace(rule_trace_edge_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), bnd, facts_to_be_applied_edge_trace[i])
+							for p1, p2 in ipl:
+								if p1==l:
+									rule_trace_edge.append((numba.types.uint16(t), numba.types.uint16(fp_cnt), comp, p2, interpretations_edge[comp].world[p2]))
+									if atom_trace:
+										_update_rule_trace(rule_trace_edge_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), interpretations_edge[comp].world[p2], facts_to_be_applied_edge_trace[i])
+								elif p2==l:
+									rule_trace_edge.append((numba.types.uint16(t), numba.types.uint16(fp_cnt), comp, p1, interpretations_edge[comp].world[p1]))
+									if atom_trace:
+										_update_rule_trace(rule_trace_edge_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), interpretations_edge[comp].world[p1], facts_to_be_applied_edge_trace[i])
+					else:
+						# Check for inconsistencies
+						if check_consistent_edge(interpretations_edge, comp, (l, bnd)):
+							mode = 'graph-attribute-fact' if graph_attribute else 'fact'
+							u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, i, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode=mode)
+
+							update = u or update
+							# Update convergence params
+							if convergence_mode=='delta_bound':
+								bound_delta = max(bound_delta, changes)
+							else:
+								changes_cnt += changes
+						# Resolve inconsistency
+						else:
+							if inconsistency_check:
+								resolve_inconsistency_edge(interpretations_edge, comp, (l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_edge, rule_trace_edge_atoms, store_interpretation_changes)
+							else:
+								mode = 'graph-attribute-fact' if graph_attribute else 'fact'
+								u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, i, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode=mode, override=True)
+
+								update = u or update
+								# Update convergence params
+								if convergence_mode=='delta_bound':
+									bound_delta = max(bound_delta, changes)
+								else:
+									changes_cnt += changes
+
+					if static:
+						facts_to_be_applied_edge_new.append((numba.types.uint16(facts_to_be_applied_edge[i][0]+1), comp, l, bnd, static, graph_attribute))
+
+				# Time doesn't match, fact to be applied later
+				else:
+					facts_to_be_applied_edge_new.append(facts_to_be_applied_edge[i])
+
+			# Update list of facts with ones that have not been applied yet (delete applied facts)
+			facts_to_be_applied_edge[:] = facts_to_be_applied_edge_new.copy()
+			facts_to_be_applied_edge_new.clear()
+
+			in_loop = True
+			while in_loop:
+				# This will become true only if delta_t = 0 for some rule, otherwise we go to the next timestep
+				in_loop = False
+
+				# Apply the rules that need to be applied at this timestep
+				# Nodes
+				rules_to_remove_idx.clear()
+				for idx, i in enumerate(rules_to_be_applied_node):
+					# If we are coming here from an immediate rule firing with delta_t=0 we have to apply that one rule. Which was just added to the list to_be_applied
+					if immediate_node_rule_fire and rules_to_be_applied_node[-1][4]:
+						i = rules_to_be_applied_node[-1]
+						idx = len(rules_to_be_applied_node) - 1
+
+					if i[0]==t:
+						comp, l, bnd, immediate, set_static = i[1], i[2], i[3], i[4], i[5]
+						sources, targets, edge_l = edges_to_be_added_node_rule[idx]
+						edges_added, changes = _add_edges(sources, targets, neighbors, reverse_neighbors, nodes, edges, edge_l, interpretations_node, interpretations_edge)
+						changes_cnt += changes
+
+						# Update bound for newly added edges. Use bnd to update all edges if label is specified, else use bnd to update normally
+						if edge_l.value!='':
+							for e in edges_added:
+								if check_consistent_edge(interpretations_edge, e, (edge_l, bnd)):
+									u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
+
+									update = u or update
+
+									# Update convergence params
+									if convergence_mode=='delta_bound':
+										bound_delta = max(bound_delta, changes)
+									else:
+										changes_cnt += changes
+								# Resolve inconsistency
+								else:
+									if inconsistency_check:
+										resolve_inconsistency_edge(interpretations_edge, e, (edge_l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_edge, rule_trace_edge_atoms, store_interpretation_changes)
+									else:
+										u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
+
+										update = u or update
+
+										# Update convergence params
+										if convergence_mode=='delta_bound':
+											bound_delta = max(bound_delta, changes)
+										else:
+											changes_cnt += changes
+						else:
+							# Check for inconsistencies
+							if check_consistent_node(interpretations_node, comp, (l, bnd)):
+								u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode='rule')
+
+								update = u or update
+								# Update convergence params
+								if convergence_mode=='delta_bound':
+									bound_delta = max(bound_delta, changes)
+								else:
+									changes_cnt += changes
+							# Resolve inconsistency
+							else:
+								if inconsistency_check:
+									resolve_inconsistency_node(interpretations_node, comp, (l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_node, rule_trace_node_atoms, store_interpretation_changes)
+								else:
+									u, changes = _update_node(interpretations_node, comp, (l, bnd), ipl, rule_trace_node, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_node_trace, idx, facts_to_be_applied_node_trace, rule_trace_node_atoms, store_interpretation_changes, mode='rule', override=True)
+
+									update = u or update
+									# Update convergence params
+									if convergence_mode=='delta_bound':
+										bound_delta = max(bound_delta, changes)
+									else:
+										changes_cnt += changes
+
+						# Delete rules that have been applied from list by adding index to list
+						rules_to_remove_idx.append(idx)
+
+						# Break out of the apply rules loop if a rule is immediate. Then we go to the fp operator and check for other applicable rules then come back
+						if immediate:
+							# If delta_t=0 we want to apply one rule and go back to the fp operator
+							# If delta_t>0 we want to come back here and apply the rest of the rules
+							if immediate_edge_rule_fire:
+								break
+							elif not immediate_edge_rule_fire and u:
+								immediate_rule_applied = True
+								break
+
+				# Remove from rules to be applied and edges to be applied lists after coming out from loop
+				rules_to_be_applied_node[:] = numba.typed.List([rules_to_be_applied_node[i] for i in range(len(rules_to_be_applied_node)) if i not in rules_to_remove_idx])
+				edges_to_be_added_node_rule[:] = numba.typed.List([edges_to_be_added_node_rule[i] for i in range(len(edges_to_be_added_node_rule)) if i not in rules_to_remove_idx])
+				if atom_trace:
+					rules_to_be_applied_node_trace[:] = numba.typed.List([rules_to_be_applied_node_trace[i] for i in range(len(rules_to_be_applied_node_trace)) if i not in rules_to_remove_idx])
+
+				# Edges
+				rules_to_remove_idx.clear()
+				for idx, i in enumerate(rules_to_be_applied_edge):
+					# If we broke from above loop to apply more rules, then break from here
+					if immediate_rule_applied and not immediate_edge_rule_fire:
+						break
+					# If we are coming here from an immediate rule firing with delta_t=0 we have to apply that one rule. Which was just added to the list to_be_applied
+					if immediate_edge_rule_fire and rules_to_be_applied_edge[-1][4]:
+						i = rules_to_be_applied_edge[-1]
+						idx = len(rules_to_be_applied_edge) - 1
+
+					if i[0]==t:
+						comp, l, bnd, immediate, set_static = i[1], i[2], i[3], i[4], i[5]
+						sources, targets, edge_l = edges_to_be_added_edge_rule[idx]
+						edges_added, changes = _add_edges(sources, targets, neighbors, reverse_neighbors, nodes, edges, edge_l, interpretations_node, interpretations_edge)
+						changes_cnt += changes
+
+						# Update bound for newly added edges. Use bnd to update all edges if label is specified, else use bnd to update normally
+						if edge_l.value!='':
+							for e in edges_added:
+								if check_consistent_edge(interpretations_edge, e, (edge_l, bnd)):
+									u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
+
+									update = u or update
+
+									# Update convergence params
+									if convergence_mode=='delta_bound':
+										bound_delta = max(bound_delta, changes)
+									else:
+										changes_cnt += changes
+								# Resolve inconsistency
+								else:
+									if inconsistency_check:
+										resolve_inconsistency_edge(interpretations_edge, e, (edge_l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_edge, rule_trace_edge_atoms, store_interpretation_changes)
+									else:
+										u, changes = _update_edge(interpretations_edge, e, (edge_l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
+
+										update = u or update
+
+										# Update convergence params
+										if convergence_mode=='delta_bound':
+											bound_delta = max(bound_delta, changes)
+										else:
+											changes_cnt += changes
+
+						else:
+							# Check for inconsistencies
+							if check_consistent_edge(interpretations_edge, comp, (l, bnd)):
+								u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule')
+
+								update = u or update
+								# Update convergence params
+								if convergence_mode=='delta_bound':
+									bound_delta = max(bound_delta, changes)
+								else:
+									changes_cnt += changes
+							# Resolve inconsistency
+							else:
+								if inconsistency_check:
+									resolve_inconsistency_edge(interpretations_edge, comp, (l, bnd), ipl, t, fp_cnt, atom_trace, rule_trace_edge, rule_trace_edge_atoms, store_interpretation_changes)
+								else:
+									u, changes = _update_edge(interpretations_edge, comp, (l, bnd), ipl, rule_trace_edge, fp_cnt, t, set_static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_edge_trace, idx, facts_to_be_applied_edge_trace, rule_trace_edge_atoms, store_interpretation_changes, mode='rule', override=True)
+
+									update = u or update
+									# Update convergence params
+									if convergence_mode=='delta_bound':
+										bound_delta = max(bound_delta, changes)
+									else:
+										changes_cnt += changes
+
+						# Delete rules that have been applied from list by adding the index to list
+						rules_to_remove_idx.append(idx)
+
+						# Break out of the apply rules loop if a rule is immediate. Then we go to the fp operator and check for other applicable rules then come back
+						if immediate:
+							# If t=0 we want to apply one rule and go back to the fp operator
+							# If t>0 we want to come back here and apply the rest of the rules
+							if immediate_edge_rule_fire:
+								break
+							elif not immediate_edge_rule_fire and u:
+								immediate_rule_applied = True
+								break
+
+				# Remove from rules to be applied and edges to be applied lists after coming out from loop
+				rules_to_be_applied_edge[:] = numba.typed.List([rules_to_be_applied_edge[i] for i in range(len(rules_to_be_applied_edge)) if i not in rules_to_remove_idx])
+				edges_to_be_added_edge_rule[:] = numba.typed.List([edges_to_be_added_edge_rule[i] for i in range(len(edges_to_be_added_edge_rule)) if i not in rules_to_remove_idx])
+				if atom_trace:
+					rules_to_be_applied_edge_trace[:] = numba.typed.List([rules_to_be_applied_edge_trace[i] for i in range(len(rules_to_be_applied_edge_trace)) if i not in rules_to_remove_idx])
+
+				# Fixed point
+				# if update or immediate_node_rule_fire or immediate_edge_rule_fire or immediate_rule_applied:
+				if update:
+					# Increase fp operator count only if not an immediate rule
+					if not (immediate_node_rule_fire or immediate_edge_rule_fire):
+						fp_cnt += 1
+
+					for i in range(len(rules)):
+						rule = rules[i]
+						immediate_rule = rule.is_immediate_rule()
+						immediate_node_rule_fire = False
+						immediate_edge_rule_fire = False
+
+						# Only go through if the rule can be applied within the given timesteps, or we're running until convergence
+						delta_t = rule.get_delta()
+						if t + delta_t <= tmax or tmax == -1 or again:
+							applicable_node_rules = _ground_node_rule(rule, interpretations_node, interpretations_edge, nodes, neighbors, reverse_neighbors, atom_trace, reverse_graph, nodes_to_skip[i])
+							applicable_edge_rules = _ground_edge_rule(rule, interpretations_node, interpretations_edge, nodes, edges, neighbors, reverse_neighbors, atom_trace, reverse_graph, edges_to_skip[i])
+
+							# Loop through applicable rules and add them to the rules to be applied for later or next fp operation
+							for applicable_rule in applicable_node_rules:
+								n, annotations, qualified_nodes, qualified_edges, edges_to_add = applicable_rule
+								# If there is an edge to add or the predicate doesn't exist or the interpretation is not static
+								if len(edges_to_add[0]) > 0 or rule.get_target() not in interpretations_node[n].world or not interpretations_node[n].world[rule.get_target()].is_static():
+									bnd = annotate(annotation_functions, rule, annotations, rule.get_weights())
+									bnd = interval.closed(bnd[0], bnd[1])
+									max_rules_time = max(max_rules_time, t + delta_t)
+									edges_to_be_added_node_rule.append(edges_to_add)
+									rules_to_be_applied_node.append((numba.types.uint16(t + delta_t), n, rule.get_target(), bnd, immediate_rule, rule.is_static_rule()))
+									if atom_trace:
+										rules_to_be_applied_node_trace.append((qualified_nodes, qualified_edges, rule.get_name()))
+
+									# We apply a rule on a node/edge only once in each timestep to prevent it from being added to the to_be_added list continuously (this will improve performance
+									nodes_to_skip[i].append(n)
+
+									# Handle loop parameters for the next (maybe) fp operation
+									# If it is a t=0 rule or an immediate rule we want to go back for another fp operation to check for new rules that may fire
+									# Next fp operation we will skip this rule on this node because anyway there won't be an update
+									if delta_t == 0:
+										in_loop = True
+										update = False
+									if immediate_rule and delta_t == 0:
+										# immediate_rule_fire becomes True because we still need to check for more eligible rules, we're not done.
+										in_loop = True
+										update = True
+										immediate_node_rule_fire = True
+										break
+
+							# Break, apply immediate rule then come back to check for more applicable rules
+							if immediate_node_rule_fire:
+								break
+
+							for applicable_rule in applicable_edge_rules:
+								e, annotations, qualified_nodes, qualified_edges, edges_to_add = applicable_rule
+								# If there is an edge to add or the predicate doesn't exist or the interpretation is not static
+								if len(edges_to_add[0]) > 0 or rule.get_target() not in interpretations_edge[e].world or not interpretations_edge[e].world[rule.get_target()].is_static():
+									bnd = annotate(annotation_functions, rule, annotations, rule.get_weights())
+									bnd = interval.closed(bnd[0], bnd[1])
+									max_rules_time = max(max_rules_time, t+delta_t)
+									edges_to_be_added_edge_rule.append(edges_to_add)
+									rules_to_be_applied_edge.append((numba.types.uint16(t+delta_t), e, rule.get_target(), bnd, immediate_rule, rule.is_static_rule()))
+									if atom_trace:
+										rules_to_be_applied_edge_trace.append((qualified_nodes, qualified_edges, rule.get_name()))
+
+									# We apply a rule on a node/edge only once in each timestep to prevent it from being added to the to_be_added list continuously (this will improve performance
+									edges_to_skip[i].append(e)
+
+									# Handle loop parameters for the next (maybe) fp operation
+									# If it is a t=0 rule or an immediate rule we want to go back for another fp operation to check for new rules that may fire
+									# Next fp operation we will skip this rule on this node because anyway there won't be an update
+									if delta_t == 0:
+										in_loop = True
+										update = False
+									if immediate_rule and delta_t == 0:
+										# immediate_rule_fire becomes True because we still need to check for more eligible rules, we're not done.
+										in_loop = True
+										update = True
+										immediate_edge_rule_fire = True
+										break
+
+							# Break, apply immediate rule then come back to check for more applicable rules
+							if immediate_edge_rule_fire:
+								break
+
+					# Go through all the rules and go back to applying the rules if we came here because of an immediate rule where delta_t>0
+					if immediate_rule_applied and not (immediate_node_rule_fire or immediate_edge_rule_fire):
+						immediate_rule_applied = False
+						in_loop = True
+						update = False
+						continue
+
+			# Check for convergence after each timestep (perfect convergence or convergence specified by user)
+			# Check number of changed interpretations or max bound change
+			# User specified convergence
+			if convergence_mode=='delta_interpretation':
+				if changes_cnt <= convergence_delta:
+					if verbose:
+						print(f'\nConverged at time: {t} with {int(changes_cnt)} changes from the previous interpretation')
+					# Be consistent with time returned when we don't converge
+					t += 1
+					break
+			elif convergence_mode=='delta_bound':
+				if bound_delta <= convergence_delta:
+					if verbose:
+						print(f'\nConverged at time: {t} with {float_to_str(bound_delta)} as the maximum bound change from the previous interpretation')
+					# Be consistent with time returned when we don't converge
+					t += 1
+					break
+			# Perfect convergence
+			# Make sure there are no rules to be applied, and no facts that will be applied in the future. We do this by checking the max time any rule/fact is applicable
+			# If no more rules/facts to be applied
+			elif convergence_mode=='perfect_convergence':
+				if t>=max_facts_time and t>=max_rules_time:
+					if verbose:
+						print(f'\nConverged at time: {t}')
+					# Be consistent with time returned when we don't converge
+					t += 1
+					break
+
+			# Increment t
+			t += 1
+
+		return fp_cnt, t
+
+	def add_edge(self, edge, l):
+		# This function is useful for pyreason gym, called externally
+		_add_edge(edge[0], edge[1], self.neighbors, self.reverse_neighbors, self.nodes, self.edges, l, self.interpretations_node, self.interpretations_edge)
+
+	def delete_edge(self, edge):
+		# This function is useful for pyreason gym, called externally
+		_delete_edge(edge, self.neighbors, self.reverse_neighbors, self.edges, self.interpretations_edge)
+
+	def get_interpretation_dict(self):
+		# This function can be called externally to retrieve a dict of the interpretation values
+		# Only values in the rule trace will be added
+
+		# Initialize interpretations for each time and node and edge
+		interpretations = {}
+		for t in range(self.tmax+1):
+			interpretations[t] = {}
+			for node in self.nodes:
+				interpretations[t][node] = InterpretationDict()
+			for edge in self.edges:
+				interpretations[t][edge] = InterpretationDict()
+
+		# Update interpretation nodes
+		for change in self.rule_trace_node:
+			time, _, node, l, bnd = change
+			interpretations[time][node][l._value] = (bnd.lower, bnd.upper)
+
+			# If canonical, update all following timesteps as well
+			if self. canonical:
+				for t in range(time+1, self.tmax+1):
+					interpretations[t][node][l._value] = (bnd.lower, bnd.upper)
+
+		# Update interpretation edges
+		for change in self.rule_trace_edge:
+			time, _, edge, l, bnd, = change
+			interpretations[time][edge][l._value] = (bnd.lower, bnd.upper)
+
+			# If canonical, update all following timesteps as well
+			if self. canonical:
+				for t in range(time+1, self.tmax+1):
+					interpretations[t][edge][l._value] = (bnd.lower, bnd.upper)
+
+		return interpretations
+
+
+@numba.njit(cache=False, parallel=True)
+def _ground_node_rule(rule, interpretations_node, interpretations_edge, nodes, neighbors, reverse_neighbors, atom_trace, reverse_graph, nodes_to_skip):
+	# Extract rule params
+	rule_type = rule.get_type()
+	clauses = rule.get_clauses()
+	thresholds = rule.get_thresholds()
+	ann_fn = rule.get_annotation_function()
+	rule_edges = rule.get_edges()
+
+	# We return a list of tuples which specify the target nodes/edges that have made the rule body true
+	applicable_rules = numba.typed.List.empty_list(node_applicable_rule_type)
+
+	# Return empty list if rule is not node rule and if we are not inferring edges
+	if rule_type != 'node' and rule_edges[0] == '':
+		return applicable_rules
+
+	# Steps
+	# 1. Loop through all nodes and evaluate each clause with that node and check the truth with the thresholds
+	# 2. Inside the clause loop it may be necessary to loop through all nodes/edges while grounding the variables
+	# 3. If the clause is true add the qualified nodes and qualified edges to the atom trace, if on. Break otherwise
+	# 4. After going through all clauses, add to the annotations list all the annotations of the specified subset. These will be passed to the annotation function
+	# 5. Finally, if there are any edges to be added, place them in the list
+
+	for piter in prange(len(nodes)):
+		target_node = nodes[piter]
+		if target_node in nodes_to_skip:
+			continue
+		# Initialize dictionary where keys are strings (x1, x2 etc.) and values are lists of qualified neighbors
+		# Keep track of qualified nodes and qualified edges
+		# If it's a node clause update (x1 or x2 etc.) qualified neighbors, if it's an edge clause update the qualified neighbors for the source and target (x1, x2)
+		subsets = numba.typed.Dict.empty(key_type=numba.types.string, value_type=list_of_nodes)
+		qualified_nodes = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+		qualified_edges = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
+		annotations = numba.typed.List.empty_list(numba.typed.List.empty_list(interval.interval_type))
+		edges_to_be_added = (numba.typed.List.empty_list(node_type), numba.typed.List.empty_list(node_type), rule_edges[-1])
+
+		satisfaction = True
+		for i, clause in enumerate(clauses):
+			# Unpack clause variables
+			clause_type = clause[0]
+			clause_label = clause[1]
+			clause_variables = clause[2]
+			clause_bnd = clause[3]
+			clause_operator = clause[4]
+
+			# Unpack thresholds
+			# This value is total/available
+			threshold_quantifier_type = thresholds[i][1][1]
+
+			# This is a node clause
+			# The groundings for node clauses are either the target node, neighbors of the target node, or an existing subset of nodes
+			if clause_type == 'node':
+				clause_var_1 = clause_variables[0]
+				subset = get_node_rule_node_clause_subset(clause_var_1, target_node, subsets, neighbors)
+
+				subsets[clause_var_1] = get_qualified_components_node_clause(interpretations_node, subset, clause_label, clause_bnd)
+
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List(subsets[clause_var_1]))
+					qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qn in subsets[clause_var_1]:
+						a.append(interpretations_node[qn].world[clause_label])
+					annotations.append(a)
+
+			# This is an edge clause
+			elif clause_type == 'edge':
+				clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+				subset_source, subset_target = get_node_rule_edge_clause_subset(clause_var_1, clause_var_2, target_node, subsets, neighbors, reverse_neighbors, nodes)
+
+				# Get qualified edges
+				qe = get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, clause_bnd, reverse_graph)
+				subsets[clause_var_1] = qe[0]
+				subsets[clause_var_2] = qe[1]
+
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List.empty_list(node_type))
+					qualified_edges.append(numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qe in numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])):
+						a.append(interpretations_edge[qe].world[clause_label])
+					annotations.append(a)
+			else:
+				# This is a comparison clause
+				# Make sure there is at least one ground atom such that pred-num(x) : [1,1] or pred-num(x,y) : [1,1]
+				# Remember that the predicate in the clause will not contain the "-num" where num is some number.
+				# We have to remove that manually while checking
+				# Steps:
+				# 1. get qualified nodes/edges as well as number associated for first predicate
+				# 2. get qualified nodes/edges as well as number associated for second predicate
+				# 3. if there's no number in steps 1 or 2 return false clause
+				# 4. do comparison with each qualified component from step 1 with each qualified component in step 2
+
+				# It's a node comparison
+				if len(clause_variables) == 2:
+					clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+					subset_1 = get_node_rule_node_clause_subset(clause_var_1, target_node, subsets, neighbors)
+					subset_2 = get_node_rule_node_clause_subset(clause_var_2, target_node, subsets, neighbors)
+
+					# 1, 2
+					qualified_nodes_for_comparison_1, numbers_1 = get_qualified_components_node_comparison_clause(interpretations_node, subset_1, clause_label, clause_bnd)
+					qualified_nodes_for_comparison_2, numbers_2 = get_qualified_components_node_comparison_clause(interpretations_node, subset_2, clause_label, clause_bnd)
+
+				# It's an edge comparison
+				elif len(clause_variables) == 4:
+					clause_var_1_source, clause_var_1_target, clause_var_2_source, clause_var_2_target = clause_variables[0], clause_variables[1], clause_variables[2], clause_variables[3]
+					subset_1_source, subset_1_target = get_node_rule_edge_clause_subset(clause_var_1_source, clause_var_1_target, target_node, subsets, neighbors, reverse_neighbors, nodes)
+					subset_2_source, subset_2_target = get_node_rule_edge_clause_subset(clause_var_2_source, clause_var_2_target, target_node, subsets, neighbors, reverse_neighbors, nodes)
+
+					# 1, 2
+					qualified_nodes_for_comparison_1_source, qualified_nodes_for_comparison_1_target, numbers_1 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_1_source, subset_1_target, clause_label, clause_bnd, reverse_graph)
+					qualified_nodes_for_comparison_2_source, qualified_nodes_for_comparison_2_target, numbers_2 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_2_source, subset_2_target, clause_label, clause_bnd, reverse_graph)
+
+			# Check if thresholds are satisfied
+			# If it's a comparison clause we just need to check if the numbers list is not empty (no threshold support)
+			if clause_type == 'comparison':
+				if len(numbers_1) == 0 or len(numbers_2) == 0:
+					satisfaction = False
+				# Node comparison. Compare stage
+				elif len(clause_variables) == 2:
+					satisfaction, qualified_nodes_1, qualified_nodes_2 = compare_numbers_node_predicate(numbers_1, numbers_2, clause_operator, qualified_nodes_for_comparison_1, qualified_nodes_for_comparison_2)
+
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1] = qualified_nodes_1
+					subsets[clause_var_2] = qualified_nodes_2
+					qualified_comparison_nodes = numba.typed.List(qualified_nodes_1)
+					qualified_comparison_nodes.extend(qualified_nodes_2)
+
+					if atom_trace:
+						qualified_nodes.append(qualified_comparison_nodes)
+						qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qn in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
+				# Edge comparison. Compare stage
+				else:
+					satisfaction, qualified_nodes_1_source, qualified_nodes_1_target, qualified_nodes_2_source, qualified_nodes_2_target = compare_numbers_edge_predicate(numbers_1, numbers_2, clause_operator,
+																																										  qualified_nodes_for_comparison_1_source,
+																																										  qualified_nodes_for_comparison_1_target,
+																																										  qualified_nodes_for_comparison_2_source,
+																																										  qualified_nodes_for_comparison_2_target)
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1_source] = qualified_nodes_1_source
+					subsets[clause_var_1_target] = qualified_nodes_1_target
+					subsets[clause_var_2_source] = qualified_nodes_2_source
+					subsets[clause_var_2_target] = qualified_nodes_2_target
+
+					qualified_comparison_nodes_1 = numba.typed.List(zip(qualified_nodes_1_source, qualified_nodes_1_target))
+					qualified_comparison_nodes_2 = numba.typed.List(zip(qualified_nodes_2_source, qualified_nodes_2_target))
+					qualified_comparison_nodes = numba.typed.List(qualified_comparison_nodes_1)
+					qualified_comparison_nodes.extend(qualified_comparison_nodes_2)
+
+					if atom_trace:
+						qualified_nodes.append(numba.typed.List.empty_list(node_type))
+						qualified_edges.append(qualified_comparison_nodes)
+
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qe in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
+
+			# Non comparison clause
+			else:
+				if threshold_quantifier_type == 'total':
+					if clause_type == 'node':
+						neigh_len = len(subset)
+					else:
+						neigh_len = sum([len(l) for l in subset_target])
+
+				# Available is all neighbors that have a particular label with bound inside [0,1]
+				elif threshold_quantifier_type == 'available':
+					if clause_type == 'node':
+						neigh_len = len(get_qualified_components_node_clause(interpretations_node, subset, clause_label, interval.closed(0,1)))
+					else:
+						neigh_len = len(get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, interval.closed(0,1), reverse_graph)[0])
+
+				qualified_neigh_len = len(subsets[clause_var_1])
+				satisfaction = _satisfies_threshold(neigh_len, qualified_neigh_len, thresholds[i]) and satisfaction
+
+			# Exit loop if even one clause is not satisfied
+			if not satisfaction:
+				break
+
+		if satisfaction:
+			# Collect edges to be added
+			source, target, _ = rule_edges
+
+			# Edges to be added
+			if source != '' and target != '':
+				# Check if edge nodes are target
+				if source == '__target':
+					edges_to_be_added[0].append(target_node)
+				elif source in subsets:
+					edges_to_be_added[0].extend(subsets[source])
+				else:
+					edges_to_be_added[0].append(source)
+
+				if target == '__target':
+					edges_to_be_added[1].append(target_node)
+				elif target in subsets:
+					edges_to_be_added[1].extend(subsets[target])
+				else:
+					edges_to_be_added[1].append(target)
+
+			# node/edge, annotations, qualified nodes, qualified edges, edges to be added
+			applicable_rules.append((target_node, annotations, qualified_nodes, qualified_edges, edges_to_be_added))
+
+	return applicable_rules
+
+
+@numba.njit(cache=False, parallel=True)
+def _ground_edge_rule(rule, interpretations_node, interpretations_edge, nodes, edges, neighbors, reverse_neighbors, atom_trace, reverse_graph, edges_to_skip):
+	# Extract rule params
+	rule_type = rule.get_type()
+	clauses = rule.get_clauses()
+	thresholds = rule.get_thresholds()
+	ann_fn = rule.get_annotation_function()
+	rule_edges = rule.get_edges()
+
+	# We return a list of tuples which specify the target nodes/edges that have made the rule body true
+	applicable_rules = numba.typed.List.empty_list(edge_applicable_rule_type)
+
+	# Return empty list if rule is not node rule
+	if rule_type != 'edge':
+		return applicable_rules
+
+	# Steps
+	# 1. Loop through all nodes and evaluate each clause with that node and check the truth with the thresholds
+	# 2. Inside the clause loop it may be necessary to loop through all nodes/edges while grounding the variables
+	# 3. If the clause is true add the qualified nodes and qualified edges to the atom trace, if on. Break otherwise
+	# 4. After going through all clauses, add to the annotations list all the annotations of the specified subset. These will be passed to the annotation function
+	# 5. Finally, if there are any edges to be added, place them in the list
+
+	for piter in prange(len(edges)):
+		target_edge = edges[piter]
+		if target_edge in edges_to_skip:
+			continue
+		# Initialize dictionary where keys are strings (x1, x2 etc.) and values are lists of qualified neighbors
+		# Keep track of qualified nodes and qualified edges
+		# If it's a node clause update (x1 or x2 etc.) qualified neighbors, if it's an edge clause update the qualified neighbors for the source and target (x1, x2)
+		subsets = numba.typed.Dict.empty(key_type=numba.types.string, value_type=list_of_nodes)
+		qualified_nodes = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+		qualified_edges = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
+		annotations = numba.typed.List.empty_list(numba.typed.List.empty_list(interval.interval_type))
+		edges_to_be_added = (numba.typed.List.empty_list(node_type), numba.typed.List.empty_list(node_type), rule_edges[-1])
+
+		satisfaction = True
+		for i, clause in enumerate(clauses):
+			# Unpack clause variables
+			clause_type = clause[0]
+			clause_label = clause[1]
+			clause_variables = clause[2]
+			clause_bnd = clause[3]
+			clause_operator = clause[4]
+
+			# Unpack thresholds
+			# This value is total/available
+			threshold_quantifier_type = thresholds[i][1][1]
+
+			# This is a node clause
+			# The groundings for node clauses are either the source, target, neighbors of the source node, or an existing subset of nodes
+			if clause_type == 'node':
+				clause_var_1 = clause_variables[0]
+				subset = get_edge_rule_node_clause_subset(clause_var_1, target_edge, subsets, neighbors)
+
+				subsets[clause_var_1] = get_qualified_components_node_clause(interpretations_node, subset, clause_label, clause_bnd)
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List(subsets[clause_var_1]))
+					qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qn in subsets[clause_var_1]:
+						a.append(interpretations_node[qn].world[clause_label])
+					annotations.append(a)
+
+			# This is an edge clause
+			elif clause_type == 'edge':
+				clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+				subset_source, subset_target = get_edge_rule_edge_clause_subset(clause_var_1, clause_var_2, target_edge, subsets, neighbors, reverse_neighbors, nodes)
+
+				# Get qualified edges
+				qe = get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, clause_bnd, reverse_graph)
+				subsets[clause_var_1] = qe[0]
+				subsets[clause_var_2] = qe[1]
+
+				if atom_trace:
+					qualified_nodes.append(numba.typed.List.empty_list(node_type))
+					qualified_edges.append(numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])))
+
+				# Add annotations if necessary
+				if ann_fn != '':
+					a = numba.typed.List.empty_list(interval.interval_type)
+					for qe in numba.typed.List(zip(subsets[clause_var_1], subsets[clause_var_2])):
+						a.append(interpretations_edge[qe].world[clause_label])
+					annotations.append(a)
+
+			else:
+				# This is a comparison clause
+				# Make sure there is at least one ground atom such that pred-num(x) : [1,1] or pred-num(x,y) : [1,1]
+				# Remember that the predicate in the clause will not contain the "-num" where num is some number.
+				# We have to remove that manually while checking
+				# Steps:
+				# 1. get qualified nodes/edges as well as number associated for first predicate
+				# 2. get qualified nodes/edges as well as number associated for second predicate
+				# 3. if there's no number in steps 1 or 2 return false clause
+				# 4. do comparison with each qualified component from step 1 with each qualified component in step 2
+
+				# It's a node comparison
+				if len(clause_variables) == 2:
+					clause_var_1, clause_var_2 = clause_variables[0], clause_variables[1]
+					subset_1 = get_edge_rule_node_clause_subset(clause_var_1, target_edge, subsets, neighbors)
+					subset_2 = get_edge_rule_node_clause_subset(clause_var_2, target_edge, subsets, neighbors)
+
+					# 1, 2
+					qualified_nodes_for_comparison_1, numbers_1 = get_qualified_components_node_comparison_clause(interpretations_node, subset_1, clause_label, clause_bnd)
+					qualified_nodes_for_comparison_2, numbers_2 = get_qualified_components_node_comparison_clause(interpretations_node, subset_2, clause_label, clause_bnd)
+
+				# It's an edge comparison
+				elif len(clause_variables) == 4:
+					clause_var_1_source, clause_var_1_target, clause_var_2_source, clause_var_2_target = clause_variables[0], clause_variables[1], clause_variables[2], clause_variables[3]
+					subset_1_source, subset_1_target = get_edge_rule_edge_clause_subset(clause_var_1_source, clause_var_1_target, target_edge, subsets, neighbors, reverse_neighbors, nodes)
+					subset_2_source, subset_2_target = get_edge_rule_edge_clause_subset(clause_var_2_source, clause_var_2_target, target_edge, subsets, neighbors, reverse_neighbors, nodes)
+
+					# 1, 2
+					qualified_nodes_for_comparison_1_source, qualified_nodes_for_comparison_1_target, numbers_1 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_1_source, subset_1_target, clause_label, clause_bnd, reverse_graph)
+					qualified_nodes_for_comparison_2_source, qualified_nodes_for_comparison_2_target, numbers_2 = get_qualified_components_edge_comparison_clause(interpretations_edge, subset_2_source, subset_2_target, clause_label, clause_bnd, reverse_graph)
+
+			# Check if thresholds are satisfied
+			# If it's a comparison clause we just need to check if the numbers list is not empty (no threshold support)
+			if clause_type == 'comparison':
+				if len(numbers_1) == 0 or len(numbers_2) == 0:
+					satisfaction = False
+				# Node comparison. Compare stage
+				elif len(clause_variables) == 2:
+					satisfaction, qualified_nodes_1, qualified_nodes_2 = compare_numbers_node_predicate(numbers_1, numbers_2, clause_operator, qualified_nodes_for_comparison_1, qualified_nodes_for_comparison_2)
+
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1] = qualified_nodes_1
+					subsets[clause_var_2] = qualified_nodes_2
+					qualified_comparison_nodes = numba.typed.List(qualified_nodes_1)
+					qualified_comparison_nodes.extend(qualified_nodes_2)
+
+					if atom_trace:
+						qualified_nodes.append(qualified_comparison_nodes)
+						qualified_edges.append(numba.typed.List.empty_list(edge_type))
+
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qn in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
+				# Edge comparison. Compare stage
+				else:
+					satisfaction, qualified_nodes_1_source, qualified_nodes_1_target, qualified_nodes_2_source, qualified_nodes_2_target = compare_numbers_edge_predicate(numbers_1, numbers_2, clause_operator,
+																																										  qualified_nodes_for_comparison_1_source,
+																																										  qualified_nodes_for_comparison_1_target,
+																																										  qualified_nodes_for_comparison_2_source,
+																																										  qualified_nodes_for_comparison_2_target)
+					# Update subsets with final qualified nodes
+					subsets[clause_var_1_source] = qualified_nodes_1_source
+					subsets[clause_var_1_target] = qualified_nodes_1_target
+					subsets[clause_var_2_source] = qualified_nodes_2_source
+					subsets[clause_var_2_target] = qualified_nodes_2_target
+
+					qualified_comparison_nodes_1 = numba.typed.List(zip(qualified_nodes_1_source, qualified_nodes_1_target))
+					qualified_comparison_nodes_2 = numba.typed.List(zip(qualified_nodes_2_source, qualified_nodes_2_target))
+					qualified_comparison_nodes = numba.typed.List(qualified_comparison_nodes_1)
+					qualified_comparison_nodes.extend(qualified_comparison_nodes_2)
+
+					if atom_trace:
+						qualified_nodes.append(numba.typed.List.empty_list(node_type))
+						qualified_edges.append(qualified_comparison_nodes)
+
+					# Add annotations for comparison clause. For now, we don't distinguish between LHS and RHS annotations
+					if ann_fn != '':
+						a = numba.typed.List.empty_list(interval.interval_type)
+						for qe in qualified_comparison_nodes:
+							a.append(interval.closed(1, 1))
+						annotations.append(a)
+			
+			# Non comparison clause
+			else:
+				if threshold_quantifier_type == 'total':
+					if clause_type == 'node':
+						neigh_len = len(subset)
+					else:
+						neigh_len = sum([len(l) for l in subset_target])
+	
+				# Available is all neighbors that have a particular label with bound inside [0,1]
+				elif threshold_quantifier_type == 'available':
+					if clause_type == 'node':
+						neigh_len = len(get_qualified_components_node_clause(interpretations_node, subset, clause_label, interval.closed(0, 1)))
+					else:
+						neigh_len = len(get_qualified_components_edge_clause(interpretations_edge, subset_source, subset_target, clause_label, interval.closed(0, 1), reverse_graph)[0])
+	
+				qualified_neigh_len = len(subsets[clause_var_1])
+				satisfaction = _satisfies_threshold(neigh_len, qualified_neigh_len, thresholds[i]) and satisfaction
+			
+			# Exit loop if even one clause is not satisfied
+			if not satisfaction:
+				break
+
+		# Here we are done going through each clause of the rule
+		# If all clauses we're satisfied, proceed to collect annotations and prepare edges to be added
+		if satisfaction:
+			# Collect edges to be added
+			source, target, _ = rule_edges
+
+			# Edges to be added
+			if source != '' and target != '':
+				# Check if edge nodes are source/target
+				if source == '__source':
+					edges_to_be_added[0].append(target_edge[0])
+				elif source == '__target':
+					edges_to_be_added[0].append(target_edge[1])
+				elif source in subsets:
+					edges_to_be_added[0].extend(subsets[source])
+				else:
+					edges_to_be_added[0].append(source)
+
+				if target == '__source':
+					edges_to_be_added[1].append(target_edge[0])
+				elif target == '__target':
+					edges_to_be_added[1].append(target_edge[1])
+				elif target in subsets:
+					edges_to_be_added[1].extend(subsets[target])
+				else:
+					edges_to_be_added[1].append(target)
+
+			# node/edge, annotations, qualified nodes, qualified edges, edges to be added
+			applicable_rules.append((target_edge, annotations, qualified_nodes, qualified_edges, edges_to_be_added))
+
+	return applicable_rules
+
+
+@numba.njit(cache=False)
+def get_node_rule_node_clause_subset(clause_var_1, target_node, subsets, neighbors):
+	# The groundings for node clauses are either the target node, neighbors of the target node, or an existing subset of nodes
+	if clause_var_1 == '__target':
+		subset = numba.typed.List([target_node])
+	else:
+		subset = neighbors[target_node] if clause_var_1 not in subsets else subsets[clause_var_1]
+	return subset
+
+
+@numba.njit(cache=False)
+def get_node_rule_edge_clause_subset(clause_var_1, clause_var_2, target_node, subsets, neighbors, reverse_neighbors, nodes):
+	# There are 5 cases for predicate(Y,Z):
+	# 1. Either one or both of Y, Z are the target node
+	# 2. Both predicate variables Y and Z have not been encountered before
+	# 3. The source variable Y has not been encountered before but the target variable Z has
+	# 4. The target variable Z has not been encountered before but the source variable Y has
+	# 5. Both predicate variables Y and Z have been encountered before
+
+	# Case 1:
+	# Check if 1st variable or 1st and 2nd variables are the target
+	if clause_var_1 == '__target':
+		subset_source = numba.typed.List([target_node])
+
+		# If both variables are the same
+		if clause_var_2 == '__target':
+			subset_target = numba.typed.List([numba.typed.List([target_node])])
+		elif clause_var_2 in subsets:
+			subset_target = numba.typed.List([subsets[clause_var_2]])
+		else:
+			subset_target = numba.typed.List([neighbors[target_node]])
+
+	# Check if 2nd variable is the target (this means 1st variable isn't the target)
+	elif clause_var_2 == '__target':
+		subset_source = reverse_neighbors[target_node] if clause_var_1 not in subsets else subsets[clause_var_1]
+		subset_target = numba.typed.List([numba.typed.List([target_node]) for _ in subset_source])
+
+	# Case 2:
+	# We replace Y by all nodes and Z by the neighbors of each of these nodes
+	elif clause_var_1 not in subsets and clause_var_2 not in subsets:
+		subset_source = numba.typed.List(nodes)
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 3:
+	# We replace Y by the sources of Z
+	elif clause_var_1 not in subsets and clause_var_2 in subsets:
+		subset_source = numba.typed.List.empty_list(node_type)
+		subset_target = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+
+		for n in subsets[clause_var_2]:
+			sources = reverse_neighbors[n]
+			for source in sources:
+				subset_source.append(source)
+				subset_target.append(numba.typed.List([n]))
+
+	# Case 4:
+	# We replace Z by the neighbors of Y
+	elif clause_var_1 in subsets and clause_var_2 not in subsets:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 5:
+	else:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([subsets[clause_var_2] for _ in subset_source])
+
+	return subset_source, subset_target
+
+
+@numba.njit(cache=False)
+def get_edge_rule_node_clause_subset(clause_var_1, target_edge, subsets, neighbors):
+	# The groundings for node clauses are either the source, target, neighbors of the source node, or an existing subset of nodes
+	if clause_var_1 == '__source':
+		subset = numba.typed.List([target_edge[0]])
+	elif clause_var_1 == '__target':
+		subset = numba.typed.List([target_edge[1]])
+	else:
+		subset = neighbors[target_edge[0]] if clause_var_1 not in subsets else subsets[clause_var_1]
+	return subset
+
+
+@numba.njit(cache=False)
+def get_edge_rule_edge_clause_subset(clause_var_1, clause_var_2, target_edge, subsets, neighbors, reverse_neighbors, nodes):
+	# There are 5 cases for predicate(Y,Z):
+	# 1. Either one or both of Y, Z are the source or target node
+	# 2. Both predicate variables Y and Z have not been encountered before
+	# 3. The source variable Y has not been encountered before but the target variable Z has
+	# 4. The target variable Z has not been encountered before but the source variable Y has
+	# 5. Both predicate variables Y and Z have been encountered before
+	# Case 1:
+	# Check if 1st variable is the source
+	if clause_var_1 == '__source':
+		subset_source = numba.typed.List([target_edge[0]])
+
+		# If 2nd variable is source/target/something else
+		if clause_var_2 == '__source':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[0]])])
+		elif clause_var_2 == '__target':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[1]])])
+		elif clause_var_2 in subsets:
+			subset_target = numba.typed.List([subsets[clause_var_2]])
+		else:
+			subset_target = numba.typed.List([neighbors[target_edge[0]]])
+
+	# if 1st variable is the target
+	elif clause_var_1 == '__target':
+		subset_source = numba.typed.List([target_edge[1]])
+
+		# if 2nd variable is source/target/something else
+		if clause_var_2 == '__source':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[0]])])
+		elif clause_var_2 == '__target':
+			subset_target = numba.typed.List([numba.typed.List([target_edge[1]])])
+		elif clause_var_2 in subsets:
+			subset_target = numba.typed.List([subsets[clause_var_2]])
+		else:
+			subset_target = numba.typed.List([neighbors[target_edge[1]]])
+
+	# Handle the cases where the 2nd variable is source/target but the 1st is something else (cannot be source/target)
+	elif clause_var_2 == '__source':
+		subset_source = reverse_neighbors[target_edge[0]] if clause_var_1 not in subsets else subsets[clause_var_1]
+		subset_target = numba.typed.List([numba.typed.List([target_edge[0]]) for _ in subset_source])
+
+	elif clause_var_2 == '__target':
+		subset_source = reverse_neighbors[target_edge[1]] if clause_var_1 not in subsets else subsets[clause_var_1]
+		subset_target = numba.typed.List([numba.typed.List([target_edge[1]]) for _ in subset_source])
+
+	# Case 2:
+	# We replace Y by all nodes and Z by the neighbors of each of these nodes
+	elif clause_var_1 not in subsets and clause_var_2 not in subsets:
+		subset_source = numba.typed.List(nodes)
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 3:
+	# We replace Y by the sources of Z
+	elif clause_var_1 not in subsets and clause_var_2 in subsets:
+		subset_source = numba.typed.List.empty_list(node_type)
+		subset_target = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+
+		for n in subsets[clause_var_2]:
+			sources = reverse_neighbors[n]
+			for source in sources:
+				subset_source.append(source)
+				subset_target.append(numba.typed.List([n]))
+
+	# Case 4:
+	# We replace Z by the neighbors of Y
+	elif clause_var_1 in subsets and clause_var_2 not in subsets:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([neighbors[n] for n in subset_source])
+
+	# Case 5:
+	else:
+		subset_source = subsets[clause_var_1]
+		subset_target = numba.typed.List([subsets[clause_var_2] for _ in subset_source])
+
+	return subset_source, subset_target
+
+
+@numba.njit(cache=False)
+def get_qualified_components_node_clause(interpretations_node, candidates, l, bnd):
+	# Get all the qualified neighbors for a particular clause
+	qualified_nodes = numba.typed.List.empty_list(node_type)
+	for n in candidates:
+		if is_satisfied_node(interpretations_node, n, (l, bnd)):
+			qualified_nodes.append(n)
+
+	return qualified_nodes
+
+
+@numba.njit(cache=False)
+def get_qualified_components_node_comparison_clause(interpretations_node, candidates, l, bnd):
+	# Get all the qualified neighbors for a particular comparison clause and return them along with the number associated
+	qualified_nodes = numba.typed.List.empty_list(node_type)
+	qualified_nodes_numbers = numba.typed.List.empty_list(numba.types.float64)
+	for n in candidates:
+		result, number = is_satisfied_node_comparison(interpretations_node, n, (l, bnd))
+		if result:
+			qualified_nodes.append(n)
+			qualified_nodes_numbers.append(number)
+
+	return qualified_nodes, qualified_nodes_numbers
+
+
+@numba.njit(cache=False)
+def get_qualified_components_edge_clause(interpretations_edge, candidates_source, candidates_target, l, bnd, reverse_graph):
+	# Get all the qualified sources and targets for a particular clause
+	qualified_nodes_source = numba.typed.List.empty_list(node_type)
+	qualified_nodes_target = numba.typed.List.empty_list(node_type)
+	for i, source in enumerate(candidates_source):
+		for target in candidates_target[i]:
+			edge = (source, target) if not reverse_graph else (target, source)
+			if is_satisfied_edge(interpretations_edge, edge, (l, bnd)):
+				qualified_nodes_source.append(source)
+				qualified_nodes_target.append(target)
+
+	return qualified_nodes_source, qualified_nodes_target
+
+
+@numba.njit(cache=False)
+def get_qualified_components_edge_comparison_clause(interpretations_edge, candidates_source, candidates_target, l, bnd, reverse_graph):
+	# Get all the qualified sources and targets for a particular clause
+	qualified_nodes_source = numba.typed.List.empty_list(node_type)
+	qualified_nodes_target = numba.typed.List.empty_list(node_type)
+	qualified_edges_numbers = numba.typed.List.empty_list(numba.types.float64)
+	for i, source in enumerate(candidates_source):
+		for target in candidates_target[i]:
+			edge = (source, target) if not reverse_graph else (target, source)
+			result, number = is_satisfied_edge_comparison(interpretations_edge, edge, (l, bnd))
+			if result:
+				qualified_nodes_source.append(source)
+				qualified_nodes_target.append(target)
+				qualified_edges_numbers.append(number)
+
+	return qualified_nodes_source, qualified_nodes_target, qualified_edges_numbers
+
+
+@numba.njit(cache=False)
+def compare_numbers_node_predicate(numbers_1, numbers_2, op, qualified_nodes_1, qualified_nodes_2):
+	result = False
+	final_qualified_nodes_1 = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_2 = numba.typed.List.empty_list(node_type)
+	for i in range(len(numbers_1)):
+		for j in range(len(numbers_2)):
+			if op == '<':
+				if numbers_1[i] < numbers_2[j]:
+					result = True
+			elif op == '<=':
+				if numbers_1[i] <= numbers_2[j]:
+					result = True
+			elif op == '>':
+				if numbers_1[i] > numbers_2[j]:
+					result = True
+			elif op == '>=':
+				if numbers_1[i] >= numbers_2[j]:
+					result = True
+			elif op == '==':
+				if numbers_1[i] == numbers_2[j]:
+					result = True
+			elif op == '!=':
+				if numbers_1[i] != numbers_2[j]:
+					result = True
+
+			if result:
+				final_qualified_nodes_1.append(qualified_nodes_1[i])
+				final_qualified_nodes_2.append(qualified_nodes_2[j])
+	return result, final_qualified_nodes_1, final_qualified_nodes_2
+
+
+@numba.njit(cache=False)
+def compare_numbers_edge_predicate(numbers_1, numbers_2, op, qualified_nodes_1a, qualified_nodes_1b, qualified_nodes_2a, qualified_nodes_2b):
+	result = False
+	final_qualified_nodes_1a = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_1b = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_2a = numba.typed.List.empty_list(node_type)
+	final_qualified_nodes_2b = numba.typed.List.empty_list(node_type)
+	for i in range(len(numbers_1)):
+		for j in range(len(numbers_2)):
+			if op == '<':
+				if numbers_1[i] < numbers_2[j]:
+					result = True
+			elif op == '<=':
+				if numbers_1[i] <= numbers_2[j]:
+					result = True
+			elif op == '>':
+				if numbers_1[i] > numbers_2[j]:
+					result = True
+			elif op == '>=':
+				if numbers_1[i] >= numbers_2[j]:
+					result = True
+			elif op == '==':
+				if numbers_1[i] == numbers_2[j]:
+					result = True
+			elif op == '!=':
+				if numbers_1[i] != numbers_2[j]:
+					result = True
+
+			if result:
+				final_qualified_nodes_1a.append(qualified_nodes_1a[i])
+				final_qualified_nodes_1b.append(qualified_nodes_1b[i])
+				final_qualified_nodes_2a.append(qualified_nodes_2a[j])
+				final_qualified_nodes_2b.append(qualified_nodes_2b[j])
+	return result, final_qualified_nodes_1a, final_qualified_nodes_1b, final_qualified_nodes_2a, final_qualified_nodes_2b
+
+
+@numba.njit(cache=False)
+def _satisfies_threshold(num_neigh, num_qualified_component, threshold):
+	# Checks if qualified neighbors satisfy threshold. This is for one clause
+	if threshold[1][0]=='number':
+		if threshold[0]=='greater_equal':
+			result = True if num_qualified_component >= threshold[2] else False
+		elif threshold[0]=='greater':
+			result = True if num_qualified_component > threshold[2] else False
+		elif threshold[0]=='less_equal':
+			result = True if num_qualified_component <= threshold[2] else False
+		elif threshold[0]=='less':
+			result = True if num_qualified_component < threshold[2] else False
+		elif threshold[0]=='equal':
+			result = True if num_qualified_component == threshold[2] else False
+
+	elif threshold[1][0]=='percent':
+		if num_neigh==0:
+			result = False
+		elif threshold[0]=='greater_equal':
+			result = True if num_qualified_component/num_neigh >= threshold[2]*0.01 else False
+		elif threshold[0]=='greater':
+			result = True if num_qualified_component/num_neigh > threshold[2]*0.01 else False
+		elif threshold[0]=='less_equal':
+			result = True if num_qualified_component/num_neigh <= threshold[2]*0.01 else False
+		elif threshold[0]=='less':
+			result = True if num_qualified_component/num_neigh < threshold[2]*0.01 else False
+		elif threshold[0]=='equal':
+			result = True if num_qualified_component/num_neigh == threshold[2]*0.01 else False
+
+	return result
+
+
+@numba.njit(cache=False)
+def _update_node(interpretations, comp, na, ipl, rule_trace, fp_cnt, t_cnt, static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_trace, idx, facts_to_be_applied_trace, rule_trace_atoms, store_interpretation_changes, mode, override=False):
+	updated = False
+	# This is to prevent a key error in case the label is a specific label
+	try:
+		world = interpretations[comp]
+		l, bnd = na
+		updated_bnds = numba.typed.List.empty_list(interval.interval_type)
+
+		# Add label to world if it is not there
+		if l not in world.world:
+			world.world[l] = interval.closed(0, 1)
+
+		# Check if update is necessary with previous bnd
+		prev_bnd = world.world[l].copy()
+
+		# override will not check for inconsistencies
+		if override:
+			world.world[l].set_lower_upper(bnd.lower, bnd.upper)
+		else:
+			world.update(l, bnd)
+		world.world[l].set_static(static)
+		if world.world[l]!=prev_bnd:
+			updated = True
+			updated_bnds.append(world.world[l])
+
+			# Add to rule trace if update happened and add to atom trace if necessary
+			if (save_graph_attributes_to_rule_trace or not mode=='graph-attribute-fact') and store_interpretation_changes:
+				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, l, world.world[l].copy()))
+				if atom_trace:
+					# Mode can be fact or rule, updation of trace will happen accordingly
+					if mode=='fact' or mode=='graph-attribute-fact':
+						qn = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+						qe = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
+						name = facts_to_be_applied_trace[idx]
+						_update_rule_trace(rule_trace_atoms, qn, qe, prev_bnd, name)
+					elif mode=='rule':
+						qn, qe, name = rules_to_be_applied_trace[idx]
+						_update_rule_trace(rule_trace_atoms, qn, qe, prev_bnd, name)
+
+		# Update complement of predicate (if exists) based on new knowledge of predicate
+		if updated:
+			ip_update_cnt = 0
+			for p1, p2 in ipl:
+				if p1==l:
+					if atom_trace:
+						_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p2], f'IPL: {l.get_value()}')
+					lower = max(world.world[p2].lower, 1 - world.world[p1].upper)
+					upper = min(world.world[p2].upper, 1 - world.world[p1].lower)
+					world.world[p2].set_lower_upper(lower, upper)
+					world.world[p2].set_static(static)
+					ip_update_cnt += 1
+					updated_bnds.append(world.world[p2])
+					if store_interpretation_changes:
+						rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p2, interval.closed(lower, upper)))
+				if p2==l:
+					if atom_trace:
+						_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p1], f'IPL: {l.get_value()}')
+					lower = max(world.world[p1].lower, 1 - world.world[p2].upper)
+					upper = min(world.world[p1].upper, 1 - world.world[p2].lower)
+					world.world[p1].set_lower_upper(lower, upper)
+					world.world[p1].set_static(static)
+					ip_update_cnt += 1
+					updated_bnds.append(world.world[p1])
+					if store_interpretation_changes:
+						rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p1, interval.closed(lower, upper)))
+
+		# Gather convergence data
+		change = 0
+		if updated:
+			# Find out if it has changed from previous interp
+			current_bnd = world.world[l]
+			prev_t_bnd = interval.closed(world.world[l].prev_lower, world.world[l].prev_upper)
+			if current_bnd != prev_t_bnd:
+				if convergence_mode=='delta_bound':
+					for i in updated_bnds:
+						lower_delta = abs(i.lower-prev_t_bnd.lower)
+						upper_delta = abs(i.upper-prev_t_bnd.upper)
+						max_delta = max(lower_delta, upper_delta)
+						change = max(change, max_delta)
+				else:
+					change = 1 + ip_update_cnt
+
+		return (updated, change)
+
+	except:
+		return (False, 0)
+
+
+@numba.njit(cache=False)
+def _update_edge(interpretations, comp, na, ipl, rule_trace, fp_cnt, t_cnt, static, convergence_mode, atom_trace, save_graph_attributes_to_rule_trace, rules_to_be_applied_trace, idx, facts_to_be_applied_trace, rule_trace_atoms, store_interpretation_changes, mode, override=False):
+	updated = False
+	# This is to prevent a key error in case the label is a specific label
+	try:
+		world = interpretations[comp]
+		l, bnd = na
+		updated_bnds = numba.typed.List.empty_list(interval.interval_type)
+
+		# Add label to world if it is not there
+		if l not in world.world:
+			world.world[l] = interval.closed(0, 1)
+
+		# Check if update is necessary with previous bnd
+		prev_bnd = world.world[l].copy()
+
+		# override will not check for inconsistencies
+		if override:
+			world.world[l].set_lower_upper(bnd.lower, bnd.upper)
+		else:
+			world.update(l, bnd)
+		world.world[l].set_static(static)
+		if world.world[l]!=prev_bnd:
+			updated = True
+			updated_bnds.append(world.world[l])
+
+			# Add to rule trace if update happened and add to atom trace if necessary
+			if (save_graph_attributes_to_rule_trace or not mode=='graph-attribute-fact') and store_interpretation_changes:
+				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, l, world.world[l].copy()))
+				if atom_trace:
+					# Mode can be fact or rule, updation of trace will happen accordingly
+					if mode=='fact' or mode=='graph-attribute-fact':
+						qn = numba.typed.List.empty_list(numba.typed.List.empty_list(node_type))
+						qe = numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type))
+						name = facts_to_be_applied_trace[idx]
+						_update_rule_trace(rule_trace_atoms, qn, qe, prev_bnd, name)
+					elif mode=='rule':
+						qn, qe, name = rules_to_be_applied_trace[idx]
+						_update_rule_trace(rule_trace_atoms, qn, qe, prev_bnd, name)
+
+		# Update complement of predicate (if exists) based on new knowledge of predicate
+		if updated:
+			ip_update_cnt = 0
+			for p1, p2 in ipl:
+				if p1==l:
+					if atom_trace:
+						_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p2], f'IPL: {l.get_value()}')
+					lower = max(world.world[p2].lower, 1 - world.world[p1].upper)
+					upper = min(world.world[p2].upper, 1 - world.world[p1].lower)
+					world.world[p2].set_lower_upper(lower, upper)
+					world.world[p2].set_static(static)
+					ip_update_cnt += 1
+					updated_bnds.append(world.world[p2])
+					if store_interpretation_changes:
+						rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p2, interval.closed(lower, upper)))
+				if p2==l:
+					if atom_trace:
+						_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p1], f'IPL: {l.get_value()}')
+					lower = max(world.world[p1].lower, 1 - world.world[p2].upper)
+					upper = min(world.world[p1].upper, 1 - world.world[p2].lower)
+					world.world[p1].set_lower_upper(lower, upper)
+					world.world[p1].set_static(static)
+					ip_update_cnt += 1
+					updated_bnds.append(world.world[p2])
+					if store_interpretation_changes:
+						rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p1, interval.closed(lower, upper)))
+	
+		# Gather convergence data
+		change = 0
+		if updated:
+			# Find out if it has changed from previous interp
+			current_bnd = world.world[l]
+			prev_t_bnd = interval.closed(world.world[l].prev_lower, world.world[l].prev_upper)
+			if current_bnd != prev_t_bnd:
+				if convergence_mode=='delta_bound':
+					for i in updated_bnds:
+						lower_delta = abs(i.lower-prev_t_bnd.lower)
+						upper_delta = abs(i.upper-prev_t_bnd.upper)
+						max_delta = max(lower_delta, upper_delta)
+						change = max(change, max_delta)
+				else:
+					change = 1 + ip_update_cnt
+		
+		return (updated, change)
+	except:
+		return (False, 0)
+
+
+@numba.njit(cache=False)
+def _update_rule_trace(rule_trace, qn, qe, prev_bnd, name):
+	rule_trace.append((qn, qe, prev_bnd.copy(), name))
+	
+
+@numba.njit(cache=False)
+def are_satisfied_node(interpretations, comp, nas):
+	result = True
+	for (label, interval) in nas:
+		result = result and is_satisfied_node(interpretations, comp, (label, interval))
+	return result
+
+
+@numba.njit(cache=False)
+def is_satisfied_node(interpretations, comp, na):
+	result = False
+	if (not (na[0] is None or na[1] is None)):
+		# This is to prevent a key error in case the label is a specific label
+		try:
+			world = interpretations[comp]
+			result = world.is_satisfied(na[0], na[1])
+		except:
+			result = False
+	else:
+		result = True
+	return result
+
+
+@numba.njit(cache=False)
+def is_satisfied_node_comparison(interpretations, comp, na):
+	result = False
+	number = 0
+	l, bnd = na
+	l_str = l.value
+
+	if not (l is None or bnd is None):
+		# This is to prevent a key error in case the label is a specific label
+		try:
+			world = interpretations[comp]
+			for world_l in world.world.keys():
+				world_l_str = world_l.value
+				if l_str in world_l_str and world_l_str[len(l_str)+1:].replace('.', '').replace('-', '').isdigit():
+					# The label is contained in the world
+					result = world.is_satisfied(world_l, na[1])
+					# Find the suffix number
+					number = str_to_float(world_l_str[len(l_str)+1:])
+					break
+
+		except:
+			result = False
+	else:
+		result = True
+	return result, number
+
+
+@numba.njit(cache=False)
+def are_satisfied_edge(interpretations, comp, nas):
+	result = True
+	for (label, interval) in nas:
+		result = result and is_satisfied_edge(interpretations, comp, (label, interval))
+	return result
+
+
+@numba.njit(cache=False)
+def is_satisfied_edge(interpretations, comp, na):
+	result = False
+	if (not (na[0] is None or na[1] is None)):
+		# This is to prevent a key error in case the label is a specific label
+		try:
+			world = interpretations[comp]
+			result = world.is_satisfied(na[0], na[1])
+		except:
+			result = False
+	else:
+		result = True
+	return result
+
+
+@numba.njit(cache=False)
+def is_satisfied_edge_comparison(interpretations, comp, na):
+	result = False
+	number = 0
+	l, bnd = na
+	l_str = l.value
+
+	if not (l is None or bnd is None):
+		# This is to prevent a key error in case the label is a specific label
+		try:
+			world = interpretations[comp]
+			for world_l in world.world.keys():
+				world_l_str = world_l.value
+				if l_str in world_l_str and world_l_str[len(l_str)+1:].replace('.', '').replace('-', '').isdigit():
+					# The label is contained in the world
+					result = world.is_satisfied(world_l, na[1])
+					# Find the suffix number
+					number = str_to_float(world_l_str[len(l_str)+1:])
+					break
+
+		except:
+			result = False
+	else:
+		result = True
+	return result, number
+
+
+@numba.njit(cache=False)
+def annotate(annotation_functions, rule, annotations, weights):
+	func_name = rule.get_annotation_function()
+	if func_name == '':
+		return rule.get_bnd().lower, rule.get_bnd().upper
+	else:
+		with numba.objmode(annotation='Tuple((float64, float64))'):
+			for func in annotation_functions:
+				if func.__name__ == func_name:
+					annotation = func(annotations, weights)
+		return annotation
+
+
+@numba.njit(cache=False)
+def check_consistent_node(interpretations, comp, na):
+	world = interpretations[comp]
+	if na[0] in world.world:
+		bnd = world.world[na[0]]
+	else:
+		bnd = interval.closed(0, 1)
+	if (na[1].lower > bnd.upper) or (bnd.lower > na[1].upper):
+		return False
+	else:
+		return True
+
+
+@numba.njit(cache=False)
+def check_consistent_edge(interpretations, comp, na):
+	world = interpretations[comp]
+	if na[0] in world.world:
+		bnd = world.world[na[0]]
+	else:
+		bnd = interval.closed(0, 1)
+	if (na[1].lower > bnd.upper) or (bnd.lower > na[1].upper):
+		return False
+	else:
+		return True
+
+
+@numba.njit(cache=False)
+def resolve_inconsistency_node(interpretations, comp, na, ipl, t_cnt, fp_cnt, atom_trace, rule_trace, rule_trace_atoms, store_interpretation_changes):
+	world = interpretations[comp]
+	if store_interpretation_changes:
+		rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, na[0], interval.closed(0,1)))
+		if atom_trace:
+			_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[na[0]], 'Inconsistency')
+	# Resolve inconsistency and set static
+	world.world[na[0]].set_lower_upper(0, 1)
+	world.world[na[0]].set_static(True)
+	for p1, p2 in ipl:
+		if p1==na[0]:
+			if atom_trace:
+				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p2], 'Inconsistency')
+			world.world[p2].set_lower_upper(0, 1)
+			world.world[p2].set_static(True)
+			if store_interpretation_changes:
+				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p2, interval.closed(0,1)))
+
+		if p2==na[0]:
+			if atom_trace:
+				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), world.world[p1], 'Inconsistency')
+			world.world[p1].set_lower_upper(0, 1)
+			world.world[p1].set_static(True)
+			if store_interpretation_changes:
+				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p1, interval.closed(0,1)))
+	# Add inconsistent predicates to a list 
+
+
+@numba.njit(cache=False)
+def resolve_inconsistency_edge(interpretations, comp, na, ipl, t_cnt, fp_cnt, atom_trace, rule_trace, rule_trace_atoms, store_interpretation_changes):
+	w = interpretations[comp]
+	if store_interpretation_changes:
+		rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, na[0], interval.closed(0,1)))
+		if atom_trace:
+			_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), w.world[na[0]], 'Inconsistency')
+	# Resolve inconsistency and set static
+	w.world[na[0]].set_lower_upper(0, 1)
+	w.world[na[0]].set_static(True)
+	for p1, p2 in ipl:
+		if p1==na[0]:
+			if atom_trace:
+				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), w.world[p2], 'Inconsistency')
+			w.world[p2].set_lower_upper(0, 1)
+			w.world[p2].set_static(True)
+			if store_interpretation_changes:
+				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p2, interval.closed(0,1)))
+
+		if p2==na[0]:
+			if atom_trace:
+				_update_rule_trace(rule_trace_atoms, numba.typed.List.empty_list(numba.typed.List.empty_list(node_type)), numba.typed.List.empty_list(numba.typed.List.empty_list(edge_type)), w.world[p1], 'Inconsistency')
+			w.world[p1].set_lower_upper(0, 1)
+			w.world[p1].set_static(True)
+			if store_interpretation_changes:
+				rule_trace.append((numba.types.uint16(t_cnt), numba.types.uint16(fp_cnt), comp, p1, interval.closed(0,1)))
+
+
+@numba.njit(cache=False)
+def _add_node(node, neighbors, reverse_neighbors, nodes, interpretations_node):
+	nodes.append(node)
+	neighbors[node] = numba.typed.List.empty_list(node_type)
+	reverse_neighbors[node] = numba.typed.List.empty_list(node_type)
+	interpretations_node[node] = world.World(numba.typed.List.empty_list(label.label_type))
+
+
+@numba.njit(cache=False)
+def _add_edge(source, target, neighbors, reverse_neighbors, nodes, edges, l, interpretations_node, interpretations_edge):
+	# If not a node, add to list of nodes and initialize neighbors
+	if source not in nodes:
+		_add_node(source, neighbors, reverse_neighbors, nodes, interpretations_node)
+
+	if target not in nodes:
+		_add_node(target, neighbors, reverse_neighbors, nodes, interpretations_node)
+
+	# Make sure edge doesn't already exist
+	# Make sure, if l=='', not to add the label
+	# Make sure, if edge exists, that we don't override the l label if it exists
+	edge = (source, target)
+	new_edge = False
+	if edge not in edges:
+		new_edge = True
+		edges.append(edge)
+		neighbors[source].append(target)
+		reverse_neighbors[target].append(source)
+		if l.value!='':
+			interpretations_edge[edge] = world.World(numba.typed.List([l]))
+		else:
+			interpretations_edge[edge] = world.World(numba.typed.List.empty_list(label.label_type))
+	else:
+		if l not in interpretations_edge[edge].world and l.value!='':
+			new_edge = True
+			interpretations_edge[edge].world[l] = interval.closed(0, 1)
+
+	return (edge, new_edge)
+
+
+@numba.njit(cache=False)
+def _add_edges(sources, targets, neighbors, reverse_neighbors, nodes, edges, l, interpretations_node, interpretations_edge):
+	changes = 0
+	edges_added = numba.typed.List.empty_list(edge_type)
+	for source in sources:
+		for target in targets:
+			edge, new_edge = _add_edge(source, target, neighbors, reverse_neighbors, nodes, edges, l, interpretations_node, interpretations_edge)
+			edges_added.append(edge)
+			changes = changes+1 if new_edge else changes
+	return edges_added, changes
+
+
+@numba.njit(cache=False)
+def _delete_edge(edge, neighbors, reverse_neighbors, edges, interpretations_edge):
+	source, target = edge
+	edges.remove(edge)
+	del interpretations_edge[edge]
+	neighbors[source].remove(target)
+	reverse_neighbors[target].remove(source)
+
+
+@numba.njit(cache=False)
+def float_to_str(value):
+	number = int(value)
+	decimal = int(value % 1 * 1000)
+	float_str = f'{number}.{decimal}'
+	return float_str
+
+
+@numba.njit(cache=False)
+def str_to_float(value):
+	decimal_pos = value.find('.')
+	if decimal_pos != -1:
+		after_decimal_len = len(value[decimal_pos+1:])
+	else:
+		after_decimal_len = 0
+	value = value.replace('.', '')
+	value = str_to_int(value)
+	value = value / 10**after_decimal_len
+	return value
+
+
+@numba.njit(cache=False)
+def str_to_int(value):
+	if value[0] == '-':
+		negative = True
+		value = value.replace('-','')
+	else:
+		negative = False
+	final_index, result = len(value) - 1, 0
+	for i, v in enumerate(value):
+		result += (ord(v) - 48) * (10 ** (final_index - i))
+	result = -result if negative else result
+	return result
diff --git a/pyreason/scripts/numba_wrapper/numba_types/rule_type.py b/pyreason/scripts/numba_wrapper/numba_types/rule_type.py
index 4eb985b5..a26eb55a 100755
--- a/pyreason/scripts/numba_wrapper/numba_types/rule_type.py
+++ b/pyreason/scripts/numba_wrapper/numba_types/rule_type.py
@@ -23,7 +23,7 @@ def __init__(self):
 rule_type = RuleType()
 
 
-# Type ann_fnerence
+# Type ann_fn
 @typeof_impl.register(Rule)
 def typeof_rule(val, c):
     return rule_type
@@ -32,67 +32,67 @@ def typeof_rule(val, c):
 # Construct object from Numba functions (Doesn't work. We don't need this currently)
 @type_callable(Rule)
 def type_rule(context):
-    def typer(name, target, tc, delta, neigh_criteria, bnd, thresholds, ann_fn, ann_label, weights, edges, immediate_rule):
-        if isinstance(name, types.UnicodeType) and isinstance(target, label.LabelType) and isinstance(tc, (types.NoneType, types.ListType)) and isinstance(delta, types.Integer) and isinstance(neigh_criteria, (types.NoneType, types.ListType)) and isinstance(bnd, interval.IntervalType) and isinstance(thresholds, types.ListType) and isinstance(ann_fn, types.UnicodeType) and isinstance(ann_label, label.LabelType) and isinstance(weights, types.Array) and isinstance(edges, types.Tuple) and isinstance(immediate_rule, types.Boolean):
+    def typer(rule_name, type, target, delta, clauses, bnd, thresholds, ann_fn, weights, edges, static, immediate_rule):
+        if isinstance(rule_name, types.UnicodeType) and isinstance(type, types.UnicodeType) and isinstance(target, label.LabelType) and isinstance(delta, types.Integer) and isinstance(clauses, (types.NoneType, types.ListType)) and isinstance(bnd, interval.IntervalType) and isinstance(thresholds, types.ListType) and isinstance(ann_fn, types.UnicodeType) and isinstance(weights, types.Array) and isinstance(edges, types.Tuple) and isinstance(static, types.Boolean) and isinstance(immediate_rule, types.Boolean):
             return rule_type
     return typer
 
 
-# Define native representation: datamodel
+# Define native representation: data-model
 @register_model(RuleType)
 class RuleModel(models.StructModel):
     def __init__(self, dmm, fe_type):
         members = [
-            ('name', types.string),
+            ('rule_name', types.string),
+            ('type', types.string),
             ('target', label.label_type),
-            ('target_criteria', types.ListType(types.Tuple((label.label_type, interval.interval_type)))),
             ('delta', types.uint16),
-            ('neigh_criteria', types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), label.label_type, interval.interval_type)))),
+            ('clauses', types.ListType(types.Tuple((types.string, label.label_type, types.ListType(types.string), interval.interval_type, types.string)))),
             ('bnd', interval.interval_type),
             ('thresholds', types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), types.float64)))),
             ('ann_fn', types.string),
-            ('ann_label', label.label_type),
             ('weights', types.float64[::1]),
             ('edges', types.Tuple((types.string, types.string, label.label_type))),
+            ('static', types.boolean),
             ('immediate_rule', types.boolean)
             ]
         models.StructModel.__init__(self, dmm, fe_type, members)
 
 
-# Expose datamodel attributes
-make_attribute_wrapper(RuleType, 'name', 'name')
+# Expose data-model attributes
+make_attribute_wrapper(RuleType, 'rule_name', 'rule_name')
+make_attribute_wrapper(RuleType, 'type', 'type')
 make_attribute_wrapper(RuleType, 'target', 'target')
-make_attribute_wrapper(RuleType, 'target_criteria', 'target_criteria')
 make_attribute_wrapper(RuleType, 'delta', 'delta')
-make_attribute_wrapper(RuleType, 'neigh_criteria', 'neigh_criteria')
+make_attribute_wrapper(RuleType, 'clauses', 'clauses')
 make_attribute_wrapper(RuleType, 'bnd', 'bnd')
 make_attribute_wrapper(RuleType, 'thresholds', 'thresholds')
 make_attribute_wrapper(RuleType, 'ann_fn', 'ann_fn')
-make_attribute_wrapper(RuleType, 'ann_label', 'ann_label')
 make_attribute_wrapper(RuleType, 'weights', 'weights')
 make_attribute_wrapper(RuleType, 'edges', 'edges')
+make_attribute_wrapper(RuleType, 'static', 'static')
 make_attribute_wrapper(RuleType, 'immediate_rule', 'immediate_rule')
 
 
 # Implement constructor
-@lower_builtin(Rule, types.string, label.label_type, types.ListType(types.Tuple((label.label_type, interval.interval_type))), types.uint16, types.ListType(types.Tuple((types.string, label.label_type, interval.interval_type))), interval.interval_type, types.ListType(types.ListType(types.Tuple((types.string, types.string, types.float64)))), types.string, label.label_type, types.float64[::1], types.Tuple((types.string, types.string, label.label_type)), types.boolean)
+@lower_builtin(Rule, types.string, types.string, label.label_type, types.uint16, types.ListType(types.Tuple((types.string, label.label_type, types.ListType(types.string), interval.interval_type, types.string))), interval.interval_type, types.ListType(types.ListType(types.Tuple((types.string, types.string, types.float64)))), types.string, types.float64[::1], types.Tuple((types.string, types.string, label.label_type)), types.boolean, types.boolean)
 def impl_rule(context, builder, sig, args):
     typ = sig.return_type
-    name, target, target_criteria, delta, neigh_criteria, bnd, thresholds, ann_fn, ann_label, weights, edges, immediate_rule = args
-    context.nrt.incref(builder, types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), label.label_type, interval.interval_type))), neigh_criteria)
+    rule_name, type, target, delta, clauses, bnd, thresholds, ann_fn, weights, edges, static, immediate_rule = args
+    context.nrt.incref(builder, types.ListType(types.Tuple((types.string, label.label_type, types.ListType(types.string), interval.interval_type, types.string))), clauses)
     context.nrt.incref(builder, types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), types.float64))), thresholds)
     rule = cgutils.create_struct_proxy(typ)(context, builder)
-    rule.name = name
+    rule.rule_name = rule_name
+    rule.type = type
     rule.target = target
-    rule.target_criteria = target_criteria
     rule.delta = delta
-    rule.neigh_criteria = neigh_criteria
+    rule.clauses = clauses
     rule.bnd = bnd
     rule.thresholds = thresholds
     rule.ann_fn = ann_fn
-    rule.ann_label = ann_label
     rule.weights = weights
     rule.edges = edges
+    rule.static = static
     rule.immediate_rule = immediate_rule
     return rule._getvalue()
 
@@ -101,21 +101,21 @@ def impl_rule(context, builder, sig, args):
 @overload_method(RuleType, "get_name")
 def get_name(rule):
     def getter(rule):
-        return rule.name
+        return rule.rule_name
     return getter
 
 
-@overload_method(RuleType, "get_target")
-def get_target(rule):
+@overload_method(RuleType, "get_type")
+def get_type(rule):
     def getter(rule):
-        return rule.target
+        return rule.type
     return getter
 
 
-@overload_method(RuleType, "get_target_criteria")
-def get_target_criteria(rule):
+@overload_method(RuleType, "get_target")
+def get_target(rule):
     def getter(rule):
-        return rule.target_criteria
+        return rule.target
     return getter
 
 
@@ -126,10 +126,10 @@ def getter(rule):
     return getter
 
 
-@overload_method(RuleType, "get_neigh_criteria")
-def get_neigh_criteria(rule):
+@overload_method(RuleType, "get_clauses")
+def get_clauses(rule):
     def getter(rule):
-        return rule.neigh_criteria
+        return rule.clauses
     return getter
 
 
@@ -154,13 +154,6 @@ def impl(rule):
     return impl
 
 
-@overload_method(RuleType, "get_annotation_label")
-def get_label(rule):
-    def impl(rule):
-        return rule.ann_label
-    return impl
-
-
 @overload_method(RuleType, "get_weights")
 def get_weights(rule):
     def impl(rule):
@@ -175,6 +168,13 @@ def impl(rule):
     return impl
 
 
+@overload_method(RuleType, "is_static_rule")
+def is_static_rule(rule):
+    def impl(rule):
+        return rule.static
+    return impl
+
+
 @overload_method(RuleType, "is_immediate_rule")
 def is_immediate_rule(rule):
     def impl(rule):
@@ -185,42 +185,42 @@ def impl(rule):
 # Tell numba how to make native
 @unbox(RuleType)
 def unbox_rule(typ, obj, c):
-    name_obj = c.pyapi.object_getattr_string(obj, "_name")
+    name_obj = c.pyapi.object_getattr_string(obj, "_rule_name")
+    type_obj = c.pyapi.object_getattr_string(obj, "_type")
     target_obj = c.pyapi.object_getattr_string(obj, "_target")
-    tc_obj = c.pyapi.object_getattr_string(obj, "_target_criteria")
     delta_obj = c.pyapi.object_getattr_string(obj, "_delta")
-    neigh_criteria_obj = c.pyapi.object_getattr_string(obj, "_neigh_criteria")
+    clauses_obj = c.pyapi.object_getattr_string(obj, "_clauses")
     bnd_obj = c.pyapi.object_getattr_string(obj, "_bnd")
     thresholds_obj = c.pyapi.object_getattr_string(obj, "_thresholds")
     ann_fn_obj = c.pyapi.object_getattr_string(obj, "_ann_fn")
-    ann_label_obj = c.pyapi.object_getattr_string(obj, "_ann_label")
     weights_obj = c.pyapi.object_getattr_string(obj, "_weights")
     edges_obj = c.pyapi.object_getattr_string(obj, "_edges")
+    static_obj = c.pyapi.object_getattr_string(obj, "_static")
     immediate_rule_obj = c.pyapi.object_getattr_string(obj, "_immediate_rule")
     rule = cgutils.create_struct_proxy(typ)(c.context, c.builder)
-    rule.name = c.unbox(types.string, name_obj).value
+    rule.rule_name = c.unbox(types.string, name_obj).value
+    rule.type = c.unbox(types.string, type_obj).value
     rule.target = c.unbox(label.label_type, target_obj).value
-    rule.target_criteria = c.unbox(types.ListType(types.Tuple((label.label_type, interval.interval_type))), tc_obj).value
     rule.delta = c.unbox(types.uint16, delta_obj).value
-    rule.neigh_criteria = c.unbox(types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), label.label_type, interval.interval_type))), neigh_criteria_obj).value
+    rule.clauses = c.unbox(types.ListType(types.Tuple((types.string, label.label_type, types.ListType(types.string), interval.interval_type, types.string))), clauses_obj).value
     rule.bnd = c.unbox(interval.interval_type, bnd_obj).value
     rule.thresholds = c.unbox(types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), types.float64))), thresholds_obj).value
     rule.ann_fn = c.unbox(types.string, ann_fn_obj).value
-    rule.ann_label = c.unbox(label.label_type, ann_label_obj).value
     rule.weights = c.unbox(types.float64[::1], weights_obj).value
     rule.edges = c.unbox(types.Tuple((types.string, types.string, label.label_type)), edges_obj).value
+    rule.static = c.unbox(types.boolean, static_obj).value
     rule.immediate_rule = c.unbox(types.boolean, immediate_rule_obj).value
     c.pyapi.decref(name_obj)
+    c.pyapi.decref(type_obj)
     c.pyapi.decref(target_obj)
-    c.pyapi.decref(tc_obj)
     c.pyapi.decref(delta_obj)
-    c.pyapi.decref(neigh_criteria_obj)
+    c.pyapi.decref(clauses_obj)
     c.pyapi.decref(bnd_obj)
     c.pyapi.decref(thresholds_obj)
     c.pyapi.decref(ann_fn_obj)
-    c.pyapi.decref(ann_label_obj)
     c.pyapi.decref(weights_obj)
     c.pyapi.decref(edges_obj)
+    c.pyapi.decref(static_obj)
     c.pyapi.decref(immediate_rule_obj)
     is_error = cgutils.is_not_null(c.builder, c.pyapi.err_occurred())
     return NativeValue(rule._getvalue(), is_error=is_error)
@@ -230,30 +230,30 @@ def unbox_rule(typ, obj, c):
 def box_rule(typ, val, c):
     rule = cgutils.create_struct_proxy(typ)(c.context, c.builder, value=val)
     class_obj = c.pyapi.unserialize(c.pyapi.serialize_object(Rule))
-    name_obj = c.box(types.string, rule.name)
+    name_obj = c.box(types.string, rule.rule_name)
+    type_obj = c.box(types.string, rule.type)
     target_obj = c.box(label.label_type, rule.target)
-    tc_obj = c.box(types.ListType(types.Tuple((label.label_type, interval.interval_type))), rule.target_criteria)
     delta_obj = c.box(types.uint16, rule.delta)
-    neigh_criteria_obj = c.box(types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), label.label_type, interval.interval_type))), rule.neigh_criteria)
+    clauses_obj = c.box(types.ListType(types.Tuple((types.string, label.label_type, types.ListType(types.string), interval.interval_type, types.string))), rule.clauses)
     bnd_obj = c.box(interval.interval_type, rule.bnd)
     thresholds_obj = c.box(types.ListType(types.Tuple((types.string, types.UniTuple(types.string, 2), types.float64))), rule.thresholds)
     ann_fn_obj = c.box(types.string, rule.ann_fn)
-    ann_label_obj = c.box(label.label_type, rule.ann_label)
     weights_obj = c.box(types.float64[::1], rule.weights)
     edges_obj = c.box(types.Tuple((types.string, types.string, label.label_type)), rule.edges)
+    static_obj = c.box(types.boolean, rule.static)
     immediate_rule_obj = c.box(types.boolean, rule.immediate_rule)
-    res = c.pyapi.call_function_objargs(class_obj, (name_obj, target_obj, tc_obj, delta_obj, neigh_criteria_obj, bnd_obj, thresholds_obj, ann_fn_obj, ann_label_obj, weights_obj, edges_obj, immediate_rule_obj))
+    res = c.pyapi.call_function_objargs(class_obj, (name_obj, type_obj, target_obj, delta_obj, clauses_obj, bnd_obj, thresholds_obj, ann_fn_obj, weights_obj, edges_obj, static_obj, immediate_rule_obj))
     c.pyapi.decref(name_obj)
+    c.pyapi.decref(type_obj)
     c.pyapi.decref(target_obj)
-    c.pyapi.decref(tc_obj)
     c.pyapi.decref(delta_obj)
-    c.pyapi.decref(neigh_criteria_obj)
+    c.pyapi.decref(clauses_obj)
     c.pyapi.decref(ann_fn_obj)
     c.pyapi.decref(bnd_obj)
     c.pyapi.decref(thresholds_obj)
-    c.pyapi.decref(ann_label_obj)
     c.pyapi.decref(weights_obj)
     c.pyapi.decref(edges_obj)
+    c.pyapi.decref(static_obj)
     c.pyapi.decref(immediate_rule_obj)
     c.pyapi.decref(class_obj)
     return res
diff --git a/pyreason/scripts/program/program.py b/pyreason/scripts/program/program.py
index 2a52c685..0889cbfc 100755
--- a/pyreason/scripts/program/program.py
+++ b/pyreason/scripts/program/program.py
@@ -1,4 +1,5 @@
-from pyreason.scripts.interpretation.interpretation import Interpretation
+from pyreason.scripts.interpretation.interpretation import Interpretation as Interpretation
+from pyreason.scripts.interpretation.interpretation import Interpretation as InterpretationParallel
 
 
 class Program:
@@ -7,18 +8,20 @@ class Program:
 	specific_node_labels = []
 	specific_edge_labels = []
 
-	def __init__(self, graph, facts_node, facts_edge, rules, ipl, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes):
+	def __init__(self, graph, facts_node, facts_edge, rules, ipl, annotation_functions, reverse_graph, atom_trace, save_graph_attributes_to_rule_trace, canonical, inconsistency_check, store_interpretation_changes, parallel_computing):
 		self._graph = graph
 		self._facts_node = facts_node
 		self._facts_edge = facts_edge
 		self._rules = rules
 		self._ipl = ipl
+		self._annotation_functions = annotation_functions
 		self._reverse_graph = reverse_graph
 		self._atom_trace = atom_trace
 		self._save_graph_attributes_to_rule_trace = save_graph_attributes_to_rule_trace
 		self._canonical = canonical
 		self._inconsistency_check = inconsistency_check
 		self._store_interpretation_changes = store_interpretation_changes
+		self._parallel_computing = parallel_computing
 		self.interp = None
 
 	def reason(self, tmax, convergence_threshold, convergence_bound_threshold, verbose=True):
@@ -29,7 +32,11 @@ def reason(self, tmax, convergence_threshold, convergence_bound_threshold, verbo
 		Interpretation.specific_node_labels = self.specific_node_labels
 		Interpretation.specific_edge_labels = self.specific_edge_labels
 
-		self.interp = Interpretation(self._graph, self._ipl, self._reverse_graph, self._atom_trace, self._save_graph_attributes_to_rule_trace, self._canonical, self._inconsistency_check, self._store_interpretation_changes)
+		# Instantiate correct interpretation class based on whether we parallelize the code or not. (We cannot parallelize with cache on)
+		if self._parallel_computing:
+			self.interp = InterpretationParallel(self._graph, self._ipl, self._annotation_functions, self._reverse_graph, self._atom_trace, self._save_graph_attributes_to_rule_trace, self._canonical, self._inconsistency_check, self._store_interpretation_changes)
+		else:
+			self.interp = Interpretation(self._graph, self._ipl, self._annotation_functions, self._reverse_graph, self._atom_trace, self._save_graph_attributes_to_rule_trace, self._canonical, self._inconsistency_check, self._store_interpretation_changes)
 		self.interp.start_fp(self._tmax, self._facts_node, self._facts_edge, self._rules, verbose, convergence_threshold, convergence_bound_threshold)
 
 		return self.interp
@@ -40,4 +47,3 @@ def reason_again(self, tmax, convergence_threshold, convergence_bound_threshold,
 		self.interp.start_fp(self._tmax, facts_node, facts_edge, self._rules, verbose, convergence_threshold, convergence_bound_threshold, again=True)
 
 		return self.interp
-		
diff --git a/pyreason/scripts/rules/rule.py b/pyreason/scripts/rules/rule.py
index 5b417c8b..69de3ca3 100755
--- a/pyreason/scripts/rules/rule.py
+++ b/pyreason/scripts/rules/rule.py
@@ -1,33 +1,33 @@
 class Rule:
 
-    def __init__(self, name, target, target_criteria, delta, neigh_criteria, bnd, thresholds, ann_fn, ann_label, weights, edges, immediate_rule):
-        self._name = name
+    def __init__(self, rule_name, rule_type, target, delta, clauses, bnd, thresholds, ann_fn, weights, edges, static, immediate_rule):
+        self._rule_name = rule_name
+        self._type = rule_type
         self._target = target
-        self._target_criteria = target_criteria
         self._delta = delta
-        self._neigh_criteria = neigh_criteria
+        self._clauses = clauses
         self._bnd = bnd
         self._thresholds = thresholds
         self._ann_fn = ann_fn
-        self._ann_label = ann_label
         self._weights = weights
         self._edges = edges
+        self._static = static
         self._immediate_rule = immediate_rule
 
-    def get_name(self):
-        return self._name
-    
+    def get_rule_name(self):
+        return self._rule_name
+
+    def get_rule_type(self):
+        return self._type
+
     def get_target(self):
         return self._target
 
-    def get_target_criteria(self):
-        return self._target_criteria
-
     def get_delta(self):
         return self._delta
 
     def get_neigh_criteria(self):
-        return self._neigh_criteria
+        return self._clauses
     
     def get_bnd(self):
         return self._bnd
@@ -37,12 +37,12 @@ def get_thresholds(self):
 
     def get_annotation_function(self):
         return self._ann_fn
-
-    def get_ann_label(self):
-        return self._ann_label
     
     def get_edges(self):
         return self._edges
 
+    def is_static(self):
+        return self._static
+
     def is_immediate_rule(self):
         return self._immediate_rule
diff --git a/pyreason/scripts/utils/filter.py b/pyreason/scripts/utils/filter.py
index 84ae7bec..8f7498a0 100755
--- a/pyreason/scripts/utils/filter.py
+++ b/pyreason/scripts/utils/filter.py
@@ -107,9 +107,11 @@ def filter_and_sort_edges(self, interpretation, labels, bound, sort_by='lower',
             dataframe = pd.DataFrame.from_dict(edges[t]).transpose()
             dataframe = dataframe.reset_index()
             if not dataframe.empty:
-                dataframe.columns = ['component', *labels]
+                dataframe.columns = ['source', 'target', *labels]
+                dataframe['component'] = dataframe[['source', 'target']].apply(tuple, axis=1)
+                dataframe.drop(columns=['source', 'target'], inplace=True)
+                dataframe = dataframe[['component', *labels]]
             else:
                 dataframe = pd.DataFrame(columns=['component', *labels])
             dataframes.append(dataframe)
         return dataframes
-
diff --git a/pyreason/scripts/utils/graphml_parser.py b/pyreason/scripts/utils/graphml_parser.py
index d51c3107..967bd4dd 100755
--- a/pyreason/scripts/utils/graphml_parser.py
+++ b/pyreason/scripts/utils/graphml_parser.py
@@ -14,18 +14,21 @@ def __init__(self):
 
     def parse_graph(self, graph_path, reverse):
         self.graph = nx.read_graphml(graph_path)
+        self.graph = nx.DiGraph(self.graph)
         if reverse:
             self.graph = self.graph.reverse()
 
         return self.graph
 
+    def load_graph(self, graph):
+        self.graph = nx.DiGraph(graph)
+        return self.graph
+
     def parse_graph_attributes(self, static_facts):
         facts_node = numba.typed.List.empty_list(fact_node.fact_type)
         facts_edge = numba.typed.List.empty_list(fact_edge.fact_type)
-        specific_node_labels = numba.typed.Dict.empty(key_type=label.label_type,
-                                                      value_type=numba.types.ListType(numba.types.string))
-        specific_edge_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(
-            numba.types.Tuple((numba.types.string, numba.types.string))))
+        specific_node_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(numba.types.string))
+        specific_edge_labels = numba.typed.Dict.empty(key_type=label.label_type, value_type=numba.types.ListType(numba.types.Tuple((numba.types.string, numba.types.string))))
         for n in self.graph.nodes:
             for key, value in self.graph.nodes[n].items():
                 # IF attribute is a float or int and it is less than 1, then make it a bound, else make it a label
@@ -54,8 +57,7 @@ def parse_graph_attributes(self, static_facts):
                 if label.Label(l) not in specific_node_labels.keys():
                     specific_node_labels[label.Label(l)] = numba.typed.List.empty_list(numba.types.string)
                 specific_node_labels[label.Label(l)].append(n)
-                f = fact_node.Fact('graph-attribute-fact', n, label.Label(l), interval.closed(l_bnd, u_bnd), 0, 0,
-                                   static=static_facts)
+                f = fact_node.Fact('graph-attribute-fact', n, label.Label(l), interval.closed(l_bnd, u_bnd), 0, 0, static=static_facts)
                 facts_node.append(f)
         for e in self.graph.edges:
             for key, value in self.graph.edges[e].items():
@@ -83,11 +85,9 @@ def parse_graph_attributes(self, static_facts):
                             pass
 
                 if label.Label(l) not in specific_edge_labels.keys():
-                    specific_edge_labels[label.Label(l)] = numba.typed.List.empty_list(
-                        numba.types.Tuple((numba.types.string, numba.types.string)))
+                    specific_edge_labels[label.Label(l)] = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, numba.types.string)))
                 specific_edge_labels[label.Label(l)].append((e[0], e[1]))
-                f = fact_edge.Fact('graph-attribute-fact', (e[0], e[1]), label.Label(l), interval.closed(l_bnd, u_bnd),
-                                   0, 0, static=static_facts)
+                f = fact_edge.Fact('graph-attribute-fact', (e[0], e[1]), label.Label(l), interval.closed(l_bnd, u_bnd), 0, 0, static=static_facts)
                 facts_edge.append(f)
 
-        return facts_node, facts_edge, specific_node_labels, specific_edge_labels                
+        return facts_node, facts_edge, specific_node_labels, specific_edge_labels
\ No newline at end of file
diff --git a/pyreason/scripts/utils/output.py b/pyreason/scripts/utils/output.py
index 075151d4..e0d12f96 100755
--- a/pyreason/scripts/utils/output.py
+++ b/pyreason/scripts/utils/output.py
@@ -1,89 +1,96 @@
-import pickle
 import csv
 import os
+import pandas as pd
+
 
 class Output:
     def __init__(self, timestamp):
         self.timestamp = timestamp
+        self.rule_trace_node = None
+        self.rule_trace_edge = None
 
-    def save_rule_trace(self, interpretation, folder='./'):
-        # Saves the rule trace in a csv, human readable format
+    def _parse_internal_rule_trace(self, interpretation):
         header_node = ['Time', 'Fixed-Point-Operation', 'Node', 'Label', 'Old Bound', 'New Bound', 'Occurred Due To']
 
         # Nodes rule trace
-        path = os.path.join(folder, f'rule_trace_nodes_{self.timestamp}.csv')
-        with open(path, 'w') as f:
-            data = []
-            max_j = -1
-
-            for i, r in enumerate(interpretation.rule_trace_node):
-                row = [r[0], r[1], r[2], r[3]._value, '-', r[4].to_str(), '-']
-                if interpretation.atom_trace:
-                    qn, qe, old_bnd, name = interpretation.rule_trace_node_atoms[i]
-                    row[4] = old_bnd.to_str()
-                    # Go through all the changes in the rule trace
-                    # len(qn) = len(qe) = num of clauses in rule that was used
-                    row[6] = name
-
-                    # Go through each clause
-                    for j in range(len(qn)):
-                        max_j = max(j, max_j)
-                        if len(qe[j])==0:
-                            # Node clause
-                            row.append(list(qn[j]))
-                        elif len(qn[j])==0:
-                            # Edge clause
-                            row.append(list(qe[j]))
-
-                data.append(row)
-
-            # Add Clause-num to header
-            if interpretation.atom_trace and max_j!=-1:
-                for i in range(1, max_j+2):
-                    header_node.append(f'Clause-{i}')
-
-            writer = csv.writer(f)
-            writer.writerow(header_node)
-            writer.writerows(data)
-        
+        data = []
+        max_j = -1
+
+        for i, r in enumerate(interpretation.rule_trace_node):
+            row = [r[0], r[1], r[2], r[3]._value, '-', r[4].to_str(), '-']
+            if interpretation.atom_trace:
+                qn, qe, old_bnd, name = interpretation.rule_trace_node_atoms[i]
+                row[4] = old_bnd.to_str()
+                # Go through all the changes in the rule trace
+                # len(qn) = len(qe) = num of clauses in rule that was used
+                row[6] = name
+
+                # Go through each clause
+                for j in range(len(qn)):
+                    max_j = max(j, max_j)
+                    if len(qe[j]) == 0:
+                        # Node clause
+                        row.append(list(qn[j]))
+                    elif len(qn[j]) == 0:
+                        # Edge clause
+                        row.append(list(qe[j]))
+
+            data.append(row)
+
+        # Add Clause-num to header
+        if interpretation.atom_trace and max_j != -1:
+            for i in range(1, max_j + 2):
+                header_node.append(f'Clause-{i}')
+
+        # Store the trace in a DataFrame
+        self.rule_trace_node = pd.DataFrame(data, columns=header_node)
+
         header_edge = ['Time', 'Fixed-Point-Operation', 'Edge', 'Label', 'Old Bound', 'New Bound', 'Occurred Due To']
 
         # Edges rule trace
-        path = os.path.join(folder, f'rule_trace_edges_{self.timestamp}.csv')
-        with open(path, 'w') as f:
-            data = []
-            max_j = -1
-
-            for i, r in enumerate(interpretation.rule_trace_edge):
-                row = [r[0], r[1], r[2], r[3]._value, '-', r[4].to_str(), '-']
-                if interpretation.atom_trace:
-                    qn, qe, old_bnd, name = interpretation.rule_trace_edge_atoms[i]
-                    row[4] = old_bnd.to_str()
-                    # Go through all the changes in the rule trace
-                    # len(qn) = num of clauses in rule that was used
-                    row[6] = name
-
-                    # Go through each clause
-                    for j in range(len(qn)):
-                        max_j = max(j, max_j)
-                        if len(qe[j])==0:
-                            # Node clause
-                            row.append(list(qn[j]))
-                        elif len(qn[j])==0:
-                            # Edge clause
-                            row.append(list(qe[j]))
-
-
-                data.append(row)
-
-            # Add Clause-num to header
-            if interpretation.atom_trace and max_j!=-1:
-                for i in range(1, max_j+2):
-                    header_edge.append(f'Clause-{i}')
-
-            writer = csv.writer(f)
-            writer.writerow(header_edge)
-            writer.writerows(data)
-            
+        data = []
+        max_j = -1
+
+        for i, r in enumerate(interpretation.rule_trace_edge):
+            row = [r[0], r[1], r[2], r[3]._value, '-', r[4].to_str(), '-']
+            if interpretation.atom_trace:
+                qn, qe, old_bnd, name = interpretation.rule_trace_edge_atoms[i]
+                row[4] = old_bnd.to_str()
+                # Go through all the changes in the rule trace
+                # len(qn) = num of clauses in rule that was used
+                row[6] = name
+
+                # Go through each clause
+                for j in range(len(qn)):
+                    max_j = max(j, max_j)
+                    if len(qe[j]) == 0:
+                        # Node clause
+                        row.append(list(qn[j]))
+                    elif len(qn[j]) == 0:
+                        # Edge clause
+                        row.append(list(qe[j]))
+
+            data.append(row)
+
+        # Add Clause-num to header
+        if interpretation.atom_trace and max_j != -1:
+            for i in range(1, max_j + 2):
+                header_edge.append(f'Clause-{i}')
+
+        # Store the trace in a DataFrame
+        self.rule_trace_edge = pd.DataFrame(data, columns=header_edge)
+
+    def save_rule_trace(self, interpretation, folder='./'):
+        if self.rule_trace_node is None and self.rule_trace_edge is None:
+            self._parse_internal_rule_trace(interpretation)
+
+        path_nodes = os.path.join(folder, f'rule_trace_nodes_{self.timestamp}.csv')
+        path_edges = os.path.join(folder, f'rule_trace_edges_{self.timestamp}.csv')
+        self.rule_trace_node.to_csv(path_nodes, index=False)
+        self.rule_trace_edge.to_csv(path_edges, index=False)
 
+    def get_rule_trace(self, interpretation):
+        if self.rule_trace_node is None and self.rule_trace_edge is None:
+            self._parse_internal_rule_trace(interpretation)
 
+        return self.rule_trace_node, self.rule_trace_edge
diff --git a/pyreason/scripts/utils/rule_parser.py b/pyreason/scripts/utils/rule_parser.py
new file mode 100644
index 00000000..fc2db4bb
--- /dev/null
+++ b/pyreason/scripts/utils/rule_parser.py
@@ -0,0 +1,220 @@
+import numba
+import numpy as np
+
+import pyreason.scripts.numba_wrapper.numba_types.rule_type as rule
+import pyreason.scripts.numba_wrapper.numba_types.label_type as label
+import pyreason.scripts.numba_wrapper.numba_types.interval_type as interval
+
+
+def parse_rule(rule_text: str, name: str, infer_edges: bool = False, set_static: bool = False, immediate_rule: bool = False) -> rule.Rule:
+    # First remove all spaces from line
+    r = rule_text.replace(' ', '')
+
+    # Separate into head and body
+    head, body = r.split('<-')
+
+    # Extract delta_t of rule if it exists else set it to 0
+    t = ''
+    is_digit = True
+    while is_digit:
+        if body[0].isdigit():
+            t += body[0]
+            body = body[1:]
+        else:
+            is_digit = False
+
+    if t == '':
+        t = 0
+    else:
+        t = int(t)
+
+    # Raw parsing steps
+    # 1. Remove whitespaces
+    # 2. replace ) by )) and ] by ]] so that we can split without damaging the string
+    # 3. Split with ), and then for each element of list, split with ], and add to new list
+    # 4. Then replace ]] with ] and )) with ) in for loop
+    # 5. Add :[1,1] to the end of each element if a bound is not specified
+    # 6. Then split each element with :
+    # 7. Transform bound strings into pr.intervals
+
+    # 2
+    body = body.replace(')', '))')
+    body = body.replace(']', ']]')
+
+    # 3
+    body = body.split('),')
+    split_body = []
+    for b in body:
+        split_body.extend(b.split('],'))
+
+    # 4
+    for i in range(len(split_body)):
+        split_body[i] = split_body[i].replace('))', ')')
+        split_body[i] = split_body[i].replace(']]', ']')
+
+    # 5
+    for i in range(len(split_body)):
+        if split_body[i][-1] != ']':
+            split_body[i] += ':[1,1]'
+
+    # 6
+    body_clauses = []
+    body_bounds = []
+    for b in split_body:
+        clause, bound = b.split(':')
+        body_clauses.append(clause)
+        body_bounds.append(bound)
+
+    # 7
+    for i in range(len(body_bounds)):
+        bound = body_bounds[i]
+        l, u = _str_bound_to_bound(bound)
+        body_bounds[i] = [l, u]
+
+    # Find the target predicate and bounds and annotation function if any.
+    # Possible heads:
+    # pred(x) : [x,y]
+    # pred(x) : f
+    # pred(x)
+
+    # This means there is no bound or annotation function specified
+    if head[-1] == ')':
+        head += ':[1,1]'
+
+    head, head_bound = head.split(':')
+    # Check if we have a bound or annotation function
+    if _is_bound(head_bound):
+        target_bound = list(_str_bound_to_bound(head_bound))
+        target_bound = interval.closed(*target_bound)
+        ann_fn = ''
+    else:
+        target_bound = interval.closed(0, 1)
+        ann_fn = head_bound
+
+    idx = head.find('(')
+    target = head[:idx]
+    target = label.Label(target)
+
+    # Variable(s) in the head of the rule
+    end_idx = head.find(')')
+    head_variables = head[idx + 1:end_idx].split(',')
+
+    # Assign type of rule
+    rule_type = 'node' if len(head_variables) == 1 else 'edge'
+
+    # Get the variables in the body
+    # If there's an operator in the body then discard anything that comes after the operator, but keep the variables
+    body_predicates = []
+    body_variables = []
+    for clause in body_clauses:
+        start_idx = clause.find('(')
+        end_idx = clause.find(')')
+        body_predicates.append(clause[:start_idx])
+
+        # Add body variables depending on whether there's an operator or not
+        variables = clause[start_idx+1:end_idx].split(',')
+        start_idx = clause.find('(', start_idx+1)
+        end_idx = clause.find(')', end_idx+1)
+        if start_idx != -1 and end_idx != -1:
+            variables += clause[start_idx+1:end_idx].split(',')
+        body_variables.append(variables)
+
+    # Change infer edge parameter if it's a node rule
+    if rule_type == 'node':
+        infer_edges = False
+
+    # Replace the variables in the body with source/target if they match the variables in the head
+    # If infer_edges is true, then we consider all rules to be node rules, we infer the 2nd variable of the target predicate from the rule body
+    # Else we consider the rule to be an edge rule and replace variables with source/target
+    # Node rules with possibility of adding edges
+    if infer_edges or len(head_variables) == 1:
+        head_source_variable = head_variables[0]
+        for i in range(len(body_variables)):
+            for j in range(len(body_variables[i])):
+                if body_variables[i][j] == head_source_variable:
+                    body_variables[i][j] = '__target'
+    # Edge rule, no edges to be added
+    elif len(head_variables) == 2:
+        for i in range(len(body_variables)):
+            for j in range(len(body_variables[i])):
+                if body_variables[i][j] == head_variables[0]:
+                    body_variables[i][j] = '__source'
+                elif body_variables[i][j] == head_variables[1]:
+                    body_variables[i][j] = '__target'
+
+    # Start setting up clauses
+    # clauses = [c1, c2, c3, c4]
+    # thresholds = [t1, t2, t3, t4]
+
+    # Array of thresholds to keep track of for each neighbor criterion. Form [(comparison, (number/percent, total/available), thresh)]
+    thresholds = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, numba.types.UniTuple(numba.types.string, 2), numba.types.float64)))
+
+    # Array to store clauses for nodes: node/edge, [subset]/[subset1, subset2], label, interval, operator
+    clauses = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, label.label_type, numba.types.ListType(numba.types.string), interval.interval_type, numba.types.string)))
+
+    # Loop though clauses
+    for body_clause, predicate, variables, bounds in zip(body_clauses, body_predicates, body_variables, body_bounds):
+        # Neigh criteria
+        clause_type = 'node' if len(variables) == 1 else 'edge'
+        op = _get_operator_from_clause(body_clause)
+        if op:
+            clause_type = 'comparison'
+
+        subset = numba.typed.List(variables)
+        l = label.Label(predicate)
+        bnd = interval.closed(bounds[0], bounds[1])
+        clauses.append((clause_type, l, subset, bnd, op))
+
+        # Threshold.
+        quantifier = 'greater_equal'
+        quantifier_type = ('number', 'total')
+        thresh = 1
+        thresholds.append((quantifier, quantifier_type, thresh))
+
+    # Assert that there are two variables in the head of the rule if we infer edges
+    # Add edges between head variables if necessary
+    if infer_edges:
+        var = '__target' if head_variables[0] == head_variables[1] else head_variables[1]
+        edges = ('__target', var, target)
+    else:
+        edges = ('', '', label.Label(''))
+
+    weights = np.ones(len(body_predicates), dtype=np.float64)
+    weights = np.append(weights, 0)
+
+    r = rule.Rule(name, rule_type, target, numba.types.uint16(t), clauses, target_bound, thresholds, ann_fn, weights, edges, set_static, immediate_rule)
+    return r
+
+
+def _str_bound_to_bound(str_bound):
+    str_bound = str_bound.replace('[', '')
+    str_bound = str_bound.replace(']', '')
+    l, u = str_bound.split(',')
+    return float(l), float(u)
+
+
+def _is_bound(str_bound):
+    str_bound = str_bound.replace('[', '')
+    str_bound = str_bound.replace(']', '')
+    try:
+        l, u = str_bound.split(',')
+        l = l.replace('.', '')
+        u = u.replace('.', '')
+        if l.isdigit() and u.isdigit():
+            result = True
+        else:
+            result = False
+    except:
+        result = False
+
+    return result
+
+
+def _get_operator_from_clause(clause):
+    operators = ['<=', '>=', '<', '>', '==', '!=']
+    for op in operators:
+        if op in clause:
+            return op
+
+    # No operator found in clause
+    return ''
diff --git a/pyreason/scripts/utils/yaml_parser.py b/pyreason/scripts/utils/yaml_parser.py
index b3d82c36..d93df9fc 100755
--- a/pyreason/scripts/utils/yaml_parser.py
+++ b/pyreason/scripts/utils/yaml_parser.py
@@ -21,12 +21,6 @@ def parse_rules(path):
         if values['target'] is not None:
             target = label.Label(values['target'])
 
-        # Set rule target criteria
-        target_criteria = numba.typed.List.empty_list(numba.types.Tuple((label.label_type, interval.interval_type)))
-        if values['target_criteria'] is not None:
-            for tc in values['target_criteria']:
-                target_criteria.append((label.Label(tc[0]), interval.closed(tc[1], tc[2])))
-
         # Set delta t
         delta_t = numba.types.uint16(values['delta_t'])
 
@@ -37,7 +31,7 @@ def parse_rules(path):
         thresholds = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, numba.types.UniTuple(numba.types.string, 2), numba.types.float64)))
 
         # Array to store clauses for nodes: node/edge, [subset]/[subset1, subset2], label, interval
-        neigh_criteria = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, numba.types.UniTuple(numba.types.string, 2), label.label_type, interval.interval_type)))
+        neigh_criteria = numba.typed.List.empty_list(numba.types.Tuple((numba.types.string, label.label_type, numba.types.UniTuple(numba.types.string, 2), interval.interval_type)))
         if values['neigh_criteria'] is not None:
             # Loop through clauses
             for clause in values['neigh_criteria']:
@@ -47,7 +41,7 @@ def parse_rules(path):
                 subset = (clause[1][0], clause[1][0]) if clause_type=='node' else (clause[1][0], clause[1][1])
                 l = label.Label(clause[2])
                 bnd = interval.closed(clause[3][0], clause[3][1])
-                neigh_criteria.append((clause_type, subset, l, bnd))
+                neigh_criteria.append((clause_type, l, subset, bnd))
 
                 # Append threshold corresponding to clause if specified in rule, else use default of greater equal to 1
                 if len(clause)>4:
@@ -98,7 +92,8 @@ def parse_rules(path):
         if 'immediate' in values and values['immediate'] is not None:
             immediate_rule = True
 
-        r = rule.Rule(rule_name, target, target_criteria, delta_t, neigh_criteria, bnd, thresholds, ann_fn, ann_label, weights, edges, immediate_rule)
+        # Dummy rule type. this file is deprecated
+        r = rule.Rule(rule_name, 'node', target, delta_t, neigh_criteria, bnd, thresholds, ann_fn, ann_label, weights, edges, immediate_rule)
 
         # Insert to beginning of list if flag for immediate rule is true
         if immediate_rule:
diff --git a/requirements.txt b/requirements.txt
index 388b4645..61b7a78a 100755
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,4 +3,5 @@ pyyaml
 pandas
 numba
 numpy
-memory_profiler
\ No newline at end of file
+memory_profiler
+pytest
\ No newline at end of file
diff --git a/setup.py b/setup.py
index ca1be29b..3bbe05b6 100644
--- a/setup.py
+++ b/setup.py
@@ -2,38 +2,39 @@
 
 # Read the contents of README file
 from pathlib import Path
+
 this_directory = Path(__file__).parent
 long_description = (this_directory / "README.md").read_text()
 
-
 setup(
-    name = 'pyreason',
-    version = '1.6.4',
-    author = 'Dyuman Aditya',
-    author_email = 'dyuman.aditya@gmail.com',
-    description = 'An explainable inference software supporting annotated, real valued, graph based and temporal logic',
-    long_description = long_description,
-    long_description_content_type = 'text/markdown',
-    url = 'https://github.com/lab-v2/pyreason',
-    license = 'BSD 3-clause',
-    project_urls = {
+    name='pyreason',
+    version='2.0.0',
+    author='Dyuman Aditya',
+    author_email='dyuman.aditya@gmail.com',
+    description='An explainable inference software supporting annotated, real valued, graph based and temporal logic',
+    long_description=long_description,
+    long_description_content_type='text/markdown',
+    url='https://github.com/lab-v2/pyreason',
+    license='BSD 3-clause',
+    project_urls={
         'Bug Tracker': 'https://github.com/lab-v2/pyreason/issues',
         'Repository': 'https://github.com/lab-v2/pyreason'
     },
-    classifiers = [
+    classifiers=[
         "Programming Language :: Python :: 3",
         "License :: OSI Approved :: BSD License",
         "Operating System :: OS Independent"
     ],
-    python_requires = '>3.6',
-    install_requires = [
+    python_requires='>3.6',
+    install_requires=[
         'networkx',
         'pyyaml',
         'pandas',
         'numba',
         'numpy',
-        'memory_profiler'
+        'memory_profiler',
+        'pytest'
     ],
-    packages = find_packages(),
+    packages=find_packages(),
     include_package_data=True
 )


Test patch:
diff --git a/.github/workflows/python-package-version-test.yml b/.github/workflows/python-package-version-test.yml
index 3ba870c5..3bc43864 100644
--- a/.github/workflows/python-package-version-test.yml
+++ b/.github/workflows/python-package-version-test.yml
@@ -35,7 +35,6 @@ jobs:
         flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
-    # Add tests later when library is made
-    # - name: Test with pytest
-    #  run: |
-    #    pytest
+    - name: Test with pytest
+      run: |
+        pytest
diff --git a/tests/friends_graph.graphml b/tests/friends_graph.graphml
new file mode 100644
index 00000000..931cd697
--- /dev/null
+++ b/tests/friends_graph.graphml
@@ -0,0 +1,33 @@
+<?xml version='1.0' encoding='utf-8'?>
+<graphml xmlns="http://graphml.graphdrawing.org/xmlns" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://graphml.graphdrawing.org/xmlns http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd">
+  <key id="owns" for="edge" attr.name="owns" attr.type="long" />
+  <key id="Friends" for="edge" attr.name="Friends" attr.type="long" />
+  <graph edgedefault="directed">
+    <node id="John" />
+    <node id="Mary" />
+    <node id="Justin" />
+    <node id="Dog" />
+    <node id="Cat" />
+    <edge source="John" target="Mary">
+      <data key="Friends">1</data>
+    </edge>
+    <edge source="John" target="Justin">
+      <data key="Friends">1</data>
+    </edge>
+    <edge source="John" target="Dog">
+      <data key="owns">1</data>
+    </edge>
+    <edge source="Mary" target="Cat">
+      <data key="owns">1</data>
+    </edge>
+    <edge source="Justin" target="Mary">
+      <data key="Friends">1</data>
+    </edge>
+    <edge source="Justin" target="Cat">
+      <data key="owns">1</data>
+    </edge>
+    <edge source="Justin" target="Dog">
+      <data key="owns">1</data>
+    </edge>
+  </graph>
+</graphml>
diff --git a/tests/test_hello_world.py b/tests/test_hello_world.py
new file mode 100644
index 00000000..b6d6f9c0
--- /dev/null
+++ b/tests/test_hello_world.py
@@ -0,0 +1,38 @@
+# Test if the simple hello world program works
+import pyreason as pr
+
+
+def test_hello_world():
+    # Modify the paths based on where you've stored the files we made above
+    graph_path = './tests/friends_graph.graphml'
+
+    # Modify pyreason settings to make verbose and to save the rule trace to a file
+    pr.settings.verbose = True     # Print info to screen
+
+    # Load all the files into pyreason
+    pr.load_graphml(graph_path)
+    pr.add_rule('popular(x) <-1 popular(y), Friends(x,y), owns(y,z), owns(x,z)', 'popular_rule')
+    pr.add_fact(pr.Fact('popular-fact', 'Mary', 'popular', [1, 1], 0, 2))
+
+    # Run the program for two timesteps to see the diffusion take place
+    interpretation = pr.reason(timesteps=2)
+
+    # Display the changes in the interpretation for each timestep
+    dataframes = pr.filter_and_sort_nodes(interpretation, ['popular'])
+    for t, df in enumerate(dataframes):
+        print(f'TIMESTEP - {t}')
+        print(df)
+        print()
+
+    assert len(dataframes[0]) == 1, 'At t=0 there should be one popular person'
+    assert len(dataframes[1]) == 2, 'At t=0 there should be two popular people'
+    assert len(dataframes[2]) == 3, 'At t=0 there should be three popular people'
+
+    assert dataframes[0].iloc[0].component == 'Mary' and dataframes[0].iloc[0].popular == [1, 1], 'Mary should have popular bounds [1,1] for t=0 timesteps'
+    assert dataframes[1].iloc[0].component == 'Mary' and dataframes[1].iloc[0].popular == [1, 1], 'Mary should have popular bounds [1,1] for t=1 timesteps'
+    assert dataframes[2].iloc[0].component == 'Mary' and dataframes[2].iloc[0].popular == [1, 1], 'Mary should have popular bounds [1,1] for t=2 timesteps'
+
+    assert dataframes[1].iloc[1].component == 'Justin' and dataframes[1].iloc[1].popular == [1, 1], 'Justin should have popular bounds [1,1] for t=1 timesteps'
+    assert dataframes[2].iloc[2].component == 'Justin' and dataframes[2].iloc[2].popular == [1, 1], 'Justin should have popular bounds [1,1] for t=2 timesteps'
+
+    assert dataframes[2].iloc[1].component == 'John' and dataframes[2].iloc[1].popular == [1, 1], 'John should have popular bounds [1,1] for t=2 timesteps'