source("ols_lfn.R")
source("simulated_setup.R")

# There are really only two distinct likelihood functions: lf1() and
# lf1.inlogs. The following test cases are desirable --
#                   BFGS      L-BFGS-B     Nelder-Mead      SANN
# lf1                 a,n         a,n            yes          yes
# lf1.inlogs           n           NO             yes          yes

# This function does measurement about the hit rate on reaching
# the correct answer, and the time taken (in ms), of alternative MLE
# algorithms.
# You send in args specific to the method or situation using "..."
measurement <- function(method, lf, trans=FALSE, ...) {
  N <- 100                              # How many tests?
  success <- 0
  costs <- system.time(
                       for (i in 1:100) {
                         seeds <- 10*runif(3) # random seeds from 0 to 10
                         if (trans) {
                           seeds[1] <- log(seeds[1])
                         }
                         optimised <- optim(seeds, lf, method=method,
                                            y=y, X=X, ...)
                         if (trans) {
                           optimised$par[1] <- exp(optimised$par[1])
                         }
                         if (SSE(optimised$par) < 0.001)
                           {success <- success + 1}
                       }
                       )
  return(list(successrate=100*success/N, cost=1000*costs[1]/N))
}

performance <- matrix(NA, nrow=9, ncol=2)
colnames(performance) <- c("Hit rate", "Cost")
rownames(performance) <- rep("Method = ", 9)
i <- 1

# Start of 9 runs --
a <- measurement(method="BFGS", lf=ols.lf1)
rownames(performance)[i] <- "BFGS, numerical"
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="BFGS", lf=ols.lf1, gr=ols.gradient)
rownames(performance)[i] <- "BFGS, analytical"
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="L-BFGS-B", lf=ols.lf1,
                 lower=c(1e-6,-Inf,-Inf), upper=rep(Inf,3))
rownames(performance)[i] <- "L-BFGS-B, numerical"
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="L-BFGS-B", lf=ols.lf1,
                 lower=c(1e-6,-Inf,-Inf), upper=rep(Inf,3), gr=ols.gradient)
rownames(performance)[i] <- "L-BFGS-B, analytical"
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="BFGS", lf=ols.lf1.inlogs, trans=TRUE)
rownames(performance)[i] <- "BFGS, trans., numerical"
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="Nelder-Mead", lf=ols.lf1)
rownames(performance)[i] <- "Nelder-Mead"
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="SANN", lf=ols.lf1)
rownames(performance)[i] <- "SANN"
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="Nelder-Mead", lf=ols.lf1.inlogs, trans=TRUE)
rownames(performance)[i] <- "Nelder-Mead, trans."
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

a <- measurement(method="SANN", lf=ols.lf1.inlogs, trans=TRUE)
rownames(performance)[i] <- "SANN, trans."
performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost
i <- i + 1

# Print out summary statistics
o <- order(performance[,2])
performance[o,]