library(VGAM)

gompertz <- function (a, b, c, t) {
  return(a * exp(b * exp(c * t)))
}

bling <- function (c, t, interval) c * interval * ceiling(t/interval)

samples <- 100            # number of days worked
ticks   <- 200            # number of time indices
minval  <- 300            # the marginal value of daily production 
shape   <- log(5)/log(4)  # 80-20 rule
#shape   <- log(10)/log(9) # 90% of everything is crap

cost <- 500 # day rate

# A and C parameters to the curve functions are random Pareto numbers
gparam  <- cbind(rpareto(samples, minval, shape),
                 -rpareto(samples, 0.05, shape))

#par(mfcol=c(4,1), mar=c(3, 6, 2, 1))
par(mar=c(3,6,2,1))
png("~/Desktop/individual-values.png", 500, 300, bg="transparent");

# first plot, individual values
plot(1, ann=TRUE, type="n",
     xlim=c(0, ticks), ylim=c(0, max(gparam[,1])),
     main="Individual Values of Created Information Artifacts",
     xlab="Days", ylab="Dollars",
     las=1, mgp=c(4.5, 0.5, 0),
     family="CorpidOffice"
     )
for (i in 1:samples) {
  x <- gompertz(gparam[i,1], -10, gparam[i,2], 1:ticks-i)
  lines(x, type="l", col=hsv(i/samples, 1, 0.5, 0.6)) #col=rgb(0, 0, 0, 0.6))
}
abline(v=samples, col=rgb(.6, 0, 0, .5))

dev.off()
par(mar=c(3,6,2,1))
png("~/Desktop/cost-function.png", 500, 300, bg="transparent");

totalcost <- unlist(lapply(1:ticks,
                           function(t) -ifelse(t<=samples, bling(cost, t, 20),
                                               cost*samples)))
# cost function
plot(# continuous version
     #unlist(lapply(1:ticks,
     #              function (t) -cost*ifelse(t<=samples, t, samples))),
     # stair-step version
     #unlist(lapply(1:ticks,
     #              function(t) -ifelse(t<=samples, bling(cost, t, 20),
     #                                  cost*samples))),
     totalcost,
     type="l", xlim=c(0, ticks), ylim=c(-cost*samples, 0),
     main="Compensation to Creative Professional",
     xlab="Days", ylab="Dollars",
     las=1, mgp=c(4.5, 0.5, 0), col=rgb(0, 0, 0, 0.5),
     family="CorpidOffice"
     )
## lines(
##      unlist(lapply(1:ticks,
##                    function(t) -ifelse(t<=samples, bling(cost, t, 20),
##                                        cost*samples))),
## )
abline(v=samples, col=rgb(.6, 0, 0, .6))
par(family="CorpidOffice")
legend("topright", c("Intervention Ends"),
       lty=1, col=rgb(.6, 0, 0, .6), bty="n", inset=0.0)


dev.off()
par(mar=c(3,6,2,1))
png("~/Desktop/cumulative-effect.png", 500, 300, bg="transparent");

# second plot, cumulative
totalval <- unlist(lapply(1:ticks,
  function (t) sum(unlist(lapply(1:samples,
    function (s) gompertz(gparam[s,1], -10, gparam[s,2], t - s))))
                        ))
                        #+ cost*ifelse(t<=samples, -t, -samples)))
                        #+ ifelse(t<=samples, -bling(cost, t, 20),
                        #         -cost*samples)))

plot(totalval, type="l", col=rgb(0, 0, 0, 0.5),
     main="Cumulative Value of Created Information Artifacts",
     xlab="Days", ylab="Dollars",
     family="CorpidOffice",
     #ylim=c(0, 10^ceiling(log(max(career), 10))),
     las=1, mgp=c(4.5, 0.5, 0),
#     yaxt="n"
     )
#axis(2, career, labels=sprintf("%.3fk", career/1000))
#lines(c2, col=rgb(0, 0, 0, 0.5))
abline(v=samples, col=rgb(.6, 0, 0, .5))

dev.off()
par(mar=c(3,6,2,1))
png("~/Desktop/net-return.png", 500, 300, bg="transparent");

career <- totalval + totalcost

## c2 <- unlist(lapply(1:ticks,
##   function (t) sum(unlist(lapply(1:samples,
##     function (s) gompertz(gparam[s,1], -10, gparam[s,2], t - s))))
##                     + cost*ifelse(t<=samples, -t, -samples)))
##                     #+ ifelse(t<=samples, -bling(cost, t, 20),
##                     #         -cost*samples)))


plot(career, type="l", col=rgb(0, 0, 0, 0.5),
     main="Net Return on Created Information Artifacts",
     xlab="Days", ylab="Dollars",
     family="CorpidOffice",
     #ylim=c(0, 10^ceiling(log(max(career), 10))),
     las=1, mgp=c(4.5, 0.5, 0),
#     yaxt="n"
     )
#axis(2, career, labels=sprintf("%.3fk", career/1000))
#lines(c2, col=rgb(0, 0, 0, 0.5))
abline(h=0, col=rgb(0, 0, 0, .6))
abline(v=samples, col=rgb(.6, 0, 0, .5))

dev.off()