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))