#########################
# Quadratic Programming #
#########################

# Goal: min -d^tx+0.5x^tDx
# Subject to: A^tx>=b
install.packages("quadprog")
library(quadprog)

# Example 1: unconstrained
Dmat = matrix(c(2,0,0,8),2,2,byrow=T)
dvec = c(8,16)
Amat = matrix(c(0,0,0,0),2,2,byrow=T)
bvec = c(0,0)
solve.QP(Dmat,dvec,Amat,bvec)

# Example 1b: unconstrained, no minimum because D is not positive definite
Dmat = matrix(c(2,4,4,8),2,2,byrow=T)
solve.QP(Dmat,dvec,Amat,bvec)

# Example 1c: unconstrained, no minimum because D is not positive definite
Dmat = matrix(c(-2,0,0,8),2,2,byrow=T)
solve.QP(Dmat,dvec,Amat,bvec)


# Example 2: non-binding constraints 
Dmat = matrix(c(2,0,0,8),2,2,byrow=T)
dvec = c(8,16)
Amat = matrix(c(1,1,1,0),2,2,byrow=T)
bvec = c(5,3)
solve.QP(Dmat,dvec,Amat,bvec)

Amat = t(matrix(c(1,1,1,0,1,0,0,1),4,2,byrow=T))
bvec = c(5,3,0,0)
solve.QP(Dmat,dvec,Amat,bvec)


# Example 3: partially binding constraints that change only one optimal value
Dmat = matrix(c(2,0,0,8),2,2,byrow=T)
dvec = c(8,16)
Amat = t(matrix(c(1,1,1,0,1,0,0,1),4,2,byrow=T))
bvec = c(5,4.5,0,0)
solve.QP(Dmat,dvec,Amat,bvec)

# Example 4: partially binding constraints that change both optimal values
Dmat = matrix(c(2,0,0,8),2,2,byrow=T)
dvec = c(8,16)
Amat = t(matrix(c(1,1,1,0,1,0,0,1),4,2,byrow=T))
bvec = c(7,3,0,0)
solve.QP(Dmat,dvec,Amat,bvec)

# Example 5: totally binding constraints
Dmat = matrix(c(2,0,0,8),2,2,byrow=T)
dvec = c(8,16)
Amat = t(matrix(c(1,1,1,0,1,0,0,1),4,2,byrow=T))
bvec = c(10,7,0,0)
solve.QP(Dmat,dvec,Amat,bvec)

##############
# Extra plot #
##############

# plot the quadratic form, i.e., sufface of y = f(x1,x2)

# unique global minimum point (positive definite Hessian)
x1 = seq(-10, 10, 1)
x2 = x1
f = function(x1, x2) {return(x1^2+4*x2^2)}
y = outer(x1, x2, f)
op = par(bg = "white")

par(mfrow=c(2,2))
persp(x1, x2, y, theta = 10, phi = 30, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 60, phi = 20, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 100, phi = 30, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 150, phi = 20, expand = 0.5, col = "lightblue")


# minumum path (Hessian is not positive definite)
x1 = seq(-10, 10, 1)
x2 = x1
f = function(x1, x2) {return(x1^2+4*x2^2+4*x1*x2)}
y = outer(x1, x2, f)
op = par(bg = "white")

par(mfrow=c(2,2))
persp(x1, x2, y, theta = 10, phi = 30, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 60, phi = 20, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 100, phi = 30, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 150, phi = 20, expand = 0.5, col = "lightblue")

# Saddle point (Hessian is indefinite)
x1 = seq(-10, 10, 1)
x2 = x1
f = function(x1, x2) {return(x1^2-4*x2^2)}
y = outer(x1, x2, f)
op = par(bg = "white")

par(mfrow=c(2,2))
persp(x1, x2, y, theta = 10, phi = 30, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 60, phi = 20, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 100, phi = 30, expand = 0.5, col = "lightblue")
persp(x1, x2, y, theta = 150, phi = 20, expand = 0.5, col = "lightblue")


##############################################
# Optional: Sequential Quadratic Programming #
##############################################
# page 13, example 5.1 from https://ecal.berkeley.edu/files/ce191/CH02-QuadraticProgramming.pdf
# https://cran.r-project.org/web/packages/quadprog/quadprog.pdf

x1 = 1; x2 = 1
i = 1
while (i < 7) {
cat("--------------------------------------", "\n")
cat("Solutions are ", c(x1,x2), "\n")
cat("f and g are ", c(exp(-x1)+(x2-2)^2,x1*x2-1), "\n")
Dmat = matrix(c(exp(-x1),0,0,2),2,2,byrow=T)
dvec = matrix(c(exp(-x1),-2*(x2-2)),2,1)
Amat = matrix(c(-x2,-x1),2,1)
bvec = c(x1*x2-1)
ou = solve.QP(Dmat,dvec,Amat,bvec)
x1 = ou$solution[1]+x1
x2 = ou$solution[2]+x2
i = i + 1
}

###############################
# Optional: OLS as QP problem #
###############################
# canned lm function
ad = "https://www.fsb.miamioh.edu/lij14/400_house.txt"
data = read.table(url(ad), header=T)  
rprice = data[,4]
age = data[,1]
area = data[,2]
lm(rprice~age+area)$coef

# QP
x = cbind(1, age,area)
y = rprice
xprimx = t(x)%*%x
xprimy = t(x)%*%y

library(quadprog)
Dmat = matrix(2*xprimx,3,3,byrow=T)
dvec = 2*xprimy
Amat = matrix(rep(0,9),3,3,byrow=T)
bvec = rep(0,3)
solve.QP(Dmat,dvec,Amat,bvec)

###########################################
# Optional: Constrained OLS as QP problem #
###########################################
# imposing the restriction that beta1 = -10*beta2
x = cbind(1, age,area)
y = rprice
xprimx = t(x)%*%x
xprimy = t(x)%*%y

library(quadprog)
Dmat = matrix(2*xprimx,3,3,byrow=T)
dvec = 2*xprimy
Amat = t(matrix(c(0,1,10),1,3,byrow=T))
bvec = c(0)
solve.QP(Dmat,dvec,Amat,bvec,meq=1)