# Examples

## Computing immanants, aka Kostant's relation

We can –-numericaly–- verify Kostant's relation. For instance, computing immanant $\lambda = (2,1,0)$ using group functions. You need to install our package Immanant.jl

  add https://github.com/davidamaro/Immanants.jl

Now you can run the test:

  julia> using RandomMatrices, Test, GroupFunctions
julia> using Immanants
julia> mat = rand(Haar(2), 3)
julia> pt_1 = GTPattern([[2,1,0], [2,0],[1]],[1])
julia> pt_2 = GTPattern([[2,1,0], [1,1],[1]],[1])
julia> suma = group_function([2,1,0], pt_1, pt_1, mat)+ group_function([2,1,0], pt_2, pt_2, mat)
julia> @test suma ≈ immanant(Partition([2,1]), mat)
>true

## Legendre's test

First, you need to define a function that generates a $\mathrm{SU}(n)$ block

function bloquesun(tamano::Int64, pos::Int64, angulos::Tuple{Float64,Float64,Float64})
@assert pos < tamano && pos > 0

base = zeros(Complex{Float64}, tamano, tamano)
α,β,γ = angulos
base[pos,pos] = exp(-im*(α+γ))*cos(β/2)
base[pos,pos+1] = -exp(-im*(α-γ))*sin(β/2)
base[pos+1,pos] = exp(im*(α-γ))*sin(β/2)
base[pos+1,pos+1] = exp(im*(α+γ))*cos(β/2)
for i in 1:tamano
if i == pos || i == pos + 1
continue
end
base[i,i] = 1.0
end

return base
end

Now you can type the folling to verify the identity between $L=0$ states and a sum of Legendre polynomials:

  julia> using GroupFunctions, Test
julia> α1,β1,γ1 = rand(Float64,3)
julia> xx=bloquesun(3,1,(α1,β1,γ1))
julia> α2,β2 = rand(Float64,2)
julia> yy=bloquesun(3,2,(α2,β2,α2))
julia> α3,β3,γ3 = rand(Float64,3)
julia> zz=bloquesun(3,1,(α3,β3,γ3))

julia> uu = xx*yy*zz

julia> edo = GTPattern([[2,1,0],[1,1], [1]], [1])
julia> @test group_function([2,1,0],edo, edo, uu) ≈ (1/4)*(1+3*cos(β2))
>true

julia> edo = GTPattern([[4,2,0],[2,2], [2]], [2])
julia> @test group_function([4,2,0],edo, edo, uu) ≈ (1/3^2)*(1+3*cos(β2)+5*(1/2)*(3*cos(β2)^2 - 1))
>true

## Sum rules

Based on the content of Sum rules in multiphoton coincidence rates by David Amaro-Alcalá, Dylan Spivak and Hubert de Guise DOI.

First, we create the unitary matrix (bloquesun is defined above):

julia> α1,β1,γ1 = rand(Float64,3)
julia> xx=bloquesun(4,1,(α1,β1,γ1))
julia> α2,β2 = rand(Float64,2)
julia> yy=bloquesun(4,2,(α2,β2,α2))
julia> α3,β3,γ3 = rand(Float64,3)
julia> zz=bloquesun(4,1,(α3,β3,γ3))
julia> α4,β4 = rand(Float64,3)
julia> xx2=bloquesun(4,3,(α4,β4,α4))
julia> α5,β5 = rand(Float64,2)
julia> yy2=bloquesun(4,2,(α5,β5,α5))
julia> α6,β6,γ6 = rand(Float64,3)
julia> zz2=bloquesun(4,1,(α6,β6,γ6))

Then, we compute irrep states:

julia> welcome = basis_states([3,0,0,0])

Now, we filter the states with first occupation number equals 1

julia> filteredstates = filter(x -> pweight(x)[1] == 1, welcome)

Next, create the actual $SU(4)$ and $SU(4) / SU(3)$ matrices

julia> mat4 = zz2*yy2*xx2*zz*yy*xx
julia> mat4quotient = zz2*yy2*xx2#*zz#*yy#*xx

Carry on the sum

julia> total = 0.0
julia> for edo in filteredstates
total+=abs( group_function([3,0,0,0], filteredstates[2 #=5 3 2 =#], edo, mat4) )^2
end

julia> totalquotient = 0.0
julia> for edo in filteredstates
totalquotient += abs( group_function([3,0,0,0], filteredstates[2 #=5 3 2=#], edo, mat4quotient) )^2
end

Finally, check that total and totalquotient are the same

julia> total ≈ totalquotient
true

### Sum rules running in parallel

Start crearing 20 workers (we assume a fresh Julia session),

using Distributed
addprocs(20)

@everywhere using GroupFunctions
using RandomMatrices, Immanants

Let's start with an example for $\mathrm{SU}6$ with $5$ photons,

α1,β1,γ1 = rand(Float64,3)
xx=bloquesun(6,1,(α1,β1,γ1))
α2,β2 = rand(Float64,2)
yy=bloquesun(6,2,(α2,β2,α2))
α3,β3,γ3 = rand(Float64,3)
zz=bloquesun(6,1,(α3,β3,γ3))

α4,β4 = rand(Float64,3)
xx2=bloquesun(6,3,(α4,β4,α4))
α5,β5 = rand(Float64,2)
yy2=bloquesun(6,2,(α5,β5,α5))
α6,β6,γ6 = rand(Float64,3)
zz2=bloquesun(6,1,(α6,β6,γ6))

α7,β7 = rand(Float64,3)
xx3=bloquesun(6,4,(α7,β7,α7))
α8,β8 = rand(Float64,2)
yy3=bloquesun(6,3,(α8,β8,α8))
α9,β9 = rand(Float64,2)
zz3=bloquesun(6,2,(α9,β9,α9))
α10,β10,γ10 = rand(Float64,3)
ww3=bloquesun(6,1,(α10,β10,γ10))

α11,β11 = rand(Float64,3)
xx4=bloquesun(6,5,(α11,β11,α11))
α12,β12 = rand(Float64,2)
yy4=bloquesun(6,4,(α12,β12,α12))
α13,β13 = rand(Float64,3)
zz4=bloquesun(6,3,(α13,β13,α13))
α14,β14 = rand(Float64,3)
ww4=bloquesun(6,2,(α14,β14,α14))
α15,β15,γ15 = rand(Float64,3)
uu4=bloquesun(6,1,(α15,β15,γ15))

mat6   = uu4*ww4*zz4*yy4*xx4*ww3*zz3*yy3*xx3*zz2*yy2*xx2*zz*yy*xx
mat6c1 = uu4*ww4*zz4*yy4*xx4#*ww3#*zz3#*yy3#*xx3#*zz2#*yy2#*xx2#*zz#*yy#*xx

Then, all the code necessary to create the states,

irrep = [5,0,0,0,0,0];
welcome = basis_states(irrep);
filteredstates = filter(x -> pweight(x)[1] == 1, welcome);
papapapa =filter(x -> maximum(pweight(x)) < 2, filteredstates)

Finally, we run both sums: using the complete and the quotient matrices

@time p1 = @distributed (+) for edo in filteredstates
abs( group_function(irrep, papapapa[2 #=5 3 2=#], edo, mat6c1) )^2
end

@time p2 = @distributed (+) for edo in filteredstates
abs( group_function(irrep, papapapa[2 #=5 3 2=#], edo, mat6) )^2
end

p2 ≈ p1