Implement ordered multiset partitions

[alauve]

An ordered multiset partition C of a multiset X is a list of subsets of X (not multisets), called the blocks of C, whose multi-union is X.

These objects appear in the work of

• [Haglund, Rhoades, Wilson] - arXiv:1509.07058
• [Haglund, Rhoades, Shimozono] - arXiv:1609.07575
• [Benkart, et al] - arXiv:1707.08709v2
• [Lauve, Mastnak] - arXiv:1803.02691

This module provides tools for manipulating ordered multiset partitions.

CC: @tscrim @darijgr @zabrocki @alauve @sagetrac-mshimo @anneschilling @saliola @amypang

Component: combinatorics

Keywords: IMA coding sprint, CHAs, sagedays@icerm

Author: Aaron Lauve, Anne Schilling

Branch: 085a93d

Reviewer: Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/25148

[alauve]

Description changed:

--- 
+++ 
@@ -1,4 +1,4 @@
-An ordered multiset partition C of a multiset C is a list of subsets of X (not multisets),
+An ordered multiset partition C of a multiset X is a list of subsets of X (not multisets),
 called the blocks of C, whose multi-union is X.

 These objects appear in the work of 

[alauve]

Commit: 8a2ded6

[alauve]

Branch: public/combinat/implement-ordered-multiset-partitions-25148

[tscrim]

Changed keywords from none to IMA coding sprint, CHAs

[anneschilling]

comment:5

Actually, I have the minimaj crystal implemented if you want to add that to this ticket!?

[tscrim]

comment:6

Yes, that would be great. Do you want to add it onto the current code or attach (or email me) the file with the implementation?

[anneschilling]

comment:7

Here is my code, but I have not recently checked it:

def multiset_partitions(k,ell,n):
    """
    Return the set of all multiset partitions with `k` parts
    of words of length `ell` on the alphabet `[1,2,...,n]`.

    EXAMPLES::

         sage: multiset_partitions(2,3,3)
         [({1, 2}, {1}),
         ({1, 2}, {2}),
         ({1, 2}, {3}),
         ({1, 3}, {1}),
         ({1, 3}, {2}),
         ({1, 3}, {3}),
         ({2, 3}, {1}),
         ({2, 3}, {2}),
         ({2, 3}, {3}),
         ({1}, {1, 2}),
         ({1}, {1, 3}),
         ({1}, {2, 3}),
         ({2}, {1, 2}),
         ({2}, {1, 3}),
         ({2}, {2, 3}),
         ({3}, {1, 2}),
         ({3}, {1, 3}),
         ({3}, {2, 3})]
    """
    E = range(1,n+1)
    return DisjointUnionEnumeratedSets(Family(Compositions(ell,length=k),
                                              lambda I: cartesian_product( [Subsets(E,i) for i in I] ))).list()

def greedy(b):
    """
    Return greedy algorithm on `b`.

    EXAMPLES::

        sage: b=multiset_partitions(2,3,3)[0]; b
        ({1, 2}, {1})
        sage: greedy(b)
        [[2, 1], [1]]
    """
    l = len(b.cartesian_factors())
    c = tuple([tuple(sorted(b[l-1]))])
    for i in range(l-2,-1,-1):
        lower = sorted([j for j in b[i] if j<=c[0][0]])
        upper = sorted([j for j in b[i] if j>c[0][0]])
        c = tuple([tuple(upper+lower)])+c
    return c

def minimaj_ordering(b):
    """
    Return minimaj ordering of `b`, where `b` is an element in decreasing order within blocks.

    EXAMPLES::

        sage: b = [[3, 1, 2], [2, 3]]
        sage: minimaj_ordering(b)
        [[3, 1, 2], [2, 3]]
        sage: b=[[2,3],[2],[1,3]]
        sage: minimaj_ordering(b)
        ((3, 2), (2,), (1, 3))
    """
    result = [tuple(sorted(b[-1]))]
    for i in range(len(b)-2,-1,-1):
        m = result[0][0]
        result = [tuple(sorted([j for j in b[i] if j>m])+sorted([j for j in b[i] if j<=m]))] + result
    return tuple(result)

def minimaj(b):
    """
    Return minimaj of `b`.

    EXAMPLES::

        sage: S = crystal_elements(2,5,3)
        sage: [minimaj(b) for b in S]
        [2, 2, 1, 1, 1, 2]
    """
    return Word(flatten(minimaj_ordering(b))).major_index()

def crystal_elements(k,ell,n):
    """
    Return list of crystal elements with `ell` letters from alphabet `\{1,2,...,n\}`
    divided into `k` blocks.

    EXAMPLES::

        sage: crystal_elements(2,3,2)
        (((2, 1), (1,)), ((1, 2), (2,)), ((1,), (1, 2)), ((2,), (1, 2)))
    """
    M = multiset_partitions(k,ell,n)
    return tuple([tuple(greedy(b)) for b in M])

def partial_sum(list):
    """
    Return partial sums of elements in `list`.

    EXAMPLES::

        sage: list = [1,3,5]
        sage: partial_sum(list)
        [0, 1, 4, 9]
    """
    result = [0]
    for i in range(len(list)):
        result += [result[-1]+list[i]]
    return result

def which_block(partial, p):
    """
    Return which block position `p` belongs to.

    EXAMPLES::

        sage: which_block([0,3,5,6],4)
        1
        sage: which_block([0,3,5,6],5)
        2
        sage: which_block([0,3,5,6],6)
        2
    """
    i=0
    while p >= partial[i+1] and i<len(partial)-2:
        i += 1
    return i

def descents(mu):
    """
    Return descents of minimaj word ``mu``.

    EXAMPLES::

        sage: mu = ((1,2,4),(4,5),(3,),(4,6,1),(2,3,1),(1,),(2,5))
        sage: descents(mu)
        [4, 7, 10]
    """
    w = flatten(mu)
    return [i for i in range(len(w)-1) if w[i]>w[i+1]]

def to_tableau(mu):
    """
    Return bijection from minimaj words to sequence of tableaux.

    EXAMPLES::

        sage: mu = ((1,2,4),(4,5),(3,),(4,6,1),(2,3,1),(1,),(2,5))
        sage: to_tableau(mu)
        [[5, 1], [3, 1], [6], [5, 4, 2], [1, 4, 3, 4, 2, 1, 2]]
    """
    b = [mu[i][0] for i in range(len(mu))]
    beginning = partial_sum([len(mu[i]) for i in range(len(mu))])
    w = flatten(mu)
    D = [0]+descents(mu)+[len(w)]
    pieces = [b]
    for i in range(len(D)-1):
        p = [w[j] for j in range(D[i]+1,D[i+1]+1) if j not in beginning]
        pieces = [p[::-1]] + pieces
    return pieces

def from_tableau(t):
    """
    Return minimaj word from sequence of tableaux.

    EXAMPLES::

        sage: all(mu == from_tableau(to_tableau(mu)) for mu in crystal_elements(3,6,3))    
        True
    """
    mu = [(i,) for i in t[-1]]
    breaks = [0]+descents(t[-1])+[len(mu)-1]
    t = [t[i][::-1] for i in range(len(t)-1)][::-1]
    for f in range(len(breaks)-1):
        for j in range(breaks[f],breaks[f+1]+1):
            mu[j] += tuple(i for i in t[f] if (mu[j][0]<i or j==breaks[f]) and (j==breaks[f+1] or i<=mu[j+1][0]))
    return tuple(mu)

def val(k,n,r):
    """
    Return val polynomials.

    """
    M = MinimajCrystal(k,n,r)
    H = [t for t in M if t.is_highest_weight()]
    Sym = SymmetricFunctions(QQ['q'])
    q = Sym.base_ring().gens()[0]
    s = Sym.schur()
    return sum((q**(minimaj(t.value))*s[[flatten(t.value).count(i) for i in range(1,r+1)]] for t in H), Sym.zero())

class MinimajCrystal(UniqueRepresentation, Parent):

    def __init__(self, k, ell, n):
        Parent.__init__(self, category = ClassicalCrystals())
        self.n = n
        self.k = k
        self.ell = ell
        self._cartan_type = CartanType(['A',n-1])
        self.module_generators = [ self(b) for b in crystal_elements(k,ell,n) ]

    def _repr_(self):
        """
        EXAMPLES::

            sage: B = MinimajCrystal(2,4,3); B
            Minimaj Crystal of type A_3 of words of length 4 into 2 blocks
        """
        return "Minimaj Crystal of type A_%s of words of length %s into %s blocks"%(self.n-1, self.ell, self.k)

    # temporary workaround while an_element is overriden by Parent
    _an_element_ = EnumeratedSets.ParentMethods._an_element_

    class Element(ElementWrapper):

        def e(self,i):
            t = to_tableau(self.value)
            B = crystals.Tableaux(['A',self.parent().n-1],shape=[1])
            T = tensor([B]*self.parent().ell)
            blocks = [len(h) for h in t]
            partial = partial_sum(blocks)
            b = T(*[B(a) for a in flatten(t)])
            if b.e(i) is None:
                return None
            b = b.e(i)
            b = [[b[a].to_tableau()[0][0] for a in range(partial[j],partial[j+1])] for j in range(len(partial)-1)]
            return self.parent()(from_tableau(b))

        def f(self,i):
            """
            Return `f_i` on ``self``.

            EXAMPLES::

                sage: B = MinimajCrystal(2,4,3)
                sage: b = B.an_element(); b
                ((2, 3, 1), (1,))
                sage: b.f(1)
                ((2, 3), (1, 2))
            """
            t = to_tableau(self.value)
            B = crystals.Tableaux(['A',self.parent().n-1],shape=[1])
            T = tensor([B]*self.parent().ell)
            blocks = [len(h) for h in t]
            partial = partial_sum(blocks)
            b = T(*[B(a) for a in flatten(t)])
            if b.f(i) is None:
                return None
            b = b.f(i)
            b = [[b[a].to_tableau()[0][0] for a in range(partial[j],partial[j+1])] for j in range(len(partial)-1)]
            return self.parent()(from_tableau(b))

[tscrim]

comment:8

Thank you.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

5069fff

rewrote for greater speed

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 8a2ded6 to 5069fff

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

af1175d

Added some new keywords for constraints on set of all ordered multiset partitions. Cleaned up documentation.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 5069fff to af1175d

[alauve]

comment:11

general query. which is better?:

-    if isinstance(n_or_OMPs, (int, Integer)):
+    if n_or_OMPs in ZZ:

The former does not catch x in the result of the following block of code.

x = 4
x = x/3
x = 6*x
type(x)
isinstance(x, (int, Integer))
x in ZZ

But perhaps there are speed reasons to choose isinstance?

[tscrim]

comment:12

The former does not include 4/2 (which is a Rational).

[alauve]

comment:13

Thanks, Anna. for review purposes, I think it's best to implement crystal structure via a new ticket. (Anyway, you take very different slices than I do in your definition of multiset_partitions, so it's not entirely straightforward for me.) Look for a new ticket and branch early next week.

Replying to @anneschilling:

Here is my code, but I have not recently checked it:

...

[anneschilling]

comment:14

Hi Aaron,

I would be happy to help to get this crystal code into sage. But you think it will be "easy" to put it on top of your code for ordered multiset partitions?

Best,

Anne

Replying to @alauve:

Thanks, Anna. for review purposes, I think it's best to implement crystal structure via a new ticket. (Anyway, you take very different slices than I do in your definition of multiset_partitions, so it's not entirely straightforward for me.) Look for a new ticket and branch early next week.

Replying to @anneschilling:

Here is my code, but I have not recently checked it:

...

[alauve]

comment:15

I wouldn't say "easy," but feasible. For my ultimate purpose (a free, graded Hopf algebra on finite subsets of NN), it would be nice to have finite dimensional slices. I couldn't think of a better way to do this than grading by total sum of entries. This makes your e and f operators not graded morphisms, but that shouldn't be an issue.

Travis suggested adding e and f operators and crystal functionality as done in shifted_primed_tableau.py specifically, ShiftedPrimedTableaux_shape.

Give me a day or two to tackle it.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from af1175d to bb57939

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

bb57939

complete overhaul, adding constraints to ordered multiset partitions

[alauve]

comment:17

Essentially starting over with this ticket...

In the use of Benkart, et al., it is natural to call e.g., OrderedMultisetPartitions([2,4], 3) or (resp.) OrderedMultisetPartitions([2,4], 3, length=2) to return the ordered multiset partitions on alphabet {2,4} with 3 total integers used (resp., and with two blocks).

In the use for a forthcoming combinatorial Hopf algebra, it is natural to call, e.g., OrderedMultisetPartitions(4) to return all ordered multiset partitions whose sum (across all blocks) is 4.

In fact, neither of these make as much sense as as the use OrderedMultisetPartitions([1,1,3,4,4,4]) or (resp.) OrderedMultisetPartitions([1,1,3,4,4,4], length=3) to return the ordered multiset partitions of the multiset {{1,1,3,4,4,4}} (resp. into two blocks).

So I have implemented all of these different uses, and added constraints as well, so one can call, e.g.,

sage: OrderedMultisetPartitions([2,3,4], 4, max_length=3)
sage: OrderedMultisetPartitions(8, alphabet=[2,4], min_length=3, max_order=4)

etc.

I still need to implement the crystal operators from Anne's comment:7.

[alauve]

comment:18

I have to say I fell down a rabbit hole and may not have worked my way out of it...

If somebody could comment on how best to handle constraints to a class (e.g., like max_length or max_block_size for SetPartitions) before I get too much farther, I'd appreciate it.

[anneschilling]

comment:19

Hi Aaron,

Perhaps it would help to see how constraints are handled in other classes? See for example Compositions or IntegerVectors. For example, you can specify various parameters for IntegerVectors

sage: IntegerVectors(7,2,min_part=2).list()
[[5, 2], [4, 3], [3, 4], [2, 5]]
sage: IntegerVectors(7,min_part=2).list()
[[7], [5, 2], [4, 3], [3, 4], [3, 2, 2], [2, 5], [2, 3, 2], [2, 2, 3]]

Best,

Anne

[alauve]

comment:20

Thanks for the IntegerVectors tip, Anne. I had looked at Compositions but didn't like how it passed things off to IntegerListsLex (and wasn't sure I understood everything happening there).

It looks like wasn't too far off on the basic structure. I should have further by end of weekend.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

4a5b2df	Moving shuffle products with overlap to combinat/shuffle.py.
4b508fa	Extending shuffle products with overlap to more general iterables.
0bd18f3	getting rid of zero
49d8364	Create the class of ordered multiset partitions.
ff02456	rewrote for greater speed
b131a25	Added some new keywords for constraints on set of all ordered multiset partitions. Cleaned up documentation.
44706f5	complete overhaul, adding constraints to ordered multiset partitions
48b5830	Merge branch 'public/combinat/implement-ordered-multiset-partitions-25148' of trac.sagemath.org:sage into public/combinat/implement-ordered-multiset-partitions-25148

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from bb57939 to 48b5830

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 48b5830 to bf0a58a

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

a0c5b82	Create the class of ordered multiset partitions.
b4dedb0	rewrote for greater speed
c9f9b02	Added some new keywords for constraints on set of all ordered multiset partitions. Cleaned up documentation.
80e04bb	complete overhaul, adding constraints to ordered multiset partitions
331d47c	Moving shuffle products with overlap to combinat/shuffle.py.
5cba60f	Create the class of ordered multiset partitions.
69c7de9	complete overhaul, adding constraints to ordered multiset partitions
bf0a58a	resolving rebase conflicts

[alauve]

comment:23

Remark: I lot of the commits above are effectively empty because I took "HEAD" each time when resolving conflicts while rebasing from current development branch. (One could safely ignore all of the above and wait for me to add e and f operators later this week before taking a look.)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

17dcd58

add crystal of ordered multiset partitions

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from bf0a58a to 17dcd58

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

e1d67d6

fixed bug in cardinality computation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 17dcd58 to e1d67d6

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from e1d67d6 to c416d02

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

c416d02

improved an_element methods

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c416d02 to cdc759b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

845fcef	cleaning up examples and tests in .soll_gt
cdc759b	deleting old code

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

2d1042b

zapping some `__iter__` bugs, though some still remain; editing examples/tests in docstrings

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from cdc759b to 2d1042b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 2d1042b to f87635b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

f87635b

added `__eq__` method for OrderedMultisetPartitions

[alauve]

comment:30

Anne, Could you have a look at this branch now? I think the crystal code is in decent shape.

Anybody, Could you explain why the following two different blocks of code behave differently?

      sage: O1 = OrderedMultisetPartitions(weight=[2,1])
      sage: list(O1) # maximum recursion depth exceded

      sage: O2 = OrderedMultisetPartitions(weight=[2,0,1])
      sage: O2.list()
      [[{1}, {1}, {3}], [{1}, {1,3}], [{1}, {3}, {1}], [{1,3}, {1}], [{3}, {1}, {1}]]
      sage: list(O2)
      [[{1}, {1}, {3}], [{1}, {1,3}], [{1}, {3}, {1}], [{1,3}, {1}], [{3}, {1}, {1}]]

[darijgr]

comment:31

I don't know the code well enough, but there is a certain code smell that might well hint at the underlying issue. Here:

+        P = OrderedMultisetPartitions(multiset)

you're using the duck-typing factory class OrderedMultisetPartitions. Are you sure your multiset will actually be understood as a multiset? Similarly, here

+            for alpha in Permutations(multiset):

you're using Permutations; the right class to use here is Permutations.mset from .permutations.

Also, I have the impression that the multiset method on your class leaks mutable internal data (namely, self._multiset which, contrary to its name, is a list).

Finally, I think "list of subsets" should be "list of nonempty subsets" everywhere in the doc, or is it not so?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

a7cf678

change type of attribute _Xlist to tuple

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from f87635b to a7cf678

[alauve]

comment:33

Replying to @darijgr:

I don't know the code well enough, but there is a certain code smell that might well hint at the underlying issue. Here:

+        P = OrderedMultisetPartitions(multiset)

you're using the duck-typing factory class OrderedMultisetPartitions. Are you sure your multiset will actually be understood as a multiset?

Pretty sure... not definite.

Similarly, here

+            for alpha in Permutations(multiset):

you're using Permutations; the right class to use here is Permutations.mset from .permutations.

Thanks, Darij

Also, I have the impression that the multiset method on your class leaks mutable internal data (namely, self._multiset which, contrary to its name, is a list).

Unfortunately, moving things to tuples instead of lists did not solve the problem. But I suppose tuples are best, so made this change.

Finally, I think "list of subsets" should be "list of nonempty subsets" everywhere in the doc, or is it not so?

Correct (and corrected).

[darijgr]

comment:34

What happens if you replace

P = OrderedMultisetPartitions(multiset)

by

P = OrderedMultisetPartitions(weight=multiset)

?

[alauve]

comment:35

calling without a keyword:

sage: P1 = OrderedMultisetPartitions([1,1,4]); P1.list()
[[{1}, {1}, {4}], [{1}, {1,4}], [{1}, {4}, {1}], [{1,4}, {1}], [{4}, {1}, {1}]]
sage: P2 = OrderedMultisetPartitions({1:2, 4:1}); P2.list()
[[{1}, {1}, {4}], [{1}, {1,4}], [{1}, {4}, {1}], [{1,4}, {1}], [{4}, {1}, {1}]]

versus calling with a keyword:

sage: P3 = OrderedMultisetPartitions(weight=[2,0,0,1]); P3.list()
[[{1}, {1}, {4}], [{1}, {1,4}], [{1}, {4}, {1}], [{1,4}, {1}], [{4}, {1}, {1}]]
sage: P4 = OrderedMultisetPartitions(weight={1:2, 4:1}); P4.list()
[[{1}, {1}, {4}], [{1}, {1,4}], [{1}, {4}, {1}], [{1,4}, {1}], [{4}, {1}, {1}]]
sage: P5 = OrderedMultisetPartitions(weight=((1,2), (4,1))); P5.list()
[[{1}, {1}, {4}], [{1}, {1,4}], [{1}, {4}, {1}], [{1,4}, {1}], [{4}, {1}, {1}]]

one exception: something that looks like a dictionary but isn't gets treated incorrectly:

sage: P6 = OrderedMultisetPartitions(((1,2), (4,1))); P6.list()
[[{(1,2)}, {(4,1)}], [{(1,2),(4,1)}], [{(4,1)}, {(1,2)}]]

(but I'm not using syntax like P6 anywhere)

[darijgr]

comment:36

I meant, what if you change the code itself and check for your recursion-depth bug again?

[alauve]

comment:37

Hmm... well something different happened (including an error!). I'll take a deeper look and report back.

[tscrim]

comment:38

Replying to @darijgr:

Here:

+        P = OrderedMultisetPartitions(multiset)

you're using the duck-typing factory class OrderedMultisetPartitions.

For the record, this is not duck-typing. This is essentially polymorphism...well, as best as Python supports it. Duck-typing is something like x.factor(), where x could be either an element of ZZ, QQ, QQ['x,y'], etc. They are completely different classes, but both quack like a duck (have a factor method), and so it is a duck (as far as that code is concerned).