# R code to accompany Mixture Design and Analysis video by Researchnology.com
library(mixexp)
library(nloptr)   # For nonlinear optimization
library(openxlsx)

##### Step1: Construct the design
design <- Xvert(4, uc=c(0.9, 0.9, 0.5, 0.9), lc=c(0.1,0,0.1,0), ndm=1, plot=TRUE)
design <- signif(design,3)
# Visualize the subset of the design
DesignPoints(design, cornerlabs=c("盐", "酸菜", "鸡汤酱", "香油"),
             axislabs=c("鸡汤酱", "酸菜", "盐", "香油"))

# If interested can fill more interior points, not for this video demo
tempDesign <- Fillv(4, design)
DesignPoints(tempDesign, cornerlabs=c("盐", "酸菜", "鸡汤酱", "香油"),
             axislabs=c("鸡汤酱", "酸菜", "盐", "香油"))

##### Step 2: Data analysis
# Load experiment data with design and result
spack <- read.csv("data/spice_mix.csv")
spack
mdl <- MixModel(spack, "y", mixcomps=c("x1", "x2", "x3", "x4"), model=2)
summary(mdl)

ModelPlot(model=mdl,
          dimensions=list(x1="x1", x2="x2", x3="x3"),
          slice=list(mix.vars = c(x4=0)),
          main =  "方便面料包顾客满意度预测",
          axislabs=c("鸡汤酱比例","酸菜比例","盐比例"),
          axislab.pars = list(),
          cornerlabs = c("鸡汤酱", "酸菜", "盐"),
          lims = c(0,0.9,0,0.9,0.1,0.5),
          constraints=TRUE, contour = TRUE, cuts=12, fill = TRUE, 
          pseudo = FALSE)

# Model effect plot
ModelEff(nfac = 4, mod = 2, dir = 2, ufunc = mdl,
         dimensions = list("x1", "x2", "x3", "x4"))

##### Step 3: Optimization and prediction
coef <- mdl$coeff
names(coef)
paste0(round(coef,3), names(coef), sep="+", collapse = "")

# The objective function
obj <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  x3 <- x[3]
  x4 <- x[4]
  fn <- (10.996*x1+13.292*x2+28.190*x3+20.094*x4+23.602*x1*x2-55.858*x1*x3-
           27.789*x1*x4-49.381*x2*x3-3.4*x2*x4-56.179*x3*x4)
}

# Equality constraints
eval_g_eq <- function(x) {
  return ( x[1] + x[2] + x[3] + x[4] - 1 )
}

# Inequality constraints
eval_g_ineq <- function (x) {
  constr <- c(x[2] + x[3] - 0.8,
              x[1] + x[4] - 0.8,
              0.1 - x[1] - x[4])
  return (constr)
}

# Lower and upper bounds
lb <- c(0.1, 0, 0.1, 0)
ub <- c(0.9,0.9,0.5,0.9)

#initial values
x0 <- c(0.3, 0.3, 0.3, 0.3)

# Set optimization options.
local_opts <- list( "algorithm" = "NLOPT_LD_MMA", "xtol_rel" = 1.0e-15 )
opts <- list("algorithm"= "NLOPT_GN_ISRES",
             "xtol_rel"= 1.0e-15,
             "maxeval"= 160000,
             "local_opts" = local_opts,
             "print_level" = 0 )

# Optimization
res <- nloptr (x0 = x0, eval_f = obj, lb = lb, ub = ub, 
               eval_g_eq = eval_g_eq, eval_g_ineq = eval_g_ineq, opts = opts)
res

# The predicted taste ranking under optimal condition
ModelPlot(model=mdl,
          dimensions=list(x1="x1", x2="x2", x3="x3"),
          slice=list(mix.vars = c(x4=0.188)),
          main =  "方便面料包满意度最佳配方",
          axislabs=c("鸡汤酱比例","酸菜比例","盐比例"),
          axislab.pars = list(),
          cornerlabs = c("鸡汤酱", "酸菜", "盐"),
          lims = c(0,0.9,0,0.9,0.1,0.5),
          constraints=TRUE, contour = TRUE, cuts=12, fill = TRUE, 
          pseudo = FALSE)