diff --git a/examples/api/sudoku_zero.py b/examples/api/sudoku_zero.py
index 341b090..330bb96 100644
--- a/examples/api/sudoku_zero.py
+++ b/examples/api/sudoku_zero.py
@@ -1,5 +1,4 @@
-# import clingo
-# from dumbo_asp.primitives import SymbolicAtom, SymbolicProgram, Model, GroundAtom, Predicate
+import time
 
 from dumbo_asp.queries import open_graph_in_xasp_navigator
 
@@ -7,16 +6,15 @@
 
 program = program_from_files(["examples/sudoku/instance4x4.lp","examples/sudoku/encoding4x4.lp"])
 
-atoms_in_answer_set = Model.of_program(program).filter(when=lambda atom: atom.predicate_name == 'assign')
+atoms_in_answer_set = Model.of_program(program)  # .filter(when=lambda atom: atom.predicate_name == 'assign')
 answer_set = tuple([tuple([a,True]) for a in atoms_in_answer_set])
 
-#print(model.as_facts)
-#assert False
+# print(atoms_in_answer_set)
+# assert False
 query_atom = GroundAtom.parse("assign((1,2),2)")
 
 herbrand_base = get_herbrand_zero(program, atoms_in_answer_set)
 
-
 # Add atoms from herbrand to answer set
 answer_set = extend_answer_set(answer_set, herbrand_base)
 
@@ -29,30 +27,38 @@
     program, answer_set, query_atom, MUS_PREDICATE
 )
 
+print(time.time())
 selectors = get_selectors(
     extended_program, MUS_PREDICATE, all_selectors
 )
+print(time.time())
 selectors_prg = SymbolicProgram.parse(Model.of_atoms(selectors).as_facts)
+print(time.time())
+# selectors_prg = SymbolicProgram.parse(Model.of_atoms(all_selectors).as_facts)  ## MALVI: REMOVE ME!
 
 # Remove __false__ from rules
 extended_program = remove_false(extended_program)
+print(time.time())
 expanded_prg_with_selectors = SymbolicProgram.of(*extended_program, *selectors_prg)
+print(time.time())
 
 # Ground program manually
-expanded_prg = expanded_prg_with_selectors.expand_global_and_local_variables()
+herbrand_base = Model.of_atoms(*herbrand_base, *all_selectors)  # MALVI: add selectors to HB (otherwise we cannot expand rules)
+print(time.time())
+expanded_prg = expanded_prg_with_selectors.expand_global_and_local_variables(herbrand_base=herbrand_base)
+print(time.time())
 
 # Serialize
-serialization_program  = get_serialization_program(expanded_prg, query_atom)
-
-# Reify
-reified_program = get_reified_program(serialization_program)
+serialization = get_serialization(expanded_prg, query_atom)
+print('ser',time.time())
 
 # Get derivation
-derivation_sequence = get_derivation_sequence_program(reified_program)
-# WARNING! Here I'm assuming that atoms are ordered according to the derivation in the solver. If it is not, we need a propagator or something different
+derivation_sequence = get_derivation_sequence(serialization)
+print(time.time())
 
 # Get graph
-graph = get_graph(derivation_sequence, serialization_program)
+graph = get_graph(derivation_sequence, serialization)
+print(time.time())
 
 # Show graph
-open_graph_in_xasp_navigator(graph)
\ No newline at end of file
+open_graph_in_xasp_navigator(graph)
diff --git a/ucorexplain/__init__.py b/ucorexplain/__init__.py
index 7c39afe..1a066b4 100644
--- a/ucorexplain/__init__.py
+++ b/ucorexplain/__init__.py
@@ -19,13 +19,13 @@
 
 
 def unpack_answer_set_element(element: AnswerSetElement) -> tuple[GroundAtom, bool]:
-    if type(element) != GroundAtom:
+    if type(element) is not GroundAtom:
         return element
     return element, True
 
 
 def answer_set_element_to_string(element: AnswerSetElement, *, flip: bool = False) -> str:
-    if type(element) == GroundAtom:
+    if type(element) is GroundAtom:
         element = (element, True)
     return f"{'' if element[1] != flip else 'not '}{element[0]}"
 
@@ -38,12 +38,21 @@ def key(element):
     return tuple(sorted(answer_set, key=key))
 
 
+def path(file: str) -> str:
+    import pathlib
+    directory = pathlib.Path(__file__).parent.resolve() / ".."
+    return str(directory / file)
+
+
+def file_to_str(file: str) -> str:
+    with open(path(file)) as f:
+        return f.read()
+
+
 def program_from_files(files) -> SymbolicProgram:
-    prg_str = ""
-    for f in files:
-        with open(f) as f:
-            prg_str += "".join(f.readlines())
-    return SymbolicProgram.parse(prg_str)
+    return SymbolicProgram.parse(
+        '\n'.join(file_to_str(file) for file in files)
+    )
 
 
 def answer_set_to_constraints(
@@ -54,7 +63,7 @@ def answer_set_to_constraints(
     """
     Produces the sequence of selecting constraints.
     """
-    if type(query_atom) == GroundAtom:
+    if type(query_atom) is GroundAtom:
         query_atom = (query_atom,)
     query_atoms = set(query_atom)
     query_literals = []
@@ -106,7 +115,7 @@ def build_extended_program_and_possible_selectors(
             + "}."
         ),
         # MALVI: WE STILL NEED SOMETHING FROM to_zero_simplification_version()
-        SymbolicRule.parse(' | '.join(str(atom) for atom in atoms_in_the_answer_set) + f" :- {false_predicate}."),
+        SymbolicRule.parse('{' + '; '.join(str(atom) for atom in atoms_in_the_answer_set) + f"}} :- {false_predicate}."),
         SymbolicRule.parse(f"{{{false_predicate}}}."),
         SymbolicRule.parse(f":- {false_predicate}."),
     )
@@ -187,7 +196,7 @@ def get_selectors(extended_program_, mus_predicate, all_selectors):
     control, selector_to_literal, literal_to_selector = build_control_and_maps(
         extended_program_, mus_predicate
     )
-    selectors = all_selectors
+    selectors = list(all_selectors)  # MALVI: make a defensive copy
     trim_selectors(
         control=control,
         selectors=selectors,
@@ -261,6 +270,7 @@ def get_herbrand_choices(program, known_atoms):
     print_with_title("HERBRAND BASE", herbrand)
     return herbrand
 
+
 def get_herbrand_zero(program, known_atoms, compact=True):
     """
     Adds choices for the known atoms and computes the herbrand base
@@ -274,12 +284,14 @@ def get_herbrand_zero(program, known_atoms, compact=True):
     print_with_title("HERBRAND BASE", herbrand)
     return herbrand
 
+
 def remove_false(extended_program):
     """
     Removes rules using __false__ added in the build_extended_program_and_possible_selectors method
     This is just the last three rules that are added.
     """
-    return SymbolicProgram.of([rule for rule in extended_program][:-3])  # MALVI: DISCARD rules with __false__
+    return SymbolicProgram.of(*([rule for rule in extended_program][:-3]))  # MALVI: DISCARD rules with __false__
+
 
 def get_reified_program(prg):
     reify_program = program_from_files(["ucorexplain/encodings/mario_reify.lp"])
@@ -287,48 +299,71 @@ def get_reified_program(prg):
     reified_program = Model.of_program(full_prg).as_facts
     return reified_program
 
-def get_serialization_program(expanded_prg, query_atom):
-    expanded_prg_facts = Model.of_atoms(expanded_prg.serialize(base64_encode=False)).as_facts
-    expanded_prg_program = SymbolicProgram.parse(expanded_prg_facts)
-    query_program = SymbolicProgram.parse(f'query({clingo.String(str(query_atom))}).')
-    serialization_program = SymbolicProgram.of(*expanded_prg_program, *query_program)
-    print_with_title("SERIALIZATION PROGRAM", serialization_program)
-
-    return serialization_program
 
-def get_derivation_sequence_program(reified_program):
-    derivation_sequence = []
-    ctl = clingo.Control(["1"])
-    ctl.load("ucorexplain/encodings/mario_meta.lp")
-    ctl.add("base",[],reified_program)
+def get_serialization(expanded_prg, query_atom) -> Model:
+    res = Model.of_atoms(
+        *expanded_prg.serialize(base64_encode=False),
+        GroundAtom.parse(f'query({clingo.String(str(query_atom))})')
+    )
+    print_with_title("SERIALIZATION", res.as_facts)
+    return res
+
+
+def get_derivation_sequence(serialization: Model) -> Model:
+    seen = set()
+    res = []
+    terminate = []
+
+    def collect(model):
+        for at in model.symbols(shown=True):
+            atom = GroundAtom(at)
+            if atom.predicate_name == "assign":  # :- atom(Atom), #count{Value, Reason : assign(Atom, Value, Reason)} > 1.
+                key = (atom.predicate_name, atom.arguments[0])
+            elif atom.predicate_name == "cannot_support":  # :- head(Rule, HeadAtom), #count{Atom : cannot_support(Rule, HeadAtom, Atom)} > 1.
+                key = (atom.predicate_name, atom.arguments[0], atom.arguments[1])
+            elif atom.predicate_name == "constraint":  # :- rule(Rule), #count{Bound : constraint(Rule, Bound)} > 1.
+                key = (atom.predicate_name, atom.arguments[0])
+            elif atom.predicate_name == "terminate":
+                terminate.append(True)
+                continue
+            else:
+                assert False
+            if key not in seen:
+                seen.add(key)
+                res.append(atom)
+
+    while not terminate:
+        previous_len = len(res)
+        ctl = clingo.Control(["1", "--solve-limit=1"])
+        ctl.add("base", [], file_to_str("ucorexplain/encodings/mario_meta.lp"))
+        ctl.add("base", [], serialization.as_facts)
+        ctl.add("base", [], '\n'.join(f"{atom}." for atom in res))
+        ctl.ground([("base", [])])
+        ctl.solve(on_model=collect)
+        assert len(res) > previous_len
+    res = Model.of_elements(*res)
+
+    print_with_title("DERIVATION SEQUENCE", res.filter(when=lambda atom: atom.predicate_name == "assign").as_facts)
+    return res
+
+
+def get_graph(derivation_sequence: Model, serialization: Model) -> Model:
+    ctl = clingo.Control(["1", "--solve-limit=1"])
+    ctl.add("base", [], file_to_str("ucorexplain/encodings/mario_graph.lp"))
+    ctl.add("base", [], serialization.as_facts)
+    ctl.add("base", [], derivation_sequence.as_facts)
     ctl.ground([("base", [])])
-    def m(atoms):
-        for atom in atoms.symbols(shown=True):
-            at = GroundAtom.parse(str(atom))
-            derivation_sequence.append(f"{at.predicate_name}({','.join(str(arg) for arg in at.arguments)},{len(derivation_sequence) + 1})")
-    ctl.solve(on_model=m)
-
-    derivation_sequence_prg = SymbolicProgram.parse(Model.of_atoms(derivation_sequence).as_facts)
-    # WARNING! Here I'm assuming that atoms are ordered according to the derivation in the solver. If it is not, we need a propagator or something different
-    print_with_title("DERIVATION SEQUENCE", derivation_sequence_prg)
-    
-    return derivation_sequence_prg
-
-def get_graph(derivation_sequence_prg, serialization_program):
-    graph_program = program_from_files(["ucorexplain/encodings/mario_graph.lp"])
-    compute_graph_prg = SymbolicProgram.of(*graph_program, *derivation_sequence_prg, *serialization_program)
-
-    graph = Model.of_program(compute_graph_prg)
-    graph = graph.filter(when=lambda atom: atom.predicate_name in ["node", "link'"])
-    graph = graph.rename(Predicate.parse("link'"), Predicate.parse("link"))
-    print_with_title("GRAPH", derivation_sequence_prg)
 
+    graph = Model.of_control(ctl)
+    print_with_title("GRAPH", graph.as_facts)
     return graph
 
+
 def print_with_title(title, value):
+    return
     console.print(f"[bold red]{title}:[/bold red]")
 
-    if type(value) == list:
+    if type(value) is list:
         for e in value:
             console.print(f"{e}")
     else:
@@ -369,7 +404,7 @@ def get_mus_program_and_selectors(
     query_atom: GroundAtom | tuple[GroundAtom, ...],
 ) -> Optional[SymbolicProgram]:
     """
-    Builds MUS program extended with  with externals and the needed selectors
+    Builds MUS program extended with externals and the needed selectors
     """
     atoms_in_answer_set = set(element[0] for element in answer_set)
     
diff --git a/ucorexplain/encodings/mario_graph.lp b/ucorexplain/encodings/mario_graph.lp
index ccbd480..bffc83a 100644
--- a/ucorexplain/encodings/mario_graph.lp
+++ b/ucorexplain/encodings/mario_graph.lp
@@ -1,55 +1,52 @@
 link(Atom, BodyAtom) :-
-  assign(Atom, _, (support, Rule), _);
+  assign(Atom, _, (support, Rule));
   pos_body(Rule, BodyAtom).
 link(Atom, BodyAtom) :-
-  assign(Atom, _, (support, Rule), _);
+  assign(Atom, _, (support, Rule));
   neg_body(Rule, BodyAtom).
-link(Atom, HeadAtom) :-
-  assign(Atom, _, (support, Rule), _);
-  head(Rule, HeadAtom), assign(HeadAtom, false, _, _).
+link(Atom, FalseHeadAtom) :-
+  assign(Atom, _, (support, Rule));
+  head(Rule, FalseHeadAtom), assign(FalseHeadAtom, false, _).
 
 link(Atom, BodyAtom) :-
-  assign(Atom, _, (choice, Rule), _);
+  assign(Atom, _, (head_upper_bound, Rule));
   pos_body(Rule, BodyAtom).
 link(Atom, BodyAtom) :-
-  assign(Atom, _, (choice, Rule), _);
+  assign(Atom, _, (head_upper_bound, Rule));
   neg_body(Rule, BodyAtom).
-link(Atom, HeadAtom) :-
-  assign(Atom, _, (choice, Rule), Index);
-  head(Rule, HeadAtom), assign(HeadAtom, true, _, Index'), Index' < Index.
-
-{link(Atom, Atom'): false_body(Rule, Atom', Index'), Index' < Index} = 1 :-
-  assign(Atom, _, (lack_of_support,), Index);
-  head(Rule, Atom).
-  % disj rule?
-
-%link(HeadAtom, Atom) :-
-%  assign(Atom, _, (last_support, Rule, Atom), Index);
-%  head(Rule, HeadAtom), assign(HeadAtom, _, _, Index'), Index' < Index.
-% handle disjunctive rules?
-link(BodyAtom, Atom) :-
-  assign(BodyAtom, _, (last_support, Rule, Atom), _).
-{link(BodyAtom, Atom'): false_body(Rule', Atom', Index'), Index' < Index} = 1 :-
-  assign(BodyAtom, _, (last_support, Rule, Atom), Index);
-  head(Rule', Atom), false_body(Rule, Atom, _).
-
-link(Atom, HeadAtom) :-
-  assign(Atom, _, (constraint, Rule), Index);
-  head(Rule, HeadAtom), assign(HeadAtom, _, _, Index'), Index' < Index.
+link(Atom, TrueHeadAtom) :-
+  assign(Atom, _, (head_upper_bound, Rule));
+  head(Rule, TrueHeadAtom), assign(TrueHeadAtom, true, _).
+
+link(Atom, BecauseOfAtom) :-
+  assign(Atom, _, (lack_of_support,));
+  head(Rule, Atom), cannot_support(Rule, Atom, BecauseOfAtom).
+
+link(Atom, AtomToSupport) :-
+  assign(Atom, _, (last_support, Rule, AtomToSupport)).
+link(Atom, BecauseOfAtom) :-
+  assign(Atom, _, (last_support, Rule, AtomToSupport));
+  cannot_support(Rule, AtomToSupport, BecauseOfAtom).
+
+link(Atom, TrueHeadAtom) :-
+  assign(Atom, _, (constraint, Rule, upper_bound));
+  head(Rule, TrueHeadAtom), assign(TrueHeadAtom, true, _).
+link(Atom, FalseHeadAtom) :-
+  assign(Atom, _, (constraint, Rule, lower_bound));
+  head(Rule, FalseHeadAtom), assign(FalseHeadAtom, false, _).
 link(Atom, BodyAtom) :-
-  assign(Atom, _, (constraint, Rule), _);
+  assign(Atom, _, (constraint, Rule, _));
   pos_body(Rule, BodyAtom), BodyAtom != Atom.
 link(Atom, BodyAtom) :-
-  assign(Atom, _, (constraint, Rule), _);
+  assign(Atom, _, (constraint, Rule, _));
   neg_body(Rule, BodyAtom), BodyAtom != Atom.
 
+
 reach(Atom) :- query(Atom).
 reach(Atom') :- reach(Atom), link(Atom, Atom'), not hide(Atom').
 
 hide(Atom) :- head(Rule, Atom); not pos_body(Rule,_); not neg_body(Rule,_).
 
-%#show.
-%#show node(X,V,R) : assign(X,V,R,_), reach(X).
-%#show link(X,Y) : link(X,Y), reach(X), reach(Y).
-node(X,V,R) :- assign(X,V,R,_), reach(X).
-link'(X,Y) :- link(X,Y), reach(X), reach(Y).
\ No newline at end of file
+#show.
+#show node(X,V,R) : assign(X,V,R), reach(X).
+#show link(X,Y) : link(X,Y), reach(X), reach(Y).
\ No newline at end of file
diff --git a/ucorexplain/encodings/mario_meta.lp b/ucorexplain/encodings/mario_meta.lp
index 986d39b..278131c 100644
--- a/ucorexplain/encodings/mario_meta.lp
+++ b/ucorexplain/encodings/mario_meta.lp
@@ -1,54 +1,76 @@
-:- atom(Atom), #count{Value, Reason : assign(Atom, Value, Reason)} > 1.
-:- query(Atom), not assign(Atom, _, _).
+head_bounds(Rule, LowerBound, UpperBound) :-
+  rule(Rule), choice(Rule, LowerBound, UpperBound).
+head_bounds(Rule, 1, Size) :-
+  rule(Rule), not choice(Rule,_,_);
+  Size = #count{HeadAtom : head(Rule, HeadAtom)}.
 
-{assign(HeadAtom, true, (support, Rule))} :-
+atom(Atom) :- head(Rule, Atom).
+atom(Atom) :- pos_body(Rule, Atom).
+atom(Atom) :- neg_body(Rule, Atom).
+
+%:- atom(Atom), #count{Value, Reason : assign(Atom, Value, Reason)} > 1.
+terminate :- assign'(Atom, _, _) : query(Atom).
+
+assign'(HeadAtom, true, (support, Rule)) :-
   rule(Rule), head_bounds(Rule, LowerBound, UpperBound);
   head(Rule, HeadAtom), #sum{1, Atom : head(Rule, Atom); -1, Atom : head(Rule, Atom), assign(Atom, false, _)} = LowerBound;
   assign(BodyAtom, true, _) : pos_body(Rule, BodyAtom);
   assign(BodyAtom, false, _) : neg_body(Rule, BodyAtom).
 
-{assign(HeadAtom, false, (choice, Rule))} :-
+assign'(HeadAtom, false, (head_upper_bound, Rule)) :-
   rule(Rule), head_bounds(Rule, LowerBound, UpperBound);
   head(Rule, HeadAtom), #count{Atom : head(Rule, Atom), assign(Atom, true, _), Atom != HeadAtom} = UpperBound;
   assign(BodyAtom, true, _) : pos_body(Rule, BodyAtom);
   assign(BodyAtom, false, _) : neg_body(Rule, BodyAtom).
 
-false_body(Rule, Atom) :-
-  rule(Rule), pos_body(Rule, Atom), assign(Atom, false, _).
-false_body(Rule, Atom) :-
-  rule(Rule), neg_body(Rule, Atom), assign(Atom, true, _).
+%:- head(Rule, HeadAtom), #count{Atom : cannot_support(Rule, HeadAtom, Atom)} > 1.
+cannot_support(Rule, HeadAtom, OtherHeadAtom) :-
+  rule(Rule), head(Rule, HeadAtom), not choice(Rule,_,_);
+  head(Rule, OtherHeadAtom), OtherHeadAtom != HeadAtom, assign(OtherHeadAtom, true, _).
+cannot_support(Rule, HeadAtom, BodyAtom) :-
+  rule(Rule), head(Rule, HeadAtom);
+  pos_body(Rule, BodyAtom), assign(BodyAtom, false, _).
+cannot_support(Rule, HeadAtom, BodyAtom) :-
+  rule(Rule), head(Rule, HeadAtom);
+  neg_body(Rule, BodyAtom), assign(BodyAtom, true, _).
 
-{assign(Atom, false, (lack_of_support,))} :-
+assign'(Atom, false, (lack_of_support,)) :-
   atom(Atom);
-  false_body(Rule, _) : head(Rule, Atom).
+  cannot_support(Rule, HeadAtom, _) : head(Rule, Atom).
 
 last_support(Rule, Atom) :-
   assign(Atom, true, _), head(Rule, Atom);
-  #count{Rule' : head(Rule', Atom), not false_body(Rule', _)} = 1.
+  cannot_support(Rule', Atom, _) : head(Rule', Atom), Rule' != Rule.
   
-{assign(BodyAtom, true, (last_support, Rule, Atom))} :-
+assign'(BodyAtom, true, (last_support, Rule, Atom)) :-
   last_support(Rule, Atom);
   pos_body(Rule, BodyAtom).
-{assign(BodyAtom, false, (last_support, Rule, Atom))} :-
+assign'(BodyAtom, false, (last_support, Rule, Atom)) :-
   last_support(Rule, Atom);
   neg_body(Rule, BodyAtom).
+assign'(HeadAtom, false, (last_support, Rule, Atom)) :-
+  last_support(Rule, Atom), not choice(Rule, _, _);
+  head(Rule, HeadAtom), HeadAtom != Atom.
 
-constraint(Rule) :-
+%:- rule(Rule), #count{Bound : constraint(Rule, Bound)} > 1.
+constraint(Rule, upper_bound) :-
   rule(Rule), head_bounds(Rule, LowerBound, UpperBound);
   #count{Atom : head(Rule, Atom), assign(Atom, true, _)} > UpperBound.
-constraint(Rule) :-
+constraint(Rule, lower_bound) :-
   rule(Rule), head_bounds(Rule, LowerBound, UpperBound);
   #sum{1, Atom : head(Rule, Atom); -1, Atom : head(Rule, Atom), assign(Atom, false, _)} < LowerBound.
 
-{assign(Atom, false, (constraint, Rule))} :-
-  constraint(Rule), pos_body(Rule, Atom);
+assign'(Atom, false, (constraint, Rule, Bound)) :-
+  constraint(Rule, Bound), pos_body(Rule, Atom);
   assign(Atom', true, _) : pos_body(Rule, Atom'), Atom' != Atom;
   assign(Atom', false, _) : neg_body(Rule,Atom').
-{assign(Atom, true, (constraint, Rule))} :-
-  constraint(Rule), neg_body(Rule, Atom);
+assign'(Atom, true, (constraint, Rule, Bound)) :-
+  constraint(Rule, Bound), neg_body(Rule, Atom);
   assign(Atom', true, _) : pos_body(Rule, Atom');
   assign(Atom', false, _) : neg_body(Rule,Atom'), Atom' != Atom.
   
 #show.
-#show assign/3.
-#show false_body/2.
\ No newline at end of file
+#show assign(Atom, Value, Reason) : assign'(Atom, Value, Reason), not assign(Atom, _, _).
+#show cannot_support/3.
+#show constraint/2.
+#show terminate/0.
\ No newline at end of file
diff --git a/ucorexplain/encodings/mario_reify.lp b/ucorexplain/encodings/mario_reify.lp
deleted file mode 100644
index 002cff3..0000000
--- a/ucorexplain/encodings/mario_reify.lp
+++ /dev/null
@@ -1,9 +0,0 @@
-head_bounds(Rule, LowerBound, UpperBound) :-
-  rule(Rule), choice(Rule, LowerBound, UpperBound).
-head_bounds(Rule, 1, Size) :- 
-  rule(Rule), not choice(Rule,_,_); 
-  Size = #count{HeadAtom : head(Rule, HeadAtom)}.
-
-atom(Atom) :- head(Rule, Atom).
-atom(Atom) :- pos_body(Rule, Atom).
-atom(Atom) :- neg_body(Rule, Atom). 
\ No newline at end of file