##################################
### call hello

rcpp_hello_world()

##################################################
## simulate data
n = 500 #sample size
beta = c(1,.5) # true (intercept, slope)

## storage for simulation results
## py1 = F(beta[1] + beta[2]x), y ~ Bern(py1)
x = rep(0,n)
y = rep(0,n)
py1 = rep(0,n)

## simulate
set.seed(34)
for(i in 1:n) {
  x[i] = rnorm(1)
  py1[i] = 1/(1 + exp(-(beta[1] + beta[2]*x[i])))
  y[i] = rbinom(1,1,py1[i])
}

############################################################
## call rcmll

junk = rcmll(x,y,beta)

###########################################################
### time against pure R version

## R version
mLL = function(x,y,beta) {
  py1 = 1/(1+exp(-(beta[1] + beta[2]*x)))
  return(-sum(ifelse(y==1,log(py1),log(1-py1))))
}

## should all be the same
cat("\n##### R and C++ mll should be the same:\n")
print(mLL(x,y,beta))
print(junk$mll)

## big grid of (b0,b1) to evaluate mll on
nval=500
b1g = seq(from=0,to=2,length.out=nval) 
b2g = seq(from=0,to=1,length.out=nval) 
## bg will be bivariate grid, all vals of b1g x all vals of b2g
bg = expand.grid(b1g,b2g)
nn = nrow(bg)
bgM = as.matrix(bg)
bgMM = cbind(bgM[,1],bgM[,2])

##  R version
llv1 = rep(0,nn)
tm1 = system.time({
  for(i in 1:nn) {
    llv1[i] = mLL(x,y,bgMM[i,])
  }
})

## C++ Rcpp version
llv2 = rep(0,nn)
tm2 = system.time({
  for(i in 1:nn) {
    llv2[i] = rcmll(x,y,bgMM[i,])$mll
  }
})

cat("\n&&&& time for vectorized R version\n")
print(tm1)

cat("\n&&&& time for C++ version\n")
print(tm2)

cat("\n$$$$$ difference in mll from R and C++\n")
print(summary(llv1-llv2))