desirability.R
In desirability: Desirabiliy Function Optimization and Ranking

### R code from vignette source 'desirability.Rnw'

###################################################
### code chunk number 1: makeDesirePlots
###################################################
library(desirability)
pdf("dMax.pdf", width = 7, height = 4)
   par(mfrow=c(1,3))
   
   plot(dMax(1, 3), nonInform = FALSE)
   
   plot(dMax(1, 3, 5),  TRUE, col = 2, nonInform = FALSE)
   text(2.73, .3, "scale = 5", col = 2) 

   plot(dMax(1, 3, 1/5),  TRUE, col = 4, nonInform = FALSE) 
   text(1.3, .8, "scale = .2", col = 4)   
   
   text(2, .62, "scale = 1", col = 1)   
   title("(a)")   

   plot(dMin(1, 3), nonInform = FALSE)
   
   plot(dMin(1, 3, 5),  TRUE, col = 2, nonInform = FALSE)
   text(1.5, 0.1, "scale = 5", col = 2)   
  
   plot(dMin(1, 3, 1/5),  TRUE, col = 4, nonInform = FALSE)  
   text(1.5, 1, "scale = .2", col = 4)  
    
   text(2, .62, "scale = 1", col = 1)   
   title("(b)")   

   plot(dTarget(1, 2, 3), nonInform = FALSE)
   
   plot(dTarget(1, 2, 3, 5),  TRUE, col = 2, nonInform = FALSE)
   text(1.9, .1, "lowScale = 5", col = 2) 

   plot(dTarget(1, 2, 3, 1/5),  TRUE, col = 4, nonInform = FALSE)      
   text(1.3, 0.9, "lowScale = .2", col = 4)   
   
   text(1.35, 0.6, "lowScale = 1", col = 1)   
   title("(c)")   
dev.off()
par(mfrow=c(1,1))


###################################################
### code chunk number 2: responseEq
###################################################
conversionPred <- function(x) 81.09 + 1.0284 * x[1] + 4.043 * x[2] + 6.2037 * x[3] - 
   1.8366 * x[1]^2 + 2.9382 * x[2]^2 - 5.1915 * x[3]^2 +
   2.2150 * x[1] * x[2] + 11.375 * x[1] * x[3] - 3.875 * x[2] * x[3]

activityPred <- function(x) 59.85 + 3.583 * x[1] + 0.2546 * x[2] + 2.2298 * x[3] + 
   0.83479 * x[1]^2 + 0.07484 * x[2]^2 + 0.05716 * x[3]^2 -
   0.3875 * x[1] * x[2] - 0.375 * x[1] * x[3] + 0.3125 * x[2] * x[3]


###################################################
### code chunk number 3: setupData
###################################################
plotGrid <- expand.grid(time = seq(-1.7, 1.7, length = 20), temperature = seq(-1.7, 1.7, length = 4), catalyst = seq(-1.7, 1.7, length = 20))
plotGrid$conversionPred <- apply(plotGrid[, 1:3], 1, conversionPred)
plotGrid$activityPred <- apply(plotGrid[, 1:3], 1, activityPred)


###################################################
### code chunk number 4: conversionPlot
###################################################
library(lattice)
print(
contourplot(conversionPred ~ time + catalyst|temperature, plotGrid, aspect = 1))


###################################################
### code chunk number 5: activityPlot
###################################################

print(
contourplot(activityPred ~ time + catalyst|temperature, plotGrid, aspect = 1))


###################################################
### code chunk number 6: makeDesire
###################################################
conversionD <- dMax(80, 97)
activityD <- dTarget(55, 57.5, 60)
predOutcomes <- c(conversionPred(c(0,0,0)), activityPred(c(0,0,0)))
print(predOutcomes)
predict(conversionD, predOutcomes[1])
predict(activityD, predOutcomes[2])


###################################################
### code chunk number 7: makeDesire
###################################################
overallD <- dOverall(conversionD, activityD)
print(overallD)
predict(overallD, predOutcomes)


###################################################
### code chunk number 8: setupDPlots
###################################################
dValues <- predict(overallD, plotGrid[, 4:5], TRUE)
names(dValues) <- gsub("Pred", "Desire", names(dValues))
plotGrid <- cbind(plotGrid, dValues)


###################################################
### code chunk number 9: conversionPlot2
###################################################
print(
contourplot(conversionDesire ~ time + catalyst|temperature, plotGrid, aspect = 1))


###################################################
### code chunk number 10: activityPlot2
###################################################
print(
contourplot(activityDesire ~ time + catalyst|temperature, plotGrid, aspect = 1))


###################################################
### code chunk number 11: overallPlot
###################################################
print(
contourplot(Overall ~ time + catalyst|temperature, plotGrid, spect = 1))


###################################################
### code chunk number 12: optFunction
###################################################
rsmOpt <- function(x, dObject, space = "square")
{
   conv <- conversionPred(x)
   acty <- activityPred(x)

   out <- do.call(
      predict, 
      list(object = dObject, newdata = c(conv, acty)))
   if(space == "circular")
   {
      if(sqrt(sum(x^2)) > 1.682) out <- 0
   } else if(space == "square") if(any(abs(x) > 1.682)) out <- 0
   out
}


###################################################
### code chunk number 13: squareRegion
###################################################
searchGrid <- expand.grid(
   time = seq(-1.5, 1.5, length = 5),
   temperature = seq(-1.5, 1.5, length = 5),
   catalyst = seq(-1.5, 1.5, length = 5))

for(i in 1:dim(searchGrid)[1])
{
   tmp <- optim(
      as.vector(searchGrid[i,]), 
      rsmOpt, 
      dObject = overallD,
      space = "square",
      control = list(fnscale = -1))
   if(i == 1)
   {
      best <- tmp
   } else {
      if(tmp$value > best$value) best <- tmp   
   }
}
print(best) 


###################################################
### code chunk number 14: circularRegion
###################################################
for(i in 1:dim(searchGrid)[1])
{
   tmp <- optim(
      as.vector(searchGrid[i,]), 
      rsmOpt, 
      space = "circular",      
      dObject = overallD,      
      control = list(fnscale = -1))
   if(i == 1)
   {
      best <- tmp
   } else {
      if(tmp$value > best$value) best <- tmp   
   }
}
print(best) 


###################################################
### code chunk number 15: logistic
###################################################
foo <- function(u) 1/(1+exp(-u))
xInput <- seq(-5,5, length = 20)
logisticD <- dArb(xInput, foo(xInput))


###################################################
### code chunk number 16: logisticPlot
###################################################
pdf("logistic.pdf", width = 5, height = 4)
   plot(logisticD)
dev.off()


###################################################
### code chunk number 17: box
###################################################
pdf("box.pdf", width = 5, height = 4)
plot(dBox(-1.682, 1.682), nonInform = FALSE)
dev.off()