# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: SageMath 8.1
#     language: ''
#     name: sagemath
# ---

# # Symmetric Superpolynomials: A Compendium of Results and Open Problems - Worksheet 
# This worksheet is intended to be used alongside the article Symmetric Superpolynomials: A Compendium of Results and Open Problems which can be obtain on ArXiv HERE(GIVELINK).
#
# This worksheet provides an introduction to the basic and more advanced features that are implemented in this library. The core of the program is based on Sage and features related to superpolynomials are provided in the following files: 
# - `sym_superfunct.py`
#     - Contains everything about superpolynomials on the basis\[spart\] representation
# - `superpartition.py`
#     - Contains everything about superpartitions
# - `superspace.py`
#     - Contains the necessary for the expansion of bases on the odd and even variables
# - `./super_cache/`
#     - Contains the pickled expression for development of different bases on the monomial superbasis.
#
# All those files can be loaded at once using the `all.py` file, this is what the following cell does, please go ahead and execute it.

#This cell has to be executed
load('all.py')
import numpy as np
from IPython.display import display, Math

# Now the library is loaded into the worksheet. As you read, make sure to execute every cell as they come by because the code of the following cell often relies on code executed in earlier cells. 
#
# The worksheet respects the structure of the compendium except for sections that are strictly about the functionalities. 
#
# Through the different sections you will get familiar with the classes and methods associated to symmetric superpolynomials. You will also see the propreties and conjectures seen in the document for yourself. Also, we tried not to hide everything in the library files and build non-trivial functions from scratch in the worksheet. These examples will help you understand how you could implement your own features. 
#
# To understand what is going on more in depth, please read the library files and the Sage documentation.

# # Superpartitions and symmetric superpolynomials

# ## Superpartitions
# ### Superpartition: the $(\Lambda^a; \Lambda^s)$ description. 
# The superpartitions are introduced in the following way. 
# #### `Superpartitions()` class
# We first introduce the complete space of superpartitions:

S = Superpartitions() # Creates an instance of the set of superpartitions
S #We print the description of the just created set. 

# Let us consider a superpartition.
#  (Note that since Python is underlying  Sage, one should never use `lambda` as variable name since it is a special function.)

my_spart = S([[3,2],[4,2,2]]) # This how you define a superpartition
my_spart # We print it

# We know that for superpartition $\Lambda$ we have
# $$
# \Lambda_1 > \Lambda_2 > \cdots > \Lambda_m \geq 0 \quad \text{and} \quad
# \Lambda_{m+1} \geq \Lambda_{m+2} \geq \cdots \geq \Lambda_{\ell} > 0
# $$
# This is enforced into the declaration of superpartitions. Trying to define a superpartition that violates these conditions will raise an error:

try: # This tries to execute the following statement
        my_bad_spart = S([[3,2,1,1],[5]]) # This generates an error
except ValueError as err: # If an error was raised, the following is executed
    print("The following error was raised: ", err)

# #### `len( a_superpartition )` and `sector()`
# Back to our superpartition, we can retrieve different information about it. 

print(my_spart)
length = len(my_spart) # Define length
print('The length of my spart is ' + str(length))
sector = my_spart.sector() 
print('My spart is in the sector ' + str(sector))

# You can access both $\Lambda^a$ and $\Lambda^s$

fermionic_part = my_spart[0]
print(fermionic_part)
bosonic_part = my_spart[1]
print(bosonic_part)

# #### `Superpartitions(n, m)`
# Previously we have defined the whole space of superpartitions. Now let say that we want to deal with only one sector. Then, we just call the `Superpartitions` with `n` (the number of boxes, or bosonic degree) and `m` (the number of circles, or fermionic degree) as arguments:

spart_nm = Superpartitions(12, 3)
spart_nm

# You can retrieve the elements of a set as a list, let's take a smaller set not to clutter the doccument. 

myset = Superpartitions(5, 2)
print(myset)
mylist = list(myset)
print(mylist)

# Note that order in the list is deterministic but somewhat arbitrary. 
#
# The set can also be used as an iterator:

an_other_set = Superpartitions(4,2)
print(an_other_set)
# We use it here as an iterator: 
for some_superpartition in an_other_set:
    print(some_superpartition)

# ### The $\Lambda^\ast, \Lambda^\circledast$ description
# #### `partition_pair()`
# We can obtain this description from a superpartition:

spart1 = S([[6,5,2,0],[4,3,3,2,1]])
print(spart1)
spart1_pair = spart1.partition_pair()
print(spart1_pair)

# #### `Superpartitions.partition_pair_to_spart( a_partition_pair )`
# We can also revert this operation (even though this operation is not really intended for end-user):

Superpartitions.partition_pair_to_spart(spart1_pair)

# ### Superdiagram
# #### `terminal_diagram()`
#
# We can obtain the representation of the superdiagram directly from a superpartition. 

S = Superpartitions()
a_spart = S([[7,5,2,0],[8,8,5,4,2,2]])
print(a_spart)
a_spart.terminal_diagram()

# ### Conjugate superpartition
# #### `conjugate()` method
# We can easily obtain the conjugate of a superpartition. 

S = Superpartitions()
a_spart = S([[9,7,5,2,0],[8,8,5,3,2,2,2,1]])
print(a_spart)
a_spart.terminal_diagram()
print('And now the conjugate superpartition')
a_spart_conjugate = a_spart.conjugate()
a_spart_conjugate.terminal_diagram()

# ### Dominance order
#
# First, let's take three superpartitions with 10 boxes and 3 circles. 

# +
S = Superpartitions(10,3)
spart_list = list(S)
spart1 = spart_list[22]
spart2 = spart_list[33]
spart3 = spart_list[35]

print('spart1')
print(spart1)
spart1.terminal_diagram()

print('spart2')
print(spart2)
spart2.terminal_diagram()

print('spart3')
print(spart3)
spart3.terminal_diagram()

# -

# #### `compare_dominance()` method
# Let us now recover their ordering relations:

print('Compare spart1 and spart2')
print(S.compare_dominance(spart1, spart2))
print('Compare spart1 and spart3')
print(S.compare_dominance(spart1, spart3))
print('Compare spart2 and spart3')
print(S.compare_dominance(spart2, spart3))

# The `compare_dominance()` method returns one of the following strings: 
#     `"=="`, 
#     `">"`, 
#     `"<"` or
#     `"Non-comparable"`. We illustrate this with an example:
#
# Let us consider the following set of superpartitions of sector $(4\,|\,1)$:

Sparts = Superpartitions(4,1)
print(Sparts)
print(list(Sparts))
for k in Sparts:
    print(k)

# We will compare the superpartitions $(2;2)$ to every superpartitions of the previous set:

# +
Sparts = Superpartitions()
our_spart = Sparts([[2],[2]]) # We define (2; 2)
Sparts_4_1 = Superpartitions(4,1)
our_spart.terminal_diagram() # Print the diagram

for an_other_spart in Sparts_4_1:
    #We store the result of the operation
    comparision = Sparts.compare_dominance(our_spart, an_other_spart)
    if comparision == ">":
        print(str(our_spart) + " dominates " + str(an_other_spart))
    elif comparision == "==":
        print(str(our_spart) + " is exactly " + str(an_other_spart))
    elif comparision == "<":
        print(str(our_spart) + " is dominated by " + str(an_other_spart))
    elif comparision == "Non-comparable":
        print(str(our_spart) + " and " + str(an_other_spart) + " are non-comparable")
# -

# Take note that the dominance ordering overloads the `>`, `<`, `>=`, ... operators. For instance

print( Sparts([[2],[2]]) <= Sparts([[3],[1]]))
print( Sparts([[2],[2]]) != Sparts([[2],[1,1]]))

# Hence, the previous example can be implemented in a much more readeable fashion: 

# +
Sparts = Superpartitions(4,1)
our_spart = Sparts([[2],[2]]) # We define (2; 2)

for an_other_spart in list(Sparts):
    if (our_spart > an_other_spart):
        print(str(our_spart) + " dominates " + str(an_other_spart))
    elif our_spart == an_other_spart:
        print(str(our_spart) + " is exactly " + str(an_other_spart))
    elif our_spart < an_other_spart:
        print(str(our_spart) + " is dominated by " + str(an_other_spart))
    else:
        print(str(our_spart) + " and " + str(an_other_spart) + " are non-comparable")
# -

# #### `get_smaller_sparts()` method
# allows one to obtain the list of every superpartitions that is dominated by the one the method is called from. For instance,

Sparts = Superpartitions()
my_spart = Sparts([[2,0],[1,1]])
print("The superpartition")
my_spart.terminal_diagram()
print("dominates the following")
dominated_sparts = my_spart.get_smaller_sparts()
for spart in dominated_sparts:
    spart.terminal_diagram()

# #### `Superpartitions.give_greatest_spart( list_of_spart )`
# this function allows one to extract the greatest superpartition out of the list:

Sparts_4_2 = list(Superpartitions(4,2))
some_list = [Sparts_4_2[x] for x in [2,4,6,8]]
print(some_list)
print("The greatest element in the list is")
greatest = Superpartitions.give_greatest_spart(some_list)
print(greatest)

# #### `Superpartitions.sort_by_dominance( list_of_spart )`
# allows to sort a list of superpartitions according to dominance. Non-comparable superpartitions are sorted in lexicographic order.

print(some_list)
sorted_list = Superpartitions.sort_by_dominance(some_list)
print("once sorted :")
print(sorted_list)

# Note that other methods are provided for evaluation of certain quantities on the diagram, such as hook length, arm length... These quantities will be introduced in the section about the normalization of Jack superpolynomials. 

# # Superspace
# We introduce here the space of supercommutative variables. Since we have Python under the hood, the index of the variables start at 0 to avoid confusion. So to write expressions in superspace, one have to create an instance of superspace:

# The syntax is superspace(N) to create a space of variables
# with x_i, theta_i for 0 <= i < N
ss = superspace(6)
ss

# Now we have created a superspace in 12 variables, the even (or bosonic) variables $x_0, ..., x_5$ and the odd (or fermionic) variables $\theta_0, ..., \theta_5$. Now we would like to be able to write expressions in these variables.
#
# Now this is a bit technical, but Sage does not have a built-in way of dealing with supercommutative variables, the superspace class serves as a interface to Singular which is doing the computation behind the scene. So when we defined ss as a superspace with $N=6$, Singular was started in the background. Also, as a side effect, the names of the variables were injected in the Sage session:

theta_1*theta_0

(theta_0*x_0 + theta_1*x_1)*theta_0

# One must be careful with operations one wants to perform because these are not sage objects:

print(type(theta_0*x_2))

# So if you are used to sage, many operations might not work. 
#
# Back to the main matter, we can look at what is the current ring in Singular. Notice that there is an extra fermionic variable dtheta, this variable is introduced as a work around for the differentiation on odd variables and this variable should not be called directly:

singular.current_ring()

# It seems that Sage parses some expressions through Singular and this could reset the `current_ring()`. So if you have problems with some calculation in superspace, make sure to check if the supercommutative relations are still there. To avoid this problem, make all computations related to superspace right after the `ss = superspace(N)` declaration.
#
#
# Differentiation is implemented with simple syntax `superspace.diff(expr, variable)`:

my_expr = theta_0*theta_1*theta_3*x_1*x_2^2*x_3^4
print(ss.diff(my_expr, theta_1))
print(ss.diff(my_expr, theta_3))
print(ss.diff(my_expr, x_3))

# The superspace instance comes with it's permutation group:

S6 = ss.SN()
S6

# These permutations can be applied with a special method `superspace.apply_permutation(expr, permutation)`, we illustrate this by printing the action of the permutations of length 12 on the following expression:

my_expr = theta_0*theta_1*x_0^2*x_5^3
for a_permutation in S6:
    if a_permutation.length() == 12:
        new_expr = ss.apply_permuation(my_expr, a_permutation)
        print(str(a_permutation) + " ---> " + str(new_expr))

# From an expression, we can extract the sector to which it belongs. Note that this only works for homogeneous superpolynomials:

my_expr = theta_0*theta_1*theta_2*theta_3*(x_1^2*x_2 - x_1*x_2^2)
the_sector = ss.get_sector(my_expr)
print(the_sector)

# We can also determine if a given polynomial is symmetric with the `superspace.is_symmetric(expr)` method. We also use this occasion to introduce a method that returns a list with the variables so that one can use loops to create polynomials instead of writting everything by hand: 

X = ss.x_vars()
Theta = ss.theta_vars()
print(X)
print(Theta)

# We create here a symmetric polynomial
my_sym_pol = sum(Theta[k]*X[k]^2 for k in range(6))
print(my_sym_pol)
print(ss.is_symmetric(my_sym_pol))

# The following is not symmetric since our superspace is defined over N=6
my_nonsym_pol = theta_1*x_1 + theta_2*x_2 + theta_3*x_3
print(ss.is_symmetric(my_nonsym_pol))

# There is a method that symmetrize your polynomial by applying every permutation of SN, it really is every permutation, not just distinct one, so use with care since you can easily crash Sage and Singular with $N>9$. For these high variables cases, special care must be taken. But here with $N=6$ it's all fine:

symmed = ss.symmetrize(my_nonsym_pol)
print(symmed)
print(ss.is_symmetric(symmed))

# We see here that a lot of permutations were computed for nothing. For cases where speed is crucial, we implemented a distinct permutation generator for the problem at hand. In `superspace.py`, if you want to take a look, see `superspace.monomial(self, spart)` where every step is well documented.
#
# For more information on how to handle these objetcs or create your own methods in superspace see:
#
# http://doc.sagemath.org/html/en/reference/interfaces/sage/interfaces/singular.html?highlight=singular#module-sage.interfaces.singular
#
# https://www.singular.uni-kl.de/index.php/singular-manual.html
#
# Reading `superspace.py` might help you with the philosphy adopted here. 

# ## Symmetric superpolynomials
# We now introduce the space of symmetric superpolynomials. We first define the space:

Sym = SymSuperfunctionsAlgebra(QQ)
Sym

# We could have declared it with another coefficient ring, such as `ZZ`, but many things would fail because  their connection with the powersum basis (which is not a basis over  $\mathbb{Z}$).
#

# ### The monomial basis
# Now this allows us to define the basis over the ring. For instance, one can introduce the monomial basis like so:

m = Sym.Monomial()
m

# We can now write an expression in the monomial basis:

Sparts = Superpartitions()
spart1 = Sparts([[3,2],[1,1]])
spart2 = Sparts([[3,0],[2,2]])
some_expr = 3*m(spart1) + 5/2*m(spart2)
print(some_expr)

# Now you see that we have defined two superpartition-type object to create our element of basis. When doing large scripts this is the way to go ( that is `m( superpartition )`). But when experimenting directly, it is rather cumbersome to define superpartitions, that's why abuse of notation is allowed. One could equally write the following:

easier_to_write = 3*m[[3,2],[1,1]] + 5/2*m[[3,0],[2,2]]
print(easier_to_write)

# Incidentally, this also allows copy-and-paste of printed results as valid sage expressions. Now back to the monomial basis, the product of two monomials is implemented:

my_pol = 3*m[[1,0],[3,2]] + 2/3*m[[2,0],[2,1,1]]
print(my_pol)
print('Result of left multiplication by m[[1],[]]')
my_other_pol = m[[1],[]]*my_pol
print(my_other_pol)

# #### Expansion on the variables
# Any expression can easily be expressed in superspace. One should take note that the number of variables must be at least the length of the longest superpartition in the expression. 

# We introduce the space of symmetric functions
sym = SymSuperfunctionsAlgebra(QQ)
m = sym.Monomial()
# We introduce the superspace in 4 variables
ss = superspace(4)

# We make an arbitrary expression with random integers [0,15] as coefficient
# Here we use the numpy package for the random integer generator
rand_integer = np.random.randint
rand_expr = rand_integer(15)*m[[2,0],[3]] + rand_integer(15)*m[[1,0],[2,2]]
print(rand_expr)
# Now we expand over the superspace ss
expr_x = rand_expr.expand(ss)
print(expr_x)

# The monomial expression can be retrieved (this works with any basis) with `basis.from_polynomial(poly, superspace)`

back_in_m = m.from_polynomial(expr_x, ss)
print(back_in_m)
print(back_in_m == rand_expr)

# Here we recall that a monomial with length $< N$ will yield 0 when expanded, the `expand` method will warn you about that:

# +
sym = SymSuperfunctionsAlgebra(QQ)
m = sym.Monomial()
ss = superspace(3)

my_expr = m[[1,0],[2]] + 3*m[[1,0],[1,1]]
my_expr_x = my_expr.expand(ss)
mono_from_x = m.from_polynomial(my_expr_x, ss)

print("")
print(my_expr)
print("")
print(mono_from_x)
# -

# ### Mutliplicative bases
# We introduce the other three multiplicative basis respectively the power-sum, elementary and completely homogeneous bases:

p = Sym.Powersum()
e = Sym.Elementary()
h = Sym.Homogeneous()
print(p); print(e); print(h)

# To define elements of this basis, it works the same as with monomials. You can either declare `p( a_superpartition )` or write directely `p[[...], [...]]`

pol_in_p = 1/4*p[[1,0],[2,2]] + 3*p[[3,2], []]
pol_in_e = 1/4*e[[1,0],[2,2]] + 3*e[[3,2], []]
pol_in_h = 1/4*h[[1,0],[2,2]] + 3*h[[3,2], []]
print(pol_in_p)
print(pol_in_e)
print(pol_in_h)

# The proper multiplication rules work out of the box:

print(p[[2],[]]*pol_in_p)
print(e[[2],[]]*pol_in_e)
print(h[[2],[]]*pol_in_h)

# It has to be noted that the basis indexed by the empty partition is considered to be unity (which differs from `p[[],[]] <-> 0` convention that we have in many formulas.)

print(p[[],[]]*pol_in_p)
print(e[[],[]]*pol_in_e)
print(h[[],[]]*pol_in_h)

# As with the monomial case, one can expand expressions in variables:

ss = superspace(3)
p_pol = p[[1,0],[2]]
p_pol_x = p_pol.expand(ss)
print(p_pol_x)
print("")
# and convert polynomials back on any basis
print(p.from_polynomial(p_pol_x, ss))
print(e.from_polynomial(p_pol_x, ss))

# ### Change of basis
# Obtaining the expression on another basis is achieved by the following syntax:
#
# `new_basis( expr )`
#
# He are some examples :

my_expr = m[[2,0],[2]] + 1/2*m[[2,1],[1]]
print(my_expr)
my_expr_p = p(my_expr)
print(my_expr_p)

# Comparision between expressions on different basis is a valid operation:

print(my_expr == my_expr_p)

# The implementation of change of basis is not lighting fast, but it seems sufficiently fast for most research purposes. 
#
# Now the fact that we can compare two expressions expressed on different bases comes from the coercion model of Sage. The coercion is some kind of automatic conversion that allows one to make operations that make sense mathematically without having to explicitely program what is meant by this operation. For instance, you can multiply two different basis:

p[[1],[1,1]]*m[[3],[2,1]]

# So for example, we can test relations quite explicitely. Let us verify
# $$
# 5 h_5 = \sum_{r=1}^5 p_r h_{5-r}
# $$

# We first evaluate the right hand side (in python range(a,b) = [a,a+1, ..., b-1])
rhs = sum(p[[],[r]]*h[[],[5-r]] for r in range(1,6))
# We convert the sum to the basis h and print the result
print(h(rhs))

# Which was expected. We now evaluate
# $$
# 6\tilde{h}_5 = \sum_{r=0}^5 p_r \tilde{h}_{n-r} + (r+1)\tilde{p}_rh_{n-r}
# $$
# Note that we split the sum in two since the program evaluates `p[[],[]]` as the identity instead of 0. 

sum1 = sum(p[[],[r]]*h[[5-r],[]] for r in range(1,6))  
sum2 = sum((r+1)*p[[r],[]]*h[[],[5-r]] for r in range(0,6))
print(h(sum1+sum2))

# as expected. Finally, we evaluate the following
# $$
# (6+1)\tilde{e}_6 = \sum_{r=0}^n (-1)^{r+1}\left[ p_r\tilde{e}_{6-r}
# - (r+1)\tilde{p}_r e_{6-r} \right]
# $$

sum1 = sum((-1)**(r+1)*p[[],[r]]*e[[6-r],[]] for r in range(1,7))  
sum2 = sum((-1)**(r+1)*(r+1)*p[[r],[]]*e[[],[6-r]] for r in range(0,7))
print(e(sum1-sum2))

# which is again the expected result.  

# ### Ring homomorphism I
# The $\hat{\omega}$ homomorphism is defined on every basis. 
# You can obtain the result of the involution with `expr.omega()`, it will return the effect of the homomorphism on the same basis you started from. Here is an exemple:

print('We have the following superpolynomial')
print(pol_in_e)
print('We apply the omega homomoprhism')
print(pol_in_e.omega())
print('We then convert the result to the homogeneous basis')
print(h(pol_in_e.omega()))

# The mechanics of `expr.omega()` converts the expression to the powersum basis, apply omega and converts back to the starting basis. So the previous operation would be faster if one converts `pol_in_e` to powersums before applying omega, then convert the result to the homogeneous basis.
#
# The operation is well defined on expression in any basis here is an example on the powersum basis:

an_expr = p[[2,0],[3,2]] + 1/2*p[[1,0],[3,2,1]]
print(an_expr)
print('Result of omega:')
an_expr.omega()

# ### Scalar product
# We now introduce the `expr.scalar_product(other_expr)` method. Here is an example of usage: 

Sparts_7_3 = Superpartitions(7,3)
sparts = list(Sparts_7_3)
print(len(sparts))
expr1 = p(sparts[0]) + 2*p(sparts[3]) + 5/2*p(sparts[6])
expr2 = p(sparts[0]) + 33*p(sparts[3]) + 13*p(sparts[5])
print(expr1)
print(expr2)
print('We now evaluate the scalar product between those two')
expr1.scalar_product(expr2)

# The method is well defined for any two basis. So for instance you can do the following:

expr_m = m[[3,0],[1,1]] + m[[3,1],[1]]+ m[[1,0],[1,1,1,1]]
expr_p = p[[5,0],[]] + p[[1,0],[2,1,1]] + 22*p[[1,0],[1,1,1,1]]
expr_m.scalar_product(expr_p)

# We add examples illustrating some properties. First, the monomial and the homogeneous are dual with respect to the scalar product (in superspace duality is defined up to a sign).

# We generate the set of superpartitions with 8 boxes and 3 circles
spart_set = Superpartitions(8,3)
# We take the scalar product between m_\Lambda and h_\Lambda for 
# every \Lambda in that set. The scalar product are store in a list.
scals = [m(spart).scalar_product(h(spart)) for spart in spart_set]
print(scals)

# It is known that the elementary basis is dual to the forgotten basis. 
# Even though the later is not explicitely programmed, we can test that property by assuming $f_\Lambda = \hat{\omega}(m_\Lambda)$. 
# Note that the following snipet of code might take a minute to execute. 
# This is because, for every partition, the monomial has to be converted to powersums in order to apply the $\hat{\omega}$ involution and then converted back into monomials and back into powersum to evaluate the scalar product with $e_\Lambda$ which itself has to be converted to powersum in the process. 
# If the forgotten basis is implemented in the future, this whole process should be much faster. 
# For now though, the implementation of the forgotten basis is not a priority. 

# Generate some set
Sparts = Superpartitions(4,2)
scals = [m(spart).omega().scalar_product(e(spart)) for spart in spart_set]
print(scals)

# ### The kernel and the Cauchy formula
# Here we do not introduce the Kernel since its usefulness lies more in analytical calculation than in a computer algebra system. However we verify properties shown in the compendium.
#
# We have
#
# $$
# \tilde{h}_n = \sum_{\Lambda \vdash (n|1)} z_\Lambda^{-1} p_\Lambda
# $$
#
# Let us verify this for $n=5$

sym = SymSuperfunctionsAlgebra(QQ.fraction_field())
p = sym.Powersum()
h = sym.Homogeneous()
sparts = Superpartitions(5,1)
rhs = sum(1/(spart.z_lambda())*p(spart) for spart in sparts)
print(rhs)
print("\n Converted to homogeneous")
print(h(rhs))

# Let us also verify 
# $$
# \tilde{e}_5 = \sum_{\Lambda \vdash (5|1)} z_\Lambda^{-1} \omega_\Lambda p_\Lambda
# $$

e = sym.Elementary()
rhs = sum(1/(spart.z_lambda())*p(spart).omega() for spart in sparts)
print(rhs)
print("\n As elementary")
print(e(rhs))

# ### The one parameter space of symmetric superfunctions
# We now want to define our superpolynomial ring over rational polynomials in one variable, namely $\alpha$. To do so, we have to introduce the ring like this:

# We first define the coefficient ring
# Note that we add .fracion_field() so that it is allowed to 
# write expressions such as (alpha+1)/(3*alpha+2)
QQalpha = QQ['alpha'].fraction_field()
print(QQalpha)
# We can now introduce the symmetric functions
Sym_alpha = SymSuperfunctionsAlgebra(QQalpha)
print(Sym_alpha)
# We inject the basis definitions
Sym_alpha.inject_shorthands()

# From there, it is now possible to define polynomials with coefficients including $\alpha$. Even though the ring knows the existence of the symbol $\alpha$, 
# we must assign the symbol to a variable to directly write expressions with it:

alpha = QQalpha.gen()
# Now the variable alpha is associated with the symbol alpha
print(alpha)
# And so we can write expressions out of it
print(alpha**2 + 2*alpha)

# And so, we are now able to introduce superpolynomials with coeffients as functions of $\alpha$: 

my_alpha_pol = m[[2,0],[1,1]] + 3*alpha/(2*alpha+2)*m[[1,0],[2,1]]
print(my_alpha_pol)

# ### One paremeter deformation of the scalar product
# Now that we have a proper ring, we can introduce a new scalar product. 

my_expr = p[[2,0],[3,1]] + p[[1,0],[3,1,1]]
my_expr.scalar_alpha(my_expr)

# Note that again, this is well defined for any classical basis:

my_mono_expr = (2*(alpha + 1)*m[[2,0],[1,1]] +
                alpha**2*m[[1,0],[3]])
my_mono_expr.scalar_alpha(my_mono_expr)

# ### The $\alpha$-deformed homogeneous basis
# We now introduce the $g_\Lambda$ basis, it is introduced as `galpha`. 
#
# So let us define, the space

QQa = QQ['alpha'].fraction_field()
Syma = SymSuperfunctionsAlgebra(QQa)
Syma.inject_shorthands()

# We now define the `galpha` basis

galpha = Syma.Galpha()
galpha

# It works just like any other basis:

MyPol = galpha[[4],[1]]
print(MyPol)
MyPol_in_m = m(MyPol)
print(MyPol_in_m)

# Now we see here that the result of the conversion to the monomial basis looks rather messy. Sage does not provide a built-in way of factoring the coefficient of monomials, but one has been implemented. We easily obtain a more readable expression by typing `expr.collect()`: 

MyPol_in_m.collect()

# which is a little bit better. It not quite what we are looking for yet, but I will make progress on the presention front later. 
#
#
# Now, let us illustrate the orthonormality property: 
# $$
# < m_\Lambda | g_\Omega >_\alpha = (-1)^{\binom{m}{2}}\delta_{\Lambda \Omega}
# $$
#
# which is shown here for the sector $n=5, m=2$: 

sparts = Superpartitions(5,2)
scals = [m(spart).scalar_alpha(galpha(spart)) for spart in sparts]
print(scals)

# ### Ring homomorphism II
# We now introduce the $\alpha$-deformation of the endomorphism $\hat{\omega}$, namely $\hat{\omega}_\alpha$. To call it, you simply use the following formulation `expr.omega_alpha()`: 

my_expr = p[[1,0],[2,1,1,1]] + 32*p[[2,1],[3]]
omega_expr = my_expr.omega_alpha()
print(omega_expr)

# To apply the operation, though, you have to make sure that you are indeed working with bases defined from the right coefficient ring, otherwise it will raise an error. 

bad_p = SymSuperfunctionsAlgebra(QQ).Powersum()
an_expr = bad_p[[1,0],[3]]
# We do the error handling here:
try:
    # It will run this 
    an_expr.omega_alpha()
    # and if it produces an error will execute instructions in except: 
except:
    print("An error has been raised")

# Now, since the $g_\Lambda$ basis is the $\alpha$-deformation of the homogeneous basis, we have the following property
# $$\hat{\omega}_\alpha(g_\Lambda) = e_\Lambda$$
# As can be seen here: 

sparts = Superpartitions(4,2)
for spart in sparts:
    expr = galpha(spart)
    print(expr)
    print("We apply omega_alpha")
    om_expr = expr.omega_alpha()
    om_expr_in_e = e(om_expr)
    print(om_expr_in_e)
    print('')

# and so the test is successful. 

# # Jack superpolynomials
# ## Combinatorial characterization
# Given the $\alpha$-scalar product and the dominance ordering on superpartitions, one can define the Jack superpolynomials. 
#
# We now introduce the implementation of those superfunctions. 

QQa = FractionField(QQ['alpha'])
alpha = QQa.gen()
Sym = SymSuperfunctionsAlgebra(QQa)
Sym.inject_shorthands()
Palpha = Sym.Jack()

# This basis works like any other basis

my_expr = Palpha[[2,0],[3,2]] + Palpha[[3,0],[2,2]]
print(my_expr)

# We know from the combinatorial definition that these superpolynomials are triangular on the monomial basis. We illustrate this fact by giving the expansion of the `Palpha`'s on the monomial basis for the sector $(5,3)$

for a_spart in Superpartitions(5,3):
    print("The Jack for " + str(a_spart) + " is ")
    print(m(Palpha(a_spart)))
    print('') #Spacing of output

# We also know that those polynomials are orthogonal:

sparts_set = Superpartitions(6,3)
sparts_pairs = Combinations(sparts_set, 2)
for a_pair in sparts_pairs:
    print(a_pair)
    spart1, spart2 = a_pair
    scalar_prod = Palpha(spart1).scalar_alpha(Palpha(spart2))
    print(scalar_prod)
    print('')

# ## Eigenfunction characterization

# We now illustrate how one can obtain the Jack polynomial from the Eigenvalue problem. We provide both the operators $D$ and $\Delta$ with the methods `superspace.Djack(expr)` and `superspace.DeltaJack(expr)`.

# +
sparts = Superpartitions()
dom_spart = sparts([[2,1],[2]])
smaller = dom_spart.get_smaller_sparts()

# Now we need a coefficient ring that will accept indeterminates
#   The parameter     The indeterminates         As many as #sparts
#         \             /                           /
params = ['alpha'] + ['b' + str(k) for k in range(len(smaller))]
coef_ring = PolynomialRing(QQ, params).fraction_field()

sym = SymSuperfunctionsAlgebra(coef_ring) 
# The ak variables 
indets = coef_ring.gens()[1:]
m = sym.Monomial()
# We get the generic expression of the Jack polynomial
expr = m(dom_spart) + sum(indets[k]*m(smaller[k]) for k in range(len(smaller)))
print(expr)
print("\n")

# We need this expression in the variables
N = max([len(spart) for spart in smaller])
ss = superspace(N, parameters = params)
expr_x = expr.expand(ss)

# We apply the two operators
# Notice that we normalize the expression so that the leading coefficient
# remains unity
Dexpr_x = ss.Djack(expr_x).normalize()
Dexpr = m.from_polynomial(Dexpr_x, ss)
Deltaexpr_x = ss.DeltaJack(expr_x).normalize()
Deltaexpr = m.from_polynomial(Deltaexpr_x, ss)
print('Result of D')
print(Dexpr)
print("\nResult of Delta")
print(Deltaexpr)
# -

# We are looking for coefficients `a0, a1, ...` such that `expr` is an eigenfunction of the two operators. This means, since we normalized the result of the action of the operators, that `expr - Operator(expr).normalize()` should be zero. By the linear independence of every monomial, each coefficient must be zero. That is how we obtain the folowing system of equations.

# +
eqn_system1 = (expr - Dexpr).monomial_coefficients().values()
eqn_system2 = (expr - Deltaexpr).monomial_coefficients().values()
eqn_system = eqn_system1 + eqn_system2

# Solving only works on the symbolic ring
# So we convert every equation to SR
SR_system = [SR(x) for x in eqn_system]
# And the indeterminates too
SR_indets = [SR(x) for x in indets]
# Now we solve. The [0] at the end is there because
# solve returns a list of all solutions (even if there is only one).
solns = solve(SR_system, SR_indets, solution_dict=True)[0]
print(solns)
print("\n")
# We convert these to the coef_ring
sub_dict = {coef_ring(x):coef_ring(solns[x]) for x in solns}

# We thus obtain
final_jack = expr.subs_coeff(sub_dict)
print(final_jack)
Palpha = sym.Jack()
print("We check if this is the right expression by comparing with the polynomial obtained from the scalar product")
print(Palpha(dom_spart) == final_jack)
# -

# We can also easily verify the eigenvalues

super_init()
Palpha = Sym.Jack()
jack = Palpha[[2,0],[2,1]]
ss = superspace(6)
jack_x = jack.expand(ss)
DJack = Palpha.from_polynomial(ss.Djack(jack_x), ss)
DeltaJack = Palpha.from_polynomial(ss.DeltaJack(jack_x), ss)
display(Math(r"\epsilon_\Lambda:" + latex(DJack)))
display(Math(r"\tilde{\epsilon}_\Lambda:" + latex(DeltaJack)))

# ## The symmetrization construction
# We now illustrate the construction of the Jack superpolynomials through the symmetrization of the non-symmetric Jack polynomial in superspace. Since the nsJacks are not implemented in Sage, we take the Macdonald non-symmetric polynomial and evaluate the limit. This limiting process is better explained in the Macdonald section. 
#
# Starting from the following definition
#
# $$
#  P^{(\alpha)}_\Lambda \propto \sum_{\omega \in S_N} \mathcal{K}_\omega \theta_0 \cdots \theta_{m-1} E_{\Lambda^R}(x)
# $$
#
# We are going to generate a Jack superpolynomial. We will use $\Lambda = (1,0;2)$ which implies a composition $\eta = \Lambda^R = (0,1,0,2)$ (We add zeroes so that the composition is of length $N$, where $N$ has to be at least equal to the length of the longest superpartition of the sector. 
#
# To obtain the polynomial, we are going to:
# 1. Initialize the space of symmetric superfunctions and define the Palpha basis.
# 2. Import the non-symmetric Macdonald polynomial and obtain a variable representation for our composition.
# 3. Take the limit case $q=t^\alpha$ $t\rightarrow 1$ to obtain the non-symmetric Jack polynomial.
# 4. Convert the resulting expression to a singular expression.
# 5. Mutliply by $\theta_0 \theta_1$.
# 6. Symmetrize the result and normalize, we thus obtain the Jack superpolynomial in terms of variables.
# 7. Convert the variable expression to the ring of symmetric function on the Palpha basis.
#
# Note that the following procedure is quite slow because of the limiting process. Implementing the non-symmetric Jack polynomial directly could speed it up significantly. We use this method here only for illustration purpose. 

super_init()
Palpha = Sym.Jack()

# +
from sage.combinat.sf.ns_macdonald import E
N = 4
mm = 2
composition = [0,1,0,2]
q, t, alpha = QQqta.gens()

# Now the non-symmetric Jack polynomial is not implemented in Sage
# But the non-symmetric Macdonald is, so we will obtain the nsJack
# as the proper limit. We have to reverse the composition and use
# a base permutation (pi = perm) because of Sage conventions not
# perfectly fitting with ours.
perm = range(1,N+1)
perm.reverse()
perm = Permutation(perm)
composition.reverse()
nsMacdo = E(composition, q, t, pi=perm)

# Now we take the limit, this process will be explained more in depth
# in the section on Macdonald polynomials
q_sr = SR(q); t_sr = SR(t); alpha_sr = SR(alpha)
assume(alpha_sr, 'integer')
srMac = factor(SR(nsMacdo).subs({q_sr:t_sr**alpha_sr}))
nsJacksr = srMac.limit(t=1)

# This polynomial is not in the right ring, so we convert it
# to a singular object. Singular accepts strings as expressions
# So we take the expr of the nsJack in the variables x0,x1, ...
# convert it to a string, change 'xk' --> 'x_k' in the string
# expression and parse the string through singular.
ss = superspace(N)  # We create the superspace
string_list = str(nsJacksr).split('x')
singular_str = 'x_'.join(string_list)
nsJack_singular = singular(singular_str).normalize()

# We now need to multiply this by the theta variables
theta = ss.theta_vars()
super_nsJack = prod(theta[:mm])*nsJack_singular

# We symmetrize the result and normalize
superJack = ss.symmetrize(super_nsJack).normalize()

# We convert this back in the Macdonald basis
superJackPalpha = Palpha.from_polynomial(superJack, ss)
print(superJackPalpha)
# -

# Which illustrates the definition from the non-symmetric polynomials. 

# ## Another scalar product
# We do not implement the analytical scalar product since it is computationally expensive and it's utility is very limited within a computer algebra system.

# ## Further properties of Jack superpolynomials
#
# ### Normalization
# We use the normalization of the Jack superpolynomials as an opportunity to introduce some methods exctracting combinatorial data off the diagram. 
#
# #### Combinatorial data
# Let us first introduce the coordinates of the cells

# +
sparts = Superpartitions()
my_spart = sparts([[7,4,2,0],[5,2,2,1]])
my_spart.terminal_diagram()

# To obtain the location of every box
coords = my_spart.cells()
print(coords)
# -

# As can be seen, the coordinates are returned as $(i,j)$ where the top left corner is $(1,1)$.
#
# Now we have two kind of cells, the bosonic and fermionic cells.
#
#
# The fermionic cells are those that have a both a circle to their right and a circle below them:

print(my_spart.fermionic_cells())

# The bosonic cells are those that are not fermionic:

print(my_spart.bosonic_cells())

# The arm length and leg length are not defined directly on the superpartitions. They are only defined on standard partitions. So these quantities are calculated either on $\Lambda^*$ or $\Lambda^{\circledast}$. Keeping the same superpartition as before:
#
# `spart.star().arm_length(i, j)`
#
# `spart.circle_star().leg_length(i, j)`

my_spart.star().terminal_diagram()
leg32 = my_spart.star().leg_length(3,2)
arm32 = my_spart.star().arm_length(3,2)
print("The leg length of (3,2) is " + str(leg32))
print("The arm length of (3,2) is " + str(arm32))

my_spart.circle_star().terminal_diagram()
legtilde32 = my_spart.circle_star().leg_length(3,2)
armtilde32 = my_spart.circle_star().arm_length(3,2)
print("The leg length of (3,2) is " + str(legtilde32))
print("The arm length of (3,2) is " + str(armtilde32))

# Therefore, we can introduce hook lengths
#
# `spart.upper_hook_length(i, j, parameter)`
#
# `spart.lower_hook_length(i, j, parameter)`

QQa = FractionField(QQ['alpha'])
alpha = QQa.gen()
print(my_spart.upper_hook_length(3, 2, alpha))
print(my_spart.lower_hook_length(3, 2, alpha))

# All the functions illustrated here allows one to program the formula for the normalization of Jack superpolynomials, which is implemented as `instance_of Jack_Basis.calc_norm( spart )`

Sym_a = SymSuperfunctionsAlgebra(QQa)
Palpha = Sym_a.Jack()
spart = Superpartitions()([[2,0],[3]])
my_jack = Palpha(spart)
# We obtain the scalar product from the actual operation
scal = my_jack.scalar_alpha(my_jack)
display(Math(r"<P^{\alpha}_{(2,0;3)} | P^\alpha_{(2,0;3)}> = " + latex(factor(SR(scal)))))
# We obtain the norm by the formula
norm_calc = Palpha.calc_norm(spart)
display(Math(r"||P^{\alpha}_{(2,0;3)}||^2 = " + latex(factor(SR(norm_calc)))))

# ### Evaluation

# The evaluation formula is already implemented with the method `expr.evaluation(N)`. We give an example. 
#
# First we introduce the space of symmetric functions and the Jack basis, note also that a handy function `super_init()` is provided to introduce the coefficient ring, the ring of symmetric superfunction and all the classical bases as well as an instance of `Superpartitions`:

super_init()
Palpha = Sym.Jack()

# Now we generate a Jack superpolynomial and obtain its evaluation obtained from the formula:

my_jack = Palpha[[2,1,0],[3,1]]
eval_jack = my_jack.evaluation(6)
print(factor(SR(eval_jack)))

# Now let us do the same in the variable representation. The code will run as follows:
#
# 1. We generate a superspace with $N=6$ and expand the Jack superpolynomial on it.
# 2. From the variables expression, we only keep the coefficient of $\theta_0 \theta_1 \theta_2$, that's what `superspace.theta_project(expr)` does
# 3. Generate the Vandermonde determinant in 3 variables
# 4. Factor out the Vandermonde determinant from the expression obtained in 2
# 5. Substitute x_i = 1
# 6. Print the result

# Superspace with N = 6
ss = superspace(6)
# Expressions of the Jack in terms of x_i, theta_i
expr = my_jack.expand(ss)
# We pick up the coefficient of theta_0 ... theta_{m-1}
p_expr = ss.theta_project(expr)
# Vandermonde determinant in 3 variables
vander = ss.vandermonde(3)
# We make the division. Here we use a Singular function for this 
# puprose since p_expr/vander will return an error. We have to do this
# because 1/x_i is not part of the ring. The result of the division is
# returned in a table hence the [1][1,1], see Singular documentation.
n_p_expr = p_expr.division(vander)[1][1,1]
# We get the list of variables
X = ss.x_vars()
# We prepare the argument for substitution, the following yields a list
# of the form [x_0, 1 x_1, 1, ...]
subs_arg = flatten([[x, 1] for x in X])
# We call the Singular subst function on our polynomial
evaluation = n_p_expr.subst(*subs_arg)
# We convert the result to a string to convert it to the symbolic ring
# without hassle, once on the symbolic ring, we factor it for formating
evaluation = factor(SR(str(evaluation)))
# We print the result
print(evaluation)

# which is what we obtained with the formula. 

# ### Integral form of the Jack superpolynomials
# Here is a simple example of expansion of the integral form of the Jack polynomial:

super_init()
Palpha = Sym.Jack()
sparts = Superpartitions()

# +
spart = sparts([[2,0],[3,1]])

jack_m = m(Palpha(spart))
print(jack_m)
# -

# Given this super-Jack, we normalize it to obtain the integral form:

# +
nn = spart.bosonic_degree()
mm = spart.fermionic_degree()
l_nm = nn - mm*(mm-1)/2

sparts_in_sector = list(Superpartitions(*spart.sector()))
lowest_spart = Superpartitions.sort_by_dominance(sparts_in_sector)[-1]
c_min = jack_m.monomial_coefficients()[lowest_spart]

Jlambda_m = factorial(l_nm)*(1/c_min)*jack_m.collect()

print(Jlambda_m)
# -

# The coefficients are seen to be $\in \mathbb{Z}(\alpha)$. 

# #### Conjecture 
# The coefficient of $m_{\Lambda_{\text{min}}}$ is given by 
# $$
#     (n-m(m-1)/2)!v_\Lambda(\alpha)^{-1}
# $$
# with $v_\Lambda(\alpha) = \prod_{s\in \mathcal{B}} h^{\text{lo}}_\Lambda(s;\alpha)$
#
# This can be easily checked and is left as an exercise. We provide a list of useful functions:
#
# `Superpartitions.sort_by_dominance( spart_list )`
#
# `a_spart.lower_hook_length(i, j, alpha)`
#
# `a_spart.bosonic_cells()`

# +
# Exercise
# Verify the previous formula for the sector (6|3)

# Initialization #
super_init()
Palpha = Sym.Jack()
alpha = QQqta.gens_dict()['alpha']

# Obtain the superpartitions of the sector (6|3) as a list

# Find the lowest superpartition

# Generate each Jack for earch superpartition

# Obtain the coefficient of \Lambda_min for every Jack

# Compute the formula for each Jack

# Compare the two results
# -

# #### Conjecture
# The coefficients in the expansion of the Jack superpolynomial in the integral from are $\in \mathbb{Z}(\alpha)$

# We now illustrate the conjecture. We are going to do the following:
# 1. Make a function which returns the integral form on the monomial basis
# 2. Generate the Jack superpolynomials of a given sector
# 3. Obtain the integral form with the function made in 1.
# 4. For each expression, we are going to extract the coefficients in a list
# 5. We then convert those coefficients to the symbolic ring in order to have a neat factorization.
# 6. We print the coefficient of every monomial belonging to every Jack in the sector

super_init()
Palpha = Sym.Jack()
alpha = QQqta.gens_dict()['alpha']


def integral_form(expr_m):
    spart_coeff = expr_m.monomial_coefficients()
    sparts = spart_coeff.keys()
    smallest_spart = Superpartitions.sort_by_dominance(sparts)[-1]
    sector = smallest_spart.sector()
    l_nm = sector[0] - sector[1]*(sector[1]-1)/2
    c_min = spart_coeff[smallest_spart]
    int_form = factorial(l_nm)*(1/c_min)*expr_m
    int_form = int_form.collect()
    return int_form


sparts = Superpartitions(5,2)
Jacks_m = [m(Palpha(spart)) for spart in sparts]
intForm_m = [integral_form(jack) for jack in Jacks_m]
Jcoeffs = [J.coefficients() for J in intForm_m]
coeffs_SR = [factor(SR(coeff)) for J in Jcoeffs for coeff in J]
for coeffs in coeffs_SR:
    display(Math(latex(coeffs)))

# We thus see that all coefficients are $\in \mathbb{Z}(\alpha)$. 

# ## Duality
# We illustrate the duality property.
#
# We first pick a superpartition and obtain the Jack superpolynomial associated with it.

sparts = Superpartitions()
spart1 = sparts([[4,1],[]])
spart1.terminal_diagram()
a_Jack = Palpha(spart1)
print(a_Jack)

# We compute the norm of this polynomial and apply $\hat{\omega}_\alpha$ on the Jack. We then divide the result by the norm.

#We get the norm from the scalar product
the_norm = a_Jack.scalar_alpha(a_Jack)
#We apply omega alpha
omega_a_Jack = a_Jack.omega_alpha()
#We then normalize
normalized_om_Jack = omega_a_Jack/the_norm
print(normalized_om_Jack.collect())

# Now, according to the duality property, this polynomial should be $P^{1/\alpha}_{(1,0;2,1,1)}$. To check this, we develop this polynomial on the monomial basis so that no $\alpha$ is hidden inside a definition and substitute $\alpha$ for $1/\alpha$. Converting back to the Jack basis, we should find $P^\alpha_{(1,0;2,1,1)}$
#
# Note that here we use the in-house method `expr.subs_coeff(dictionnary)` to make the substitution $\alpha \longrightarrow 1/\alpha$

# We first get the monomial expansion of the sJack 
# to unravel the hidden alpha dependence 
m_expr = m(normalized_om_Jack)
# alpha -> 1/alpha
invert_alpha_m = m_expr.subs_coeff({alpha:1/alpha})
# We convert the resulting expression back to the sJack basis. 
resulting_jack = Palpha(invert_alpha_m)
print(resulting_jack)
new_spart = resulting_jack.monomial_coefficients().keys()[0].terminal_diagram()

# Which is indeed the expected result. 

# ## Particular expressions

# ### $P^{\alpha}$ vs $g^\alpha$
# Let us illustrate the following properties:

# $$
# P^{(\alpha)}_{(;n)} = \frac{\alpha^n n!}{(1 + \alpha(n-1))_n} g_n
# \qquad \text{and}\qquad 
# P^{(\alpha)}_{(n;)} = \frac{\alpha^{n+1} n!}{(1 + \alpha n)_n} \tilde{g}_n
# $$
# For $\Lambda = (;\, 5)$ we have

# We initialize the space
super_init()
Palpha = Sym.Jack()
galpha = Sym.Galpha()

jack_g = galpha(Palpha[[],[5]]).collect()
print(jack_g)
display(Math(latex(factor(SR(jack_g.coefficients()[0])))))

# and for $\Lambda = (5;\,)$:

jack_g2 = galpha(Palpha[[5],[]]).collect()
print(jack_g2)
display(Math(latex(factor(SR(jack_g2.coefficients()[0])))))

# $\lim_{\alpha \rightarrow \infty} P^{(\alpha)}_{\Lambda} = m_\Lambda$

aJack = Palpha[[2,0],[3,1]]
print(aJack)
aJack_m = m(aJack)
# Here the expr.subs_coeff(dictionnary) is not appropriate
# for limiting process. So we use the more involved 
# expr.map_coefficients(func) provided by SAGE. 
# Note also that we transfer the coefficients to the symbolic ring
# in order to evaluate the limit. 
lim_aJack = aJack_m.map_coefficients(
    lambda x: QQa(SR(x).limit(alpha = infinity)))
print("alpha --> \\infty")
print(lim_aJack)

# $
# P^{(0)}_\Lambda = (-1)^{\binom{m}{2}} e_{\Lambda^\prime}
# $

alpha_zero_Jack = aJack_m.subs_coeff({alpha:0})
print(aJack)
print("alpha --> 0")
print(e(alpha_zero_Jack))

# $P^{(1)}_{(;n)} = h_n$

nJack_m = m(Palpha[[],[5]])
alpha_one = nJack_m.subs_coeff({alpha:1})
print(h(alpha_one))

# $P^{(1)}_{(n;)} = \frac{1}{n+1} \tilde{h}_n$

nJack_m = m(Palpha[[5],[]])
alpha_one = nJack_m.subs_coeff({alpha:1})
print(h(alpha_one))


# ## Negative fractional parameter
# ### Admissible Jack superpolynomials
# We first build a function that tells us if a superpartition is $k,r$ admissible:

def is_krN_admissible(spart, k, r, N):
    assert(k>=1)
    assert(r>=2)
    assert(N>=k)
    # We check the coprime condition
    assert(gcd(k+1,r-1)==1)
    L_star = list(spart.star())
    add_zeros = [0 for k in range(N-len(L_star))]
    L_star = L_star + add_zeros
    L_cstar = list(spart.circle_star())
    add_zeros = [0 for k in range(N-len(L_cstar))]
    L_cstar = L_cstar + add_zeros
    for i in range(N-k):
        diffL = L_cstar[i] - L_star[i+k]
        if  diffL >= r:
            pass
        else:
            return False
    return True


# Now we set $k,\,r$ and $N$ a find admissible superpartitions, in the present case we will find one.

admissible = []
k, r, N= (1, 2, 3)
sectors = Sym._Jack_m_cache.keys()
sparts = Superpartitions(6,2)
for spart in sparts:
    if is_krN_admissible(spart, k, r, N):
        admissible.append(spart)
print(admissible)

# Given this superpartition, we will generate the sJack associated to it and set $\alpha = -\frac{k+1}{r-1}$:

ad_spart = admissible[0]
super_init()
alpha = QQqta.gens_dict()['alpha']
Palpha = Sym.Jack()
jack_m = m(Palpha(ad_spart))
alpha_spec = -(k+1)/(r-1)
jack_akr = jack_m.subs_coeff({alpha:alpha_spec})
print("\n")
print(jack_akr)

# So we see that this function does not have poles. We now expand on the variables to see the vanishing and clustering property:

ss = superspace(N)
in_var = jack_akr.expand(ss)
print(in_var)

# The Jack polynomial in $\alpha_{k,r}$ with admissible partition vanishes whenever $k+1$ of its variable are equal:

for x,y in [[x_0,x_1], [x_0, x_2], [x_1,x_2]]:
    k1eq = in_var.subst(x, y)
    print(k1eq)

# We finally illustrate the clustering property

projected = ss.diff(ss.diff(in_var, theta_0), theta_1)
sr_expr = factor(SR(str(projected)))
print(sr_expr)

# Another clustering exemple:
#
# $N=4, k=2, r=2$. 
#
# With those parameters, we select a set of superpartitions with $m=1$ that are admissible. For each of those superpartitions, we obtain the Jack polynomial with the proper $\alpha$ specialization. Then we expand the superpolynomials in variables and set $x_3 = x_2$. We then see that each expression goes to 0 as $(x_1 - x_2)^2$ for $x_1 \rightarrow x_2$. 

admissible = []
k, r, N= (2, 2, 4)
sectors = Sym._Jack_m_cache.keys()
for sector in sectors:
    sparts = Superpartitions(*sector)
    for spart in sparts:
        if len(spart) <= N:
            if is_krN_admissible(spart, k, r, N):
                if spart.fermionic_degree() == 1:
                    admissible.append(spart)

jacks_m = [m(Palpha(spart)) for spart in admissible]
alpha_spec = -(k+1)/(r-1)
jack_akr = [jack.subs_coeff({alpha:alpha_spec}) for jack in jacks_m]
ss = superspace(N)
jacks_x = [jack.expand(ss) for jack in jack_akr]
projected_jacks = [ss.diff(in_var, theta_0) for in_var in jacks_x]
sr_expr = [factor(SR(str(projected.subst(x_3,x_2)))) for projected in projected_jacks]
for expr in sr_expr:
    display(Math(latex(expr)))

# The previous strategy is very computationally expensive for large $N$. A better way approach goes through the non-symmetric Jack polynomial. The final step where we project the expression on $\theta_0 \cdots \theta_{m-1}$ amounts to reduce the whole superpolynomial to the non-symmetric Jack with prescribed symmetry (antisymmetrized on the first $m$ variables). If one desires to explore the conjecture further, this would be the way to go. 

# # Macdonald superpolynomials

QQqt = FractionField(QQ['q', 't'])
q, t = QQqt.gens()
symqt = SymSuperfunctionsAlgebra(QQqt)
symqt.inject_shorthands()
Pqt = symqt.Macdonald()

# As usual the basis works either by invoking `Pqt(spart)` or `Pqt[[...],[...]]`

sparts = Superpartitions()
my_spart = sparts([[3,2,1],[3,2]])
print(Pqt(my_spart))
print(Pqt[[4,1],[3,2]])

# ## Particular cases
# ### Jack limit
# $q=t^\alpha, t\rightarrow 1$
#

# We would like to test the limit from the Macdonald to the Jack case. So we must declare the ring of coefficients with all three parameters:

QQqta = FractionField(QQ['q','t','alpha'])
Symqta = SymSuperfunctionsAlgebra(QQqta)
Symqta.inject_shorthands()
Pqt = Symqta.Macdonald()
Palpha = Symqta.Jack()

# So we are in a space of symmetric superfunctions where both bases are well defined as we illustrate here:

MyJack = Palpha[[1,0],[2]]
MyJack_m = m(MyJack)
print(MyJack_m)
MyMacdo = Pqt[[1,0],[2]]
MyMacdo_m = m(MyMacdo)
print(MyMacdo_m)


# Let us now test the Jack limiting case. Now we confronted with the problem that for Sage $t^{\alpha}$ is not part of the coefficient ring. Instead of simply providing a function for the evaluation of this limit, we will here define a Sage function from scratch. This will serve as an introduction to how one would implement one's own features and present subtle details that are not apparent otherwise. 
#
# So we have to design a function that takes a coefficient and does the proper substitution and take the proper limit. More precisely it has to do the following: 
#
# 1. Save the original coefficient ring.
# 2. Convert the coefficient and q,t,alpha to the symbolic ring.
# 3. Make the substitution `q = t^alpha`.
# 4. Take the limit `t=1`.
# 5. Convert the result back to the original coefficient ring.

def mac_jack_coeff(coeff):
    # First, if the coefficient is not a polynomial in either
    # q or t we don't have to do anything
    if coeff in QQ:
        return coeff
    
    # Now we must extract the coefficient ring
    BR = coeff.parent()
    # This is a dictionnary with q, t and alpha
    parameters = BR.gens_dict()
    q = parameters['q']
    t = parameters['t']
    alpha = parameters['alpha']
    
    # Now we need to convert our coefficient to the symbolic ring
    coeff_SR = SR(coeff)
    # We also need a version of the parameters on the symbolic ring
    q_SR = SR(q)
    t_SR = SR(t)
    alpha_SR = SR(alpha)
    
    # We can now do the substitution since it is allowed on the 
    # symbolic ring
    coeff_SR_qtalpha = coeff_SR.subs({q_SR:t_SR**alpha_SR})

    # To take the limit, we have to prepare the underlying 
    # engine (GAP) that will compute it. First, we tell
    # SAGE that alpha is in ZZ
    assume(alpha_SR, 'integer')
    # And here, one must understand that the limit is computed
    # by GAP and that the limit argument is sent by sage as a
    # string. So instead of writing t_SR = 1, we must write 
    # t = 1 since t_SR is represented as the string 't' 
    # in the equations
    coeff_jack = coeff_SR_qtalpha.limit(t = 1)
    # We simplify the resulting coefficient
    coeff_jack = factor(coeff_jack)
    
    # Now we must convert the result of our computation back to
    # the coefficient ring. To do so, we have to substitute each
    # variable since Sage does not understand automaticaly that,
    # for instance, q_SR from SR is related to q taken from BR. 
    # So we substitute with the right variable and apply BR()
    # to convert it to the right ring. 
    BR_coeff = BR(coeff_jack.subs({alpha_SR : alpha}))
    
    return BR_coeff  
    


# We now have a function that maps the coefficients. So we just have to go on and apply it to every coefficient of a polynomial to verify the limiting behaviour. To procede, we are going to:
#
# 1. Generate a Macdonald superpolynomial and expand it on the monomial basis.
# 2. Use the `expr.map_coefficients(FUNCTION)` method to apply our brand new function to the coefficient.
# 3. Convert the resulting expression to Jack superpolynomials.

# +
# Step 1
mac = Pqt[[3,0],[2]]
mac_m = m(mac)
print("Starting from " + str(mac) + "\n")

# Step 2
print("We apply the limit resulting in")
# Note here that all we need to provide to map_coefficients
# is the function name
mac_m_apply_limit = mac_m.map_coefficients(mac_jack_coeff)
print(mac_m_apply_limit)
print("\n")
#
resulting_jack = Palpha(mac_m_apply_limit)
print("Converting to Jack superpolynomials we have")
print(resulting_jack)
# -

# ### Monomial limit
# $t=1$
#
# $P_{\Lambda}(q, 1) = m_\Lambda$

a_mac = Pqt[[2,1],[2]]
print("Starting from " + str(a_mac) + "\n")
a_mac_m = m(a_mac)
mono = a_mac_m.subs_coeff({t:1})
print("We substitute t=1 and obtain :")
print(mono)

# ### Elementary limit
# $q=1$
#
# $P_\Lambda(1,t) = (-1)^{\binom{m}{2}}e_{\Lambda^\prime}$

# +
a_mac = Pqt[[1,0],[1,1]]
a_mac_e = e(a_mac)
print("Starting from " + str(a_mac) + "\n")
# On the following line we take the expression a_mac
# expr.monomial_coefficients() will extract a dict of the form
# {spart:coeff, spart:coeff, ...}, using keys() we obtain
# a list of the superpartition, then [0] exctract the first
# (and only) superpartition
spart = a_mac.monomial_coefficients().keys()[0]
spart.terminal_diagram()

# Now for the substitution
elem = a_mac_e.subs_coeff({q:1})
print("We expand on the e basis and evaluate q=1 obtaining:")
print(elem)
elem.monomial_coefficients().keys()[0].terminal_diagram()
# -

# ### Schur
# The super-Schur function is presented in its own section.

# ## Cauchy kernel
#
# ### $g_\Lambda(q,t)$
# We now introduce the 2 parameter deformation of the homogeneous basis. 

# +
QQqt = FractionField(QQ['q','t'])
Sym = SymSuperfunctionsAlgebra(QQqt)
Sym.inject_shorthands()
Sparts = Superpartitions()

print('\nAs with any other basis')
gqt = Sym.Gqt()
my_expr = gqt[[3,0],[1]]
print(my_expr)
print("\n we can convert it to other bases")
print(m(my_expr).collect())
# -

# Now a defining property of this basis is that it is orthonormal to the monomial basis, so let us illustrate that in for superpartitions in the sector $(6|3)$:

spart_set = list(Superpartitions(6,3))
ortho_pairs = Combinations(spart_set, 2)
ortho_scal = [gqt(spart1).scalar_qt(m(spart2))
              for spart1, spart2 in ortho_pairs]
unit_scal = [gqt(spart).scalar_qt(m(spart))
             for spart in spart_set]
print("Lambda =/= Omega")
print(ortho_scal)
print("\nLambda == Omega")
print(unit_scal)

# ## Macdonald superpolynomials from an eigenvalue problem

# ## Further properties 
# ### One part superpartition
# $P_{(n;)} = \left( \sum_{i=0}^n q^{i-n} \chi_i \right)^{-1} \tilde{g}_n(q,t)$

# +
QQqt = FractionField(QQ['q','t'])
Sym = SymSuperfunctionsAlgebra(QQqt)
Pqt = Sym.Macdonald()
gqt = Sym.Gqt()

one_part_mac = Pqt[[5],[]]
print(gqt(one_part_mac).collect())
# -

# $P_{(;\, n)} = \chi_n^{-1} g_n(q,t)$ 

one_part_mac2 = Pqt[[],[5]]
print(gqt(one_part_mac2))

# ### Duality
# We now illustrate the duality property. 
#
# $$\hat{\omega}_{q,t} P_\Lambda (q, t) = 
# \frac{(q/t)^{|\Lambda|}}{<P_\Lambda | P_\Lambda>_{q,t}}P_{\Lambda^\prime} (t^{-1}, q^{-1})$$
#
# To illustrate this, we will follow these steps: 
#
# 1. We generate a superMacdonald for a superpartition
# 2. We apply the omega_qt automorphism
# 3. We normalize the result (we devide by the factor on the rhs)
# 4. We expand the expression on the monomial basis to unravel the q, t dependence.
# 5. We substitute q=1/t and t=1/q
# 6. Convert the result back to Macdonald, this should result in the Macdonald of the conjugate partition.

# +
# We generate a macdonald polynomial
spart1 = Superpartitions()([[2,0],[1,1]])
a_macdo = Pqt(spart1)
print(a_macdo)
spart1.terminal_diagram()

# We apply omega
w_a_macdo = a_macdo.omega_qt()

# We compute the norm 
norm = a_macdo.scalar_qt(a_macdo)

# We normalize the result (DOES NOT FIT WITH Compendium document)
new_mac = 1/norm*w_a_macdo

# Expand on monomial basis
new_mac_m = m(new_mac)
# Substitution
new_mac_m = new_mac_m.subs_coeff({q:1/t, t:1/q})

# Convert back to Macdonald
new_mac = Pqt(new_mac_m)
print(new_mac.collect())
new_mac.monomial_coefficients().keys()[0].terminal_diagram()
# -

# ### The second duality
# The second duality is expressed through the following automorphism:
# $$
#     \hat{\rho}_{q,t}(p_r) = (-1)^{r-1} \frac{1-q^r}{1-t^r} p_r, \quad
#     \hat{\rho}_{q,t}(\tilde{p}_r) =
#     \sum_{\Lambda \vdash (r|1)} \frac{\omega_\Lambda}{z_\Lambda}
#     \prod_{i}(1-q^{\Lambda^s_i})p_\Lambda
# $$
# with $\omega_\Lambda = (-1)^{|\Lambda| - \ell (\Lambda^{s})}$. This automorphism is provided as a method on expressions as `expr.rho_qt()` :

super_init()
expr_p = p[[2,0],[3,1]].rho_qt()
print(expr_p)
print('\n')
expr_m = m[[2],[1]].rho_qt()
print(expr_m)

# We now illustrate the duality through the following steps:
# 1. We choose a superpartition and get the associated Macdonald superpolynomial
# 2. We could apply $\hat{\rho}_{q,t}$ right away, but since it goes through the powersum basis and we'll need to go back to this basis we convert our Macdonald to powersum and apply $\hat{\rho}_{q,t}$
# 3. We make the substitution $q\leftrightarrow  t$
# 4. Convert the expression in the powersums back to the Macdonald basis
# 5. We print the result, the diagram and the coefficient.

Pqt = Sym.Macdonald()
q, t, alpha = QQqta.gens()
# We generate our partition and print the diagram
my_spart = Sparts([[2,0],[2,1,1]])
my_spart.terminal_diagram()
# We now take the Macdonald associated to the spart and
# apply rho, we will convert it to the powersum basis to
# speed up the following steps
mac = Pqt(my_spart)
rho_mac_p = p(mac).rho_qt()
# Now rho_mac should be Q of conjugate superpartition
# with t <--> q, so we develop on the m basis, exchange q and t
# and finally we put it back on the Jack basis
rho_mac_p_tq = rho_mac_p.subs_coeff({q:t, t:q})
rho_mac_conj = Pqt(rho_mac_p_tq)
print("")
print(rho_mac_conj)
new_spart, coeff = rho_mac_conj.monomial_coefficients().items()[0]
new_spart.terminal_diagram()
print("The coefficient is")
display(Math(latex(factor(SR(coeff)))))

# ### Normalization
# $$ < P_\Lambda(q,t) | P_\Lambda(q,t)>_{q,t} =
# (-1)^{\binom{m}{2}} q^{| \Lambda^a |}
# \prod_{s \in \mathcal{B} \Lambda}
# \frac{1 - q^{a(s) + 1}t^{\tilde{l}(s)}}
# {1 - q^{\tilde{a}(s)}t^{l(s) + 1}}
# $$
#
# We verify some examples of the conjectured normalization formula here. The value of the given by the formula can be obtained with
# `Pqt.calc_norm(spart)`, where `Pqt` must be an instance of `Sym.Macdonald()`

Sparts = Superpartitions(5,3)
Pqt = Sym.Macdonald()
for a_spart in Sparts:
    a_Mac = Pqt(a_spart)
    scalar_prod = a_Mac.scalar_qt(a_Mac)
    calc_norm = Pqt.calc_norm(a_spart)
    display(Math(str(a_spart)+':\,' +latex(factor(SR(calc_norm)))))
    display(Math(latex(factor(SR(scalar_prod)))))
    print("\n")

# ### Specialization
# We now illustrate the specialization conjecture. We present an example by going through the following steps:
# 1. Create a Macdonald polynomial
# 2. Expand it on an appropritate superspace
# 3. Project the expression to obtain only the coefficient of $\theta_0\cdots \theta_{m-1}$ and factor out the Vandermonde determinant
# 4. Obtain the $\delta^k$ staircase and use it to create the list of substitution
# 5. Make the substitutions and print the result

super_init()
Pqt = Sym.Macdonald()
q,t,_ = QQqta.gens()
sparts = Superpartitions()
# We use the Macdonald with the following superpartition
spart = sparts([[3,1,0],[2,1]])
a_mac = Pqt(spart)
# We use a number of variable equivalent to the longest 
# superpartition in the sector
N = max([len(x) for x in Superpartitions(*spart.sector())])
ss = superspace(N)
# Expand the macdonald polynomial
mac_x = a_mac.expand(ss)
# Projection
projected = ss.theta_project(mac_x)
# Now we must divide by the Vandermonde determinant
projected = projected.division(ss.vandermonde(spart.fermionic_degree()))[1][1,1]
# Staircase
delta_m = list(sparts.stair(spart.fermionic_degree()-1))
delta_m = delta_m + [0 for k in range(N - len(delta_m))]
X = ss.x_vars()
# Substitution list
sub_list = flatten([[X[k], q**(-delta_m[k])*t**(k)] for k in range(N)])
# Substitution
mac_spec = projected.subst(*sub_list)
# Use the symbolic ring to print the result nicely
spec_SR = factor(SR(str(mac_spec)))
display(Math(latex(spec_SR)))
conjectural_spec = a_mac.specialize(N)
print(factor(SR(str(mac_spec)))/factor(SR(str(conjectural_spec))))


# ### Symmetry
# To illustrate the symmetry conjecture we will define a function that applies the specialization. If you read the specialization section you should understand what is going on. 

def spec_macdo(Pqt, mac_spart, spec_spart, ss):
    BR = Pqt.base_ring()
    params = BR.gens_dict()
    q = params['q']
    t = params['t']
    N = ss._N
    # Obtain the Macdonald
    mac_x = Pqt(mac_spart).expand(ss)
    # Get the coefficient of theta_0...theta_{m-1}
    projected = ss.theta_project(mac_x)
    ferm_degree = mac_spart.fermionic_degree()
    projected = projected.division(
        ss.vandermonde(ferm_degree))[1][1,1]
    
    X = ss.x_vars()
    # Substitution list 
    spec_composition = list(spec_spart[0]) + list(spec_spart[1]) 
    spec_composition = (spec_composition +
                        [0
                         for k
                         in range(N - len(spec_composition))])
    index_part = sorted(
        enumerate(spec_composition, 0),
        key=lambda x: x[1], reverse=True)
    sub_list = flatten([
        [X[index_part[k][0]], q**(-index_part[k][1])*t**k]
        for k in range(N)]) 
    # Substitute
    spec = projected.subst(*sub_list)
    return spec


# Now we test our function, if we set the `spec_spart` to $(2,1,0;\,)$ and the `mac_spart` to $(3,1,0;2,1)$ we should recover the result from the last section.

super_init(); Pqt = Sym.Macdonald()
spec_spart = Sparts([[2,1,0],[]])
mac_spart = Sparts([[3,1,0],[2,1]])
ss = superspace(7)
res = spec_macdo(Pqt, mac_spart, spec_spart, ss)
display(Math(latex(factor(SR(str(res))))))

# So it's all fine. With that in hand, we are now in position of testing the conjecture. We are going to take every pair of superpartitions of the sector $(6|3)$, compute their respective $\tilde{P}_\Lambda(u_\Omega)$ and compare the results.

# +
sparts = list(Superpartitions(6,3))
delta = Sparts([[2,1,0],[]])
ss = superspace(max([len(x) for x in sparts]))

spart_pairs = Combinations(sparts, 2)
for spart_pair in spart_pairs:
    Lambda, Omega = spart_pair
    
    PLa_stair = spec_macdo(Pqt, Lambda, delta, ss)
    PLa_Om = spec_macdo(Pqt, Lambda, Omega, ss)
    PtildeLa_Om = PLa_Om.division(PLa_stair)[1][1,1]
    
    POm_stair = spec_macdo(Pqt, Omega, delta, ss)
    POm_La = spec_macdo(Pqt, Omega, Lambda, ss)
    PtildeOm_La = POm_La.division(POm_stair)[1][1,1]
    
    print(PtildeOm_La == PtildeLa_Om)
# -

# ### Integral form and the positivity
# The integral form is given by 
# $$
#     J_\Lambda(q,t) = h^{\text{lo}}_\Lambda(q,t) P_\Lambda(q,t)
# $$
# The $h^{\text{lo}}_\Lambda$ quantity can be obtained with `pqt_expr.hlo_Lambda(q,t,spart)`

super_init()
Pqt = Sym.Macdonald()
q, t, alpha = QQqta.gens()
spart = Sparts([[2,0],[3,1]])
expr = Pqt(spart)
hlo = expr.hlo_Lambda(q,t,spart)
print(factor(SR(hlo)))

# It can be seen that this normalization gives coefficients $\in \mathbb{Z}(q,t)$ when expanded on the monomial basis:

Jmono = m(hlo*expr)
coeffs = Jmono.monomial_coefficients().values()
for coef in coeffs:
    display(Math(latex(factor(SR(coef)))))

# So we introduce the automorphism $\phi_t$

expr = p[[2,0],[3]] + p[[1,0],[3,1]]
print(expr.phi_t())
print('')
expr_m = m[[1,0],[3]]
print(expr_m.phi_t())

# With this in hand we can obtain the modified Macdonald superpolynomials

spart = Sparts([[2,0],[2,1]])
mac = Pqt(spart)
q, t, alpha = QQqta.gens()
mod_mac = mac.hlo_Lambda(q, t, spart)*mac.phi_t()
print(mod_mac)

# Converting this expression to the Schur basis allows us to illustrate the conjecture, the super-Schur function is presented in a later section.
#
# Here we convert the expression to Schur, then we for each superpartition, we print the coefficient.

s = Sym.Schur()
mod_schur = s(mod_mac)
spart_Kostka = mod_schur.monomial_coefficients().items()
for spart, coeff in spart_Kostka:
    display(Math(str(spart) + ":" + latex(factor(SR(coeff)))))

# And so we see that the coefficients are all $\in \mathbb{N}(q,t)$.

# ### Antoher scalar product
# This scalar product is more interesting from the analytical point of view, thus it is not implemented.

# ## Constructing the Macdonald superpolynomials from the non-symmetric Macdonald polynomials

# ### Cherednik algebra
# The Demazure-Lusztig operators are part of Sage libraries but not as operators on polynomials. So we have our own implementation. For instance:

# +
ss = superspace(3)
Ti = ss.T_i

my_expr = x_0 + x_1 + x_2
for k in range(2):
    print(Ti(my_expr, k))
# -

# It is also implemented for a serie of elementary permutations $s_i$:

# +
S3 = ss.SN()
as_elementary = [omega.reduced_word() for omega in S3]
Tw = ss.Tw

expr = (x_0-x_1)*x_2^2
apply_Tw = [Tw(expr, w) for w in as_elementary]
print(apply_Tw)

# -

# ### The non-symmetric Macdonald polynomials
# Now using the $T_i$ operators, we are going to obtain Macdonald superpolynomials from the non-symmetric Macdonald polynomials.
#
# We have are going to use the following simple form 
# $$
# P_\Lambda (q,t) \propto \sum_{\sigma \in S_N} \sigma 
# \theta_0 \cdots \theta_{m-1}
# U^+_{m,...,N-1}
# E_{\Lambda^R}(q,t)
# $$

# The procedure is explained as comments in the following code:

# We get an instance of Symmetric superpolynomials
super_init()
Pqt = Sym.Macdonald()

# +
# First, we need the non-symmetric Macdonald polynomial
# We import it
from sage.combinat.sf.ns_macdonald import E
# We define the q,t parameter and set the superspace with N=5
q, t = QQ['q','t'].fraction_field().gens()
N = 5
ss = superspace(N)


# We define the composition Lambda^R, let's say we want
# Lambda = (2,0;1,1), this yields the following composition
# and sector (mm stands for m)
composition = [0,2,0,1,1]
mm = 2

# Now we need the nsMacdonald, the Sage implementation does not
# follow the same convention we use. So we have to reverse the
# composition and also reverse the variables.
composition.reverse()
perm_vars = range(1,N+1)
perm_vars.reverse()
perm = Permutation(perm_vars)
# Now we get the non-symmetric macdonald polynomial
nsMacdo = E(composition, q, t, pi=perm)

# This polynomial is not in the right ring, so we convert it
# to a singular object
string_list = str(nsMacdo).split('x')
singular_str = 'x_'.join(string_list)
nsMacdo_singular = singular(singular_str).normalize()
# We now need to t-symmetrize on the variables mm...N-1
SNslash = SymmetricGroup(range(mm, N))
# we express permutations as chains of s_i
elem_perms = [omega.reduced_word() for omega in SNslash]
# We can now apply the Tw on the nsMacdonald and sum the result
# to get a Macdonald with prescribed symmetry
terms = [ss.Tw(nsMacdo_singular, w) for w in elem_perms]
# We also normalize the leading terms to limit proliferation of 
# q,t factors 
prescMac = sum(terms).normalize()

# We now need to multiply this by the theta variables
theta = ss.theta_vars()
super_prescMac = prod(theta[:mm])*prescMac

# We symmetrize the result and normalize
macdo = ss.symmetrize(super_prescMac).normalize()

# We convert this back in the Macdonald basis
macdo2 = Pqt.from_polynomial(macdo, ss)
print(macdo2)
# -

# Which is the expected result. You can try with your own compostion. The composition must be exactly of length $N$ to get sensible results. 

# # The double Macdonald polynomials

# # The super-Schur and their Pieri rules 
# The two types of Schur superpolynomials can be obtained from `sym.Schur()` and `sym.SchurBar()`. Note that your ring of symmetric superfunctions must have coefficients in $q,t$ since the expansion on other basis is obtained from the limiting process. 

super_init();
s = Sym.Schur()
sbar = Sym.SchurBar()

# We also have their dual superfunctions, respectively:

sStar = Sym.SchurStar()
sbarStar = Sym.SchurBarStar()

# We see that these are indeed dual:

sparts = Superpartitions(5,2)
scals = [s(spart).scalar_product(sStar(spart)) for spart in sparts]
display(Math('<s_\Lambda|s^*_\Lambda>'))
print(scals)

sparts = Superpartitions(5,2)
scals = [sbar(spart).scalar_product(sbarStar(spart)) for spart in sparts]
display(Math('<\overline{s}_\Lambda|\overline{s}^*_\Lambda>'))
print(scals)

# We now illustrate some standard operations that work out of the box. Note that the first two expression are indeed triangular on the monomial basis.

spart = sparts([[2,0],[2,2]])
print(m(s(spart)))
print('')
print(m(sbar(spart)))
print('')
print(m(sStar(spart)))
print('')
print(m(sbarStar(spart)))
print('')
print(s(sbar(spart)))
print('')
print(sStar(sbar(spart)))

# ## Pieri rules
# We first illustrate the one-row and one-column super-Schur identities:

display(Math("s_\Lambda: "))
s_row_ferm = p(s[[3],[]])
s_row_bos = h(s[[],[3]])
s_col_ferm = e(s[[0],[1,1,1]])
s_col_bos = e(s[[],[1,1,1]])
display(Math(latex(s_row_ferm)))
display(Math(latex(s_row_bos)))
display(Math(latex(s_col_ferm)))
display(Math(latex(s_col_bos)))

display(Math("s^*_\Lambda: "))
s_row_ferm = h(sStar[[3],[]])
s_row_bos = h(sStar[[],[3]])
s_col_ferm = e(sStar[[0],[1,1,1]])
s_col_bos = e(sStar[[],[1,1,1]])
display(Math(latex(s_row_ferm)))
display(Math(latex(s_row_bos)))
display(Math(latex(s_col_ferm)))
display(Math(latex(s_col_bos)))

# We illustrate the Pieri rules with some examples. 
# ### A-

s[[2,0],[1]]*s[[3],[]]

s[[2,0],[1]]*s[[],[3]]

# ### B-

s[[2,0],[1]]*s[[0],[1,1,1]]

s[[2,0],[1]]*s[[],[1,1,1]]

# ### C-

sStar[[1],[1]]*sStar[[2],[]]

sStar[[1],[1]]*sStar[[],[2]]

# ### D-

sStar[[1],[1]]*sStar[[0],[1,1,1]]

sStar[[1,0],[1]]*sStar[[],[1,1]]

StarcoderdataJupyterScriptsDedupFiltered