findNonTrans <- function()
{
  diceSet <- unique(t(apply(expand.grid(1:4, 1:4, 1:4, 1:4), 1, sort))) 
  winningDice <- function(X, Y) 
  {
    comparisonTable <- data.frame(X = rep(X, each = length(X)), Y = rep(Y, times = length(Y)))
    rowWinner <- ifelse(comparisonTable["X"] > comparisonTable["Y"], "X",
                        ifelse(comparisonTable["X"] == comparisonTable["Y"], "-", "Y"))
    netXWins <- sum(rowWinner == "X") - sum(rowWinner == "Y")
    if(netXWins > 0) "X" else if(netXWins == 0) "Draw" else "Y"
  }
  rows <- seq_len(nrow(diceSet)) 
  XvsAllY <- function(X) sapply(rows, function(i) winningDice(X, diceSet[i, ]))
  winners <- as.data.frame(lapply(rows, function(i) XvsAllY(diceSet[i, ])),
                row.names = paste("Y=Dice", rows), col.names = paste("X=Dice", rows), check.names = FALSE)
  solutionCount <- 0
  for(S in rows)
  {
    beatsS <- which(winners[paste("X=Dice", S)] == "Y") 
    beatsT <- lapply(beatsS, function(X) which(winners[paste("X=Dice", X)] == "Y")) 
    beatenByS <- which(winners[paste("X=Dice", S)] == "X") 
    potentialU <- lapply(beatsT, function(X) intersect(X, beatenByS)) 
    nonEmptyIndices <- lengths(potentialU) != 0 
    if(any(nonEmptyIndices)) 
    {
      solutionCount <- solutionCount + 1
      diceUIndex <- potentialU[nonEmptyIndices][[1]] 
      diceTIndex <- beatsS[nonEmptyIndices][[1]] 
      cat("Solution", solutionCount, "is:\n")
      output <- rbind(S = diceSet[S,], T = diceSet[diceTIndex,], U = diceSet[diceUIndex,])
      colnames(output) <- paste("Dice", 1:4)
      print(output)
    }
  }
}
findNonTrans()