################################################################################
# Kostka Polynomials
#
# Copyright (C) 2020 Ulrich Thiel, ulthiel.com/math
################################################################################

export kostka_polynomial, charge


"""
	kostka_polynomial(Î»::Partition{T}, Î¼::Partition{T})
	kostka_polynomial(Î»::Array{Int,1}, Î¼::Array{Int,1})

The (one-variable) **Kostka polymial** ``K_{Î»Î¼}(t)`` associated to partitions Î» and Î¼ can be defined as

```math
K_{Î»Î¼}(t) = âˆ‘_{T âˆˆ SSYT(Î»,Î¼)} t^{charge(T)} âˆˆ â„•[t]	\\;,
```

where SSYT(Î»,Î¼) is the set of all semistandard Young tableaux of shape Î» and weight Î¼, and charge(T) is the charge of the tableau T. The Kostka polynomials relate the Hallâ€“Littlewood polynomials ``P_Î¼``
to the Schur polynomials ``s_Î»`` via
```math
s_Î»(x_1,â€¦,x_n) = \\sum_Î¼ K_{Î»Î¼}(t)P_Î¼(x_1,â€¦,x_n;t)
```
This function returns the Kostka polynomial ``K_{Î»Î¼}(t)`` as an fmpz_poly over ZZ in t.

# Example
```julia-repl
julia> kostka_polynomial([4,2,1],[3,2,1,1])
t^3 + 2*t^2 + t
```

# Algorithm
The computation here is not based on the above formula but on an explicit description due to Kirillovâ€“Reshetikhin in "The Bethe ansatz and the combinatorics of Young tableaux", J. Sov. Math. 41 (1988) 925., which is summarized by Doreyâ€“Tongaâ€“Turner in "[A matrix model for WZW](https://inspirehep.net/files/ab8568896dcce9b3115b76dc9d096da4)" (Apendix B). Namely:

```math
K_{Î»,Î¼}(t)=âˆ‘_{\\{v\\}}âˆ�_{K=1}^{l(Î»)-1}âˆ�_{nâ‰¥1}
\\begin{bmatrix}
â„™_n^{(K)}+m_n(v^{(K)})\\\\
m_n(v^{(K)})
\\end{bmatrix}_t \\;,
```

where the sum is over all admissible configurations ``\\{v\\}`` i.e. sequences of partitions ``v^{(K)}`` with

```math
v^{(0)}=Î¼ \\text{ ,\\hspace{2mm} } |v^{(K)}|=âˆ‘_{jâ‰¥K+1}Î»_j \\text{\\hspace{5mm} and \\hspace{5mm}} â„™_n^{(K)}â‰¥0	\\text{\\hspace{2mm} for all \\hspace{1mm}} n>0, K=0,1,â€¦l(Î»)
```

```math
\\begin{aligned}
& â„™_n^{(K)} := âˆ‘_{jâ‰¥1}\\left[\\min\\left(n,v_j^{(K+1)}\\right) - 2â‹…\\min\\left(n,v_j^{(K)}\\right) + \\min\\left(n,v_j^{(K-1)}\\right)\\right]
\\\\
& c(v) := âˆ‘_{iâ‰¥1}(i-1)Î¼_i + âˆ‘_{K=1}^{l(Î»)-1}\\left(ğ�•„\\left[v^{(K)}, v^{(K)}\\right] - ğ�•„\\left[v^{(K)}, v^{(K-1)}\\right]\\right)
\\\\
& ğ�•„[Ï�,Îº] := âˆ‘_{i,jâ‰¥1} \\min(Ï�_i,Îº_j)
\\end{aligned}
```

Here, ``\\left[\\genfrac{[}{]}{0pt}{0}{m}{n} \\right]_t`` is the [Gaussian binomial coefficient](https://en.wikipedia.org/wiki/Gaussian_binomial_coefficient).
"""
function kostka_polynomial(Î»::Partition{T}, Î¼::Partition{T}) where T<:Integer
	sum(Î») == sum(Î¼) || throw(ArgumentError("Î» and Î¼ have to be Partitions of the same Integer"))

	R,t = PolynomialRing(ZZ, "t")
	kos_poly = R()

	if length(Î»)==0 || !dominates(Î»,Î¼)
	return R()
	end

	#Here we apply fact 2.12 from https://arxiv.org/pdf/1608.01775.pdf	(AN ITERATIVE FORMULA FOR THE KOSTKA-FOULKES POLYNOMIALS - TIMOTHEE W. BRYAN AND NAIHUAN JING)
	copied = false
	i = 1
	while i<=length(Î») && i<=length(Î¼) && Î»[i]==Î¼[i]
	i += 1
	end
	if i != 1
	if Î» == Î¼
		return R(1)
	end
	Î»_copy = Î»
	Î¼_copy = Î¼
	Î» = Partition{T}(Î»[i:end])
	Î¼ = Partition{T}(Î¼[i:end])
	copied = true
	end

	len_Î» = length(Î») #note that length(v)=Length(Î»)
	len_Î¼ = length(Î¼)

	#calculate the sizes of the partitions from the admissible configurations
	size_v = Partition{T}(zeros(T,len_Î»))
	for k = 1:len_Î»
	for j = k:len_Î»
		size_v[k] += Î»[j]
	end
	end

	#compute the Arrays of available Partitions
	parts = Array{Array{Partition{T},1},1}([[] for i = 1:len_Î»])
	for i = 2:len_Î»
	parts[i] = partitions(size_v[i])
	end
	parts[1] = [Î¼]


	index_max = [length(parts[i]) for i=1:len_Î»] #gives an upper bound to how big index can get
	index = ones(Int, len_Î»)	#will be our iterator with which we will look at all admissible configurations
	empty_partition = Partition{T}([])	#this will act as a placeholder during the Algorithm
	pointer = 2	#this will index the part of the configuration we are currently looking at
	v = Array{Partition{T},1}([empty_partition for i=1:len_Î»])
	v[1] = Î¼

	#returns P_n^(k)
	function vacancy(n,k)
	res = 0

	for vj in v[k-1]
		res += min(n,vj)
	end

	if k < len_Î»
		for vj in v[k+1]
		res += min(n,vj)
		end
	end

	for vj in v[k]
		res -= 2*min(n,vj)
	end

	return res
	end

	#Here the actual Algorithm begins
	while pointer > 1

	if pointer == len_Î» + 1
		summand = t^charge(v)
		divisor = one(R)
		for k = 2:len_Î»
		m = partition_to_partcount(v[k])
		for n = 1:length(m)
			if m[n] > 0
			P_nk = vacancy(n,k)
			if P_nk > 0
				#multiply by the q-binomial coefficient [P_nk+m[n] ; m[n]]
				high = P_nk + m[n]
				for i = 0:m[n]-1
				mul!(summand, summand, (1-t^(high-i)))
				mul!(divisor, divisor, (1-t^(i+1)))
				end
			end
			end
		end
		end
		addeq!(kos_poly, divexact(summand, divisor))
		pointer -= 1
		index[pointer] += 1

	else
		if index[pointer] > index_max[pointer]
		v[pointer] = empty_partition
		pointer -= 1
		index[pointer] += 1
		else
		v[pointer] = parts[pointer][index[pointer]]
		n = 1
		if pointer > 2
			max_n = max(v[pointer-2][1], v[pointer-1][1], v[pointer][1])
		else
			max_n = max(v[pointer-1][1], v[pointer][1])
		end
		while n <= max_n
			if pointer > 2 && vacancy(n,pointer-1) < 0 || vacancy(n,pointer) < -getindex_safe(size_v,pointer+1)
			index[pointer] += 1
			if index[pointer] > index_max[pointer]
				break
			else
				v[pointer] = parts[pointer][index[pointer]]
			end
			n = 0
			end
			n += 1
		end
		if n == max_n + 1
			pointer += 1
			if pointer <= len_Î»
			index[pointer] = 1
			end
		end
		end
	end
	end

	if copied
	Î» = Î»_copy
	Î¼ = Î¼_copy
	end

	return kos_poly
end


function kostka_polynomial(Î»::Array{Int,1}, Î¼::Array{Int,1})
	 return kostka_polynomial(Partition(Î»), Partition(Î¼))
end


"""
	charge(config::Array{Partition{T},1})

The **charge** ``c`` of an admissible configuration ``v:=`` ```config```, is defined by:

```math
\\begin{aligned}
& c(v) := âˆ‘_{iâ‰¥1}(i-1)v^{(0)}_i + âˆ‘_{K=1}^{l(Î»)-1}\\left(ğ�•„\\left[v^{(K)}, v^{(K)}\\right] - ğ�•„\\left[v^{(K)}, v^{(K-1)}\\right]\\right)
\\\\
& ğ�•„[Ï�,Îº] := âˆ‘_{i,jâ‰¥1} \\min(Ï�_i,Îº_j)
\\end{aligned}
```
"""
function charge(config::Array{Partition{T},1}) where T<:Integer
	function M(p::Partition{T}, k::Partition{T})
	res = 0
	for i in p
		for j in k
		if i > j
			res += j
		else
			res += i
		end
		end
	end
	return res
	end

	#M(p) := M(p,p)
	function M(p::Partition{T})
	res = 0
	for i = 1:length(p)
		res+= p[i]*(2*i-1)
	end
	return res
	end

	c = 0
	#n[Î¼]
	for i = 2:length(config[1])
	c += (i-1)*config[1][i]
	end

	for k = 2:length(config)
	c += M(config[k])-M(config[k],config[k-1])	 #first call of M(c,c) could be improved/specialized
	end
	return c
end



"""
	charge(Tab::Tableau)

returns the **charge** of a Tableau ```Tab``` which is defined by the charge of itÂ´s **reading word**.
```
charge(Tab):=charge(reading_word(Tab))
```
"""
function charge(Tab::Tableau)
	return charge(reading_word(Tab))
end


"""
	charge(word::Array{Int,1},standard=false::Bool)

This returns the **charge** of the Tableau corresponding to the reading word ```word```.

We call a word **standard**, if each of its letters or numbers only appears once in it. If you are shure that your word is standard, you can set this argument to ``true``, to improve the efficiency.

This Algorithm is based on the Algorithm following Example 7.3 in "[Hall-Littlewood Functions and Kostka-Foulkes Polynomials in Representation Theory](https://www.emis.de/journals/SLC/opapers/s32leclerc.pdf)", J. DÃ©sarmÃ©nien, B. Leclerc and J.-Y. Thibon.
"""
function charge(word::Array{Int,1},standard=false::Bool)
	c = 0
	if !standard
	baseword = copy(word)
	while length(baseword) > 0
		indices = Array{Int,1}()
		a = 1
		maxim = maximum(baseword)
		while a <= maxim
		emptyrun = true
		for i = length(baseword):-1:1
			if baseword[i] == a
			push!(indices,i)
			a += 1
			emptyrun = false
			end
		end
		if emptyrun
			a += 1
		end
		end
		sort!(indices)
		w = baseword[indices]		 #splice!(baseword,indices)
		deleteat!(baseword,indices) #
		c += charge(w,true)
	end
	else
	index = 0
	pointer = argmin(word)
	maxim = maximum(word)
	while word[pointer]!=maxim
		nextmin = word[pointer]+1
		while findfirst(==(nextmin),word) == nothing
		nextmin += 1
		end
		if findfirst(==(nextmin),word) > pointer
		index += 1
		end
		c += index
		pointer = findfirst(==(nextmin),word)
	end
	end
	return c
end
