LiF_Decoy_Feasibility_Test.R

library(Cairo)
library(ggplot2)
library(animation)

############################## configuration ##############################


# set working directory

setwd("C:/Users/robin/Desktop/Monopolist")

# do animated plots?

plots <- FALSE

# do spline? if false, we'll do line

is_spline <- TRUE

# if do_spline, which method?

spline_method <- "monoH.FC"  # straight lines and bends through points
spline_method <- "periodic"  # bends, but not necessarly through points, flattens toward sine
spline_method <- "fmm"       # default, always bends, through points

# if not do_spline, which method?

nspline_method <- "linear"   #lines
nspline_method <- "constant"  #blocks

# Number of customers %% 200

customers <- 1000

# LiF settings

amount_features_start <- 45

TT <- customers
T <- 8
w <- (2*pi)/T
A <- 10
gamma <- 0.9
window.length <- 1


############################## helper functions  #########################


push <- function(A, value){
  if(any(is.na(A))){
    i <- sum(!is.na(A))
    A[i+1] <- value
  } else {
    A <- c(A,value)
    A <- A[-1]
  }
  return(A)
}

slideFunct <- function(data, window, step){
  total <- length(data)
  spots <- seq(from=1, to=(total-window), by=step)
  result <- vector(length = length(spots))
  for(i in 1:length(spots)){
    result[i] <- mean(data[spots[i]:(spots[i]+window)])
  }
  return(result)
}


############################## sketch some graphs  #####################


x_decoy_features <- c(30,60)  # A has 30 features B has 60 features
y_decoy_price <- c(30,60)     # A costs 30 B costs 60

plot(x_decoy_features, y_decoy_price, xlim = c(0, 90), ylim = c(90,0), col="red")
text(x_decoy_features+0.2, y_decoy_price-4, labels = c("A","B"), col="red")
abline(h=60,lty = 2)
text(25, 65, labels = "Decoy values of C")
text(25, 71, labels = "along dashed line")


x1   <- c(  0, 30, 60, 90   )  
f1   <- c( .4, .4, .2, .2   )  

x2   <- c(   0, 30, 60, 90  ) 
f2   <- c(   .4, .4, .7, .0 )

x3   <- c(   0, 30, 60, 90  )
f3   <- c(  .2, .2, .1, .8  )


if (is_spline) {
  
  s1 <- splinefun(x1, f1,method=spline_method)
  plot(x1, f1, ylim = c(0, 1))
  curve(s1(x), add = TRUE, col = 2, n = 1001) 
  
  s2 <- splinefun(x2, f2,method=spline_method)
  plot(x2, f2, ylim = c(0, 1))
  curve(s2(x), add = TRUE, col = 2, n = 1001) 
  
  s3 <- function(x){
    return(  1 - s1(x) - s2(x)  )
  }
  plot(x3, f3, ylim = c(0, 1))
  curve(s3(x), add = TRUE, col = 2, n = 1001) 
  
} else {
  
  s1 <- approxfun(x1, f1,method = nspline_method)
  plot(x1, f1, ylim = c(0, 1))
  curve(s1(x), add = TRUE, col = 2, n = 1001) 
  
  s2 <- approxfun(x2, f2,method = nspline_method)
  plot(x2, f2, ylim = c(0, 1))
  curve(s2(x), add = TRUE, col = 2, n = 1001) 
  
  s3 <- function(x){
    return(  1 - s1(x) - s2(x)  )
  }
  plot(x3, f3, ylim = c(0, 1))
  curve(s3(x), add = TRUE, col = 2, n = 1001) 
  
} 

# Plot functions, see if conform sketch

p <- ggplot(data.frame(x=c(0,90)), aes(x)) 
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + stat_function(fun=function(x){s1(x)}, geom="line", colour="red")
p <- p + stat_function(fun=function(x){s1(x)+s2(x)}, geom="line", colour="blue")
p <- p + stat_function(fun=function(x){s1(x)+s2(x)+s3(x)}, geom="line", colour="green")
print (p)

# So, maximize sales of *B* by changing how many features *C* (both at same price)

p <- ggplot(data.frame(x=c(0,90)), aes(x)) 
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + stat_function(fun=function(x){s2(x)}, geom="line", colour="blue")
print (p)


############################## run a test of probs  #####################


# set amount of features to 60

feat = 45

choices_sampling <- sample( c("A","B","C"), customers, replace=TRUE, prob=c(s1(feat), s2(feat), s3(feat)))

# plot barchart of prob

choicetable = prop.table(table(choices_sampling))

barplot(choicetable, col = heat.colors(3), ylim=c(0,1))


################################ Do lif! ################################


# now see if LIF can find the optimum C 
# which optimizes B, where B is most profitable

yws <- rep(NA, T)

ytm <- .3

length_features <- 90

profit_per_sale <- rep(NA,length_features) 

arr.features <- rep(NA, TT)
arr.amount_features_start <- rep(NA,TT)
arr.doesbuy <- rep(NA, TT)
arr.sampleChoices <- rep(NA, TT)

############################## Lif sim loop ################################


for(i in 1:customers){
  
  # Oscillate around features
  features <- amount_features_start + A*cos(w*i)
  
  # This is only necessary for non-spline constant
  if(features > 90) features <- 90
  if(features < 0) features <- 0
  
  sampleChoice <- sample( c("A","B","C"), 1, replace=TRUE, prob=c(s1(features), s2(features), s3(features)))
  
  if (sampleChoice == "A") {
    yt <- 30*0.05     # 5% profit on $30 of A
  } else if (sampleChoice == "B"){
    yt <- 60*0.06     # 6% profit of $60 for B
  } else {
    yt <- 60*0.04     # 4% profit on $60 for C
  }
  
  
  yws <- push(yws, yt * cos(w*i))
  
  # reset amount_features_start using lif
  if(i > T){
    yw <- sum(yws) / T
    amount_features_start <- amount_features_start + (gamma / T) * yw
  }
  
  # Log me!
  profit_per_sale[i] <- yt
  arr.sampleChoices[i] <- sampleChoice
  arr.features[i] <- features 
  arr.amount_features_start[i] <- amount_features_start
  
  # Want animations?
  if (i %% 200 == 0 && plots) { 
    Cairo(file = paste0("S", i/200, ".png"), bg = "white", width = 1300,height = 900)  
    par(mfrow=c(2,2))
    plot(arr.features, type="l",ylim=c(0,90))
    plot(arr.amount_features_start, type="l", xlab="Time in stream", ylab="features",ylim=c(0,90))
    lines(arr.amount_features_start, col="red",ylim=c(0,90))
    dev.off()  
    
  }
}

############################## Model profit ###################################

profit_plot <- rep(NA,length_features) 
for (j in 1:length_features) profit_plot[j] <- s1(j)*30*0.05 + s2(j)*60*0.06  +s3(j)*60*0.04 


############################## Result #########################################

max_profit_model = which(profit_plot==max(profit_plot))

print(paste0("Maximum profit at feature X of profit function: ", max_profit_model))

print(paste0("Maximum profit at feature X found by LiF ", amount_features_start))


############################## do some plotting  ##############################

# Do some simple plots



plot(arr.features, type="l",ylim=c(0,90))

plot(profit_plot, type="l")
lines(c(max_profit_model[1],max_profit_model[1]),c(-5,max(profit_plot)), col="red")

plot(slideFunct(profit_per_sale,100,50), type="l")

plot(arr.amount_features_start, type="l", xlab="Time in stream", ylab="features",ylim=c(0,90))
lines(arr.amount_features_start, col="red",ylim=c(0,90))


# Do some animated plots
if (plots) {
  ani.options(convert = shQuote("convert.exe"),ani.width = 1300, ani.height = 900, interval = 0.05, ani.dev = "png",ani.type = "png")
  bm.files = sprintf("ani/S%i.png", 1:(customers/200) )
  im.convert(bm.files, output = paste0("ani",format(Sys.time(), "%d%b%H%M%S"),".gif"), extra.opts = "-dispose Background", clean = T)  
}