Littlewood-Richardson rule for skew Schur

LittlewoodRichardsonRule00.R

insertWith <- function(f, mp, key, value) {
  if(key %in% names(mp)) {
    mp[[key]] <- f(value, mp[[key]])
  } else {
    mp[[key]] <- value
  }
  mp
}

zip3 <- function(v1, v2, v3) {
  lapply(1L:length(v1), function(i) {
    c(v1[i], v2[i], v3[i])
  })
}
# _lrRule :: Partition -> Partition -> Map Partition Int
# _lrRule plam@(Partition lam) pmu@(Partition mu0) = 
#   if not (pmu `isSubPartitionOf` plam) 
#     then Map.empty
#     else foldl' f Map.empty [ nu | (nu,_) <- fillings n diagram ]
#   where
#     f old nu = Map.insertWith (+) (Partition nu) 1 old
#     diagram  = [ (i,j) | (i,a,b) <- reverse (zip3 [1..] lam mu) , j <- [b+1..a] ]    
#     mu       = mu0 ++ repeat 0
#     n        = sum' $ zipWith (-) lam mu 

lrRule <- function(lambda, mu) {
  f <- function(old, nu) {
    insertWith(`+`, old, toString(nu), 1L)
  }
  l <- length(lambda)
  mu <- c(mu, rep(0L, l - length(mu)))
  Liab <- rev(zip3(seq_len(l), lambda, mu))
  diagram <- do.call(rbind, do.call(c, lapply(Liab, function(iab) {
    i <- iab[1L]
    a <- iab[2L]
    b <- iab[3L]
    jvec <- if(b + 1L <= a) (b + 1L):a else integer(0L)
    lapply(jvec, function(j) {
      c(i, j)
    })
  })))
  n <- sum(lambda - mu)
  Lnu <- lapply(fillings(n, diagram), `[[`, 1L)
  Reduce(f, Lnu, init = integer(0L))
}

# fillings :: Int -> Diagram -> [Filling]
# fillings _ []                   = [ ([],[]) ]
# fillings n diagram@((x,y):rest) = concatMap (nextLetter lower upper) (fillings (n-1) rest) where
#   upper = case findIndex (==(x  ,y+1)) diagram of { Just j -> n-j ; Nothing -> 0 }
#   lower = case findIndex (==(x-1,y  )) diagram of { Just j -> n-j ; Nothing -> 0 }

fillings <- function(n, diagram) {
  if(nrow(diagram) == 0L) {
    list(list(integer(0L), integer(0L)))
  } else {
    xy <- diagram[1L, ]
    x <- xy[1L]
    y <- xy[2L]
    rest <- diagram[-1L, , drop = FALSE]
    diagram <- apply(diagram, 1L, toString)
    upper <- n - (match(toString(c(x, y + 1L)), diagram, nomatch = n + 1L) - 1L)
    lower <- n - (match(toString(c(x - 1L, y)), diagram, nomatch = n + 1L) - 1L)
    L <- lapply(fillings(n - 1L, rest), function(filling) {
      nextLetter(lower, upper, filling)
    })
    do.call(c, L)
  }
}

# nextLetter :: Int -> Int -> Filling -> [Filling]
# nextLetter lower upper (nu,lpart) = stuff where
#   stuff = [ ( incr i shape , lpart++[i] ) | i<-nlist ] 
#   shape = nu ++ [0] 
#   lb = if lower>0
#     then lpart !! (lower-1)
#     else 0
#   ub = if upper>0 
#     then min (length shape) (lpart !! (upper-1))  
#     else      length shape
# 
#   nlist = filter (>0) $ map f [lb+1..ub] 
#   f j   = if j==1 || shape!!(j-2) > shape!!(j-1) then j else 0

nextLetter <- function(lower, upper, filling) {
  nu <- filling[[1L]]
  lpart <- filling[[2L]]
  shape <- c(nu, 0L)
  lb <- if(lower > 0L) lpart[lower] else 0L
  ub <- if(upper > 0L) min(length(shape), lpart[upper]) else length(shape)
  f <- function(j) {
    if(j == 1L || shape[j-1L] > shape[j]) j else 0L
  }
  v <- vapply((lb+1L):ub, f, integer(1L))
  nlist <- v[v > 0L]
  lapply(nlist, function(i) {
    list(incr(i, shape), c(lpart, i))
  })
}

#   incr :: Int -> [Int] -> [Int]
#   incr i (x:xs) = case i of
#     0 ->         finish (x:xs)
#     1 -> (x+1) : finish xs
#     _ -> x     : incr (i-1) xs
#   incr _ []     = []
incr <- function(i, xxs) {
  if(length(xxs) == 0L) {
    integer(0L)
  } else {
    if(i == 0L) {
      finish(xxs)
    } else if(i == 1L) {
      c(xxs[1L] + 1L, finish(xxs[-1L]))
    } else {
      c(xxs[1L], incr(i - 1L, xxs[-1L]))
    }
  }
}

#   finish :: [Int] -> [Int]
#   finish (x:xs) = if x>0 then x : finish xs else []    
#   finish []     = [] 
finish <- function(xxs) {
  if(length(xxs) == 0L) {
    integer(0L)
  } else {
    x <- xxs[1L]
    if(x > 0L) {
      c(x, finish(xxs[-1L]))
    } else {
      integer(0L)
    }
  }
}