Simplify this grid such that each row and column has 1 value

Question

Example code here:

> temp2
  a b c d e f g h
i 1 1 0 0 0 1 0 1
j 0 1 0 0 0 1 0 1
k 0 1 1 0 0 1 1 1
l 0 0 0 0 1 0 0 1
m 0 0 1 1 0 0 1 1
n 0 0 1 1 0 0 1 1
o 0 0 0 1 0 0 1 1
p 0 0 0 0 1 0 0 1

> dput(temp2)
structure(list(a = c(1, 0, 0, 0, 0, 0, 0, 0), b = c(1, 1, 1, 
0, 0, 0, 0, 0), c = c(0, 0, 1, 0, 1, 1, 0, 0), d = c(0, 0, 0, 
0, 1, 1, 1, 0), e = c(0, 0, 0, 1, 0, 0, 0, 1), f = c(1, 1, 1, 
0, 0, 0, 0, 0), g = c(0, 0, 1, 0, 1, 1, 1, 0), h = c(1, 1, 1, 
1, 1, 1, 1, 1)), .Names = c("a", "b", "c", "d", "e", "f", "g", 
"h"), class = "data.frame", row.names = c("i", "j", "k", "l", 
"m", "n", "o", "p"))

I have this 8x8 grid of 1s and 0s. I need to solve for some grid where each row and each column has exactly one 1 and the rest 0s, but the 1 has to be in a place where the original grid has a 1. It's almost like a sudoku question but not exactly. Any thoughts on how to get started?

I would need some function that can do this for a general grid, not simply this specific one. We can assume that there's always a solution grid, given some starting grid.

Thanks!

Edit: a valid solution

> temp3
  a b c d e f g h
i 1 0 0 0 0 0 0 0
j 0 1 0 0 0 0 0 0
k 0 0 0 0 0 1 0 0
l 0 0 0 0 1 0 0 0
m 0 0 0 1 0 0 0 0
n 0 0 1 0 0 0 0 0
o 0 0 0 0 0 0 1 0
p 0 0 0 0 0 0 0 1

EDIT2: given that there's only 8! unique solutions for any grid, i may attempt a brute force / matching approach.

Answer 1

This can be solved as a transportation problem or as an integer programming problem. We also show a one-line solution using only base R which generates random matrices for which each row and each columns column sums to 1 filtering out and returning the ones satisfying the additional constraints that each element of the solution matrix be less than or equal to the corresponding element of temp2.

1) transportation problem Using lp.transport in lpSolve we can solve it in one statement:

library(lpSolve)

res <- lp.transport(as.matrix(temp2), "max", 
  rep("=", 8), rep(1, 8), rep("=", 8), rep(1, 8), integers = 0:1)

res
## Success: the objective function is 8

soln <- array(res$solution, dim(temp2))

# verify

all(colSums(soln)==1) && all(rowSums(soln)==1) && all(temp2>=soln) && all(soln %in% 0:1)
## [1] TRUE

2) integer programming

If X is the solution we have specified the row and column constraints but have not specified the X <= temp2 constraints since they will be satisfied automatically as no solution putting a 1 where a temp2 0 is can have the maximum objective of 8.

library(lpSolve)

n <- nrow(temp2)
obj <- unlist(temp2)
const_row <- t(sapply(1:n, function(i) c(row(temp2)) == i)) # each row sums to 1
const_col <- t(sapply(1:n, function(i) c(col(temp2)) == i)) # each col sums to 1
const.mat <- rbind(const_row, const_col)
res <- lp("max", obj, const.mat, "=", 1, all.bin = TRUE)
res
## Success: the objective function is 8

soln <- array(res$solution, dim(temp2))

# verify

all(colSums(soln)==1) && all(rowSums(soln)==1) && all(temp2>=soln) && all(soln %in% 0:1)
## [1] TRUE

(Note that by the same argument we could have relaxed the problem to a linear programming problem provided we add 0 <= soln[i, j] <= 1 constraints since by the same argument that allowed us to omit the soln[i, j] <= temp2[i, j] constraints the maximization will force the soln elements to be 0 or 1 anyways.)

2a) This approach is longer but does spell out the X <= temp2 constraints explicitly:

n <- nrow(temp2)
obj <- numeric(n*n)
const1 <- diag(n*n) # soln[i,j] <= temp2[i,j]
const2 <- t(sapply(1:n, function(i) c(row(temp2)) == i)) # each row sums to 1
const3 <- t(sapply(1:n, function(i) c(col(temp2)) == i)) # each col sums to 1
const.mat <- rbind(const1, const2, const3)
const.dir <- rep(c("<=", "="), c(n*n, 2*n))
const.rhs <- c(unlist(temp2), rep(1, 2*n))

res <- lp("max", obj, const.mat, const.dir, const.rhs, all.bin = TRUE)
res
## Success: the objective function is 0

soln <- array(res$solution, dim(temp2))

# verify

all(colSums(soln)==1) && all(rowSums(soln)==1) && all(temp2>=soln) && all(soln %in% 0:1)
## [1] TRUE

2b) Note that if X is the solution matrix then in X <= temp2 only the positions of X corresponding to zeros in temp2 actually constrain so we could eliminate any constraint corresponding to a 1 in temp2 in the (2a) solution. With this change all constraints become equality constraints.

n <- nrow(temp2)
obj <- numeric(n*n)
const1 <- diag(n*n)[unlist(temp2) == 0, ]
const2 <- t(sapply(1:n, function(i) c(row(temp2)) == i)) # each row sums to 1
const3 <- t(sapply(1:n, function(i) c(col(temp2)) == i)) # each col sums to 1
const.mat <- rbind(const1, const2, const3)
const.dir <- "="
const.rhs <- c(numeric(nrow(const1)), rep(1, 2*n))

res <- lp("max", obj, const.mat, const.dir, const.rhs, all.bin = TRUE)
res
## Success: the objective function is 0

soln <- array(res$solution, dim(temp2))

# verify

all(colSums(soln)==1) && all(rowSums(soln)==1) && all(temp2>=soln) && all(soln %in% 0:1)
## [1] TRUE

In fact, we could go further and remove the variables that correspond to zero elements of temp2.

3) r2dtable Here we use rd2table to generate 10,000 8x8 tables whose rows and columns sum to 1 and then filter them to pick out only those satisfying the X < temp2 constrainsts. Withtemp2` from the question and the random seed shown has found 3 solutions. If with different inputs it finds no solutions then try generating a higher number of random proposals. This approach does not use any packages.

set.seed(123) # for reproducibility
Filter(function(x) all(x <= temp2), r2dtable(10000, rep(1, 8), rep(1, 8)))

giving:

[[1]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    0    0    0    0    0    0    0
[2,]    0    0    0    0    0    1    0    0
[3,]    0    1    0    0    0    0    0    0
[4,]    0    0    0    0    0    0    0    1
[5,]    0    0    0    0    0    0    1    0
[6,]    0    0    1    0    0    0    0    0
[7,]    0    0    0    1    0    0    0    0
[8,]    0    0    0    0    1    0    0    0

[[2]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    0    0    0    0    0    0    0
[2,]    0    0    0    0    0    1    0    0
[3,]    0    1    0    0    0    0    0    0
[4,]    0    0    0    0    1    0    0    0
[5,]    0    0    0    1    0    0    0    0
[6,]    0    0    1    0    0    0    0    0
[7,]    0    0    0    0    0    0    1    0
[8,]    0    0    0    0    0    0    0    1

[[3]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    0    0    0    0    0    0    0
[2,]    0    1    0    0    0    0    0    0
[3,]    0    0    0    0    0    1    0    0
[4,]    0    0    0    0    1    0    0    0
[5,]    0    0    1    0    0    0    0    0
[6,]    0    0    0    0    0    0    1    0
[7,]    0    0    0    1    0    0    0    0
[8,]    0    0    0    0    0    0    0    1

Answer 2

A brute-force way:

m = as.matrix(temp2)
w = data.frame(which(m == 1, arr.ind = TRUE))
combos = as.matrix(do.call(expand.grid, with(w, split(col, row))))
combos[ apply(combos, 1, function(x) !anyDuplicated(x)), ]


      1 2 3 4 5 6 7 8
 [1,] 1 6 2 8 7 3 4 5
 [2,] 1 2 6 8 7 3 4 5
 [3,] 1 6 2 8 3 7 4 5
 [4,] 1 2 6 8 3 7 4 5
 [5,] 1 6 2 8 4 3 7 5
 [6,] 1 2 6 8 4 3 7 5
 [7,] 1 6 2 8 3 4 7 5
 [8,] 1 2 6 8 3 4 7 5
 [9,] 1 6 2 5 7 3 4 8
[10,] 1 2 6 5 7 3 4 8
[11,] 1 6 2 5 3 7 4 8
[12,] 1 2 6 5 3 7 4 8
[13,] 1 6 2 5 4 3 7 8
[14,] 1 2 6 5 4 3 7 8
[15,] 1 6 2 5 3 4 7 8
[16,] 1 2 6 5 3 4 7 8

OP claims to only ever need to handle an 8x8 grid, so I guess this performs well enough. Each row of the result is a solution. The first row says that (1,1), (2,6), (3,2) ... is a solution.

---

A variation using data.table:

library(data.table)
m = as.matrix(temp2)
comboDT = setDT(melt(m))[ value == 1, do.call(CJ, split(Var2, Var1)) ][, 
  rid := .I ][, melt(.SD, id="rid", variable.name="row", value.name="col")]

setkey(comboDT, rid)
comboDT[ .( comboDT[, !anyDuplicated(col), by=rid][(V1), rid]) ]

Answer 3

this works. Let grid be my grid (temp2 from above). then this will return a grid that works

# create random sufficient grid
  counter = 0
  while(2 > 1) {
    counter = counter + 1
    if(counter == 10000) {
      break
    }
    rand_grid = matrix(0, nrow = 8, ncol = 8)
    indices_avail = seq(1,8,by=1)
    for(i in 1:8) {
      k = sample(indices_avail, 1)
      rand_grid[i, k] = 1
      indices_avail = indices_avail[indices_avail != k]
    }
    if(sum(grid[which(rand_grid == 1)]) == 8) {
      break
    }
    print(counter)
  }

Answer 4

This approach will return all valid combinations. First find all matrix row combinations. Then search through exhaustively. This method would have to be improved if your matrix size increased. One simple improvement would be to run the diag test in parallel.

st<-as.matrix(temp2) # make sure we are working with matrices
## This method will return all possible matrices of combinations
## in essence if you have diag(matr) = width matrix than you have
## a valid choice


## Helper function to build all combinations, there may be better way to
## do this but it gets the job done
allCombinationsAux<-function(z,nreg,x){
    if(sum(nreg)>1){
        innerLoop<-do.call(rbind,lapply(x[nreg&(z!=x)], test1,nreg&(z!=x),x))
        ret<-cbind(z,innerLoop )
    }
    else{
        ret<-x[nreg]
    }
    ret
}

## Build all of the combinations of possible matrices
combs<-do.call(rbind,lapply(x,function(y) allCombinationsAux(y,y!=x,x)))

## iterate through all the possible combinations of matrices, to find out
## which ones have 1s throughout the  diag
inds<-which(apply(combs,1,function(x) sum(diag(st[x,]))==8))

lapply(inds,function(x) st[combs[x,],])