diff --git a/qtensor/CircuitComposer.py b/qtensor/CircuitComposer.py
index 387f9b01..022d084e 100644
--- a/qtensor/CircuitComposer.py
+++ b/qtensor/CircuitComposer.py
@@ -51,7 +51,7 @@ def layer_of_Hadamards(self):
         for q in self.qubits:
             self.apply_gate(self.operators.H, q)
 
-    def expectation(self, operator, *qubits, **params):
+    def expectation(self, operator, *qubits, do_simplify=True, **params):
         """
         Args:
             operator: an element from OpFactory
@@ -70,14 +70,14 @@ def expectation(self, operator, *qubits, **params):
         second_part.inverse()
 
         circ = first_part.circuit + second_part.circuit
-        do_simplify = len(circ) < 100
         if do_simplify:
             try:
                 circ = qtensor.simplify_circuit.simplify_qtree_circuit(circ)
             except Exception as e:
+                print('failed to simplify:', type(e), e)
+                raise
 
                 pass
-                #print('failed to simplify:', type(e), e)
 
         return circ
 
diff --git a/qtensor/simplify_circuit/simplify_circuit_api.py b/qtensor/simplify_circuit/simplify_circuit_api.py
index 8fee4e6f..5b1cccdd 100644
--- a/qtensor/simplify_circuit/simplify_circuit_api.py
+++ b/qtensor/simplify_circuit/simplify_circuit_api.py
@@ -101,7 +101,14 @@ def simplify_qtree_circuit(qtreeCircuit):
     for gate in simplified:
         GateClass = GATE_MAP[gate.__class__]
         if issubclass(GateClass, ParametricGate):
-            qtree_gate = GateClass(*gate.index, alpha=gate.angle)
+            try:
+                qtree_gate = GateClass(*gate.index, alpha=gate.angle, beta=gate.angle)
+            except:
+                pass
+            try:
+                qtree_gate = GateClass(*gate.index, alpha=gate.angle)
+            except:
+                pass
         else:
             qtree_gate = GateClass(*gate.index)
         qtree_circuit.append(qtree_gate)
diff --git a/qtensor/toolbox.py b/qtensor/toolbox.py
index 26958f72..c8d26bb9 100644
--- a/qtensor/toolbox.py
+++ b/qtensor/toolbox.py
@@ -15,18 +15,20 @@
 import qtensor
 
 def bethe_graph(p, degree):
-    def add_two_nodes_to_leafs(graph):
+    def add_nodes_to_leafs(graph, deg=3):
         """ Works in-place """
-        leaves = [n for n in graph.nodes() if graph.degree(n) <= degree-2]
+        leaves = [n for n in graph.nodes() if graph.degree(n) <= 1]
         n = graph.number_of_nodes()
         for leaf in leaves:
-            next_edges = [(leaf, n+x) for x in range(1, degree)]
+            next_edges = [(leaf, n+x) for x in range(1, deg)]
             graph.add_edges_from(next_edges)
-            n += 2
+            n += deg-1
     graph = nx.Graph()
     graph.add_edges_from([(0,1)])
     for i in range(p):
-        add_two_nodes_to_leafs(graph)
+        add_nodes_to_leafs(graph, deg=degree)
+    n = graph.number_of_nodes()
+    assert n == 2*sum([(degree-1)**y for y in range(p+1)])
     return graph
 
 def random_graph(nodes, type='random', **kwargs):
diff --git a/qtensor/tools/qaoa_ordering/__init__.py b/qtensor/tools/qaoa_ordering/__init__.py
new file mode 100644
index 00000000..e3339ee4
--- /dev/null
+++ b/qtensor/tools/qaoa_ordering/__init__.py
@@ -0,0 +1,114 @@
+from networkx.classes.function import number_of_nodes
+import networkx as nx
+import qtensor
+
+def get_graph_order(graph: nx.Graph, opt_key:str):
+    opt = qtensor.toolbox.get_ordering_algo(opt_key)
+    peo, path = opt._get_ordering_ints(graph) #
+    return peo
+
+
+def circ2gvars(qubit_count, circuit, pdict={},
+               omit_terminals=True):
+    """
+    Constructs a graph from a circuit in the form of a
+    list of lists.
+
+    Parameters
+    ----------
+    qubit_count : int
+            number of qubits in the circuit
+    circuit : list of lists
+            quantum circuit as returned by
+            :py:meth:`operators.read_circuit_file`
+    pdict : dict
+            Dictionary with placeholders if any parameteric gates
+            were unresolved
+    max_depth : int, default None
+            Maximal depth of the circuit which should be used
+    omit_terminals : bool, default True
+            If terminal nodes should be excluded from the final
+            graph.
+
+    Returns
+    -------
+        gvars: dict(int: list)
+    """
+    import functools, itertools
+    import qtree.operators as ops
+    from qtree.optimizer import Var
+
+    # The circuit is built from left to right, as it operates
+    # on the ket ( |0> ) from the left. We thus first place
+    # the bra ( <x| ) and then put gates in the reverse order
+
+    # Fill the variable `frame`
+    layer_variables = list(range(qubit_count))
+    result = {i: [v] for i, v in enumerate(layer_variables)}
+    current_var_idx = qubit_count
+
+    # Populate nodes and save variables of the bra
+    bra_variables = []
+    for var in layer_variables:
+        bra_variables.append(Var(var, name=f"o_{var}"))
+
+    # Place safeguard measurement circuits before and after
+    # the circuit
+    measurement_circ = [[ops.M(qubit) for qubit in range(qubit_count)]]
+
+    combined_circ = functools.reduce(
+        lambda x, y: itertools.chain(x, y),
+        [measurement_circ, reversed(circuit)])
+
+    # Start building the graph in reverse order
+    for layer in combined_circ:
+        for op in layer:
+            # Update current variable frame
+            for qubit in op.changed_qubits:
+                new_var = Var(current_var_idx, name=f"o_{current_var_idx}")
+                result[qubit].append(new_var)
+                current_var_idx += 1
+
+    for qubit in range(qubit_count):
+        var = layer_variables[qubit]
+        new_var = Var(current_var_idx, name=f'i_{qubit}', size=2)
+        # update graph and variable `frame`
+        result[qubit].append(new_var)
+        current_var_idx += 1
+
+    if omit_terminals:
+        for k in result:
+            # first and last elements will always be terminals
+            result[k] = result[k][1:-1]
+
+    return result
+
+
+def build_connectivity_graph(circuit):
+    g = nx.Graph()
+    for gate in circuit:
+        if len(gate.qubits) == 2:
+            # will add muliple edges several times, but who cares
+            g.add_edge(gate.qubits[0], gate.qubits[1])
+    return g
+
+
+def get_qaoa_exp_ordering(circuit: list, algo='greedy') -> list:
+    graph = build_connectivity_graph(circuit)
+    cverts = circ2gvars(graph.number_of_nodes(), [circuit])
+    super_order = get_graph_order(graph, algo)
+    order = []
+    for i in super_order:
+        order += cverts[i]
+    return order
+
+
+class QAOAEnergyOptimizer(qtensor.optimisation.GreedyOptimizer):
+    # debt: this class is only usable for a single instance
+    def __init__(self, circuit, *args, algo='greedy', **kwargs):
+        super().__init__(*args, **kwargs)
+        self.circuit = circuit
+        self.algo = algo
+
+    def _get_ordering_ints(self, old_graph, free_vars):
+        return get_qaoa_exp_ordering(self.circuit, algo=self.algo)
\ No newline at end of file
diff --git a/qtensor/tools/qaoa_ordering/test_qaoa_ordering.py b/qtensor/tools/qaoa_ordering/test_qaoa_ordering.py
new file mode 100644
index 00000000..dac69846
--- /dev/null
+++ b/qtensor/tools/qaoa_ordering/test_qaoa_ordering.py
@@ -0,0 +1,65 @@
+import qtensor
+import qtree
+import networkx as nx
+
+
+def test_column_verts():
+    from qtensor.tools.qaoa_ordering import circ2gvars
+    g = nx.Graph()
+    g.add_edge(0, 1)
+    p = 2
+    comp = qtensor.DefaultQAOAComposer(g, gamma=[1]*p, beta=[2]*p)
+    comp.energy_expectation_lightcone((0, 1))
+    cverts = circ2gvars(g.number_of_nodes(), [comp.circuit])
+    print('cverts', cverts)
+    assert len(cverts.keys()) == g.number_of_nodes()
+    assert len(cverts[0]) == 2*p + 2 + 1
+
+    tn = qtensor.optimisation.QtreeTensorNet.from_qtree_gates(comp.circuit)
+    lg = tn.get_line_graph()
+    verts = set(lg.nodes)
+    cvert_merge = sum(cverts.values(), [])
+    cvert_merge = [int(x) for x in cvert_merge]
+    assert len(cvert_merge) == len(verts)
+    assert set(cvert_merge) == verts
+
+
+def test_qaoa_ordering():
+    from qtensor.tools.qaoa_ordering import get_qaoa_exp_ordering
+    p = 3
+    N = 50
+    g = nx.path_graph(N)  
+    comp = qtensor.DefaultQAOAComposer(g, gamma=[1]*p, beta=[2]*p)
+    comp.energy_expectation_lightcone((0, 1))
+    order = get_qaoa_exp_ordering(comp.circuit, algo='greedy')
+    tn = qtensor.optimisation.QtreeTensorNet.from_qtree_gates(comp.circuit)
+    lg = tn.get_line_graph()
+    nodes, path = qtensor.utils.get_neighbors_path(lg, peo=order)
+    width = max(path)
+    
+    assert width == 2*p + 1
+
+def test_qaoa_ordering_tree():
+    import time
+    from qtensor.tools.qaoa_ordering import get_qaoa_exp_ordering
+    p = 8
+    degree = 5
+    start = time.time()
+    g = qtensor.toolbox.bethe_graph(p, degree=degree)
+    print('P = ', p)
+    print('D = ', degree)
+    print('Bethe num of nodes', g.number_of_nodes())
+    comp = qtensor.DefaultQAOAComposer(g, gamma=[1]*p, beta=[2]*p)
+    comp.energy_expectation_lightcone((0, 1))
+    order = get_qaoa_exp_ordering(comp.circuit, algo='greedy')
+    print('Got ordering')
+    tn = qtensor.optimisation.QtreeTensorNet.from_qtree_gates(comp.circuit)
+    lg = tn.get_line_graph()
+    assert lg.number_of_nodes() == len(order)
+    print('Bethe LineGraph # of nodes', len(order))
+    nodes, path = qtensor.utils.get_neighbors_path(lg, peo=order)
+    width = max(path)
+    print('Synthetic width', width)
+    print('Total time', time.time() - start)
+    
+    assert width == 2*p + 1
\ No newline at end of file