---
title: "Introduction to gMOIP"
author: "Lars Relund Nielsen"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to gMOIP}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r setup, include = FALSE}
library(knitr)
library(rgl)
library(ggsci)
library(magrittr)
rgl::setupKnitr()
options(rgl.useNULL=TRUE)
rgl::par3d("family" = "serif")
opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning=FALSE, message=FALSE, include = TRUE, 
  out.width = "99%", fig.width = 8, fig.align = "center", fig.asp = 0.62
)

if (!requireNamespace("rmarkdown", quietly = TRUE) || !rmarkdown::pandoc_available("1.14")) {
   warning(call. = FALSE, "These vignettes assume rmarkdown and pandoc version 1.14 (or higher). These were not found. Older versions will not work.")
   knitr::knit_exit()
}
```

The `gMOIP` package can be used to make 2D and 3D plots of linear programming (LP), integer linear
programming (ILP), or mixed integer linear programming (MILP) models with up to three objectives.
This include the polytope, integer points, ranges and iso profit curve. Plots of both the solution
and criterion space are possible. For instance the non-dominated (Pareto) set for bi-objective
LP/ILP/MILP programming models.

The package also include an `inHull` function for checking if a set of points is inside/at/outside
the convex hull of a set of vertices (for arbitrary dimension).

Finally, the package also contains functions for generating (non-dominated) points in $\mathbb{R}_n$
and classifying non-dominated points as supported extreme, supported non-extreme and unsupported.


## Installation

Install the latest stable release from CRAN:

``` r
install.packages("gMOIP")
```

Alternatively, install the latest development version from GitHub (recommended):

``` r
install.packages("devtools")
devtools::install_github("relund/gMOIP")
library(gMOIP)
```

```{r, cache=FALSE, include=FALSE}
library(gMOIP)
```


## Single criterion models

We define the model \(\max\{cx | Ax \leq b\}\) (could also be minimized) using matrix `A` and
vectors `b` and `obj`:

```{r}
A <- matrix(c(-3,2,2,4,9,10), ncol = 2, byrow = TRUE)
b <- c(3,27,90)
obj <- c(7.75, 10)  # coefficients c
```

The polytope of the LP model with non-negative continuous variables
(\(x \geq 0\)):

```{r}
plotPolytope(
   A,
   b,
   obj,
   type = rep("c", ncol(A)),
   crit = "max",
   faces = rep("c", ncol(A)),
   plotFaces = TRUE,
   plotFeasible = TRUE,
   plotOptimum = TRUE,
   labels = "coord"
)
```

The polytope of the ILP model with LP faces (\(x\in \mathbb{Z}_0\)):

```{r}
plotPolytope(
   A,
   b,
   obj,
   type = rep("i", ncol(A)),
   crit = "max",
   faces = rep("c", ncol(A)),
   plotFaces = TRUE,
   plotFeasible = TRUE,
   plotOptimum = TRUE,
   labels = "coord"
)
```

The polytope of the MILP model (first variable integer) with LP faces:

```{r}
plotPolytope(
   A,
   b,
   obj,
   type = c("i", "c"),
   crit = "max",
   faces = c("c", "c"),
   plotFaces = TRUE,
   plotFeasible = TRUE,
   plotOptimum = TRUE,
   labels = "coord"
)
```

### Three variables

You can do the same with three variables:

```{r}
A <- matrix( c(
   3, 2, 5,
   2, 1, 1,
   1, 1, 3,
   5, 2, 4
), nc = 3, byrow = TRUE)
b <- c(55, 26, 30, 57)
obj <- c(20, 10, 15)
```

LP model:

```{r, webgl=TRUE}
view <- matrix( c(-0.412063330411911, -0.228006735444069, 0.882166087627411, 0, 0.910147845745087,
                  -0.0574885793030262, 0.410274744033813, 0, -0.042830865830183, 0.97196090221405,
                  0.231208890676498, 0, 0, 0, 0, 1), nc = 4)   
loadView(v = view)  # set view angle
plotPolytope(A,
             b,
             obj,
             plotOptimum = TRUE,
             labels = "n")
```

ILP model (here with ILP faces):

```{r, webgl=TRUE}
view <- matrix( c(-0.412063330411911, -0.228006735444069, 0.882166087627411, 0, 0.910147845745087,
                  -0.0574885793030262, 0.410274744033813, 0, -0.042830865830183, 0.97196090221405,
                  0.231208890676498, 0, 0, 0, 0, 1), nc = 4)   
loadView(v = view)  # set view angle
plotPolytope(A,
             b,
             obj, 
             type = c("i", "i", "i"), 
             plotOptimum = TRUE, 
             labels = "n")
```

MILP model (here with continuous faces):

```{r, webgl=TRUE, cache=FALSE}
view <- matrix( c(-0.412063330411911, -0.228006735444069, 0.882166087627411, 0, 0.910147845745087,
                  -0.0574885793030262, 0.410274744033813, 0, -0.042830865830183, 0.97196090221405,
                  0.231208890676498, 0, 0, 0, 0, 1), nc = 4)   
loadView(v = view)  # set view angle
plotPolytope(A,
             b,
             obj, 
             type = c("i", "i", "c"), 
             faces = c("c", "c", "c"), 
             plotOptimum = TRUE, 
             # plotFaces = FALSE,
             labels = "n")
```


## Bi-objective models

With `gMOIP` you can also make plots of the criterion space for bi-objective models (linear
programming (LP), integer linear programming (ILP), or mixed integer linear programming (MILP)).

First let us have a look at bi-objective model with two variables. We define a function for
grouping plots of the solution and criterion space:

```{r 2DPlotFunction}
plotBiObj2D <- function(A, b, obj,
   type = rep("c", ncol(A)),
   crit = "max",
   faces = rep("c", ncol(A)),
   plotFaces = TRUE,
   plotFeasible = TRUE,
   plotOptimum = FALSE,
   labels = "numb",
   addTriangles = TRUE,
   addHull = TRUE)
{
   p1 <- plotPolytope(A, b, type = type, crit = crit, faces = faces, plotFaces = plotFaces,
                      plotFeasible = plotFeasible, plotOptimum = plotOptimum, labels = labels) + 
      ggplot2::ggtitle("Solution space")
   p2 <- plotCriterion2D(A, b, obj, type = type, crit = crit, addTriangles = addTriangles,
                         addHull = addHull, plotFeasible = plotFeasible, labels = labels) +
      ggplot2::ggtitle("Criterion space")
   gridExtra::grid.arrange(p1, p2, nrow = 1) 
}
```

Let us define the constraints:

```{r 2DConstraints}
A <- matrix(c(-3,2,2,4,9,10), ncol = 2, byrow = TRUE)
b <- c(3,27,90)
```

First let us have a look at a LP model (maximize):

```{r 2DLP}
obj <- matrix(
   c(7, -10, # first criterion
     -10, -10), # second criterion
   nrow = 2)
plotBiObj2D(A, b, obj, addTriangles = FALSE)
```

Note the non-dominated (Pareto) set consists of all supported extreme non-dominated points
(illustrated with triangles) and the line segments between them.

ILP model (maximize):

```{r 2DILPMax}
obj <- matrix(c(3, -1, -2, 2), nrow = 2)
plotBiObj2D(A, b, obj, type = rep("i", ncol(A)))
```

Note the non-dominated set consists of all points in black (with shape supported extreme:triangle,
supported non-extreme:round, unsupported:round (not on the border)). The triangles drawn using the
extreme non-dominated points illustrate areas where unsupported non-dominated points may be found.
A point in the solution space is identified in the criterion space using the same number.

MILP model ($x_2$ integer) (minimize):

```{r 2DMILPMin}
obj <- matrix(c(3, -1, -2, 2), nrow = 2)
plotBiObj2D(A, b, obj, type = c("c", "i"), crit = "min")
```

Note the solution space now consists to segments and hence the non-dominated set may consist of
points and segments (open and closed). Note these segments are not highlighted in the current 
version of `gMOIP`.


## Bi-objective models with three variables

We define functions for plotting the solution and criterion space:

```{r 3DPlotFunction}
plotSol <- function(A, b, type = rep("c", ncol(A)),
                        faces = rep("c", ncol(A)),
                        plotFaces = TRUE, labels = "numb")
{
   loadView(v = view, close = F, zoom = 0.75)
   plotPolytope(A, b, type = type, faces = faces, labels = labels, plotFaces = plotFaces, 
                argsTitle3d = list(main = "Solution space"))
}

plotCrit <- function(A, b, obj, crit = "min", type = rep("c", ncol(A)), addTriangles = TRUE, 
                     labels = "numb") 
{
    plotCriterion2D(A, b, obj, type = type, crit = crit, addTriangles = addTriangles, 
                   labels = labels) + 
      ggplot2::ggtitle("Criterion space")
}
```

Four examples are given. A few plots of the solution space are made interactive to illustrate the
functionality.

### Example 1

We define the model  $\max \{cx | Ax \leq b\}$ (could also be minimized) with three variables:

```{r Ex1Model}
Ab <- matrix( c(
   1, 1, 2, 5,
   2, -1, 0, 3,
   -1, 2, 1, 3,
   0, -3, 5, 2
), nc = 4, byrow = TRUE)
A <- Ab[,1:3]
b <- Ab[,4]
obj <- matrix(c(1, -6, 3, -4, 1, 6), nrow = 2)
```

We load the preferred view angle for the RGL window:

```{r Ex1View}
view <- matrix( c(-0.452365815639496, -0.446501553058624, 0.77201122045517, 0, 0.886364221572876,
                  -0.320795893669128, 0.333835482597351, 0, 0.0986008867621422, 0.835299551486969,
                  0.540881276130676, 0, 0, 0, 0, 1), nc = 4)
loadView(v = view)
```

LP model (solution space):

```{r Ex1LPSol, webgl=TRUE, cache=FALSE}
plotSol(A, b)
```

LP model (criterion space):

```{r Ex1LPCrit}
plotCrit(A, b, obj, addTriangles = FALSE) 
```

ILP model (solution space):

```{r Ex1ILPSol, webgl=TRUE, cache=FALSE}
plotSol(A, b, type = c("i","i","i"))
```

ILP model (criterion space):

```{r Ex1ILPCrit}
plotCrit(A, b, obj, type = c("i","i","i"))
```

## LaTeX support

You may create a TikZ file of the 2D plots for LaTeX using

``` r
library(tikzDevice)
tikz(file = "plot_polytope.tex", standAlone=F, width = 7, height = 6)
plotPolytope(
   A,
   b,
   obj,
   type = rep("i", ncol(A)),
   crit = "max",
   faces = rep("c", ncol(A)),
   plotFaces = TRUE,
   plotFeasible = TRUE,
   plotOptimum = TRUE,
   labels = "coord"
)
dev.off()
```

## Further examples

For further examples see the [documentation](https://relund.github.io/gMOIP/), the vignettes and the [articles](https://relund.github.io/gMOIP/).

``` r
library(gMOIP)
browseVignettes('gMOIP')
example("gMOIP-package")
```

```{r, include=F}
rm(list = ls(all.names = TRUE))
```

