code: def _elimination_trees(theta, decision_variables):
      
      # auxiliary variables are any variables that are not decision
      auxiliary_variables = set(n for n in theta.linear if n not in decision_variables)
  
      # get the adjacency of the auxiliary subgraph
      adj = {v: {u for u in theta.adj[v] if u in auxiliary_variables}
             for v in theta.adj if v in auxiliary_variables}
  
      # get the elimination order that minimizes treewidth
      tw, order = dnx.treewidth_branch_and_bound(adj)
  
      ancestors = {}
      for n in order:
          ancestors[n] = set(adj[n])
  
          # now make v simplicial by making its neighborhood a clique, then
          # continue
          neighbors = adj[n]
          for u, v in itertools.combinations(neighbors, 2):
              adj[u].add(v)
              adj[v].add(u)
          for v in neighbors:
              adj[v].discard(n)
          del adj[n]
  
      roots = {}
      nodes = {v: {} for v in ancestors}
      for vidx in range(len(order) - 1, -1, -1):
          v = order[vidx]
  
          if ancestors[v]:
              for u in order[vidx + 1:]:
                  if u in ancestors[v]:
                      # v is a child of u
                      nodes[u][v] = nodes[v]  # nodes[u][v] = children of v
                      break
          else:
              roots[v] = nodes[v]  # roots[v] = children of v
  
      return roots, ancestors
docs: From Theta and the decision variables, determine the elimination order and the induced
      trees.