options(amelia.parallel="no",
        amelia.ncpus=1)
library(matrixStats) # for numerically stable logsumexps
source("../misclassificationmodels/R/likelihoods.R")
source("../misclassificationmodels/R/glm_fixit.R")

## ... are the parts of the formular.
run_simulation <-  function(..., family=gaussian(),
                      proxy_family = binomial(link='logit'),
                      truth_family = binomial(link='logit'),
                      data,
                      data2,
                      result
                      ){

    accuracy <- df[,mean(w_pred==x)]
    accuracy.x.y0 <- df[y<=0,mean(w_pred==x)]
    accuracy.x.y1 <- df[y>=0,mean(w_pred==x)]
    accuracy.z.y0 <- df[y<=0,mean(a_pred==z)]
    accuracy.z.y1 <- df[y>=0,mean(a_pred==z)]
    cor.y.xi <- cor(df$x - df$w_pred, df$y)
    cor.y.zi <- cor(df$z - df$a_pred, df$y)

    fnr <- df[w_pred==0,mean(w_pred!=x)]
    fnr.y0 <- df[(w_pred==0) & (y<=0),mean(w_pred!=x)]
    fnr.y1 <- df[(w_pred==0) & (y>=0),mean(w_pred!=x)]

    fpr <- df[w_pred==1,mean(w_pred!=x)]
    fpr.y0 <- df[(w_pred==1) & (y<=0),mean(w_pred!=x)]
    fpr.y1 <- df[(w_pred==1) & (y>=0),mean(w_pred!=x)]
    cor.resid.w_pred <- cor(resid(lm(y~x+z,df)),df$w_pred)

    result <- append(result, list(accuracy=accuracy,
                                  accuracy.x.y0=accuracy.x.y0,
                                  accuracy.x.y1=accuracy.x.y1,
                                  accuracy.z.y0=accuracy.z.y0,
                                  accuracy.z.y1=accuracy.z.y1,
                                  cor.y.xi=cor.y.xi,
                                  fnr=fnr,
                                  fnr.y0=fnr.y0,
                                  fnr.y1=fnr.y1,
                                  fpr=fpr,
                                  fpr.y0=fpr.y0,
                                  fpr.y1=fpr.y1,
                                  cor.resid.w_pred=cor.resid.w_pred
                                  ))

    result <- append(result, list(cor.xz=cor(df$x,df$z)))
    (model.true <- lm(y ~ x + z, data=df))
    true.ci.Bxy <- confint(model.true)['x',]
    true.ci.Bzy <- confint(model.true)['z',]

    result <- append(result, list(Bxy.est.true=coef(model.true)['x'],
                                  Bzy.est.true=coef(model.true)['z'],
                                  Bxy.ci.upper.true = true.ci.Bxy[2],
                                  Bxy.ci.lower.true = true.ci.Bxy[1],
                                  Bzy.ci.upper.true = true.ci.Bzy[2],
                                  Bzy.ci.lower.true = true.ci.Bzy[1]))
                                  
  ## mle_result <- list(
  ##   Bxy.est.naive = NULL,
  ##   Bxy.ci.upper.naive = NULL,
  ##   Bxy.ci.lower.naive = NULL,
  ##   Bzy.est.naive = NULL,
  ##   Bzy.ci.upper.naive = NULL,
  ##   Bzy.ci.lower.naive = NULL,
  ##   Bxy.est.feasible = NULL,
  ##   Bxy.ci.upper.feasible = NULL,
  ##   Bxy.ci.lower.feasible = NULL,
  ##   Bzy.est.feasible = NULL,
  ##   Bzy.ci.upper.feasible = NULL,
  ##   Bzy.ci.lower.feasible = NULL,
  ##   Bxy.est.mle = NULL,
  ##   Bxy.ci.upper.mle = NULL,
  ##   Bxy.ci.lower.mle = NULL,
  ##   Bzy.est.mle = NULL,
  ##   Bzy.ci.upper.mle = NULL,
  ##   Bzy.ci.lower.mle = NULL
  ## )
  ## tryCatch({
    mod <- glm_fixit(..., family=family, data=data, data2=data2)
  
    coef_corrected <- coef(mod, which_model = "corrected")
    ci_corrected <- confint(mod, which_model = "corrected")
    
    coef_feasible <- coef(mod, which_model = "feasible")
    ci_feasible <- confint(mod, which_model = "feasible")
    
    coef_naive <- coef(mod, which_model = "naive")
    ci_naive <- confint(mod, which_model = "naive")
    
    mle_result <- list(
      Bxy.est.mle = coef_corrected["x"],
      Bxy.ci.upper.mle = ci_corrected["x", 2],
      Bxy.ci.lower.mle = ci_corrected["x", 1],
      Bzy.est.mle = coef_corrected["z"],
      Bzy.ci.upper.mle = ci_corrected["z", 2],
      Bzy.ci.lower.mle = ci_corrected["z", 1],
      Bxy.est.feasible = coef_feasible["x"],
      Bxy.ci.upper.feasible = ci_feasible["x", 2],
      Bxy.ci.lower.feasible = ci_feasible["x",1],
      Bzy.est.feasible = coef_feasible["z"],
      Bzy.ci.upper.feasible = ci_feasible["z",2],
      Bzy.ci.lower.feasible = ci_feasible["z",1],
      Bxy.est.naive = coef_naive["w"],
      Bxy.ci.upper.naive = ci_naive["w",2],
      Bxy.ci.lower.naive = ci_naive["w",1],
      Bzy.est.naive = coef_naive["a"],
      Bzy.ci.upper.naive = ci_naive["a",2],
      Bzy.ci.lower.naive = ci_naive["a",1]
    )
  ##     print(mle_result)
  ## },
  ## error=function(e) {result[['error']] <- as.character(e)
  ## })
  result <- append(result, mle_result)
  return(result)
}