From d9b6bf3316af5a79999449b1879a852534c5f947 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Wed, 21 Aug 2024 16:26:39 +0200
Subject: [PATCH 01/31] provide is_constant

---
 src/sage/rings/species.py | 231 +++++++++++++++++++-------------------
 1 file changed, 113 insertions(+), 118 deletions(-)

diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index 4e07ef076be..41eea7498f7 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -1,4 +1,5 @@
 from itertools import accumulate, chain
+
 from sage.categories.graded_algebras_with_basis import GradedAlgebrasWithBasis
 from sage.categories.infinite_enumerated_sets import InfiniteEnumeratedSets
 from sage.categories.monoids import Monoids
@@ -10,12 +11,14 @@
 from sage.groups.perm_gps.permgroup_named import SymmetricGroup
 from sage.libs.gap.libgap import libgap
 from sage.misc.cachefunc import cached_method
-from sage.monoids.indexed_free_monoid import IndexedFreeAbelianMonoid, IndexedFreeAbelianMonoidElement, IndexedMonoid
+from sage.monoids.indexed_free_monoid import (IndexedFreeAbelianMonoid,
+                                              IndexedFreeAbelianMonoidElement,
+                                              IndexedMonoid)
 from sage.rings.integer_ring import ZZ
+from sage.sets.finite_enumerated_set import FiniteEnumeratedSet
 from sage.structure.element import Element, parent
 from sage.structure.parent import Parent
 from sage.structure.unique_representation import UniqueRepresentation
-from sage.sets.finite_enumerated_set import FiniteEnumeratedSet
 
 GAP_FAIL = libgap.eval('fail')
 
@@ -43,25 +46,6 @@ def __init__(self):
         """
         self._cache = dict()
 
-    def clear_cache(self):
-        r"""
-        Clear the cache.
-
-        EXAMPLES::
-
-            sage: A = AtomicSpecies(1)
-            sage: A(SymmetricGroup(1))
-            {((),): ((1,),)}
-            sage: A(SymmetricGroup(0))
-            {(): ((),)}
-            sage: A._cache
-            {((0,), ()): [{(): ((),)}], ((1,), ((),)): [{((),): ((1,),)}]}
-            sage: A.clear_cache()
-            sage: A._cache
-            {}
-        """
-        self._cache = dict()
-
     def _cache_get(self, elm):
         r"""
         Return the cached element, or add it
@@ -225,18 +209,12 @@ def _element_constructor_(self, x):
 
             sage: from sage.rings.species import ConjugacyClassesOfDirectlyIndecomposableSubgroups
             sage: C = ConjugacyClassesOfDirectlyIndecomposableSubgroups()
-            sage: C.clear_cache()
             sage: C(PermutationGroup([("a", "b", "c")]))
             ((1,2,3),)
-            sage: C._cache
-            {(3, 3, (3,)): [((1,2,3),)]}
             sage: C(PermutationGroup([(1, 3, 5)], domain=[1,3,5]))
             ((1,2,3),)
             sage: C(PermutationGroup([[(1,3),(4,7)],[(2,5),(6,8)], [(1,4),(2,5),(3,7)]]))
             ((5,6)(7,8), (1,2)(3,4), (1,3)(2,4)(5,6))
-            sage: C._cache
-            {(3, 3, (3,)): [((1,2,3),)],
-            (8, 8, (2, 2, 4)): [((5,6)(7,8), (1,2)(3,4), (1,3)(2,4)(5,6))]}
         """
         if parent(x) == self:
             return x
@@ -329,7 +307,7 @@ def __init__(self, parent, dis, domain_partition):
 
         TESTS::
 
-            sage: A = AtomicSpecies(1)
+            sage: A = AtomicSpecies("X")
             sage: G = PermutationGroup([[(1,3),(4,7)], [(2,5),(6,8)], [(1,4),(2,5),(3,7)]])
             sage: TestSuite(A(G)).run()
         """
@@ -357,7 +335,7 @@ def _canonicalize(self):
 
         EXAMPLES::
 
-            sage: At = AtomicSpecies(2)
+            sage: At = AtomicSpecies("X, Y")
             sage: G = PermutationGroup([[(1,2),(3,4),(5,6),(7,8,9,10)]]); G
             Permutation Group with generators [(1,2)(3,4)(5,6)(7,8,9,10)]
             sage: H = PermutationGroup([[(1,2,3,4),(5,6),(7,8),(9,10)]]); H
@@ -378,7 +356,7 @@ def __hash__(self):
 
         TESTS::
 
-            sage: At = AtomicSpecies(2)
+            sage: At = AtomicSpecies("X, Y")
             sage: G = PermutationGroup([[(1,2),(3,4),(5,6),(7,8,9,10)]]); G
             Permutation Group with generators [(1,2)(3,4)(5,6)(7,8,9,10)]
             sage: H = PermutationGroup([[(1,2,3,4),(5,6),(7,8),(9,10)]]); H
@@ -403,7 +381,7 @@ def __eq__(self, other):
 
         TESTS::
 
-            sage: At = AtomicSpecies(2)
+            sage: At = AtomicSpecies("X, Y")
             sage: G = PermutationGroup([[(1,2),(3,4),(5,6),(7,8,9,10)]]); G
             Permutation Group with generators [(1,2)(3,4)(5,6)(7,8,9,10)]
             sage: H = PermutationGroup([[(1,2,3,4),(5,6),(7,8),(9,10)]]); H
@@ -440,7 +418,7 @@ def _repr_(self):
 
         TESTS::
 
-            sage: At = AtomicSpecies(2)
+            sage: At = AtomicSpecies("X, Y")
             sage: G = PermutationGroup([[(1,2),(3,4),(5,6),(7,8,9,10)]]); G
             Permutation Group with generators [(1,2)(3,4)(5,6)(7,8,9,10)]
             sage: A = At(G, {1: [1,2,3,4], 2: [5,6,7,8,9,10]}); A
@@ -459,13 +437,16 @@ def __classcall__(cls, names):
 
             sage: A1 = AtomicSpecies("X")
             sage: A2 = AtomicSpecies("Y")
+            sage: A3 = AtomicSpecies("X, Y")
+            sage: A4 = AtomicSpecies(["X", "Y"])
             sage: A1 is A2
             False
+            sage: A3 is A4
+            True
         """
-        if all(isinstance(X, str) for X in names):
-            raise ValueError(f"singleton_names must be a tuple of {k} strings")
-
-        return super().__classcall__(cls, k, singleton_names)
+        from sage.structure.category_object import normalize_names
+        names = normalize_names(-1, names)
+        return super().__classcall__(cls, names)
 
     def __init__(self, names):
         r"""
@@ -483,11 +464,10 @@ def __init__(self, names):
             sage: TestSuite(At2).run(skip="_test_graded_components")
         """
         category = SetsWithGrading().Infinite()
-        Parent.__init__(self, category=category)
+        Parent.__init__(self, names=names, category=category)
         ElementCache.__init__(self)
-        self._k = k
-        self._grading_set = IntegerVectors(length=k)
-        self._singleton_names = singleton_names
+        self._k = len(names)
+        self._grading_set = IntegerVectors(length=self._k)
         self._renamed = set()  # the degrees that have been renamed already
 
     @cached_method
@@ -497,10 +477,10 @@ def an_element(self):
 
         TESTS::
 
-            sage: At1 = AtomicSpecies(1)
-            sage: At2 = AtomicSpecies(2)
+            sage: At1 = AtomicSpecies("X")
+            sage: At2 = AtomicSpecies("X, Y")
             sage: At1.an_element()
-            {((1,2),): ((1, 2),)}
+            E_2
             sage: At2.an_element()
             {((1,2)(3,4),): ((1, 2), (3, 4))}
         """
@@ -548,27 +528,27 @@ def _element_constructor_(self, G, pi=None):
         for k, v in pi.items():
             dpart[k - 1] = tuple(mapping2(mapping[x]) for x in v)
         elm = self._cache_get(self.element_class(self, dis_elm, tuple(dpart)))
-        if self._singleton_names and elm._tc not in self._renamed:
+        if elm._tc not in self._renamed:
             self._rename(elm._tc)
         return elm
 
     def _rename(self, n):
         from sage.groups.perm_gps.permgroup import PermutationGroup
-        from sage.groups.perm_gps.permgroup_named import (SymmetricGroup,
-                                                            CyclicPermutationGroup,
-                                                            DihedralGroup,
-                                                            AlternatingGroup)
+        from sage.groups.perm_gps.permgroup_named import (AlternatingGroup,
+                                                          CyclicPermutationGroup,
+                                                          DihedralGroup,
+                                                          SymmetricGroup)
 
         # prevent infinite recursion in self._element_constructor_
         self._renamed.add(n)
         for i in range(self._k):
             if n == 1:
-                self(SymmetricGroup(1), {i+1: [1]}).rename(self._singleton_names[i])
+                self(SymmetricGroup(1), {i+1: [1]}).rename(self._names[i])
 
             if self._k == 1:
                 sort = ""
             else:
-                sort = f"({self._singleton_names[i]})"
+                sort = f"({self._names[i]})"
 
             if n >= 2:
                 self(SymmetricGroup(n),
@@ -625,12 +605,14 @@ def _repr_(self):
 
         TESTS::
 
-            sage: At1 = AtomicSpecies(1); At1
-            Infinite set of 1-variate atomic species
-            sage: At2 = AtomicSpecies(2); At2
-            Infinite set of 2-variate atomic species
+            sage: AtomicSpecies("X")
+            Atomic species in X
+            sage: AtomicSpecies("X, Y")
+            Atomic species in X, Y
         """
-        return f"Atomic species in sorts "
+        if len(self._names) == 1:
+            return f"Atomic species in {self._names[0]}"
+        return f"Atomic species in {', '.join(self._names)}"
 
     Element = AtomicSpeciesElement
 
@@ -704,7 +686,7 @@ def one(self):
 
         EXAMPLES::
 
-            sage: P = PolynomialSpecies(2)
+            sage: P = PolynomialSpecies(ZZ, "X, Y")
             sage: P._indices.one()
             1
         """
@@ -853,7 +835,7 @@ def _canonicalize(self):
 
             EXAMPLES::
 
-                sage: P = PolynomialSpecies(2)
+                sage: P = PolynomialSpecies(ZZ, "X, Y")
                 sage: G = PermutationGroup([[(1,2),(3,4),(5,6),(7,8,9,10)]]); G
                 Permutation Group with generators [(1,2)(3,4)(5,6)(7,8,9,10)]
                 sage: H = PermutationGroup([[(1,2,3,4),(5,6),(7,8),(9,10)]]); H
@@ -892,7 +874,7 @@ def hadamard_product(self, other):
 
             Exercise 2.1.9 from the BLL book::
 
-                sage: P = PolynomialSpecies(["X"])
+                sage: P = PolynomialSpecies(ZZ, ["X"])
                 sage: M = P._indices
                 sage: C3 = M(CyclicPermutationGroup(3))
                 sage: X = M(SymmetricGroup(1))
@@ -904,16 +886,11 @@ def hadamard_product(self, other):
                 sage: (X*E2).hadamard_product(X*E2)
                 X*E_2 + X^3
             """
-            if not isinstance(other, MolecularSpecies.Element):
-                raise ValueError("other must be a molecular species")
-            if self.parent()._k != other.parent()._k:
-                raise ValueError("other must have same arity")
-
             P = self.parent()
-            if not P is other.parent():
+            if P is not other.parent():
                 raise ValueError("the factors of a Hadamard product must be the same.")
- 
-            Pn = PolynomialSpecies(names=P._indices.names)
+            Pn = PolynomialSpecies(ZZ, P._indices._names)
+
             if self._mc != other._mc:
                 return Pn.zero()
             # create S
@@ -934,31 +911,40 @@ def hadamard_product(self, other):
                 res += Pn(PermutationGroup(gap_group=F, domain=self.domain()), dpart)
             return res
 
-        def inner_sum(self, *args):
+        def inner_sum(self, base_ring, names, *args):
             r"""
             Compute the inner sum of exercise 2.6.16 of BLL book.
 
-            args are the compositions (in Compositions) each of which
-            sum to the corresponding cardinality of ``self``. The number
-            of args is equal to the arity of ``self``.
+            INPUT:
+
+                - ``base_ring``, the base ring of the result
+
+                - ``names``, the names of the result
+
+                - ``args``, the sequence of compositions, each of
+                  which sums to the corresponding cardinality of
+                  ``self``. The number of ``args`` is equal to the
+                  arity of ``self``.
 
             EXAMPLES::
 
-                sage: P = PolynomialSpecies(1)
+                sage: P = PolynomialSpecies(ZZ, "X")
                 sage: M = P._indices
                 sage: C4 = M(CyclicPermutationGroup(4))
-                sage: C4.inner_sum([2, 2]) # X^2Y^2 + C2(XY)
-                {((),): ((1,), ())}^2*{((),): ((), (1,))}^2 + {((1,2)(3,4),): ((1, 2), (3, 4))}
+                sage: C4.inner_sum(ZZ, "X, Y", [2, 2]) # X^2Y^2 + C2(XY)
+                {((1,2)(3,4),): ((1, 2), (3, 4))} + X^2*Y^2
+
             """
             # TODO: No checks are performed right now, must be added.
             # Checks: all args in compositions, sums must match cardinalities.
 
-            res = 0
+            # Create group of the composition
+            Pn = PolynomialSpecies(base_ring, names)
+            res = Pn.zero()
             # conjugate self._group so that [1..k] is sort 1, [k+1,..] is sort 2, so on
             conj = PermutationGroupElement(list(chain.from_iterable(self._dompart))).inverse()
             G = libgap.ConjugateGroup(self._group, conj)
-            # Create group of the composition
-            Pn = PolynomialSpecies(len(args[0]))
+
             comp = list(chain.from_iterable(args))
             # Create the double coset representatives.
             S_down = SymmetricGroup(sum(comp)).young_subgroup(comp)
@@ -980,7 +966,7 @@ def __call__(self, *args):
 
             EXAMPLES::
 
-                sage: P = PolynomialSpecies(["X"])
+                sage: P = PolynomialSpecies(ZZ, ["X"])
                 sage: M = P._indices
                 sage: X = M(SymmetricGroup(1))
                 sage: E2 = M(SymmetricGroup(2))
@@ -1029,32 +1015,23 @@ def __call__(self, *args):
 
 
 class PolynomialSpecies(CombinatorialFreeModule):
-    def __classcall__(cls, names, base_ring=ZZ):
+    def __classcall__(cls, base_ring, names):
         r"""
         Normalize the arguments.
 
         TESTS::
 
-            sage: P1 = PolynomialSpecies("X", ZZ)
-            sage: P2 = PolynomialSpecies(1, base_ring=ZZ)
-            sage: P3 = PolynomialSpecies(["X", "Y"], base_ring=ZZ)
-            sage: P4 = PolynomialSpecies(2, ["X", "Y"])
+            sage: P1 = PolynomialSpecies(ZZ, "X, Y")
+            sage: P2 = PolynomialSpecies(ZZ, "X, Y")
+            sage: P3 = PolynomialSpecies(ZZ, ["X", "Z"])
             sage: P1 is P2
-            False
-            sage: P3 is P4
             True
-
-        .. TODO::
-
-            Reconsider the order of the arguments.
-            ``singleton_names`` are not generators, and may even be
-            omitted, and ``base_ring`` will usually be ``ZZ``, but
-            maybe it would be nice to allow ``P.<X,Y,Z> =
-            PolynomialSpecies(ZZ)`` anyway.
-
+            sage: P1 is P3
+            False
         """
-        A = AtomicSpecies(k, singleton_names=singleton_names)
-        return super().__classcall__(cls, A._k, A._singleton_names, base_ring)
+        from sage.structure.category_object import normalize_names
+        names = normalize_names(-1, names)
+        return super().__classcall__(cls, base_ring, names)
 
     def __init__(self, base_ring, names):
         r"""
@@ -1062,21 +1039,21 @@ def __init__(self, base_ring, names):
 
         TESTS::
 
-            sage: P = PolynomialSpecies(["X"])
+            sage: P = PolynomialSpecies(ZZ, "X")
             sage: TestSuite(P).run()
-            sage: P2 = PolynomialSpecies(["X", "Y"])
+            sage: P2 = PolynomialSpecies(ZZ, "X, Y")
             sage: TestSuite(P2).run()
         """
         # should we pass a category to basis_keys?
-        A = AtomicSpecies(names=names)
+        A = AtomicSpecies(names)
         basis_keys = MolecularSpecies(A, prefix='', bracket=False)
         category = GradedAlgebrasWithBasis(base_ring).Commutative()
         CombinatorialFreeModule.__init__(self, base_ring,
-                                        basis_keys=basis_keys,
-                                        category=category,
-                                        element_class=self.Element,
-                                        prefix='', bracket=False)
-        self._k = k
+                                         basis_keys=basis_keys,
+                                         category=category,
+                                         element_class=self.Element,
+                                         prefix='', bracket=False)
+        self._k = len(names)
         self._atomic_basis = basis_keys.indices()
 
     def degree_on_basis(self, m):
@@ -1114,10 +1091,10 @@ def one_basis(self):
 
         EXAMPLES::
 
-            sage: P = PolynomialSpecies(1)
+            sage: P = PolynomialSpecies(ZZ, "X")
             sage: P.one_basis()
             1
-            sage: P2 = PolynomialSpecies(2)
+            sage: P2 = PolynomialSpecies(ZZ, "X, Y")
             sage: P2.one_basis()
             1
         """
@@ -1128,10 +1105,10 @@ def an_element(self):
         """
         Return an element of ``self``.
 
-            sage: P = PolynomialSpecies(1)
+            sage: P = PolynomialSpecies(ZZ, "X")
             sage: P.an_element()
             1
-            sage: P2 = PolynomialSpecies(2)
+            sage: P2 = PolynomialSpecies(ZZ, "X, Y")
             sage: P2.an_element()
             1
         """
@@ -1143,7 +1120,7 @@ def product_on_basis(self, H, K):
 
         EXAMPLES::
 
-            sage: P = PolynomialSpecies("X")
+            sage: P = PolynomialSpecies(ZZ, "X")
             sage: L1 = [P(H) for H in SymmetricGroup(3).conjugacy_classes_subgroups()]
             sage: L2 = [P(H) for H in SymmetricGroup(2).conjugacy_classes_subgroups()]
             sage: matrix([[F * G for F in L1] for G in L2])
@@ -1158,23 +1135,41 @@ def _repr_(self):
 
         EXAMPLES::
 
-            sage: P = PolynomialSpecies(1)
-            sage: P
-            Ring of 1-variate virtual species
-            sage: P2 = PolynomialSpecies(2)
-            sage: P2
-            Ring of 2-variate virtual species
+            sage: PolynomialSpecies(ZZ, "X")
+            Polynomial species in X over Integer Ring
+            sage: PolynomialSpecies(ZZ, "X, Y")
+            Polynomial species in X, Y over Integer Ring
         """
-        return f"Ring of {self._k}-variate virtual species"
+        names = self._indices._indices._names
+        if len(names) == 1:
+            return f"Polynomial species in {names[0]} over {self.base_ring()}"
+        return f"Polynomial species in {', '.join(names)} over {self.base_ring()}"
 
     class Element(CombinatorialFreeModule.Element):
+        def is_constant(self):
+            """
+            Return ``True`` if this is a constant polynomial species.
+    
+            EXAMPLES::
+    
+                sage: P = PolynomialSpecies(ZZ, ["X", "Y"])
+                sage: X = P(SymmetricGroup(1), {1: [1]})
+                sage: X.is_constant()
+                False
+                sage: (3*P.one()).is_constant()
+                True
+                sage: P(0).is_constant()
+                True
+            """
+            return self.is_zero() or not self.degree()
+        
         def is_virtual(self):
             r"""
             Return if ``self`` is a virtual species.
 
             TESTS::
 
-                sage: P = PolynomialSpecies(["X", "Y"])
+                sage: P = PolynomialSpecies(ZZ, ["X", "Y"])
                 sage: X = P(SymmetricGroup(1), {1: [1]})
                 sage: Y = P(SymmetricGroup(1), {2: [1]})
                 sage: V = 2 * X - 3 * Y; V
@@ -1192,7 +1187,7 @@ def is_molecular(self):
 
             TESTS::
 
-                sage: P = PolynomialSpecies(["X", "Y"])
+                sage: P = PolynomialSpecies(ZZ, ["X", "Y"])
                 sage: X = P(SymmetricGroup(1), {1: [1]})
                 sage: Y = P(SymmetricGroup(1), {2: [1]})
                 sage: V = 2 * X - 3 * Y; V
@@ -1212,7 +1207,7 @@ def is_atomic(self):
 
             TESTS::
 
-                sage: P = PolynomialSpecies(["X", "Y"])
+                sage: P = PolynomialSpecies(ZZ, ["X", "Y"])
                 sage: X = P(SymmetricGroup(1), {1: [1]})
                 sage: Y = P(SymmetricGroup(1), {2: [1]})
                 sage: V = 2 * X - 3 * Y; V

From 9e31e27711d2fc4eecf81bd01d53e58d7d9d7be9 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Wed, 21 Aug 2024 16:26:50 +0200
Subject: [PATCH 02/31] basic design, non-working __call__

---
 src/sage/rings/lazy_species.py | 344 +++++++++++++++++++++++++++++++++
 1 file changed, 344 insertions(+)
 create mode 100644 src/sage/rings/lazy_species.py

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
new file mode 100644
index 00000000000..79378278ab4
--- /dev/null
+++ b/src/sage/rings/lazy_species.py
@@ -0,0 +1,344 @@
+from sage.rings.integer_ring import ZZ
+from sage.rings.lazy_series import LazyCompletionGradedAlgebraElement, LazyModuleElement
+from sage.rings.lazy_series_ring import LazyCompletionGradedAlgebra
+from sage.data_structures.stream import (
+    Stream_zero,
+    Stream_exact,
+    Stream_function,
+    Stream_cauchy_compose,
+)
+from sage.rings.species import AtomicSpecies, PolynomialSpecies
+from sage.libs.gap.libgap import libgap
+from sage.categories.sets_cat import cartesian_product
+from sage.combinat.set_partition import SetPartitions
+from sage.structure.element import parent
+import itertools
+
+
+def weighted_compositions(n, d, weights):
+    r"""
+    Return all compositions of `n` of weight `d`.
+
+    The weight of a composition `n_1, n_2, \dots` is `\sum_i w_i
+    n_i`.
+
+    ``weights`` is assumed to be a weakly increasing list of positive
+    integers.
+
+    EXAMPLES::
+
+        sage: from sage.rings.lazy_species import weighted_compositions
+        sage: list(weighted_compositions(1, 1, [1,1,2]))
+        [[0, 1], [1]]
+
+        sage: list(weighted_compositions(2, 1, [1,1,2]))
+        []
+
+        sage: list(weighted_compositions(1, 2, [1,1,2,3]))
+        [[0, 0, 1]]
+
+        sage: list(weighted_compositions(3, 4, [1,1,2,2,3,3,4,4,5]))
+        [[0, 2, 0, 1], [0, 2, 1], [1, 1, 0, 1], [1, 1, 1], [2, 0, 0, 1], [2, 0, 1]]
+
+    """
+    n = int(n)
+    d = int(d)
+    # the empty composition exists if and only if n == d == 0
+    if not n:
+        if not d:
+            yield []
+        return
+    if not d:
+        return
+
+    # otherwise we iterate over the possibilities for the first part
+    w0 = weights[0]
+    if w0 > d:
+        return
+    wr = weights[1:]
+    for i in range(min(n, d // w0) + 1):
+        for c in weighted_compositions(n - i, d - i * w0, wr):
+            yield [i] + c
+
+
+######################################################################
+
+class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
+    r"""
+    Compute the molecular expansion of `E(-X)`::
+
+        sage: from sage.rings.lazy_species import LazySpecies
+        sage: L = LazySpecies(ZZ, "X")
+        sage: E = L(lambda n: SymmetricGroup(n))
+        sage: 1 / E
+        1 + (-X) + (-E_2+X^2) + (-E_3+2*X*E_2-X^3)
+          + (-E_4+2*X*E_3+E_2^2-3*X^2*E_2+X^4)
+          + (-E_5+2*X*E_4+2*E_2*E_3-3*X^2*E_3-3*X*E_2^2+4*X^3*E_2-X^5)
+          + (-E_6+2*X*E_5+2*E_2*E_4-3*X^2*E_4+E_3^2-6*X*E_2*E_3+4*X^3*E_3-E_2^3+6*X^2*E_2^2-5*X^4*E_2+X^6)
+          + O^7
+
+    Compare with the explicit formula::
+
+        sage: def coefficient(m):
+        ....:     return sum((-1)^len(la) * multinomial((n := la.to_exp())) * prod(E[i]^ni for i, ni in enumerate(n, 1)) for la in Partitions(m))
+
+        sage: all(coefficient(m) == (1/E)[m] for m in range(10))
+        True
+    """
+    def generating_series(self):
+        pass
+
+    def isotype_generating_series(self):
+        pass
+
+    def cycle_index_series(self):
+        pass
+
+    def _add_(self, other):
+        r"""
+        Return the sum of ``self`` and ``other``.
+
+        In particular, the method to obtain the structures is
+        adapted.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: F = L(lambda n: SymmetricGroup(n))
+            sage: list(E.structures([1,2,3]))
+            [((1, 2, 3), E_3)]
+            sage: list((E+F).structures([1,2,3]))
+            [((1, 2, 3), E_3), ((1, 2, 3), E_3)]
+
+        """
+        def add_structures(labels):
+            yield from self.structures(labels)
+            yield from other.structures(labels)
+
+        result = super()._add_(other)
+        result.structures = add_structures
+        return result
+
+    def _mul_(self, other):
+        """
+        Return the product of this series with ``other``.
+
+        In particular, the method to obtain the structures is
+        adapted.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: list((E^2).structures([1,2,3]))
+            [(((), 1), ((1, 2, 3), E_3)),
+             (((1,), X), ((2, 3), E_2)),
+             (((2,), X), ((1, 3), E_2)),
+             (((3,), X), ((1, 2), E_2)),
+             (((1, 2), E_2), ((3,), X)),
+             (((1, 3), E_2), ((2,), X)),
+             (((2, 3), E_2), ((1,), X))]
+        """
+        def mul_structures(labels):
+            n = len(labels)
+            l = set(labels)
+            assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
+            for k in range(n):
+                for U in itertools.combinations(l, k):
+                    V = l.difference(U)
+                    yield from itertools.product(self.structures(U),
+                                                 other.structures(V))
+
+        result = super()._mul_(other)
+        result.structures = mul_structures
+        return result
+
+    def structures(self, labels):
+        r"""
+        Iterate over the structures on the given set of labels.
+
+        Generically, this yields a list of relabelled representatives
+        of the cosets of corresponding groups.
+
+        The relabelling is such that the first few labels correspond
+        to the first factor in the atomic decomposition, etc.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: list(E.structures([1,2,3]))
+            [((1, 2, 3), E_3)]
+
+            sage: P = L(lambda n: CyclicPermutationGroup(n))
+            sage: list(P.structures([1,2,3]))
+            [((1, 2, 3), C_3), ((2, 1, 3), C_3)]
+
+            sage: F = 1/(2-E)
+            sage: list(F.structures([1,2,3]))
+            [((1, 2, 3), E_3),
+             ((1, 2, 3), X*E_2, 0),
+             ((3, 2, 1), X*E_2, 0),
+             ((2, 3, 1), X*E_2, 0),
+             ((1, 2, 3), X*E_2, 1),
+             ((3, 2, 1), X*E_2, 1),
+             ((2, 3, 1), X*E_2, 1),
+             ((1, 2, 3), X^3),
+             ((3, 2, 1), X^3),
+             ((3, 1, 2), X^3),
+             ((2, 1, 3), X^3),
+             ((2, 3, 1), X^3),
+             ((1, 3, 2), X^3)]
+        """
+        from sage.groups.perm_gps.permgroup_named import SymmetricGroup
+        n = len(labels)
+        l = set(labels)
+        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
+        S = SymmetricGroup(n)
+        F = self[n]
+        l = list(l)[::-1]
+        for M, c in F.monomial_coefficients().items():
+            if c < 0:
+                raise NotImplementedError("only implemented for proper non-virtual species")
+            types = [tuple(S(rep)._act_on_list_on_position(l))[::-1]
+                     for rep in libgap.RightTransversal(S, M._group)]
+            if c == 1:
+                for s in types:
+                    yield s, M
+            else:
+                for e, s in cartesian_product([range(c), types]):
+                    yield s, M, e
+
+    def isotypes(self, labels):
+        pass
+
+    def polynomial(self, degree=None, names=None):
+        r"""
+        Return ``self`` as a polynomial if ``self`` is actually so.
+
+        INPUT:
+
+        - ``degree`` -- ``None`` or an integer
+        - ``names`` -- names of the variables; if it is ``None``, the name of
+          the variables of the series is used
+
+        OUTPUT:
+
+        If ``degree`` is not ``None``, the terms of the series of
+        degree greater than ``degree`` are first truncated.  If
+        ``degree`` is ``None`` and the series is not a polynomial
+        polynomial, a ``ValueError`` is raised.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: E.polynomial(3)
+            1 + X + E_2 + E_3
+        """
+        S = self.parent()
+        R = S._laurent_poly_ring
+
+        if isinstance(self._coeff_stream, Stream_zero):
+            return R.zero()
+
+        if degree is None:
+            if (isinstance(self._coeff_stream, Stream_exact)
+                and not self._coeff_stream._constant):
+                m = self._coeff_stream._degree
+            else:
+                raise ValueError("not a polynomial species")
+        else:
+            m = degree + 1
+
+        return R.sum(self[:m])
+
+    def __call__(self, *g):
+        """
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: L(SymmetricGroup(2))(L(SymmetricGroup(2)))
+            P_4
+        """
+        fP = parent(self)
+        if len(g) != fP._arity:
+            raise ValueError("arity of must be equal to the number of arguments provided")
+
+        # Find a good parent for the result
+        from sage.structure.element import get_coercion_model
+        cm = get_coercion_model()
+        P = cm.common_parent(self.base_ring(), *[parent(h) for h in g])
+
+        # f = 0
+        if isinstance(self._coeff_stream, Stream_zero):
+            return P.zero()
+
+        # g = (0, ..., 0)
+        if all((not isinstance(h, LazyModuleElement) and not h)
+               or (isinstance(h, LazyModuleElement)
+                   and isinstance(h._coeff_stream, Stream_zero))
+               for h in g):
+            return P(self[0])
+
+        # f has finite length and f != 0
+        if (isinstance(self._coeff_stream, Stream_exact)
+            and not self._coeff_stream._constant):
+            # constant polynomial
+            poly = self.polynomial()
+            if poly.is_constant():
+                return P(poly)
+            return P(poly(g))
+
+        g = [P(h) for h in g]
+        R = P._internal_poly_ring.base_ring()
+
+        for h in g:
+            if h._coeff_stream._approximate_order == 0:
+                if not h._coeff_stream.is_uninitialized() and h[0]:
+                    raise ValueError("can only compose with a positive valuation series")
+                h._coeff_stream._approximate_order = 1
+
+        sorder = self._coeff_stream._approximate_order
+        gv = min(h._coeff_stream._approximate_order for h in g)
+
+        def coefficient(n):
+            r = R.zero()
+            for i in range(n // gv + 1):
+                # Make sure the element returned from the composition is in P
+                r += P(self[i](g))[n]
+            return r
+        coeff_stream = Stream_function(coefficient, P._sparse, sorder * gv)
+        return P.element_class(P, coeff_stream)
+
+
+class LazySpecies(LazyCompletionGradedAlgebra):
+    """
+    The ring of combinatorial species.
+    """
+    Element = LazySpeciesElement
+
+    @staticmethod
+    def __classcall_private__(cls, base_ring, names, sparse=True):
+        """
+        Normalize input to ensure a unique representation.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: LazySpecies(QQ, "X")
+            Lazy completion of Polynomial species in X over Rational Field
+        """
+        from sage.structure.category_object import normalize_names
+        names = normalize_names(-1, names)
+        return super().__classcall__(cls, base_ring, names, sparse)
+
+    def __init__(self, base_ring, names, sparse):
+        self._arity = len(names)
+        super().__init__(PolynomialSpecies(base_ring, names))

From 041dcad01f48961635960e8fed1932f1493877cb Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Thu, 29 Aug 2024 13:22:52 +0200
Subject: [PATCH 03/31] lazy composition possibly working, at least univariate

---
 src/sage/rings/lazy_species.py | 152 ++++++++++++++++++++++++---------
 src/sage/rings/species.py      |  49 ++++++++++-
 2 files changed, 160 insertions(+), 41 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 79378278ab4..f056f1e425f 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -11,11 +11,11 @@
 from sage.libs.gap.libgap import libgap
 from sage.categories.sets_cat import cartesian_product
 from sage.combinat.set_partition import SetPartitions
+from sage.combinat.integer_vector import IntegerVectors
 from sage.structure.element import parent
 import itertools
-
-
-def weighted_compositions(n, d, weights):
+from collections import defaultdict
+def weighted_compositions(n, d, weights, offset=0):
     r"""
     Return all compositions of `n` of weight `d`.
 
@@ -41,8 +41,6 @@ def weighted_compositions(n, d, weights):
         [[0, 2, 0, 1], [0, 2, 1], [1, 1, 0, 1], [1, 1, 1], [2, 0, 0, 1], [2, 0, 1]]
 
     """
-    n = int(n)
-    d = int(d)
     # the empty composition exists if and only if n == d == 0
     if not n:
         if not d:
@@ -52,14 +50,49 @@ def weighted_compositions(n, d, weights):
         return
 
     # otherwise we iterate over the possibilities for the first part
-    w0 = weights[0]
+    if offset < len(weights):
+        w0 = weights[offset]
+    else:
+        return
     if w0 > d:
         return
-    wr = weights[1:]
     for i in range(min(n, d // w0) + 1):
-        for c in weighted_compositions(n - i, d - i * w0, wr):
+        for c in weighted_compositions(n - i, d - i * w0, weights, offset=offset+1):
             yield [i] + c
 
+def weighted_vector_compositions(n_vec, d, weights_vec):
+    r"""
+    Return all compositions of the vector `n` of weight `d`.
+
+    INPUT:
+
+    - ``n_vec``, a `k`-tuple of non-negative integers.
+
+    - ``d``, a non-negative integer.
+
+    - ``weights_vec``, `k`-tuple of weakly increasing lists of
+      positive integers.
+
+    EXAMPLES::
+
+        sage: list(weighted_vector_compositions([1,1], 2, [[1,1,2,3], [1,2,3]]))
+        [([0, 1], [1]), ([1], [1])]
+
+        sage: list(weighted_vector_compositions([3,1], 4, [[1,1,2,5], [1,1,2,5]]))
+        [([0, 3], [0, 1]),
+         ([0, 3], [1]),
+         ([1, 2], [0, 1]),
+         ([1, 2], [1]),
+         ([2, 1], [0, 1]),
+         ([2, 1], [1]),
+         ([3], [0, 1]),
+         ([3], [1])]
+
+    """
+    k = len(n_vec)
+    for d_vec in IntegerVectors(d, length=k):
+        yield from itertools.product(*map(weighted_compositions, n_vec, d_vec, weights_vec))
+
 
 ######################################################################
 
@@ -258,62 +291,105 @@ def polynomial(self, degree=None, names=None):
 
         return R.sum(self[:m])
 
-    def __call__(self, *g):
+    def __call__(self, *args):
         """
         EXAMPLES::
 
             sage: from sage.rings.lazy_species import LazySpecies
-            sage: L = LazySpecies(ZZ, "X")
-            sage: L(SymmetricGroup(2))(L(SymmetricGroup(2)))
-            P_4
+            sage: L = LazySpecies(QQ, "X")
+            sage: E2 = L(SymmetricGroup(2))
+            sage: E2(E2)
+            P_4 + O^11
+
+            sage: P = PolynomialSpecies(QQ, "X")
+            sage: Gc = L(lambda n: sum(P(G.automorphism_group()) for G in graphs(n) if G.is_connected()) if n else 0)
+            sage: E = L(lambda n: SymmetricGroup(n))
+
+            sage: G = L(lambda n: sum(P(G.automorphism_group()) for G in graphs(n)))
+            sage: (E-1)(Gc) - G
+            (-1) + O^7
+
+
+            sage: A = L.undefined(1)
+            sage: X = L(SymmetricGroup(1))
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: A.define(X*(E-1)(A) + X)
+            sage: A
+
+        TESTS::
+
+            sage: X = L(SymmetricGroup(1))
+            sage: X(X + E2)
+            X + E_2 + O^8
+            sage: E2(X + E2)
+            E_2 + X*E_2 + P_4 + O^9
+
+            sage: (1+E2)(E2)
         """
         fP = parent(self)
-        if len(g) != fP._arity:
+        if len(args) != fP._arity:
             raise ValueError("arity of must be equal to the number of arguments provided")
 
         # Find a good parent for the result
         from sage.structure.element import get_coercion_model
         cm = get_coercion_model()
-        P = cm.common_parent(self.base_ring(), *[parent(h) for h in g])
+        P = cm.common_parent(self.base_ring(), *[parent(g) for g in args])
 
         # f = 0
         if isinstance(self._coeff_stream, Stream_zero):
             return P.zero()
 
-        # g = (0, ..., 0)
-        if all((not isinstance(h, LazyModuleElement) and not h)
-               or (isinstance(h, LazyModuleElement)
-                   and isinstance(h._coeff_stream, Stream_zero))
-               for h in g):
+        # args = (0, ..., 0)
+        if all((not isinstance(g, LazyModuleElement) and not g)
+               or (isinstance(g, LazyModuleElement)
+                   and isinstance(g._coeff_stream, Stream_zero))
+               for g in args):
             return P(self[0])
 
-        # f has finite length and f != 0
+        # f is a constant polynomial
         if (isinstance(self._coeff_stream, Stream_exact)
-            and not self._coeff_stream._constant):
-            # constant polynomial
-            poly = self.polynomial()
-            if poly.is_constant():
-                return P(poly)
-            return P(poly(g))
-
-        g = [P(h) for h in g]
-        R = P._internal_poly_ring.base_ring()
+            and not self._coeff_stream._constant
+            and self.polynomial().is_constant()):
+            return P(self.polynomial())
 
-        for h in g:
-            if h._coeff_stream._approximate_order == 0:
-                if not h._coeff_stream.is_uninitialized() and h[0]:
+        args = [P(g) for g in args]
+
+        for g in args:
+            if g._coeff_stream._approximate_order == 0:
+                if not g._coeff_stream.is_uninitialized() and g[0]:
                     raise ValueError("can only compose with a positive valuation series")
-                h._coeff_stream._approximate_order = 1
+                g._coeff_stream._approximate_order = 1
 
         sorder = self._coeff_stream._approximate_order
-        gv = min(h._coeff_stream._approximate_order for h in g)
+        gv = min(g._coeff_stream._approximate_order for g in args)
+        R = P._internal_poly_ring.base_ring()
 
         def coefficient(n):
-            r = R.zero()
+            args_flat = [[(M, c) for i in range(n+1) for M, c in g[i]]
+                         for g in args]
+            weights = [[M._tc for i in range(n+1) for M, _ in g[i]]
+                       for g in args]
+            result = R.zero()
             for i in range(n // gv + 1):
-                # Make sure the element returned from the composition is in P
-                r += P(self[i](g))[n]
-            return r
+                # compute homogeneous components
+                lF = defaultdict(R)
+                for M, c in self[i]:
+                    lF[M._mc] += R._from_dict({M: c})
+                for mc, F in lF.items():
+                    for degrees in weighted_vector_compositions(mc, n, weights):
+                        multiplicities = [c for alpha, g_flat in zip(degrees, args_flat)
+                                          for d, (_, c) in zip(alpha, g_flat) if d]
+                        molecules = [M for alpha, g_flat in zip(degrees, args_flat)
+                                          for d, (M, _) in zip(alpha, g_flat) if d]
+                        non_zero_degrees = [[d for d in alpha if d] for alpha in degrees]
+                        names = ["X%s" % i for i in range(len(molecules))]
+                        FX = F._compose_with_weighted_singletons(names,
+                                                                 multiplicities,
+                                                                 non_zero_degrees)
+                        FG = [(M(*molecules), c) for M, c in FX]
+                        result += R.sum_of_terms(FG)
+            return result
+
         coeff_stream = Stream_function(coefficient, P._sparse, sorder * gv)
         return P.element_class(P, coeff_stream)
 
diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index dcd7ac4460d..866a60059ab 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -1469,10 +1469,29 @@ def exponential(self, multiplicities, degrees):
 
         TESTS::
 
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(QQ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: P = PolynomialSpecies(QQ, ["X"])
+
+            sage: c = 3/2; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            True
+
+            sage: c = -5/3; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            True
+
+            sage: c = 0; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            True
+
+            sage: c = 1; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            True
+
+            sage: c = -1; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            True
+
             sage: P = PolynomialSpecies(QQ, ["X"])
             sage: P.exponential([1], [0]).parent()
             Polynomial species in X over Rational Field
-
         """
         def stretch(c, k):
             r"""
@@ -1513,8 +1532,32 @@ def is_constant(self):
                 True
                 sage: P(0).is_constant()
                 True
+                sage: (1 + X).is_constant()
+            """
+            return self.is_zero() or not self.maximal_degree()
+
+        def homogeneous_degree(self):
             """
-            return self.is_zero() or not self.degree()
+
+            ..TODO::
+
+               This implementation should not be necessary.
+
+            EXAMPLES::
+
+                sage: P = PolynomialSpecies(ZZ, ["X"])
+                sage: C3 = P(CyclicPermutationGroup(3))
+                sage: X = P(SymmetricGroup(1))
+                sage: E2 = P(SymmetricGroup(2))
+                sage: (E2*X + C3).homogeneous_degree()
+                3
+            """
+            if not self.support():
+                raise ValueError("the zero element does not have a well-defined degree")
+            if not self.is_homogeneous():
+                raise ValueError("element is not homogeneous")
+            return self.parent().degree_on_basis(self.support()[0])
+
 
         def is_virtual(self):
             r"""
@@ -1735,6 +1778,7 @@ def __call__(self, *args):
             multiplicities = list(chain.from_iterable([[c for _, c in g] for g in args]))
             molecules = list(chain.from_iterable([[M for M, _ in g] for g in args]))
             F_degrees = sorted(set(M._mc for M, _ in self))
+            names = ["X%s" % i for i in range(sum(len(arg) for arg in args))]
 
             result = P0.zero()
             for n in F_degrees:
@@ -1742,7 +1786,6 @@ def __call__(self, *args):
                 for degrees in cartesian_product([IntegerVectors(n_i, length=len(arg))
                                                   for n_i, arg in zip(n, args)]):
                     # each degree is a weak composition of the degree of F in sort i
-                    names = ["X%s" % i for i in range(sum(len(arg) for arg in args))]
                     FX = F._compose_with_weighted_singletons(names,
                                                              multiplicities,
                                                              degrees)

From a1d42ef63a976e94d8f06c75e187e14504b4fd9f Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Thu, 29 Aug 2024 17:37:57 +0200
Subject: [PATCH 04/31] fix multisort bugs

---
 src/sage/rings/lazy_species.py | 61 +++++++++++++++++++++++++---------
 src/sage/rings/species.py      | 44 +++++++++++++++---------
 2 files changed, 75 insertions(+), 30 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index f056f1e425f..55990f015f0 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -1,20 +1,17 @@
-from sage.rings.integer_ring import ZZ
 from sage.rings.lazy_series import LazyCompletionGradedAlgebraElement, LazyModuleElement
 from sage.rings.lazy_series_ring import LazyCompletionGradedAlgebra
-from sage.data_structures.stream import (
-    Stream_zero,
-    Stream_exact,
-    Stream_function,
-    Stream_cauchy_compose,
-)
-from sage.rings.species import AtomicSpecies, PolynomialSpecies
+from sage.data_structures.stream import (Stream_zero,
+                                         Stream_exact,
+                                         Stream_function)
+from sage.rings.species import PolynomialSpecies
 from sage.libs.gap.libgap import libgap
 from sage.categories.sets_cat import cartesian_product
-from sage.combinat.set_partition import SetPartitions
 from sage.combinat.integer_vector import IntegerVectors
 from sage.structure.element import parent
 import itertools
 from collections import defaultdict
+
+
 def weighted_compositions(n, d, weights, offset=0):
     r"""
     Return all compositions of `n` of weight `d`.
@@ -60,6 +57,7 @@ def weighted_compositions(n, d, weights, offset=0):
         for c in weighted_compositions(n - i, d - i * w0, weights, offset=offset+1):
             yield [i] + c
 
+
 def weighted_vector_compositions(n_vec, d, weights_vec):
     r"""
     Return all compositions of the vector `n` of weight `d`.
@@ -75,6 +73,7 @@ def weighted_vector_compositions(n_vec, d, weights_vec):
 
     EXAMPLES::
 
+        sage: from sage.rings.lazy_species import weighted_vector_compositions
         sage: list(weighted_vector_compositions([1,1], 2, [[1,1,2,3], [1,2,3]]))
         [([0, 1], [1]), ([1], [1])]
 
@@ -306,25 +305,53 @@ def __call__(self, *args):
             sage: E = L(lambda n: SymmetricGroup(n))
 
             sage: G = L(lambda n: sum(P(G.automorphism_group()) for G in graphs(n)))
-            sage: (E-1)(Gc) - G
-            (-1) + O^7
-
+            sage: E(Gc) - G
+            O^7
 
             sage: A = L.undefined(1)
             sage: X = L(SymmetricGroup(1))
             sage: E = L(lambda n: SymmetricGroup(n))
-            sage: A.define(X*(E-1)(A) + X)
+            sage: A.define(X*E(A))
             sage: A
+            X + X^2 + (X^3+X*E_2) + (X^2*E_2+2*X^4+X*E_3)
+             + (X^2*E_3+3*X^5+3*X^3*E_2+X*{((1,2)(3,4),):({1,2,3,4})}+X*E_4)
+             + (X^2*E_4+2*X^2*{((1,2)(3,4),):({1,2,3,4})}+6*X^4*E_2+6*X^6
+                +3*X^3*E_3+X^2*E_2^2+X*E_5)
+             + (X^2*E_5+2*X^3*E_2^2+6*X^4*E_3+12*X^7+14*X^5*E_2
+                +3*X^3*{((1,2)(3,4),):({1,2,3,4})}+3*X^3*E_4
+                +X*{((3,4),(1,2),(1,3)(2,4)(5,6)):({1,2,3,4,5,6})}
+                +X*{((1,2)(3,4)(5,6),):({1,2,3,4,5,6})}
+                +X*{((2,3)(4,5),(1,3)(5,6)):({1,2,3,4,5,6})}+2*X^2*E_2*E_3
+                +E_2*{((1,2)(3,4),):({1,2,3,4})}*X+X*E_6)
+             + O^8
 
         TESTS::
 
+            sage: L = LazySpecies(QQ, "X")
             sage: X = L(SymmetricGroup(1))
+            sage: E2 = L(SymmetricGroup(2))
             sage: X(X + E2)
             X + E_2 + O^8
             sage: E2(X + E2)
             E_2 + X*E_2 + P_4 + O^9
 
-            sage: (1+E2)(E2)
+            sage: (1+E2)(X)
+            1 + E_2 + O^7
+
+            sage: L = LazySpecies(QQ, "X, Y")
+            sage: P = PolynomialSpecies(QQ, "X, Y")
+            sage: X = L(P(SymmetricGroup(1), {1: [1]}))
+            sage: Y = L(P(SymmetricGroup(1), {2: [1]}))
+            sage: X(Y, 0)
+            Y + O^8
+
+            sage: L1 = LazySpecies(QQ, "X")
+            sage: L = LazySpecies(QQ, "X, Y")
+            sage: P = PolynomialSpecies(QQ, "X, Y")
+            sage: X = L(P(SymmetricGroup(1), {1:[1]}))
+            sage: E = L1(lambda n: SymmetricGroup(n))
+            sage: E(X)
+            1 + X + E_2(X) + E_3(X) + E_4(X) + E_5(X) + E_6(X) + O^7
         """
         fP = parent(self)
         if len(args) != fP._arity:
@@ -365,6 +392,10 @@ def __call__(self, *args):
         R = P._internal_poly_ring.base_ring()
 
         def coefficient(n):
+            if not n:
+                if self[0]:
+                    return R(list(self[0])[0][1])
+                return R.zero()
             args_flat = [[(M, c) for i in range(n+1) for M, c in g[i]]
                          for g in args]
             weights = [[M._tc for i in range(n+1) for M, _ in g[i]]
@@ -416,5 +447,5 @@ def __classcall_private__(cls, base_ring, names, sparse=True):
         return super().__classcall__(cls, base_ring, names, sparse)
 
     def __init__(self, base_ring, names, sparse):
-        self._arity = len(names)
         super().__init__(PolynomialSpecies(base_ring, names))
+        self._arity = len(names)
diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index 866a60059ab..e0be0b6d7c4 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -1,6 +1,4 @@
 from itertools import accumulate, chain
-
-from sage.arith.misc import multinomial
 from sage.categories.graded_algebras_with_basis import GradedAlgebrasWithBasis
 from sage.categories.infinite_enumerated_sets import InfiniteEnumeratedSets
 from sage.categories.monoids import Monoids
@@ -14,11 +12,9 @@
 from sage.groups.perm_gps.permgroup_named import SymmetricGroup
 from sage.libs.gap.libgap import libgap
 from sage.misc.cachefunc import cached_method
-from sage.misc.misc_c import prod
 from sage.monoids.indexed_free_monoid import (IndexedFreeAbelianMonoid,
                                               IndexedFreeAbelianMonoidElement,
                                               IndexedMonoid)
-from sage.functions.other import binomial
 from sage.rings.integer_ring import ZZ
 from sage.sets.finite_enumerated_set import FiniteEnumeratedSet
 from sage.structure.element import Element, parent
@@ -1166,8 +1162,14 @@ def _compose_with_singletons(self, base_ring, names, args):
 
         def __call__(self, *args):
             r"""
-            Substitute M_1...M_k into self.
-            M_i must all have same arity and must be molecular.
+            Substitute `M_1,\dots, M_k` into ``self``.
+
+            The arguments must all have the same parent and must all
+            be molecular.  The number of arguments must be equal to
+            the arity of ``self``.
+
+            The result is a molecular species, whose parent is the
+            same as those of the arguments.
 
             EXAMPLES::
 
@@ -1198,13 +1200,21 @@ def __call__(self, *args):
                 sage: Y = M2(SymmetricGroup(1), {2: [1]})
                 sage: C3(X*Y)
                 {((1,2,3)(4,5,6),): ({1, 2, 3}, {4, 5, 6})}
-            """
+
+            TESTS::
+
+                sage: P = PolynomialSpecies(QQ, ["X"])
+                sage: P.one()()
+                Traceback (most recent call last):
+                ...
+                ValueError: number of args must match arity of self
+                """
             if len(args) != self.parent()._k:
                 raise ValueError("number of args must match arity of self")
+            if len(set(arg.parent() for arg in args)) > 1:
+                raise ValueError("all args must have the same parent")
             if not all(isinstance(arg, MolecularSpecies.Element) for arg in args):
                 raise ValueError("all args must be molecular species")
-            if len(set(arg.parent()._k for arg in args)) > 1:
-                raise ValueError("all args must have same arity")
 
             gens = []
 
@@ -1533,6 +1543,7 @@ def is_constant(self):
                 sage: P(0).is_constant()
                 True
                 sage: (1 + X).is_constant()
+                False
             """
             return self.is_zero() or not self.maximal_degree()
 
@@ -1558,7 +1569,6 @@ def homogeneous_degree(self):
                 raise ValueError("element is not homogeneous")
             return self.parent().degree_on_basis(self.support()[0])
 
-
         def is_virtual(self):
             r"""
             Return if ``self`` is a virtual species.
@@ -1769,12 +1779,16 @@ def __call__(self, *args):
 
             """
             P = self.parent()
-            if not self.support():
-                return P.zero()
+            if len(args) != P._k:
+                raise ValueError("number of args must match arity of self")
+            if len(set(arg.parent() for arg in args)) > 1:
+                raise ValueError("all args must have the same parent")
+
             P0 = args[0].parent()
-            assert all(P0 == arg.parent() for arg in args), "all parents must be the same"
-            args = [sorted(g, key=lambda x: x[0]._mc)
-                    for g in args]
+            if not self.support():
+                return P0.zero()
+
+            args = [sorted(g, key=lambda x: x[0]._mc) for g in args]
             multiplicities = list(chain.from_iterable([[c for _, c in g] for g in args]))
             molecules = list(chain.from_iterable([[M for M, _ in g] for g in args]))
             F_degrees = sorted(set(M._mc for M, _ in self))

From 049642205bb7705ddac05731aace52e443dfbf22 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Thu, 29 Aug 2024 18:01:41 +0200
Subject: [PATCH 05/31] add weighted example

---
 src/sage/rings/lazy_species.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 55990f015f0..dddc146c7aa 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -325,6 +325,22 @@ def __call__(self, *args):
                 +E_2*{((1,2)(3,4),):({1,2,3,4})}*X+X*E_6)
              + O^8
 
+            sage: R.<q> = QQ[]
+            sage: L = LazySpecies(R, "X")
+            sage: P = PolynomialSpecies(QQ, "X, Y")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: E1 = L(lambda n: SymmetricGroup(n) if n else 0)
+            sage: E(q*E1)
+            1 + q*X + ((q^2+q)*E_2) + ((q^3+q)*E_3+q^2*X*E_2)
+             + ((q^4+q)*E_4+q^2*P_4+q^2*X*E_3+q^3*E_2^2)
+             + ((q^5+q)*E_5+(q^4+q^3+q^2)*E_3*E_2+q^2*X*E_4+q^3*P_4*X)
+             + ((q^6+q)*E_6+q^2*{((1,2,3)(4,6),(1,4)(2,5)(3,6)):({1,2,3,4,5,6})}
+                +(q^5+q^3+q^2)*E_4*E_2+q^2*X*E_5
+                +q^3*{((1,4)(2,3),(1,6,3,2,5,4)):({1,2,3,4,5,6})}
+                +q^3*E_3*E_2*X+q^4*E_2*P_4+q^4*E_3^2)
+             + O^7
+
+
         TESTS::
 
             sage: L = LazySpecies(QQ, "X")

From 3c423645a33d622c7fd0685c4980ac03825b1cb2 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Fri, 30 Aug 2024 15:53:58 +0200
Subject: [PATCH 06/31] simplify unisort _repr_

---
 src/sage/rings/species.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index e0be0b6d7c4..d574e2823c7 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -467,6 +467,8 @@ def _repr_(self):
             sage: A = At(G, {2: [1,2,3,4,5,6,7,8,9,10]}); A
             {((1,2,3,4)(5,6)(7,8)(9,10),): ({}, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10})}
         """
+        if self.parent()._k == 1:
+            return "{" + f"{self._dis}" + "}"
         dompart = ', '.join("{" + repr(sorted(b))[1:-1] + "}"
                            for b in self._dompart)
         return "{" + f"{self._dis}: ({dompart})" + "}"

From f07046314753c2c3d08dfaf75b13841216da88eb Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Fri, 30 Aug 2024 15:54:20 +0200
Subject: [PATCH 07/31] polish examples

---
 src/sage/rings/lazy_species.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index dddc146c7aa..8c3fce740ce 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -308,9 +308,10 @@ def __call__(self, *args):
             sage: E(Gc) - G
             O^7
 
-            sage: A = L.undefined(1)
+            sage: L = LazySpecies(QQ, "X")
             sage: X = L(SymmetricGroup(1))
             sage: E = L(lambda n: SymmetricGroup(n))
+            sage: A = L.undefined(1)
             sage: A.define(X*E(A))
             sage: A
             X + X^2 + (X^3+X*E_2) + (X^2*E_2+2*X^4+X*E_3)
@@ -325,9 +326,16 @@ def __call__(self, *args):
                 +E_2*{((1,2)(3,4),):({1,2,3,4})}*X+X*E_6)
              + O^8
 
+            sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
+            sage: F = E(C(A))
+            sage: [sum(F[n].monomial_coefficients().values()) for n in range(1, 7)]
+            [1, 3, 7, 19, 47, 130]
+            sage: oeis(_)
+            0: A001372: Number of unlabeled mappings (or mapping patterns) from n points to themselves; number of unlabeled endofunctions.
+
+
             sage: R.<q> = QQ[]
             sage: L = LazySpecies(R, "X")
-            sage: P = PolynomialSpecies(QQ, "X, Y")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: E1 = L(lambda n: SymmetricGroup(n) if n else 0)
             sage: E(q*E1)

From 259f844d227fda7fcf435cd654c8dd40aa68d297 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Fri, 30 Aug 2024 20:52:55 +0200
Subject: [PATCH 08/31] add a doctest

---
 src/sage/rings/lazy_species.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 0c62a366c1b..0123498a59c 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -470,5 +470,17 @@ def __classcall_private__(cls, base_ring, names, sparse=True):
         return super().__classcall__(cls, base_ring, names, sparse)
 
     def __init__(self, base_ring, names, sparse):
+        """
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: LazySpecies(QQ, "X, Y")
+            Lazy completion of Polynomial species in X, Y over Rational Field
+
+        TESTS::
+
+            sage: LazySpecies(QQ, "X, Y, Z")._arity
+            3
+        """
         super().__init__(PolynomialSpecies(base_ring, names))
         self._arity = len(names)

From cc32428cd671f4168e3584fe194acdda7c6a11b3 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Fri, 30 Aug 2024 20:53:00 +0200
Subject: [PATCH 09/31] add a doctest

---
 src/sage/rings/species.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index 12de97ab907..f7dec40b26c 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -845,6 +845,14 @@ def gen(self, x):
             sage: type(m)
             <class 'sage.rings.species.MolecularSpecies_with_category.element_class'>
 
+            sage: P = PolynomialSpecies(ZZ, ["X"])
+            sage: M = P._indices
+            sage: A = AtomicSpecies("X")
+            sage: a = A(CyclicPermutationGroup(4))
+            sage: M.gen(a)
+            C_4
+
+            sage: P = PolynomialSpecies(ZZ, ["X", "Y"])
             sage: M = P._indices
             sage: m = M(SymmetricGroup(6).young_subgroup([2, 2, 2]), {1: [1,2], 2: [3,4,5,6]})
             sage: list(m)

From 3c2a3d13ecc4354b77fbc4c877745fd7522944dc Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sat, 31 Aug 2024 20:36:35 +0200
Subject: [PATCH 10/31] implement generating series, delegate sum and product
 to separate classes

---
 src/sage/rings/lazy_species.py | 178 ++++++++++++++++++++++++++++-----
 1 file changed, 153 insertions(+), 25 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 0123498a59c..aa1cdd595eb 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -10,7 +10,8 @@
 from sage.structure.element import parent
 import itertools
 from collections import defaultdict
-
+from sage.rings.lazy_series_ring import LazyPowerSeriesRing, LazySymmetricFunctions
+from sage.combinat.sf.sf import SymmetricFunctions
 
 def weighted_compositions(n, d, weights, offset=0):
     r"""
@@ -97,6 +98,9 @@ def weighted_vector_compositions(n_vec, d, weights_vec):
 
 class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
     r"""
+
+    EXAMPLES:
+
     Compute the molecular expansion of `E(-X)`::
 
         sage: from sage.rings.lazy_species import LazySpecies
@@ -117,14 +121,129 @@ class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
         sage: all(coefficient(m) == (1/E)[m] for m in range(10))
         True
     """
+    def isotype_generating_series(self):
+        r"""
+        Return the isotype generating series of ``self``.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(QQ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: E.isotype_generating_series()
+            1 + X + X^2 + X^3 + X^4 + X^5 + X^6 + O(X^7)
+
+            sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
+            sage: E(C).isotype_generating_series()
+            1 + X + 2*X^2 + 3*X^3 + 5*X^4 + 7*X^5 + 11*X^6 + O(X^7)
+
+            sage: L2 = LazySpecies(QQ, "X, Y")
+            sage: P2 = PolynomialSpecies(QQ, "X, Y")
+            sage: X = L2(P2(SymmetricGroup(1), {1: [1]}))
+            sage: Y = L2(P2(SymmetricGroup(1), {2: [1]}))
+            sage: E(X + Y).isotype_generating_series()
+            1 + (X+Y) + (X^2+X*Y+Y^2) + (X^3+X^2*Y+X*Y^2+Y^3)
+            + (X^4+X^3*Y+X^2*Y^2+X*Y^3+Y^4)
+            + (X^5+X^4*Y+X^3*Y^2+X^2*Y^3+X*Y^4+Y^5)
+            + (X^6+X^5*Y+X^4*Y^2+X^3*Y^3+X^2*Y^4+X*Y^5+Y^6)
+            + O(X,Y)^7
+
+            sage: C(X + Y).isotype_generating_series()
+            (X+Y) + (X^2+X*Y+Y^2) + (X^3+X^2*Y+X*Y^2+Y^3)
+            + (X^4+X^3*Y+2*X^2*Y^2+X*Y^3+Y^4)
+            + (X^5+X^4*Y+2*X^3*Y^2+2*X^2*Y^3+X*Y^4+Y^5)
+            + (X^6+X^5*Y+3*X^4*Y^2+4*X^3*Y^3+3*X^2*Y^4+X*Y^5+Y^6)
+            + O(X,Y)^7
+        """
+        P = self.parent()
+        L = LazyPowerSeriesRing(P.base_ring().fraction_field(),
+                                P._laurent_poly_ring._indices._indices.variable_names())
+        if P._arity == 1:
+            def coefficient(n):
+                return sum(c for M, c in self[n].monomial_coefficients().items())
+        else:
+            def coefficient(n):
+                return sum(c * P.base_ring().prod(v ** d for v, d in zip(L.gens(), M.grade()))
+                           for M, c in self[n].monomial_coefficients().items())
+        return L(coefficient)
+
     def generating_series(self):
-        pass
+        r"""
+        Return the (exponential) generating series of ``self``.
 
-    def isotype_generating_series(self):
-        pass
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: E.generating_series()
+            1 + X + 1/2*X^2 + 1/6*X^3 + 1/24*X^4 + 1/120*X^5 + 1/720*X^6 + O(X^7)
+
+            sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
+            sage: C.generating_series()
+            X + 1/2*X^2 + 1/3*X^3 + 1/4*X^4 + 1/5*X^5 + 1/6*X^6 + O(X^7)
+
+            sage: L2 = LazySpecies(QQ, "X, Y")
+            sage: P2 = PolynomialSpecies(QQ, "X, Y")
+            sage: X = L2(P2(SymmetricGroup(1), {1: [1]}))
+            sage: Y = L2(P2(SymmetricGroup(1), {2: [1]}))
+            sage: E(X + Y).generating_series()
+            1 + (X+Y) + (1/2*X^2+X*Y+1/2*Y^2)
+            + (1/6*X^3+1/2*X^2*Y+1/2*X*Y^2+1/6*Y^3)
+            + (1/24*X^4+1/6*X^3*Y+1/4*X^2*Y^2+1/6*X*Y^3+1/24*Y^4)
+            + (1/120*X^5+1/24*X^4*Y+1/12*X^3*Y^2+1/12*X^2*Y^3+1/24*X*Y^4+1/120*Y^5)
+            + (1/720*X^6+1/120*X^5*Y+1/48*X^4*Y^2+1/36*X^3*Y^3+1/48*X^2*Y^4+1/120*X*Y^5+1/720*Y^6)
+            + O(X,Y)^7
+
+            sage: C(X + Y).generating_series()
+            (X+Y) + (1/2*X^2+X*Y+1/2*Y^2) + (1/3*X^3+X^2*Y+X*Y^2+1/3*Y^3)
+            + (1/4*X^4+X^3*Y+3/2*X^2*Y^2+X*Y^3+1/4*Y^4)
+            + (1/5*X^5+X^4*Y+2*X^3*Y^2+2*X^2*Y^3+X*Y^4+1/5*Y^5)
+            + (1/6*X^6+X^5*Y+3*X^4*Y^2+4*X^3*Y^3+3*X^2*Y^4+X*Y^5+1/6*Y^6)
+            + O(X,Y)^7
+        """
+        P = self.parent()
+        L = LazyPowerSeriesRing(P.base_ring().fraction_field(),
+                                P._laurent_poly_ring._indices._indices.variable_names())
+        if P._arity == 1:
+            def coefficient(n):
+                return sum(c / M._group.cardinality()
+                           for M, c in self[n].monomial_coefficients().items())
+        else:
+            def coefficient(n):
+                return sum(c / M._group.cardinality()
+                           * P.base_ring().prod(v ** d for v, d in zip(L.gens(), M.grade()))
+                           for M, c in self[n].monomial_coefficients().items())
+        return L(coefficient)
 
     def cycle_index_series(self):
-        pass
+        r"""
+        Return the cycle index series for this species.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: h = SymmetricFunctions(QQ).h()
+            sage: LazySymmetricFunctions(h)(E.cycle_index_series())
+            h[] + h[1] + h[2] + h[3] + h[4] + h[5] + h[6] + O^7
+
+            sage: s = SymmetricFunctions(QQ).s()
+            sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
+            sage: s(C.cycle_index_series()[5])
+            s[1, 1, 1, 1, 1] + s[2, 2, 1] + 2*s[3, 1, 1] + s[3, 2] + s[5]
+        """
+        P = self.parent()
+        p = SymmetricFunctions(P.base_ring().fraction_field()).p()
+        if P._arity == 1:
+            L = LazySymmetricFunctions(p)
+            def coefficient(n):
+                return sum(c * M._group.cycle_index()
+                           for M, c in self[n].monomial_coefficients().items())
+        else:
+            raise NotImplementedError
+        return L(coefficient)
 
     def _add_(self, other):
         r"""
@@ -145,13 +264,7 @@ def _add_(self, other):
             [((1, 2, 3), E_3), ((1, 2, 3), E_3)]
 
         """
-        def add_structures(labels):
-            yield from self.structures(labels)
-            yield from other.structures(labels)
-
-        result = super()._add_(other)
-        result.structures = add_structures
-        return result
+        return SumSpeciesElement(self, other)
 
     def _mul_(self, other):
         """
@@ -174,19 +287,7 @@ def _mul_(self, other):
              (((1, 3), E_2), ((2,), X)),
              (((2, 3), E_2), ((1,), X))]
         """
-        def mul_structures(labels):
-            n = len(labels)
-            l = set(labels)
-            assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
-            for k in range(n):
-                for U in itertools.combinations(l, k):
-                    V = l.difference(U)
-                    yield from itertools.product(self.structures(U),
-                                                 other.structures(V))
-
-        result = super()._mul_(other)
-        result.structures = mul_structures
-        return result
+        return ProductSpeciesElement(self, other)
 
     def structures(self, labels):
         r"""
@@ -447,6 +548,33 @@ def coefficient(n):
         coeff_stream = Stream_function(coefficient, P._sparse, sorder * gv)
         return P.element_class(P, coeff_stream)
 
+class SumSpeciesElement(LazySpeciesElement):
+    def __init__(self, left, right):
+        self._left = left
+        self._right = right
+        add = super(LazySpeciesElement, type(left))._add_(left, right)
+        super().__init__(add.parent(), add._coeff_stream)
+
+    def structures(self, labels):
+        yield from self._left.structures(labels)
+        yield from self._right.structures(labels)
+
+class ProductSpeciesElement(LazySpeciesElement):
+    def __init__(self, left, right):
+        self._left = left
+        self._right = right
+        add = super(LazySpeciesElement, type(left))._mul_(left, right)
+        super().__init__(add.parent(), add._coeff_stream)
+
+    def structures(self, labels):
+        n = len(labels)
+        l = set(labels)
+        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
+        for k in range(n):
+            for U in itertools.combinations(l, k):
+                V = l.difference(U)
+                yield from itertools.product(self._left.structures(U),
+                                             self._right.structures(V))
 
 class LazySpecies(LazyCompletionGradedAlgebra):
     """

From 3ea68d1432ecc015e9c9478cac1e7cdb66575ef2 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Tue, 3 Sep 2024 08:38:22 +0200
Subject: [PATCH 11/31] start to fix multisort structures

---
 src/sage/rings/lazy_species.py | 55 +++++++++++++++++++++-------------
 1 file changed, 35 insertions(+), 20 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index aa1cdd595eb..65e67b34ccc 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -12,6 +12,7 @@
 from collections import defaultdict
 from sage.rings.lazy_series_ring import LazyPowerSeriesRing, LazySymmetricFunctions
 from sage.combinat.sf.sf import SymmetricFunctions
+from sage.groups.perm_gps.permgroup_named import SymmetricGroup
 
 def weighted_compositions(n, d, weights, offset=0):
     r"""
@@ -289,7 +290,7 @@ def _mul_(self, other):
         """
         return ProductSpeciesElement(self, other)
 
-    def structures(self, labels):
+    def structures(self, *labels):
         r"""
         Iterate over the structures on the given set of labels.
 
@@ -326,14 +327,28 @@ def structures(self, labels):
              ((2, 1, 3), X^3),
              ((2, 3, 1), X^3),
              ((1, 3, 2), X^3)]
+
+            sage: L2 = LazySpecies(QQ, "X, Y")
+            sage: P2 = PolynomialSpecies(QQ, "X, Y")
+            sage: X = L2(P2(SymmetricGroup(1), {1: [1]}))
+            sage: Y = L2(P2(SymmetricGroup(1), {2: [1]}))
+            sage: list((X*Y).structures([1],[2]))
+            [((1, 2), X*Y)]
+
+            sage: list(E(X*Y).structures([1,2],[3,4]))
+            [((3, 4, 1, 2), {((1,2)(3,4),): ({1, 2}, {3, 4})}),
+             ((4, 3, 1, 2), {((1,2)(3,4),): ({1, 2}, {3, 4})})]
+
         """
-        from sage.groups.perm_gps.permgroup_named import SymmetricGroup
-        n = len(labels)
-        l = set(labels)
-        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
-        S = SymmetricGroup(n)
-        F = self[n]
-        l = list(l)[::-1]
+        fP = parent(self)
+        if len(labels) != fP._arity:
+            raise ValueError("arity of must be equal to the number of arguments provided")
+        assert all(len(U) == len(set(U)) for U in labels), f"The argument labels must be a set, but {labels} has duplicates"
+
+        n = tuple([len(U) for U in labels])
+        S = SymmetricGroup(sum(n)).young_subgroup(n)
+        F = self[sum(n)]
+        l = [e for l in labels for e in l][::-1]
         for M, c in F.monomial_coefficients().items():
             if c < 0:
                 raise NotImplementedError("only implemented for proper non-virtual species")
@@ -555,9 +570,9 @@ def __init__(self, left, right):
         add = super(LazySpeciesElement, type(left))._add_(left, right)
         super().__init__(add.parent(), add._coeff_stream)
 
-    def structures(self, labels):
-        yield from self._left.structures(labels)
-        yield from self._right.structures(labels)
+    def structures(self, *labels):
+        yield from self._left.structures(*labels)
+        yield from self._right.structures(*labels)
 
 class ProductSpeciesElement(LazySpeciesElement):
     def __init__(self, left, right):
@@ -566,15 +581,15 @@ def __init__(self, left, right):
         add = super(LazySpeciesElement, type(left))._mul_(left, right)
         super().__init__(add.parent(), add._coeff_stream)
 
-    def structures(self, labels):
-        n = len(labels)
-        l = set(labels)
-        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
-        for k in range(n):
-            for U in itertools.combinations(l, k):
-                V = l.difference(U)
-                yield from itertools.product(self._left.structures(U),
-                                             self._right.structures(V))
+#    def structures(self, labels):
+#        n = len(labels)
+#        l = set(labels)
+#        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
+#        for k in range(n):
+#            for U in itertools.combinations(l, k):
+#                V = l.difference(U)
+#                yield from itertools.product(self._left.structures(U),
+#                                             self._right.structures(V))
 
 class LazySpecies(LazyCompletionGradedAlgebra):
     """

From 7c8e8bf203c18f44b736a988401b3e5901e82b17 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Thu, 12 Sep 2024 08:25:43 +0200
Subject: [PATCH 12/31] alternative way to compute group

---
 src/sage/rings/species.py | 27 +++++++++++++++++++++++++--
 1 file changed, 25 insertions(+), 2 deletions(-)

diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index 6b228fb084a..619146df152 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -960,11 +960,34 @@ def group(self):
             EXAMPLES::
 
                 sage: from sage.rings.species import MolecularSpecies
-                sage: M = MolecularSpecies(AtomicSpecies("X,Y"), prefix='', bracket=False)
+                sage: M = MolecularSpecies("X,Y")
                 sage: G = PermutationGroup([[(1,2),(3,4)], [(5,6)]])
                 sage: F = M(G, {1: [5,6], 2: [1,2,3,4]})
+                sage: F.group()
+                Permutation Group with generators [(1,2)(3,4), (5,6)]
             """
-            pass
+            factors = list(self)
+            if not factors:
+                return SymmetricGroup(0)
+
+            if len(factors) == 1:
+                A, n = factors[0]
+                if n == 1:
+                    return list(A._monomial)[0]._dis._C
+
+                if n % 2 == 1:
+                    f_gap = libgap.DirectProduct(list(A._monomial)[0]._dis._C,
+                                                 (A ** (n-1)).group())
+                else:
+                    f1 = (A ** (n // 2)).group()
+                    f_gap = libgap.DirectProduct(f1, f1)
+
+                f = PermutationGroup(gap_group=f_gap)
+                return f
+
+            f_gap = libgap.DirectProduct(*[(A ** n).group() for A, n in factors])
+            f = PermutationGroup(gap_group=f_gap)
+            return f
 
         def _assign_group_info(self, other):
             r"""

From f7bfd48c4b4d479e5ec5bd5a8de8c21ca311b7ce Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sat, 14 Sep 2024 22:52:29 +0200
Subject: [PATCH 13/31] adapt doctests

---
 src/sage/rings/lazy_species.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index a28c3bc724e..c578276588e 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -140,8 +140,8 @@ def isotype_generating_series(self):
 
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L2(P2(SymmetricGroup(1), {1: [1]}))
-            sage: Y = L2(P2(SymmetricGroup(1), {2: [1]}))
+            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
+            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
             sage: E(X + Y).isotype_generating_series()
             1 + (X+Y) + (X^2+X*Y+Y^2) + (X^3+X^2*Y+X*Y^2+Y^3)
             + (X^4+X^3*Y+X^2*Y^2+X*Y^3+Y^4)
@@ -186,8 +186,8 @@ def generating_series(self):
 
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L2(P2(SymmetricGroup(1), {1: [1]}))
-            sage: Y = L2(P2(SymmetricGroup(1), {2: [1]}))
+            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
+            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
             sage: E(X + Y).generating_series()
             1 + (X+Y) + (1/2*X^2+X*Y+1/2*Y^2)
             + (1/6*X^3+1/2*X^2*Y+1/2*X*Y^2+1/6*Y^3)
@@ -330,8 +330,8 @@ def structures(self, *labels):
 
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L2(P2(SymmetricGroup(1), {1: [1]}))
-            sage: Y = L2(P2(SymmetricGroup(1), {2: [1]}))
+            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
+            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
             sage: list((X*Y).structures([1],[2]))
             [((1, 2), X*Y)]
 
@@ -479,15 +479,15 @@ def __call__(self, *args):
 
             sage: L = LazySpecies(QQ, "X, Y")
             sage: P = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L(P(SymmetricGroup(1), {1: [1]}))
-            sage: Y = L(P(SymmetricGroup(1), {2: [1]}))
+            sage: X = L(P(SymmetricGroup(1), {0: [1]}))
+            sage: Y = L(P(SymmetricGroup(1), {1: [1]}))
             sage: X(Y, 0)
             Y + O^8
 
             sage: L1 = LazySpecies(QQ, "X")
             sage: L = LazySpecies(QQ, "X, Y")
             sage: P = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L(P(SymmetricGroup(1), {1:[1]}))
+            sage: X = L(P(SymmetricGroup(1), {0: [1]}))
             sage: E = L1(lambda n: SymmetricGroup(n))
             sage: E(X)
             1 + X + E_2(X) + E_3(X) + E_4(X) + E_5(X) + E_6(X) + O^7

From 48fb1168c788b750fc208b9d81b886c57e1d25ad Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Mon, 23 Sep 2024 12:14:43 +0200
Subject: [PATCH 14/31] generic cycle_index_series for multisort species

---
 src/sage/rings/lazy_species.py | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index c578276588e..2445b1d28b4 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -6,6 +6,7 @@
 from sage.rings.species import PolynomialSpecies
 from sage.libs.gap.libgap import libgap
 from sage.categories.sets_cat import cartesian_product
+from sage.categories.tensor import tensor
 from sage.combinat.integer_vector import IntegerVectors
 from sage.structure.element import parent
 import itertools
@@ -234,6 +235,18 @@ def cycle_index_series(self):
             sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
             sage: s(C.cycle_index_series()[5])
             s[1, 1, 1, 1, 1] + s[2, 2, 1] + 2*s[3, 1, 1] + s[3, 2] + s[5]
+
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: L2 = LazySpecies(QQ, "X, Y")
+            sage: P2 = PolynomialSpecies(QQ, "X, Y")
+            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
+            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
+            sage: E(X + Y).cycle_index_series()[3]
+            1/6*p[] # p[1, 1, 1] + 1/2*p[] # p[2, 1] + 1/3*p[] # p[3]
+            + 1/2*p[1] # p[1, 1] + 1/2*p[1] # p[2] + 1/2*p[1, 1] # p[1]
+            + 1/6*p[1, 1, 1] # p[] + 1/2*p[2] # p[1] + 1/2*p[2, 1] # p[]
+            + 1/3*p[3] # p[]
         """
         P = self.parent()
         p = SymmetricFunctions(P.base_ring().fraction_field()).p()
@@ -243,7 +256,10 @@ def coefficient(n):
                 return sum(c * M.group_and_partition()[0].cycle_index()
                            for M, c in self[n].monomial_coefficients().items())
         else:
-            raise NotImplementedError
+            L = LazySymmetricFunctions(tensor([p for _ in range(P._arity)]))
+            def coefficient(n):
+                return sum(c * M.cycle_index()
+                           for M, c in self[n].monomial_coefficients().items())
         return L(coefficient)
 
     def _add_(self, other):

From 7c0b2bb9211ef8b4925c39c135aff1700f767591 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Wed, 25 Sep 2024 23:14:13 +0200
Subject: [PATCH 15/31] fix a doctest

---
 src/sage/rings/lazy_species.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 2445b1d28b4..09910918b65 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -201,7 +201,7 @@ def generating_series(self):
             (X+Y) + (1/2*X^2+X*Y+1/2*Y^2) + (1/3*X^3+X^2*Y+X*Y^2+1/3*Y^3)
             + (1/4*X^4+X^3*Y+3/2*X^2*Y^2+X*Y^3+1/4*Y^4)
             + (1/5*X^5+X^4*Y+2*X^3*Y^2+2*X^2*Y^3+X*Y^4+1/5*Y^5)
-            + (1/6*X^6+X^5*Y+3*X^4*Y^2+4*X^3*Y^3+3*X^2*Y^4+X*Y^5+1/6*Y^6)
+            + (1/6*X^6+X^5*Y+5/2*X^4*Y^2+10/3*X^3*Y^3+5/2*X^2*Y^4+X*Y^5+1/6*Y^6)
             + O(X,Y)^7
         """
         P = self.parent()
@@ -448,14 +448,14 @@ def __call__(self, *args):
             sage: A
             X + X^2 + (X^3+X*E_2) + (X^2*E_2+2*X^4+X*E_3)
              + (X^2*E_3+3*X^5+3*X^3*E_2+X*{((1,2)(3,4),)}+X*E_4)
-             + (X^2*E_4+2*X^2*{((1,2)(3,4),))}+6*X^4*E_2+6*X^6
-                +3*X^3*E_3+X^2*E_2^2+X*E_5)
-             + (X^2*E_5+2*X^3*E_2^2+6*X^4*E_3+12*X^7+14*X^5*E_2
-                +3*X^3*{((1,2)(3,4),))}+3*X^3*E_4
-                +X*{((3,4),(1,2),(1,3)(2,4)(5,6)))}
+             + (X^2*E_4+2*X^2*{((1,2)(3,4),)}+6*X^4*E_2+6*X^6
+                +3*X^3*E_3+E_2^2*X^2+X*E_5)
+             + (X^2*E_5+2*E_2^2*X^3+6*X^4*E_3+12*X^7+14*X^5*E_2
+                +3*X^3*{((1,2)(3,4),)}+3*X^3*E_4
+                +X*{((3,4),(1,2),(1,3)(2,4)(5,6))}
                 +X*{((1,2)(3,4)(5,6),)}
-                +X*{((2,3)(4,5),(1,3)(5,6))}+2*X^2*E_2*E_3
-                +E_2*{((1,2)(3,4),)}*X+X*E_6)
+                +X*{((1,2)(4,5),(1,2,3)(4,5,6))}
+                +2*E_2*X^2*E_3+E_2*{((1,2)(3,4),)}*X+X*E_6)
              + O^8
 
             sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
@@ -473,11 +473,11 @@ def __call__(self, *args):
             sage: E(q*E1)
             1 + q*X + ((q^2+q)*E_2) + ((q^3+q)*E_3+q^2*X*E_2)
              + ((q^4+q)*E_4+q^2*P_4+q^2*X*E_3+q^3*E_2^2)
-             + ((q^5+q)*E_5+(q^4+q^3+q^2)*E_3*E_2+q^2*X*E_4+q^3*P_4*X)
+             + ((q^5+q)*E_5+(q^4+q^3+q^2)*E_2*E_3+q^2*X*E_4+q^3*X*P_4)
              + ((q^6+q)*E_6+q^2*{((1,2,3)(4,6),(1,4)(2,5)(3,6))}
-                +(q^5+q^3+q^2)*E_4*E_2+q^2*X*E_5
+                +(q^5+q^3+q^2)*E_2*E_4+q^2*X*E_5
                 +q^3*{((1,4)(2,3),(1,6,3,2,5,4))}
-                +q^3*E_3*E_2*X+q^4*E_2*P_4+q^4*E_3^2)
+                +q^3*X*E_2*E_3+q^4*E_2*P_4+q^4*E_3^2)
              + O^7
 
         TESTS::

From c9819823b416619e67a70015476a67e533b73780 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Thu, 24 Oct 2024 23:04:17 +0200
Subject: [PATCH 16/31] import PolynomialSpecies in doctests

---
 src/sage/rings/lazy_species.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 09910918b65..7ec71fd0b11 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -139,6 +139,7 @@ def isotype_generating_series(self):
             sage: E(C).isotype_generating_series()
             1 + X + 2*X^2 + 3*X^3 + 5*X^4 + 7*X^5 + 11*X^6 + O(X^7)
 
+            sage: from sage.rings.species import PolynomialSpecies
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
             sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
@@ -185,6 +186,7 @@ def generating_series(self):
             sage: C.generating_series()
             X + 1/2*X^2 + 1/3*X^3 + 1/4*X^4 + 1/5*X^5 + 1/6*X^6 + O(X^7)
 
+            sage: from sage.rings.species import PolynomialSpecies
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
             sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
@@ -236,6 +238,7 @@ def cycle_index_series(self):
             sage: s(C.cycle_index_series()[5])
             s[1, 1, 1, 1, 1] + s[2, 2, 1] + 2*s[3, 1, 1] + s[3, 2] + s[5]
 
+            sage: from sage.rings.species import PolynomialSpecies
             sage: L = LazySpecies(ZZ, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: L2 = LazySpecies(QQ, "X, Y")
@@ -344,6 +347,7 @@ def structures(self, *labels):
              ((2, 3, 1), X^3),
              ((1, 3, 2), X^3)]
 
+            sage: from sage.rings.species import PolynomialSpecies
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
             sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
@@ -432,6 +436,7 @@ def __call__(self, *args):
             sage: E2(E2)
             P_4 + O^11
 
+            sage: from sage.rings.species import PolynomialSpecies
             sage: P = PolynomialSpecies(QQ, "X")
             sage: Gc = L(lambda n: sum(P(G.automorphism_group()) for G in graphs(n) if G.is_connected()) if n else 0)
             sage: E = L(lambda n: SymmetricGroup(n))

From 47c78142d1a2357da77e3372cea59f973260c1ec Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Tue, 29 Oct 2024 20:16:50 +0100
Subject: [PATCH 17/31] adapt to changes in species.py

---
 src/sage/rings/lazy_species.py | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 7ec71fd0b11..1302ece1248 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -177,7 +177,7 @@ def generating_series(self):
         EXAMPLES::
 
             sage: from sage.rings.lazy_species import LazySpecies
-            sage: L = LazySpecies(ZZ, "X")
+            sage: L = LazySpecies(QQ, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: E.generating_series()
             1 + X + 1/2*X^2 + 1/6*X^3 + 1/24*X^4 + 1/120*X^5 + 1/720*X^6 + O(X^7)
@@ -211,11 +211,11 @@ def generating_series(self):
                                 P._laurent_poly_ring._indices._indices.variable_names())
         if P._arity == 1:
             def coefficient(n):
-                return sum(c / M.group_and_partition()[0].cardinality()
+                return sum(c / M.permutation_group()[0].cardinality()
                            for M, c in self[n].monomial_coefficients().items())
         else:
             def coefficient(n):
-                return sum(c / M.group_and_partition()[0].cardinality()
+                return sum(c / M.permutation_group()[0].cardinality()
                            * P.base_ring().prod(v ** d for v, d in zip(L.gens(), M.grade()))
                            for M, c in self[n].monomial_coefficients().items())
         return L(coefficient)
@@ -239,7 +239,7 @@ def cycle_index_series(self):
             s[1, 1, 1, 1, 1] + s[2, 2, 1] + 2*s[3, 1, 1] + s[3, 2] + s[5]
 
             sage: from sage.rings.species import PolynomialSpecies
-            sage: L = LazySpecies(ZZ, "X")
+            sage: L = LazySpecies(QQ, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
@@ -256,7 +256,7 @@ def cycle_index_series(self):
         if P._arity == 1:
             L = LazySymmetricFunctions(p)
             def coefficient(n):
-                return sum(c * M.group_and_partition()[0].cycle_index()
+                return sum(c * M.permutation_group()[0].cycle_index()
                            for M, c in self[n].monomial_coefficients().items())
         else:
             L = LazySymmetricFunctions(tensor([p for _ in range(P._arity)]))
@@ -322,7 +322,7 @@ def structures(self, *labels):
         EXAMPLES::
 
             sage: from sage.rings.lazy_species import LazySpecies
-            sage: L = LazySpecies(ZZ, "X")
+            sage: L = LazySpecies(QQ, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: list(E.structures([1,2,3]))
             [((1, 2, 3), E_3)]
@@ -373,7 +373,7 @@ def structures(self, *labels):
             if c < 0:
                 raise NotImplementedError("only implemented for proper non-virtual species")
             types = [tuple(S(rep)._act_on_list_on_position(l))[::-1]
-                     for rep in libgap.RightTransversal(S, M.group_and_partition()[0])]
+                     for rep in libgap.RightTransversal(S, M.permutation_group()[0])]
             if c == 1:
                 for s in types:
                     yield s, M
@@ -550,6 +550,7 @@ def __call__(self, *args):
         sorder = self._coeff_stream._approximate_order
         gv = min(g._coeff_stream._approximate_order for g in args)
         R = P._internal_poly_ring.base_ring()
+        L = fP._internal_poly_ring.base_ring()
 
         def coefficient(n):
             if not n:
@@ -563,9 +564,9 @@ def coefficient(n):
             result = R.zero()
             for i in range(n // gv + 1):
                 # compute homogeneous components
-                lF = defaultdict(R)
+                lF = defaultdict(L)
                 for M, c in self[i]:
-                    lF[M.grade()] += R._from_dict({M: c})
+                    lF[M.grade()] += L._from_dict({M: c})
                 for mc, F in lF.items():
                     for degrees in weighted_vector_compositions(mc, n, weights):
                         multiplicities = [c for alpha, g_flat in zip(degrees, args_flat)

From 1d91d1e498c4cc028aea7afe13ad9f7d9bd825e5 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Mon, 11 Nov 2024 19:39:03 +0100
Subject: [PATCH 18/31] extract __call__ into separate class, proof of concept
 for special species

---
 src/sage/rings/lazy_species.py | 144 ++++++++++++++++++++++++---------
 1 file changed, 104 insertions(+), 40 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 1302ece1248..65c3c140ab6 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -513,17 +513,51 @@ def __call__(self, *args):
             sage: E(X)
             1 + X + E_2(X) + E_3(X) + E_4(X) + E_5(X) + E_6(X) + O^7
         """
-        fP = parent(self)
+        return CompositionSpeciesElement(self, *args)
+
+
+class SumSpeciesElement(LazySpeciesElement):
+    def __init__(self, left, right):
+        self._left = left
+        self._right = right
+        F = super(LazySpeciesElement, type(left))._add_(left, right)
+        super().__init__(F.parent(), F._coeff_stream)
+
+    def structures(self, *labels):
+        yield from self._left.structures(*labels)
+        yield from self._right.structures(*labels)
+
+class ProductSpeciesElement(LazySpeciesElement):
+    def __init__(self, left, right):
+        self._left = left
+        self._right = right
+        F = super(LazySpeciesElement, type(left))._mul_(left, right)
+        super().__init__(F.parent(), F._coeff_stream)
+
+#    def structures(self, labels):
+#        n = len(labels)
+#        l = set(labels)
+#        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
+#        for k in range(n+1):
+#            for U in itertools.combinations(l, k):
+#                V = l.difference(U)
+#                yield from itertools.product(self._left.structures(U),
+#                                             self._right.structures(V))
+
+
+class CompositionSpeciesElement(LazySpeciesElement):
+    def __init__(self, left, *args):
+        fP = parent(left)
         if len(args) != fP._arity:
             raise ValueError("arity of must be equal to the number of arguments provided")
 
         # Find a good parent for the result
         from sage.structure.element import get_coercion_model
         cm = get_coercion_model()
-        P = cm.common_parent(self.base_ring(), *[parent(g) for g in args])
+        P = cm.common_parent(left.base_ring(), *[parent(g) for g in args])
 
         # f = 0
-        if isinstance(self._coeff_stream, Stream_zero):
+        if isinstance(left._coeff_stream, Stream_zero):
             return P.zero()
 
         # args = (0, ..., 0)
@@ -531,13 +565,13 @@ def __call__(self, *args):
                or (isinstance(g, LazyModuleElement)
                    and isinstance(g._coeff_stream, Stream_zero))
                for g in args):
-            return P(self[0])
+            return P(left[0])
 
         # f is a constant polynomial
-        if (isinstance(self._coeff_stream, Stream_exact)
-            and not self._coeff_stream._constant
-            and self.polynomial().is_constant()):
-            return P(self.polynomial())
+        if (isinstance(left._coeff_stream, Stream_exact)
+            and not left._coeff_stream._constant
+            and left.polynomial().is_constant()):
+            return P(left.polynomial())
 
         args = [P(g) for g in args]
 
@@ -547,15 +581,15 @@ def __call__(self, *args):
                     raise ValueError("can only compose with a positive valuation series")
                 g._coeff_stream._approximate_order = 1
 
-        sorder = self._coeff_stream._approximate_order
+        sorder = left._coeff_stream._approximate_order
         gv = min(g._coeff_stream._approximate_order for g in args)
         R = P._internal_poly_ring.base_ring()
         L = fP._internal_poly_ring.base_ring()
 
         def coefficient(n):
             if not n:
-                if self[0]:
-                    return R(list(self[0])[0][1])
+                if left[0]:
+                    return R(list(left[0])[0][1])
                 return R.zero()
             args_flat = [[(M, c) for i in range(n+1) for M, c in g[i]]
                          for g in args]
@@ -565,7 +599,7 @@ def coefficient(n):
             for i in range(n // gv + 1):
                 # compute homogeneous components
                 lF = defaultdict(L)
-                for M, c in self[i]:
+                for M, c in left[i]:
                     lF[M.grade()] += L._from_dict({M: c})
                 for mc, F in lF.items():
                     for degrees in weighted_vector_compositions(mc, n, weights):
@@ -583,35 +617,8 @@ def coefficient(n):
             return result
 
         coeff_stream = Stream_function(coefficient, P._sparse, sorder * gv)
-        return P.element_class(P, coeff_stream)
+        super().__init__(P, coeff_stream)
 
-class SumSpeciesElement(LazySpeciesElement):
-    def __init__(self, left, right):
-        self._left = left
-        self._right = right
-        add = super(LazySpeciesElement, type(left))._add_(left, right)
-        super().__init__(add.parent(), add._coeff_stream)
-
-    def structures(self, *labels):
-        yield from self._left.structures(*labels)
-        yield from self._right.structures(*labels)
-
-class ProductSpeciesElement(LazySpeciesElement):
-    def __init__(self, left, right):
-        self._left = left
-        self._right = right
-        add = super(LazySpeciesElement, type(left))._mul_(left, right)
-        super().__init__(add.parent(), add._coeff_stream)
-
-#    def structures(self, labels):
-#        n = len(labels)
-#        l = set(labels)
-#        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
-#        for k in range(n):
-#            for U in itertools.combinations(l, k):
-#                V = l.difference(U)
-#                yield from itertools.product(self._left.structures(U),
-#                                             self._right.structures(V))
 
 class LazySpecies(LazyCompletionGradedAlgebra):
     """
@@ -642,6 +649,12 @@ def __init__(self, base_ring, names, sparse):
             sage: LazySpecies(QQ, "X, Y")
             Lazy completion of Polynomial species in X, Y over Rational Field
 
+            sage: L = LazySpecies(QQ, "X")
+            sage: G = L.Graphs()
+            sage: P = L.SetPartitions()
+            sage: S = L.Sets()
+            sage: C = L.Cycles()
+
         TESTS::
 
             sage: LazySpecies(QQ, "X, Y, Z")._arity
@@ -649,3 +662,54 @@ def __init__(self, base_ring, names, sparse):
         """
         super().__init__(PolynomialSpecies(base_ring, names))
         self._arity = len(names)
+        if self._arity == 1:
+            self.Graphs = lambda: GraphSpecies(self)
+            self.SetPartitions = lambda: SetPartitionSpecies(self)
+            self.Sets = lambda: SetSpecies(self)
+            self.Cycles = lambda: CycleSpecies(self)
+
+
+from sage.groups.perm_gps.permgroup_named import SymmetricGroup
+class SetSpecies(LazySpeciesElement):
+    def __init__(self, parent):
+        P = parent._laurent_poly_ring
+        S = parent(SymmetricGroup)
+        super().__init__(parent, S._coeff_stream)
+
+from sage.groups.perm_gps.permgroup_named import CyclicPermutationGroup
+class CycleSpecies(LazySpeciesElement):
+    def __init__(self, parent):
+        P = parent._laurent_poly_ring
+        S = parent(lambda n: CyclicPermutationGroup(n) if n else 0)
+        super().__init__(parent, S._coeff_stream)
+
+
+from sage.graphs.graph_generators import graphs
+from sage.rings.integer_ring import ZZ
+class GraphSpecies(LazySpeciesElement):
+    def __init__(self, parent):
+        P = parent._laurent_poly_ring
+        S = parent(lambda n: sum(P(G.automorphism_group()) for G in graphs(n)))
+        super().__init__(parent, S._coeff_stream)
+
+    def isotypes(self, labels):
+        if labels in ZZ:
+            yield from graphs(labels)
+
+
+from sage.combinat.partition import Partitions
+from sage.groups.perm_gps.permgroup_named import SymmetricGroup
+from sage.combinat.set_partition import SetPartitions
+class SetPartitionSpecies(LazySpeciesElement):
+    def __init__(self, parent):
+        P = parent._laurent_poly_ring
+        E = parent(SymmetricGroup)
+        E1 = parent(lambda n: SymmetricGroup(n) if n else 0)
+        super().__init__(parent, E(E1)._coeff_stream)
+
+    def isotypes(self, labels):
+        if labels in ZZ:
+            yield from Partitions(labels)
+
+    def structures(self, labels):
+        yield from SetPartitions(labels)

From cbf8b1d91ef73c6a938901b661e885061aae3f2e Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Wed, 13 Nov 2024 14:06:44 +0100
Subject: [PATCH 19/31] fix bug in LazySpeciesElement.structures, add
 ProductSpeciesElement.structures

---
 src/sage/rings/lazy_species.py | 133 ++++++++++++++++++++++++++++-----
 1 file changed, 113 insertions(+), 20 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 65c3c140ab6..9678aa833b1 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -96,6 +96,22 @@ def weighted_vector_compositions(n_vec, d, weights_vec):
         yield from itertools.product(*map(weighted_compositions, n_vec, d_vec, weights_vec))
 
 
+def label_sets(arity, labels):
+    r"""
+    Return labels as a list of sets.
+
+    INPUT:
+
+    - ``arity`` -- the arity of the species
+    - ``labels`` -- an iterable of iterables
+    """
+    if len(labels) != arity:
+        raise ValueError("arity of must be equal to the number of arguments provided")
+
+    label_sets = [set(U) for U in labels]
+    assert all(len(U) == len(V) for U, V in zip(labels, label_sets)), f"The argument labels must be a set, but {labels} has duplicates"
+    return label_sets
+
 ######################################################################
 
 class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
@@ -350,30 +366,30 @@ def structures(self, *labels):
             sage: from sage.rings.species import PolynomialSpecies
             sage: L2 = LazySpecies(QQ, "X, Y")
             sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
-            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
-            sage: list((X*Y).structures([1],[2]))
+            sage: XY = L2(P2(PermutationGroup([], domain=[1, 2]), {0: [1], 1: [2]}))
+            sage: list((XY).structures([1], [2]))
             [((1, 2), X*Y)]
 
-            sage: list(E(X*Y).structures([1,2],[3,4]))
+            sage: list(E(XY).structures([1,2],[3,4]))
             [((3, 4, 1, 2), {((1,2)(3,4),): ({1, 2}, {3, 4})}),
              ((4, 3, 1, 2), {((1,2)(3,4),): ({1, 2}, {3, 4})})]
 
+            sage: list(XY.structures([], [1, 2]))
+            []
         """
-        fP = parent(self)
-        if len(labels) != fP._arity:
-            raise ValueError("arity of must be equal to the number of arguments provided")
-        assert all(len(U) == len(set(U)) for U in labels), f"The argument labels must be a set, but {labels} has duplicates"
-
+        labels = label_sets(self.parent()._arity, labels)
         n = tuple([len(U) for U in labels])
         S = SymmetricGroup(sum(n)).young_subgroup(n)
         F = self[sum(n)]
         l = [e for l in labels for e in l][::-1]
         for M, c in F.monomial_coefficients().items():
-            if c < 0:
+            if c not in ZZ or c < 0:
                 raise NotImplementedError("only implemented for proper non-virtual species")
+            G, dompart = M.permutation_group()
+            if not tuple(len(b) for b in dompart) == n:
+                continue
             types = [tuple(S(rep)._act_on_list_on_position(l))[::-1]
-                     for rep in libgap.RightTransversal(S, M.permutation_group()[0])]
+                     for rep in libgap.RightTransversal(S, G)]
             if c == 1:
                 for s in types:
                     yield s, M
@@ -524,9 +540,12 @@ def __init__(self, left, right):
         super().__init__(F.parent(), F._coeff_stream)
 
     def structures(self, *labels):
+        labels = label_sets(self.parent()._arity, labels)
         yield from self._left.structures(*labels)
         yield from self._right.structures(*labels)
 
+from itertools import chain, product
+from sage.combinat.subset import subsets
 class ProductSpeciesElement(LazySpeciesElement):
     def __init__(self, left, right):
         self._left = left
@@ -534,16 +553,35 @@ def __init__(self, left, right):
         F = super(LazySpeciesElement, type(left))._mul_(left, right)
         super().__init__(F.parent(), F._coeff_stream)
 
-#    def structures(self, labels):
-#        n = len(labels)
-#        l = set(labels)
-#        assert len(l) == n, f"The argument labels must be a set, but {labels} has duplicates"
-#        for k in range(n+1):
-#            for U in itertools.combinations(l, k):
-#                V = l.difference(U)
-#                yield from itertools.product(self._left.structures(U),
-#                                             self._right.structures(V))
+    def structures(self, *labels):
+        """
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L.Sets()
+            sage: C = L.Cycles()
+            sage: P = E * C
+            sage: list(P.structures([1,2]))
+            [(set(), [1, 2]), ({1}, [2]), ({2}, [1])]
+
+            sage: P = E * E
+            sage: list(P.structures([1,2]))
+            [(set(), {1, 2}), ({1}, {2}), ({2}, {1}), ({1, 2}, set())]
+
+            sage: L.<X, Y> = LazySpecies(QQ)
+            sage: list((X*Y).structures([1], [2]))
+            [(((1,), X), ((2,), Y))]
+        """
+        def dissections(s):
+            for subset in subsets(s):
+                subset = set(subset)
+                yield (subset, s - subset)
 
+        labels = label_sets(self.parent()._arity, labels)
+        for d in product(*[dissections(u) for u in labels]):
+            yield from product(self._left.structures(*[U for U, _ in d]),
+                               self._right.structures(*[V for _, V in d]))
 
 class CompositionSpeciesElement(LazySpeciesElement):
     def __init__(self, left, *args):
@@ -641,6 +679,26 @@ def __classcall_private__(cls, base_ring, names, sparse=True):
         names = normalize_names(-1, names)
         return super().__classcall__(cls, base_ring, names, sparse)
 
+    def _first_ngens(self, n):
+        r"""
+        Used by the preparser for ``F.<x> = ...``.
+
+        We do not use the generic implementation of
+        :class:`sage.combinat.CombinatorialFreeModule`, because we do
+        not want to implement `gens`.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: P.<X, Y> = LazySpecies(QQ)  # indirect doctest
+            sage: 1/(1-X-Y)
+            1 + (X+Y) + (X^2+2*X*Y+Y^2) + (X^3+3*X^2*Y+3*X*Y^2+Y^3)
+             + (X^4+4*X^3*Y+6*X^2*Y^2+4*X*Y^3+Y^4)
+             + (X^5+5*X^4*Y+10*X^3*Y^2+10*X^2*Y^3+5*X*Y^4+Y^5)
+             + (X^6+6*X^5*Y+15*X^4*Y^2+20*X^3*Y^3+15*X^2*Y^4+6*X*Y^5+Y^6) + O^7
+        """
+        return tuple([self(g) for g in self._laurent_poly_ring._first_ngens(n)])
+
     def __init__(self, base_ring, names, sparse):
         """
         EXAMPLES::
@@ -676,13 +734,47 @@ def __init__(self, parent):
         S = parent(SymmetricGroup)
         super().__init__(parent, S._coeff_stream)
 
+    def structures(self, *labels):
+        """
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: E = L.Sets()
+            sage: list(E.structures([1,2,3]))
+            [{1, 2, 3}]
+        """
+        labels = label_sets(self.parent()._arity, labels)
+        yield labels[0]
+
 from sage.groups.perm_gps.permgroup_named import CyclicPermutationGroup
+from sage.combinat.permutation import CyclicPermutations
 class CycleSpecies(LazySpeciesElement):
     def __init__(self, parent):
         P = parent._laurent_poly_ring
         S = parent(lambda n: CyclicPermutationGroup(n) if n else 0)
         super().__init__(parent, S._coeff_stream)
 
+    def structures(self, *labels):
+        """
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(ZZ, "X")
+            sage: C = L.Cycles()
+            sage: list(C.structures([]))
+            []
+            sage: list(C.structures([1]))
+            [[1]]
+            sage: list(C.structures([1,2]))
+            [[1, 2]]
+            sage: list(C.structures([1,2,3]))
+            [[1, 2, 3], [1, 3, 2]]
+        """
+        labels = label_sets(self.parent()._arity, labels)
+        yield from CyclicPermutations(labels[0])
 
 from sage.graphs.graph_generators import graphs
 from sage.rings.integer_ring import ZZ
@@ -712,4 +804,5 @@ def isotypes(self, labels):
             yield from Partitions(labels)
 
     def structures(self, labels):
+        labels = label_sets(self.parent()._arity, labels)
         yield from SetPartitions(labels)

From f4849eecd703674ff5aee7caf53a15cadc3ce1e1 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Wed, 13 Nov 2024 15:18:39 +0100
Subject: [PATCH 20/31] cosmetics

---
 src/sage/rings/lazy_species.py | 104 +++++++++++++++------------------
 1 file changed, 47 insertions(+), 57 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 9678aa833b1..ca9fb675c89 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -1,19 +1,27 @@
+from sage.rings.integer_ring import ZZ
 from sage.rings.lazy_series import LazyCompletionGradedAlgebraElement, LazyModuleElement
-from sage.rings.lazy_series_ring import LazyCompletionGradedAlgebra
+from sage.rings.lazy_series_ring import (LazyCompletionGradedAlgebra,
+                                         LazyPowerSeriesRing,
+                                         LazySymmetricFunctions)
+from sage.rings.species import PolynomialSpecies
 from sage.data_structures.stream import (Stream_zero,
                                          Stream_exact,
                                          Stream_function)
-from sage.rings.species import PolynomialSpecies
-from sage.libs.gap.libgap import libgap
 from sage.categories.sets_cat import cartesian_product
 from sage.categories.tensor import tensor
 from sage.combinat.integer_vector import IntegerVectors
+from sage.combinat.subset import subsets
+from sage.combinat.sf.sf import SymmetricFunctions
+from sage.combinat.partition import Partitions
+from sage.combinat.permutation import CyclicPermutations
+from sage.combinat.set_partition import SetPartitions
+from sage.graphs.graph_generators import graphs
+from sage.groups.perm_gps.permgroup_named import SymmetricGroup, CyclicPermutationGroup
+from sage.libs.gap.libgap import libgap
 from sage.structure.element import parent
 import itertools
 from collections import defaultdict
-from sage.rings.lazy_series_ring import LazyPowerSeriesRing, LazySymmetricFunctions
-from sage.combinat.sf.sf import SymmetricFunctions
-from sage.groups.perm_gps.permgroup_named import SymmetricGroup
+
 
 def weighted_compositions(n, d, weights, offset=0):
     r"""
@@ -114,6 +122,7 @@ def label_sets(arity, labels):
 
 ######################################################################
 
+
 class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
     r"""
 
@@ -271,14 +280,17 @@ def cycle_index_series(self):
         p = SymmetricFunctions(P.base_ring().fraction_field()).p()
         if P._arity == 1:
             L = LazySymmetricFunctions(p)
+
             def coefficient(n):
                 return sum(c * M.permutation_group()[0].cycle_index()
                            for M, c in self[n].monomial_coefficients().items())
         else:
             L = LazySymmetricFunctions(tensor([p for _ in range(P._arity)]))
+
             def coefficient(n):
                 return sum(c * M.cycle_index()
                            for M, c in self[n].monomial_coefficients().items())
+
         return L(coefficient)
 
     def _add_(self, other):
@@ -314,14 +326,15 @@ def _mul_(self, other):
             sage: from sage.rings.lazy_species import LazySpecies
             sage: L = LazySpecies(ZZ, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
-            sage: list((E^2).structures([1,2,3]))
+            sage: sorted((E^2).structures([1,2,3]))
             [(((), 1), ((1, 2, 3), E_3)),
              (((1,), X), ((2, 3), E_2)),
-             (((2,), X), ((1, 3), E_2)),
-             (((3,), X), ((1, 2), E_2)),
              (((1, 2), E_2), ((3,), X)),
+             (((1, 2, 3), E_3), ((), 1)),
              (((1, 3), E_2), ((2,), X)),
-             (((2, 3), E_2), ((1,), X))]
+             (((2,), X), ((1, 3), E_2)),
+             (((2, 3), E_2), ((1,), X)),
+             (((3,), X), ((1, 2), E_2))]
         """
         return ProductSpeciesElement(self, other)
 
@@ -348,31 +361,31 @@ def structures(self, *labels):
             [((1, 2, 3), C_3), ((2, 1, 3), C_3)]
 
             sage: F = 1/(2-E)
-            sage: list(F.structures([1,2,3]))
+            sage: sorted(F.structures([1,2,3]))
             [((1, 2, 3), E_3),
              ((1, 2, 3), X*E_2, 0),
-             ((3, 2, 1), X*E_2, 0),
-             ((2, 3, 1), X*E_2, 0),
              ((1, 2, 3), X*E_2, 1),
-             ((3, 2, 1), X*E_2, 1),
-             ((2, 3, 1), X*E_2, 1),
              ((1, 2, 3), X^3),
-             ((3, 2, 1), X^3),
-             ((3, 1, 2), X^3),
+             ((1, 3, 2), X^3),
              ((2, 1, 3), X^3),
+             ((2, 3, 1), X*E_2, 0),
+             ((2, 3, 1), X*E_2, 1),
              ((2, 3, 1), X^3),
-             ((1, 3, 2), X^3)]
+             ((3, 1, 2), X*E_2, 0),
+             ((3, 1, 2), X*E_2, 1),
+             ((3, 1, 2), X^3),
+             ((3, 2, 1), X^3)]
 
             sage: from sage.rings.species import PolynomialSpecies
-            sage: L2 = LazySpecies(QQ, "X, Y")
-            sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: XY = L2(P2(PermutationGroup([], domain=[1, 2]), {0: [1], 1: [2]}))
+            sage: L = LazySpecies(QQ, "X, Y")
+            sage: P = PolynomialSpecies(QQ, "X, Y")
+            sage: XY = L(P(PermutationGroup([], domain=[1, 2]), {0: [1], 1: [2]}))
             sage: list((XY).structures([1], [2]))
             [((1, 2), X*Y)]
 
             sage: list(E(XY).structures([1,2],[3,4]))
-            [((3, 4, 1, 2), {((1,2)(3,4),): ({1, 2}, {3, 4})}),
-             ((4, 3, 1, 2), {((1,2)(3,4),): ({1, 2}, {3, 4})})]
+            [((1, 2, 3, 4), {((1,2)(3,4),): ({1, 2}, {3, 4})}),
+             ((2, 1, 3, 4), {((1,2)(3,4),): ({1, 2}, {3, 4})})]
 
             sage: list(XY.structures([], [1, 2]))
             []
@@ -466,18 +479,8 @@ def __call__(self, *args):
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: A = L.undefined(1)
             sage: A.define(X*E(A))
-            sage: A
-            X + X^2 + (X^3+X*E_2) + (X^2*E_2+2*X^4+X*E_3)
-             + (X^2*E_3+3*X^5+3*X^3*E_2+X*{((1,2)(3,4),)}+X*E_4)
-             + (X^2*E_4+2*X^2*{((1,2)(3,4),)}+6*X^4*E_2+6*X^6
-                +3*X^3*E_3+E_2^2*X^2+X*E_5)
-             + (X^2*E_5+2*E_2^2*X^3+6*X^4*E_3+12*X^7+14*X^5*E_2
-                +3*X^3*{((1,2)(3,4),)}+3*X^3*E_4
-                +X*{((3,4),(1,2),(1,3)(2,4)(5,6))}
-                +X*{((1,2)(3,4)(5,6),)}
-                +X*{((1,2)(4,5),(1,2,3)(4,5,6))}
-                +2*E_2*X^2*E_3+E_2*{((1,2)(3,4),)}*X+X*E_6)
-             + O^8
+            sage: A[5]
+            X*E_4 + X^2*E_3 + 3*X^3*E_2 + X*{((1,2)(3,4),)} + 3*X^5
 
             sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
             sage: F = E(C(A))
@@ -491,15 +494,8 @@ def __call__(self, *args):
             sage: L = LazySpecies(R, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: E1 = L(lambda n: SymmetricGroup(n) if n else 0)
-            sage: E(q*E1)
-            1 + q*X + ((q^2+q)*E_2) + ((q^3+q)*E_3+q^2*X*E_2)
-             + ((q^4+q)*E_4+q^2*P_4+q^2*X*E_3+q^3*E_2^2)
-             + ((q^5+q)*E_5+(q^4+q^3+q^2)*E_2*E_3+q^2*X*E_4+q^3*X*P_4)
-             + ((q^6+q)*E_6+q^2*{((1,2,3)(4,6),(1,4)(2,5)(3,6))}
-                +(q^5+q^3+q^2)*E_2*E_4+q^2*X*E_5
-                +q^3*{((1,4)(2,3),(1,6,3,2,5,4))}
-                +q^3*X*E_2*E_3+q^4*E_2*P_4+q^4*E_3^2)
-             + O^7
+            sage: E(q*E1)[4]
+            (q^4+q)*E_4 + q^2*P_4 + q^2*X*E_3 + q^3*E_2^2
 
         TESTS::
 
@@ -544,8 +540,7 @@ def structures(self, *labels):
         yield from self._left.structures(*labels)
         yield from self._right.structures(*labels)
 
-from itertools import chain, product
-from sage.combinat.subset import subsets
+
 class ProductSpeciesElement(LazySpeciesElement):
     def __init__(self, left, right):
         self._left = left
@@ -579,9 +574,10 @@ def dissections(s):
                 yield (subset, s - subset)
 
         labels = label_sets(self.parent()._arity, labels)
-        for d in product(*[dissections(u) for u in labels]):
-            yield from product(self._left.structures(*[U for U, _ in d]),
-                               self._right.structures(*[V for _, V in d]))
+        for d in itertools.product(*[dissections(u) for u in labels]):
+            yield from itertools.product(self._left.structures(*[U for U, _ in d]),
+                                         self._right.structures(*[V for _, V in d]))
+
 
 class CompositionSpeciesElement(LazySpeciesElement):
     def __init__(self, left, *args):
@@ -727,7 +723,6 @@ def __init__(self, base_ring, names, sparse):
             self.Cycles = lambda: CycleSpecies(self)
 
 
-from sage.groups.perm_gps.permgroup_named import SymmetricGroup
 class SetSpecies(LazySpeciesElement):
     def __init__(self, parent):
         P = parent._laurent_poly_ring
@@ -748,8 +743,7 @@ def structures(self, *labels):
         labels = label_sets(self.parent()._arity, labels)
         yield labels[0]
 
-from sage.groups.perm_gps.permgroup_named import CyclicPermutationGroup
-from sage.combinat.permutation import CyclicPermutations
+
 class CycleSpecies(LazySpeciesElement):
     def __init__(self, parent):
         P = parent._laurent_poly_ring
@@ -776,8 +770,7 @@ def structures(self, *labels):
         labels = label_sets(self.parent()._arity, labels)
         yield from CyclicPermutations(labels[0])
 
-from sage.graphs.graph_generators import graphs
-from sage.rings.integer_ring import ZZ
+
 class GraphSpecies(LazySpeciesElement):
     def __init__(self, parent):
         P = parent._laurent_poly_ring
@@ -789,9 +782,6 @@ def isotypes(self, labels):
             yield from graphs(labels)
 
 
-from sage.combinat.partition import Partitions
-from sage.groups.perm_gps.permgroup_named import SymmetricGroup
-from sage.combinat.set_partition import SetPartitions
 class SetPartitionSpecies(LazySpeciesElement):
     def __init__(self, parent):
         P = parent._laurent_poly_ring

From 92ba51fe937b51ff56eeb9775574865e92d2f7a9 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Wed, 13 Nov 2024 21:37:00 +0100
Subject: [PATCH 21/31] adapt doctest

---
 src/sage/rings/species.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index a43f41d2a75..7ccbd288439 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -2378,19 +2378,19 @@ def _exponential(self, multiplicities, degrees):
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: P = PolynomialSpecies(QQ, ["X"])
 
-            sage: c = 3/2; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            sage: c = 3/2; all((E^c)[i] == P._exponential([c], [i]) for i in range(6))
             True
 
-            sage: c = -5/3; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            sage: c = -5/3; all((E^c)[i] == P._exponential([c], [i]) for i in range(6))
             True
 
-            sage: c = 0; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            sage: c = 0; all((E^c)[i] == P._exponential([c], [i]) for i in range(6))
             True
 
-            sage: c = 1; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            sage: c = 1; all((E^c)[i] == P._exponential([c], [i]) for i in range(6))
             True
 
-            sage: c = -1; all((E^c)[i] == P.exponential([c], [i]) for i in range(6))
+            sage: c = -1; all((E^c)[i] == P._exponential([c], [i]) for i in range(6))
             True
 
             sage: P = PolynomialSpecies(QQ, ["X"])

From 84086df2adbf8f02658ccaceaffe86dbdeca6803 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Thu, 14 Nov 2024 23:21:41 +0100
Subject: [PATCH 22/31] preliminary support for generating structures of a
 composition

---
 src/sage/combinat/set_partition.py |  26 +++++--
 src/sage/rings/lazy_species.py     | 106 +++++++++++++++++++++--------
 2 files changed, 100 insertions(+), 32 deletions(-)

diff --git a/src/sage/combinat/set_partition.py b/src/sage/combinat/set_partition.py
index 092e2f9dbb7..5078003c6f2 100644
--- a/src/sage/combinat/set_partition.py
+++ b/src/sage/combinat/set_partition.py
@@ -73,7 +73,11 @@ def _repr_(self):
             sage: S([[1,3],[2,4]])
             {{1, 3}, {2, 4}}
         """
-        return '{' + ', '.join('{' + repr(sorted(x))[1:-1] + '}' for x in self) + '}'
+        try:
+            s = [sorted(x) for x in self]
+        except TypeError:
+            s = [sorted(x, key=str) for x in self]
+        return '{' + ', '.join('{' + repr(x)[1:-1] + '}' for x in s) + '}'
 
     def __hash__(self):
         """
@@ -532,7 +536,7 @@ def pre_conjugate(sp):
 
 
 class SetPartition(AbstractSetPartition,
-        metaclass=InheritComparisonClasscallMetaclass):
+                   metaclass=InheritComparisonClasscallMetaclass):
     r"""
     A partition of a set.
 
@@ -620,7 +624,11 @@ def __init__(self, parent, s, check=True):
             {}
         """
         self._latex_options = {}
-        ClonableArray.__init__(self, parent, sorted(map(frozenset, s), key=min), check=check)
+        try:
+            s = sorted(map(frozenset, s), key=min)
+        except TypeError:
+            s = sorted(map(frozenset, s), key=lambda b: min(str(b)))
+        ClonableArray.__init__(self, parent, s, check=check)
 
     def check(self):
         """
@@ -2821,7 +2829,11 @@ def __iter__(self):
             sage: SetPartitions(["a", "b"]).list()
             [{{'a', 'b'}}, {{'a'}, {'b'}}]
         """
-        for sp in set_partition_iterator(sorted(self._set)):
+        try:
+            s = sorted(self._set)
+        except TypeError:
+            s = sorted(self._set, key=str)
+        for sp in set_partition_iterator(s):
             yield self.element_class(self, sp, check=False)
 
     def base_set(self):
@@ -3179,7 +3191,11 @@ def __iter__(self):
             sage: SetPartitions(["a", "b", "c"], 2).list()
             [{{'a', 'c'}, {'b'}}, {{'a'}, {'b', 'c'}}, {{'a', 'b'}, {'c'}}]
         """
-        for sp in set_partition_iterator_blocks(sorted(self._set), self._k):
+        try:
+            s = sorted(self._set)
+        except TypeError:
+            s = sorted(self._set, key=str)
+        for sp in set_partition_iterator_blocks(s, self._k):
             yield self.element_class(self, sp, check=False)
 
     def __contains__(self, x):
diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index ca9fb675c89..33004a29f91 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -468,15 +468,13 @@ def __call__(self, *args):
             sage: from sage.rings.species import PolynomialSpecies
             sage: P = PolynomialSpecies(QQ, "X")
             sage: Gc = L(lambda n: sum(P(G.automorphism_group()) for G in graphs(n) if G.is_connected()) if n else 0)
-            sage: E = L(lambda n: SymmetricGroup(n))
-
-            sage: G = L(lambda n: sum(P(G.automorphism_group()) for G in graphs(n)))
+            sage: E = L.Sets()
+            sage: G = L.Graphs()
             sage: E(Gc) - G
             O^7
 
-            sage: L = LazySpecies(QQ, "X")
-            sage: X = L(SymmetricGroup(1))
-            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: L.<X> = LazySpecies(QQ)
+            sage: E = L.Sets()
             sage: A = L.undefined(1)
             sage: A.define(X*E(A))
             sage: A[5]
@@ -487,8 +485,8 @@ def __call__(self, *args):
             sage: [sum(F[n].monomial_coefficients().values()) for n in range(1, 7)]
             [1, 3, 7, 19, 47, 130]
             sage: oeis(_)
-            0: A001372: Number of unlabeled mappings (or mapping patterns) from n points to themselves; number of unlabeled endofunctions.
-
+            0: A001372: Number of unlabeled mappings (or mapping patterns)
+             from n points to themselves; number of unlabeled endofunctions.
 
             sage: R.<q> = QQ[]
             sage: L = LazySpecies(R, "X")
@@ -499,8 +497,7 @@ def __call__(self, *args):
 
         TESTS::
 
-            sage: L = LazySpecies(QQ, "X")
-            sage: X = L(SymmetricGroup(1))
+            sage: L.<X> = LazySpecies(QQ)
             sage: E2 = L(SymmetricGroup(2))
             sage: X(X + E2)
             X + E_2 + O^8
@@ -510,18 +507,13 @@ def __call__(self, *args):
             sage: (1+E2)(X)
             1 + E_2 + O^7
 
-            sage: L = LazySpecies(QQ, "X, Y")
-            sage: P = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L(P(SymmetricGroup(1), {0: [1]}))
-            sage: Y = L(P(SymmetricGroup(1), {1: [1]}))
+            sage: L.<X,Y> = LazySpecies(QQ)
             sage: X(Y, 0)
             Y + O^8
 
             sage: L1 = LazySpecies(QQ, "X")
-            sage: L = LazySpecies(QQ, "X, Y")
-            sage: P = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L(P(SymmetricGroup(1), {0: [1]}))
-            sage: E = L1(lambda n: SymmetricGroup(n))
+            sage: E = L1.Sets()
+            sage: L.<X,Y> = LazySpecies(QQ)
             sage: E(X)
             1 + X + E_2(X) + E_3(X) + E_4(X) + E_5(X) + E_6(X) + O^7
         """
@@ -530,10 +522,10 @@ def __call__(self, *args):
 
 class SumSpeciesElement(LazySpeciesElement):
     def __init__(self, left, right):
-        self._left = left
-        self._right = right
         F = super(LazySpeciesElement, type(left))._add_(left, right)
         super().__init__(F.parent(), F._coeff_stream)
+        self._left = left
+        self._right = right
 
     def structures(self, *labels):
         labels = label_sets(self.parent()._arity, labels)
@@ -543,10 +535,10 @@ def structures(self, *labels):
 
 class ProductSpeciesElement(LazySpeciesElement):
     def __init__(self, left, right):
-        self._left = left
-        self._right = right
         F = super(LazySpeciesElement, type(left))._mul_(left, right)
         super().__init__(F.parent(), F._coeff_stream)
+        self._left = left
+        self._right = right
 
     def structures(self, *labels):
         """
@@ -652,6 +644,66 @@ def coefficient(n):
 
         coeff_stream = Stream_function(coefficient, P._sparse, sorder * gv)
         super().__init__(P, coeff_stream)
+        self._left = left
+        self._args = args
+
+    def structures(self, *labels):
+        r"""
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L = LazySpecies(QQ, "X")
+            sage: E = L.Sets()
+            sage: E1 = L.Sets(min=1)
+            sage: list(E(E1).structures([1,2,3]))
+            [({((1, 0), (2, 0), (3, 0))}, ({1, 2, 3},)),
+             ({((1, 0), (2, 0)), ((3, 0),)}, ({1, 2}, {3})),
+             ({((1, 0), (3, 0)), ((2, 0),)}, ({1, 3}, {2})),
+             ({((1, 0),), ((2, 0), (3, 0))}, ({1}, {2, 3})),
+             ({((1, 0),), ((2, 0),), ((3, 0),)}, ({1}, {2}, {3}))]
+
+            sage: C = L.Cycles()
+            sage: L.<X, Y> = LazySpecies(QQ)
+            sage: sum(1 for s in C(X*Y).structures([1,2,3], [1,2,3]))
+
+            sage: C(X*Y).generating_series()[6]
+            1/3*X^3*Y^3
+
+            sage: sum(1 for s in E(X*Y).structures([1,2,3], ["a", "b", "c"]))
+            6
+        """
+        F = self._left
+        G = self._args
+        m = len(G)  # == F.parent()._arity
+        k = self.parent()._arity  # == G[i].parent()._arity
+        labels = label_sets(k, labels)
+        # make label sets disjoint
+        U = [(e, i) for i, l in enumerate(labels) for e in l]
+
+        def split_set(C):
+            C_split = defaultdict(list)
+            for e, i in C:
+                C_split[i].append(e)
+            return [C_split[i] for i in range(k)]
+
+        Par_U = SetPartitions(U)
+        for pi in Par_U:
+            # Fix an arbitrary order of the blocks
+            pi_list = list(pi)
+            # Generate all functions chi from pi to {0, ..., m-1}
+            for chi in itertools.product(range(m), repeat=len(pi_list)):
+                chi_inv = defaultdict(list)
+                for b, i in zip(pi_list, chi):
+                    chi_inv[i].append(b)
+
+                # The set of structures is the Cartesian product of
+                # the structures in F[chi_inv[i] for i in range(m)]
+                # and for each set C in chi_inv[i] the set of
+                # structures in G_i[C]
+                F_s = F.structures(*[[tuple(b) for b in chi_inv[i]] for i in range(m)])
+                G_s = [G[i].structures(*split_set(C)) for i in range(m) for C in chi_inv[i]]
+                yield from itertools.product(F_s, itertools.product(*G_s))
 
 
 class LazySpecies(LazyCompletionGradedAlgebra):
@@ -719,14 +771,14 @@ def __init__(self, base_ring, names, sparse):
         if self._arity == 1:
             self.Graphs = lambda: GraphSpecies(self)
             self.SetPartitions = lambda: SetPartitionSpecies(self)
-            self.Sets = lambda: SetSpecies(self)
+            self.Sets = lambda min=0: SetSpecies(self, min)
             self.Cycles = lambda: CycleSpecies(self)
 
 
 class SetSpecies(LazySpeciesElement):
-    def __init__(self, parent):
+    def __init__(self, parent, min):
         P = parent._laurent_poly_ring
-        S = parent(SymmetricGroup)
+        S = parent(lambda n: SymmetricGroup(n) if n >= min else 0)
         super().__init__(parent, S._coeff_stream)
 
     def structures(self, *labels):
@@ -785,8 +837,8 @@ def isotypes(self, labels):
 class SetPartitionSpecies(LazySpeciesElement):
     def __init__(self, parent):
         P = parent._laurent_poly_ring
-        E = parent(SymmetricGroup)
-        E1 = parent(lambda n: SymmetricGroup(n) if n else 0)
+        E = parent.Sets()
+        E1 = parent.Sets(min=1)
         super().__init__(parent, E(E1)._coeff_stream)
 
     def isotypes(self, labels):

From b91806690740f293fc5687fd1377ed0975724361 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sat, 16 Nov 2024 11:02:08 +0100
Subject: [PATCH 23/31] provide structures for atomic and molecular species

---
 src/sage/rings/species.py | 97 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 97 insertions(+)

diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index f0dcf95291b..cb038fb67f8 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -43,6 +43,7 @@
 from sage.combinat.free_module import CombinatorialFreeModule
 from sage.combinat.integer_vector import IntegerVectors
 from sage.combinat.partition import Partitions, _Partitions
+from sage.combinat.set_partition_ordered import OrderedSetPartitions
 from sage.combinat.sf.sf import SymmetricFunctions
 from sage.groups.perm_gps.constructor import PermutationGroupElement
 from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic
@@ -69,6 +70,23 @@
 GAP_FAIL = libgap.eval('fail')
 
 
+def label_sets(arity, labels):
+    r"""
+    Return labels as a list of sets.
+
+    INPUT:
+
+    - ``arity`` -- the arity of the species
+    - ``labels`` -- an iterable of iterables
+    """
+    if len(labels) != arity:
+        raise ValueError("arity of must be equal to the number of arguments provided")
+
+    label_sets = [set(U) for U in labels]
+    assert all(len(U) == len(V) for U, V in zip(labels, label_sets)), f"The argument labels must be a set, but {labels} has duplicates"
+    return label_sets
+
+
 class AtomicSpeciesElement(WithEqualityById,
                            Element,
                            WithPicklingByInitArgs,
@@ -345,6 +363,46 @@ def __le__(self, other):
         """
         return self is other or self < other
 
+    def structures(self, *labels):
+        r"""
+        Iterate over the structures on the given set of labels.
+
+        Generically, this yields a list of relabelled representatives
+        of the cosets of corresponding groups.
+
+        EXAMPLES::
+
+            sage: from sage.rings.species import AtomicSpecies
+            sage: A = AtomicSpecies("X, Y")
+            sage: G = PermutationGroup([[("a", "b", "c", "d"), ("e", "f")]])
+            sage: a = A(G, {0: "abcd", 1: "ef"})
+            sage: sorted(a.structures([1, 2, 3, 4], ["a", "b"]))
+            [(1, 2, 3, 4, 'a', 'b'),
+             (1, 2, 3, 4, 'b', 'a'),
+             (1, 2, 4, 3, 'a', 'b'),
+             (1, 2, 4, 3, 'b', 'a'),
+             (1, 3, 2, 4, 'a', 'b'),
+             (1, 3, 2, 4, 'b', 'a'),
+             (1, 3, 4, 2, 'a', 'b'),
+             (1, 3, 4, 2, 'b', 'a'),
+             (1, 4, 2, 3, 'a', 'b'),
+             (1, 4, 2, 3, 'b', 'a'),
+             (1, 4, 3, 2, 'a', 'b'),
+             (1, 4, 3, 2, 'b', 'a')]
+
+            sage: G = PermutationGroup([[(2,3),(4,5)]], domain=[2,3,4,5])
+            sage: a = A(G, {0: [2, 3], 1: [4, 5]})
+            sage: sorted(a.structures([1, 2],["a", "b"]))
+            [(1, 2, 'a', 'b'), (1, 2, 'b', 'a')]
+        """
+        labels = label_sets(self.parent()._arity, labels)
+        n = tuple([len(U) for U in labels])
+        S = SymmetricGroup(sum(n)).young_subgroup(n)
+        l = [e for l in labels for e in l]
+        if self._mc == n:
+            for rep in libgap.RightTransversal(S, self._dis):
+                yield tuple(S(rep)._act_on_list_on_position(l))
+
 
 class AtomicSpecies(UniqueRepresentation, Parent):
     r"""
@@ -1520,6 +1578,45 @@ def __call__(self, *args):
             return P(PermutationGroup(gens, domain=range(1, starts[-1] + 1)),
                      pi, check=False)
 
+        def structures(self, *labels):
+            r"""
+            Iterate over the structures on the given set of labels.
+
+                sage: from sage.rings.species import MolecularSpecies
+                sage: M = MolecularSpecies("X,Y")
+                sage: a = M(PermutationGroup([(3,4),(5,)]), {0:[1,3,4], 1:[2,5]})
+                sage: a
+                X*Y^2*E_2(X)
+                sage: list(a.structures([1, 2, 3], ["a", "b"]))
+                [((1,), ('a',), ('b',), (2, 3)),
+                 ((1,), ('b',), ('a',), (2, 3)),
+                 ((2,), ('a',), ('b',), (1, 3)),
+                 ((2,), ('b',), ('a',), (1, 3)),
+                 ((3,), ('a',), ('b',), (1, 2)),
+                 ((3,), ('b',), ('a',), (1, 2))]
+
+                sage: G = PermutationGroup([[(2,3),(4,5)]])
+                sage: a = M(G, {0: [1, 2, 3], 1: [4, 5]})
+                sage: a
+                X*E_2(XY)
+                sage: list(a.structures([1, 2, 3], ["a", "b"]))
+                [((1,), (2, 3, 'a', 'b')),
+                 ((1,), (2, 3, 'b', 'a')),
+                 ((2,), (1, 3, 'a', 'b')),
+                 ((2,), (1, 3, 'b', 'a')),
+                 ((3,), (1, 2, 'a', 'b')),
+                 ((3,), (1, 2, 'b', 'a'))]
+            """
+            k = self.parent()._arity
+            labels = label_sets(k, labels)
+            atoms = [a for a, n in self._monomial.items() for _ in range(n)]
+            sizes = [a._mc for a in atoms]
+            dissections = product(*[OrderedSetPartitions(l, [mc[i] for mc in sizes])
+                                    for i, l in enumerate(labels)])
+            for d in dissections:
+                yield from product(*[a.structures(*[l[i] for l in d])
+                                     for i, a in enumerate(atoms)])
+
 
 class PolynomialSpeciesElement(CombinatorialFreeModule.Element):
     r"""

From 7af0ee84c90f72075051cb761569ea5f228b423b Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sat, 16 Nov 2024 16:16:04 +0100
Subject: [PATCH 24/31] simplify and correct structures

---
 src/sage/rings/lazy_species.py | 106 +++++++++++++++------------------
 src/sage/rings/species.py      |  16 +++--
 2 files changed, 60 insertions(+), 62 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 33004a29f91..ec60079b06c 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -17,7 +17,6 @@
 from sage.combinat.set_partition import SetPartitions
 from sage.graphs.graph_generators import graphs
 from sage.groups.perm_gps.permgroup_named import SymmetricGroup, CyclicPermutationGroup
-from sage.libs.gap.libgap import libgap
 from sage.structure.element import parent
 import itertools
 from collections import defaultdict
@@ -307,9 +306,9 @@ def _add_(self, other):
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: F = L(lambda n: SymmetricGroup(n))
             sage: list(E.structures([1,2,3]))
-            [((1, 2, 3), E_3)]
+            [(E_3, ((1, 2, 3),))]
             sage: list((E+F).structures([1,2,3]))
-            [((1, 2, 3), E_3), ((1, 2, 3), E_3)]
+            [(E_3, ((1, 2, 3),)), (E_3, ((1, 2, 3),))]
 
         """
         return SumSpeciesElement(self, other)
@@ -327,14 +326,14 @@ def _mul_(self, other):
             sage: L = LazySpecies(ZZ, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: sorted((E^2).structures([1,2,3]))
-            [(((), 1), ((1, 2, 3), E_3)),
-             (((1,), X), ((2, 3), E_2)),
-             (((1, 2), E_2), ((3,), X)),
-             (((1, 2, 3), E_3), ((), 1)),
-             (((1, 3), E_2), ((2,), X)),
-             (((2,), X), ((1, 3), E_2)),
-             (((2, 3), E_2), ((1,), X)),
-             (((3,), X), ((1, 2), E_2))]
+            [((1, ()), (E_3, ((1, 2, 3),))),
+             ((X, ((1,),)), (E_2, ((2, 3),))),
+             ((X, ((2,),)), (E_2, ((1, 3),))),
+             ((X, ((3,),)), (E_2, ((1, 2),))),
+             ((E_2, ((1, 2),)), (X, ((3,),))),
+             ((E_2, ((1, 3),)), (X, ((2,),))),
+             ((E_2, ((2, 3),)), (X, ((1,),))),
+             ((E_3, ((1, 2, 3),)), (1, ()))]
         """
         return ProductSpeciesElement(self, other)
 
@@ -354,61 +353,55 @@ def structures(self, *labels):
             sage: L = LazySpecies(QQ, "X")
             sage: E = L(lambda n: SymmetricGroup(n))
             sage: list(E.structures([1,2,3]))
-            [((1, 2, 3), E_3)]
+            [(E_3, ((1, 2, 3),))]
 
             sage: P = L(lambda n: CyclicPermutationGroup(n))
             sage: list(P.structures([1,2,3]))
-            [((1, 2, 3), C_3), ((2, 1, 3), C_3)]
+            [(C_3, ((1, 2, 3),)), (C_3, ((1, 3, 2),))]
 
             sage: F = 1/(2-E)
             sage: sorted(F.structures([1,2,3]))
-            [((1, 2, 3), E_3),
-             ((1, 2, 3), X*E_2, 0),
-             ((1, 2, 3), X*E_2, 1),
-             ((1, 2, 3), X^3),
-             ((1, 3, 2), X^3),
-             ((2, 1, 3), X^3),
-             ((2, 3, 1), X*E_2, 0),
-             ((2, 3, 1), X*E_2, 1),
-             ((2, 3, 1), X^3),
-             ((3, 1, 2), X*E_2, 0),
-             ((3, 1, 2), X*E_2, 1),
-             ((3, 1, 2), X^3),
-             ((3, 2, 1), X^3)]
+            [(E_3, ((1, 2, 3),)),
+             (X*E_2, ((1,), (2, 3)), 0),
+             (X*E_2, ((1,), (2, 3)), 1),
+             (X*E_2, ((2,), (1, 3)), 0),
+             (X*E_2, ((2,), (1, 3)), 1),
+             (X*E_2, ((3,), (1, 2)), 0),
+             (X*E_2, ((3,), (1, 2)), 1),
+             (X^3, ((1,), (2,), (3,))),
+             (X^3, ((1,), (3,), (2,))),
+             (X^3, ((2,), (1,), (3,))),
+             (X^3, ((2,), (3,), (1,))),
+             (X^3, ((3,), (1,), (2,))),
+             (X^3, ((3,), (2,), (1,)))]
 
             sage: from sage.rings.species import PolynomialSpecies
             sage: L = LazySpecies(QQ, "X, Y")
             sage: P = PolynomialSpecies(QQ, "X, Y")
             sage: XY = L(P(PermutationGroup([], domain=[1, 2]), {0: [1], 1: [2]}))
-            sage: list((XY).structures([1], [2]))
-            [((1, 2), X*Y)]
+            sage: list((XY).structures([1], ["a"]))
+            [(X*Y, ((1,), ('a',)))]
 
-            sage: list(E(XY).structures([1,2],[3,4]))
-            [((1, 2, 3, 4), {((1,2)(3,4),): ({1, 2}, {3, 4})}),
-             ((2, 1, 3, 4), {((1,2)(3,4),): ({1, 2}, {3, 4})})]
+            sage: sorted(E(XY).structures([1,2], [3, 4]))  # random
+            [((E_2, ((((2, 'X'), (3, 'Y')), ((1, 'X'), (4, 'Y'))),)),
+              ((X*Y, ((1,), (4,))), (X*Y, ((2,), (3,))))),
+             ((E_2, ((((3, 'Y'), (1, 'X')), ((2, 'X'), (4, 'Y'))),)),
+              ((X*Y, ((1,), (3,))), (X*Y, ((2,), (4,)))))]
 
             sage: list(XY.structures([], [1, 2]))
             []
         """
         labels = label_sets(self.parent()._arity, labels)
-        n = tuple([len(U) for U in labels])
-        S = SymmetricGroup(sum(n)).young_subgroup(n)
-        F = self[sum(n)]
-        l = [e for l in labels for e in l][::-1]
+        F = self[sum(map(len, labels))]
         for M, c in F.monomial_coefficients().items():
             if c not in ZZ or c < 0:
                 raise NotImplementedError("only implemented for proper non-virtual species")
-            G, dompart = M.permutation_group()
-            if not tuple(len(b) for b in dompart) == n:
-                continue
-            types = [tuple(S(rep)._act_on_list_on_position(l))[::-1]
-                     for rep in libgap.RightTransversal(S, G)]
             if c == 1:
-                for s in types:
-                    yield s, M
+                for s in M.structures(*labels):
+                    yield M, s
             else:
-                for e, s in cartesian_product([range(c), types]):
-                    yield s, M, e
+                for e, s in cartesian_product([range(c), M.structures(*labels)]):
+                    yield M, s, e
 
     def isotypes(self, labels):
         pass
@@ -484,7 +477,7 @@ def __call__(self, *args):
             sage: F = E(C(A))
             sage: [sum(F[n].monomial_coefficients().values()) for n in range(1, 7)]
             [1, 3, 7, 19, 47, 130]
-            sage: oeis(_)
+            sage: oeis(_)  # optional -- internet
             0: A001372: Number of unlabeled mappings (or mapping patterns)
              from n points to themselves; number of unlabeled endofunctions.
 
@@ -558,7 +551,7 @@ def structures(self, *labels):
 
             sage: L.<X, Y> = LazySpecies(QQ)
             sage: list((X*Y).structures([1], [2]))
-            [(((1,), X), ((2,), Y))]
+            [((X, ((1,),)), (Y, ((2,),)))]
         """
         def dissections(s):
             for subset in subsets(s):
@@ -656,16 +649,17 @@ def structures(self, *labels):
             sage: L = LazySpecies(QQ, "X")
             sage: E = L.Sets()
             sage: E1 = L.Sets(min=1)
-            sage: list(E(E1).structures([1,2,3]))
-            [({((1, 0), (2, 0), (3, 0))}, ({1, 2, 3},)),
-             ({((1, 0), (2, 0)), ((3, 0),)}, ({1, 2}, {3})),
-             ({((1, 0), (3, 0)), ((2, 0),)}, ({1, 3}, {2})),
-             ({((1, 0),), ((2, 0), (3, 0))}, ({1}, {2, 3})),
-             ({((1, 0),), ((2, 0),), ((3, 0),)}, ({1}, {2}, {3}))]
+            sage: sorted(E(E1).structures([1,2,3]))  # random
+            [({((2, 'X'), (3, 'X'), (1, 'X'))}, ({1, 2, 3},)),
+             ({((2, 'X'), (1, 'X')), ((3, 'X'),)}, ({1, 2}, {3})),
+             ({((2, 'X'),), ((3, 'X'), (1, 'X'))}, ({1, 3}, {2})),
+             ({((1, 'X'),), ((2, 'X'), (3, 'X'))}, ({1}, {2, 3})),
+             ({((1, 'X'),), ((2, 'X'),), ((3, 'X'),)}, ({1}, {2}, {3}))]
 
             sage: C = L.Cycles()
             sage: L.<X, Y> = LazySpecies(QQ)
             sage: sum(1 for s in C(X*Y).structures([1,2,3], [1,2,3]))
+            12
 
             sage: C(X*Y).generating_series()[6]
             1/3*X^3*Y^3
@@ -677,15 +671,16 @@ def structures(self, *labels):
         G = self._args
         m = len(G)  # == F.parent()._arity
         k = self.parent()._arity  # == G[i].parent()._arity
+        names = self.parent()._laurent_poly_ring._indices._indices._names
         labels = label_sets(k, labels)
         # make label sets disjoint
-        U = [(e, i) for i, l in enumerate(labels) for e in l]
+        U = [(e, i) for l, i in zip(labels, names) for e in l]
 
         def split_set(C):
             C_split = defaultdict(list)
             for e, i in C:
                 C_split[i].append(e)
-            return [C_split[i] for i in range(k)]
+            return [C_split[i] for i in names]
 
         Par_U = SetPartitions(U)
         for pi in Par_U:
@@ -777,7 +772,6 @@ def __init__(self, base_ring, names, sparse):
 
 class SetSpecies(LazySpeciesElement):
     def __init__(self, parent, min):
-        P = parent._laurent_poly_ring
         S = parent(lambda n: SymmetricGroup(n) if n >= min else 0)
         super().__init__(parent, S._coeff_stream)
 
@@ -798,7 +792,6 @@ def structures(self, *labels):
 
 class CycleSpecies(LazySpeciesElement):
     def __init__(self, parent):
-        P = parent._laurent_poly_ring
         S = parent(lambda n: CyclicPermutationGroup(n) if n else 0)
         super().__init__(parent, S._coeff_stream)
 
@@ -836,7 +829,6 @@ def isotypes(self, labels):
 
 class SetPartitionSpecies(LazySpeciesElement):
     def __init__(self, parent):
-        P = parent._laurent_poly_ring
         E = parent.Sets()
         E1 = parent.Sets(min=1)
         super().__init__(parent, E(E1)._coeff_stream)
diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index cb038fb67f8..a9dc8236ee5 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -1587,7 +1587,7 @@ def structures(self, *labels):
                 sage: a = M(PermutationGroup([(3,4),(5,)]), {0:[1,3,4], 1:[2,5]})
                 sage: a
                 X*Y^2*E_2(X)
-                sage: list(a.structures([1, 2, 3], ["a", "b"]))
+                sage: sorted(a.structures([1, 2, 3], ["a", "b"]))
                 [((1,), ('a',), ('b',), (2, 3)),
                  ((1,), ('b',), ('a',), (2, 3)),
                  ((2,), ('a',), ('b',), (1, 3)),
@@ -1599,7 +1599,7 @@ def structures(self, *labels):
                 sage: a = M(G, {0: [1, 2, 3], 1: [4, 5]})
                 sage: a
                 X*E_2(XY)
-                sage: list(a.structures([1, 2, 3], ["a", "b"]))
+                sage: sorted(a.structures([1, 2, 3], ["a", "b"]))
                 [((1,), (2, 3, 'a', 'b')),
                  ((1,), (2, 3, 'b', 'a')),
                  ((2,), (1, 3, 'a', 'b')),
@@ -1611,9 +1611,15 @@ def structures(self, *labels):
             labels = label_sets(k, labels)
             atoms = [a for a, n in self._monomial.items() for _ in range(n)]
             sizes = [a._mc for a in atoms]
-            dissections = product(*[OrderedSetPartitions(l, [mc[i] for mc in sizes])
-                                    for i, l in enumerate(labels)])
-            for d in dissections:
+            try:
+                # TODO: maybe OrderedSetPartitions should not raise
+                # an error if the second argument is a composition,
+                # but not of the right size
+                dissections = [OrderedSetPartitions(l, [mc[i] for mc in sizes])
+                               for i, l in enumerate(labels)]
+            except ValueError:
+                return
+            for d in product(*dissections):
                 yield from product(*[a.structures(*[l[i] for l in d])
                                      for i, a in enumerate(atoms)])
 

From 10c90e0af8f2de08201e75ad63cc8603be4ba24a Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sat, 16 Nov 2024 17:54:14 +0100
Subject: [PATCH 25/31] try to make structures more predictable

---
 src/sage/rings/lazy_species.py | 71 ++++++++++++++--------------------
 src/sage/rings/species.py      | 23 ++++++-----
 2 files changed, 43 insertions(+), 51 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index ec60079b06c..284ce81c521 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -3,7 +3,7 @@
 from sage.rings.lazy_series_ring import (LazyCompletionGradedAlgebra,
                                          LazyPowerSeriesRing,
                                          LazySymmetricFunctions)
-from sage.rings.species import PolynomialSpecies
+from sage.rings.species import PolynomialSpecies, _label_sets
 from sage.data_structures.stream import (Stream_zero,
                                          Stream_exact,
                                          Stream_function)
@@ -102,23 +102,6 @@ def weighted_vector_compositions(n_vec, d, weights_vec):
     for d_vec in IntegerVectors(d, length=k):
         yield from itertools.product(*map(weighted_compositions, n_vec, d_vec, weights_vec))
 
-
-def label_sets(arity, labels):
-    r"""
-    Return labels as a list of sets.
-
-    INPUT:
-
-    - ``arity`` -- the arity of the species
-    - ``labels`` -- an iterable of iterables
-    """
-    if len(labels) != arity:
-        raise ValueError("arity of must be equal to the number of arguments provided")
-
-    label_sets = [set(U) for U in labels]
-    assert all(len(U) == len(V) for U, V in zip(labels, label_sets)), f"The argument labels must be a set, but {labels} has duplicates"
-    return label_sets
-
 ######################################################################
 
 
@@ -382,16 +365,16 @@ def structures(self, *labels):
             sage: list((XY).structures([1], ["a"]))
             [(X*Y, ((1,), ('a',)))]
 
-            sage: sorted(E(XY).structures([1,2], [3, 4]))  # random
-            [((E_2, ((((2, 'X'), (3, 'Y')), ((1, 'X'), (4, 'Y'))),)),
-              ((X*Y, ((1,), (4,))), (X*Y, ((2,), (3,))))),
-             ((E_2, ((((3, 'Y'), (1, 'X')), ((2, 'X'), (4, 'Y'))),)),
-              ((X*Y, ((1,), (3,))), (X*Y, ((2,), (4,)))))]
+            sage: sorted(E(XY).structures([1,2], [3, 4]))
+            [((E_2, ((((1, 'X'), (3, 'Y')), ((2, 'X'), (4, 'Y'))),)),
+              ((X*Y, ((1,), (3,))), (X*Y, ((2,), (4,))))),
+             ((E_2, ((((1, 'X'), (4, 'Y')), ((2, 'X'), (3, 'Y'))),)),
+              ((X*Y, ((1,), (4,))), (X*Y, ((2,), (3,)))))]
 
             sage: list(XY.structures([], [1, 2]))
             []
         """
-        labels = label_sets(self.parent()._arity, labels)
+        labels = _label_sets(self.parent()._arity, labels)
         F = self[sum(map(len, labels))]
         for M, c in F.monomial_coefficients().items():
             if c not in ZZ or c < 0:
@@ -521,7 +504,7 @@ def __init__(self, left, right):
         self._right = right
 
     def structures(self, *labels):
-        labels = label_sets(self.parent()._arity, labels)
+        labels = _label_sets(self.parent()._arity, labels)
         yield from self._left.structures(*labels)
         yield from self._right.structures(*labels)
 
@@ -543,11 +526,11 @@ def structures(self, *labels):
             sage: C = L.Cycles()
             sage: P = E * C
             sage: list(P.structures([1,2]))
-            [(set(), [1, 2]), ({1}, [2]), ({2}, [1])]
+            [((), [1, 2]), ((1,), [2]), ((2,), [1])]
 
             sage: P = E * E
             sage: list(P.structures([1,2]))
-            [(set(), {1, 2}), ({1}, {2}), ({2}, {1}), ({1, 2}, set())]
+            [((), (1, 2)), ((1,), (2,)), ((2,), (1,)), ((1, 2), ())]
 
             sage: L.<X, Y> = LazySpecies(QQ)
             sage: list((X*Y).structures([1], [2]))
@@ -555,10 +538,10 @@ def structures(self, *labels):
         """
         def dissections(s):
             for subset in subsets(s):
-                subset = set(subset)
-                yield (subset, s - subset)
+                subset_set = set(subset)
+                yield (subset, tuple([e for e in s if e not in subset_set]))
 
-        labels = label_sets(self.parent()._arity, labels)
+        labels = _label_sets(self.parent()._arity, labels)
         for d in itertools.product(*[dissections(u) for u in labels]):
             yield from itertools.product(self._left.structures(*[U for U, _ in d]),
                                          self._right.structures(*[V for _, V in d]))
@@ -649,12 +632,12 @@ def structures(self, *labels):
             sage: L = LazySpecies(QQ, "X")
             sage: E = L.Sets()
             sage: E1 = L.Sets(min=1)
-            sage: sorted(E(E1).structures([1,2,3]))  # random
-            [({((2, 'X'), (3, 'X'), (1, 'X'))}, ({1, 2, 3},)),
-             ({((2, 'X'), (1, 'X')), ((3, 'X'),)}, ({1, 2}, {3})),
-             ({((2, 'X'),), ((3, 'X'), (1, 'X'))}, ({1, 3}, {2})),
-             ({((1, 'X'),), ((2, 'X'), (3, 'X'))}, ({1}, {2, 3})),
-             ({((1, 'X'),), ((2, 'X'),), ((3, 'X'),)}, ({1}, {2}, {3}))]
+            sage: sorted(E(E1).structures([1,2,3]))
+            [((((1, 'X'),), ((2, 'X'),), ((3, 'X'),)), ((1,), (2,), (3,))),
+             ((((1, 'X'),), ((2, 'X'), (3, 'X'))), ((1,), (2, 3))),
+             ((((1, 'X'), (2, 'X')), ((3, 'X'),)), ((1, 2), (3,))),
+             ((((1, 'X'), (2, 'X'), (3, 'X')),), ((1, 2, 3),)),
+             ((((1, 'X'), (3, 'X')), ((2, 'X'),)), ((1, 3), (2,)))]
 
             sage: C = L.Cycles()
             sage: L.<X, Y> = LazySpecies(QQ)
@@ -672,7 +655,7 @@ def structures(self, *labels):
         m = len(G)  # == F.parent()._arity
         k = self.parent()._arity  # == G[i].parent()._arity
         names = self.parent()._laurent_poly_ring._indices._indices._names
-        labels = label_sets(k, labels)
+        labels = _label_sets(k, labels)
         # make label sets disjoint
         U = [(e, i) for l, i in zip(labels, names) for e in l]
 
@@ -685,7 +668,11 @@ def split_set(C):
         Par_U = SetPartitions(U)
         for pi in Par_U:
             # Fix an arbitrary order of the blocks
-            pi_list = list(pi)
+            try:
+                pi_list = sorted([sorted(b) for b in pi])
+            except TypeError:
+                pi_list = sorted([sorted(b, key=str) for b in pi], key=str)
+
             # Generate all functions chi from pi to {0, ..., m-1}
             for chi in itertools.product(range(m), repeat=len(pi_list)):
                 chi_inv = defaultdict(list)
@@ -784,9 +771,9 @@ def structures(self, *labels):
             sage: L = LazySpecies(ZZ, "X")
             sage: E = L.Sets()
             sage: list(E.structures([1,2,3]))
-            [{1, 2, 3}]
+            [(1, 2, 3)]
         """
-        labels = label_sets(self.parent()._arity, labels)
+        labels = _label_sets(self.parent()._arity, labels)
         yield labels[0]
 
 
@@ -812,7 +799,7 @@ def structures(self, *labels):
             sage: list(C.structures([1,2,3]))
             [[1, 2, 3], [1, 3, 2]]
         """
-        labels = label_sets(self.parent()._arity, labels)
+        labels = _label_sets(self.parent()._arity, labels)
         yield from CyclicPermutations(labels[0])
 
 
@@ -838,5 +825,5 @@ def isotypes(self, labels):
             yield from Partitions(labels)
 
     def structures(self, labels):
-        labels = label_sets(self.parent()._arity, labels)
+        labels = _label_sets(self.parent()._arity, labels)
         yield from SetPartitions(labels)
diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index a9dc8236ee5..fa5a7e9da0e 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -70,9 +70,9 @@
 GAP_FAIL = libgap.eval('fail')
 
 
-def label_sets(arity, labels):
+def _label_sets(arity, labels):
     r"""
-    Return labels as a list of sets.
+    Return labels as a list of tuples.
 
     INPUT:
 
@@ -84,7 +84,12 @@ def label_sets(arity, labels):
 
     label_sets = [set(U) for U in labels]
     assert all(len(U) == len(V) for U, V in zip(labels, label_sets)), f"The argument labels must be a set, but {labels} has duplicates"
-    return label_sets
+    try:
+        label_sets_sorted = [sorted(x) for x in label_sets]
+    except TypeError:
+        label_sets_sorted = [sorted(x, key=str) for x in label_sets]
+
+    return [tuple(U) for U in label_sets_sorted]
 
 
 class AtomicSpeciesElement(WithEqualityById,
@@ -376,7 +381,7 @@ def structures(self, *labels):
             sage: A = AtomicSpecies("X, Y")
             sage: G = PermutationGroup([[("a", "b", "c", "d"), ("e", "f")]])
             sage: a = A(G, {0: "abcd", 1: "ef"})
-            sage: sorted(a.structures([1, 2, 3, 4], ["a", "b"]))
+            sage: list(a.structures([1, 2, 3, 4], ["a", "b"]))
             [(1, 2, 3, 4, 'a', 'b'),
              (1, 2, 3, 4, 'b', 'a'),
              (1, 2, 4, 3, 'a', 'b'),
@@ -392,10 +397,10 @@ def structures(self, *labels):
 
             sage: G = PermutationGroup([[(2,3),(4,5)]], domain=[2,3,4,5])
             sage: a = A(G, {0: [2, 3], 1: [4, 5]})
-            sage: sorted(a.structures([1, 2],["a", "b"]))
+            sage: list(a.structures([1, 2],["a", "b"]))
             [(1, 2, 'a', 'b'), (1, 2, 'b', 'a')]
         """
-        labels = label_sets(self.parent()._arity, labels)
+        labels = _label_sets(self.parent()._arity, labels)
         n = tuple([len(U) for U in labels])
         S = SymmetricGroup(sum(n)).young_subgroup(n)
         l = [e for l in labels for e in l]
@@ -1587,7 +1592,7 @@ def structures(self, *labels):
                 sage: a = M(PermutationGroup([(3,4),(5,)]), {0:[1,3,4], 1:[2,5]})
                 sage: a
                 X*Y^2*E_2(X)
-                sage: sorted(a.structures([1, 2, 3], ["a", "b"]))
+                sage: list(a.structures([1, 2, 3], ["a", "b"]))
                 [((1,), ('a',), ('b',), (2, 3)),
                  ((1,), ('b',), ('a',), (2, 3)),
                  ((2,), ('a',), ('b',), (1, 3)),
@@ -1599,7 +1604,7 @@ def structures(self, *labels):
                 sage: a = M(G, {0: [1, 2, 3], 1: [4, 5]})
                 sage: a
                 X*E_2(XY)
-                sage: sorted(a.structures([1, 2, 3], ["a", "b"]))
+                sage: list(a.structures([1, 2, 3], ["a", "b"]))
                 [((1,), (2, 3, 'a', 'b')),
                  ((1,), (2, 3, 'b', 'a')),
                  ((2,), (1, 3, 'a', 'b')),
@@ -1608,7 +1613,7 @@ def structures(self, *labels):
                  ((3,), (1, 2, 'b', 'a'))]
             """
             k = self.parent()._arity
-            labels = label_sets(k, labels)
+            labels = _label_sets(k, labels)
             atoms = [a for a, n in self._monomial.items() for _ in range(n)]
             sizes = [a._mc for a in atoms]
             try:

From 568df5feb46835e3cb9e5fd8ec0a5d4fc7466ce5 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sun, 17 Nov 2024 00:55:38 +0100
Subject: [PATCH 26/31] add missing 'EXAMPLES::'

---
 src/sage/rings/species.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/sage/rings/species.py b/src/sage/rings/species.py
index fa5a7e9da0e..4903cdba2e4 100644
--- a/src/sage/rings/species.py
+++ b/src/sage/rings/species.py
@@ -1587,6 +1587,8 @@ def structures(self, *labels):
             r"""
             Iterate over the structures on the given set of labels.
 
+            EXAMPLES::
+
                 sage: from sage.rings.species import MolecularSpecies
                 sage: M = MolecularSpecies("X,Y")
                 sage: a = M(PermutationGroup([(3,4),(5,)]), {0:[1,3,4], 1:[2,5]})

From a3525df107cc0ff4242a8fbb8df3b44ffce60607 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Mon, 16 Dec 2024 07:45:40 +0100
Subject: [PATCH 27/31] add (currently failing) tests

---
 src/sage/rings/lazy_species.py | 20 +++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 284ce81c521..31ebe005328 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -492,6 +492,12 @@ def __call__(self, *args):
             sage: L.<X,Y> = LazySpecies(QQ)
             sage: E(X)
             1 + X + E_2(X) + E_3(X) + E_4(X) + E_5(X) + E_6(X) + O^7
+
+            sage: L.<X> = LazySpecies(QQ)
+            sage: E1 = L.Sets(min=1)
+            sage: Omega = L.undefined(1)
+            sage: L.define_implicitly([Omega], [E1(Omega) - X])
+            sage: Omega[1]
         """
         return CompositionSpeciesElement(self, *args)
 
@@ -549,6 +555,14 @@ def dissections(s):
 
 class CompositionSpeciesElement(LazySpeciesElement):
     def __init__(self, left, *args):
+        """
+        TESTS::
+
+            sage: P.<X> = LazySpecies(QQ)
+            sage: X(P.zero())
+
+            sage: (1+X)(P.zero())
+        """
         fP = parent(left)
         if len(args) != fP._arity:
             raise ValueError("arity of must be equal to the number of arguments provided")
@@ -560,7 +574,11 @@ def __init__(self, left, *args):
 
         # f = 0
         if isinstance(left._coeff_stream, Stream_zero):
-            return P.zero()
+            coeff_stream = Stream_zero()
+            super().__init__(P, coeff_stream)
+            self._left = left
+            self._args = args
+            return
 
         # args = (0, ..., 0)
         if all((not isinstance(g, LazyModuleElement) and not g)

From 6b8f13950a829561754b4b9d614db3d67425f2e3 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Mon, 16 Dec 2024 11:24:11 +0100
Subject: [PATCH 28/31] add forgotten import

---
 src/sage/rings/lazy_species.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 31ebe005328..e08d5a82932 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -558,6 +558,7 @@ def __init__(self, left, *args):
         """
         TESTS::
 
+            sage: from sage.rings.lazy_species import LazySpecies
             sage: P.<X> = LazySpecies(QQ)
             sage: X(P.zero())
 

From 50d3082c22863c72bf264f2eef33fe784c3e92d3 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Thu, 26 Dec 2024 20:33:46 +0100
Subject: [PATCH 29/31] provide reversion, switch to multiplicities in
 weighted_compositions

---
 src/sage/rings/lazy_species.py | 257 ++++++++++++++++++++++-----------
 1 file changed, 171 insertions(+), 86 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index e08d5a82932..94437f62573 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -1,3 +1,4 @@
+from sage.misc.lazy_list import lazy_list
 from sage.rings.integer_ring import ZZ
 from sage.rings.lazy_series import LazyCompletionGradedAlgebraElement, LazyModuleElement
 from sage.rings.lazy_series_ring import (LazyCompletionGradedAlgebra,
@@ -22,31 +23,40 @@
 from collections import defaultdict
 
 
-def weighted_compositions(n, d, weights, offset=0):
+def weighted_compositions(n, d, weight_multiplicities, _w0=0):
     r"""
     Return all compositions of `n` of weight `d`.
 
-    The weight of a composition `n_1, n_2, \dots` is `\sum_i w_i
-    n_i`.
+    The weight of a composition `n_1, n_2, \dots` is `\sum_i w_i n_i`.
 
-    ``weights`` is assumed to be a weakly increasing list of positive
-    integers.
+    INPUT:
+
+
+    - ``n`` -- a nonnegative integer, the sum of the parts
+    - ``d`` -- a nonnegative integer, the total weight
+    - ``weight_multiplicities`` -- an iterable,
+      ``weight_multiplicities[i]`` is the number of positions with
+      weight `i+1`.
 
     EXAMPLES::
 
         sage: from sage.rings.lazy_species import weighted_compositions
-        sage: list(weighted_compositions(1, 1, [1,1,2]))
-        [[0, 1], [1]]
+        sage: list(weighted_compositions(1, 1, [2,1]))
+        [[1, 0], [0, 1]]
 
-        sage: list(weighted_compositions(2, 1, [1,1,2]))
+        sage: list(weighted_compositions(2, 1, [2,1]))
         []
 
-        sage: list(weighted_compositions(1, 2, [1,1,2,3]))
+        sage: list(weighted_compositions(1, 2, [2,1,1]))
         [[0, 0, 1]]
 
-        sage: list(weighted_compositions(3, 4, [1,1,2,2,3,3,4,4,5]))
-        [[0, 2, 0, 1], [0, 2, 1], [1, 1, 0, 1], [1, 1, 1], [2, 0, 0, 1], [2, 0, 1]]
-
+        sage: list(weighted_compositions(3, 4, [2,2]))
+        [[2, 0, 1, 0],
+         [1, 1, 1, 0],
+         [0, 2, 1, 0],
+         [2, 0, 0, 1],
+         [1, 1, 0, 1],
+         [0, 2, 0, 1]]
     """
     # the empty composition exists if and only if n == d == 0
     if not n:
@@ -56,51 +66,57 @@ def weighted_compositions(n, d, weights, offset=0):
     if not d:
         return
 
-    # otherwise we iterate over the possibilities for the first part
-    if offset < len(weights):
-        w0 = weights[offset]
-    else:
-        return
-    if w0 > d:
+    # otherwise we iterate over the possibilities for the first
+    # weight_multiplicities[_w0] parts
+    try:
+        if _w0 >= len(weight_multiplicities):
+            return
+    except TypeError:
+        pass
+    if _w0 > d:
         return
-    for i in range(min(n, d // w0) + 1):
-        for c in weighted_compositions(n - i, d - i * w0, weights, offset=offset+1):
-            yield [i] + c
+    for s in range(n + 1):
+        for c in weighted_compositions(n - s, d - s * (_w0 + 1), weight_multiplicities, _w0=_w0+1):
+            m = weight_multiplicities[_w0]
+            for v in map(list, IntegerVectors(s, length=m)):
+                yield v + c
 
-
-def weighted_vector_compositions(n_vec, d, weights_vec):
+def weighted_vector_compositions(n_vec, d, weight_multiplicities_vec):
     r"""
     Return all compositions of the vector `n` of weight `d`.
 
     INPUT:
 
-    - ``n_vec``, a `k`-tuple of non-negative integers.
+    - ``n_vec`` -- a `k`-tuple of non-negative integers.
 
-    - ``d``, a non-negative integer.
+    - ``d`` -- a non-negative integer, the total sum of the parts in
+      all components
 
-    - ``weights_vec``, `k`-tuple of weakly increasing lists of
-      positive integers.
+    - ``weight_multiplicities_vec`` -- `k`-tuple of iterables, an
+      iterable, ``weight_multiplicities_vec[j][i]`` is the number of
+      positions with weight `i+1` in the `j`-th component.
 
     EXAMPLES::
 
         sage: from sage.rings.lazy_species import weighted_vector_compositions
-        sage: list(weighted_vector_compositions([1,1], 2, [[1,1,2,3], [1,2,3]]))
-        [([0, 1], [1]), ([1], [1])]
+        sage: list(weighted_vector_compositions([1,1], 2, [[2,1,1], [1,1,1]]))
+        [([1, 0], [1]), ([0, 1], [1])]
 
-        sage: list(weighted_vector_compositions([3,1], 4, [[1,1,2,5], [1,1,2,5]]))
-        [([0, 3], [0, 1]),
-         ([0, 3], [1]),
-         ([1, 2], [0, 1]),
-         ([1, 2], [1]),
+        sage: list(weighted_vector_compositions([3,1], 4, [[2,1,0,0,1], [2,1,0,0,1]]))
+        [([3, 0], [1, 0]),
+         ([3, 0], [0, 1]),
+         ([2, 1], [1, 0]),
          ([2, 1], [0, 1]),
-         ([2, 1], [1]),
-         ([3], [0, 1]),
-         ([3], [1])]
-
+         ([1, 2], [1, 0]),
+         ([1, 2], [0, 1]),
+         ([0, 3], [1, 0]),
+         ([0, 3], [0, 1])]
     """
     k = len(n_vec)
     for d_vec in IntegerVectors(d, length=k):
-        yield from itertools.product(*map(weighted_compositions, n_vec, d_vec, weights_vec))
+        yield from itertools.product(*map(weighted_compositions,
+                                          n_vec, d_vec,
+                                          weight_multiplicities_vec))
 
 ######################################################################
 
@@ -115,7 +131,8 @@ class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
         sage: from sage.rings.lazy_species import LazySpecies
         sage: L = LazySpecies(ZZ, "X")
         sage: E = L(lambda n: SymmetricGroup(n))
-        sage: 1 / E
+        sage: E_inv = 1 / E
+        sage: E_inv
         1 + (-X) + (-E_2+X^2) + (-E_3+2*X*E_2-X^3)
           + (-E_4+2*X*E_3+E_2^2-3*X^2*E_2+X^4)
           + (-E_5+2*X*E_4+2*E_2*E_3-3*X^2*E_3-3*X*E_2^2+4*X^3*E_2-X^5)
@@ -127,7 +144,7 @@ class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
         sage: def coefficient(m):
         ....:     return sum((-1)^len(la) * multinomial((n := la.to_exp())) * prod(E[i]^ni for i, ni in enumerate(n, 1)) for la in Partitions(m))
 
-        sage: all(coefficient(m) == (1/E)[m] for m in range(10))
+        sage: all(coefficient(m) == E_inv[m] for m in range(10))
         True
     """
     def isotype_generating_series(self):
@@ -146,11 +163,7 @@ def isotype_generating_series(self):
             sage: E(C).isotype_generating_series()
             1 + X + 2*X^2 + 3*X^3 + 5*X^4 + 7*X^5 + 11*X^6 + O(X^7)
 
-            sage: from sage.rings.species import PolynomialSpecies
-            sage: L2 = LazySpecies(QQ, "X, Y")
-            sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
-            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
+            sage: L2.<X, Y> = LazySpecies(QQ)
             sage: E(X + Y).isotype_generating_series()
             1 + (X+Y) + (X^2+X*Y+Y^2) + (X^3+X^2*Y+X*Y^2+Y^3)
             + (X^4+X^3*Y+X^2*Y^2+X*Y^3+Y^4)
@@ -185,7 +198,7 @@ def generating_series(self):
 
             sage: from sage.rings.lazy_species import LazySpecies
             sage: L = LazySpecies(QQ, "X")
-            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: E = L.Sets()
             sage: E.generating_series()
             1 + X + 1/2*X^2 + 1/6*X^3 + 1/24*X^4 + 1/120*X^5 + 1/720*X^6 + O(X^7)
 
@@ -193,11 +206,7 @@ def generating_series(self):
             sage: C.generating_series()
             X + 1/2*X^2 + 1/3*X^3 + 1/4*X^4 + 1/5*X^5 + 1/6*X^6 + O(X^7)
 
-            sage: from sage.rings.species import PolynomialSpecies
-            sage: L2 = LazySpecies(QQ, "X, Y")
-            sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
-            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
+            sage: L2.<X, Y> = LazySpecies(QQ)
             sage: E(X + Y).generating_series()
             1 + (X+Y) + (1/2*X^2+X*Y+1/2*Y^2)
             + (1/6*X^3+1/2*X^2*Y+1/2*X*Y^2+1/6*Y^3)
@@ -235,7 +244,7 @@ def cycle_index_series(self):
 
             sage: from sage.rings.lazy_species import LazySpecies
             sage: L = LazySpecies(ZZ, "X")
-            sage: E = L(lambda n: SymmetricGroup(n))
+            sage: E = L.Sets()
             sage: h = SymmetricFunctions(QQ).h()
             sage: LazySymmetricFunctions(h)(E.cycle_index_series())
             h[] + h[1] + h[2] + h[3] + h[4] + h[5] + h[6] + O^7
@@ -245,13 +254,9 @@ def cycle_index_series(self):
             sage: s(C.cycle_index_series()[5])
             s[1, 1, 1, 1, 1] + s[2, 2, 1] + 2*s[3, 1, 1] + s[3, 2] + s[5]
 
-            sage: from sage.rings.species import PolynomialSpecies
             sage: L = LazySpecies(QQ, "X")
-            sage: E = L(lambda n: SymmetricGroup(n))
-            sage: L2 = LazySpecies(QQ, "X, Y")
-            sage: P2 = PolynomialSpecies(QQ, "X, Y")
-            sage: X = L2(P2(SymmetricGroup(1), {0: [1]}))
-            sage: Y = L2(P2(SymmetricGroup(1), {1: [1]}))
+            sage: E = L.Sets()
+            sage: L2.<X, Y> = LazySpecies(QQ)
             sage: E(X + Y).cycle_index_series()[3]
             1/6*p[] # p[1, 1, 1] + 1/2*p[] # p[2, 1] + 1/3*p[] # p[3]
             + 1/2*p[1] # p[1, 1] + 1/2*p[1] # p[2] + 1/2*p[1, 1] # p[1]
@@ -433,6 +438,7 @@ def polynomial(self, degree=None, names=None):
 
     def __call__(self, *args):
         """
+
         EXAMPLES::
 
             sage: from sage.rings.lazy_species import LazySpecies
@@ -493,14 +499,95 @@ def __call__(self, *args):
             sage: E(X)
             1 + X + E_2(X) + E_3(X) + E_4(X) + E_5(X) + E_6(X) + O^7
 
+        It would be extremely nice to allow the following, but this
+        poses theoretical problems::
+
             sage: L.<X> = LazySpecies(QQ)
             sage: E1 = L.Sets(min=1)
             sage: Omega = L.undefined(1)
             sage: L.define_implicitly([Omega], [E1(Omega) - X])
-            sage: Omega[1]
+            sage: Omega[1]  # not tested
         """
+        fP = self.parent()
+        if len(args) != fP._arity:
+            raise ValueError("arity of must be equal to the number of arguments provided")
+        # Find a good parent for the result
+        from sage.structure.element import get_coercion_model
+        cm = get_coercion_model()
+        P = cm.common_parent(self.base_ring(), *[parent(g) for g in args])
+        # f = 0
+        if isinstance(self._coeff_stream, Stream_zero):
+            return P.zero()
+
+        # args = (0, ..., 0)
+        if all((not isinstance(g, LazyModuleElement) and not g)
+               or (isinstance(g, LazyModuleElement)
+                   and isinstance(g._coeff_stream, Stream_zero))
+               for g in args):
+            return P(self[0])
+
         return CompositionSpeciesElement(self, *args)
 
+    def revert(self):
+        r"""
+        Return the compositional inverse of ``self``.
+
+        EXAMPLES::
+
+            sage: from sage.rings.lazy_species import LazySpecies
+            sage: L.<X> = LazySpecies(QQ)
+            sage: E1 = L.Sets(1)
+            sage: g = E1.revert()
+            sage: g[:5]
+            [X, -E_2, -E_3 + X*E_2, -E_4 + P_4 + X*E_3 - X^2*E_2]
+
+            sage: E = L.Sets()
+            sage: P = E(X*E1(-X))*(1+X) - 1
+            sage: P.revert()[:5]
+            [X, X^2, X*E_2 + 2*X^3, X*E_3 + 2*X^2*E_2 + {((1,2)(3,4),)} + 5*X^4]
+        """
+        P = self.parent()
+        if P._arity != 1:
+            raise ValueError("arity must be equal to 1")
+        coeff_stream = self._coeff_stream
+        if isinstance(coeff_stream, Stream_zero):
+            raise ValueError("compositional inverse does not exist")
+        R = P._laurent_poly_ring
+        if (isinstance(coeff_stream, Stream_exact)
+            and coeff_stream.order() >= 0
+            and coeff_stream._degree == 2):
+            # self = a + b * p_1; self.revert() = -a/b + 1/b * p_1
+            a = coeff_stream[0]
+            b = coeff_stream[1][Partition([1])]
+            X = R(Partition([1]))
+            coeff_stream = Stream_exact((-a/b, 1/b * X),
+                                        order=0)
+            return P.element_class(P, coeff_stream)
+
+        # TODO: coefficients should not be checked here, it prevents
+        # us from using self.define in some cases!
+        if coeff_stream[0]:
+            raise ValueError("cannot determine whether the compositional inverse exists")
+
+        X_mol = P._laurent_poly_ring._indices.subset(1)[0]  # as a molecular species
+        X = P(SymmetricGroup(1))  # as a lazy species
+
+        def coefficient(n):
+            if n:
+                return 0
+            c = coeff_stream[1].coefficient(X_mol)
+            if c.is_unit():
+                return ~c
+            raise ValueError("compositional inverse does not exist")
+
+        b = P(lambda n: 0 if n else coeff_stream[1].coefficient(X_mol))  # TODO: we want a lazy version of Stream_exact
+        b_inv = P(coefficient)  # TODO: we want a lazy version of Stream_exact
+        g = P.undefined(valuation=1)
+        g.define(b_inv * (X - (self - b * X)(g)))
+        return g
+
+    compositional_inverse = revert
+
 
 class SumSpeciesElement(LazySpeciesElement):
     def __init__(self, left, right):
@@ -560,34 +647,19 @@ def __init__(self, left, *args):
 
             sage: from sage.rings.lazy_species import LazySpecies
             sage: P.<X> = LazySpecies(QQ)
+            sage: P.zero()(X)
+            0
             sage: X(P.zero())
-
+            0
             sage: (1+X)(P.zero())
+            1
         """
-        fP = parent(left)
-        if len(args) != fP._arity:
-            raise ValueError("arity of must be equal to the number of arguments provided")
-
+        fP = left.parent()
         # Find a good parent for the result
         from sage.structure.element import get_coercion_model
         cm = get_coercion_model()
         P = cm.common_parent(left.base_ring(), *[parent(g) for g in args])
 
-        # f = 0
-        if isinstance(left._coeff_stream, Stream_zero):
-            coeff_stream = Stream_zero()
-            super().__init__(P, coeff_stream)
-            self._left = left
-            self._args = args
-            return
-
-        # args = (0, ..., 0)
-        if all((not isinstance(g, LazyModuleElement) and not g)
-               or (isinstance(g, LazyModuleElement)
-                   and isinstance(g._coeff_stream, Stream_zero))
-               for g in args):
-            return P(left[0])
-
         # f is a constant polynomial
         if (isinstance(left._coeff_stream, Stream_exact)
             and not left._coeff_stream._constant
@@ -607,23 +679,36 @@ def __init__(self, left, *args):
         R = P._internal_poly_ring.base_ring()
         L = fP._internal_poly_ring.base_ring()
 
+        def coeff(g, i):
+            c = g._coeff_stream[i]
+            if not isinstance(c, PolynomialSpecies.Element):
+                return R(c)
+            return c
+
+        # args_flat and weights contain one list for each g
+        weight_exp = [lazy_list(lambda j, g=g: len(coeff(g, j+1)))
+                      for g in args]
+        # work around python's scoping rules
+        def flat(g):
+            return itertools.chain.from_iterable((coeff(g, j) for j in itertools.count()))
+        args_flat1 = [lazy_list(flat(g)) for g in args]
+
         def coefficient(n):
             if not n:
                 if left[0]:
                     return R(list(left[0])[0][1])
                 return R.zero()
-            args_flat = [[(M, c) for i in range(n+1) for M, c in g[i]]
-                         for g in args]
-            weights = [[sum(M.grade()) for i in range(n+1) for M, _ in g[i]]
-                       for g in args]
             result = R.zero()
-            for i in range(n // gv + 1):
+            for i in range(1, n // gv + 1):
+                # skip i=0 because it produces a term only for n=0
+
                 # compute homogeneous components
                 lF = defaultdict(L)
                 for M, c in left[i]:
                     lF[M.grade()] += L._from_dict({M: c})
                 for mc, F in lF.items():
-                    for degrees in weighted_vector_compositions(mc, n, weights):
+                    for degrees in weighted_vector_compositions(mc, n, weight_exp):
+                        args_flat = [list(a[0:len(degrees[j])]) for j, a in enumerate(args_flat1)]
                         multiplicities = [c for alpha, g_flat in zip(degrees, args_flat)
                                           for d, (_, c) in zip(alpha, g_flat) if d]
                         molecules = [M for alpha, g_flat in zip(degrees, args_flat)

From 5340aacc173003f3ebbaa62dfdbb05de472ebbf7 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sun, 5 Jan 2025 11:31:09 +0100
Subject: [PATCH 30/31] provide non-trivial examples for molecular
 decompositions

---
 src/doc/en/reference/references/index.rst | 10 +++-
 src/sage/rings/lazy_species.py            | 69 ++++++++++++++++++++++-
 2 files changed, 75 insertions(+), 4 deletions(-)

diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst
index ca81cda3f75..ef0e3df53df 100644
--- a/src/doc/en/reference/references/index.rst
+++ b/src/doc/en/reference/references/index.rst
@@ -872,7 +872,7 @@ REFERENCES:
 .. [BrHu2019] Petter Brändén, June Huh. *Lorentzian polynomials*.
               Ann. Math. (2) 192, No. 3, 821-891 (2020).
               :arxiv:`1902.03719`, :doi:`10.4007/annals.2020.192.3.4`.
-        
+
 .. [Bru2014] Erwan Brugalle and Kristin Shaw. *A bit of tropical geometry*.
              Amer. Math. Monthly, 121(7):563-589, 2014.
 
@@ -3107,6 +3107,10 @@ REFERENCES:
             descent algebra.* Adv. Math. **77** (1989).
             http://www.lacim.uqam.ca/~christo/Publi%C3%A9s/1989/Decomposition%20Solomon.pdf
 
+.. [GL2011] Ira M. Gessel, Ji Li.
+            *Enumeration of point-determining graphs*. Journal of
+            Combinatorial Theory, Series A, 118 (2011), pp. 591--612.
+
 .. [GR1993] Ira M. Gessel, Christophe Reutenauer.
             *Counting Permutations with Given Cycle Structure and Descent Set*. Journal of
             Combinatorial Theory, Series A, 64 (1993), pp. 189--215.
@@ -5243,7 +5247,7 @@ REFERENCES:
 .. [NT2007] Serguei Norine and Robin Thomas. *Minimally Non-Pfaffian Graphs*.
             Combinatorica, vol. 27, no. 5, pages: 587 -- 600, Springer. 2007.
             :doi:`10.1016/j.jctb.2007.12.005`.
-            
+
 .. [Nur2004] K. Nurmela. *Upper bounds for covering arrays by tabu search*.
              Discrete Applied Math., 138 (2004), 143-152.
 
@@ -6193,7 +6197,7 @@ REFERENCES:
              of the chromatic polynomial of a graph*, Adv. Math., ***111***
              no.1 (1995), 166-194. :doi:`10.1006/aima.1995.1020`.
 
-.. [Sta1998] \R. P. Stanley, *Graph colorings and related symmetric functions: 
+.. [Sta1998] \R. P. Stanley, *Graph colorings and related symmetric functions:
              ideas and applications A description of results, interesting applications,
              & notable open problems*, Discrete Mathematics, 193, no.1-3, (1998),
              267-286. :doi:`10.1016/S0012-365X(98)00146-0`.
diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 94437f62573..478905017e4 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -1,3 +1,71 @@
+"""r
+Lazy Combinatorial Species
+
+We regard a combinatorial species as a sequence of group actions of
+the symmetric groups `\mathfrak S_n`, for `n\in\NN`.
+
+Coefficients of lazy species are computed on demand.  They have
+infinite precision, although equality can only be decided in special
+cases.
+
+AUTHORS:
+
+- Mainak Roy, Martin Rubey, Travis Scrimshaw (2024-2025)
+
+EXAMPLES::
+
+    We can reproduce the molecular expansions from Appendix B in
+    [GL2011]_ with little effort.  The molecular expansion of the
+    species of point determining graphs can be computed as the
+    species of graphs composed with the compositional inverse of the
+    species of non-empty sets.  To make the result more readable, we
+    provide a name for `E_2(X^2)`::
+
+        sage: from sage.rings.lazy_species import LazySpecies
+        sage: L.<X> = LazySpecies(QQ)
+        sage: E = L.Sets()
+        sage: E_2 = L(SymmetricGroup(2))
+        sage: Ep = L.Sets(1)
+        sage: G = L.Graphs()
+        sage: E_2(X^2)[4].support()[0].support()[0].rename("E_2(X^2)")
+
+    The molecular decomposition begins with::
+
+        sage: P = G(Ep.revert())
+        sage: P.truncate(6)
+        1 + X + E_2 + (E_3+X*E_2) + (E_4+X*E_3+P_4+X^2*E_2+E_2(X^2))
+        + (E_5+E_2*E_3+X*E_4+X*E_2^2+X^2*E_3+2*X*P_4+P_5+5*X*E_2(X^2)+3*X^3*E_2)
+
+    Note that [GL2011]_ write `E_2(E_2)` instead of `P_4` and `D_5`
+    instead of `P_5`, and there is apparently a misprint: `X*P_4 + 4
+    X^3 E_2` should be `2 X P_4 + 3 X^3 E_2`.
+
+    To compute the molecular decomposition of the species of
+    connected graphs with no endpoints, we use Equation (3.3) in
+    [GL2011]_.  Before that we need to define the species of
+    connected graphs::
+
+        sage: Gc = Ep.revert()(G-1)
+        sage: Mc = Gc(X*E(-X)) + E_2(-X)
+        sage: E(Mc).truncate(5)
+        1 + X + E_2 + 2*E_3 + (2*E_4+P_4+E_2^2+X*E_3)
+
+    Note that [GL2011]_ apparently contains a misprint: `2 X E_3`
+    should be `X E_3 + E_2^2`.  Indeed, the graphs on four vertices
+    without endpoints are the complete graph and the empty graph, the
+    square, the diamond graph and the triangle with an extra isolated
+    vertex.
+
+
+    To compute the molecular decomposition of the species of
+    bi-point-determining graphs we use Corollary (4.6) in
+    [GL2011]_::
+
+        sage: B = G(2*Ep.revert() - X)
+        sage: B.truncate(6)
+        1 + X + E_2(X^2) + (P_5+5*X*E_2(X^2))
+"""
+
 from sage.misc.lazy_list import lazy_list
 from sage.rings.integer_ring import ZZ
 from sage.rings.lazy_series import LazyCompletionGradedAlgebraElement, LazyModuleElement
@@ -123,7 +191,6 @@ def weighted_vector_compositions(n_vec, d, weight_multiplicities_vec):
 
 class LazySpeciesElement(LazyCompletionGradedAlgebraElement):
     r"""
-
     EXAMPLES:
 
     Compute the molecular expansion of `E(-X)`::

From f48e209996308eb50082c310a200fb1d387efb98 Mon Sep 17 00:00:00 2001
From: Martin Rubey <martin.rubey@tuwien.ac.at>
Date: Sun, 5 Jan 2025 11:36:10 +0100
Subject: [PATCH 31/31] fix doctests and syntax error

---
 src/sage/rings/lazy_species.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/sage/rings/lazy_species.py b/src/sage/rings/lazy_species.py
index 478905017e4..55e56f9bb2a 100644
--- a/src/sage/rings/lazy_species.py
+++ b/src/sage/rings/lazy_species.py
@@ -1,4 +1,4 @@
-"""r
+r"""
 Lazy Combinatorial Species
 
 We regard a combinatorial species as a sequence of group actions of
@@ -527,7 +527,7 @@ def __call__(self, *args):
             sage: A = L.undefined(1)
             sage: A.define(X*E(A))
             sage: A[5]
-            X*E_4 + X^2*E_3 + 3*X^3*E_2 + X*{((1,2)(3,4),)} + 3*X^5
+            X*E_4 + X^2*E_3 + 3*X^3*E_2 + X*E_2(X^2) + 3*X^5
 
             sage: C = L(lambda n: CyclicPermutationGroup(n) if n else 0)
             sage: F = E(C(A))
@@ -611,7 +611,7 @@ def revert(self):
             sage: E = L.Sets()
             sage: P = E(X*E1(-X))*(1+X) - 1
             sage: P.revert()[:5]
-            [X, X^2, X*E_2 + 2*X^3, X*E_3 + 2*X^2*E_2 + {((1,2)(3,4),)} + 5*X^4]
+            [X, X^2, X*E_2 + 2*X^3, X*E_3 + 2*X^2*E_2 + E_2(X^2) + 5*X^4]
         """
         P = self.parent()
         if P._arity != 1: