Braid groups

ASCII00.R

# hSepString :: HSep -> String
# hSepString hsep = case hsep of
#   HSepEmpty    -> ""
#   HSepSpaces k -> replicate k ' '
#   HSepString s -> s
hSepString <- function(hsep) {
  sep <- attr(hsep, "sep")
  switch(
    sep,
    "empty" = "",
    "spaces" = paste0(rep(" ", hsep), collapse = ""),
    "string" = hsep
  )
}

vSepString <- function(vsep) {
  sep <- attr(vsep, "sep")
  switch(
    sep,
    "empty" = character(0L),
    "spaces" = paste0(rep(" ", vsep), collapse = ""),
    "string" = vsep
  )
}

vSepSpaces <- function(k) {
  out <- k
  attr(out, "sep") <- "spaces"
  out
}

vSepSize <- function(vsep) {
  nchar(vSepString(vsep))
}

HSepEmpty <- function() {
  out <- "."
  attr(out, "sep") <- "empty"
  out
}

ASCII <- function(x, y, lines) {
  list("x" = x, "y" = y, "lines" = lines)
}

asciiLines <- function(ascii) {
  ascii[["lines"]]
}

# -- | Extends an ASCII figure with spaces vertically to the given height.
# -- Note: the alignment is the alignment of the original picture in the new bigger picture!
# vExtendTo :: VAlign -> Int -> ASCII -> ASCII
# vExtendTo valign n0 rect@(ASCII (x,y) ls) = vExtendWith valign (max n0 y - y) rect
vExtendTo <- function(valign, n0, rect) {
  y <- rect[["y"]]
  vExtendWith(valign, max(n0, y) - y, rect)
}


# -- | Extend vertically with the given number of empty lines.
# vExtendWith :: VAlign -> Int -> ASCII -> ASCII
# vExtendWith valign d (ASCII (x,y) ls) = ASCII (x,y+d) (f ls) where
#   f ls = case valign of
#     VTop     -> ls ++ replicate d emptyline
#     VBottom  -> replicate d emptyline ++ ls
#     VCenter  -> replicate a emptyline ++ ls ++ replicate (d-a) emptyline
#   a = div d 2
#   emptyline = replicate x ' '
vExtendWith <- function(valign, d, rect) {
  x <- rect[["x"]]
  y <- rect[["y"]]
  lines <- rect[["lines"]]
  emptyLine <- paste0(rep(" ", x), collapse = "")
  f <- function(ls) {
    a <- d %/% 2L
    switch(
      valign,
      "Vtop" = c(ls, rep(emptyLine, d)),
      "Vbottom" = c(rep(emptyLine, d), ls),
      "Vcenter" = c(rep(emptyLine, a), ls, rep(emptyLine, d - a))
    )
  }
  ASCII(x, y + d, f(lines))
}

# -- | Extends an ASCII figure with spaces horizontally to the given width.
# -- Note: the alignment is the alignment of the original picture in the new bigger picture!
# hExtendTo :: HAlign -> Int -> ASCII -> ASCII
# hExtendTo halign n0 rect@(ASCII (x,y) ls) = hExtendWith halign (max n0 x - x) rect
hExtendTo <- function(halign, n0, rect) {
  x <- rect[["x"]]
  hExtendWith(halign, max(n0, x) - x, rect)
}

# -- | Extend horizontally with the given number of spaces.
# hExtendWith :: HAlign -> Int -> ASCII -> ASCII
# hExtendWith alignment d (ASCII (x,y) ls) = ASCII (x+d,y) (map f ls) where
#   f l = case alignment of
#     HLeft   -> l ++ replicate d ' '
#     HRight  -> replicate d ' ' ++ l
#     HCenter -> replicate a ' ' ++ l ++ replicate (d-a) ' '
#   a = div d 2
hExtendWith <- function(halign, d, rect) {
  x <- rect[["x"]]
  y <- rect[["y"]]
  lines <- rect[["lines"]]
  f <- function(l) {
    a <- d %/% 2L
    switch(
      halign,
      "Hleft" = paste0(c(l, rep(" ", d)), collapse = ""),
      "Hright" = paste0(c(rep(" ", d), l), collapse = ""),
      "Hcenter" = paste0(c(rep(" ", a), l, rep(" ", d - a)), collapse = "")
    )
  }
  ASCII(x + d, y, vapply(lines, f, character(1L)))
}


# asciiFromLines :: [String] -> ASCII
# asciiFromLines ls = ASCII (x,y) (map f ls) where
#   y   = length ls
#   x   = maximum (map length ls)
#   f l = l ++ replicate (x - length l) ' '
asciiFromLines <- function(ls) {
  y <- length(ls)
  x <- max(vapply(ls, nchar, integer(1L)))
  f <- function(l) {
    paste0(c(l, rep(" ", x - nchar(l))), collapse = "")
  }
  ASCII(x, y, vapply(ls, f, character(1L)))
}

# -- | Horizontal concatenation, top-aligned, no separation
# hCatTop :: [ASCII] -> ASCII
# hCatTop = hCatWith VTop HSepEmpty
#
hCatTop <- function(asciis) {
  hCatWith("Vtop", HSepEmpty(), asciis)
}

# -- | General horizontal concatenation
# hCatWith :: VAlign -> HSep -> [ASCII] -> ASCII
# hCatWith valign hsep rects = ASCII (x',maxy) final where
#   n    = length rects
#   maxy = maximum [ y | ASCII (_,y) _ <- rects ]
#   xsz  =         [ x | ASCII (x,_) _ <- rects ]
#   sep   = hSepString hsep
#   sepx  = length sep
#   rects1 = map (vExtendTo valign maxy) rects
#   x' = sum' xsz + (n-1)*sepx
#   final = map (intercalate sep) $ transpose (map asciiLines rects1)
hCatWith <- function(valign, hsep, rects) {
  n <- length(rects)
  maxy <- max(vapply(rects, `[[`, integer(1L), "y"))
  xsz <- vapply(rects, `[[`, integer(1L), "x")
  sep <- hSepString(hsep)
  sepx <- length(sep)
  rects1 <- lapply(rects, function(rect) {
    vExtendTo(valign, maxy, rect)
  })
  x2 <- sum(xsz) + (n - 1L) * sepx
  M <- do.call(rbind, lapply(rects1, asciiLines))
  final <- apply(M, 2L, paste0, collapse = sep)
  ASCII(x2, maxy, final)
}

# intercalate : List A → List (List A) → List A
# intercalate xs []         = []
# intercalate xs (ys ∷ [])  = ys
# intercalate xs (ys ∷ yss) = ys ++ xs ++ intercalate xs yss
intercalate <- function(sep, x) {
  if(length(x) == 0L) {
    list()
  } else if(length(x) == 1L) {
    x[1L]
  } else {
    unlist(c(x[1L], list(sep), intercalate(sep, x[-1L])))
  }
}

#intercalate(c("a", "b"), list(c("xxx", "yyy"), c("zzz", "ooo"), c("uuu", "vvv")))

# -- | General vertical concatenation
# vCatWith :: HAlign -> VSep -> [ASCII] -> ASCII
# vCatWith halign vsep rects = ASCII (maxx,y') final where
#   n    = length rects
#   maxx = maximum [ x | ASCII (x,_) _ <- rects ]
#   ysz  =         [ y | ASCII (_,y) _ <- rects ]
#   sepy    = vSepSize vsep
#   fullsep = transpose (replicate maxx $ vSepString vsep) :: [String]
#   rects1  = map (hExtendTo halign maxx) rects
#   y'    = sum' ysz + (n-1)*sepy
#   final = intercalate fullsep $ map asciiLines rects1
vCatWith <- function(halign, vsep, rects) {
  n <- length(rects)
  maxx <- max(vapply(rects, `[[`, integer(1L), "x"))
  ysz <- vapply(rects, `[[`, integer(1L), "y")
  sepy <- vSepSize(vsep)
  vsepstring <- vSepString(vsep)
  fullsep <- apply(
    rbind(
      vapply(strsplit(vsepstring, "")[[1L]], rep, character(maxx), times = maxx)
    ),
    2L, paste0, collapse = ""
  )
  rects1 <- lapply(rects, function(rect) {
    hExtendTo(halign, maxx, rect)
  })
  y2 <- sum(ysz) + (n - 1L) * sepy
  final <- intercalate(fullsep, lapply(rects1, asciiLines))
  ASCII(maxx, y2, final)
}

# -- | A box simply filled with the given character
# filledBox :: Char -> (Int,Int) -> ASCII
# filledBox c (x0,y0) = asciiFromLines $ replicate y (replicate x c) where
#   x = max 0 x0
#   y = max 0 y0
#
# -- | A box of spaces
# transparentBox :: (Int,Int) -> ASCII
# transparentBox = filledBox ' '
filledBox <- function(c, x0, y0) {
  x <- max(0L, x0)
  y <- max(0L, y0)
  asciiFromLines(rep(paste0(rep(c, x), collapse = ""), y))
}

transparentBox <- function(x0, y0) {
  filledBox(" ", x0, y0)
}

asciiShow <- function(x) {
  asciiFromLines(x)
}
# horizBraidASCII' :: KnownNat n => Bool -> Braid n -> ASCII
# horizBraidASCII' flipped braid@(Braid gens) = final where
#
#   n = numberOfStrands braid
#
#   final        = vExtendWith VTop 1 $ hCatTop allBlocks
#   allBlocks    = prelude ++ middleBlocks ++ epilogue
#   prelude      = [ numberBlock   , spaceBlock , beginEndBlock ]
#   epilogue     = [ beginEndBlock , spaceBlock , numberBlock'  ]
#   middleBlocks = map block gens
#
#   block g = case g of
#     Sigma    i -> block' i $ if flipped then over  else under
#     SigmaInv i -> block' i $ if flipped then under else over
#
#   block' i middle = asciiFromLines $ drop 2 $ concat
#                   $ replicate a horiz ++ [space3, middle] ++ replicate b horiz
#     where
#       (a,b) = if flipped then (n-i-1,i-1) else (i-1,n-i-1)
#
#   spaceBlock    = transparentBox (1,n*3-2)
#   beginEndBlock = asciiFromLines $ drop 2 $ concat $ replicate n horiz
#   numberBlock   = mkNumbers [1..n]
#   numberBlock'  = mkNumbers $ P.fromPermutation $ braidPermutation braid
#
#   mkNumbers :: [Int] -> ASCII
#   mkNumbers list = vCatWith HRight (VSepSpaces 2) $ map asciiShow
#                  $ (if flipped then reverse else id) $ list
#
#   under  = [ "\\ /" , " / "  , "/ \\" ]
#   over   = [ "\\ /" , " \\ " , "/ \\" ]
#   horiz  = [ "   "  , "   "  , "___"  ]
#   space3 = [ "   "  , "   "  , "   "  ]
horizBraidASCII <- function(flipped, braid) {
  under  = c("\\ /" , " / "  , "/ \\")
  over   = c("\\ /" , " \\ " , "/ \\")
  horiz  = c("   "  , "   "  , "___")
  space3 = c("   "  , "   "  , "   ")
  n <- numberOfStrands(braid)
  block2 <- function(i, middle) {
    if(flipped) {
      a <- n - i - 1L
      b <- i - 1L
    } else {
      a <- i - 1L
      b <- n - i - 1L
    }
    x <- c(rep(horiz, a), c(space3, middle), rep(horiz, b))
    asciiFromLines(x[-c(1L, 2L)])
  }
  block <- function(g) {
    i <- g[1L]
    if(g[2L] == 1L) {
      block2(i, if(flipped) over else under)
    } else {
      block2(i, if(flipped) under else over)
    }
  }
  mkNumbers <- function(x) {
    if(flipped) x <- rev(x)
    vCatWith(
      "Hright",
      vSepSpaces(2L),
      lapply(x, asciiShow)
    )
  }
  #   spaceBlock    = transparentBox (1,n*3-2)
  #   beginEndBlock = asciiFromLines $ drop 2 $ concat $ replicate n horiz
  #   numberBlock   = mkNumbers [1..n]
  #   numberBlock'  = mkNumbers $ P.fromPermutation $ braidPermutation braid
  spaceBlock <- transparentBox(1L, 3L*n - 2L)
  beginEndBlock <- asciiFromLines(rep(horiz, n)[-c(1L, 2L)])
  numberBlock <- mkNumbers(1L:n)
  numberBlock2 <- mkNumbers(braidPermutation(braid))
  prelude <- list(numberBlock, spaceBlock, beginEndBlock)
  epilogue <- list(beginEndBlock, spaceBlock, numberBlock2)
  middleBlocks <- lapply(braid, block)
  allBlocks <- c(prelude, middleBlocks, epilogue)
  print(allBlocks)
  final <- vExtendWith("Vtop", 1L, hCatTop(allBlocks))
}

braidGroups00.R

library(maybe)

isPositiveInteger <- function(n) {
  length(n) == 1L && is.numeric(n) && !is.na(n) && n != 0 && floor(n) == n
}

#' @title Artin generators
#' @description A standard Artin generator of a braid: \code{Sigma(i)} 
#'   represents twisting the neighbour strands \code{i} and \code{i+1}, 
#'   such that strand \code{i} goes \emph{under} strand \code{i+1}.
#' 
#' @param i index of the strand, a positive integer
#' 
#' @returns A vector of two integers.
#' @export
#' @name ArtinGenerators
#' @rdname ArtinGenerators
Sigma <- function(i) {
  stopifnot(isPositiveInteger(i))
  c(as.integer(i), 1L)
}

#' @export
#' @rdname ArtinGenerators
SigmaInv <- function(i) {
  stopifnot(isPositiveInteger(i))
  c(as.integer(i), -1L)
}

mkBraid <- function(n, brgens) {
  stopifnot(isPositiveInteger(n))
  out <- brgens
  idx <- vapply(brgens, `[[`, integer(1L), 1L)
  if(any(idx) >= n) {
    stop("Found a generator with a too large index.")
  }
  attr(out, "n") <- as.integer(n)
  class(out) <- "braid"
  out
}

#' @exportS3Method print braid
print.braid <- function(x, ...) {
  idx <- vapply(x, `[[`, integer(1L), 1L)
  signs <- vapply(x, `[[`, integer(1L), 2L)
  print(paste0(vapply(seq_along(idx), function(i) {
    if(signs[i] == 1L) {
      sprintf("sigma_%d", idx[i])
    } else {
      sprintf("sigmaInv_%d", idx[i])
    }
  }, character(1L)), collapse = " "))
  invisible()
}

numberOfStrands <- function(braid) {
  attr(braid, "n")
}

# freeReduceBraidWord :: Braid n -> Braid n
# freeReduceBraidWord (Braid orig) = Braid (loop orig) where
# 
#   loop w = case reduceStep w of
#     Nothing -> w
#     Just w' -> loop w'
#   
#   reduceStep :: [BrGen] -> Maybe [BrGen]
#   reduceStep = go False where    
#     go !changed w = case w of
#       (Sigma    x : SigmaInv y : rest) | x==y   -> go True rest
#       (SigmaInv x : Sigma    y : rest) | x==y   -> go True rest
#       (this : rest)                             -> liftM (this:) $ go changed rest
#       _                                         -> if changed then Just w else Nothing
freeReduceBraidWord <- function(braid) {
  reduceStep <- function(gens) {
    go <- function(changed, w) {
      if(length(w) >= 2L) {
        w1 <- w[[1L]]
        w2 <- w[[2L]]
        x <- w1[1L]
        y <- w2[1L]
        s1 <- w1[2L]
        s2 <- w2[2L]
        if(x == y && ((s1 == 1L && s2 == -1L) || (s1 == -1L && s2 == 1L))) {
          go(TRUE, w[-c(1L, 2L)])
        } else {
          gg <- go(changed, w[-1L])
          if(is_nothing(gg)) {
            nothing()
          } else {
            just(c(list(w1), from_just(gg)))
          }
        }
      } else if(length(w) == 1L) {
        if(changed) {
          just(c(list(w[[1L]]), w))
        } else {
          nothing()
        }
      } else {
        if(changed) {
          just(w)
        } else {
          nothing()
        }
      }
    }
    go(FALSE, gens)
  }
  loop <- function(w) {
    x <- reduceStep(w)
    if(is_nothing(x)) {
      w
    } else {
      loop(from_just(x))
    }
  }
  mkBraid(numberOfStrands(braid), loop(braid))
}

# -- | This is an untyped version of 'braidPermutation'
# _braidPermutation :: Int -> [Int] -> Permutation
# _braidPermutation n idxs = P.uarrayToPermutationUnsafe (runSTUArray action) where
# 
#   action :: forall s. ST s (STUArray s Int Int) 
#   action = do 
#     arr <- newArray_ (1,n) 
#     forM_ [1..n] $ \i -> writeArray arr i i
#     worker arr idxs
#     return arr
#     
#   worker arr = go where
#     go []     = return arr 
#     go (i:is) = do
#       a <- readArray arr  i
#       b <- readArray arr (i+1)
#       writeArray arr  i    b
#       writeArray arr (i+1) a
#       go is
braidPermutation <- function(brain) {
  n <- numberOfStrands(brain)
  idxs <- vapply(brain, `[[`, integer(1L), 1L)
  worker <- function(arr, idxs) {
    if(length(idxs) == 0L) {
      arr
    } else {
      i <- idxs[1L]
      is <- idxs[-1L]
      a <- arr[i]
      b <- arr[i + 1L]
      arr[i] <- b
      arr[i + 1L] <- a
      worker(arr, is)
    }
  }
  worker(1L:n, idxs)
}

# doubleSigma :: KnownNat n => Int -> Int -> Braid (n :: Nat)
# doubleSigma s t = braid where
#   n = numberOfStrands braid
#   braid
#     | s < 1 || s > n   = error "doubleSigma: s index out of range"
#     | t < 1 || t > n   = error "doubleSigma: t index out of range"
#     | s >= t           = error "doubleSigma: s >= t"
#     | otherwise        = Braid $
#        [ Sigma i | i<-[t-1,t-2..s] ] ++ [ SigmaInv i | i<-[s+1..t-1] ]
doubleSigma <- function(n, s, t) {
  stopifnot(isPositiveInteger(n), isPositiveInteger(s), isPositiveInteger(t))
  if(s > n) {
    stop("The `s` index is out of range.")
  }
  if(t > n) {
    stop("The `t` index is out of range.")
  }
  if(s >= t) {
    stop("`s` must be strictly smaller than `t`.")
  }
  if(t - 1L <= s) {
    gens <- lapply((t-1L):s, Sigma)
  } else {
    gens <- lapply((s+1L):(t-1L), SigmaInv)
  }
  mkBraid(n, gens)
}

# -- | The (positive) half-twist of all the braid strands, usually denoted by @Delta@.
# halfTwist :: KnownNat n => Braid n
# halfTwist = braid where
#   braid = Braid $ map Sigma $ _halfTwist n 
#   n     = numberOfStrands braid
# 
# -- | The untyped version of 'halfTwist'
# _halfTwist :: Int -> [Int]
# _halfTwist n = gens where
#   gens  = concat [ sub k | k<-[1..n-1] ]
#   sub k = [ j | j<-[n-1,n-2..k] ]
halfTwist <- function(n) {
  stopifnot(isPositiveInteger(n))
  if(n == 1L) {
    gens <- list()
  } else {
    subs <- lapply(1L:(n-1L), function(k) {
      (n-1L):k
    })
    gens <- lapply(do.call(c, subs), Sigma)
  }
  mkBraid(n, gens)
}

# tau :: KnownNat n => Braid n -> Braid n
# tau :: forall (n :: Nat). KnownNat n => Braid n -> Braid n
# tau braid@(Braid gens) = Braid (map f gens) where
#   n = numberOfStrands braid
#   f (Sigma    i) = Sigma    (n-i)
#   f (SigmaInv i) = SigmaInv (n-i)
tau <- function(braid) {
  n <- numberOfStrands(braid)
  gens <- lapply(braid, function(gen) {
    i <- gen[2L]
    if(i == 1L) {
      Sigma(n - i)
    } else {
      SigmaInv(n - i)
    }
  })
  mkBraid(n, gens)
}

# tauPerm :: Permutation -> Permutation
# tauPerm :: Permutation -> Permutation
# tauPerm perm = P.toPermutationUnsafeN n [ (n+1) - perm !!! (n-i) | i<-[0..n-1] ] where
#   n = P.permutationSize perm
tauPerm <- function(perm) {
  n <- length(perm)
  # ??
}

# -- | The inverse of a braid. Note: we do not perform reduction here,
# -- as a word is reduced if and only if its inverse is reduced.
# inverse :: Braid n -> Braid n
# inverse = Braid . reverse . map invBrGen . braidWord
inverseBraid <- function(braid) {
  n <- numberOfStrands(braid)
  invgens <- lapply(braid, function(gen) {
    c(gen[1L], -gen[2L])
  })
  mkBraid(n, rev(invgens))
}

# -- | Composes two braids, doing free reduction on the result 
# -- (that is, removing @(sigma_k * sigma_k^-1)@ pairs@)
# compose :: Braid n -> Braid n -> Braid n
# compose (Braid gs) (Braid hs) = freeReduceBraidWord $ Braid (gs++hs)
composeTwoBraids <- function(braid1, braid2) {
  n <- numberOfStrands(braid1) 
  if(n != numberOfStrands(braid2)) {
    stop("Unequal numbers of strands.")
  }
  freeReduceBraidWord(mkBraid(n, c(braid1, braid2)))
}

composeManyBraids <- function(braids) {
  ns <- vapply(braids, numberOfStrands, integer(1L))
  n <- ns[1L]
  if(any(ns != n)) {
    stop("Unequal numbers of strands.")
  }
  freeReduceBraidWord(mkBraid(n, do.call(c, braids)))
}

isPureBraid <- function(braid) {
  identical(braidPermutation(braid), seq_len(numberOfStrands(braid)))
}

# -- | A positive braid word contains only positive (@Sigma@) generators.
# isPositiveBraidWord :: KnownNat n => Braid n -> Bool
# isPositiveBraidWord (Braid gs) = all (isPlus . brGenSign) gs 
isPositiveBraidWord <- function(braid) {
  signs <- vapply(braid, `[[`, integer(1L), 2L)
  all(signs == 1L)
}

# -- | We compute the linking numbers between all pairs of strands:
# --
# -- > linkingMatrix braid ! (i,j) == strandLinking braid i j 
# --
# linkingMatrix :: KnownNat n => Braid n -> UArray (Int,Int) Int
# linkingMatrix braid@(Braid gens) = _linkingMatrix (numberOfStrands braid) gens where
# 
# -- | Untyped version of 'linkingMatrix'
# _linkingMatrix :: Int -> [BrGen] -> UArray (Int,Int) Int
# _linkingMatrix n gens = runSTUArray action where
# 
#   action :: forall s. ST s (STUArray s (Int,Int) Int)
#   action = do
#     perm <- newArray_ (1,n) :: ST s (STUArray s Int Int)
#     forM_ [1..n] $ \i -> writeArray perm i i
#     let doSwap :: Int -> ST s ()
#         doSwap i = do
#           a <- readArray perm  i
#           b <- readArray perm (i+1)
#           writeArray perm  i    b
#           writeArray perm (i+1) a
#                
#     mat <- newArray ((1,1),(n,n)) 0 :: ST s (STUArray s (Int,Int) Int)
#     let doAdd :: Int -> Int -> Int -> ST s ()
#         doAdd i j pm1 = do
#           x <- readArray mat (i,j)
#           writeArray mat (i,j) (x+pm1) 
#           writeArray mat (j,i) (x+pm1)
#        
#     forM_ gens $ \g -> do
#       let (sgn,k) = brGenSignIdx g
#       u <- readArray perm  k 
#       v <- readArray perm (k+1)
#       doAdd u v (signValue sgn)
#       doSwap k 
#         
#     return mat
linkingMatrix <- function(braid) {
  n <- numberOfStrands(braid)
  perm <- 1L:n
  doSwap <- function(i) {
    a <- perm[i]
    b <- perm[i+1L]
    perm[i] <<- b
    perm[i+1L] <<- a
    invisible()
  }
  mat <- matrix(0L, nrow = n, ncol = n)
  doAdd <- function(i, j, pm1) {
    x <- mat[i, j]
    mat[i, j] <<- mat[j, i] <<- x + pm1
    invisible()
  }
  for(gen in braid) {
    k <- gen[1L]
    u <- perm[k]
    v <- perm[k+1L]
    doAdd(u, v, gen[2L])
    doSwap(k)
  }
  mat
}

braid <- mkBraid(4, list(Sigma(2), SigmaInv(3)))
linkingMatrix(braid)

# -- | A /permutation braid/ is a positive braid where any two strands cross
# -- at most one, and /positively/. 
# --
# isPermutationBraid :: KnownNat n => Braid n -> Bool
# isPermutationBraid braid = isPositiveBraidWord braid && crosses where
#   crosses     = and [ check i j | i<-[1..n-1], j<-[i+1..n] ] 
#   check i j   = zeroOrOne (lkMatrix ! (i,j)) 
#   zeroOrOne a = (a==1 || a==0)
#   lkMatrix    = linkingMatrix   braid
#   n           = numberOfStrands braid
isPermutationBraid <- function(braid) {
  if(isPositiveBraidWord(braid)) {
    lkMatrix <- linkingMatrix(braid)
    all(lkMatrix[upper.tri(lkMatrix)] %in% c(0L, 1L))
  } else {
    FALSE
  }
}

braidGroups01.R

library(maybe)

isPositiveInteger <- function(n) {
  length(n) == 1L && is.numeric(n) && !is.na(n) && n != 0 && floor(n) == n
}

#' @title Artin generators
#' @description A standard Artin generator of a braid: \code{Sigma(i)} 
#'   represents twisting the neighbour strands \code{i} and \code{i+1}, 
#'   such that strand \code{i} goes \emph{under} strand \code{i+1}.
#' 
#' @param i index of the strand, a positive integer
#' 
#' @returns A vector of two integers.
#' @export
#' @name ArtinGenerators
#' @rdname ArtinGenerators
Sigma <- function(i) {
  stopifnot(isPositiveInteger(i))
  c(as.integer(i), 1L)
}

#' @export
#' @rdname ArtinGenerators
SigmaInv <- function(i) {
  stopifnot(isPositiveInteger(i))
  c(as.integer(i), -1L)
}

#' @title Make a braid
#' @description Make a braid.
#'
#' @param n number of strands, an integer, at least 2
#' @param brgens list of generators obtained with \code{\link{Sigma}} or 
#'   \code{\link{SigmaInv}}
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' mkBraid(4, list(Sigma(2), SigmaInv(3)))
mkBraid <- function(n, brgens) {
  stopifnot(isPositiveInteger(n), n >= 2)
  out <- brgens
  idx <- vapply(brgens, `[[`, integer(1L), 1L)
  if(any(idx) >= n) {
    stop("Found a generator with a too large index.")
  }
  attr(out, "n") <- as.integer(n)
  class(out) <- "braid"
  out
}

#' @exportS3Method print braid
print.braid <- function(x, ...) {
  idx <- vapply(x, `[[`, integer(1L), 1L)
  signs <- vapply(x, `[[`, integer(1L), 2L)
  print(paste0(vapply(seq_along(idx), function(i) {
    if(signs[i] == 1L) {
      sprintf("sigma_%d", idx[i])
    } else {
      sprintf("sigmaInv_%d", idx[i])
    }
  }, character(1L)), collapse = " "))
  invisible()
}

#' @title Number of strands
#' @description The number of strands of a braid.
#'
#' @param braid a \code{braid} object created with \code{\link{mkBraid}}
#'
#' @return An integer.
#' @export
#'
#' @examples
numberOfStrands <- function(braid) {
  attr(braid, "n")
}

#' @title Free reduction of a braid
#' @description Applies free reduction to a braid, i.e. removes pairs of 
#'   consecutive generators inverse of each other.
#'
#' @param braid a \code{braid} object created with \code{\link{mkBraid}}
#'
#' @return A \code{braid} object.
#' @export
#' @importFrom maybe just nothing is_nothing from_just
#'
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' freeReduceBraidWord(braid)
freeReduceBraidWord <- function(braid) {
  reduceStep <- function(gens) {
    go <- function(changed, w) {
      if(length(w) >= 2L) {
        w1 <- w[[1L]]
        w2 <- w[[2L]]
        x <- w1[1L]
        y <- w2[1L]
        s1 <- w1[2L]
        s2 <- w2[2L]
        if(x == y && ((s1 == 1L && s2 == -1L) || (s1 == -1L && s2 == 1L))) {
          go(TRUE, w[-c(1L, 2L)])
        } else {
          gg <- go(changed, w[-1L])
          if(is_nothing(gg)) {
            nothing()
          } else {
            just(c(list(w1), from_just(gg)))
          }
        }
      } else if(length(w) == 1L) {
        if(changed) {
          just(c(list(w[[1L]]), w))
        } else {
          nothing()
        }
      } else {
        if(changed) {
          just(w)
        } else {
          nothing()
        }
      }
    }
    go(FALSE, gens)
  }
  loop <- function(w) {
    x <- reduceStep(w)
    if(is_nothing(x)) {
      w
    } else {
      loop(from_just(x))
    }
  }
  mkBraid(numberOfStrands(braid), loop(braid))
}

#' @title Braid permutation
#' @description Returns the left-to-right permutation associated to a braid. 
#'
#' @param braid a \code{braid} object created with \code{\link{mkBraid}}
#'
#' @return A permutation.
#' @export
#'
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' braidPermutation(braid)
braidPermutation <- function(braid) {
  n <- numberOfStrands(braid)
  idxs <- vapply(braid, `[[`, integer(1L), 1L)
  worker <- function(arr, idxs) {
    if(length(idxs) == 0L) {
      arr
    } else {
      i <- idxs[1L]
      is <- idxs[-1L]
      a <- arr[i]
      b <- arr[i + 1L]
      arr[i] <- b
      arr[i + 1L] <- a
      worker(arr, is)
    }
  }
  worker(1L:n, idxs)
}

#' @title Double generator
#' @description Generator \code{sigma_{s,t}} in the Birman-Ko-Lee new 
#'   presentation. It twistes the strands \code{s} and \code{t} whie going over 
#'   all other strands (for \code{t=s+1}, this is \code{sigma_s}).
#' 
#' @param n number of strands, integer \code{>= 2}
#' @param s,t indices of two strands, \code{s < t} 
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' doubleSigma(5, 1, 3)
doubleSigma <- function(n, s, t) {
  stopifnot(isPositiveInteger(n), isPositiveInteger(s), isPositiveInteger(t))
  stopifnot(n >= 2)
  if(s > n) {
    stop("The `s` index is out of range.")
  }
  if(t > n) {
    stop("The `t` index is out of range.")
  }
  if(s >= t) {
    stop("`s` must be strictly smaller than `t`.")
  }
  if(t - 1L <= s) {
    gens <- lapply((t-1L):s, Sigma)
  } else {
    gens <- lapply((s+1L):(t-1L), SigmaInv)
  }
  mkBraid(n, gens)
}

#' @title Half-twist
#' @description The (positive) half-twist of all the braid strands, usually 
#'   denoted by \eqn{\Delta}.
#' 
#' @param n number of strands, integer \code{>= 2}
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' halfTwist(4)
halfTwist <- function(n) {
  stopifnot(isPositiveInteger(n), n >= 2)
  # if(n == 1L) {
  #   gens <- list()
  # } else {
    subs <- lapply(1L:(n-1L), function(k) {
      (n-1L):k
    })
    gens <- lapply(do.call(c, subs), Sigma)
  # }
  mkBraid(n, gens)
}

#' @title Inner automorphism
#' @description The inner automorphism defined by 
#'   \eqn{\tau X = \Delta^{-1} X \Delta}, where \eqn{\Delta} is the 
#'   positive half-twist; it send each generator \eqn{\sigma_j} to 
#'   \eqn{\sigma_{n-j}}.
#'
#' @param braid a \code{braid} object created with \code{\link{mkBraid}}
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' tau(braid)
tau <- function(braid) {
  n <- numberOfStrands(braid)
  gens <- lapply(braid, function(gen) {
    i <- gen[2L]
    if(i == 1L) {
      Sigma(n - i)
    } else {
      SigmaInv(n - i)
    }
  })
  mkBraid(n, gens)
}

#' @title Inverse braid
#' @description The inverse of a braid (without performing reduction).
#' 
#' @param braid a \code{braid} object created with \code{\link{mkBraid}}
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' ibraid <- inverseBraid(braid)
#' composeTwoBraids(braid, ibraid)
inverseBraid <- function(braid) {
  n <- numberOfStrands(braid)
  invgens <- lapply(braid, function(gen) {
    c(gen[1L], -gen[2L])
  })
  mkBraid(n, rev(invgens))
}

#' @title Composition of two braids
#' @description Composes two braids, doing free reduction on the result.
#' 
#' @param braid1,braid2 \code{braid} objects with the same number of strands
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' composeTwoBraids(braid, braid)
composeTwoBraids <- function(braid1, braid2) {
  n <- numberOfStrands(braid1) 
  if(n != numberOfStrands(braid2)) {
    stop("Unequal numbers of strands.")
  }
  freeReduceBraidWord(mkBraid(n, c(braid1, braid2)))
}

#' @title Composition of many braids.
#' @description Composes many braids, doing free reduction on the result.
#' 
#' @param braids list of \code{braid} objects with the same number of strands
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' composeManyBraids(list(braid, braid, braid))
composeManyBraids <- function(braids) {
  ns <- vapply(braids, numberOfStrands, integer(1L))
  n <- ns[1L]
  if(any(ns != n)) {
    stop("Unequal numbers of strands.")
  }
  freeReduceBraidWord(mkBraid(n, do.call(c, braids)))
}

#' @title Whether a braid is pure
#' @description Checks whether a braid is pure, i.e. its permutation is trivial.
#' 
#' @param braid a \code{braid} object
#'
#' @return A Boolean value.
#' @export
#'
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' isPureBraid(braid)
isPureBraid <- function(braid) {
  identical(braidPermutation(braid), seq_len(numberOfStrands(braid)))
}

#' @title Whether a braid is positive
#' @description Checks whether a braid has only positive Artin generators.
#' 
#' @param braid a \code{braid} object
#'
#' @return A Boolean value.
#' @export
#' 
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' isPositiveBraidWord(braid)
isPositiveBraidWord <- function(braid) {
  signs <- vapply(braid, `[[`, integer(1L), 2L)
  all(signs == 1L)
}

#' @title Linking matrix
#' @description Linking numbers between all pairs of strands of a braid.
#' 
#' @param braid a \code{braid} object
#'
#' @return A matrix.
#' @export
#' 
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' linkingMatrix(braid)
linkingMatrix <- function(braid) {
  n <- numberOfStrands(braid)
  perm <- 1L:n
  doSwap <- function(i) {
    a <- perm[i]
    b <- perm[i+1L]
    perm[i] <<- b
    perm[i+1L] <<- a
    invisible()
  }
  mat <- matrix(0L, nrow = n, ncol = n)
  doAdd <- function(i, j, pm1) {
    x <- mat[i, j]
    mat[i, j] <<- mat[j, i] <<- x + pm1
    invisible()
  }
  for(gen in braid) {
    k <- gen[1L]
    u <- perm[k]
    v <- perm[k+1L]
    doAdd(u, v, gen[2L])
    doSwap(k)
  }
  mat
}

#' @title Whether a braid is a permutation braid
#' @description Checks whether a braid is a permutation braid, that is, 
#'   a positive braid where any two strands cross at most one, and positively.
#' 
#' @param braid a \code{braid} object
#'
#' @return A Boolean value.
#' @export
#' 
#' @examples
#' braid <- mkBraid(4, list(Sigma(2), SigmaInv(3), Sigma(3)))
#' isPermutationBraid(braid)
isPermutationBraid <- function(braid) {
  if(isPositiveBraidWord(braid)) {
    lkMatrix <- linkingMatrix(braid)
    all(lkMatrix[upper.tri(lkMatrix)] %in% c(0L, 1L))
  } else {
    FALSE
  }
}

isPermutation <- function(x) {
  setequal(x, seq_along(x))
}

.permutationBraid <- function(perm) {
  n <- length(perm)
  cfwd <- cinv <- 1L:n
  doSwap <- function(i) {
    a <- cinv[i]
    b <- cinv[i+1L]
    cinv[i] <<- b
    cinv[i+1L] <<- a
    u <- cfwd[a]
    v <- cfwd[b]
    cfwd[a] <<- v
    cfwd[b] <<- u
    invisible()
  }
  worker <- function(phase) {
    if(phase >= n) {
      list()
    } else {
      tgt <- perm[phase]
      src <- cfwd[tgt]
      this <- (src-1L):phase
      lapply(this, doSwap)
      rest <- worker(phase + 1L)
      c(list(this), rest)
    }
  }
  worker(1L)
}

#' @title Permutation braid
#' @description Makes a permutation braid from a permutation.
#'
#' @param perm a permutation
#'
#' @return A \code{braid} object.
#' @export
#'
#' @examples
#' perm <- c(3, 1, 4, 2)
#' braid <- permutationBraid(perm)
#' isPermutationBraid(braid)
#' braidPermutation(braid)
permutationBraid <- function(perm) {
  stopifnot(isPermutation(perm), length(perm) >= 2)
  gens <- lapply(do.call(c, .permutationBraid(perm)), Sigma)
  mkBraid(length(perm), gens)
}

.allPositiveBraidWords <- function(n, l) {
  go <- function(k) {
    if(k == 0L) {
      list(list())
    } else {
      do.call(c, lapply(1L:(n-1L), function(i) {
        lapply(go(k - 1L), function(rest) {
          c(list(Sigma(i)), rest)
        })
      }))
    }
  }
  go(l)
}

#' @title Positive braid words of given length
#' @description All positive braid words of the given length.
#'
#' @param n number of strands, positive integer \code{>= 2}
#' @param l length of the words
#'
#' @return A list of \code{braid} objects.
#' @export
#'
#' @examples
#' allPositiveBraidWords(3, 4)
allPositiveBraidWords <- function(n, l) {
  stopifnot(isPositiveInteger(n), n >= 2)
  lapply(.allPositiveBraidWords(n, l), function(gens) {
    mkBraid(n, gens)
  })
}

.allBraidWords <- function(n, l) {
  go <- function(k) {
    if(k == 0L) {
      list(list())
    } else {
      gens <- do.call(c, lapply(1L:(n-1L), function(i) {
        c(list(Sigma(i)), list(SigmaInv(i)))
      }))
      do.call(c, lapply(go(k - 1L), function(rest) {
        lapply(gens, function(gen) {
          c(list(gen), rest)
        })
      }))
    }
  }
  go(l)
}

#' @title Braid words of given length
#' @description All braid words of the given length.
#'
#' @param n number of strands, positive integer \code{>= 2}
#' @param l length of the words
#'
#' @return A list of \code{braid} objects.
#' @export
#'
#' @examples
#' allPositiveBraidWords(3, 4)
allBraidWords <- function(n, l) {
  stopifnot(isPositiveInteger(n), n >= 2)
  lapply(.allBraidWords(n, l), function(gens) {
    mkBraid(n, gens)
  })
}