docs: def closure(self):
          """
          Returns a new `Independencies()`-object that additionally contains those `IndependenceAssertions`
          that are implied by the the current independencies (using with the `semi-graphoid axioms
          <https://en.wikipedia.org/w/index.php?title=Conditional_independence&oldid=708760689#Rules_of_conditional_independence>`_;
          see (Pearl, 1989, `Conditional Independence and its representations
          <http://www.cs.technion.ac.il/~dang/journal_papers/pearl1989conditional.pdf>`_)).
  
          Might be very slow if more than six variables are involved.
  
          Examples
          --------
          >>> from pgmpy.independencies import Independencies
          >>> ind1 = Independencies(('A', ['B', 'C'], 'D'))
          >>> ind1.closure()
          (A _|_ B | D, C)
          (A _|_ B, C | D)
          (A _|_ B | D)
          (A _|_ C | D, B)
          (A _|_ C | D)
  
          >>> ind2 = Independencies(('W', ['X', 'Y', 'Z']))
          >>> ind2.closure()
          (W _|_ Y)
          (W _|_ Y | X)
          (W _|_ Z | Y)
          (W _|_ Z, X, Y)
          (W _|_ Z)
          (W _|_ Z, X)
          (W _|_ X, Y)
          (W _|_ Z | X)
          (W _|_ Z, Y | X)
          [..]
          """
  
          def single_var(var):
              "Checks if var represents a single variable"
              if not hasattr(var, '__iter__'):
                  return True
              else:
                  return len(var) == 1
  
          def sg0(ind):
              "Symmetry rule: 'X ⟂ Y | Z' -> 'Y ⟂ X | Z'"
              return IndependenceAssertion(ind.event2, ind.event1, ind.event3)
  
          # since X⟂Y|Z == Y⟂X|Z in pgmpy, sg0 (symmetry) is not used as an axiom/rule.
          # instead we use a decorator for the other axioms to apply them on both sides
          def apply_left_and_right(func):
              def symmetric_func(*args):
                  if len(args) == 1:
                      return func(args[0]) + func(sg0(args[0]))
                  if len(args) == 2:
                      return (func(*args) + func(args[0], sg0(args[1])) +
                              func(sg0(args[0]), args[1]) + func(sg0(args[0]), sg0(args[1])))
              return symmetric_func
  
          @apply_left_and_right
          def sg1(ind):
              "Decomposition rule: 'X ⟂ Y,W | Z' -> 'X ⟂ Y | Z', 'X ⟂ W | Z'"
              if single_var(ind.event2):
                  return []
              else:
                  return [IndependenceAssertion(ind.event1, ind.event2 - {elem}, ind.event3)
                          for elem in ind.event2]
  
          @apply_left_and_right
          def sg2(ind):
              "Weak Union rule: 'X ⟂ Y,W | Z' -> 'X ⟂ Y | W,Z', 'X ⟂ W | Y,Z' "
              if single_var(ind.event2):
                  return []
              else:
                  return [IndependenceAssertion(ind.event1, ind.event2 - {elem}, {elem} | ind.event3)
                          for elem in ind.event2]
  
          @apply_left_and_right
          def sg3(ind1, ind2):
              "Contraction rule: 'X ⟂ W | Y,Z' & 'X ⟂ Y | Z' -> 'X ⟂ W,Y | Z'"
              if ind1.event1 != ind2.event1:
                  return []
  
              Y = ind2.event2
              Z = ind2.event3
              Y_Z = ind1.event3
              if Y < Y_Z and Z < Y_Z and Y.isdisjoint(Z):
                  return [IndependenceAssertion(ind1.event1, ind1.event2 | Y, Z)]
              else:
                  return []
  
          # apply semi-graphoid axioms as long as new independencies are found.
          all_independencies = set()
          new_inds = set(self.independencies)
  
          while new_inds:
              new_pairs = (set(itertools.permutations(new_inds, 2)) |
                           set(itertools.product(new_inds, all_independencies)) |
                           set(itertools.product(all_independencies, new_inds)))
  
              all_independencies |= new_inds
              new_inds = set(sum([sg1(ind) for ind in new_inds] +
                                 [sg2(ind) for ind in new_inds] +
                                 [sg3(*inds) for inds in new_pairs], []))
              new_inds -= all_independencies
  
          return Independencies(*list(all_independencies))