Write vignette to be included with the package

[dirkschumacher]

Write a vignette that does not use a solver, otherwise we would have a circular dependency (or is that ok to depend on ompr.roi in Suggests?) Best candidate would the modelling vignette...

[CleKraus]

Hi Dirk,

first of all thanks a lot for this amazing package!

I'm currently playing around with the MILP backend. However, when I follow your vignette I get "Error: This function only works with models of type'optimization_model' from package ompr" when trying to solve the model. Reason seems to be that MILPModel() is of type "milp_model" rather than "optimization_model".

Do you have any suggestions on how to solve this?

Thanks a lot, Cle

[dirkschumacher]

Heyho @CleKraus ,

please also use the current github version of ompr.roi. I implemented this uneccessary check in the CRAN version 🙃. Will be fixed in the next version.

[CleKraus]

Hi Dirk,

cool, that works! Thanks a lot :-)

One other thing I came across though: When trying to slightly modify your MILP example by changing the constraint, it seems like the constraint was not picked up as e.g. sum over j for weights[1,j]*x[1,j] = 22 > 20 in my solution (see example below) Oddly, it works when I change to the MIP backend.

library(ROI) library(ROI.plugin.glpk) library(ompr.roi) library(dplyr)

set.seed(1) n <- 5 weights <- matrix(rpois(n * n, 5), ncol = n, nrow = n)

w <- function(i, j) { vapply(seq_along(i), function(k) weights[i[k], j[k]], numeric(1L)) }

result <- MILPModel() %>% add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>% set_objective(sum_expr(colwise(w(i, j)) x[i, j], i = 1:n, j = 1:n)) %>% add_constraint(sum_expr(colwise(w(i,j))x[i, j], j = 1:n) <= 20, i = 1:n) %>% solve_model(with_ROI("glpk", verbose = TRUE))

result <- get_solution(result, x[i, j]) %>% filter(value == 1) %>% rowwise() %>% mutate(weightsCol = weights[i,j]) %>% arrange(i,j) %>% select(i,j,weightsCol) print(result)

i j weightsCol 1 2 8 1 4 5 1 5 9 2 1 4 2 2 9 2 4 6 3 1 5 3 3 6 3 4 11 3 5 6 4 1 8 4 3 4 4 4 4 4 5 3 5 1 3 5 3 7 5 4 7 5 5 4

[dirkschumacher]

Can you please send me an email - the semantics of the milp backend are slightly different. This is unrelated to this issue. @CleKraus PS: really excited that you try the backend 🎉

[rickyars]

I think I'm getting the same error when using sum_expr in add_constraint with the new back-end. I'll try to make reproducible example. In the meantime, here's the error message:

Error in as.data.frame.default(x[[i]], optional = TRUE) : 
  cannot coerce class ‘structure("LinearVariableSum", package = "ompr")’ to a data.frame

[dirkschumacher]

@rickyars That woudl be great. Thanks! Then I can fix the issue before pushing the package to CRAN

[rickyars]

@CleKraus did you ever figure out how to use a matrix in the constraint? I'm having a hard time figuring out how to get this working on either the left-hand side or right-hand side.

[rickyars]

If you save the MIP and MILP models and extract the constraint matrices using extract_constraints, you'll see that there's a problem:

> mipc$matrix
5 x 25 sparse Matrix of class "dgCMatrix"

[1,] 4 . . . . 8 . . . . 3 . . . . 5 .  . . . 9 . . . .
[2,] . 4 . . . . 9 . . . . 3 . . . . 6  . . . . 3 . . .
[3,] . . 5 . . . . 6 . . . . 6 . . . . 11 . . . . 6 . .
[4,] . . . 8 . . . . 6 . . . . 4 . . .  . 4 . . . . 3 .
[5,] . . . . 3 . . . . 2 . . . . 7 . .  . . 7 . . . . 4

> milpc$matrix
5 x 25 sparse Matrix of class "dgCMatrix"

[1,] 4 . . . . 9 . . . . 6 . . . . 4 . . . . 4 . . . .
[2,] . 4 . . . . 9 . . . . 6 . . . . 4 . . . . 4 . . .
[3,] . . 4 . . . . 9 . . . . 6 . . . . 4 . . . . 4 . .
[4,] . . . 4 . . . . 9 . . . . 6 . . . . 4 . . . . 4 .
[5,] . . . . 4 . . . . 9 . . . . 6 . . . . 4 . . . . 4

[dirkschumacher]

@rickyars is this the example from above?

[rickyars]

@dirkschumacher yes, this is using the example shared by @CleKraus above.

I'm trying to wrap my head around this colwise/vectorized change for MILP and I just can't get anything to work. So I thought I was start with a simple example.

[dirkschumacher]

Does this help? https://dirkschumacher.github.io/ompr/dev/articles/milp-modelling.html#vectorized-semantics-revisited

[rickyars]

@dirkschumacher i've read through all the issues where colwise appears and all the documentation, but i just can't get this to work. i've tried the seq_length and matrix functions described by @hugolarzabal. i've even tried to use array_branch from the purrr package. for whatever reason, it just doesn't work. have you tried the example above?

[hugolarzabal]

@rickyars the w function in the exemple above doesn't seems to work as needed by the milp model.

w(1,1:2) should return c(4, 8) but it returns 4

Try this for exemple (I think it gives the right result):

w <- list(weights[1,],weights[2,],weights[3,],weights[4,],weights[5,])
  result <- MILPModel() %>%
    add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>%
    set_objective(sum_expr(colwise(w[i]) * x[i, j], i = 1:n, j = 1:n)) %>%
    add_constraint(sum_expr(colwise(w[i])*x[i, j], j = 1:n) <= 20, i = 1:n) %>%
  solve_model(with_ROI("glpk", verbose = TRUE))

What colwise does is make the sum_expr function understand that it should take the next element of w for every i and match the values of w[i] with every j.

[rickyars]

@hugolarzabal I can't even get that to run.

[hugolarzabal]

I looked at that problem a little further today and there is indeed something that doesn t seems to work when using sum_expr with colwise.

I created a simple exemple using the sum_fun function from

library(ompr)
library(dplyr)
library(magrittr)

sum_fun <- function(i, j) {
  n_rows <- length(i)
  n_col <- length(j)
  colwise(vapply(seq_len(n_rows * n_col), function(x) {
    x # 1, 2, 3, 4, 5, 6
    # where 1, 2, 3 are for cols in the first row
    # 4, 5, 6 are the cols in the second row and so on
  }, numeric(1L)))
}

n <- 5
result <- MILPModel()
  result %<>% add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") 
  result %<>% set_objective(sum_expr(sum_fun(i, j) * x[i, j], i = 1:n, j = 1:n)) 
  result %<>% add_constraint(sum_expr( sum_fun(i,j) * x[i, j], j = 1:n) <= 20, i = 1:n) 

I always get the following error message: Error in[[<-.data.frame(tmp, "coef", value = c(1, 2, 3, 4, 5, 6, : replacement has 625 rows, data has 25

First the messaged appeared when executing the last line
result %<>% add_constraint(sum_expr( sum_fun(i,j) * x[i, j], j = 1:n) <= 20, i = 1:n)

Then, I installed the latest version of ompr and the last line worked fine but the error appeared in another line:
result %<>% set_objective(sum_expr(sum_fun(i, j) * x[i, j], i = 1:n, j = 1:n))

[hugolarzabal]

@dirkschumacher : I do think that the problem is caused by the use of expand.grid in sum_expr.

[dirkschumacher]

Thanks @hugolarzabal . Will try to look into it. I am currenty super busy and ompr is sadly just a hobby project that I do after work.

[dirkschumacher]

Regarding your example:

n <- 5
result <- MILPModel() %>% 
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>% 
  set_objective(sum_expr(colwise(1:25) * x[i, j], i = 1:n, j = 1:n)) %>% 
  add_constraint(sum_expr(colwise(1:25) * x[i, j], j = 1:5) <= 20, i = 1:5) 
ompr::extract_constraints(result)

This is what your summation function needs to produce to have the desired effect.

I really need to think harder to make this more userfriendly :)

[rickyars]

I found something that works!!

Try purrr::array_branch():

w_i <- purrr::array_branch(weights) %>% as.numeric()
w_j <- purrr::array_branch(t(weights)) %>% as.numeric()

milp <- MILPModel() %>%
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>%
  set_objective(sum_expr(colwise(w_i) * x[i, j], i = 1:n, j = 1:n)) %>%
  add_constraint(sum_expr(colwise(w_j) * x[i, j], j = 1:n) <= 20, i = 1:n) 

However, I don't understand why I have to use two versions of the the weight matrix. I guess it has to do with the indices used in the sum_expr call. I'm also not sure how I'd replicate this for multi-dimensional matrices. I have a some weight matrices with 3 indices.

[dirkschumacher]

sum_expr(x[i, j], i = 1:n, j = 1:n)

Defines variables in one row with n*n many columns. Thus you need a colwise vector of length 25, where each component is one coefficient on that row.

sum_expr(x[i, j], j = 1:n) <= 20, i = 1:n) 

Defines variables over n rows with n columns each. Thus again you need a colwise vector of length 25, but ordered differently. The first n components go in the first row, the second n components in the second one and so forth.

[rickyars]

I just can't figure this out. Here's a piece of my code that works with MIP:

model <- model %>%
  # satisfy all passenger demand
  add_constraint(
    sum_expr(capacity[k] * mm[i, k, c], k = 1:K) >=
      sum_expr(pax.dem[i, c] * sm[i, k, c], k = 1:K),
    i = 1:I, c = 1:C)

How do I format pax.dem[i,c] for this constraint using MILP?

I realize MILP offers big performance gains over MIP but at least MIP works intuitively.

[dirkdegel]

I think there is a small error in the documentation: https://dirkschumacher.github.io/ompr/dev/articles/milp-modelling.html#vectorized-semantics-revisited

In this example:

n <- 10L
m <- 20L

model <- MILPModel() %>% 
  add_variable(x[i, j], i = 1:m, j = 1:n) %>% 
  #add_constraint(sum_expr(x[i, j], i = 1:m) == 1, j = 1:n) #%>% 
  #add_constraint(x[1:m, colwise(1:n)] == 1) # this this equivalent # -> wrong
  # should be this
  add_constraint(x[colwise(1:m), 1:n] == 1) #%>% 

extract_constraints(model)

In the example in the vignette only n <- 10 is used.

I still find it very difficult and think it should be possible to use the MILP backend in the same way as the MIP backend.

I am also still not really able to implement a set cover model and find the results very confusing.

[dirkschumacher]

Yeah, the MILP backend is a bit unintuitive. I will try to improve it once I find some time.

@dirkdegel thanks. Could you please make a seperate issue with the bug you found?

[dirkschumacher]

Can we find another way to communicate? This issue as evolved into a discuassion board :)

[rickyars]

I agree this has become a bit of a combination discussion board and Stackover flow for the package. What's the best way to communicate. I'm also happy to help with some of the coding, but looking at this part of the code I realize I won't be much help until I get more experience with how this works and package development.

[hugolarzabal]

I am also ready to help.

I think that the ultimate solution to all those problems with sum_expr and colwise would be to have a standardized object in ompr to represent parameters.

In my mind, it could be something that works like that: add_parameter(p[i,j], i=1:n, j=1:m, values = 1:(n*m))

Then p could be used just like a variable in every constraint.

Obviously, I realise that it would require quite some work. I am just giving that as an idea for future improvement.

[dirkdegel]

I am happy to help as well, but have very limited R skills. I think the best would be to have the option to pass parameters in exactly the same way as in the mathematical notation, e.g. set_objective(sum_expr(weight[i, j] * x[i, j], i = 1:n, j = 1:m)) I am not sure if that is doable, but it would be the best (most natural) syntax.