## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(diffcp)
library(Matrix)

## ----lp-data------------------------------------------------------------------
A <- sparseMatrix(i = c(1, 2, 1, 3, 1, 4),
                  j = c(1, 1, 2, 2, 3, 3),
                  x = c(1, -1, 1, -1, 1, -1),
                  dims = c(4, 3))
b <- c(1, 0, 0, 0)
c <- c(1, 2, 3)
cone_dict <- list(z = 1L, l = 3L)

res <- solve_and_derivative(A, b, c, cone_dict, mode = "lsqr",
                            tol_gap_abs = 1e-10,
                            tol_gap_rel = 1e-10,
                            tol_feas    = 1e-10)
res$x

## ----lp-D---------------------------------------------------------------------
dA <- sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                   dims = c(4, 3))
db <- numeric(4)
dc <- c(0.001, 0, 0)
d <- res$D(dA, db, dc)
d$dx

## ----lp-fd--------------------------------------------------------------------
res_pert <- solve_only(A, b, c + dc, cone_dict,
                       tol_gap_abs = 1e-10,
                       tol_gap_rel = 1e-10,
                       tol_feas    = 1e-10)
max(abs(d$dx - (res_pert$x - res$x)))

## ----soc-data-----------------------------------------------------------------
A <- sparseMatrix(
  i = c(2, 1, 3, 1, 4, 1, 5),
  j = c(1, 2, 2, 3, 3, 4, 4),
  x = c(-1, 1, -1, 1, -1, 1, -1),
  dims = c(5, 4))
b <- c(1, 0, 0, 0, 0)
c <- c(1, 0, 0, 0)
cone_dict <- list(z = 1L, q = 4L)

res <- solve_and_derivative(A, b, c, cone_dict, mode = "dense",
                            tol_gap_abs = 1e-10,
                            tol_gap_rel = 1e-10,
                            tol_feas    = 1e-10)
res$x

## ----sdp-data-----------------------------------------------------------------
DIM <- 5L
n   <- DIM * (DIM + 1L) %/% 2L          # 15

trace_row <- numeric(n)
pos <- 1L
for (col in seq_len(DIM)) {
  trace_row[pos] <- 1.0
  pos <- pos + (DIM - col + 1L)         # next diagonal in column-major lower-tri
}

A <- rbind(
  Matrix(matrix(trace_row, nrow = 1L), sparse = TRUE),
  -Diagonal(n))
A <- as(A, "CsparseMatrix")
b <- c(1.0, numeric(n))

C <- matrix(0, DIM, DIM)
C[1, 1] <- 1; C[DIM, DIM] <- -1
c_vec <- vec_symm(C)

cone_dict <- list(z = 1L, s = DIM)

res <- solve_only(A, b, c_vec, cone_dict,
                  tol_gap_abs = 1e-10,
                  tol_gap_rel = 1e-10,
                  tol_feas    = 1e-10)
X <- unvec_symm(res$x, DIM)
sum(diag(X))
range(eigen(X, symmetric = TRUE, only.values = TRUE)$values)

## ----shadow-data--------------------------------------------------------------
A_lp <- sparseMatrix(i = c(1, 2, 3,  1, 2, 4),
                     j = c(1, 1, 1,  2, 2, 2),
                     x = c(2, 1, -1, 1, 3, -1),
                     dims = c(4, 2))
b_lp <- c(100, 90, 0, 0)
c_lp <- c(-3, -5)
cones_lp <- list(l = 4L)

res <- solve_and_derivative(A_lp, b_lp, c_lp, cones_lp, mode = "dense",
                            tol_gap_abs = 1e-10, tol_gap_rel = 1e-10,
                            tol_feas    = 1e-10)
res$x      # optimal production plan
-sum(c_lp * res$x)   # max profit

## ----shadow-D-----------------------------------------------------------------
dA <- sparseMatrix(i = integer(0), j = integer(0), x = numeric(0),
                   dims = c(4, 2))
db <- c(1, 0, 0, 0)
dc <- numeric(2)
d  <- res$D(dA, db, dc)
d$dx                       # change in optimal production
-sum(c_lp * d$dx)          # change in profit

## ----shadow-DT----------------------------------------------------------------
## Suppose dL/dx* is given (e.g. from autodiff in a larger model).
dL_dx <- c(1.0, -0.5)
dL_dy <- numeric(length(res$y))
dL_ds <- numeric(length(res$s))
da    <- res$DT(dL_dx, dL_dy, dL_ds)
da$db                 # gradient of L w.r.t. each entry of b
da$dc                 # gradient of L w.r.t. each entry of c