CODEWITHSUNDEEP
rggplot2plotinequality
Adrian
Adrian

Reputation: 9793

R function to plot inequalities with shading

Suppose I have a set of inequalities:

  -2x + y <= -3
1.25x + y <= 2.5
        y >= -3

And I can summarize the information as follows:

mat <- matrix(c(-2, 1, 1.25, 1, 0, 1), nrow = 3, byrow = TRUE)
dir <- c("<=", "<=", ">=")
rhs <- c(-3, 2.5, -3)

I wrote the following function to plot the inequalities:

plot(0, 0, xlim = c(-1, 5), ylim = c(-4, 1))
plot_ineq <- function(mat, dir, rhs, xlow, xhigh){
  line <- list()
  for(i in 1:nrow(mat)){
    if(mat[i, 2] > 0){
      line[[i]] <- sapply(seq(xlow, xhigh, 0.1), function(x) (rhs[i] - mat[i, 1] * x)/mat[i, 2])
    }else if(mat[i, 2] < 0){
      line[[i]] <- sapply(seq(xlow, xhigh, 0.1), function(x) (rhs[i] - mat[i, 1] * x)/mat[i, 2])
      if(dir[i] == ">="){
        dir[i] = "<="
      }else dir[i] = ">="
    }
    lines(seq(xlow, xhigh, 0.1), line[[i]])
  }
}

plot_ineq(mat = mat, dir = dir, rhs = rhs, xlow = -1, xhigh = 5)

enter image description here

I have two questions: (1) how can I have a blank plot without having the (0, 0) point there? and (2) how can I shade the corresponding region according to dir? Should I try ggplot2?

I'm simply looking to shade the area that is described by the set of inequalities above. Not where (0, 0) lies.

Upvotes: 3

Views: 349

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 269461

1) Change the last inequality to be the same direction as the others and then use plotPolytope in gMOIP.

library(gMOIP)
mat <- matrix(c(-2, 1, 1.25, 1, 0, -1), nrow = 3, byrow = TRUE)
rhs <- c(-3, 2.5, 3)

argsFaces <- list(argsGeom_polygon = list(fill = "blue"))
plotPolytope(mat, rhs, argsFaces = argsFaces)

giving (continued after image)

screenshot

2) The above uses ggplot2 graphics but if you prefer classic graphics then using mat and rhs from above:

library(gMOIP)

cp <- cornerPoints(mat, rhs)
cp <- cp[chull(cp), ]  # chull gives indices of convex hull in order
plot(cp, type = "n")
polygon(cp, col = "blue")

# not shown but to add lines run this too
for(i in 1:nrow(cp)) {
  ix <- if (i < nrow(cp)) i + 0:1 else c(i, 1)
  b <- diff(cp[ix, 2]) / (d <- diff(cp[ix, 1]))
  if (abs(d) < 1e-5) abline(v = a <- cp[i, 1])
  else abline(a = a <- cp[i, 2] - b * cp[i, 1], b = b)
}

giving (continued after image)

screenshot

3) Note that there is an archived package named intpoint on CRAN and it could be used to draw the boundary of the feasible region and lines. It does have the limitation that it is hard coded to show X and Y axes between -1 and 5 although it might not be hard to generalize it. It is used like this (output not shown) where mat, rhs and cp are from above.

library(intpoint)

intpoint:::show2d(mat, rhs, c = numeric(2))
polygon(cp, col = "blue")

Upvotes: 3

Related Questions