docs: def _query(self, variables, operation, evidence=None, joint=True):
          """
          This is a generalized query method that can be used for both query and map query.
  
          Parameters
          ----------
          variables: list
              list of variables for which you want to compute the probability
          operation: str ('marginalize' | 'maximize')
              The operation to do for passing messages between nodes.
          evidence: dict
              a dict key, value pair as {var: state_of_var_observed}
              None if no evidence
  
          Examples
          --------
          >>> from pgmpy.inference import BeliefPropagation
          >>> from pgmpy.models import BayesianModel
          >>> import numpy as np
          >>> import pandas as pd
          >>> values = pd.DataFrame(np.random.randint(low=0, high=2, size=(1000, 5)),
          ...                       columns=['A', 'B', 'C', 'D', 'E'])
          >>> model = BayesianModel([('A', 'B'), ('C', 'B'), ('C', 'D'), ('B', 'E')])
          >>> model.fit(values)
          >>> inference = BeliefPropagation(model)
          >>> phi_query = inference.query(['A', 'B'])
  
          References
          ----------
          Algorithm 10.4 Out-of-clique inference in clique tree
          Probabilistic Graphical Models: Principles and Techniques Daphne Koller and Nir Friedman.
          """
  
          is_calibrated = self._is_converged(operation=operation)
          # Calibrate the junction tree if not calibrated
          if not is_calibrated:
              self.calibrate()
  
          if not isinstance(variables, (list, tuple, set)):
              query_variables = [variables]
          else:
              query_variables = list(variables)
          query_variables.extend(evidence.keys() if evidence else [])
  
          # Find a tree T' such that query_variables are a subset of scope(T')
          nodes_with_query_variables = set()
          for var in query_variables:
              nodes_with_query_variables.update(filter(lambda x: var in x, self.junction_tree.nodes()))
          subtree_nodes = nodes_with_query_variables
  
          # Conversion of set to tuple just for indexing
          nodes_with_query_variables = tuple(nodes_with_query_variables)
          # As junction tree is a tree, that means that there would be only path between any two nodes in the tree
          # thus we can just take the path between any two nodes; no matter there order is
          for i in range(len(nodes_with_query_variables) - 1):
              subtree_nodes.update(nx.shortest_path(self.junction_tree, nodes_with_query_variables[i],
                                                    nodes_with_query_variables[i + 1]))
          subtree_undirected_graph = self.junction_tree.subgraph(subtree_nodes)
          # Converting subtree into a junction tree
          if len(subtree_nodes) == 1:
              subtree = JunctionTree()
              subtree.add_node(subtree_nodes.pop())
          else:
              subtree = JunctionTree(subtree_undirected_graph.edges())
  
          # Selecting a node is root node. Root node would be having only one neighbor
          if len(subtree.nodes()) == 1:
              root_node = list(subtree.nodes())[0]
          else:
              root_node = tuple(filter(lambda x: len(list(subtree.neighbors(x))) == 1, subtree.nodes()))[0]
          clique_potential_list = [self.clique_beliefs[root_node]]
  
          # For other nodes in the subtree compute the clique potentials as follows
          # As all the nodes are nothing but tuples so simple set(root_node) won't work at it would update the set with'
          # all the elements of the tuple; instead use set([root_node]) as it would include only the tuple not the
          # internal elements within it.
          parent_nodes = set([root_node])
          nodes_traversed = set()
          while parent_nodes:
              parent_node = parent_nodes.pop()
              for child_node in set(subtree.neighbors(parent_node)) - nodes_traversed:
                  clique_potential_list.append(self.clique_beliefs[child_node] /
                                               self.sepset_beliefs[frozenset([parent_node, child_node])])
                  parent_nodes.update([child_node])
              nodes_traversed.update([parent_node])
  
          # Add factors to the corresponding junction tree
          subtree.add_factors(*clique_potential_list)
  
          # Sum product variable elimination on the subtree
          variable_elimination = VariableElimination(subtree)
          if operation == 'marginalize':
              return variable_elimination.query(variables=variables, evidence=evidence, joint=joint)
          elif operation == 'maximize':
              return variable_elimination.map_query(variables=variables, evidence=evidence)