R/npsf-internal.R

In npsf: Nonparametric and Stochastic Efficiency and Productivity Analysis

.prepareYX <- function(formula, data, subset, rts = c("C", "NI", "V"),
                       base = c("output", "input"), ref = NULL,
                       data.ref = NULL, subset.ref = NULL, print.level = 1,
                       type = "RM", winw = 50, rts.subst = NULL, sysnframe = 1, ...)
{
 # get y and x matrices

 # mf0 <- match.call(expand.dots = FALSE, call = sys.call(which = 1))
 # needed.frame <- sys.nframe() - 1
 needed.frame <- sys.nframe() - sysnframe
 # cat("sys.nframe() = ",sys.nframe(),"\n")
 # cat("needed.frame = ",needed.frame,"\n")
 # cat("sysnframe = ",sysnframe,"\n")
 # this is the entire call
 # this one is from previous call
 # mf0 <- match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = needed.frame)))
 # print(mf0)
 # this one is from the initial call (where tenonradial might have been called)
 mf0 <- match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = needed.frame)))
 # print(mf0)
 # print( match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = sysnframe+2))) )
 # for(i in 1:1){
 #  cat("i=",i,"\n")
 #  mf0 <- match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = i)))
 #  mf <- mf0
 #  m <- match(c("formula", "data", "subset"), names(mf), 0L)
 #  mf <- mf[c(1L, m)]
 #  mf[[1L]] <- as.name("model.frame")
 #  # mf$formula <- Formula( formula )
 #  # mf <- eval(mf, parent.frame(n = 1))
 #  print(mf)
 # }
 # cat("end\n")

 # if(print.level >= 1) print(mf0)

 # check if it is a matrix
 datasupplied <- !(match("data", names(mf0), 0) == 0)
 # subsetsupplied <- !(match("subset", names(mf0), 0) == 0)
 # print(subsetsupplied)
 # cat("1\n")

 # if data are supplied
 if(datasupplied){
  # begin get a logical vector equal TRUE if !missing

  # first using data and subset to get x without NA
  mf <- mf0
  m <- match(c("formula", "data", "subset"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  # mf$formula <- Formula( formula )
  # print(mf)
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  # mf <- eval(mf, parent.frame())
  mt <- attr(mf, "terms")
  x <- as.matrix(model.matrix(mt, mf))
  # print(x)
  # now get the names in the entire data
  # esample <- seq_len( nrow(data) ) %in% as.numeric(rownames(x))
  esample <- rownames(data) %in% rownames(x)

  # print(rownames(data) )
  # print(rownames(x))
  #
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing

  # get the data
  y <- tryCatch( as.matrix(model.part(Formula(formula), data = data[esample,], lhs = 1)), error = function(e) e )
  if(inherits(y, "error")) stop("Most probably you supplied both matricies and data")
  x <- as.matrix(model.matrix(Formula(formula), data = data[esample,], rhs = 1)[,-1])

  # print(y)
  # print(x)
 }
 # if data are not supplied
 else {
  mf <- mf0
  subsetsupplied <- !(match("subset", names(mf), 0) == 0)
  if(subsetsupplied) stop("Subset with matrices is not allowed")
  m <- match(c("formula"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  # mf <- eval(mf, parent.frame())
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  y <- as.matrix(model.response(mf))
  # cat.print(y)
  x <- as.matrix(model.matrix(mt, mf)[,-1, drop = FALSE])
  # cat.print(x)
  # cat.print(as.matrix(model.matrix(mt, mf)))
  # get a logical vector equal TRUE if !missing
  with.na <- model.frame(mt, na.action = na.pass)
  esample <- rownames(with.na) %in% rownames(x)
  # print(table(esample))
  # print(y)
  # print(x)
 }

 if( !is.numeric(y) ) stop("Some of the outputs are not numeric.")
 if( !is.numeric(x) ) stop("Some of the inputs are not numeric.")
 
 if(type == "RM"){
  # check negative
  x.negative <- apply(x, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
  # cat.print(x.negative)
  # cat.print(x)
  if(sum(x.negative) == nrow(x)){
   stop("Negative values in at least one of the inputs for each data point", call. = FALSE)
  }
  y.negative <- apply(y, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
  # cat.print(y.negative)
  # cat.print(y)
  if(sum(y.negative) == nrow(y)){
   stop("Negative values in at least one of the outputs for each data point", call. = FALSE)
  }
  # deal with negative
  if(sum(x.negative) > 0){
   warning(paste0("There are negative values in inputs for ",sum(x.negative)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  if(sum(y.negative) > 0){
   warning(paste0("There are negative values in outputs for ",sum(y.negative)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  x <- x[!x.negative & !y.negative,,drop = FALSE]
  y <- y[!x.negative & !y.negative,,drop = FALSE]
 } else {
  # check nonpositive
  x.nonpositive <- apply(x, MARGIN = 1, FUN = function(q) sum(q<=0)) >= 1
  if(sum(x.nonpositive) == nrow(x)){
   stop("Nonpositive values in at least one of the inputs for each data point", call. = FALSE)
  }
  y.nonpositive <- apply(y, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
  if(sum(y.nonpositive) == nrow(y)){
   stop("Nonpositive values in at least one of the outputs for each data point", call. = FALSE)
  }
  # deal with negative
  if(sum(x.nonpositive) > 0){
   warning(paste0("There are nonpositive values in inputs for ",sum(x.nonpositive)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  if(sum(y.nonpositive) > 0){
   warning(paste0("There are nonpositive values in outputs for ",sum(y.nonpositive)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  x <- x[(!x.nonpositive & !y.nonpositive),,drop = FALSE]
  y <- y[(!x.nonpositive & !y.nonpositive),,drop = FALSE]
 }
 
 rts <- rts[1]
 base <- base[1]

 rts1 <- tolower(substr(rts, 1,1 ))
 base1 <- tolower(substr(base, 1,1 ))

 if(is.null(rts.subst)){
  if(rts1 == "c" | rts == 3){
   myrts <- 3
   myrts1 <- "CRS"
  } else if (rts1 == "n" | rts == 2){
   myrts <- 2
   myrts1 <- "NIRS"
  } else if (rts1 == "v" | rts == 1){
   myrts <- 1
   myrts1 <- "VRS"
  } else {
   stop("invalid 'rts'; 'rts' must be 'CRS', 'NIRS', or 'VRS'")
  }
 }
 else {
  myrts1 <- rts.subst
 }

 if (base1 == "o" | base == 2){
  mybase <- 2
  mybase1 <- "output"
 } else if (base1 == "i" | base == 1){
  mybase <- 1
  mybase1 <- "input"
 } else {
  stop("invalid 'base'; 'base' must be 'output' or 'input'")
 }

 # for printing
 if(print.level >= 1 & winw > 50){
  mymesage <- paste("\n",ifelse(type == "RM", "Nonradial (Russell)", "Radial (Debrue-Farrell)")," ",mybase1,"-based measures of technical efficiency under assumption of ",myrts1," technology are computed for the following data:\n", sep = "")
  cat("",unlist(strsplit(mymesage, " ")),"", sep = " ", fill = winw-10 )
  # 		cat("\n ",ifelse(type == "RM", "Nonradial (Russell)", "Radial Debrue-Farrell")," ",mybase1,"-based measures of technical efficiency \n", sep = "")
  # 		cat(" under assumption of ",myrts1," technology are computed for the \n", sep = "")
  # 		cat(" following data:\n\n", sep = "")
  cat("  Number of data points (K) = ",nrow(y),"\n", sep = "")
  cat("  Number of outputs     (M) = ",ncol(y),"\n", sep = "")
  cat("  Number of inputs      (N) = ",ncol(x),"\n", sep = "")
 }

 # get y_ref and x_ref matrices

 if(!is.null(ref)){
  mf <- mf0

  # check if it is a matrix
  datasupplied <- !(match("data.ref", names(mf), 0) == 0)

  # if data are supplied
  if(datasupplied){
   N_all_ref <- nrow(data.ref)
   if(N_all_ref == 0) warning("Provided data for reference set does not have a single data point", call. = FALSE)

   # begin get a logical vector equal TRUE if !missing

   # first using data and subset to get x without NA
   mf <- mf0
   # cat.print(mf)
   m <- match(c("ref", "data.ref", "subset.ref"), names(mf), 0L)
   # cat.print(m)
   mf <- mf[c(1L, m)]
   # cat.print(mf)
   # change the names for eval
   names(mf)[which(names(mf) == "ref")] <- "formula"
   names(mf)[which(names(mf) == "data.ref")] <- "data"
   names(mf)[which(names(mf) == "subset.ref")] <- "subset"
   mf[[1L]] <- as.name("model.frame")
   # mf$formula <- Formula( formula )
   mf <- eval(mf, parent.frame())
   # cat.print(mf)
   mt <- attr(mf, "terms")
   # cat.print(mt)
   x_ref <- as.matrix(model.matrix(mt, mf))
   # cat.print(x_ref)
   if(length(x_ref) == 0){
    warning(" Given 'subset' reference set contains zero data points", call. = FALSE)
    # warning(" Reference set will be based on observations,")
    # warning(" for which efficiency is to be computed")
    y_ref <- y
    x_ref <- x
    esample_ref <- NULL
   } else {
    # now get the names in the entire data
     esample_ref <- rownames(data.ref) %in% rownames(x_ref)
    # esample_ref <- seq_len( N_all_ref ) %in% as.numeric(rownames(x_ref))
    # print(table(esample_ref))
    # end get a logical vector equal TRUE if !missing
    # cat.print(ref)
    # cat.print(formula)
    # cat.print(data.ref[esample_ref,])
    # cat.print(data[esample,])
    y_ref <- as.matrix(model.part(Formula(ref), data = data.ref[esample_ref,], lhs = 1))
    # y_ref <- as.matrix(model.part(Formula(formula), data = data[esample,], lhs = 1))
    x_ref <- as.matrix(model.matrix(Formula(ref), data = data.ref[esample_ref,], rhs = 1)[,-1])
   }
   # cat.print(y_ref)
   # cat.print(x_ref)
  }
  # if data are not supplied
  else {
   mf <- mf0
   subsetsupplied <- !(match("subset.ref", names(mf), 0) == 0)
   if(subsetsupplied) stop("Subset for reference set with matrices is not allowed. \n   Or: \nOption 'data.ref' must be provided if option 'ref' is provided", call. = FALSE)
   m <- match(c("ref", "data.ref", "subset.ref"), names(mf), 0L)
   mf <- mf[c(1L, m)]
   names(mf)[which(names(mf) == "ref")] <- "formula"
   mf[[1L]] <- as.name("model.frame")
   mf <- eval(mf, parent.frame())
   mt <- attr(mf, "terms")
   y_ref <- as.matrix(model.response(mf))
   x_ref <- as.matrix(model.matrix(mt, mf)[,-1])
   # print(y_ref)
   # print(x_ref)
   # get a logical vector equal TRUE if !missing
   with.na <- model.frame(mt, na.action = na.pass)
   esample_ref <- rownames(with.na) %in% rownames(x_ref)
   # print(table(esample_ref))
  }
  # print(x_ref)
  # print(y_ref)
  # if reference set not provided
  # for printing
  # 		if(print.level >= 1){
  # 			if(!is.null(esample_ref)){
  # 				cat("\n Reference set is composed as follows:\n\n", sep = "")
  # 				cat(" Number of observations (K) = ",nrow(y_ref),"\n", sep = "")
  # 				cat(" Number of outputs      (M) = ",ncol(y_ref),"\n", sep = "")
  # 				cat(" Number of inputs       (N) = ",ncol(x_ref),"\n\n", sep = "")
  # 			}
  # 		}
  if(!is.null(y_ref)){
   if(ncol(y_ref) != ncol(y)) stop("Number of outputs for data points and reference set must be the same.")
   if(ncol(x_ref) != ncol(x)) stop("Number of inputs for data points and reference set must be the same.")
  }
  
  if(type == "RM"){
   # check negative
   x.negative <- apply(x_ref, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
   if(sum(x.negative) == nrow(x_ref)){
    stop("Negative values in at least one of the reference inputs for each data point", call. = FALSE)
   }
   y.negative <- apply(y_ref, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
   if(sum(y.negative) == nrow(y_ref)){
    stop("Negative values in at least one of the reference outputs for each data point", call. = FALSE)
   }
   # deal with negative
   if(sum(x.negative) > 0){
    warning(paste0("There are negative values in reference inputs for ",sum(x.negative)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
    cat("\n")
   }
   if(sum(y.negative) > 0){
    warning(paste0("There are negative values in reference outputs for ",sum(y.negative)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
    cat("\n")
   }
   x_ref <- x_ref[!x.negative & !y.negative,,drop = FALSE]
   y_ref <- y_ref[!x.negative & !y.negative,,drop = FALSE]
  } else {
   # check nonpositive
   x.nonpositive <- apply(x_ref, MARGIN = 1, FUN = function(q) sum(q<=0)) >= 1
   if(sum(x.nonpositive) == nrow(x)){
    stop("Nonpositive values in at least one of the reference inputs for each data point", call. = FALSE)
   }
   y.nonpositive <- apply(y_ref, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
   if(sum(y.nonpositive) == nrow(y_ref)){
     print(y_ref)
    stop("Nonpositive values in at least one of the reference outputs for each data point", call. = FALSE)
   }
   # deal with negative
   if(sum(x.nonpositive) > 0){
    warning(paste0("There are nonpositive values in reference inputs for ",sum(x.nonpositive)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
    cat("\n")
   }
   if(sum(y.nonpositive) > 0){
    warning(paste0("There are nonpositive values in outputs for ",sum(y.nonpositive)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
    cat("\n")
   }
   x_ref <- x_ref[!x.nonpositive & !y.nonpositive,,drop = FALSE]
   y_ref <- y_ref[!x.nonpositive & !y.nonpositive,,drop = FALSE]
  }
  
 } else {
  y_ref <- y
  x_ref <- x
  esample_ref <- NULL
 }

 if(print.level >= 1  & winw > 50){
  if(is.null(esample_ref)){
   mymesage <- paste("\nData for reference set are not provided. Reference set is formed by ",nrow(y)," data ",ngettext(nrow(y), "point", "point(s)"), " for which measures of technical efficiency are computed", sep = "")
   cat("",unlist(strsplit(mymesage, " ")),"", sep = " ", fill = winw-10 )
   #     cat("\n Data for reference set are not provided. Reference set is formed by ",nrow(y),"\n", sep = "")
   #     cat(" data points, for which measures of technical efficiency are computed.\n\n", sep = "")
  }
  else {
   mymesage <- paste("\nReference set is formed by ",nrow(y_ref)," provided reference data ",ngettext(nrow(y_ref), "point", "point(s)"), "", sep = "")
   cat("",unlist(strsplit(mymesage, " ")),"", sep = " ", fill = winw-10 )
   # cat("\n Reference set is formed by ",nrow(y_ref)," provided reference data points.\n\n", sep = "")
  }
 }

 tymch <- list(y = y, x = x, y.ref = y_ref, x.ref = x_ref, esample = esample, esample.ref = esample_ref, myrts = myrts, rts.string = myrts1, mybase = mybase, base.string = mybase1)
 class(tymch) <- "npsf"
 return(tymch)
}

.prepareYXnoRef <- function(formula, data, subset,
                            rts = c("C", "NI", "V"), rts.subst = NULL,
                            base = c("output", "input"), print.level = 1,
                            type = "RM", winw = 50, ...)
{
 # get y and x matrices

 # mf0 <- match.call(expand.dots = FALSE, call = sys.call(which = 1))
 needed.frame <- sys.nframe() - 1
 mf0 <- match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = needed.frame)))

 # if(print.level >= 1) print(mf0)

 # check if it is a matrix
 datasupplied <- !(match("data", names(mf0), 0) == 0)
 # subsetsupplied <- !(match("subset", names(mf0), 0) == 0)
 # print(subsetsupplied)
 # cat("1\n")

 # if data are supplied
 if(datasupplied){
  # begin get a logical vector equal TRUE if !missing

  # first using data and subset to get x without NA
  mf <- mf0
  m <- match(c("formula", "data", "subset"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  # mf$formula <- Formula( formula )
  # mf <- eval(mf, parent.frame(n = 1))
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  x <- as.matrix(model.matrix(mt, mf))
  # now get the names in the entire data
  esample <- seq_len( nrow(data) ) %in% as.numeric(rownames(x))
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing

  # get the data
  y <- as.matrix(model.part(Formula(formula), data = data[esample,], lhs = 1))
  x <- as.matrix(model.matrix(Formula(formula), data = data[esample,], rhs = 1)[,-1])

  # print(y)
  # print(x)
 }
 # if data are not supplied
 else {
  mf <- mf0
  subsetsupplied <- !(match("subset", names(mf), 0) == 0)
  if(subsetsupplied) stop("Subset with matrices is not allowed")
  m <- match(c("formula", "data", "subset"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  # mf <- eval(mf, parent.frame())
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  y <- as.matrix(model.response(mf))
  x <- as.matrix(model.matrix(mt, mf)[,-1])
  # get a logical vector equal TRUE if !missing
  with.na <- model.frame(mt, na.action = na.pass)
  esample <- rownames(with.na) %in% rownames(x)
  # print(table(esample))
  # print(y)
  # print(x)
 }

 if( !is.numeric(y) ) stop("Some of the outputs are not numeric.")
 if( !is.numeric(x) ) stop("Some of the inputs are not numeric.")
 
 if(type == "RM"){
  # check negative
  x.negative <- apply(x, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
  if(sum(x.negative) == nrow(x)){
   stop("Negative values in at least one of the inputs for each data point", call. = FALSE)
  }
  y.negative <- apply(y, MARGIN = 1, FUN = function(q) sum(q<0)) >= 1
  if(sum(y.negative) == nrow(y)){
   stop("Negative values in at least one of the outputs for each data point", call. = FALSE)
  }
  # deal with negative
  if(sum(x.negative) > 0){
   warning(paste0("There are negative values in inputs for ",sum(x.negative)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  if(sum(y.negative) > 0){
   warning(paste0("There are negative values in outputs for ",sum(y.negative)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  x <- x[!x.negative & !y.negative,,drop = FALSE]
  y <- y[!x.negative & !y.negative,,drop = FALSE]
 } else {
  # check nonpositive
  x.nonpositive <- apply(x, MARGIN = 1, FUN = function(q) sum(q<=0)) >= 1
  if(sum(x.nonpositive) == nrow(x)){
   stop("Nonpositive values in at least one of the inputs for each data point", call. = FALSE)
  }
  y.nonpositive <- apply(y, MARGIN = 1, FUN = function(q) sum(q<00)) >= 1
  if(sum(y.nonpositive) == nrow(y)){
   stop("Nonpositive values in at least one of the outputs for each data point", call. = FALSE)
  }
  # deal with negative
  if(sum(x.nonpositive) > 0){
   warning(paste0("There are nonpositive values in inputs for ",sum(x.nonpositive)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  if(sum(y.nonpositive) > 0){
   warning(paste0("There are nonpositive values in outputs for ",sum(y.nonpositive)," data points; these data points will not be considered"), call. = FALSE, immediate. = TRUE)
   cat("\n")
  }
  x <- x[!x.nonpositive & !y.nonpositive,,drop = FALSE]
  y <- y[!x.nonpositive & !y.nonpositive,,drop = FALSE]
 }

 rts <- rts[1]
 base <- base[1]

 rts1 <- tolower(substr(rts, 1,1 ))
 base1 <- tolower(substr(base, 1,1 ))

 if(is.null(rts.subst)){
  if(rts1 == "c" | rts == 1){
   myrts <- 3
   myrts1 <- "CRS"
  } else if (rts1 == "n" | rts == 2){
   myrts <- 2
   myrts1 <- "NIRS"
  } else if (rts1 == "v" | rts == 3){
   myrts <- 1
   myrts1 <- "VRS"
  } else {
   stop("invalid 'rts'; 'rts' must be 'CRS', 'NIRS', or 'VRS'")
  }
 }
 else {
  myrts1 <- rts.subst
 }

 if (base1 == "o" | base == 2){
  mybase <- 2
  mybase1 <- "output"
 } else if (base1 == "i" | base == 1){
  mybase <- 1
  mybase1 <- "input"
 } else {
  stop("invalid 'base'; 'base' must be 'output' or 'input'")
 }

 # for printing
 #  if(print.level >= 1){
 #   cat("\n Radial (Debreu-Farrell) ",mybase1,"-based measures of technical efficiency \n", sep = "")
 #   cat(" under assumption of CRS, NIRS, and VRS technology are computed for the\n", sep = "")
 #   cat(" following data:\n\n", sep = "")
 #   cat("  Number of data points (K) = ",nrow(y),"\n", sep = "")
 #   cat("  Number of outputs     (M) = ",ncol(y),"\n", sep = "")
 #   cat("  Number of inputs      (N) = ",ncol(x),"\n", sep = "")
 #  }

 if(print.level >= 1 & winw > 50){
  mymesage <- paste("\n",ifelse(type == "RM", "Nonradial (Russell)", "Radial (Debrue-Farrell)")," ",mybase1,"-based measures of technical efficiency under assumption of ",myrts1," technology are computed for the following data:\n", sep = "")
  cat("",unlist(strsplit(mymesage, " ")),"", sep = " ", fill = winw-10 )
  # 		cat("\n ",ifelse(type == "RM", "Nonradial (Russell)", "Radial Debrue-Farrell")," ",mybase1,"-based measures of technical efficiency \n", sep = "")
  # 		cat(" under assumption of ",myrts1," technology are computed for the \n", sep = "")
  # 		cat(" following data:\n\n", sep = "")
  cat("  Number of data points (K) = ",nrow(y),"\n", sep = "")
  cat("  Number of outputs     (M) = ",ncol(y),"\n", sep = "")
  cat("  Number of inputs      (N) = ",ncol(x),"\n", sep = "")

  mymesage <- paste("\nReference set is formed by ",nrow(y)," data point(s), for which measures of technical efficiency are computed", sep = "")
  cat("",unlist(strsplit(mymesage, " ")),"", sep = " ", fill = winw-10 )
 }

 #  if(print.level >= 1){
 #    cat("\n Reference set is formed by ",nrow(y)," data points,\n", sep = "")
 #    cat(" for which measures of technical efficiency are computed.\n\n", sep = "")
 #  }

 tymch <- list(y = y, x = x, esample = esample, mybase = mybase, base.string = mybase1)
 if(is.null(rts.subst)){
  tymch$myrts <- myrts
  tymch$rts.string <- myrts1
 }
 class(tymch) <- "npsf"
 return(tymch)
}

# .su <- function(x, mat.var.in.col = TRUE, digits = 4, probs = c(0.1, 0.25, 0.5, 0.75, 0.9), print = FALSE){
# 
#  xvec2 <- xvec1 <- FALSE
# 
#  if(is.matrix(x)){
#   if(min(dim(x)) == 1){
#    # xvec1 <- TRUE
#    if(which(dim(x) == 1) == 2){
#     mynames <- colnames(x)
#    } else {
#     mynames <- rownames(x)
#     x <- t(x)
#    }
#    # x <- as.vector(x)
#   } else {
#    if(!mat.var.in.col){
#     x <- t(x)
#    }
#    mynames <- colnames(x)
#   }
#   # print(mynames)
#   if(is.null(mynames)) mynames <- paste("Var", seq_len(ncol(x)), sep = "")
#   x <- as.data.frame(x)
#   # print(x)
#   # mynames <- colnames(x)
#  } # end if matrix
# 
#  if(is.vector(x)){
#   xvec2 <- TRUE
#   mynames <- deparse(substitute(x))
#   x <- data.frame(Var1 = x)
#  } # end if vector
# 
#  # cat("nymanes", sep ="")
#  # print(mynames)
# 
#  if(!is.vector(x) & !is.matrix(x) & !is.data.frame(x)){
#   stop("Provide vector, matrix, or data.frame")
#  } else {
#   t1 <- apply(x, 2, function(x) c(Obs = length(x), NAs = sum(is.na(x)), Mean = mean(x, na.rm = TRUE), StDev = sd(x, na.rm = TRUE), IQR = IQR(x, na.rm = TRUE), Min = min(x, na.rm = TRUE), quantile( x, probs = probs, na.rm = TRUE ), Max = max(x, na.rm = TRUE)))
#   # print(t1)
#   # print(mynames)
#   # print(class(t1))
#   # print(dim(t1))
#   if(xvec2 & !xvec1) colnames(t1) <- mynames
#   if(print){
#    tymch <- formatC(t(t1), digits = 4, format = "f", width = 4+1)
#    tymch <- gsub(".0000","",tymch, fixed = TRUE)
#    print(tymch, quote = FALSE)
#   }
#  }
#  class(t1) <- "npsf"
#  return(t(t1))
# }

.my.prettyNum <- function(xx2){
 n <- length(xx2)
 xx2.pn <- prettyNum(xx2)
 my.integers <- !(1:n %in%  grep(".", xx2.pn, fixed = TRUE))
 xx3 <- ifelse(abs(xx2) < 1e-04 & xx2 != 0, formatC(xx2, digits = 1, format = "e", width = 5), ifelse(abs(xx2) >= 1e-04 & abs(xx2) < 10, formatC(xx2, format = "f", digits = 4), ifelse(abs(xx2) >= 10 & abs(xx2) < 100, formatC(xx2, format = "f", digits = 3), ifelse(abs(xx2) >= 100 & abs(xx2) < 1000, formatC(xx2, format = "f", digits = 2), ifelse(abs(xx2) >= 1000, formatC(xx2, format = "e", digits = 1), formatC(xx2, format = "f", digits = 1))))))
 xx3[my.integers] <- xx2.pn[my.integers]
 xx3
}

.su <- function(x, mat.var.in.col = TRUE, digits = 4, probs = c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9), print = FALSE, names = NULL, ...){
 
 xvec2 <- xvec1 <- FALSE
 # print(class(x))
 
 if(is.list(x)){
  # cat("this is a list\n")
  mynames <- paste0("Var",1:length(x))
 }
 
 if(is.matrix(x)){
  # cat("this is a matrix\n")
  if(min(dim(x)) == 1){
   # xvec1 <- TRUE
   if(which(dim(x) == 1) == 2){
    mynames <- colnames(x)
    
   } else {
    mynames <- rownames(x)
    x <- t(x)
   }
   # x <- as.vector(x)
  } else {
   if(!mat.var.in.col){
    x <- t(x)
   }
   mynames <- colnames(x)
  }
  
  # print(mynames)
  if(is.null(mynames)) mynames <- paste("Var", seq_len(ncol(x)), sep = "")
  x <- as.data.frame(x)
  # print(x)
  # mynames <- colnames(x)
 } # end if matrix
 
 if(is.vector(x) & !is.list(x)){
  # cat("this is a vector\n")
  xvec2 <- TRUE
  mynames <- deparse(substitute(x))
  x <- data.frame(Var1 = x)
 } # end if vector
 
 # print(mynames)
 # print(names)
 
 if(!is.null(names)){
  if(length(mynames) == length(names)){
   mynames <- names
  }
 }
 
 # cat("nynames\n", sep ="")
 # print(mynames)
 if(is.list(x)){
  t1 <- sapply(x, function(x) c(Obs = length(x), NAs = sum(is.na(x)), Mean = mean(x, na.rm = TRUE), StDev = sd(x, na.rm = TRUE), IQR = IQR(x, na.rm = TRUE), Min = min(x, na.rm = TRUE), quantile( x, probs = probs, na.rm = TRUE ), Max = max(x, na.rm = TRUE)))
  # print(t1)
  # print(mynames)
  # print(class(t1))
  # print(is.integer(t1))
  # print( formatC(t1, format = "f") == formatC(t1, format = "fg") )
  # print(dim(t1))
  # if(xvec2 & !xvec1) colnames(t1) <- mynames
  colnames(t1) <- mynames
  # if(print){
  #  # t2 <- rbind(
  #  #   formatC(t1[c(1:2),,drop = FALSE], digits = digits, format = "fg"),
  #  #   formatC(t1[-c(1:2),,drop = FALSE], digits = digits, format = "f")
  #  # )
  #  print(data.frame(.my.prettyNum(t1)))
  # }
  if(print){
   t2 <- data.frame(.my.prettyNum(t1))
   colnames(t2) <- mynames
   print(t2)
  } 
 } else if(!is.vector(x) & !is.matrix(x) & !is.data.frame(x)){
  stop("Provide list, vector, matrix, or data.frame")
 } else {
  t1 <- apply(x, 2, function(x) c(Obs = length(x), NAs = sum(is.na(x)), Mean = mean(x, na.rm = TRUE), StDev = sd(x, na.rm = TRUE), IQR = IQR(x, na.rm = TRUE), Min = min(x, na.rm = TRUE), quantile( x, probs = probs, na.rm = TRUE ), Max = max(x, na.rm = TRUE)))
  # print(t1)
  # print(mynames)
  print(is.integer(t1))
  print(class(t1))
  # print(dim(t1))
  if(xvec2 & !xvec1) colnames(t1) <- mynames
  if(print) print(data.frame(.my.prettyNum(t1)))#t1 <- data.frame(formatC(t1, digits = digits, ...)); print(t1)
  if(print){
   t2 <- data.frame(.my.prettyNum(t1))
   colnames(t2) <- mynames
   print(t2)
  } 
 }
 tymch <- t(t1)
 class(tymch) <- "sfgtrehet"
 return(tymch)
}

.rownames.Intercept.change <- function(x){
 # x <- gsub("(Intercept)","intercept",x, fixed = TRUE)
 x
}

.teRad <- function(Y,X,M,N,K,
                   Yr,Xr,Kref,
                   rts,base,ifqh,
                   print.level=0){
 .C("radial",
    as.double(Y),
    as.double(X),
    as.integer(M),
    as.integer(N),
    as.integer(K),
    as.double(Yr),
    as.double(Xr),
    as.integer(Kref),
    as.integer(rts),
    as.integer(base),
    as.integer(ifqh),
    as.integer(print.level),
    te = double(K),
    ifsintensities = as.integer(0),
    intensities = double(Kref*K) ) $te
}

.teRad1 <- function(Y,X,M,N,K,
                   Yr,Xr,Kref,
                   rts,base,ifqh,
                   print.level=0){
 t1 <- .C("radial",
    Y = as.double(Y),
    X = as.double(X),
    as.integer(M),
    as.integer(N),
    as.integer(K),
    Yr = as.double(Yr),
    Xr = as.double(Xr),
    as.integer(Kref),
    as.integer(rts),
    as.integer(base),
    as.integer(ifqh),
    as.integer(print.level),
    te = double(K),
    ifsintensities = as.integer(1),
    intensities = double(Kref*K) )
 tymch <- list(te = t1$te, 
               Y = t1$Y,X = t1$X,
               Yr = t1$Yr,Xr = t1$Xr, 
               intensity = matrix(t1$intensities, nrow = K, byrow = TRUE))
 tymch$te <- ifelse(tymch$te < -998, NA, tymch$te)
 tymch$intensity[is.na(tymch$te),] <- NA
 return(tymch)
}

# .teNonrad <- function(Y,X,M,N,K,
#                       Yr,Xr,Kref,
#                       rts,base,ifqh,
#                       print.level=0, only.eff = TRUE){
#  .C("nonradial",
#     as.double(Y),
#     as.double(X),
#     as.integer(M),
#     as.integer(N),
#     as.integer(K),
#     as.double(Yr),
#     as.double(Xr),
#     as.integer(Kref),
#     as.integer(rts),
#     as.integer(base),
#     as.integer(ifqh),
#     as.integer(print.level),
#     te = double(K), 
#     te.all = as.double(1), 
#     as.integer(0) )$te
# }


.teNonrad2 <- function( Y,X,M,N,K,
                        Yr,Xr,Kref, 
                        lmdConstr = TRUE,
                        full.solution = TRUE,
                        rts, base,ifqh, 
                        print.level = 0 )
{
  ncol1 <- ifelse(base == 1, N, M)
  # if(base == 1){
  #   sol.full <- K*N
  # } else {
  #   sol.full <- K*M
  # }
  sol.full <- ncol1*K
  sol.z    <- Kref*K
  te <- .C("nonradial",
             as.double(Y),
             as.double(X),
             as.integer(M),
             as.integer(N),
             as.integer(K),
             as.double(Yr),
             as.double(Xr),
             as.integer(Kref),
             as.integer(rts),
             as.integer(base),
             as.integer(ifqh),
             as.integer(print.level),
             te    = double(K),
             teall = double(sol.full),
             zall  = double(sol.z),
             as.integer(full.solution),
             as.integer(lmdConstr))
  if(full.solution){
    rez <- list(te = te$te, 
                te.detail = t(matrix(te$teall, nrow = ncol1)),
                intensity  = matrix(te$zall, nrow = K, byrow = TRUE))
    rez$te <- ifelse(rez$te < -998, NA, rez$te)
    rez$te.detail[is.na(rez$te),] <- NA
    rez$intensity[is.na(rez$te),] <- NA
  } else {
    rez <- ifelse(te$te < -998, NA, te$te)
  }
  return(rez)
}

.dots <- function(nrep, message = NULL, width = 50, character = "."){
 if (!is.numeric(width)) {
  stop("'width' should be numeric")
 }
 if (width != 50 & width != 60 & width != 70 & width != 80 & width != 90  & width != 100){
  stop("'width' should be 50, 60, 70, 80, 90, or 100")
 }
 if (nrep == 0){
  if (!is.null(message)){
   cat("",message,"\n", sep = "")
  }
  cat("____|___ 1 ___|___ 2 ___|___ 3 ___|___ 4 ___|___ 5", sep = "")
  if (width == 60) {
   cat(" ___|___ 6", sep = "")
  }
  if (width == 70) {
   cat(" ___|___ 6 ___|___ 7", sep = "")
  }
  if (width == 80) {
   cat(" ___|___ 6 ___|___ 7 ___|___ 8", sep = "")
  }
  if (width == 90) {
   cat(" ___|___ 6 ___|___ 7 ___|___ 8 ___|___ 9", sep = "")
  }
  if (width == 100) {
   cat(" ___|___ 6 ___|___ 7 ___|___ 8 ___|___ 9 ___|___ 10", sep = "")
  }
  cat("\n")
 }
 else {
  if (nrep/width != floor(nrep/width)){
   cat("",character,"", sep = "")
  }
  else {
   cat("",character," ",nrep,"\n", sep = "")
  }
 }
}

.smplHom <- function(teRef,terfl,Kr,mybw,scVarHom,YorX,ba){
 # sample with replacement from Farrell measures
 # (including reflected ones)
 newSample <- sample( seq_len(2*Kr), Kr, replace = TRUE)
 bStar     <- terfl[newSample]
 # generate the sequence randomly
 epsStar   <- rnorm(Kr)
 tStar     <- bStar + mybw * epsStar
 tStar     <- ifelse(tStar >= 1, tStar, 2-tStar)
 # correct the sequence
 tStar     <- scVarHom * (tStar - mean(bStar)) + mean(bStar)
 # get new xobs if input-based, new yobs if output-based
 if (ba == 1){
  MatStar  <- YorX * ( teRef * tStar )
 }
 else {
  MatStar  <- YorX * ( teRef / tStar )
 }
 return(MatStar)
}

.smplHomTE <- function(terfl,Kr,mybw,scVarHom,ba){
 # sample with replacement from Farrell measures
 # (including reflected ones)
 newSample <- sample( seq_len(2*Kr), Kr, replace = TRUE)
 bStar     <- terfl[newSample]
 # generate the sequence randomly
 epsStar   <- rnorm(Kr)
 tStar     <- bStar + mybw * epsStar
 tStar     <- ifelse(tStar >= 1, tStar, 2-tStar)
 # correct the sequence
 tStar     <- scVarHom * (tStar - mean(bStar)) + mean(bStar)
 # inverse if input-based
 if (ba == 1){
  tStar    <- 1/tStar
 }
 return(tStar)
}

.smplHet1 <- function(Zt,L1,L2,M1,cons1,cons2,onesN,Kr,nZ,ba,M){
 ZStar1 <- NULL
 # print(class(Zt))
 if (ba == 1){
  ymax <- apply(Zt[, 1:M, drop = FALSE], MARGIN = 2, FUN = max)
  # print(ymax)
 } 	else {
  xmin <- apply(Zt[, M:(nZ-1), drop = FALSE], MARGIN = 2, FUN = min)
  # print(xmin)
 }
 cons3 <- cons1[1,] * cons2[1,]
 for(ww in seq_len(Kr)){
  flag <- TRUE
  while( flag ){
   # step 5
   newSample <- sample( seq_len(2*Kr), 1, replace = TRUE)
   ZStar <- Zt[newSample,]
   # step 6
   epsStar <- rnorm(nZ)
   if(newSample <= Kr){
    epsStar <- epsStar %*% L1
   } else {
    epsStar <- epsStar %*% L2
   }
   # step 7
   ZStar <- ZStar + as.vector(epsStar) * cons3
   if(ba == 1){
    flag0 <- max( ZStar[1:M] - ymax ) > 0
   }
   else {
    flag0 <- min( ZStar[M:(nZ-1)] - xmin ) < 0
   }
   flag <- min(ZStar) < 0 || flag0
  }
  ZStar1 <- rbind(ZStar1, ZStar)
 }

 ZStar <- ZStar1

 #  # toss possible negative values
 #  nonNegV <- apply(ZStar, MARGIN = 1, FUN = function(x) ifelse(min(x) < 0, FALSE, TRUE))
 #  # nonNegV <- rep(TRUE, Kr)
 #  ZStar <- ZStar[nonNegV,]
 nonNegV <- rep(TRUE, Kr)

 # step 8
 ZStar[,nZ] <- ifelse(ZStar[,nZ] > 1, ZStar[,nZ], 2 - ZStar[,nZ])

 # step 9
 if (ba == 1){
  teb <- 1 / ZStar[,nZ]
  yrb <- ZStar[, 1:M]
  xrb <- cbind( 1, tan(ZStar[,(M+1):(nZ-1)]) )
 } 	else {
  teb <- ZStar[,nZ]
  yrb <- cbind( 1, tan(ZStar[,1:(M-1)]) )
  xrb <- ZStar[, M:(nZ-1)]
 }
 return( list(Yrb = yrb, Xrb = xrb, teb = teb, Krb = sum(nonNegV)) )
}

.smplHet <- function(Zt,L1,L2,M1,cons1,cons2,onesN,Kr,nZ,ba,M,N){
 MinZ <- ifelse(M == 1, M, M - 1)
 if (ba == 1){
  ymax <- apply(Zt[, 1:M, drop = FALSE], MARGIN = 2, FUN = max)
 } 	else {
  xmin <- apply(Zt[, (MinZ+1):(nZ-1), drop = FALSE], MARGIN = 2, FUN = min)
 }
 # first try
 # step 5
 newSample <- sample( seq_len(2*Kr), Kr, replace = TRUE)
 ZStar <- Zt[newSample,]
 zBarStar <- matrix( colMeans( ZStar ), nrow  = 1)

 # step 6
 epsStar <- matrix(rnorm(Kr * nZ), nrow = Kr, ncol = nZ)
 for(ww in seq_len(Kr)){
  if(newSample[ww] <= Kr){
   epsStar[ww,] <- epsStar[ww, , drop = F] %*% L1
  } else {
   epsStar[ww,] <- epsStar[ww, , drop = F] %*% L2
  }
 }
 # step 7
 ZStar <- kronecker (onesN, zBarStar) + cons1 * ( M1 %*% ZStar + cons2 * epsStar)

 # toss values outside the frontier (including possible negative values)
 nonNeg <- apply(ZStar, MARGIN = 1, FUN = function(x) ifelse(min(x) < 0, FALSE, TRUE))
 if(ba == 1){
  withinFr <- apply(ZStar[,1:M, drop = FALSE], MARGIN = 1, FUN = function(z) min(ymax-z) > 0)
 }
 else {
  withinFr <- apply(ZStar[,(MinZ+1):(nZ-1), drop = FALSE], MARGIN = 1, FUN = function(z) max(xmin-z) < 0)
 }
 Zgood <- withinFr & nonNeg
 ZStar0 <- ZStar[Zgood,]

 # complete Zstar if not full
 if( nrow(ZStar0) < Kr){
  while (nrow(ZStar0) < Kr){
   # step 5
   newSample <- sample( seq_len(2*Kr), Kr, replace = TRUE)
   ZStar <- Zt[newSample,]
   zBarStar <- matrix( colMeans( ZStar ), nrow  = 1)
   # step 6
   epsStar <- matrix(rnorm(Kr * nZ), nrow = Kr, ncol = nZ)
   for(ww in seq_len(Kr)){
    if(newSample[ww] <= Kr){
     epsStar[ww,] <- epsStar[ww, , drop = F] %*% L1
    } else {
     epsStar[ww,] <- epsStar[ww, , drop = F] %*% L2
    }
   }
   # step 7
   ZStar <- kronecker (onesN, zBarStar) + cons1 * ( M1 %*% ZStar + cons2 * epsStar)
   # toss values outside the frontier (including possible negative values)
   nonNeg <- apply(ZStar, MARGIN = 1, FUN = function(x) ifelse(min(x) < 0, FALSE, TRUE))
   if(ba == 1){
    withinFr <- apply(ZStar[,1:M, drop = FALSE], MARGIN = 1, FUN = function(z) min(ymax-z) > 0)
   }
   else {
    withinFr <- apply(ZStar[,(MinZ+1):(nZ-1), drop = FALSE], MARGIN = 1, FUN = function(z) max(xmin-z) < 0)
   }
   Zgood <- withinFr & nonNeg
   ZStar0 <- rbind(ZStar0, ZStar[Zgood,])
  }
 }
 ZStar <- ZStar0[seq_len(Kr),]
 # print(nrow(ZStar))

 # step 8
 ZStar[,nZ] <- ifelse(ZStar[,nZ] > 1, ZStar[,nZ], 2 - ZStar[,nZ])

 # step 9
 if (ba == 1){
  teb <- 1 / ZStar[,nZ]
  if ( N == 1){
   xrb <- ZStar[, (M+1), drop = FALSE]
  }
  else {
   xrb <- cbind( 1, tan(ZStar[,(M+1):(nZ-1)]) )
  }
  yrb <- ZStar[, 1:M, drop = FALSE]
 }
 else {
  teb <- ZStar[,nZ]
  if ( M == 1){
   yrb <- ZStar[, 1, drop = FALSE]
  }
  else {
   yrb <- cbind( 1, tan(ZStar[,1:MinZ]) )
  }
  xrb <- ZStar[, (MinZ+1):(nZ-1), drop = FALSE]
 }
 return( list(Yrb = yrb, Xrb = xrb, teb = teb, Krb = nrow(ZStar)) )
}

.biasAndCI <- function(te, teboot, msub  = 0, K, Kr, M, N,
                       level, smoothed, forceLargerOne){

 if ( forceLargerOne ){
  teff   <- 1/te
  teboot <- 1/teboot
 }
 else {
  teff <- te
 }

 if (smoothed) {
  con1 <- 1
  con2 <- 1
 }
 else {
  con1 <- msub ^ ( 2 / (M + N + 1) )
  con2 <- Kr ^ (-2 / (M + N + 1) )
 }
 con3 <- con1 * con2

 # 	# bias-correction
 # 	tebm <- colMeans(teboot, na.rm = TRUE)
 # 	bias <- con3 * (tebm - teff)
 # 	tebc <- teff - bias
 # 	vari <- apply(teboot, 2, var, na.rm = TRUE)
 # 	BovV <- bias^2 / vari * 3
 #
 # 	# CI
 # 	teOverTEhat <- sweep(teboot, 2, FUN = "/", teff)
 # 	teOverTEhat <- (teOverTEhat - 1) * con1
 # 	quans <- apply(teOverTEhat, 2, quantile, probs = (100 + c(-level, level))/200, na.rm = T)
 # 	LL <- teff / (1 + quans[2] * con2)
 # 	UU <- teff / (1 + quans[1] * con2)
 #
 # 	# drop if inference is based on less than 100
 # 	reaB <- apply(teboot, 2, function(x) sum( !is.na(x) ) )
 #
 # 	bias <- ifelse(reaB < 100, NA, bias)
 # 	vari <- ifelse(reaB < 100, NA, vari)
 # 	tebc <- ifelse(reaB < 100, NA, tebc)
 # 	LL <- ifelse(reaB < 100, NA, LL)
 # 	UU <- ifelse(reaB < 100, NA, UU)

 bias <- vari <- reaB <- BovV <- tebc <- LL <- UU <- numeric(K)

 for (i in seq_len(K) ) {
  TEi <- teboot[,i]
  TEi <- TEi[!is.na(TEi)]
  reaB[i] <- length(TEi)
  if( reaB[i] < 100 ){
   bias[i] = NA
   vari[i] = NA
   BovV[i] = NA
   tebc[i] = NA
   LL[i]   = NA
   UU[i]   = NA
  }
  else {
   bias[i] = con3 * ( mean(TEi) - teff[i] )
   vari[i] = var(TEi)
   BovV[i] = (bias[i])^2 / vari[i] * 3
   tebc[i] = teff[i] - bias[i]
   teOverTEhat = (TEi / teff[i] - 1) * con1
   quans = quantile(teOverTEhat, probs = (100 + c(-level, level))/200, na.rm = TRUE)
   LL[i] = teff[i] / (1 + quans[2] * con2)
   UU[i] = teff[i] / (1 + quans[1] * con2)
  }
 }

 return(list(reaB=reaB, bias=bias, vari=vari, BovV=BovV, tebc=tebc, LL=LL, UU=UU))

}


.pvalsTestOne <- function(seCrs, seCrsMean, te1boot, te2boot, ba){
 seCrsB <- te1boot / te2boot
 seCrsMeanB <- rowMeans(te1boot, na.rm = TRUE) / rowMeans(te2boot, na.rm = TRUE)
 seCrsMeanB <- na.omit( seCrsMeanB )
 if(ba ==1){
  # global CRS
  pgCRS <- mean(seCrsMeanB <= seCrsMean)
  # local CRS
  plCRS <- rowMeans( apply(seCrsB, MARGIN = 1, FUN = function(z) z <= seCrs), na.rm = TRUE )
 }
 else {
  pgCRS <- mean(seCrsMeanB >= seCrsMean)
  # local CRS
  plCRS <- rowMeans( apply(seCrsB, MARGIN = 1, FUN = function(z) z >= seCrs), na.rm = TRUE )
 }

 # count those that are not NA
 nonNa <- apply(seCrsB, MARGIN = 2, FUN = function(z) length(na.omit(z)))
 plCRS <- ifelse(nonNa >= 100, plCRS, NA)

 return(list(pgCRS = pgCRS, plCRS = plCRS, nonNa = nonNa, reaB = length(seCrsMeanB)))
}

.pvalsTestTwo <- function(seNrs, seNrsMean, te1boot, te2boot, ba, performGlobal, s.inefficient){
 seNrsB <- te1boot / te2boot
 # cat("ncol(te1boot) = ",ncol(te1boot),"\n")
 # if(ncol(te1boot)==1) seNrsB <- matrix(seNrsB, ncol = 1)
 if ( performGlobal ){
  seNrsMeanB <- rowMeans(te1boot, na.rm = TRUE) / rowMeans(te2boot, na.rm = TRUE)
  seNrsMeanB <- na.omit( seNrsMeanB )

  if(ba ==1){
   # global NRS
   pgNRS <- mean(seNrsMeanB <= seNrsMean)
   # local NRS
   plNRS <- rowMeans( apply(seNrsB, MARGIN = 1, FUN = function(z) z <= seNrs), na.rm = TRUE )
  }
  else {
   # global NRS
   pgNRS <- mean(seNrsMeanB >= seNrsMean)
   # local NRS
   plNRS <- rowMeans( apply(seNrsB, MARGIN = 1, FUN = function(z) z >= seNrs), na.rm = TRUE )
  }
 }
 else {
  seNrs1 <- seNrs[s.inefficient]
  if(ba ==1){
   # local NRS
   if(ncol(te1boot) == 1){
    plNRS <- mean( apply(seNrsB, MARGIN = 1, FUN = function(z) z <= seNrs1), na.rm = TRUE )
   }
   else {
    plNRS <- rowMeans( apply(seNrsB, MARGIN = 1, FUN = function(z) z <= seNrs1), na.rm = TRUE )
   }
  }
  else {
   # local NRS
   if(ncol(te1boot) == 1){
    plNRS <- mean( apply(seNrsB, MARGIN = 1, FUN = function(z) z >= seNrs1), na.rm = TRUE )
   }
   else {
    plNRS <- rowMeans( apply(seNrsB, MARGIN = 1, FUN = function(z) z >= seNrs1), na.rm = TRUE )
   }
  }
 }
 # count those that are not NA
 nonNa <- apply(seNrsB, MARGIN = 2, FUN = function(z) length(na.omit(z)))
 plNRS <- ifelse(nonNa >= 100, plNRS, NA)

 # write into appropriate cells
 pineffdrs <- rep(NA, length(s.inefficient))

 if (performGlobal){
  pineffdrs <- ifelse(s.inefficient, plNRS, pineffdrs)
 }
 else {
  pineffdrs[s.inefficient] <- plNRS
 }
 if ( performGlobal ){
  tymch <- list(pgNRS = ifelse(performGlobal, pgNRS, NA), pineffdrs = pineffdrs, nonNa = nonNa, reaB = length(seNrsMeanB))
 }
 else {
  tymch <- list(pgNRS = ifelse(performGlobal, pgNRS, NA), pineffdrs = pineffdrs, nonNa = nonNa)
 }
 return(tymch)
}

# empirical distribution function
.edf <- function(Y,y){
 if(length(y) == 1){
  f1 <- Y<=y
  # f1 <- y-Y > 1e-6
 }
 else {
  f3 <- cbind(y, t(Y))
  f2 <- apply(f3, MARGIN = 1, FUN = function(z) z[-1] <= z[1])
  # f2 <- apply(f3, MARGIN = 1, FUN = function(z)  z[1]-z[-1]> 1e-6)
  f1 <- rowMeans(f2)==1
 }
 return(mean(f1))
}

# empirical joint distribution function
.ejdf <- function(Y, X, y, x){
 # 1
 if(length(y) == 1){
  f1 <- Y<=y
 }
 else {
  f3 <- cbind(y, t(Y))
  f2 <- apply(f3, MARGIN = 1, FUN = function(z) z[-1] <= z[1])
  f1 <- rowMeans(f2)==1
 }
 # 2
 if(length(x) == 1){
  g1 <- X<=x
 }
 else {
  g3 <- cbind(x, t(X))
  g2 <- apply(g3, MARGIN = 1, FUN = function(z) z[-1] <= z[1])
  g1 <- rowMeans(g2)==1
 }
 sum(f1 & g1) / length(f1)
}

# empirical joint distribution function
.ejdfedf1 <- function(Y, X, y, x){
 # 1
 if(length(y) == 1){
  f1 <- Y<=y
 }
 else {
  f3 <- cbind(y, t(Y))
  f2 <- apply(f3, MARGIN = 1, FUN = function(z) z[-1] <= z[1])
  f1 <- rowMeans(f2)==1
 }
 # 2
 if(length(x) == 1){
  g1 <- X<=x
 }
 else {
  g3 <- cbind(x, t(X))
  g2 <- apply(g3, MARGIN = 1, FUN = function(z) z[-1] <= z[1])
  g1 <- rowMeans(g2)==1
 }
 # mean(f1 & g1) - mean(f1)*mean(g1)
 mean(f1)
}

.ejdfedf <- function(Y, X, y, x){
 # 1
 f1 <- Y <= y
 # 2
 g2 <- X
 for(i in seq_len(ncol(X))){
  g2[,i] <- X[,i] <= x[i]
 }
 g1 <- ( rowMeans(g2) ==1 )
 mean(f1 & g1) - mean(f1)*mean(g1)
 # mean(f1)
}

.t4n <- function(w,d, print=FALSE){
 n <- length(d)
 tt <- numeric(n)
 for(i in 1:n){
  tt[i] <- .ejdfedf(Y = d, y = d[i], X = w, x = w[i,]) # - .edf(Y = d, y = d[i]) * .edf(Y = w, y = w[i,]) #
 }
 if(print) print(cbind(tt,d))
 return(sum(tt^2))
}

# begin parallel computing
.run.boots.hom.rts <- function(inp, ...){
 # Read the arguments
 args <- list(...)[[1]]
 Y <- args$Y
 X <- args$X
 ba <- args$ba
 teRef <- args$teRef
 terfl <- args$terfl
 mybw <- args$mybw
 scVarHom <- args$scVarHom

 M  <- ncol(Y)
 N  <- ncol(X)
 K  <- nrow(Y)
 # Prepare output structure
 te1boot <- rep(NA, K)
 te2boot <- rep(NA, K)

 for(ii in seq_len(K)){
  # begin homogeneous
  # print(1)
  if (ba == 2) {
   # step 1: sampling
   Yb <- .smplHom(teRef,terfl,Kr=K,mybw,scVarHom,Y,ba)
   # step 2: applying DEA
   # CRS or NiRS
   te1boot[ii] <- .teRad(Y[ii,],X[ii,],M,N,1,t(Yb),t(X),K,args$rtsHo,ba,
                         0,print.level=0)
   # VRS
   te2boot[ii] <- .teRad(Y[ii,],X[ii,],M,N,1,t(Yb),t(X),K,1,ba,
                         0,print.level=0)
  }
  else {
   # step 1: sampling
   Xb <- .smplHom(teRef,terfl,Kr=K,mybw,scVarHom,X,ba)
   # step 2: applying DEA
   # CRS or NiRS
   te1boot[ii] <- .teRad(Y[ii,],X[ii,],M,N,1,t(Y),t(Xb),K,args$rtsHo,ba,
                         0,print.level=0)
   # VRS
   te2boot[ii] <- .teRad(Y[ii,],X[ii,],M,N,1,t(Y),t(Xb),K,1,ba,
                         0,print.level=0)
  }
  # end homogeneous
  # end for ii
 }
 return (c(te1boot, te2boot, K, 0))
}

.run.boots.hom.bc <- function(inp, ...){
 # Read the arguments
 args <- list(...)[[1]]
 Y <- args$Y
 X <- args$X
 Yr <- args$Yr
 Xr <- args$Xr
 rt <- args$rt
 ba <- args$ba
 teRef <- args$teRef
 terfl <- args$terfl
 mybw <- args$mybw
 scVarHom <- args$scVarHom

 M  <- ncol(Y)
 N  <- ncol(X)
 K  <- nrow(Y)
 Kr <- nrow(Yr)

 # begin homogeneous
 # print(1)
 # teB <- numeric(K)
 if (ba == 2) {
  # step 1: sampling
  Yrb <- .smplHom(teRef,terfl,Kr,mybw,scVarHom,Yr,ba)
  # step 2: applying DEA
  teB <- .teRad(t(Y),t(X),M,N,K,t(Yrb),t(Xr),Kr,rt,ba,
                0,print.level=0)
 }
 else {
  # step 1: sampling
  Xrb <- .smplHom(teRef,terfl,Kr,mybw,scVarHom,Xr,ba)
  # step 2: applying DEA
  teB <- .teRad(t(Y),t(X),M,N,K,t(Yr),t(Xrb),Kr,rt,ba,
                0,print.level=0)
 }
 return (c(teB, K, 1))
}

.run.boots.het.rts <- function(inp, ...){
 # Read the arguments
 args <- list(...)[[1]]
 Y <- args$Y
 X <- args$X
 Yr <- args$Yr
 Xr <- args$Xr
 Zt <- args$Zt
 L1 <- args$L1
 L2 <- args$L2
 M1 <- args$M1
 cons1 <- args$cons1
 cons2 <- args$cons2
 onesN <- args$onesN
 ba <- args$ba

 M  <- ncol(Y)
 N  <- ncol(X)
 K  <- nrow(Y)
 Kr <- nrow(Yr)
 nZ <- ncol(Zt)
 # Prepare output structure
 te1boot <- rep(NA, K)
 te2boot <- rep(NA, K)

 for(ii in seq_len(K)){
  # begin heterogeneous
  # print(2)
  # step 1: sampling
  smplHet <- .smplHet(Zt,L1,L2,M1,cons1,cons2,onesN,Kr,nZ,ba,M,N)
  # step 2: get efficiency under H0
  teRefB  <- .teRad(t(smplHet$Yrb),t(smplHet$Xrb),M,N,smplHet$Krb,
                    t(Yr),t(Xr),Kr,rts=3,ba,0,print.level=0)
  # step 3: get efficient xRef or yRef
  if (ba == 1) {
   Yrb <- smplHet$Yrb
   Xrb <- smplHet$Xrb * teRefB / smplHet$teb
  }
  else {
   Yrb <- smplHet$Yrb * teRefB / smplHet$teb
   Xrb <- smplHet$Xrb
  }
  # handle infeasible
  teRefB.good <- !is.na(teRefB)
  # modify the existing matrices by tossing the missing values
  Yrb <- Yrb[teRefB.good, , drop = FALSE]
  Xrb <- Xrb[teRefB.good, , drop = FALSE]
  # new number of obs in reference
  KrB <- sum(teRefB.good)
  cat("my sample is ",KrB,"\n", sep = "")
  # step 4: get the bootstrapped TE
  # CRS or NiRS
  te1boot[ii] <- .teRad(Y[ii,],X[ii,],M,N,1,t(Yrb),t(Xrb),KrB,args$rtsHo,ba,
                        0,print.level=0)
  # VRS
  te2boot[ii] <- .teRad(Y[ii,],X[ii,],M,N,1,t(Yrb),t(Xrb),KrB,1,ba,
                        0,print.level=0)
  # end heterogeneous
  # end for ii
 }
 return (c(te1boot, te2boot, K, 0))
}

.run.boots.het.bc <- function(inp, ...){
 # Read the arguments
 args <- list(...)[[1]]
 Y <- args$Y
 X <- args$X
 Yr <- args$Yr
 Xr <- args$Xr
 rt <- args$rt
 Zt <- args$Zt
 L1 <- args$L1
 L2 <- args$L2
 M1 <- args$M1
 cons1 <- args$cons1
 cons2 <- args$cons2
 onesN <- args$onesN
 ba <- args$ba

 M  <- ncol(Y)
 N  <- ncol(X)
 K  <- nrow(Y)
 Kr <- nrow(Yr)
 nZ <- ncol(Zt)

 # begin heterogeneous
 # print(2)
 # step 1: sampling
 smplHet <- .smplHet(Zt,L1,L2,M1,cons1,cons2,onesN,Kr,nZ,ba,M,N)
 # step 2: get efficiency under H0
 teRefB  <- .teRad(t(smplHet$Yrb),t(smplHet$Xrb),M,N,smplHet$Krb,
                   t(Yr),t(Xr),Kr,rt,ba, 0,print.level=0)
 # step 3: get efficient xRef or yRef
 if (ba == 1) {
  Yrb <- smplHet$Yrb
  Xrb <- smplHet$Xrb * teRefB / smplHet$teb
 }
 else {
  Yrb <- smplHet$Yrb * teRefB / smplHet$teb
  Xrb <- smplHet$Xrb
 }
 # handle infeasible
 teRefB.good <- !is.na(teRefB)
 # modify the existing matrices by tossing the missing values
 Yrb <- Yrb[teRefB.good, , drop = FALSE]
 Xrb <- Xrb[teRefB.good, , drop = FALSE]
 # new number of obs in reference
 KrB <- sum(teRefB.good)
 # step 4: get the bootstrapped TE
 teB <- .teRad(t(Y),t(X),M,N,K,t(Yrb),t(Xrb),KrB,rt,ba,
               0,print.level=0)
 # end heterogeneous
 return (c(teB, KrB, 1))
}

.run.dots <- function(xs, nrep, args){
 width <- args$width
 character <- "."
 #  mylen <- length(xs[[nrep]])
 #  if (mylen > 0){
 #   if(xs[[nrep]][ mylen ] == 1){
 #    # boot for BC
 #    K <- mylen - 2
 #   }
 #   else {
 #    # boot for RTS test
 #    K <- (mylen - 2 )/2
 #   }
 #   # Pre-Last value in each output is 'KrB'
 #   KrB <- xs[[nrep]][ mylen-1 ]
 #   # cat("krB = ",KrB," K = ", K, " \n", sep = "")
 #   if (KrB/K < 0.80){
 #    character <- "?"
 #   }
 #   else if (KrB/K < 0.95){
 #    character <- "@"
 #   }
 #   else if (KrB/K < 1){
 #    character <- "x"
 #   }
 #  }

 #  if (!is.numeric(width)) {
 #   stop("'width' should be numeric")
 #  }
 #  if (width != 50 & width != 60 & width != 70 &
 #      width != 80 & width != 90  & width != 100){
 #   stop("'width' should be 50, 60, 70, 80, 90, or 100")
 #  }
 #  if (nrep == 0){
 #   if (!is.null(message)){
 #    cat("",message,"\n", sep = "")
 #   }
 #   cat("____|___ 1 ___|___ 2 ___|___ 3 ___|___ 4 ___|___ 5", sep = "")
 #   if (width == 60) {
 #    cat(" ___|___ 6", sep = "")
 #   }
 #   if (width == 70) {
 #    cat(" ___|___ 6 ___|___ 7", sep = "")
 #   }
 #   if (width == 80) {
 #    cat(" ___|___ 6 ___|___ 7 ___|___ 8", sep = "")
 #   }
 #   if (width == 90) {
 #    cat(" ___|___ 6 ___|___ 7 ___|___ 8 ___|___ 9", sep = "")
 #   }
 #   if (width == 100) {
 #    cat(" ___|___ 6 ___|___ 7 ___|___ 8 ___|___ 9 ___|___ 10", sep = "")
 #   }
 #   cat("\n")
 #  }
 #  else {
 if (nrep/width != floor(nrep/width)){
  cat("",character,"", sep = "")
 }
 else {
  cat("",character," ",nrep,"\n", sep = "")
 }
 # }
}

.run.boots.subs.bc <- function(inp, ...){
 # Read the arguments
 args <- list(...)[[1]]
 Y <- args$Y
 X <- args$X
 Yr <- args$Yr
 Xr <- args$Xr
 rt <- args$rt
 ba <- args$ba
 msub <- args$msub

 M  <- ncol(Y)
 N  <- ncol(X)
 K  <- nrow(Y)
 Kr <- nrow(Yr)

 # begin subsampling
 # print(3)
 # step 1: sub-sampling
 newSample <- sample( seq_len(Kr), msub, replace = TRUE)
 ystar <- Yr[newSample, , drop = FALSE]
 xstar <- Xr[newSample, , drop = FALSE]
 # step 2: applying DEA
 teB <- .teRad(t(Y),t(X),M,N,K,t(ystar),t(xstar),msub,rt,ba,
               0,print.level=0)
 mychar <- "."
 # end subsampling
 return (c(teB, K, 1))
}

# end parallel computing

.findInterval <- function(x, vec){
  inte <- NULL
  for(i in seq_len(length(x))){
    inte <- c(inte, which.min( ifelse(vec - x[i] < 0, +Inf, vec - x[i]) ))
  }
  return(inte)
}

.prepareYXZ.cs <- function(formula, ln.var.u.0i = NULL, ln.var.v.0i = NULL, mean.u.0i = NULL, data, subset, sysnframe = 1, ...) {
 # needed.frame <- sys.nframe() - 1
 needed.frame <- sys.nframe() - sysnframe
 # cat("sys.nframe() = ",sys.nframe(),"\n")
 # cat("needed.frame = ",needed.frame,"\n")
 mf0 <- match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = needed.frame)))
 if ( is.null(ln.var.u.0i) ){
  ln.var.u.0i <- ~ 1
  ku <- 1
 } else {
  ku <- 17
 }
 if ( is.null(ln.var.v.0i) ){
  ln.var.v.0i <- ~ 1
  kv <- 1
 } else {
  kv <- 17
 }
 if ( is.null(mean.u.0i) ){
  mean.u.0i <- ~ 1
  kdel <- 1
 } else {
  kdel <- 17
 }
 
 form1 <- Formula(as.formula(paste("",deparse(formula, width.cutoff = 500L)," | ",ln.var.u.0i[2]," | ",ln.var.v.0i[2]," | ",mean.u.0i[2],"", sep = "")))
 
 form <- Formula(as.formula(paste("",deparse(formula, width.cutoff = 500L)," + ",ln.var.u.0i[2]," + ",ln.var.v.0i[2]," + ",mean.u.0i[2],"", sep = "")))
 
 # check if it is a matrix
 datasupplied <- !(match("data", names(mf0), 0) == 0)
 
 
 if(datasupplied){
  # begin get a logical vector equal TRUE if !missing
  # first using data and subset to get x without NA
  mf <- mf0
  mf$formula <- formula( form )
  m <- match(c("formula", "data", "subset"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  X <- as.matrix(model.matrix(mt, mf))
  # now get the names in the entire data
  # esample <- seq_len( nrow(data) ) %in% as.numeric(rownames(X))
  esample <- rownames(data) %in% rownames(X)
  if(sum(esample)==0){
   stop("no valid data points")
  }
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing
  
  # get the data
  # print(form1)
  dataesample <- model.frame(form1, data = data[esample,])
  # print(dataesample)
  Y <- as.matrix(model.part(form1, data = dataesample, lhs = 1, drop = FALSE))
  # Y <- as.matrix( model.matrix(formula(form1, lhs = 1, rhs = 0), data = data[esample,]))
  # print(Y)
  X <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 1), data = dataesample))
  # print(X)
  n <- nrow(Y)
  k <- ncol(X)
  if(ku == 1){
   Zu <- matrix(1, nrow = sum(esample), ncol = 1)
   colnames(Zu) <- "(Intercept)"
  }
  else {
   Zu <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 2), data = dataesample))
   ku <- ncol(Zu)
  }
  if(kv == 1){
   Zv <- matrix(1, nrow = sum(esample), ncol = 1)
   colnames(Zv) <- "(Intercept)"
  }
  else {
   Zv <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 3), data = dataesample))
   kv <- ncol(Zv)
  }
  if(kdel == 1){
   Zdel <- matrix(1, nrow = sum(esample), ncol = 1)
   colnames(Zdel) <- "(Intercept)"
  }
  else {
   Zdel <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 4), data = dataesample))
   kdel <- ncol(Zdel)
  }
  #   print(head(Y))
  #   print(head(X))
  #   print(head(Zu))
  #   print(head(Zv))
  #   print(head(Zdel))
 }
 # if data are not supplied
 else {
  # begin get a logical vector equal TRUE if !missing
  
  # first using data and subset to get XZ without NA
  mf <- mf0
  mf$formula <- formula( form )
  subsetsupplied <- !(match("subset", names(mf), 0) == 0)
  if(subsetsupplied) stop("Subset with matrices is not allowed")
  m <- match(c("formula"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  # mf <- eval(mf, parent.frame())
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  Y <- as.matrix(model.response(mf))
  n <- nrow(Y)
  XZ <- as.matrix(model.matrix(mt, mf))
  # get a logical vector equal TRUE if !missing
  with.na <- model.frame(mt, na.action = na.pass)
  esample <- rownames(with.na) %in% rownames(XZ)
  if(sum(esample)==0){
   stop("no valid data points")
  }
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing
  
  # get the data
  #   print(head(Y))
  #   print(head(XZ))
  
  # get X
  mf <- mf0
  m <- match(c("formula"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  X <- as.matrix(model.matrix(mt, mf))
  # print(head(X))
  k <- ncol(X)
  
  # get Zu
  if(ku == 1){
   Zu <- matrix(1, nrow = sum(esample), ncol = 1)
   colnames(Zu) <- "(Intercept)"
  }
  else {
   mf <- mf0
   mf$formula <- formula( Formula(as.formula(paste("",deparse(formula)," + ",ln.var.u.0i[2],"", sep = ""))))
   m <- match(c("formula"), names(mf), 0L)
   mf <- mf[c(1L, m)]
   mf[[1L]] <- as.name("model.frame")
   mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
   mt <- attr(mf, "terms")
   XZu <- as.matrix(model.matrix(mt, mf))
   Zu <- XZu[, -(2:k), drop = FALSE]
   ku <- ncol(Zu)
   # print(head(Zu))
  }
  
  # get Zv
  if(kv == 1){
   Zv <- matrix(1, nrow = sum(esample), ncol = 1)
   colnames(Zv) <- "(Intercept)"
  }
  else {
   mf <- mf0
   mf$formula <- formula( Formula(as.formula(paste("",deparse(formula)," + ",ln.var.u.0i[2]," + ",ln.var.v.0i[2],"", sep = ""))))
   m <- match(c("formula"), names(mf), 0L)
   mf <- mf[c(1L, m)]
   mf[[1L]] <- as.name("model.frame")
   mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
   mt <- attr(mf, "terms")
   XZuZv <- as.matrix(model.matrix(mt, mf))
   Zv <- XZuZv[, -(2:(k+ku-1)), drop = FALSE]
   kv <- ncol(Zv)
   # print(head(Zv))
  }
  
  # get Zdel
  if(kdel == 1){
   Zdel <- matrix(1, nrow = sum(esample), ncol = 1)
   colnames(Zdel) <- "(Intercept)"
  }
  else {
   Zdel <- XZ[, -(2:(k+ku-1+kv-1)), drop = FALSE]
   # print(head(Zdel))
  }
  
  # print(head(Zu))
  # print(head(Zv))
  # print(head(Zdel))
 }
 # colnames(X)[1] <- colnames(Zu)[1] <- colnames(Zv)[1] <- colnames(Zdel)[1] <- "Intercept"
 
 # if(sum(colnames(X) == "(Intercept)") > 0){
 #  # print(which(colnames(X) == "(Intercept)"))
 #  colnames(X)[which(colnames(X) == "(Intercept)")] <- "Intercept"
 # }
 # if(sum(colnames(Zu) == "(Intercept)") > 0){
 #  # print(which(colnames(Zu) == "(Intercept)"))
 #  colnames(Zu)[which(colnames(Zu) == "(Intercept)")] <- "Intercept"
 # }
 # if(sum(colnames(Zv) == "(Intercept)") > 0){
 #  # print(which(colnames(Zv) == "(Intercept)"))
 #  colnames(Zv)[which(colnames(Zv) == "(Intercept)")] <- "Intercept"
 # }
 # if(sum(colnames(Zdel) == "(Intercept)") > 0){
 #  # print(which(colnames(Zv) == "(Intercept)"))
 #  colnames(Zdel)[which(colnames(Zdel) == "(Intercept)")] <- "Intercept"
 # }
 
 colnames(X) <- .rownames.Intercept.change(colnames(X))
 colnames(Zu) <- .rownames.Intercept.change(colnames(Zu))
 colnames(Zv) <- .rownames.Intercept.change(colnames(Zv))
 colnames(Zdel) <- .rownames.Intercept.change(colnames(Zdel))
 
 # print(colnames(Zdel))
 
 tymch <- list(Y = Y, X = X, Zu = Zu, Zv = Zv, Zdel = Zdel, n = n, k = k, ku = ku, kv = kv, kdel = kdel, esample = esample)
 class(tymch) <- "npsf"
 return(tymch)
 
}

# Half-normal model

# Log-likelihood
.ll.hn <- function(theta, prod, k, kv, ku, kdel = NULL, y, Zv, Zu, Zdel = NULL, X) {
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[-c(1:k, (k+1):(k+kv))]
 e <- y-X%*%beta
 sig <- sqrt(exp(Zv%*%gv) + exp(Zu%*%gu))
 lmd <- sqrt(exp(Zu%*%gu)/exp(Zv%*%gv))

 # Log-likelihood
 llf <- sum(log(2) - log(sqrt(2*pi)) - log(sig) + pnorm(s*e*lmd/sig, log.p = TRUE) - 0.5*(e/sig)^2)
 return(llf)
}

# Gradient
.g.hn <- function(theta, prod, k, kv, ku, kdel = NULL, y, Zv, Zu, Zdel = NULL, X) {
 if(prod == TRUE){ s = -1 } else {s = 1}
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[-c(1:k, (k+1):(k+kv))]
 e <- y - X%*%beta
 expu <- exp(Zu%*%gu); expv <- exp(Zv%*%gv)
 sig <- sqrt(expv + expu); lmd <- sqrt(expu/expv)
 ls <- lmd/sig; g <- dnorm(s*e*ls)/pnorm(s*e*ls)
 gels <- g*e*ls;

 # Gradient
 gb <- t(X)%*%(-s*g*ls + e/sig^2)
 ggv <- 0.5*t(Zv)%*%(expv/sig^2 * ((e/sig)^2 - 1) - s*gels*(1 + expv/sig^2))
 ggu <- 0.5*t(Zu)%*%(expu/sig^2 * ((e/sig)^2 - 1) - s*gels*(expu/sig^2 - 1))
 grad <- rbind(gb, ggv, ggu)
 return(grad)
}

# Hessian
.hess.hn <- function(theta, prod, k, kv, ku, kdel = NULL, y, Zv, Zu, Zdel = NULL, X) {
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[-c(1:k, (k+1):(k+kv))]
 e <- y - X%*%beta
 expu <- exp(Zu%*%gu); expv <- exp(Zv%*%gv)
 sig <- sqrt(expv + expu); lmd <- sqrt(expu/expv)
 ls <- lmd/sig; els <- e*ls;
 g <- dnorm(s*e*ls)/pnorm(s*e*ls); gels <- g*els

 # Hessian
 Hb <- t(as.numeric(-s*g*(els + s*g)*ls^2 - sig^(-2))*X)%*%X
 Hbgu <- t(as.numeric(0.5*g*ls*(1 - expu/sig^2)*((g + s*els)*els - s*1) - e*expu/sig^4)*X)%*%Zu
 Hbgv <- t(as.numeric(-s*0.5*g*ls*(1 + expv/sig^2)*((els + s*g)*els - 1) - e*expv/sig^4)*X)%*%Zv
 Hgvu <- t(0.5*as.numeric(-s*0.5*gels*(1 - expu/sig^2)*(1 + expv/sig^2)*(1 - s*(g + s*els)*els) - expv*expu*(2*(e/sig)^2 - s*gels - 1)/sig^4)*Zv)%*%Zu
 Hgu <- t(0.5*as.numeric((1 - expu/sig^2)*(expu/sig^2 * ((e/sig)^2 - s*gels - 1) + s*0.5*gels*(1 - expu/sig^2)*(1 - s*(g + s*els)*els)) - (expu*e/sig^3)^2)*Zu)%*%Zu
 Hgv <- t(0.5*as.numeric(expv/sig^2 * ((1 - expv/sig^2)*((e/sig)^2 - s*gels - 1) - e^2*expv/sig^4) - s*0.5*gels*(1 + expv/sig^2)^2 * ((els + s*g)*els - 1))*Zv)%*%Zv
 H <- cbind(rbind(Hb, t(Hbgv), t(Hbgu)), rbind(Hbgv, Hgv, t(Hgvu)), rbind(Hbgu, Hgvu, Hgu))
 return(H)
}

# Truncated-normal model

#Log-likelihood
.ll.tn <- function(theta, prod, k, kv, ku, kdel, y, Zv, Zu, X, Zdel) {
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[(k+kv+1):(k+kv+ku)]
 delta <- theta[(k+kv+ku+1):(k+kv+ku+kdel)]
 e <- y-X%*%beta
 sig <- sqrt(exp(Zv%*%gv) + exp(Zu%*%gu))
 lmd <- sqrt(exp(Zu%*%gu)/exp(Zv%*%gv))
 mu <- Zdel%*%delta

 # Log-likelihood
 llf <- sum(-log(sqrt(2*pi)) - log(sig) - pnorm(mu/sqrt(exp(Zu%*%gu)), log.p = TRUE) + pnorm(mu/(sig*lmd) + s*e*lmd/sig, log.p = TRUE) - 0.5*(((e - s*mu)^2)/sig^2))
 return(llf)
}

#Gradient

.g.tn <- function(theta, prod, k, kv, ku, kdel, y, Zv, Zu, X, Zdel) {
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[(k+kv+1):(k+kv+ku)]
 delta <- theta[(k+kv+ku+1):(k+kv+ku+kdel)]
 e <- y-X%*%beta
 expu <- exp(Zu%*%gu); expv <- exp(Zv%*%gv)
 sig <- sqrt(exp(Zv%*%gv) + exp(Zu%*%gu))
 lmd <- sqrt(exp(Zu%*%gu)/exp(Zv%*%gv))
 mu <- Zdel%*%delta; ls <- lmd/sig
 g1 <- dnorm(ls*(mu/lmd^2 + s*e))/pnorm(ls*(mu/lmd^2 + s*e))
 g2 <- dnorm(mu/sqrt(expu))/pnorm(mu/sqrt(expu)); ml = mu/lmd^2
 d1 <- 0.5*mu/(sig*lmd)*(1 - expv/sig^2)
 d2 <- -0.5*s*e*ls * (1 + expv/sig^2)
 d3 = -0.5*mu/(sig*lmd)*(1 + expu/sig^2)
 d4 = 0.5*s*e*ls * (1 - expu/sig^2)

 # Gradient
 gb <- t(X)%*%((e - s*mu)/sig^2 - s*g1*ls)
 gdel <- t(Zdel)%*%(-g2/sqrt(expu) + g1/(sig*lmd) + s*(e - s*mu)/sig^2)
 gv <- t(Zv)%*%(-0.5*expv/sig^2 + g1*(d1 + d2) + 0.5*expv*(e - s*mu)^2/sig^4)
 gu  <- t(Zu)%*%(-0.5*expu/sig^2 + 0.5*g2*mu/sqrt(expu) + g1*(d3 + d4) + 0.5*expu*(e - s*mu)^2/sig^4)
 grad <- rbind(gb, gv, gu, gdel)
 return(grad)
}

#Hessian
.hess.tn <- function(theta, prod, k, kv, ku, kdel, y, Zv, Zu, X, Zdel) {
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[(k+kv+1):(k+kv+ku)]
 delta <- theta[(k+kv+ku+1):(k+kv+ku+kdel)]
 e <- y-X%*%beta
 expu <- exp(Zu%*%gu); expv <- exp(Zv%*%gv)
 sig <- sqrt(exp(Zv%*%gv) + exp(Zu%*%gu))
 lmd <- sqrt(exp(Zu%*%gu)/exp(Zv%*%gv))
 mu <- Zdel%*%delta; ls <- lmd/sig
 g1 <- dnorm(ls*(mu/lmd^2 + s*e))/pnorm(ls*(mu/lmd^2 + s*e))
 gr1u <- ls*(0.5*(1 - expu/sig^2)*(mu/lmd^2 + s*e) - mu/lmd^2)*g1*(-s*ls*(e + s*mu/lmd^2) - g1)
 gr1v <- ls*(mu/lmd^2 - 0.5*(1 + expv/sig^2)*(mu/lmd^2 + s*e))*g1*(-s*ls*(e + s*mu/lmd^2) - g1)
 g2 <- dnorm(mu/sqrt(expu))/pnorm(mu/sqrt(expu)); ml = mu/lmd^2

 # Hessian
 Hb <- t(as.numeric(-1/sig^2 - s*g1*ls^2*(ls*(e + s*ml) + s*g1))*X)%*%X
 Hbdel <- t(as.numeric(-s/sig^2 - s*g1*(ls/lmd)^2*(-s*ls*(e + s*ml) - g1))*X)%*%Zdel
 Hbgv <- t(as.numeric(-(e - s*mu)*expv/sig^4 + s*0.5*ls*(1 + expv/sig^2)*g1 - s*ls*gr1v)*X)%*%Zv
 Hbgu <- t(as.numeric(-(e - s*mu)*expu/sig^4 - s*0.5*ls*(1 - expu/sig^2)*g1 - s*ls*gr1u)*X)%*%Zu
 Hdel <- t(as.numeric(g2/expu*(mu/sqrt(expu) + g2) - g1*ls/(lmd^3*sig) * (ls*(ml + s*e) + g1) - sig^(-2))*Zdel)%*%Zdel
 Hdelgv <- t(as.numeric(1/(lmd*sig) * (0.5*(1 - expv/sig^2)*g1 + gr1v) - s*(e - s*mu)*expv/sig^4)*Zdel)%*%Zv
 Hdelgu <- t(as.numeric(-s*(e - s*mu)*expu/sig^4 + 1/(sig*lmd)*(gr1u - 0.5*g1*(1 + expu/sig^2)) - 0.5*g2/sqrt(expu)*(-1 + mu/sqrt(expu)*(mu/sqrt(expu) + g2)))*Zdel)%*%Zu
 Hgv <- t(as.numeric(0.5*expv/sig^2 * ((1 - expv/sig^2)*((e - s*mu)^2/sig^2 - 1) - (e - s*mu)^2*expv/sig^4) + (ml - 0.5*(1 + expv/sig^2)*(ml + s*e))*(gr1v*ls - 0.5*g1*ls*(1 + expv/sig^2)) + g1*ls*(ml - 0.5*(expv/sig^2 * (1 - expv/sig^2)*(ml + s*e) + (1 + expv/sig^2)*ml)))*Zv)%*%Zv
 Hgu <- t(as.numeric((0.5*(1 - expu/sig^2)*(ml + s*e) - ml)*ls*(gr1u + 0.5*g1*(1 - expu/sig^2)) + g1*ls*(0.5*(expu/sig^2 - 1)*(expu/sig^2*(ml + s*e) + ml) + ml) + 0.5*(expu/sig^2 * ((1 - expu/sig^2)*((e - s*mu)^2/sig^2 - 1) - expu*(e - s*mu)^2/sig^4) + 0.5*g2*mu/sqrt(expu)*(-1 + mu/sqrt(expu)*(mu/sqrt(expu) + g2))))*Zu)%*%Zu
 Hgvgu <- t(as.numeric(0.5*expv*expu/sig^4 * (1 - 2*((e - s*mu)/sig)^2) + (ml - 0.5*(1 + expv/sig^2)*(ml + s*e))*ls*(gr1u + 0.5*(1 - expu/sig^2)*g1) + g1*ls*(-ml + 0.5*(expv*expu/sig^4*(ml + s*e) + (1 + expv/sig^2)*ml)))*Zv)%*%Zu

 H <- cbind(rbind(Hb, t(Hbgv),  t(Hbgu),t(Hbdel)), rbind(Hbgv, Hgv, t(Hgvgu), Hdelgv),rbind(Hbgu, Hgvgu, Hgu, Hdelgu), rbind(Hbdel, t(Hdelgv), t(Hdelgu), Hdel ))
 return(H)
}

# :: N-Exp ----------------------------------------------------------------

.ll.en <- function(theta, prod, y, X, Zv, Zu, n.ids, k, kv, ku, ...){
  myprod <- ifelse(prod, 1, -1)
  # k:    betas
  # kv:  noise, scale/variance
  # ku:  inefficiency: scale/variance
  
  # eps0  <- y - X %*%    theta[1                    : k]
  # # cat.print(2)
  # sv    <- sqrt( exp( Zv %*% theta[(k + 1)        : (k + kv)]) )
  # # cat.print(sv2)
  # su    <- sqrt( exp( Zu %*% theta[(k + kv  + 1) : (k + kv + ku)] ) )
  # # cat.print(su2)
  # 
  # ll    <- -log(su) + myprod*eps0/su + sv^2/2/su^2 + 
  #   pnorm(-sv/su - myprod*eps0/sv, log.p = TRUE)
  e  <- y - X %*%    theta[1                    : k]
  # cat.print(2)
  lnsigv2    <- Zv %*% theta[(k + 1)        : (k + kv)]
  # cat.print(sv2)
  lnsigu2    <- Zu %*% theta[(k + kv  + 1) : (k + kv + ku)]
  # cat.print(su2)
  
  z  <- -myprod * e * exp(-lnsigv2/2) - exp((lnsigv2-lnsigu2)/2)
  
  ll <- -1/2 * lnsigu2 + 0.5 * exp(lnsigv2-lnsigu2) + pnorm(z, log.p = TRUE) + myprod * e *exp(-lnsigu2/2) 
  
  return( sum(ll, na.rm = FALSE) )
}

# grad --------------------------------------------------------------------

.g.en <- function(theta, prod, y, X, Zv, Zu, n.ids, k, kv, ku, ...) {
  myprod <- ifelse(prod, 1, -1)
  e  <- y - X %*%    theta[1                    : k]
  # cat.print(2)
  lnsigv2    <- Zv %*% theta[(k + 1)        : (k + kv)]
  # cat.print(sv2)
  lnsigu2    <- Zu %*% theta[(k + kv  + 1) : (k + kv + ku)]
  # cat.print(su2)
  
  z  <- -myprod * e * exp(-lnsigv2/2) - exp((lnsigv2-lnsigu2)/2)
  
  zr = dnorm(z)/pnorm(z)
  # cat.print(length(zr))
  ztd1 = myprod*exp(-lnsigv2/2)
  # cat.print(length(ztd1))
  ztd2 = 1/2*myprod*e*exp(-lnsigv2/2)- 1/2*exp((lnsigv2-lnsigu2)/2)
  # cat.print(length(ztd2))
  ztd3 = 0.5*exp((lnsigv2-lnsigu2)/2)
  # cat.print(length(ztd3))
  g1 = zr * ztd1 - myprod*exp(-lnsigu2/2)
  # cat.print(dim(g1))
  g2 = 1/2*exp(lnsigv2-lnsigu2) + zr*ztd2
  g3 = -1/2-1/2*exp(lnsigv2-lnsigu2) + zr*ztd3-1/2*myprod*e*exp(-lnsigu2/2)
  # cat.print(dim(x))
  g  = cbind(X*as.vector(g1), Zv*as.vector(g2), Zu*as.vector(g3))
  return(colSums(g))
}

# hess --------------------------------------------------------------------

.hess.en <- function(theta, prod, y, X, Zv, Zu, n.ids, k, kv, ku, ...) {
  myprod <- ifelse(prod, 1, -1)
  e  <- y - X %*%    theta[1                    : k]
  # cat.print(2)
  lnsigv2    <- Zv %*% theta[(k + 1)        : (k + kv)]
  # cat.print(sv2)
  lnsigu2    <- Zu %*% theta[(k + kv  + 1) : (k + kv + ku)]
  # cat.print(su2)
  
  z  <- -myprod * e * exp(-lnsigv2/2) - exp((lnsigv2-lnsigu2)/2)
  
  zr = dnorm(z)/pnorm(z)
  # cat.print(length(zr))
  ztd1 = myprod*exp(-lnsigv2/2)
  # cat.print(length(ztd1))
  ztd2 = 1/2*myprod*e*exp(-lnsigv2/2)- 1/2*exp((lnsigv2-lnsigu2)/2)
  # cat.print(length(ztd2))
  ztd3 = 0.5*exp((lnsigv2-lnsigu2)/2)
  
  ztd = zr*(zr+z)
  d11 = ztd*ztd1*ztd1#, eq(1)
  # cat.print(dim(d11))
  d12 = ztd*ztd2*ztd1  - zr*(-1/2*myprod*exp(-lnsigv2/2))#, eq(1,2)
  d13 = ztd*ztd3*ztd1  - 1/2*myprod*exp(-lnsigu2/2)#, eq(1,3)
  d22 = ztd*ztd2*ztd2  - zr*(-1/4*myprod*e*exp(-lnsigv2/2) - 1/4*exp((lnsigv2-lnsigu2)/2))  - 1/2*exp(lnsigv2-lnsigu2)#, eq(2)
  d23 = ztd*ztd2*ztd3 - zr*(1/4*exp((lnsigv2-lnsigu2)/2))  + 1/2*exp(lnsigv2-lnsigu2)#, eq(2,3)
  d33 = ztd*ztd3*ztd3  - zr*(-1/4*exp((lnsigv2-lnsigu2)/2)) - 1/2*exp(lnsigv2-lnsigu2) - 1/4*myprod*e*exp(-lnsigu2/2)#, eq(3)
  
  d12 <- as.vector(d12)
  d13 <- as.vector(d13)
  
  h0 <- matrix(0,k + kv + ku,k + kv + ku)
  for(q in 1:length(y)){
    # cat.print(tcrossprod(X[q,])*d11[q])
    # cat.print(X[q,]*d12[q])
    # cat.print(X[q,]*d13[q])
    # cat.print(c( X[q,]*d12[q], d22[q], d23[q] ))
    # cat.print(c( X[q,]*d13[q], d22[q], d33[q] ))
    # cat.print(cbind( tcrossprod(X[q,])*d11[q], X[q,]*d12[q], X[q,]*d13[q] ))
    # hoq <- rbind()
    # cat.print(crossprod(X[q,,drop=FALSE],Zv[q,,drop=FALSE])*d12[q])
    # cat.print(crossprod(X[q,,drop=FALSE],Zu[q,,drop=FALSE])*d13[q])
    dd11 <- tcrossprod(X[q,])*d11[q]
    dd12 <- crossprod(X[q,,drop=FALSE],Zv[q,,drop=FALSE])*d12[q]
    dd13 <- crossprod(X[q,,drop=FALSE],Zu[q,,drop=FALSE])*d13[q]
    dd22 <- tcrossprod(Zv[q,])*d22[q]
    dd33 <- tcrossprod(Zu[q,])*d33[q]
    dd23 <- crossprod(Zv[q,,drop=FALSE],Zu[q,,drop=FALSE])*d23[q]
    # cat.print(dd11)
    # cat.print(dd12)
    # cat.print(dd13)
    # cat.print(dd22)
    # cat.print(dd23)
    # cat.print(dd33)
    h0 <- h0 - 
      rbind( 
        cbind(dd11, dd12, dd13),
        cbind(t(dd12), dd22, dd23),
        cbind(t(dd13), t(dd23), dd33)
      )
  }
  # cat.print(h0)
  return(h0)
  
  # negH = (d11, d12, d13 \ d12, d22, d23 \ d13, d23, d33)
  
  # h0 <- hessian1(func = .ll.n.exp, at = theta, myprod=myprod, y=y, x=x, Zv=Zv, Zu=Zu, n.ids=n.ids, k=k, kv=kv, ku=ku)
  # return(ifelse(is.na(h0) | is.infinite(h0), 0, h0) )
  
}

# grad + hess -------------------------------------------------------------

.gr.hess.en <- function(theta, prod, y, X, Zv, Zu, n.ids, k, kv, ku, ...) {
  myprod <- ifelse(prod, 1, -1)
  e  <- y - X %*%    theta[1                    : k]
  # cat.print(2)
  lnsigv2    <- Zv %*% theta[(k + 1)        : (k + kv)]
  # cat.print(sv2)
  lnsigu2    <- Zu %*% theta[(k + kv  + 1) : (k + kv + ku)]
  # cat.print(su2)
  
  z  <- -myprod * e * exp(-lnsigv2/2) - exp((lnsigv2-lnsigu2)/2)
  
  zr = dnorm(z)/pnorm(z)
  # cat.print(length(zr))
  ztd1 = myprod*exp(-lnsigv2/2)
  # cat.print(length(ztd1))
  ztd2 = 1/2*myprod*e*exp(-lnsigv2/2)- 1/2*exp((lnsigv2-lnsigu2)/2)
  # cat.print(length(ztd2))
  ztd3 = 0.5*exp((lnsigv2-lnsigu2)/2)
  # cat.print(length(ztd3))
  g1 = zr * ztd1 - myprod*exp(-lnsigu2/2)
  # cat.print(dim(g1))
  g2 = 1/2*exp(lnsigv2-lnsigu2) + zr*ztd2
  g3 = -1/2-1/2*exp(lnsigv2-lnsigu2) + zr*ztd3-1/2*myprod*e*exp(-lnsigu2/2)
  # cat.print(dim(x))
  g  = cbind(X*as.vector(g1), Zv*as.vector(g2), Zu*as.vector(g3))
  g0 <- colSums(g)
  
  ztd=zr*(zr+z)
  d11=ztd*ztd1*ztd1#, eq(1)
  # cat.print(dim(d11))
  d12=ztd*ztd2*ztd1  - zr*(-1/2*myprod*exp(-lnsigv2/2))#, eq(1,2)
  d13=ztd*ztd3*ztd1  - 1/2*myprod*exp(-lnsigu2/2)#, eq(1,3)
  d22=ztd*ztd2*ztd2  - zr*(-1/4*myprod*e*exp(-lnsigv2/2) - 1/4*exp((lnsigv2-lnsigu2)/2))  - 1/2*exp(lnsigv2-lnsigu2)#, eq(2)
  d23=ztd*ztd2*ztd3 - zr*(1/4*exp((lnsigv2-lnsigu2)/2))  + 1/2*exp(lnsigv2-lnsigu2)#, eq(2,3)
  d33=ztd*ztd3*ztd3  - zr*(-1/4*exp((lnsigv2-lnsigu2)/2)) - 1/2*exp(lnsigv2-lnsigu2) - 1/4*myprod*e*exp(-lnsigu2/2)#, eq(3)
  
  d12 <- as.vector(d12)
  d13 <- as.vector(d13)
  
  h0 <- matrix(0,k + kv + ku,k + kv + ku)
  # for(q in 1:length(y)){
  #   h0 <- h0 - rbind( cbind( tcrossprod(X[q,])*d11[q], X[q,]*d12[q], X[q,]*d13[q] ),
  #                     c( X[q,]*d12[q], d22[q], d23[q] ),
  #                     c( X[q,]*d13[q], d23[q], d33[q] ))
  # }
  for(q in 1:length(y)){
    # cat.print(tcrossprod(X[q,])*d11[q])
    # cat.print(X[q,]*d12[q])
    # cat.print(X[q,]*d13[q])
    # cat.print(c( X[q,]*d12[q], d22[q], d23[q] ))
    # cat.print(c( X[q,]*d13[q], d22[q], d33[q] ))
    # cat.print(cbind( tcrossprod(X[q,])*d11[q], X[q,]*d12[q], X[q,]*d13[q] ))
    # hoq <- rbind()
    # cat.print(crossprod(X[q,,drop=FALSE],Zv[q,,drop=FALSE])*d12[q])
    # cat.print(crossprod(X[q,,drop=FALSE],Zu[q,,drop=FALSE])*d13[q])
    dd11 <- tcrossprod(X[q,])*d11[q]
    dd12 <- crossprod(X[q,,drop=FALSE],Zv[q,,drop=FALSE])*d12[q]
    dd13 <- crossprod(X[q,,drop=FALSE],Zu[q,,drop=FALSE])*d13[q]
    dd22 <- tcrossprod(Zv[q,])*d22[q]
    dd33 <- tcrossprod(Zu[q,])*d33[q]
    dd23 <- crossprod(Zv[q,,drop=FALSE],Zu[q,,drop=FALSE])*d23[q]
    # cat.print(dd11)
    # cat.print(dd12)
    # cat.print(dd13)
    # cat.print(dd22)
    # cat.print(dd23)
    # cat.print(dd33)
    h0 <- h0 - 
      rbind( 
        cbind(dd11, dd12, dd13),
        cbind(t(dd12), dd22, dd23),
        cbind(t(dd13), t(dd23), dd33)
      )
  }
  
  
  list( grad = ifelse(is.na(g0) | is.infinite(g0), 0, g0),  hessian1 = ifelse(is.na(h0) | is.infinite(h0), 0, h0) )
  
}

# Technical efficiencies and prediction intervals

.u2efftnm <- function( e, su, sv, mu, alpha = 0.05, prod, cost.eff.less.one) {
 # if(prod){sn = -1} else {sn = 1}
 sn <- sn1 <- ifelse(prod, -1, 1)
 if(!prod & cost.eff.less.one) sn1 <- -1
 s  <- sqrt(su^2 + sv^2)
 m1 <- (sn*su^2 * e + mu*sv^2)/s^2
 s1 <- su * sv / s
 z  <- m1 / s1
 point.est.mean <- m1 + s1 * dnorm(-z) / pnorm(z)
 point.est.mode <- ifelse( m1 >= 0, m1, 0 )
 te_jlms_mean <- exp( sn1*point.est.mean)
 te_jlms_mode <- exp( sn1*point.est.mode)
 zl    <- qnorm( 1 - alpha / 2 * pnorm(z) )
 zu    <- qnorm( 1 - ( 1 - alpha/2 ) * pnorm(z) )
 te_l  <- exp( sn1*m1 + sn1* zl * s1 )
 te_u  <- exp( sn1*m1 + sn1* zu * s1 )
 te_bc <- exp( sn1*m1 + .5 * s1^2) * pnorm( sn1*s1 + z ) / pnorm(z)
 tymch <- data.frame(te_l, te_jlms_mean, te_jlms_mode,te_bc,te_u)
 if(!prod & !cost.eff.less.one) tymch <- data.frame(te_u, te_jlms_mean, te_jlms_mode,te_bc,te_l)
 colnames(tymch) <- c("Lower bound","JLMS", "Mode", "BC","Upper bound" )
 return(tymch)
}



# Marginal effects
.me = function(theta, Zu, Zdel, ku, kdel, n, dist = c("h", "t")){

 mat.equal <- function(x, y) is.matrix(x) && is.matrix(y) && ncol(x) == ncol(y) && all(colnames(x) == colnames(y))

 gu = theta[1:ku,1,drop = F]
 expu = exp(Zu%*%gu)
 if(dist == "t"){
  delta = theta[-c(1:ku),1,drop = F]
  mu <- Zdel%*%delta
  if(length(delta) == 1) delta = rep(0, max(ku, kdel))
 }
 if(length(gu) == 1) gu = rep(0, max(ku, kdel))
 meff = matrix(NA, ncol = max(ku, kdel) - 1, nrow = n)

 if(dist == "h"){
  if(ncol(Zu) == 1) {
   warning("Marginal effects are not returned: scale of half-normal distribution of inefficiency term is not expressed as function of any exogenous variables", call. = FALSE)
   meff = NULL
  } else {
   arg = sqrt(1/(2*pi)) * sqrt(expu)
   for(i in 2:ku){
    meff[,i-1] = arg*gu[i]
    meff = round(meff, digits = 5)
   }
   colnames(meff) = rownames(gu)[-1]
  }} else if(ncol(Zu) == 1 & ncol(Zdel) == 1){
   warning("Marginal effects are not returned: mean or variance of (pre-)truncated normal distribution of inefficiency term is not expressed as function of any exogenous variables", call. = FALSE)
   meff = NULL
  }else  if(dist == "t" & (mat.equal(Zu, Zdel)| all(delta == 0) | all(gu == 0))){
   arg = mu/sqrt(expu)
   g = dnorm(arg)/pnorm(arg)
   arg1 = (1 - arg*g - g^2); arg2 = sqrt(expu)/2 * ((1 + arg^2)*g + arg*g^2)
   for(i in 2:max(ku, kdel)){
    meff[,i-1] = arg1*delta[i] + arg2*gu[i]
    meff = round(meff, digits = 5)
   }
   if(all(gu == 0)){colnames(meff) = rownames(delta)[-1]} else {colnames(meff) = rownames(gu)[-1]}
  }  else {
   warning("Marginal effects are not returned: mean and variance of (pre-)truncated normal distribution of inefficiency term are expressed as functions not of the same exogenous variables", call. = FALSE)
   meff = NULL
  }
 return( meff)
}


.skewness <- function (x, na.rm = FALSE, type = 3)
{
 if (any(ina <- is.na(x))) {
  if (na.rm)
   x <- x[!ina]
  else return(NA)
 }
 if (!(type %in% (1:3)))
  stop("Invalid 'type' argument.")
 n <- length(x)
 # x <- x - mean(x)
 y <- sqrt(n) * sum(x^3)/(sum(x^2)^(3/2))
 if (type == 2) {
  if (n < 3)
   stop("Need at least 3 complete observations.")
  y <- y * sqrt(n * (n - 1))/(n - 2)
 }
 else if (type == 3)
  y <- y * ((1 - 1/n))^(3/2)
 y
}

# Print the estimation results

.printoutcs = function(x, digits, k, kv, ku, kdel, na.print = "NA", dist, max.name.length, mycutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), mysymbols = c("***", "**", "*", ".", " ")){

 Cf = cbind(ifelse(x[,1, drop = FALSE]> 999, formatC(x[,1, drop = FALSE], digits = 1, format = "e",width = 10), formatC(x[,1, drop = FALSE], digits = digits, format = "f", width = 10)),
            ifelse(x[,2, drop = FALSE]>999, formatC(x[,2, drop = FALSE], digits = 1, format = "e", width = 10), formatC(x[,2, drop = FALSE], digits = digits, format = "f", width = 10)),
            ifelse(x[,3, drop = FALSE]>999, formatC(x[,3, drop = FALSE], digits = 1, format = "e", width = 7), formatC(x[,3, drop = FALSE], digits = 2, format = "f", width = 7)),
            ifelse(x[,4, drop = FALSE]>999, formatC(x[,4, drop = FALSE], digits = 1, format = "e", width = 10), formatC(x[,4, drop = FALSE], digits = digits, format = "f", width = 10)),
            formatC(mysymbols[findInterval(x = x[,4], vec = mycutpoints)], flag = "-"))

 # mycutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1)
 # mysymbols = c("***", "**", "*", ".", " ")
 #
 # pvals <- m0$table[,4]
 #
 # findInterval(x = pvals, vec = mycutpoints)
 #
 # cbind(pvals,mysymbols[findInterval(x = pvals, vec = cutpoints)])
 #
 # pval_sym <- mysymbols[findInterval(x = x[,4], vec = cutpoints)]


 # cat("               Coef.        SE       z       P>|z|\n", sep = "")
 row.names(Cf) <- formatC(row.names(Cf), width = max(nchar(row.names(Cf))), flag = "-")
 cat("",rep(" ", max.name.length+6),"Coef.        SE       z       P>|z|\n", sep = "")
 dimnames(Cf)[[2]] <- rep("", dim(Cf)[[2]])
 cat("",rep("_", max.name.length+42-1),"", "\n", "Frontier", "\n", sep = "")
 print.default(Cf[1:k,,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 cat("",rep("-", max.name.length+42-1),"", "\n", "Random noise component: log(sigma_v^2)", "\n", sep = "")
 # dimnames(Cf)[[2]] <- rep("", dim(Cf)[[2]])
 print.default(Cf[(k+1):(k+kv),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 cat("",rep("-", max.name.length+42-1),"", "\n", "Inefficiency component: log(sigma_u^2)", "\n", sep = "")
 print.default(Cf[(k+kv+1):(k+kv+ku),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 if(dist == "t"){
  cat("",rep("-", max.name.length+42-1),"", "\n", "Mu (location parameter)", "\n", sep = "")
  print.default(Cf[(k+kv+ku+1):(k+kv+ku+kdel),,drop=F], quote = FALSE, right = TRUE, na.print = na.print)
 }
 if(nrow(Cf[-c(1:(k+kv+ku+kdel)),,drop=FALSE]) >= 1){
  cat("",rep("-", max.name.length+42-1),"", "\n", "Parameters of compound error distribution", "\n", sep = "")
  print.default(Cf[-c(1:(k+kv+ku+kdel)),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 }
 cat("",rep("_", max.name.length+42-1),"", "\n", sep = "")
 cat("Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n")
 invisible(x)
}



.timing <- function(x, wording = "Time elapsed is")
{
 if(x < 60){
  # cat("\n")
  cat("",wording,"",x," seconds\n", sep = "")
  # cat("\n")
 } else {
  if(x >= 60 & x < 60*60){
   minutes <- floor(x/60)
   seconds <- round(x - minutes * 60,1)
   # cat("\n")
   cat("",wording,"",minutes," minute(s) and ",seconds," second(s)\n", sep = "")
   # cat("\n")
  } else {
   if(x >= 60*60 & x < 60*60*24){
    hours   <- floor(x / 60 / 60)
    minutes <- round( (x - hours * 60 *60) / 60, 1)
    seconds <- floor(x - hours * 60 *60 - minutes * 60)
    # cat("\n")
    cat("",wording,"",hours," hour(s) and ",minutes," minute(s) \n", sep = "")
    # cat("\n")
   } else {
    if(x >= 60*60*24){
     days    <- floor(x / 60 / 60 / 24)
     hours   <- round( (x - days * 60 * 60 * 24) / 60 /60 ,1)
     minutes <- floor( (x - days * 60 * 60 * 24 - hours * 60 *60) / 60)
     seconds <- floor(x - days * 60 * 60 * 24 - hours * 60 *60 - minutes * 60)
     # cat("\n")
     cat("",wording,"",days," day(s) and ",hours," hour(s)\n", sep = "")
     # cat("\n")
    }
   }
  }
 }
}

# library(matrixcalc)


.mlmaximize <- function(theta0, ll, gr = NULL, hess = NULL, alternate = NULL, BHHH = F, level = 0.99, step.back = .Machine$double.eps, reltol = .Machine$double.eps, lmtol = sqrt(.Machine$double.eps), steptol = sqrt(.Machine$double.eps), digits = 4, when.backedup = sqrt(.Machine$double.eps), max.backedup = 17, print.level = 6, only.maximize = FALSE, maxit = 150, n = 100, ...){

 theta00 <- theta0

 k4 <- length(theta0)

 if(print.level >= 6){
  cat("\n=================")
  cat(" Initial values:\n\n", sep = "")
  print(theta0)
 }

#  if(print.level >= 2){
#   cat("\n=================")
#   cat(" Maximization:\n\n", sep = "")
#  }

 # step.back = 2^-217

 ll0 <- ll(theta0, ...)
 ltol <- reltol * (abs(ll0) + reltol)
 typf <- ll0
 theta1 <- theta0

 iter <- iter.total <- backedup <- backedups <- wasconcave <- wasconcaves <- 0

 if( is.na(ll0) | ll0 == -Inf ){
  if(print.level >= 2){
   cat("Could not compute ll at starting values: trying something else\n")
  }
  iter1 <- backedups
  repeat{
   iter1 <- iter1 + 1
   theta0 <- theta0 * runif(length(theta0), 0.96, 1.05) # not sure what to do here
   ll0 <- ll( theta0, ... )
   if(!is.na(ll0)) break
   if(iter1 == 55){
    stop("it's not gonna happen...")
   }
  }
  # backedups <- iter1
 }

 delta1 <- gHg <- s1 <- 1
 h1 <- tryCatch( 2, error = function(e) e )
 cant.invert.hess <- FALSE

 if(print.level >= 2){
  cat(paste("Iteration ",formatC(iter, width = 3)," (at starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
 }

 repeat{
  iter.total <- iter.total + 1
  # cat("backedup = ",backedup," backedups = ",backedups,"\n", sep = "")
  if(print.level >= 6){
   print(theta0)
  }

  # cumulate how many times did it backed-up in a row
  if(s1 < when.backedup){
   backedup <- backedup + 1
  } else {
   backedup <- 0
  }
  # print(s1)
  # cumulate how many times was concave
  if( inherits(h1, "error") ){
   wasconcave <- wasconcave + 1
  } else {
   wasconcave <- 0
  }

  # try different values if was concave more than @@@ times
  if(wasconcave == max.backedup){
   # start over
   if(print.level >= 2){
    cat("Not concave ",max.backedup," times in a row: trying something else (not concave ",wasconcaves+1," times in total)\n")
   }
   iter <- wasconcave <- backedup <- 0
   wasconcaves <- wasconcaves + 1
   theta0 <- theta0*runif(length(theta0), 0.96, 1.05) # not sure what to do here
   ll0 <- ll( theta0, ... )
   # if(print.theta) print(theta0)
   if(print.level >= 2){
    cat(paste("Iteration ",formatC(iter, width = 3),"  (at slightly perturbed starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
   }
  }

  # try different values if backed-up more than @@@ times
  if(backedup == max.backedup){
   # start over
   if(print.level >= 2){
    cat("Backed-up ",max.backedup," times in a row: trying something else (backup-up ",backedups+1," times in total)\n", sep = "")
   }
   iter <- backedup <- wasconcave <- 0
   backedups <- backedups + 1
   wasconcaves <- wasconcaves + 1
   theta0 <- theta0*runif(length(theta0), 0.96, 1.04) # not sure what to do here
   ll0 <- ll( theta0, ... )
   if(print.level >= 6){
    print(theta0)
   }
   if(print.level >= 2){
    cat(paste("Iteration ",formatC(iter, width = 3)," (at slightly perturbed starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
   }
  }

  # see if it calculated ll
  if( is.na(ll0) | ll0 == -Inf | ll0 == Inf | ll0 == 0 ){
   if(print.level >= 2){
    cat("Could not compute ll: trying something else\n")
   }
   iter1 <- backedups
   repeat{
    iter1 <- iter1 + 1
    # theta0 <- c( cons0, beta0, mu = 0, eta = 0, lnsv2 = -1*iter1/2, lnsu2 = -1*iter1/2)
    theta0 <- theta00*runif(length(theta0), 0.96, 1.04) # not sure what to do here
    ll0 <- ll( theta0, ... )
    if(!is.na(ll0) & ll0 != 0) break
    if(iter1 == 15){
     stop("it's not gonna happen... could not compute at initial and find feasible values")
    }
   }
   iter <- backedup <- 0
   backedups <- iter1
   if(print.level >= 2){
    cat(paste("Iteration ",formatC(iter, width = 3)," (at slightly perturbed starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
   }
  }
  if(backedups == 5){
   stop("it's not gonna happen... backuped 5 times")
  }
  if(wasconcaves == 5){
   stop("it's not gonna happen... not concave 5 times")
  }
  iter <- iter + 1
  delta3 <- 1
  # step 2: direction vector

  # previous theta

  if(iter.total > 1) theta1 <- theta0 - s1 * d0

  # BHHH (faster, but different):
  # The Hessian is approximated as the negative
  # of the sum of the outer products of the gradients
  # of individual observations,
  # -t(gradient) %*% gradient = - crossprod( gradient )
  g1 <- gr(theta0, ...)
  # print(g1)
  if(!is.null(alternate)) BHHH <- floor(iter/alternate) %% 2 == 1
  h0 <- hess(theta0,  ...)
  # print(h0)
  # check if negative definite
  eigen1 <- eigen( h0 )
  # eigen.tol <- k4 * max(abs(eigen1$values)) * .Machine$double.eps # this is for positive definiteness
  eigen.val <- ifelse(eigen1$values < .Machine$double.eps^.1, 0, eigen1$values)
  # hess.pos.def <- sum(eigen1$values > eigen.tol) == k4
  hess.neg.def <- !any(eigen.val >= 0)
  # 1. replace negative with small ones
  # eigen.val <- ifelse(eigen1$values < 0, .0001, eigen.val)
  # 2. replace negative with absolut values
  eigen.val <- abs(eigen1$values)
  # print(hess.neg.def)
  # make it negative definite if it is not already
  if(!hess.neg.def){
   h0_ <- matrix(0, k4, k4)
   # eigen1 <- eigen( h0 )
   for( i in seq_len( k4 ) ){
    # h0_ <- h0_ - abs(eigen1$values[i]) * eigen1$vectors[,i] %*% t(eigen1$vectors[,i])
    h0_ <- h0_ - eigen.val[i] * eigen1$vectors[,i] %*% t(eigen1$vectors[,i])
   }
   h0.old <- h0
   h0 <- h0_
  }
  # print( is.negative.definite(h0) )
  # remember hessian and negative of its inverse from previous iter that could have been inverted
  if( !cant.invert.hess ){
   h0_previous <- h0
   h1_previous <- h1
  }
  # easier to invert positive definite matrix
  h1 <- tryCatch( qr.solve(-h0), error = function(e) e )
  # check if it can be inverted
  cant.invert.hess <- FALSE
  cant.invert.hess <- inherits(h1, "error")
  if( cant.invert.hess ){
   # print(h1)
   if(print.level >= 2){
    cat(paste("cannot invert Hessian, using eigenvalues\n", sep = ""), sep = "")
   }
   # this was just to get the uninvertable hessian
   # return(list(hess = h0, grad = g1))
   # stop("this")
   # @14@ this
   eig1 <- eigen( -h0_previous )
   d0 <- rep(0, length(theta0))
   # eig2 <- ifelse(eig1$values < eps1, 1, eig1$values)
   for (i in 1:length(eig1$values)){
    d0 <- d0 + (g1%*%eig1$vectors[,i])%*%eig1$vectors[,i] / eig1$values[i]
   }
   # @14@ could be done easier
   # d0 <- qr.solve(-h0, g1, tol = 1e-10)
   gHg <- sum( g1 * d0)
   # in the part of the ortogonal subspace where the eigenvalues
   # are negative or small positive numbers, use steepest ascent
   # in other subspace use NR step
   # d0 <- ifelse(eigen(-h0, only.values = TRUE)$values < reltol, g1, d0)
   gg <- sqrt( crossprod(g1) )
   gHg <- gg
   # d0 <- g1
   # d0
  } else {
   d0 <- as.vector( h1 %*% g1 )
   gg <- sqrt( crossprod(g1) )
   # h1.old <- solve(-h0.old)
   gHg <- as.vector( t(g1) %*% h1 %*% g1 )
  }
  # gg_scaled <- gg * max( crossprod(theta0), crossprod(theta1) ) / max( abs(ll0), abs(typf))
  # theta_rel <- max( abs(theta0 - theta1) / apply( cbind( abs(theta0),abs(theta1) ), 1, max) )
  theta_rel <- max( abs(theta0 - theta1) / (abs(theta1)+1) )


  # begin stopping criteria calculated using new values of g1 and h1
  if(s1 > when.backedup*10^-100 & delta1 != 17.17){ # if(s1 > when.backedup*10^-100 & !cant.invert.hess){
   if(abs(gHg) < lmtol & iter.total > 1){
    conv.crite <- 1
    if(print.level >= 2){
     cat("\nConvergence given g inv(H) g' = ",abs(gHg)," < lmtol\n", sep = "")
    }
    break
   }
   if(theta_rel < steptol & iter.total > 2){
    conv.crite <- 3
    # print(theta_rel)
    if(print.level >= 2){
     cat("\nConvergence given relative change in parameters = ",theta_rel," < steptol\n", sep = "")
    }
    break
   }
  }
  # end stopping criteria
  # use steepest ascent when backed-up
  if(s1 < when.backedup*10^-3){
   # eig1 <- eigen( -h0 )
   d0 <- g1
   # d0 <- ifelse(eig1$values < reltol, g1, d0)
   # theta0 <- theta0 - 1 * d0
  }
  # print(d0)
  # step 3: new guess
  # a: s = 1
  # b: funct(theta0 + d0) > funct(theta0)
  s1 <- 1
  theta1 <- theta0 + s1 * d0
  # print(12)
  # print(theta1)
  ll1 <- ll( theta1, ... )
  # print(13)
  delta2 <- ll1 - ll0
  flag <- (!is.na(delta2) & delta2 != -Inf & delta2 > 0)
  # begin Cases
  if( flag ){
   # begin Case 1: f(theta1) > f(theta0)
   ll.temp <- ll0
   # check if s1 = 2, 3, ... increases f even more
   while( flag ){
    if(print.level >= 6){
     cat(paste("\t\tCase 1: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
    }
    ll0 <- ll1
    s1 <- s1 + 1
    theta1 <- theta0 + s1 * d0
    ll1 <- ll( theta1, ... )
    delta2 <- ll1 - ll0
    flag <- (!is.na(delta2) & delta2 != -Inf & delta2 > 0)
   }
   # overall delta
   delta1 <- ll0 - ll.temp
   delta_rel <- abs(delta1 / ll.temp)
   # print(delta_rel)
   s1 <- s1 - 1
   # overwrite the values
   theta0 <- theta0 + s1 * d0
   # end Case 1: f(theta1) > f(theta0)
  } else {
   # begin Case 2: f(theta1) < f(theta0)
   # check only if s1=1/2 increases f
   s1 <- 0.5
   theta1 <- theta0 + s1 * d0
   ll1 <- ll( theta1, ... )
   # cat(" ll1 = ",ll1,", ll0 = ",ll0,", s = ",s1,"\n", sep = "")
   delta2 <- ll1 - ll0
   if(print.level >= 6){
    cat(paste("\t\tCase 2: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
   }
   flag2 <- (!is.na(delta2) & delta2 != -Inf & delta2 > 0)
   # end Case 2: f(theta1) < f(theta0)
   if( flag2 ){
    # begin Case 2a: f(theta1) > f(theta0)
    ll.temp <- ll0
    # check if s1=1/2^2,1/2^3,... increases f even more
    while( flag2 ){
     if(print.level >= 6){
      cat(paste("\t\t\tCase 2a: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
     }
     ll0 <- ll1
     s1 <- 0.5 * s1
     theta1 <- theta0 + s1 * d0
     ll1 <- ll( theta1, ... )
     # cat(" ll1 = ",ll1,", ll0 = ",ll0,", s = ",s1,"\n", sep = "")
     delta2 <- ll1 - ll0
     flag2 <- (!is.na(delta2) & delta2 != -Inf & delta2 > ltol)
    }
    # overall delta
    delta1 <- ll0 - ll.temp
    delta_rel <- abs(delta1 / ll.temp)
    # print(delta_rel)
    s1 <- 2 * s1
    # overwrite the values
    theta0 <- theta0 + s1 * d0
    # end Case 2a: f(theta1) > f(theta0)
   } else {
    # begin Case 2b: f(theta1) < f(theta0)
    ll.temp <- ll0
    # try s1=1/2^2,1/2^3,... so that f(theta1) > f(theta0)
    while ( !flag2 & s1 > step.back ){
     s1 <- 0.5 * s1
     theta1 <- theta0 + s1 * d0
     ll1 <- ll( theta1, ... )
     delta2 <- ll1 - ll0
     if(print.level >= 6){
      cat(paste("\t\t\tCase 2b: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
     }
     flag2 <- (!is.na(delta2) & delta2 != -Inf & delta2 > 0)
    }
    if( !flag2 | s1 < step.back ){
     # stop("provide different starting values")
     delta1 <- 17.17
    } else {
     # overwrite the values
     delta1 <- delta2
     delta_rel <- abs(delta1 / ll.temp)
     ll0 <- ll1
     theta0 <- theta0 + s1 * d0
    }
    # end Case 2b: f(theta1) < f(theta0)
   }
  }


  if(print.level >= 2){
   if( cant.invert.hess ){
    cat(paste("Iteration ",formatC(iter, width = 3)," (hessian is ",ifelse(BHHH, "BHHH", "analytical"),", ",formatC(iter.total, width = 3)," in total):   log likelihood = ",format(ll0, digits = 13)," (not concave)\n", sep = ""), sep = "")
   } else if (s1 <= when.backedup) {
    cat(paste("Iteration ",formatC(iter, width = 3)," (hessian is ",ifelse(BHHH, "BHHH", "analytical"),", ",formatC(iter.total, width = 3)," in total):   log likelihood = ",format(ll0, digits = 13)," (backed up)\n", sep = ""), sep = "")
   } else {
    cat(paste("Iteration ",formatC(iter, width = 3)," (hessian is ",ifelse(BHHH, "BHHH", "analytical"),", ",formatC(iter.total, width = 3)," in total):   log likelihood = ",format(ll0, digits = 13),"\n", sep = ""), sep = "")
   }
  }

  # printing criteria
  if(print.level >= 5){
   if( cant.invert.hess ){
    cat(paste(" (in iter ",formatC(iter, width = 3),": delta = ",format(delta1, digits = 6),"; s = ",format(s1, digits = 6),"; quasi-gHg = ",format(gHg, digits = 6),"; theta_rel_change = ",format(theta_rel, digits = 6),")\n\n", sep = ""), sep = "")
    if(print.level >= 5.5){
     print(theta0)
     cat("\n")
    }
   } else {
    cat(paste(" (in iter ",formatC(iter, width = 3),": delta = ",format(delta1, digits = 6),"; s = ",format(s1, digits = 6),"; gHg = ",format(gHg, digits = 6),"; theta_rel_change = ",format(theta_rel, digits = 6),")\n\n", sep = ""), sep = "")
    if(print.level >= 5.5){
     print(theta0)
     cat("\n")
    }
   }
  }
  # print(s1)
  if(s1 > when.backedup^2 & !cant.invert.hess){ # if(s1 > when.backedup^2 & !cant.invert.hess){
   # ltol <- reltol * (abs(ll0) + reltol)
   # print(cant.invert.hess)
   if(delta1 > 0 & !is.na(delta_rel) & delta_rel < ltol & iter.total > 1){
    conv.crite <- 2
    if(print.level >= 2){
     cat("\nConvergence given relative change in log likelihood = ",delta_rel," < ltol\n", sep = "")
    }
    break
   }
  }
  if(iter.total > maxit){
   cat("\n Maximum number of iterations (",maxit,") reached without convergence\n", sep = "")
   break
  }
 } # end repeat

 if( !only.maximize & !cant.invert.hess){
  names(ll0) <- NULL
  colnames(h1) <- rownames(h1) <- names(g1) <- names(theta0)

  # sqrt(crossprod(g1))

  b0 <- theta0
  sd0 <- sqrt( diag( h1 ) )
  t0 <- b0 / sd0
  p0 <- pt(abs(t0), n-length(b0), lower.tail = FALSE) * 2
  t10 <- qt((1-0.01*level)/2, n-length(b0), lower.tail = FALSE)
  t17 <- cbind( b0, sd0, t0, p0, b0 - t10*sd0, b0 + t10*sd0)
  # t17 <- cbind( b0, sd0, t0, p0)
  colnames(t17) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)", paste("",level,"_CI_LB", sep = ""), paste("",level,"_CI_UB", sep = ""))
  # colnames(t17) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
  # t17

  if(print.level >= 2){
   cat(paste("\nFinal log likelihood = ",formatC(ll0,digits=7,format="f"),"\n\n", sep = ""), sep = "")
  }
  # cat(paste("Stoc. frontier normal/",distribution,"\n", sep = ""), sep = "")
  if(print.level >= 5.5){
   cat("\nCoefficients:\n\n", sep = "")
   printCoefmat(t17[,1:4], digits = digits)
  }

  return(list(par = theta0, table = t17, gradient = g1, vcov = h1, ll = ll0, gg = gg, gHg = gHg, delta_rel = delta_rel, ltol = ltol, theta_rel_ch = theta_rel, conv.crite = conv.crite))
 } else {
  return(list(par = theta0, gradient = g1, vcov = h1, ll = ll0, gg = gg, gHg = gHg, delta_rel = delta_rel, ltol = ltol, theta_rel_ch = theta_rel, conv.crite = conv.crite))
 }


}
.su1 <- function(x, mat.var.in.col = TRUE, digits = 4, probs = c(0.1, 0.25, 0.5, 0.75, 0.9), print = TRUE, transpose = FALSE){

 xvec2 <- xvec1 <- FALSE

 if(is.matrix(x)){
  if(min(dim(x)) == 1){
   xvec1 <- TRUE
   # mynames <- deparse(substitute(x))
   x <- as.vector(x)
  } else {
   if(!mat.var.in.col){
    x <- t(x)
   }
   mynames <- paste("Var", seq_len(ncol(x)), sep = "")
  }
  # print(x)
  # mynames <- colnames(x)
 } # end if matrix

 if(is.vector(x)){
  xvec2 <- TRUE
  mynames <- deparse(substitute(x))
  x <- data.frame(Var1 = x)
 } # end if vector

 # cat("nymanes", sep ="")
 # print(mynames)

 if(!is.vector(x) & !is.matrix(x) & !is.data.frame(x)){
  stop("Provide vector, matrix, or data.frame")
 } else {
  t1 <- apply(x, 2, function(x) c(Obs = length(x), NAs = sum(is.na(x)), Mean = mean(x, na.rm = TRUE), StDev = sd(x, na.rm = TRUE), IQR = IQR(x, na.rm = TRUE), Min = min(x, na.rm = TRUE), quantile( x, probs = probs, na.rm = TRUE ), Max = max(x, na.rm = TRUE)))
  # print(class(t1))
  # print(dim(t1))
  if(xvec2 & !xvec1) colnames(t1) <- mynames

  if(print){
   if(transpose){
    tymch <- formatC(t1, digits = 4, format = "f", width = 4+1)
   } else {
    tymch <- formatC(t(t1), digits = 4, format = "f", width = 4+1)
   }
   tymch <- gsub(".0000","",tymch, fixed = TRUE)
   print(tymch, quote = FALSE)
  }
 }
 return(t1)
}

.prepareYXZ.panel.simple <- function(formula, data, it, subset,
                                     ln.var.v.it = ln.var.v.it, 
                                     ln.var.u.0i = ln.var.u.0i, 
                                     mean.u.0i = mean.u.0i,
                                     print.level = 1, ...) {
 needed.frame <- sys.nframe() - 1
 mf0 <- match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = 1)))
 # print(mf0)
 
 # check if it is a matrix
 datasupplied <- !(match("data", names(mf0), 0) == 0)
 subssupplied <- !(match("subset", names(mf0), 0) == 0)
 
 if(subssupplied & !datasupplied){
  stop("Cannot specify 'subset' without specifying 'data'\n")
 }
 
 itsupplied <- !(match("it", names(mf0), 0) == 0)
 if(!itsupplied){
  stop("Panel structure must be specified by using 'it' argument\n")
 }
 
 if(length(it) != 2){
  stop("Invalid panel structure in 'it' argument\n")
 }
 
 if ( is.null(ln.var.v.it) ){
  ln.var.v.it <- ~ 1
  kvi <- 1
 } else {
  kvi <- 17
 }
 if ( is.null(ln.var.u.0i) ){
  ln.var.u.0i <- ~ 1
  ku0 <- 1
 } else {
  ku0 <- 17
 }
 if ( is.null(mean.u.0i) ){
  mean.u.0i <- ~ 1
  kdeli <- 1
 } else {
  kdeli <- 17
 }
 
 form3 <- as.formula(paste0("~", paste0(it, collapse = "+")))
 form1 <- as.Formula(formula,ln.var.v.it,ln.var.u.0i,mean.u.0i,form3)
 
 # print(form2)
 
 if(datasupplied){
  data.order <- as.numeric(rownames(data)) # seq_len(nrow(data))
  
  # print(rownames(data)[1:100])
  # print(as.numeric(rownames(data))[1:100])
  
  # begin get a logical vector equal TRUE if !missing
  # first using data and subset to get x without NA
  mf <- mf0
  mf$formula <- form1 #formula( form )
  m <- match(c("formula", "data", "subset"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  #   cat("Print mt\n")
  #   print(mt)
  #   cat("Print mf\n")
  #   print(mf)
  X <- model.matrix(mt, mf)
  # Y <- as.matrix(model.response(mf, "numeric"))
  # print(dim(X))
  # X <- model.frame(mt, mf)
  #   cat("Print X\n")
  #   print(X)
  #   cat("End print X\n")
  # now get the names in the entire data
  esample.nu <- as.numeric(rownames(X))
  # print(esample.nu[1:100])
  # print(old.order)
  esample <- data.order %in% esample.nu
  # print(length(esample))
  # print(table(esample))
  if(sum(esample)==0){
   stop("no valid data points")
  }
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing
  
  # get the data
  # print(form1)
  #   dataesample <- model.frame(form1, data = data)
  #   dataesample <- dataesample[esample.nu,]
  #   print(as.numeric(rownames(dataesample)))
  # esample1 <- esample
  #   esample <- as.numeric(rownames(dataesample)) %in% esample.nu
  #   # print(table(esample))
  #   print(c(length(esample), nrow(dataesample)))
  # dataesample <- dataesample[as.numeric(rownames(dataesample)) %in% esample.nu,]
  # print(dim(dataesample))
  # print(dataesample)
  Y <- as.matrix(model.part(form1, data = mf, lhs = 1, drop = FALSE))
  # print(length(Y))
  # print(length(model.response(dataesample)))
  # Y <- as.matrix( model.matrix(formula(form1, lhs = 1, rhs = 0), data = data[esample,]))
  # print(Y)
  X <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 1), data = mf))
  # print(X)
  nt <- nrow(Y)
  k <- ncol(X)
  # Zvi
  if(is.null(ln.var.v.it)){
   Zvi <- matrix(1, nrow = sum(esample), ncol = 1)
  }
  else {
   Zvi <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 2), data = mf))
   kvi <- ncol(Zvi)
  }
  # Zu0
  if(is.null(ln.var.u.0i)){
   Zu0 <- matrix(1, nrow = sum(esample), ncol = 1)
  }
  else {
   Zu0 <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 3), data = mf))
   ku0 <- ncol(Zu0)
  }
  # kdeli
  if(is.null(mean.u.0i)){
   zdeli <- matrix(1, nrow = sum(esample), ncol = 1)
  }
  else {
   zdeli <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 4), data = mf))
   kdeli <- ncol(zdeli)
  }
  # IT
  itvar <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 5), data = mf))[,-1]
  sorted <- order(itvar[,1],itvar[,2])
  old.order <- seq(nrow(itvar))
  old.sorted <- order(old.order[sorted])
  itvar <- itvar[sorted,]
  idvar <- itvar[,1]
  timevar <- itvar[,2]
  ii <- unique(idvar)
  n <- length(ii)
  t.years <- length(unique(itvar[,2]))
  
  # The size of each panel
  
  t0 <- as.vector( by(data = itvar[,1], INDICES = itvar[,1], function(x) sum(!is.na(x))) )
  t1 <- summary( t0 )
  
  # summary of the panel data
  
  # if(print.level >= 1){
  #
  #  cat("\n=================")
  #  cat(" Summary of the panel data: \n\n", sep = "")
  #  cat("   Number of obs       (NT) =",nt,"", "\n")
  #  cat("   Number of groups     (N) =",n,"", "\n")
  #  cat("   Obs per group: (T_i) min =",t1["Min."],"", "\n")
  #  cat("                        avg =",t1["Mean"],"", "\n")
  #  cat("                        max =",t1["Max."],"", "\n\n")
  # }
  
  ids <- sort( unique(idvar) )
  idlenmax <- t1["Max."]
  dat.descr <- c(NT = nt, N = n, t1["Min."], t1["Mean"], t1["Max."])
  #   print(head(Y))
  #   print(head(X))
  #   print(head(Zu))
  #   print(head(Zv))
  #   print(head(Zdel))
 }
 # if data are not supplied
 else {
  # begin get a logical vector equal TRUE if !missing
  
  # first using data and subset to get XZ without NA
  mf <- mf0
  mf$formula <- formula( formula )
  subsetsupplied <- !(match("subset", names(mf), 0) == 0)
  if(subsetsupplied) stop("Subset with matrices is not allowed")
  m <- match(c("formula"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  # mf <- eval(mf, parent.frame())
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  Y <- as.matrix(model.response(mf))
  n <- nrow(Y)
  XZ <- as.matrix(model.matrix(mt, mf))
  # get a logical vector equal TRUE if !missing
  with.na <- model.frame(mt, na.action = na.pass)
  esample.nu <- as.numeric(rownames(XZ))
  esample <- rownames(with.na) %in% esample.nu
  if(sum(esample)==0){
   stop("no valid data points")
  }
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing
  
  # get the data
  #   print(head(Y))
  #   print(head(XZ))
  
  # get X
  mf <- mf0
  m <- match(c("formula"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  X <- as.matrix(model.matrix(mt, mf))
  # print(head(X))
  k <- ncol(X)
  
  # get Zvi
  if(is.null(ln.var.v.it)){
   Zvi <- matrix(1, nrow = sum(esample), ncol = 1)
  }
  else {
   mf <- mf0
   mf$formula <- formula( Formula(as.formula(paste("",deparse(formula)," + ",ln.var.v.it[2],"", sep = ""))))
   m <- match(c("formula"), names(mf), 0L)
   mf <- mf[c(1L, m)]
   mf[[1L]] <- as.name("model.frame")
   mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
   mt <- attr(mf, "terms")
   XZs <- as.matrix(model.matrix(mt, mf))
   Zvi <- XZs[, -(2:k), drop = FALSE]
   kvi <- ncol(Zvi)
   # print(head(Zu))
  }
  # get Zu0
  if(is.null(ln.var.u.0i)){
   Zu0 <- matrix(1, nrow = sum(esample), ncol = 1)
  }
  else {
   mf <- mf0
   mf$formula <- formula( Formula(as.formula(paste("",deparse(formula)," + ",ln.var.v.it[2]," + ",ln.var.u.0i[2],"", sep = ""))))
   m <- match(c("formula"), names(mf), 0L)
   mf <- mf[c(1L, m)]
   mf[[1L]] <- as.name("model.frame")
   mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
   mt <- attr(mf, "terms")
   XZs <- as.matrix(model.matrix(mt, mf))
   Zu0 <- XZs[, -(2:(k+kvi-1)), drop = FALSE]
   ku0 <- ncol(Zu0)
   # print(head(Zv))
  }
  # get zdeli
  if(is.null(mean.u.0i)){
   zdeli <- matrix(1, nrow = sum(esample), ncol = 1)
  }
  else {
   mf <- mf0
   mf$formula <- formula( Formula(as.formula(paste("",deparse(formula)," + ",ln.var.v.it[2]," + ",ln.var.u.0i[2]," + ",mean.u.0i[2],"", sep = ""))))
   m <- match(c("formula"), names(mf), 0L)
   mf <- mf[c(1L, m)]
   mf[[1L]] <- as.name("model.frame")
   mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
   mt <- attr(mf, "terms")
   XZs <- as.matrix(model.matrix(mt, mf))
   zdeli <- XZs[, -(2:(k+kvi-1+ku0-1)), drop = FALSE]
   kdeli <- ncol(zdeli)
   # print(head(Zv))
  }
  # IT
  itvar <- XZ[, -(2:(k+kvi-1+ku0-1+kvi-1+kdeli-1)), drop = FALSE]
  sorted <- order(itvar[,1],itvar[,2])
  itvar <- itvar[sorted,]
  idvar <- itvar[,1]
  timevar <- itvar[,2]
  ii <- unique(idvar)
  n <- length(ii)
  t.years <- length(unique(itvar[,2]))
  
  # The size of each panel
  
  t0 <- as.vector( by(data = itvar[,1], INDICES = itvar[,1], function(x) sum(!is.na(x))) )
  t1 <- summary( t0 )
  
  # summary of the panel data
  
  # if(print.level >= 1){
  #  cat("\n=================")
  #  cat(" Summary of the panel data: \n\n", sep = "")
  #  cat("   Number of obs       (NT) =",nt,"", "\n")
  #  cat("   Number of groups     (N) =",n,"", "\n")
  #  cat("   Obs per group: (T_i) min =",t1["Min."],"", "\n")
  #  cat("                        avg =",t1["Mean"],"", "\n")
  #  cat("                        max =",t1["Max."],"", "\n")
  # }
  
  ids <- sort( unique(idvar) )
  idlenmax <- t1["Max."]
  dat.descr <- c(NT = nt, N = n, t1["Min."], t1["Mean"], t1["Max."])
  
  #   print(head(Zu))
  #   print(head(Zv))
  #   print(head(Zdel))
 }
 
 # sort all
 Y <- Y[sorted,, drop = FALSE]
 X <- X[sorted,, drop = FALSE]
 
 Zu0 <- Zu0[sorted,, drop = FALSE]
 Zvi <- Zvi[sorted,, drop = FALSE]
 zdeli <- zdeli[sorted,, drop = FALSE]
 # print(head(zdeli))
 
 abbr.length <- 121
 names_x  <- abbreviate(colnames(X), abbr.length, strict = TRUE, dot = FALSE, method = "both.sides")
 names_u0 <- abbreviate(colnames(Zu0), abbr.length, strict = TRUE, dot = FALSE, method = "both.sides")
 names_vi <- abbreviate(colnames(Zvi), abbr.length, strict = TRUE, dot = FALSE, method = "both.sides")
 names_udeli <- abbreviate(colnames(zdeli), abbr.length, strict = TRUE, dot = FALSE, method = "both.sides")
 # print(names_udeli)
 max.name.length <- max(nchar(c(names_x,names_u0,names_vi,names_udeli)))
 
 names_u0 <- colnames(Zu0)
 names_udeli <- colnames(zdeli)
 
 # print(names_udeli)
 
 # print(c(dim(Zv0), length(idvar)))
 
 # print(as.vector(unlist( by(data = Zv0[,-1,drop = F], INDICES = idvar, function(x) apply(x, 2, sd)) )))
 
 if(!is.null(mean.u.0i)){
  # print(as.vector(unlist(by(data = zdeli[,-1, drop = F], INDICES = idvar, FUN = function(x) apply(x, 2, sd)))))
  if(sum(as.vector(unlist(by(data = zdeli[,-1, drop = F], INDICES = idvar, FUN = function(x) apply(x, 2, sd)))), na.rm = TRUE) > 0){
   stop("Time-varying variable(s) in 'mean.u.0i'", call. = FALSE)
  } 
 }
 zdeli <- matrix(unlist(by(zdeli, idvar, colMeans)), nrow = n, byrow = TRUE)
 colnames(zdeli) <- names_udeli
 # print(head(zdeli))
 
 
 
 if(!is.null(ln.var.u.0i)){
  # print(as.vector(unlist(by(data = Zu0[,-1, drop = F], INDICES = idvar, FUN = function(x) apply(x, 2, sd)))))
  if(sum(as.vector(unlist(by(data = Zu0[,-1, drop = F], INDICES = idvar, FUN = function(x) apply(x, 2, sd)))), na.rm = TRUE) > 0){
   stop("Time-varying variable(s) in 'ln.var.u.0i'", call. = FALSE)
  }
 }
 Zu0 <- matrix(unlist(by(Zu0, idvar, colMeans)), nrow = n, byrow = TRUE)
 colnames(Zu0) <- names_u0
 # print(head(Zu0))
 
 
 # esample <- esample[sorted,, drop = FALSE]
 # esample.nu <- esample.nu[sorted]
 
 # colnames(X)[1] <- colnames(Zvi)[1] <- colnames(Zu0)[1] <- colnames(zdeli)[1] <- "(Intercept)"
 
 
 
 tymch <- list(Formula = form1, Y = Y, X = X,
               nt = nt, n = n, Kb = k,
               idvar = idvar, timevar = timevar, ii = ii, 
               ids = ids, idlenmax = idlenmax,
               t0 = t0, itvar = itvar,
               dat.descr = dat.descr, t_i = t1, old.sorted = old.sorted,
               Zvi = Zvi, Zu0 = Zu0, Zdeli = zdeli,
               Kvi = kvi, Ku0 = ku0, Kdeli = kdeli,
               esample = esample, esample.nu = esample.nu)
 
 # cat("3\n")
 # print(head(zdeli))
 
 if(print.level >= 1){
  est.spd.left <- floor( (max.name.length+42-29) / 2 )
  est.spd.right <- max.name.length+42-29 - est.spd.left
  # cat("\n------------")
  # cat(" Summary of the panel data: ------------\n\n", sep = "")
  # cat("\n------------ Summary of the panel data: ----------\n\n", sep = "")
  cat("\n",rep("-", est.spd.left)," Summary of the panel data: ",rep("-", est.spd.right),"\n\n", sep ="")
  cat("   Number of obs       (NT) =",nt,"", "\n")
  cat("   Number of groups     (N) =",n,"", "\n")
  cat("   Obs per group: (T_i) min =",dat.descr[3],"", "\n")
  cat("                        avg =",dat.descr[4],"", "\n")
  cat("                        max =",dat.descr[5],"", "\n")
 }
 # cat("0")
 class(tymch) <- "npsf"
 return(tymch)
 
}

# Log-likelihood

.ll.panel <- function(theta, prod, coef.fixed, my.n, Ktheta, k, kv, ku, kdel, yit, zvit, zui, xit, zdeli, model, timevar, maxT, ids, idvar, t0, eff.time.invariant = FALSE, mean.u.0i.zero = FALSE) {
 
 # print(theta)

 # theta <- theta[!coef.fixed]
 # print(coef.fixed)
 s <- ifelse(prod, -1, 1)
 
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[(k+kv+1):(k+kv+ku)]
 
 eit <- as.vector( yit - xit%*%beta )
 # print(summary(as.vector(eit)))
 sigu2i <- as.vector( exp(zui%*%gu) )
 sigv2it <- as.vector( exp(zvit%*%gv) )
 
 if(!mean.u.0i.zero){
  delta <- theta[(k+kv+ku+1):(k+kv+ku+kdel)]
  mui <- as.vector( zdeli%*%delta )
 } else {
  mui <- rep(0, my.n)
 }
 # print(k+kv+ku+kdel)
 # print(kdel)
 # print(Ktheta)
 
 maxT_all <- maxT
 
 if(eff.time.invariant){
  Git <- rep(1, length(timevar))
 } else {
  if(model == "BC1992"){
   # print(theta)
   # print(c( Ktheta) )
   # print(length(theta))
   # print(theta[c( Ktheta) ])
   # eta <- theta[12 ]
   # cat("eta = ", theta[11 ],"\n")
   # print(eta)
   eta <- theta[c( Ktheta )]
   Git = exp(-eta*(timevar - maxT_all))
  } else if(model == "K1990modified"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = 1 + eta[1]*(timevar - maxT_all) + eta[2]*(timevar - maxT_all)^2
  } else if(model == "K1990"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = (1 + exp(eta[1]*timevar + eta[2]*timevar^2))^(-1)
  } else {
   stop("Unknown model\n")
  }
  # cat("eta\n")
  # print(eta)
 }
 
 # define indices
 m.end <- cumsum(t0)
 m.begin <- c(1,(m.end+1)[-my.n])

 # Log-likelihood
 llf = 0
 for(i in 1:length(ids)){
  # cat("i = ", i, "ll = ",llf,"\n")
  sample.i <- (m.begin[i]):(m.end[i])
  ei = eit[sample.i]
  sigv2i = sigv2it[sample.i]
  Gi = Git[sample.i]
  Ti = length(ei)
  lmdi = mui[i]/sigu2i[i]
  sstart = (1/sigu2i[i] + sum(Gi^2/sigv2i))
  
  llf = llf - 1/2*log(sstart) - 1/2*(mui[i]*lmdi + sum(ei^2/sigv2i)) + 1/2*((lmdi + s*sum(ei*Gi/sigv2i))^2)/sstart - Ti/2 * log(2*pi) - 1/2*sum(log(sigv2i)) - 1/2*log(sigu2i[i]) - pnorm(mui[i]/sqrt(sigu2i[i]), log.p = TRUE) + pnorm((lmdi + s*sum(ei*Gi/sigv2i))/sqrt(sstart), log.p = TRUE)
  
  # if (i == 2){
   # cat("\n")
   # print(Gi)
   # cat("sigu2i[i]\n")
   # print(sigu2i[i] )
   # print(sigv2i)
   # cat("1/2*log(sstart)\n")
   # print(1/2*log(sstart))
   # cat("1/2*(mui[i]*lmdi + sum(ei^2/sigv2i))\n")
   # print(1/2*(mui[i]*lmdi + sum(ei^2/sigv2i)))
   # cat("1/2*((lmdi + s*sum(ei*Gi/sigv2i))^2)/sstart\n")
   # print(1/2*((lmdi + s*sum(ei*Gi/sigv2i))^2)/sstart)
   # cat("1/2*sum(log(sigv2i))\n")
   # print(1/2*sum(log(sigv2i)))
   # cat("1/2*log(sigu2i[i])\n")
   # print(1/2*log(sigu2i[i]))
   # cat("pnorm(mui[i]/sqrt(sigu2i[i]), log.p = TRUE)\n")
   # print(pnorm(mui[i]/sqrt(sigu2i[i]), log.p = TRUE))
   # cat("pnorm((lmdi + s*sum(ei*Gi/sigv2i))/sqrt(sstart), log.p = TRUE)\n")
   # print(pnorm((lmdi + s*sum(ei*Gi/sigv2i))/sqrt(sstart), log.p = TRUE))
   # cat("\n")
  # }
 }
 
 return(llf)
}

.gr.hess.panel <- function(theta, prod, coef.fixed, my.n, Ktheta, k, kv, ku, kdel, yit, zvit, zui, xit, zdeli, model, timevar, maxT, ids, idvar, t0, eff.time.invariant = FALSE, mean.u.0i.zero = FALSE) {
 
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[(k+kv+1):(k+kv+ku)]
 
 eit <- as.vector( yit - xit%*%beta )
 sigu2i <- as.vector( exp(zui%*%gu) )
 sigv2it <- as.vector( exp(zvit%*%gv) )
 
 if(!mean.u.0i.zero){
  delta <- theta[(k+kv+ku+1):(k+kv+ku+kdel)]
  mui <- as.vector( zdeli%*%delta )
  lmdi = mui/sigu2i
 } else {
  mui <- lmdi <- rep(0, my.n)
 }
 
 
 maxT_all <- maxT
 
 if(eff.time.invariant){
  Git <- rep(1, length(timevar))
 } else {
  if(model == "BC1992"){
   # print(theta)
   # print(c( Ktheta) )
   # print(length(theta))
   # print(theta[c( Ktheta) ])
   # eta <- theta[12 ]
   # cat("eta = ", theta[11 ],"\n")
   # print(eta)
   eta <- theta[c( Ktheta )]
   Git = exp(-eta*(timevar - maxT_all))
  } else if(model == "K1990modified"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = 1 + eta[1]*(timevar - maxT_all) + eta[2]*(timevar - maxT_all)^2
  } else if(model == "K1990"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = (1 + exp(eta[1]*timevar + eta[2]*timevar^2))^(-1)
  } else {
   stop("Unknown model\n")
  }
  # cat("eta\n")
  # print(eta)
 }
 # Keta <- length(coef.fixed) - Ktheta
 # cat.print(Keta)
 # cat.print(length(coef.fixed))
 # cat.print(Ktheta)
 # # cat("Git\n")
 # print(Git)
 
 # define indices
 m.end <- cumsum(t0)
 m.begin <- c(1,(m.end+1)[-my.n])
 
 gb = ggu = ggv = gdel = geta = 0
 
 Hb = Hbgv = Hbgu = Hbdel = Hbeta = Hgvgu = Hgvdel = Hgv = Hdel = Hdelgu = Hdeleta = Hetagu = Hgu = Hetagv = Heta = 0;
 # Hb = Hbgv = Hbgu = Hbdel = Hbeta = Hgvgu = Hgvdel = Hgv = Hdel = Hdelgu = Hdeleta = Hetagu = Hgu = Hetagv = Heta = 0;
 # if(eff.time.invariant) Hbeta <- matrix(0, nrow = Keta, ncol = k)
 
 for(i in 1:length(ids)){

  sample.i <- (m.begin[i]):(m.end[i])
  ei = eit[sample.i]
  sigv2i = sigv2it[sample.i]
  Gi = Git[sample.i]
  
  sstart = (1/sigu2i[i] + sum(Gi^2/sigv2i))
  
  xi = xit[sample.i, , drop = FALSE]; 
  zvi = zvit[sample.i, , drop = FALSE];
  Ti = length(ei)
  timevari = timevar[sample.i]
  c1 = (1/sigu2i[i] + sum(Gi^2/sigv2i))^(-1);
  c2 = ei*Gi/sigv2i; 
  c2s = sum(c2);
  c3 = t(xi)%*%(Gi/sigv2i)
  a1i = mui[i]/sqrt(sigu2i[i]); 
  a2i = (lmdi[i] + s*sum(ei*Gi/sigv2i))*sqrt(c1)
  d1 = dnorm(a1i)/pnorm(a1i); 
  d2 = dnorm(a2i)/pnorm(a2i)
  d3 = (1 - d2*(a2i + d2)) # this is A
  d4 = d1*(a1i + d1)
  
  if(!eff.time.invariant){
   if(model == "BC1992"){
    tBC1 = (timevari - maxT_all[sample.i]);
   } else if (model == "K1990modified"){
    tBC2 = cbind((timevari - maxT_all[sample.i]), (timevari - maxT_all[sample.i])^2);
   } else if (model == "K1990"){
    Gk = (Gi^(-1)-1)
    tSK = cbind(timevari, timevari^2)
   } else {
    stop("Unknown model\n")
   }
  }
  
  # if (i == 2){
  #  cat("\n")
  #  cat("Gi\n")
  #  print(Gi)
  #  cat(" c1*c3%*%t(c3)*d3\n")
  #  print( c1*c3%*%t(c3)*d3)
  #  cat.print(c1)
  #  cat.print(c3)
  #  cat.print(d3)
  #  cat.print(a2i)
  #  cat.print(d2)
  #  cat("head(xi)\n")
  #  print(head(xi))
  #  cat("sigu2i[i]\n")
  #  print(sigu2i[i] )
  #  print(sigv2i)
  #  cat("1/2*log(sstart)\n")
  #  print(1/2*log(sstart))
  #  cat("1/2*(mui[i]*lmdi + sum(ei^2/sigv2i))\n")
  #  print(1/2*(mui[i]*lmdi + sum(ei^2/sigv2i)))
  #  cat("1/2*((lmdi + s*sum(ei*Gi/sigv2i))^2)/sstart\n")
  #  print(1/2*((lmdi + s*sum(ei*Gi/sigv2i))^2)/sstart)
  #  cat("1/2*sum(log(sigv2i))\n")
  #  print(1/2*sum(log(sigv2i)))
  #  cat("1/2*log(sigu2i[i])\n")
  #  print(1/2*log(sigu2i[i]))
  #  cat("pnorm(mui[i]/sqrt(sigu2i[i]), log.p = TRUE)\n")
  #  print(pnorm(mui[i]/sqrt(sigu2i[i]), log.p = TRUE))
  #  cat("pnorm((lmdi + s*sum(ei*Gi/sigv2i))/sqrt(sstart), log.p = TRUE)\n")
  #  print(pnorm((lmdi + s*sum(ei*Gi/sigv2i))/sqrt(sstart), log.p = TRUE))
  #  cat("\n")
  # }
  
  # Gradient
  
  gb = gb + t(xi)%*%((ei/sigv2i)- s*(Gi/sigv2i)*(lmdi[i]*c1 + s*c2s*c1 + sqrt(c1)*d2))
  ggu = ggu + zui[i,,drop=FALSE]/sigu2i[i]*(0.5*c1 + mui[i]*(0.5*mui[i] + 0.5*d1*sqrt(sigu2i[i]) - sqrt(c1)*d2) + a2i*(c1*a2i/2 - sqrt(c1)*mui[i] + c1*d2/2)) - zui[i,,drop=FALSE]*0.5
  ggv = ggv + t(zvi)%*%(0.5*c1*(Gi^2/sigv2i) + 0.5*(ei^2/sigv2i) - 0.5*rep.int(1, Ti) + a2i*sqrt(c1)*((a2i + d2)*0.5*sqrt(c1)*(Gi^2/sigv2i) - s*(ei*Gi/sigv2i)) - s*sqrt(c1)*d2*(ei*Gi/sigv2i))
  gdel = gdel + zdeli[i,,drop=FALSE]/sigu2i[i] *(mui[i]*(c1/sigu2i[i] - 1) + s*c2s*c1 - d1 *sqrt(sigu2i[i]) + d2*sqrt(c1))
  
  # Hessian
  
  Hb = Hb - t(xi)%*%sweep(xi, 1, sigv2i, "/") + c1*c3%*%t(c3)*d3
  
  Hbgu = Hbgu + s*c3%*%zui[i,,drop=FALSE]*c1*(lmdi[i]*d3 - c1/sigu2i[i]*(lmdi[i] + s*sum(c2)) + sqrt(c1)/(2*sigu2i[i])*d2*(a2i*(a2i + d2) - 1))
  Hbdel = Hbdel - s*c3%*%(zdeli[i,]/sigu2i[i]*d3*c1)
  Hbgv = Hbgv + t(xi)%*%(-sweep(zvi, 1,ei/sigv2i, "*") + s*sweep(zvi, 1, Gi/sigv2i, "*")*(sqrt(c1)*(a2i + d2))) - s*(t(xi)%*%(Gi/sigv2i))%*%t(t(zvi)%*%(-s*c1*c2*d3 + (c1^1.5*(Gi^2/sigv2i))*(a2i + 0.5*d2*(1 - a2i*(a2i + d2)))))
  Hgvgu = Hgvgu + t(zvi)%*%(c1*s*c2*(mui[i]*d3 - sqrt(c1)*a2i/2*(1+d3)) + c1^1.5*a2i*0.5*(Gi^2/sigv2i)*(sqrt(c1)*a2i*0.5*(3+d3) - mui[i]*(1+d3)) + 0.5*c1^1.5*(Gi^2/sigv2i)*(sqrt(c1) + d2*(1.5*sqrt(c1)*a2i - mui[i])) - 0.5*c1^1.5*s*d2*c2)%*%zui[i,,drop=FALSE]/sigu2i[i]
  Hgvdel = Hgvdel + t(zvi)%*%(-s*c1*c2*d3 + (Gi^2/sigv2i)*(0.5*c1^1.5*(a2i*(1+d3) + d2)))%*%zdeli[i,]/sigu2i[i]
  Hgv = Hgv +  t(zvi)%*%sweep(zvi,1,Gi^2/sigv2i, "*")*0.5*c1*(-d2*a2i - a2i^2 -1) + t(zvi)%*%(Gi^2/sigv2i)%*%t(t(zvi)%*%(Gi^2/sigv2i))*0.5*c1^2*(1 + 0.5*a2i^2*(3+d3) + d2*1.5*a2i) - 0.5*t(zvi)%*%sweep(zvi,1,ei^2/sigv2i, "*") + t(zvi)%*%sweep(zvi,1,c2, "*")*s*sqrt(c1)*(d2 + a2i) - (d2+a2i*(1+d3))*0.5*c1^1.5*s*(t(zvi)%*%(c2)%*%t(t(zvi)%*%(Gi^2/sigv2i)) + t(zvi)%*%(Gi^2/sigv2i)%*%t(t(zvi)%*%(c2))) + t(zvi)%*%(c2)%*%t(t(zvi)%*%(c2))*c1*d3
  Hdel = Hdel + t(zdeli[i,,drop=FALSE]/sigu2i[i] *( - 1 + d3*c1/sigu2i[i] + d4))%*%zdeli[i,,drop=FALSE]
  Hdelgu =  Hdelgu + t(zdeli[i,,drop=FALSE])%*%zui[i,,drop=FALSE]/sigu2i[i]*(c1^1.5/sigu2i[i]*0.5*a2i*(1+d3) - mui[i]*c1/sigu2i[i]*(1+d3) + d2*sqrt(c1)*(0.5*c1/sigu2i[i] - 1) - 0.5*d1*(mui[i]*(a1i + d1) - sqrt(sigu2i[i])) + mui[i] - s*c1*sum(c2))
  Hgu = Hgu + (c1/sigu2i[i]*(0.5*c1 + (sqrt(c1)*a2i - mui[i])^2) - 0.5*(c1 + (sqrt(c1)*a2i - mui[i])^2) + d1*mui[i]/4*(mui[i]*(a1i + d1) - sqrt(sigu2i[i])) + d2*(-d2*2*c1/sigu2i[i]*(mui[i]/sqrt(2) - sqrt(c1/8)*a2i)^2 - c1*a2i/sigu2i[i]*(mui[i] - sqrt(c1)*a2i/2)^2 + sqrt(c1)*mui[i]*(1 - c1/sigu2i[i]) - c1*a2i/2*(1 - 1.5*c1/sigu2i[i])))*t(zui[i,,drop=FALSE])%*%zui[i,,drop=FALSE]/sigu2i[i]
  
  if(!eff.time.invariant){
   if(model == "BC1992"){
    geta = geta + c1*sum(Gi^2*tBC1/sigv2i) - s*sqrt(c1)*d2*sum(ei*Gi*tBC1/sigv2i) + sqrt(c1)*a2i*(sqrt(c1)*(a2i + d2)*sum(Gi^2*tBC1/sigv2i) - s*sum(ei*Gi*tBC1/sigv2i))
    Hbeta = Hbeta + s*t(xi)%*%(Gi*tBC1/sigv2i)*sqrt(c1)*(a2i + d2) - s*t(xi)%*%(Gi/sigv2i)*c1*(-s*sum(c2*tBC1)*d3 + sqrt(c1)*a2i*(1 + d3)*sum(Gi^2*tBC1/sigv2i) + sqrt(c1)*d2*sum(Gi^2*tBC1/sigv2i))
    Hdeleta = Hdeleta + t((sum(Gi^2*tBC1/sigv2i)*c1^1.5*(a2i*(1+d3) + d2)-s*c1*sum(c2*tBC1)*d3)*zdeli[i,,drop=FALSE]/sigu2i[i])
    Hetagu = Hetagu + t((sum(Gi^2*tBC1/sigv2i)*c1^1.5/sigu2i[i]*(sqrt(c1) -a2i*mui[i]*(1+d3) + sqrt(c1)*a2i^2*0.5*(3+d3) - d2*(mui[i] - 1.5*a2i*sqrt(c1))) + s*c1/sigu2i[i]*sum(c2*tBC1)*(mui[i]*d3 - d2*0.5*sqrt(c1) - 0.5*sqrt(c1)*a2i*(1+d3)))*zui[i,,drop=FALSE])
    Hetagv = Hetagv - t(zvi)%*%(Gi^2*tBC1/sigv2i)*c1*(1 + a2i^2 + a2i*d2) + sum(Gi^2*tBC1/sigv2i)*t(zvi)%*%(Gi^2/sigv2i)*c1^2*(1 + 0.5*a2i^2*(3 + d3) + 1.5*a2i*d2) + t(zvi)%*%(c2*tBC1)*s*sqrt(c1)*(d2 + a2i) + sum(c2*tBC1)*t(zvi)%*%(c2)*c1*d3 - (a2i*(1+d3) + d2)*s*c1^1.5*(0.5*sum(c2*tBC1)*t(zvi)%*%(Gi^2/sigv2i) + sum(Gi^2*tBC1/sigv2i)*t(zvi)%*%(c2))
    Heta = Heta - 2*c1*sum(Gi^2*tBC1^2/sigv2i)*(1 + a2i^2 + d2*a2i) + c1^2*sum(Gi^2*tBC1/sigv2i)^2*(2 + a2i^2*(3 + d3) + 3*d2*a2i) + sum(c2*tBC1^2)*s*sqrt(c1)*(d2 + a2i) - 2*s*c1^1.5*sum(c2*tBC1)*sum(Gi^2*tBC1/sigv2i)*(d2 + a2i*(1 + d3)) + sum(c2*tBC1)^2*c1*d3
    
   } else if(model == "K1990modified"){
    geta =  geta - c1*t(Gi/sigv2i)%*%tBC2 + 
     s*sqrt(c1)*d2*t(ei/sigv2i)%*%tBC2 - 
     sqrt(c1)*a2i*(sqrt(c1)*(a2i + d2)*t(Gi/sigv2i)%*%tBC2 - s*t(ei/sigv2i)%*%tBC2)
    Hbeta = Hbeta - s*t(xi)%*%sweep(tBC2, 1, sigv2i, "/")*sqrt(c1)*(a2i + d2) - s*(t(xi)%*%(Gi/sigv2i))%*%t(t(tBC2)%*%(c1*(s*(ei/sigv2i)*d3 - sqrt(c1)*a2i*(1 + d3)*(Gi/sigv2i) - sqrt(c1)*d2*(Gi/sigv2i))))
    Hdeleta = Hdeleta + t(zdeli[i,,drop=FALSE])%*%(t(t(tBC2)%*%(s*c1*(ei/sigv2i)*d3 - c1^1.5*(a2i*(1+d3) + d2)*(Gi/sigv2i)))/sigu2i[i])
    Hetagu = Hetagu + t(zui[i,,drop=FALSE])%*%(- t(Gi/sigv2i)%*%tBC2*c1^1.5/sigu2i[i]*(sqrt(c1) -a2i*mui[i]*(1+d3) + sqrt(c1)*a2i^2*0.5*(3+d3) - d2*(mui[i] - 1.5*a2i*sqrt(c1))) - t(ei/sigv2i)%*%tBC2*s*c1/sigu2i[i]*(mui[i]*d3 - d2*0.5*sqrt(c1) - 0.5*sqrt(c1)*a2i*(1+d3)))
    Hetagv = Hetagv + t(zvi)%*%sweep(tBC2, 1, Gi/sigv2i, "*")*c1*(1 + a2i^2 + a2i*d2) - (t(zvi)%*%(Gi^2/sigv2i))%*%(t(Gi/sigv2i)%*%tBC2)*c1^2*(1 + 0.5*a2i^2*(3 + d3) + 1.5*a2i*d2) - t(zvi)%*%sweep(tBC2, 1, ei/sigv2i, "*")*s*sqrt(c1)*(d2 + a2i) - (t(zvi)%*%(c2))%*%(t(ei/sigv2i)%*%tBC2)*c1*d3 + (a2i*(1+d3) + d2)*s*c1^1.5*(0.5*(t(zvi)%*%(Gi^2/sigv2i))%*%(t(ei/sigv2i)%*%tBC2) + (t(zvi)%*%(c2))%*%(t(Gi/sigv2i)%*%tBC2))
    Heta = Heta - c1*(1 + a2i^2 + a2i*d2)*t(tBC2)%*%sweep(tBC2, 1, sigv2i, "/") + c1^2*(2 + a2i^2*(3 + d3) + d2*a2i*3)*t(t(Gi/sigv2i)%*%tBC2)%*%(t(Gi/sigv2i)%*%tBC2) - c1^1.5*s*(d2+a2i*(1+d3))*t(t(ei/sigv2i)%*%tBC2)%*%(t(Gi/sigv2i)%*%tBC2) + c1*d3*t(t(ei/sigv2i)%*%tBC2)%*%(t(ei/sigv2i)%*%tBC2) - c1^1.5*s*(d2 + a2i*(1+d3))*t(t(Gi/sigv2i)%*%tBC2)%*%(t(ei/sigv2i)%*%tBC2)
    
   } else if(model == "K1990"){
    geta =  geta + c1*t(Gi^3*(Gi^(-1)-1)/sigv2i)%*%tSK -  s*sqrt(c1)*d2*t(ei*Gi^2*(Gi^(-1)-1)/sigv2i)%*%tSK + sqrt(c1)*a2i*(sqrt(c1)*(a2i + d2)*t(Gi^3*(Gi^(-1)-1)/sigv2i)%*%tSK - s*t(ei*Gi^2*(Gi^(-1)-1)/sigv2i)%*%tSK)
    Hbeta = Hbeta +  s*t(xi)%*%sweep(tSK, 1, Gi^2*Gk/sigv2i, "*")*sqrt(c1)*(a2i + d2) - s*(t(xi)%*%(Gi/sigv2i))%*%t(t(tSK)%*%(c1*(-s*(ei*Gi^2*Gk/sigv2i)*d3 + sqrt(c1)*a2i*(1 + d3)*(Gi^3*Gk/sigv2i) + sqrt(c1)*d2*(Gi^3*Gk/sigv2i))))
    
    Hdeleta = Hdeleta + t(zdeli[i,,drop=FALSE])%*%(s*t(t(tSK)%*%(-s*c1*(c2*Gi*Gk)*d3 + c1^1.5*(a2i*(1+d3) + d2)*(Gi^3*Gk/sigv2i)))/sigu2i[i])
    
    Hetagu = Hetagu + t(zui[i,,drop=FALSE])%*%(t(Gi^3*Gk/sigv2i)%*%tSK*c1^1.5/sigu2i[i]*(sqrt(c1) -a2i*mui[i]*(1+d3) + sqrt(c1)*a2i^2*0.5*(3+d3) - d2*(mui[i] - 1.5*a2i*sqrt(c1))) + t(c2*Gi*Gk)%*%tSK*s*c1/sigu2i[i]*(mui[i]*d3 - d2*0.5*sqrt(c1) - 0.5*sqrt(c1)*a2i*(1+d3)))
    
    Hetagv = Hetagv - t(zvi)%*%sweep(tSK, 1, Gi^3*Gk/sigv2i, "*")*c1*(1 + a2i^2 + a2i*d2) + (t(zvi)%*%(Gi^2/sigv2i))%*%(t(Gi^3*Gk/sigv2i)%*%tSK)*c1^2*(1 + 0.5*a2i^2*(3 + d3) + 1.5*a2i*d2) + t(zvi)%*%sweep(tSK, 1, c2*Gi*Gk, "*")*s*sqrt(c1)*(d2 + a2i) + (t(zvi)%*%(c2))%*%(t(c2*Gi*Gk)%*%tSK)*c1*d3 - (a2i*(1+d3) + d2)*s*c1^1.5*(0.5*(t(zvi)%*%(Gi^2/sigv2i))%*%(t(c2*Gi*Gk)%*%tSK) + (t(zvi)%*%(c2))%*%(t(Gi^3*Gk/sigv2i)%*%tSK))
    
    Heta = Heta + c1*(1 + a2i^2 + a2i*d2)*t(tSK)%*%sweep(tSK, 1, Gi^3*Gk*(1- 3*Gi*Gk)/sigv2i, "*") + c1^2*(2 + a2i^2*(3 + d3) + d2*a2i*3)*t(t(Gi^3*Gk/sigv2i)%*%tSK)%*%(t(Gi^3*Gk/sigv2i)%*%tSK) - s*sqrt(c1)*(a2i + d2)*t(tSK)%*%sweep(tSK, 1, ei*Gi^2*Gk*(1 - 2*Gi*Gk)/sigv2i, "*") - c1^1.5*s*(d2 + a2i*(1 + d3))*(t(t(ei*Gi^2*Gk/sigv2i)%*%tSK)%*%(t(Gi^3*Gk/sigv2i)%*%tSK) + t(t(Gi^3*Gk/sigv2i)%*%tSK)%*%(t(ei*Gi^2*Gk/sigv2i)%*%tSK)) + c1*d3*t(t(ei*Gi^2*Gk/sigv2i)%*%tSK)%*%(t(ei*Gi^2*Gk/sigv2i)%*%tSK)
    # print(geta)
   }
   else {
    stop("Unknown model\n")
   }
  }
  
  # if(i == 2){
  #  cat("Hb\n")
  #  print(Hb)
  # }

 }
 
 
 
 if(mean.u.0i.zero){
  
  if(eff.time.invariant){
   H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu)),
              rbind(Hbgv, 	Hgv,		t(Hgvgu)),
              rbind(Hbgu, 	Hgvgu,		Hgu))
  } else {
   H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),	t(Hbeta)),
              rbind(Hbgv, 	Hgv,		t(Hgvgu),	t(Hetagv)),
              rbind(Hbgu, 	Hgvgu,		Hgu,		 	t(Hetagu)),
              rbind(Hbeta,	Hetagv,		Hetagu,		Heta ))
  }
 } else {
  
  if(eff.time.invariant){
   H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),		t(Hbdel)),
              rbind(Hbgv, 	Hgv,		t(Hgvgu),		t(Hgvdel)),
              rbind(Hbgu, 	Hgvgu,		Hgu,			Hdelgu),
              rbind(Hbdel,	Hgvdel,		t(Hdelgu),		Hdel))
   
  } else {
   H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),		t(Hbdel),	t(Hbeta)),
              rbind(Hbgv, 	Hgv,		t(Hgvgu),		t(Hgvdel),	t(Hetagv)),
              rbind(Hbgu, 	Hgvgu,		Hgu,			Hdelgu, 	t(Hetagu)),
              rbind(Hbdel,	Hgvdel,		t(Hdelgu),		Hdel,		t(Hdeleta) ),
              rbind(Hbeta,	Hetagv,		Hetagu,			Hdeleta,	Heta ))
  }
  
 }
 
 # if(eff.time.invariant){
 #  grad <- rbind(gb, as.matrix(ggv),t(as.matrix(ggu)), t(as.matrix(gdel)))
 #  H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),		t(Hbdel)),
 #             rbind(Hbgv, 	Hgv,		t(Hgvgu),		t(Hgvdel)),
 #             rbind(Hbgu, 	Hgvgu,		Hgu,			Hdelgu),
 #             rbind(Hbdel,	Hgvdel,		t(Hdelgu),		Hdel))
 # } else {
 #  grad <- rbind(gb, as.matrix(ggv),t(as.matrix(ggu)), t(as.matrix(gdel)), t(geta))
 #  H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),		t(Hbdel),	t(Hbeta)),
 #             rbind(Hbgv, 	Hgv,		t(Hgvgu),		t(Hgvdel),	t(Hetagv)),
 #             rbind(Hbgu, 	Hgvgu,		Hgu,			Hdelgu, 	t(Hetagu)),
 #             rbind(Hbdel,	Hgvdel,		t(Hdelgu),		Hdel,		t(Hdeleta) ),
 #             rbind(Hbeta,	Hetagv,		Hetagu,			Hdeleta,	Heta ))
 # }
 
 grad <- rbind(gb, as.matrix(ggv),t(as.matrix(ggu)), t(as.matrix(gdel)), t(geta))
 grad <- as.vector(grad)
 
 # cat.print(Hbeta)
 # 
 # H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),		t(Hbdel),	t(Hbeta)), 
 #            rbind(Hbgv, 	Hgv,		t(Hgvgu),		t(Hgvdel),	t(Hetagv)),
 #            rbind(Hbgu, 	Hgvgu,		Hgu,			Hdelgu, 	t(Hetagu)), 
 #            rbind(Hbdel,	Hgvdel,		t(Hdelgu),		Hdel,		t(Hdeleta) ),
 #            rbind(Hbeta,	Hetagv,		Hetagu,			Hdeleta,	Heta ))
 
 
 grad <- grad[!coef.fixed]
 # print(H)
 # H <- H[!coef.fixed,!coef.fixed]
 # print(H)
 
 # if(mean.u.0i.zero){
 #  grad <- grad[-(k+kv+ku+kdel)]
 #  H <- H[-(k+kv+ku+kdel),-(k+kv+ku+kdel)]
 # } 
 
 # print(length(grad))
 # print(dim(H))
 # print(length(theta))
 
 
 names(grad) <- colnames(H) <- rownames(H) <- names(theta)
 return(list(grad = grad, hessian1 = H))
 
}

# Gradient

.gr.panel <- function(theta, prod, coef.fixed, my.n, Ktheta, k, kv, ku, kdel, yit, zvit, zui, xit, zdeli, model, timevar, maxT, ids, idvar, t0, eff.time.invariant = FALSE, mean.u.0i.zero = FALSE) {
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[(k+kv+1):(k+kv+ku)]
 
 eit <- as.vector( yit - xit%*%beta )
 sigu2i <- as.vector( exp(zui%*%gu) )
 sigv2it <- as.vector( exp(zvit%*%gv) )
 
 if(!mean.u.0i.zero){
  delta <- theta[(k+kv+ku+1):(k+kv+ku+kdel)]
  mui <- as.vector( zdeli%*%delta )
  lmdi = mui/sigu2i
 } else {
  mui <- lmdi <- rep(0, my.n)
 }
 
 maxT_all <- maxT
 
 if(eff.time.invariant){
  Git <- rep(1, length(timevar))
 } else {
  if(model == "BC1992"){
   # print(theta)
   # print(c( Ktheta) )
   # print(length(theta))
   # print(theta[c( Ktheta) ])
   # eta <- theta[12 ]
   # cat("eta = ", theta[11 ],"\n")
   # print(eta)
   eta <- theta[c( Ktheta )]
   Git = exp(-eta*(timevar - maxT_all))
  } else if(model == "K1990modified"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = 1 + eta[1]*(timevar - maxT_all) + eta[2]*(timevar - maxT_all)^2
  } else if(model == "K1990"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = (1 + exp(eta[1]*timevar + eta[2]*timevar^2))^(-1)
  } else {
   stop("Unknown model\n")
  }
  # cat("eta\n")
  # print(eta)
 }
 
 # define indices
 m.end <- cumsum(t0)
 m.begin <- c(1,(m.end+1)[-my.n])
 
 # Gradient
 gb = 0; ggu = 0; ggv = 0; gdel = 0; geta = 0
 for(i in 1:length(ids)){
  # maxT <- t0[i]
  sample.i <- (m.begin[i]):(m.end[i])
  ei = eit[sample.i]
  sigv2i = sigv2it[sample.i]
  Gi = Git[sample.i]
  
  xi = xit[sample.i, , drop = FALSE]; 
  zvi = zvit[sample.i, , drop = FALSE];
  Ti = length(ei)
  timevari = timevar[sample.i]
  c1 = (1/sigu2i[i] + sum(Gi^2/sigv2i))^(-1); 
  c2 = sum(ei*Gi/sigv2i);
  a1i = mui[i]/sqrt(sigu2i[i]); 
  a2i = (lmdi[i] + s*sum(ei*Gi/sigv2i))*sqrt(c1)
  d1 = dnorm(a1i)/pnorm(a1i); 
  d2 = dnorm(a2i)/pnorm(a2i)
  
  if(!eff.time.invariant){
   if(model == "BC1992"){
    tBC1 = (timevari - maxT_all[sample.i]);
   } else if (model == "K1990modified"){
    tBC2 = cbind((timevari - maxT_all[sample.i]), (timevari - maxT_all[sample.i])^2);
   } else if (model == "K1990"){
    tSK = cbind(timevari, timevari^2)
   } else {
    stop("Unknown model\n")
   }
  }
  
  gb = gb + t(xi)%*%((ei/sigv2i)- s*(Gi/sigv2i)*(lmdi[i]*c1 + s*c2*c1 + sqrt(c1)*d2))
  ggu = ggu + zui[i,,drop=FALSE]/sigu2i[i]*(0.5*c1 + mui[i]*(0.5*mui[i] + 0.5*d1*sqrt(sigu2i[i]) - sqrt(c1)*d2) + a2i*(c1*a2i/2 - sqrt(c1)*mui[i] + c1*d2/2)) - zui[i,,drop=FALSE]*0.5
  ggv = ggv + t(zvi)%*%(0.5*c1*(Gi^2/sigv2i) + 0.5*(ei^2/sigv2i) - 0.5*rep.int(1, Ti) + a2i*sqrt(c1)*((a2i + d2)*0.5*sqrt(c1)*(Gi^2/sigv2i) - s*(ei*Gi/sigv2i)) - s*sqrt(c1)*d2*(ei*Gi/sigv2i))
  gdel = gdel + zdeli[i,,drop=FALSE]/sigu2i[i] *(mui[i]*(c1/sigu2i[i] - 1) + s*c2*c1 - d1 *sqrt(sigu2i[i]) + d2*sqrt(c1))
  
  if(!eff.time.invariant){
   if(model == "BC1992"){
    geta = geta + c1*sum(Gi^2*tBC1/sigv2i) - s*sqrt(c1)*d2*sum(ei*Gi*tBC1/sigv2i) + sqrt(c1)*a2i*(sqrt(c1)*(a2i + d2)*sum(Gi^2*tBC1/sigv2i) - s*sum(ei*Gi*tBC1/sigv2i))
   } else if(model == "K1990modified"){
    geta =  geta - c1*t(Gi/sigv2i)%*%tBC2 + 
     s*sqrt(c1)*d2*t(ei/sigv2i)%*%tBC2 - 
     sqrt(c1)*a2i*(sqrt(c1)*(a2i + d2)*t(Gi/sigv2i)%*%tBC2 - s*t(ei/sigv2i)%*%tBC2)
   } else if(model == "K1990"){
    geta =  geta + c1*t(Gi^3*(Gi^(-1)-1)/sigv2i)%*%tSK -  s*sqrt(c1)*d2*t(ei*Gi^2*(Gi^(-1)-1)/sigv2i)%*%tSK + sqrt(c1)*a2i*(sqrt(c1)*(a2i + d2)*t(Gi^3*(Gi^(-1)-1)/sigv2i)%*%tSK - s*t(ei*Gi^2*(Gi^(-1)-1)/sigv2i)%*%tSK)
    # print(geta)
   }
   else {
    stop("Unknown model\n")
   }
  }
 }
 grad <- rbind(gb, as.matrix(ggv),t(as.matrix(ggu)), t(as.matrix(gdel)), t(geta))
 grad <- grad[!coef.fixed]
 grad <- as.vector(grad)
 names(grad) <- names(theta)
 return(grad)
}

# .gr.panel <- function(theta, prod, coef.fixed, my.n, Ktheta, k, kv, ku, kdel, yit, zvit, zui, xit, zdeli, model, timevar, maxT, ids, idvar, t0, eff.time.invariant = FALSE, mean.u.0i.zero = FALSE){
#  numDeriv::grad(func = .ll.panel, x = theta, prod = prod, my.n = my.n, k = k, kv = kv, ku = ku, kdel = kdel, yit = yit, zvit = zvit, zui = zui, xit = xit, zdeli = zdeli, eff.time.invariant = eff.time.invariant, model = model, timevar = timevar, maxT = maxT, ids = ids, idvar = idvar, t0 = t0)
# }

# Hessian

.hess.panel <- function(theta, prod, coef.fixed, my.n, Ktheta, k, kv, ku, kdel, yit, zvit, zui, xit, zdeli, model, timevar, maxT, ids, idvar, t0, eff.time.invariant = FALSE, mean.u.0i.zero = FALSE) {
 s <- ifelse(prod, -1, 1)
 beta <- theta[1:k]
 gv <- theta[(k+1):(k+kv)]
 gu <- theta[(k+kv+1):(k+kv+ku)]
 
 eit <- as.vector( yit - xit%*%beta )
 sigu2i <- as.vector( exp(zui%*%gu) )
 sigv2it <- as.vector( exp(zvit%*%gv) )
 
 if(!mean.u.0i.zero){
  delta <- theta[(k+kv+ku+1):(k+kv+ku+kdel)]
  mui <- as.vector( zdeli%*%delta )
  lmdi = mui/sigu2i
 } else {
  mui <- lmdi <- rep(0, my.n)
 }

 maxT_all <- maxT
 
 if(eff.time.invariant){
  Git <- rep(1, length(timevar))
 } else {
  if(model == "BC1992"){
   # print(theta)
   # print(c( Ktheta) )
   # print(length(theta))
   # print(theta[c( Ktheta) ])
   # eta <- theta[12 ]
   # cat("eta = ", theta[11 ],"\n")
   # print(eta)
   eta <- theta[c( Ktheta )]
   Git = exp(-eta*(timevar - maxT_all))
  } else if(model == "K1990modified"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = 1 + eta[1]*(timevar - maxT_all) + eta[2]*(timevar - maxT_all)^2
  } else if(model == "K1990"){
   eta <- theta[c( (Ktheta-1) : (Ktheta) )]
   Git = (1 + exp(eta[1]*timevar + eta[2]*timevar^2))^(-1)
  } else {
   stop("Unknown model\n")
  }
  # cat("eta\n")
  # print(eta)
 }
 
 # define indices
 m.end <- cumsum(t0)
 m.begin <- c(1,(m.end+1)[-my.n])
 
 # Hessian
 Hb = Hbgv = Hbgu = Hbdel = Hbeta = Hgvgu = Hgvdel = Hgv = Hdel = Hdelgu = Hdeleta = Hetagu = Hgu = Hetagv = Heta = 0;
 
 for(i in 1:length(ids)){
  sample.i <- (m.begin[i]):(m.end[i])
  # maxT <- maxT_all[sample.i]
  xi = xit[sample.i, , drop=FALSE]; 
  zvi = zvit[sample.i, , drop=FALSE];
  ei = eit[sample.i]; 
  sigv2i = sigv2it[sample.i]
  Gi = Git[sample.i]; 
  Ti = length(ei); 
  timevari = timevar[sample.i];
  
  if(!eff.time.invariant){
   if(model == "BC1992"){
    tBC1 = (timevari - maxT_all[sample.i]);
   } else if (model == "K1990modified"){
    tBC2 = cbind((timevari - maxT_all[sample.i]), (timevari - maxT_all[sample.i])^2);
   } else if (model == "K1990"){
    Gk = (Gi^(-1)-1)
    tSK = cbind(timevari, timevari^2)
   } else {
    stop("Unknown model\n")
   }
  }

  c1 = (1/sigu2i[i] + sum(Gi^2/sigv2i))^(-1); 
  c2 = ei*Gi/sigv2i; 
  c3 = t(xi)%*%(Gi/sigv2i)
  a1i = mui[i]/sqrt(sigu2i[i]); 
  a2i = (lmdi[i] + s*sum(ei*Gi/sigv2i))*sqrt(c1)
  d1 = dnorm(a1i)/pnorm(a1i); 
  d2 = dnorm(a2i)/pnorm(a2i); 
  d3 = (1 - d2*(a2i + d2)); # this is A
  d4 = d1*(a1i + d1)
  
  Hb = Hb - t(xi)%*%sweep(xi, 1, sigv2i, "/") + c1*c3%*%t(c3)*d3
  
  Hbgu = Hbgu + s*c3%*%zui[i,,drop=FALSE]*c1*(lmdi[i]*d3 - c1/sigu2i[i]*(lmdi[i] + s*sum(c2)) + sqrt(c1)/(2*sigu2i[i])*d2*(a2i*(a2i + d2) - 1))
  Hbdel = Hbdel - s*c3%*%(zdeli[i,]/sigu2i[i]*d3*c1)
  Hbgv = Hbgv + t(xi)%*%(-sweep(zvi, 1,ei/sigv2i, "*") + s*sweep(zvi, 1, Gi/sigv2i, "*")*(sqrt(c1)*(a2i + d2))) - s*(t(xi)%*%(Gi/sigv2i))%*%t(t(zvi)%*%(-s*c1*c2*d3 + (c1^1.5*(Gi^2/sigv2i))*(a2i + 0.5*d2*(1 - a2i*(a2i + d2)))))
  Hgvgu = Hgvgu + t(zvi)%*%(c1*s*c2*(mui[i]*d3 - sqrt(c1)*a2i/2*(1+d3)) + c1^1.5*a2i*0.5*(Gi^2/sigv2i)*(sqrt(c1)*a2i*0.5*(3+d3) - mui[i]*(1+d3)) + 0.5*c1^1.5*(Gi^2/sigv2i)*(sqrt(c1) + d2*(1.5*sqrt(c1)*a2i - mui[i])) - 0.5*c1^1.5*s*d2*c2)%*%zui[i,,drop=FALSE]/sigu2i[i]
  Hgvdel = Hgvdel + t(zvi)%*%(-s*c1*c2*d3 + (Gi^2/sigv2i)*(0.5*c1^1.5*(a2i*(1+d3) + d2)))%*%zdeli[i,]/sigu2i[i]
  Hgv = Hgv +  t(zvi)%*%sweep(zvi,1,Gi^2/sigv2i, "*")*0.5*c1*(-d2*a2i - a2i^2 -1) + t(zvi)%*%(Gi^2/sigv2i)%*%t(t(zvi)%*%(Gi^2/sigv2i))*0.5*c1^2*(1 + 0.5*a2i^2*(3+d3) + d2*1.5*a2i) - 0.5*t(zvi)%*%sweep(zvi,1,ei^2/sigv2i, "*") + t(zvi)%*%sweep(zvi,1,c2, "*")*s*sqrt(c1)*(d2 + a2i) - (d2+a2i*(1+d3))*0.5*c1^1.5*s*(t(zvi)%*%(c2)%*%t(t(zvi)%*%(Gi^2/sigv2i)) + t(zvi)%*%(Gi^2/sigv2i)%*%t(t(zvi)%*%(c2))) + t(zvi)%*%(c2)%*%t(t(zvi)%*%(c2))*c1*d3
  Hdel = Hdel + t(zdeli[i,,drop=FALSE]/sigu2i[i] *( - 1 + d3*c1/sigu2i[i] + d4))%*%zdeli[i,,drop=FALSE]
  Hdelgu =  Hdelgu + t(zdeli[i,,drop=FALSE])%*%zui[i,,drop=FALSE]/sigu2i[i]*(c1^1.5/sigu2i[i]*0.5*a2i*(1+d3) - mui[i]*c1/sigu2i[i]*(1+d3) + d2*sqrt(c1)*(0.5*c1/sigu2i[i] - 1) - 0.5*d1*(mui[i]*(a1i + d1) - sqrt(sigu2i[i])) + mui[i] - s*c1*sum(c2))
  Hgu = Hgu + (c1/sigu2i[i]*(0.5*c1 + (sqrt(c1)*a2i - mui[i])^2) - 0.5*(c1 + (sqrt(c1)*a2i - mui[i])^2) + d1*mui[i]/4*(mui[i]*(a1i + d1) - sqrt(sigu2i[i])) + d2*(-d2*2*c1/sigu2i[i]*(mui[i]/sqrt(2) - sqrt(c1/8)*a2i)^2 - c1*a2i/sigu2i[i]*(mui[i] - sqrt(c1)*a2i/2)^2 + sqrt(c1)*mui[i]*(1 - c1/sigu2i[i]) - c1*a2i/2*(1 - 1.5*c1/sigu2i[i])))*t(zui[i,,drop=FALSE])%*%zui[i,,drop=FALSE]/sigu2i[i]
  
  if(!eff.time.invariant){
   if(model == "BC1992"){
    Hbeta = Hbeta + s*t(xi)%*%(Gi*tBC1/sigv2i)*sqrt(c1)*(a2i + d2) - s*t(xi)%*%(Gi/sigv2i)*c1*(-s*sum(c2*tBC1)*d3 + sqrt(c1)*a2i*(1 + d3)*sum(Gi^2*tBC1/sigv2i) + sqrt(c1)*d2*sum(Gi^2*tBC1/sigv2i))
    Hdeleta = Hdeleta + t((sum(Gi^2*tBC1/sigv2i)*c1^1.5*(a2i*(1+d3) + d2)-s*c1*sum(c2*tBC1)*d3)*zdeli[i,,drop=FALSE]/sigu2i[i])
    Hetagu = Hetagu + t((sum(Gi^2*tBC1/sigv2i)*c1^1.5/sigu2i[i]*(sqrt(c1) -a2i*mui[i]*(1+d3) + sqrt(c1)*a2i^2*0.5*(3+d3) - d2*(mui[i] - 1.5*a2i*sqrt(c1))) + s*c1/sigu2i[i]*sum(c2*tBC1)*(mui[i]*d3 - d2*0.5*sqrt(c1) - 0.5*sqrt(c1)*a2i*(1+d3)))*zui[i,,drop=FALSE])
    Hetagv = Hetagv - t(zvi)%*%(Gi^2*tBC1/sigv2i)*c1*(1 + a2i^2 + a2i*d2) + sum(Gi^2*tBC1/sigv2i)*t(zvi)%*%(Gi^2/sigv2i)*c1^2*(1 + 0.5*a2i^2*(3 + d3) + 1.5*a2i*d2) + t(zvi)%*%(c2*tBC1)*s*sqrt(c1)*(d2 + a2i) + sum(c2*tBC1)*t(zvi)%*%(c2)*c1*d3 - (a2i*(1+d3) + d2)*s*c1^1.5*(0.5*sum(c2*tBC1)*t(zvi)%*%(Gi^2/sigv2i) + sum(Gi^2*tBC1/sigv2i)*t(zvi)%*%(c2))
    Heta = Heta - 2*c1*sum(Gi^2*tBC1^2/sigv2i)*(1 + a2i^2 + d2*a2i) + c1^2*sum(Gi^2*tBC1/sigv2i)^2*(2 + a2i^2*(3 + d3) + 3*d2*a2i) + sum(c2*tBC1^2)*s*sqrt(c1)*(d2 + a2i) - 2*s*c1^1.5*sum(c2*tBC1)*sum(Gi^2*tBC1/sigv2i)*(d2 + a2i*(1 + d3)) + sum(c2*tBC1)^2*c1*d3
    
   } else if (model == "K1990modified"){
    Hbeta = Hbeta - s*t(xi)%*%sweep(tBC2, 1, sigv2i, "/")*sqrt(c1)*(a2i + d2) - s*(t(xi)%*%(Gi/sigv2i))%*%t(t(tBC2)%*%(c1*(s*(ei/sigv2i)*d3 - sqrt(c1)*a2i*(1 + d3)*(Gi/sigv2i) - sqrt(c1)*d2*(Gi/sigv2i))))
    Hdeleta = Hdeleta + t(zdeli[i,,drop=FALSE])%*%(t(t(tBC2)%*%(s*c1*(ei/sigv2i)*d3 - c1^1.5*(a2i*(1+d3) + d2)*(Gi/sigv2i)))/sigu2i[i])
    Hetagu = Hetagu + t(zui[i,,drop=FALSE])%*%(- t(Gi/sigv2i)%*%tBC2*c1^1.5/sigu2i[i]*(sqrt(c1) -a2i*mui[i]*(1+d3) + sqrt(c1)*a2i^2*0.5*(3+d3) - d2*(mui[i] - 1.5*a2i*sqrt(c1))) - t(ei/sigv2i)%*%tBC2*s*c1/sigu2i[i]*(mui[i]*d3 - d2*0.5*sqrt(c1) - 0.5*sqrt(c1)*a2i*(1+d3)))
    Hetagv = Hetagv + t(zvi)%*%sweep(tBC2, 1, Gi/sigv2i, "*")*c1*(1 + a2i^2 + a2i*d2) - (t(zvi)%*%(Gi^2/sigv2i))%*%(t(Gi/sigv2i)%*%tBC2)*c1^2*(1 + 0.5*a2i^2*(3 + d3) + 1.5*a2i*d2) - t(zvi)%*%sweep(tBC2, 1, ei/sigv2i, "*")*s*sqrt(c1)*(d2 + a2i) - (t(zvi)%*%(c2))%*%(t(ei/sigv2i)%*%tBC2)*c1*d3 + (a2i*(1+d3) + d2)*s*c1^1.5*(0.5*(t(zvi)%*%(Gi^2/sigv2i))%*%(t(ei/sigv2i)%*%tBC2) + (t(zvi)%*%(c2))%*%(t(Gi/sigv2i)%*%tBC2))
    Heta = Heta - c1*(1 + a2i^2 + a2i*d2)*t(tBC2)%*%sweep(tBC2, 1, sigv2i, "/") + c1^2*(2 + a2i^2*(3 + d3) + d2*a2i*3)*t(t(Gi/sigv2i)%*%tBC2)%*%(t(Gi/sigv2i)%*%tBC2) - c1^1.5*s*(d2+a2i*(1+d3))*t(t(ei/sigv2i)%*%tBC2)%*%(t(Gi/sigv2i)%*%tBC2) + c1*d3*t(t(ei/sigv2i)%*%tBC2)%*%(t(ei/sigv2i)%*%tBC2) - c1^1.5*s*(d2 + a2i*(1+d3))*t(t(Gi/sigv2i)%*%tBC2)%*%(t(ei/sigv2i)%*%tBC2)
    
   } else if (model == "K1990"){
    Hbeta = Hbeta +  s*t(xi)%*%sweep(tSK, 1, Gi^2*Gk/sigv2i, "*")*sqrt(c1)*(a2i + d2) - s*(t(xi)%*%(Gi/sigv2i))%*%t(t(tSK)%*%(c1*(-s*(ei*Gi^2*Gk/sigv2i)*d3 + sqrt(c1)*a2i*(1 + d3)*(Gi^3*Gk/sigv2i) + sqrt(c1)*d2*(Gi^3*Gk/sigv2i))))
    
    Hdeleta = Hdeleta + t(zdeli[i,,drop=FALSE])%*%(s*t(t(tSK)%*%(-s*c1*(c2*Gi*Gk)*d3 + c1^1.5*(a2i*(1+d3) + d2)*(Gi^3*Gk/sigv2i)))/sigu2i[i])
    
    Hetagu = Hetagu + t(zui[i,,drop=FALSE])%*%(t(Gi^3*Gk/sigv2i)%*%tSK*c1^1.5/sigu2i[i]*(sqrt(c1) -a2i*mui[i]*(1+d3) + sqrt(c1)*a2i^2*0.5*(3+d3) - d2*(mui[i] - 1.5*a2i*sqrt(c1))) + t(c2*Gi*Gk)%*%tSK*s*c1/sigu2i[i]*(mui[i]*d3 - d2*0.5*sqrt(c1) - 0.5*sqrt(c1)*a2i*(1+d3)))
    
    Hetagv = Hetagv - t(zvi)%*%sweep(tSK, 1, Gi^3*Gk/sigv2i, "*")*c1*(1 + a2i^2 + a2i*d2) + (t(zvi)%*%(Gi^2/sigv2i))%*%(t(Gi^3*Gk/sigv2i)%*%tSK)*c1^2*(1 + 0.5*a2i^2*(3 + d3) + 1.5*a2i*d2) + t(zvi)%*%sweep(tSK, 1, c2*Gi*Gk, "*")*s*sqrt(c1)*(d2 + a2i) + (t(zvi)%*%(c2))%*%(t(c2*Gi*Gk)%*%tSK)*c1*d3 - (a2i*(1+d3) + d2)*s*c1^1.5*(0.5*(t(zvi)%*%(Gi^2/sigv2i))%*%(t(c2*Gi*Gk)%*%tSK) + (t(zvi)%*%(c2))%*%(t(Gi^3*Gk/sigv2i)%*%tSK))
    
    Heta = Heta + c1*(1 + a2i^2 + a2i*d2)*t(tSK)%*%sweep(tSK, 1, Gi^3*Gk*(1- 3*Gi*Gk)/sigv2i, "*") + c1^2*(2 + a2i^2*(3 + d3) + d2*a2i*3)*t(t(Gi^3*Gk/sigv2i)%*%tSK)%*%(t(Gi^3*Gk/sigv2i)%*%tSK) - s*sqrt(c1)*(a2i + d2)*t(tSK)%*%sweep(tSK, 1, ei*Gi^2*Gk*(1 - 2*Gi*Gk)/sigv2i, "*") - c1^1.5*s*(d2 + a2i*(1 + d3))*(t(t(ei*Gi^2*Gk/sigv2i)%*%tSK)%*%(t(Gi^3*Gk/sigv2i)%*%tSK) + t(t(Gi^3*Gk/sigv2i)%*%tSK)%*%(t(ei*Gi^2*Gk/sigv2i)%*%tSK)) + c1*d3*t(t(ei*Gi^2*Gk/sigv2i)%*%tSK)%*%(t(ei*Gi^2*Gk/sigv2i)%*%tSK)
   } else {
    stop("Unknown model\n")
   }
  }
 }
 
 # if(eff.time.invariant){
 #  H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),		t(Hbdel)), 
 #             rbind(Hbgv, 	Hgv,		t(Hgvgu),		t(Hgvdel)),
 #             rbind(Hbgu, 	Hgvgu,		Hgu,			Hdelgu), 
 #             rbind(Hbdel,	Hgvdel,		t(Hdelgu),		Hdel))
 # } else {
 #  H <- cbind(rbind(Hb,	t(Hbgv),	t(Hbgu),		t(Hbdel),	t(Hbeta)), 
 #             rbind(Hbgv, 	Hgv,		t(Hgvgu),		t(Hgvdel),	t(Hetagv)),
 #             rbind(Hbgu, 	Hgvgu,		Hgu,			Hdelgu, 	t(Hetagu)), 
 #             rbind(Hbdel,	Hgvdel,		t(Hdelgu),		Hdel,		t(Hdeleta) ),
 #             rbind(Hbeta,	Hetagv,		Hetagu,			Hdeleta,	Heta ))
 # }

  H <- H[!coef.fixed,!coef.fixed]
 # print(H)
 
 # if(mean.u.0i.zero){
 #  grad <- grad[-(k+kv+ku+kdel)]
 #  H <- H[-(k+kv+ku+kdel),-(k+kv+ku+kdel)]
 # } 
 
 # print(length(grad))
 # print(dim(H))
 # print(length(theta))
 
 
 colnames(H) <- rownames(H) <- names(theta)
 return(H)
}

# .hess.panel <- function(theta, prod, my.n, k, kv, ku, kdel, yit, zvit, zui, xit, zdeli, model, timevar, maxT, ids, idvar, t0, eff.time.invariant = FALSE){
#  numDeriv::hessian(func = .ll.panel, x = theta, prod = prod, my.n = my.n, k = k, kv = kv, ku = ku, kdel = kdel, yit = yit, zvit = zvit, zui = zui, xit = xit, zdeli = zdeli, eff.time.invariant = eff.time.invariant, model = model, timevar = timevar, maxT = maxT, ids = ids, idvar = idvar, t0 = t0)
# }
# 
# .hess.panel <- function(theta, prod, my.n, k, kv, ku, kdel, yit, zvit, zui, xit, zdeli, model, timevar, maxT, ids, idvar, t0, eff.time.invariant = FALSE){
#  hessian1(func = .ll.panel, at = theta, prod = prod, my.n = my.n, k = k, kv = kv, ku = ku, kdel = kdel, yit = yit, zvit = zvit, zui = zui, xit = xit, zdeli = zdeli, eff.time.invariant = eff.time.invariant, model = model, timevar = timevar, maxT = maxT, ids = ids, idvar = idvar, t0 = t0)
# }

.is.negative.definite <- function (x, tol = 1e-16)
{
 # if (!is.square.matrix(x))
 # stop("argument x is not a square matrix")
 # if (!is.symmetric.matrix(x))
 # stop("argument x is not a symmetric matrix")
 # if (!is.numeric(x))
 # stop("argument x is not a numeric matrix")
 eigenvalues <- eigen(x, only.values = TRUE, symmetric = TRUE)$values
 # cat("Eigenvalues\n")
 # print(eigenvalues)
 n <- nrow(x)
 for (i in 1:n) {
  if (abs(eigenvalues[i]) < tol) {
   eigenvalues[i] <- 0
  }
 }
 if (any(eigenvalues >= 0)) {
  return(FALSE)
 }
 return(TRUE)
}

.make.neg.def <- function(hess){
 # K4 <- ncol(hess)
 # h0_ <- matrix(0, K4, K4)
 eigen1 <- eigen( hess )
 # eigen.val <- abs(eigen1$values)
 #   for( i in seq_len( K4 ) ){
 #     h0_ <- h0_ - abs(eigen1$values[i]) * eigen1$vectors[,i] %*% t(eigen1$vectors[,i])
 #     # h0_ <- h0_ - eigen.val[i] * eigen1$vectors[,i] %*% t(eigen1$vectors[,i])
 #   }
 return(-crossprod(t(eigen1$vectors)*abs(eigen1$values), t(eigen1$vectors)))
}

.make.invertible <- function(hess){
 K <- ncol(hess)
 ok <- FALSE
 adj <- sqrt(.Machine$double.eps)
 # tr0 <- sum( 1/eigen(H,only.values = T)$values )
 i <- 0
 while(!ok & i < 1000){
  i <- 1 + i
  hess <- hess + adj*diag(rep(1, K))
  mytry <- tryCatch( solve(-hess), error = function(e) e )
  ok <- !inherits( mytry, "error")
  # print(mytry)
  # ok <- all(diag(H) != 0)
  adj <- adj * 2
  # print(adj)
 }
 # print(c(i, adj))
 return(hess)
}

.make.zero.one <- function(qq){
 for (i in 1:length(qq)) {
  qqq <- qq[i]
  if(abs(qqq) > 1){
   while(abs(qqq) > 1){
    qqq <- qqq / 10
   }
   qq[i] <- qqq
  }
 }
 return(qq)
}

.mlmaximize.panel <- function(theta0, ll, gr = NULL, hess = NULL, gr.hess = NULL, alternate = NULL, BHHH = FALSE, level = 0.99, step.back = .Machine$double.eps^.5, reltol =  sqrt(.Machine$double.eps), lmtol =  sqrt(.Machine$double.eps), steptol =  .Machine$double.eps, digits = 4, when.backedup = sqrt(.Machine$double.eps), max.backedup = 17, print.level = 6, only.maximize = FALSE, maxit = 150, n = 100, ...){
 
 theta00 <- theta0
 
 # cat.print(lmtol)
 
 k4 <- length(theta0)
 
 # if(print.level >= 6){
 #  cat("\n=================")
 #  cat(" Initial values:\n\n", sep = "")
 #  print(theta0)
 # }
 # 
 # if(print.level >= 2){
 #  cat("\n=================")
 #  cat(" Maximization:\n\n", sep = "")
 # }
 
 # step.back = 2^-217
 
 # cat.print(n)
 
 ll0 <- ll(theta0, ...)
 ltol <- reltol * (abs(ll0) + reltol)
 # print(ltol)
 typf <- ll0
 theta1 <- theta0
 
 iter <- iter.total <- backedup <- backedups <- wasconcave <- wasconcaves <- 0
 
 if( is.na(ll0) | ll0 == -Inf ){
  if(print.level >= 2){
   cat("Could not compute ll at starting values: trying something else\n")
  }
  iter1 <- backedups
  repeat{
   iter1 <- iter1 + 1
   theta0 <- theta0 * runif(length(theta0), 0.98, 1.02) # not sure what to do here
   ll0 <- ll( theta0, ... )
   if(!is.na(ll0)) break
   if(iter1 == 55){
    stop("it's not gonna happen...")
   }
  }
  # backedups <- iter1
 }
 
 delta1 <- gHg <- s1 <- 1
 h1 <- tryCatch( 2, error = function(e) e )
 cant.invert.hess <- FALSE
 
 if(print.level >= 2){
  cat(paste("Iteration ",formatC(iter, width = 3)," (at starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
 }
 
 repeat{
  iter.total <- iter.total + 1
  # cat("backedup = ",backedup," backedups = ",backedups,"\n", sep = "")
  if(print.level >= 6){
   print(theta0)
  }
  
  # cumulate how many times did it backed-up in a row
  if(s1 < when.backedup){
   backedup <- backedup + 1
  } else {
   backedup <- 0
  }
  # print(s1)
  # cumulate how many times was concave
  if( inherits(h1, "error") ){
   wasconcave <- wasconcave + 1
  } else {
   wasconcave <- 0
  }
  
  # try different values if was concave more than @@@ times
  if(wasconcave == max.backedup){
   # start over
   if(print.level >= 2){
    cat("Not concave ",max.backedup," times in a row: trying something else (not concave ",wasconcaves+1," times in total)\n")
   }
   iter <- wasconcave <- backedup <- 0
   wasconcaves <- wasconcaves + 1
   theta0 <- theta0*runif(length(theta0), 0.9999, 1.0001) # not sure what to do here
   ll0 <- ll( theta0, ... )
   # if(print.theta) print(theta0)
   if(print.level >= 2){
    cat(paste("Iteration ",formatC(iter, width = 3),"  (at slightly perturbed starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
   }
  }
  
  # try different values if backed-up more than @@@ times
  if(backedup == max.backedup){
   # start over
   if(print.level >= 2){
    cat("Backed-up ",max.backedup," times in a row: trying something else (backup-up ",backedups+1," times in total)\n", sep = "")
   }
   iter <- backedup <- wasconcave <- 0
   backedups <- backedups + 1
   wasconcaves <- wasconcaves + 1
   theta0 <- theta0*runif(length(theta0), 0.9999, 1.0001) # not sure what to do here
   ll0 <- ll( theta0, ... )
   if(print.level >= 6){
    print(theta0)
   }
   if(print.level >= 2){
    cat(paste("Iteration ",formatC(iter, width = 3)," (at slightly perturbed starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
   }
  }
  
  # see if it calculated ll
  if( is.na(ll0) | ll0 == -Inf | ll0 == Inf | ll0 == 0 ){
   if(print.level >= 2){
    cat("Could not compute ll: trying something else\n")
   }
   iter1 <- backedups
   repeat{
    iter1 <- iter1 + 1
    # theta0 <- c( cons0, beta0, mu = 0, eta = 0, lnsv2 = -1*iter1/2, lnsu2 = -1*iter1/2)
    theta0 <- theta00*runif(length(theta0), 0.9999, 1.0001) # not sure what to do here
    ll0 <- ll( theta0, ... )
    if(!is.na(ll0) & ll0 != 0) break
    if(iter1 == 15){
     stop("it's not gonna happen... could not compute at initial and find feasible values")
    }
   }
   iter <- backedup <- 0
   backedups <- iter1
   if(print.level >= 2){
    cat(paste("Iteration ",formatC(iter, width = 3)," (at slightly perturbed starting values):       log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
   }
  }
  if(backedups == 5){
   stop("it's not gonna happen... backuped 5 times")
  }
  if(wasconcaves == 5){
   stop("it's not gonna happen... not concave 5 times")
  }
  iter <- iter + 1
  delta3 <- 1
  # step 2: direction vector
  
  # previous theta
  
  if(iter.total > 1) theta1 <- theta0 - s1 * d0
  
  # BHHH (faster, but different):
  # The Hessian is approximated as the negative
  # of the sum of the outer products of the gradients
  # of individual observations,
  # -t(gradient) %*% gradient = - crossprod( gradient )

  # gradient and hessian together
  g1_h0 <- gr.hess(theta0, ...)
  g1 <- g1_h0$grad
  h0 <- g1_h0$hessian1
  
  
  # gradient and hessian separately
  # g1 <- gr(theta0, ...)
  # # print(g1)
  # if(!is.null(alternate)) BHHH <- floor(iter/alternate) %% 2 == 1
  # h0 <- hess(theta0,  ...)
  
  # print(h0)
  
  
  # check if negative definite

  
  # eigen1 <- eigen( h0 )
  # # eigen.tol <- k4 * max(abs(eigen1$values)) * .Machine$double.eps # this is for positive definiteness
  # eigen.val <- ifelse(eigen1$values < .Machine$double.eps^.1, 0, eigen1$values)
  # eigen.val <- eigen1$values
  # # hess.pos.def <- sum(eigen1$values > eigen.tol) == k4
  # hess.neg.def <- !any(eigen.val >= 0)
  # 1. replace negative with small ones
  # eigen.val <- ifelse(eigen1$values < 0, .0001, eigen.val)
  # 2. replace negative with absolut values
  # eigen.val <- abs(eigen1$values)

  
  hess.neg.def <- .is.negative.definite(h0)
  h00 <- h0
  if(!hess.neg.def) h0 <- .make.neg.def(h0)
  # print(hess.neg.def)
  # make it negative definite if it is not already
  # if(!hess.neg.def){
  #  # cat(" !hess.neg.def; k4 =",k4,"length(eigen1$vectors[,i]) = ",length(eigen1$vectors[,1]),"\n")
  #  h0_ <- matrix(0, k4, k4)
  #  # eigen1 <- eigen( h0 )
  #  for( i in seq_len( k4 ) ){
  #   # h0_ <- h0_ - abs(eigen1$values[i]) * eigen1$vectors[,i] %*% t(eigen1$vectors[,i])
  #   h0_ <- h0_ - eigen.val[i] * crossprod(t(eigen1$vectors[,i]))
  #  }
  #  h0.old <- h0
  #  h0 <- h0_
  # }
  # print( is.negative.definite(h0) )
  # remember hessian and negative of its inverse from previous iter that could have been inverted
  # if( !cant.invert.hess ){
  #  h0_previous <- h0
  #  h1_previous <- h1
  # }
  # easier to invert positive definite matrix
  # h1 <- tryCatch( qr.solve(-h0, tol = 1e-16), error = function(e) e )
  h1 <- tryCatch( solve(-h0), error = function(e) e )
  if(inherits(h1, "error")) h0 <- .make.invertible(h0)
  h1 <- solve(-h0)
  # print(diag(h1))
  # check if it can be inverted
  cant.invert.hess <- FALSE
  cant.invert.hess <- inherits(h1, "error")
  # print(cant.invert.hess)
  if( cant.invert.hess ){
   # # print(h1)
   # if(print.level >= 2){
   #  cat(paste("cannot invert Hessian, using eigenvalues\n", sep = ""), sep = "")
   # }
   # # this was just to get the uninvertable hessian
   # # return(list(hess = h0, grad = g1))
   # # stop("this")
   # # @14@ this
   # eig1 <- eigen( -h0_previous )
   # d0 <- rep(0, length(theta0))
   # # eig2 <- ifelse(eig1$values < eps1, 1, eig1$values)
   # for (i in 1:length(eig1$values)){
   #  d0 <- d0 + (g1%*%eig1$vectors[,i])%*%eig1$vectors[,i] / eig1$values[i]
   # }
   # # @14@ could be done easier
   # # d0 <- qr.solve(-h0, g1, tol = 1e-134)
   # # print(g1)
   # # print(dim(g1))
   # # print(d0)
   # # print(dim(d0))
   # gHg <- sum( as.vector(g1) * d0)
   # # in the part of the ortogonal subspace where the eigenvalues
   # # are negative or small positive numbers, use steepest ascent
   # # in other subspace use NR step
   # # d0 <- ifelse(eigen(-h0, only.values = TRUE)$values < reltol, g1, d0)
   # gg <- sqrt( crossprod(g1) )
   # gHg <- gg
   # # d0 <- g1
   # # d0
  } else {
   # print(h1)
   # cat("dim H_inv",dim(h1)," length g1 ", length(g1), "\n")
   d0 <- as.vector( h1 %*% g1 )
   # d0 <- g1 # steepest
   gg <- sqrt( crossprod(g1) )
   # h1.old <- solve(-h0.old)
   gHg <- as.vector( t(g1) %*% h1 %*% g1 )
  }
  # gg_scaled <- gg * max( crossprod(theta0), crossprod(theta1) ) / max( abs(ll0), abs(typf))
  # theta_rel <- max( abs(theta0 - theta1) / apply( cbind( abs(theta0),abs(theta1) ), 1, max) )
  theta_rel <- max( abs(theta0 - theta1) / (abs(theta1)+1) )
  
  
  # begin stopping criteria calculated using new values of g1 and h1
  # if(s1 > when.backedup*10^-100 & delta1 != 17.17){ # if(s1 > when.backedup*10^-100 & !cant.invert.hess){
  if(s1 > when.backedup*1e-10 & delta1 != 17.17){ # if(s1 > when.backedup*10^-100 & !cant.invert.hess){
   # cat.print(abs(gHg))
   # cat.print(lmtol)
   if(abs(gHg) < lmtol & iter.total > 1){
    if(print.level >= 2){
     cat("\nConvergence given abs(gHg) = ",abs(gHg)," < lmtol\n", sep = "")
    }
    break
   }
   if(theta_rel < steptol*1e-100 & iter.total > 2){
    # print(theta_rel)
    if(print.level >= 2){
     cat("\nConvergence given relative change in theta = ",theta_rel," < steptol\n", sep = "")
    }
    break
   }
  }
  # end stopping criteria
  
  
  # use steepest ascent when backed-up
  if(s1 < when.backedup | delta1 == 17.17){ # was 1e-3
   # print("try this\n\n")
   ll0 <- ll( theta1, ... )
   eig1 <- eigen( -h0 )
   # cat("Gradient\n")
   # print(g1)
   # cat("Eigenvalues of the Hessian\n")
   # print(eig1$values)
   # cat("step before transformation 0\n")
   # print(d0)
   d0 <- as.vector( solve(-.make.neg.def(h0)) %*% g1 )
   # cat("step before transformation 1\n")
   # print(d0)
   # d0 <- ifelse(eig1$values < reltol, 0, d0)
   d0 <- ifelse(eig1$values < -1e-2, 0, d0)
   # d0 <- ifelse(abs(d0) > 1, .make.zero.one(d0), d0)
   # cat("step after transformation\n")
   # print(d0)
   # cat("\n\n")
   # theta0 <- theta0 - 1 * d0
  }
  
  
  
  # print(d0)
  # step 3: new guess
  # a: s = 1
  # b: funct(theta0 + d0) > funct(theta0)
  s1 <- 1
  # s1_ <- Inf
  # try.back.up.n <- 30
  # s1s <- 2^(3:-11)#seq(-5-1e-4, 5-1e-4, length.out = try.back.up.n)
  # my.loc.ll.max <- -1e100
  # for(qq in 1:15){
  #  theta1_ <- theta0 + s1s[qq] * d0
  #  ll1_ <- ll( theta1_, ... )
  #  if(is.finite(ll1_) & ll1_> my.loc.ll.max & ll1_ - my.loc.ll.max < 1e6){
  #   my.loc.ll.max <- ll1_
  #   s1_ <- s1s[qq]
  #  }
  # }
  # if(is.finite(s1_)){
  #  s1 <- s1_
  #  cat("opt s1 =", s1_,"opt ll ==",my.loc.ll.max, "ll before =" ,ll0,"\n")
  # } 
  # print(s1)
  # print(d0)
  # print(dim(d0))
  # print(theta0)
  # print(dim(theta0))
  theta1 <- theta0 + s1 * d0
  # print(12)
  # print(theta1)
  ll1 <- ll( theta1, ... )
  # print(13)
  delta2 <- ll1 - ll0
  flag <- (!is.na(delta2) & is.finite(delta2) & delta2 > 0 & delta2 < 1e+7)
  # begin Cases
  if( flag ){
   # begin Case 1: f(theta1) > f(theta0)
   ll.temp <- ll0
   # check if s1 = 2, 3, ... increases f even more
   while( flag ){
    if(print.level >= 6){
     cat(paste("\t\tCase 1: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
    }
    # cat("ll1 from before = ",ll1,"\n")
    ll0 <- ll1
    s1 <- s1 + 1
    theta1 <- theta0 + s1 * d0
    ll1 <- ll( theta1, ... )
    delta2 <- ll1 - ll0
    flag <- (!is.na(delta2) & is.finite(delta2) & delta2 > 0)
    if(is.infinite(ll1)){
     cat("s1 = ",s1,"\n")
     cat("ll1 = ",ll1,"\n")
     cat("ll0 = ",ll0,"\n")
     cat("ll.temp = ",ll.temp,"\n")
     cat("ll( theta0, ... ) = ",ll( theta0, ... ),"\n")
     cat("ll( theta1, ... ) = ",ll( theta1, ... ),"\n")
     # s1 <- s1 - 1
    }
   }
   # overall delta
   delta1 <- ll0 - ll.temp
   delta_rel <- abs(delta1 / ll.temp)
   # print(delta_rel)
   s1 <- s1 - 1
   # overwrite the values
   theta0 <- theta0 + s1 * d0
   # end Case 1: f(theta1) > f(theta0)
  } else {
   # begin Case 2: f(theta1) < f(theta0)
   # check only if s1=1/2 increases f
   s1 <- 0.5
   # s1_ <- Inf
   # try.back.up.n <- 30
   # s1s <- s1s <- 2^(3:-26)#seq(1e-10, 1-1e-4, length.out = try.back.up.n)
   # my.loc.ll.max <- -1e100
   # for(qq in 1:try.back.up.n){
   #  theta1_ <- theta0 + s1s[qq] * d0
   #  ll1_ <- ll( theta1_, ... )
   #  if(is.finite(ll1_) & ll1_> my.loc.ll.max & ll1_ - my.loc.ll.max  < 1e7){
   #   my.loc.ll.max <- ll1_
   #   s1_ <- s1s[qq]
   #  }
   # }
   # if(is.finite(s1_)){
   #  s1 <- s1_
   #  cat("opt s1 =", s1_,"opt ll ==",my.loc.ll.max, "ll before =" ,ll0,"\n")
   # }
   theta1 <- theta0 + s1 * d0
   ll1 <- ll( theta1, ... )
   # cat(" ll1 = ",ll1,", ll0 = ",ll0,", s = ",s1,"\n", sep = "")
   delta2 <- ll1 - ll0
   if(print.level >= 6){
    cat(paste("\t\tCase 2: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
   }
   flag2 <- (!is.na(delta2) & is.finite(delta2) & delta2 > 0)
   # end Case 2: f(theta1) < f(theta0)
   if( flag2 ){
    # begin Case 2a: f(theta1) > f(theta0)
    ll.temp <- ll0
    # check if s1=1/2^2,1/2^3,... increases f even more
    while( flag2 ){
     if(print.level >= 6){
      cat(paste("\t\t\tCase 2a: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
     }
     ll0 <- ll1
     s1 <- 0.5 * s1
     theta1 <- theta0 + s1 * d0
     ll1 <- ll( theta1, ... )
     # cat(" ll1 = ",ll1,", ll0 = ",ll0,", s = ",s1,"\n", sep = "")
     delta2 <- ll1 - ll0
     flag2 <- (!is.na(delta2) & is.finite(delta2) & delta2 > ltol)
    }
    # overall delta
    delta1 <- ll0 - ll.temp
    delta_rel <- abs(delta1 / ll.temp)
    # print(delta_rel)
    s1 <- 2 * s1
    # overwrite the values
    theta0 <- theta0 + s1 * d0
    # end Case 2a: f(theta1) > f(theta0)
   } else {
    # begin Case 2b: f(theta1) < f(theta0)
    ll.temp <- ll0
    # try s1=1/2^2,1/2^3,... so that f(theta1) > f(theta0)
    while ( !flag2 & s1 > step.back^2 ){
     s1 <- 0.5 * s1
     theta1 <- theta0 + s1 * d0
     ll1 <- ll( theta1, ... )
     delta2 <- ll1 - ll0
     if(print.level >= 6){
      cat(paste("\t\t\tCase 2b: s = ",s1,", delta = ",delta2,"\n", sep = ""), sep = "")
     }
     flag2 <- (!is.na(delta2) & is.finite(delta2) & delta2 > 0)
    }
    if( !flag2 | s1 < step.back^2 ){
     # stop("provide different starting values")
     delta1 <- 17.17
    } else {
     # overwrite the values
     delta1 <- delta2
     delta_rel <- abs(delta1 / ll.temp)
     ll0 <- ll1
     theta0 <- theta0 + s1 * d0
    }
    # end Case 2b: f(theta1) < f(theta0)
   }
  }
  
  
  if(print.level >= 2){
   if( cant.invert.hess ){
    cat(paste("Iteration ",formatC(iter, width = 3)," (hessian is ",ifelse(BHHH, "BHHH", "analytical"),", ",formatC(iter.total, width = 3)," in total):   log likelihood = ",format(ll0, digits = 13)," (not concave)\n", sep = ""), sep = "")
   } else if (s1 <= when.backedup) {
    cat(paste("Iteration ",formatC(iter, width = 3)," (hessian is ",ifelse(BHHH, "BHHH", "analytical"),", ",formatC(iter.total, width = 3)," in total):   log likelihood = ",format(ll0, digits = 13)," (backed up)\n", sep = ""), sep = "")
   } else {
    cat(paste("Iteration ",formatC(iter, width = 3)," (hessian is ",ifelse(BHHH, "BHHH", "analytical"),", ",formatC(iter.total, width = 3)," in total):   log likelihood = ",format(ll0, digits = 13),"\n", sep = ""), sep = "")
   }
  }
  
  # printing criteria
  if(print.level >= 5){
   if( cant.invert.hess ){
    cat(paste(" (in iter ",formatC(iter, width = 3),": delta = ",format(delta1, digits = 6),"; s = ",format(s1, digits = 6),"; quasi-gHg = ",format(gHg, digits = 6),"; theta_rel_change = ",format(theta_rel, digits = 6),")\n\n", sep = ""), sep = "")
    if(print.level >= 6.5){
     print(theta0)
    }
   } else {
    cat(paste(" (in iter ",formatC(iter, width = 3),": delta = ",format(delta1, digits = 6),"; s = ",format(s1, digits = 6),"; gHg = ",format(abs(gHg), digits = 6),"; theta_rel_change = ",format(theta_rel, digits = 6),")\n\n", sep = ""), sep = "")
    if(print.level >= 5.5){
     print(theta0)
    }
   }
   cat("Par\n")
   print(theta0)
   cat("Gradient\n")
   print(g1)
   cat("\n")
  }
  # print(s1)
  if(s1 > when.backedup & !cant.invert.hess){ # if(s1 > when.backedup^2 & !cant.invert.hess){
   # ltol <- reltol * (abs(ll0) + reltol)
   # print(cant.invert.hess)
   if(delta1 > 0 & !is.na(delta_rel) & delta_rel < ltol*1e-16 & iter.total > 1){
    if(print.level >= 2){
     cat("\nConvergence given delta = ",delta_rel," < dtol\n", sep = "")
    }
    g1_h0 <- gr.hess(theta0, ...)
    g1 <- g1_h0$grad
    h0 <- g1_h0$hessian
    
    # g1 <- gr(theta0, ...)
    # h0 <- hess(theta0,  ...)
    # h1 <- tryCatch( qr.solve(-h0), error = function(e) e )
    break
   }
  }
  if(iter.total >= maxit){
   cat("\n Maximum number of iterations (",iter.total,") reached without convergence\n", sep = "")
   cat(" Only 'theta' from iteration ",iter," will be returned\n\n", sep = "")
   print(theta0)
   return(list(par = theta0, converged = 0))
   break
   cat(" Only 'theta' from iteration ",iter," will be returned\n", sep = "")
  }
 } # end repeat
 
 if( !only.maximize & !cant.invert.hess){
  names(ll0) <- NULL
  colnames(h1) <- rownames(h1) <- names(g1) <- names(theta0)
  
  # sqrt(crossprod(g1))
  
  b0 <- theta0
  suppressWarnings( sd0 <- sqrt( diag( h1 ) ) )
  sd.nas <- is.na(sd0)
  # cat.print(sum(sd.nas))
  if(sum(sd.nas) > 0){
   sd1 <- sqrt(diag( solve(-.make.neg.def(h00), tol = 1e-268 ) ))
   sd0[sd.nas] <- sd1[sd.nas]
  }
  # cat.print(diag( solve(-.make.neg.def(h00) ) ))
  # cat("1\n")
  # cat.print(diag( solve(-.make.invertible(h00) ) ))
  # cat.print(diag( solve(-.make.neg.def(h00), tol =1e-68 ) ))
  t0 <- b0 / sd0
  p0 <- pt(abs(t0), n-length(b0), lower.tail = FALSE) * 2
  t10 <- qt((1-0.01*level)/2, n-length(b0), lower.tail = FALSE)
  t17 <- cbind( b0, sd0, t0, p0, b0 - t10*sd0, b0 + t10*sd0)
  # t17 <- cbind( b0, sd0, t0, p0)
  colnames(t17) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)", paste("",level,"_CI_LB", sep = ""), paste("",level,"_CI_UB", sep = ""))
  # colnames(t17) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
  # t17
  
  if(print.level >= 2){
   cat(paste("\nFinal log likelihood = ",format(ll0, digits = 13),"\n\n", sep = ""), sep = "")
  }
  # cat(paste("Stoc. frontier normal/",distribution,"\n", sep = ""), sep = "")
  if(print.level >= 7){
   cat("\nCoefficients:\n\n", sep = "")
   printCoefmat(t17[,1:4], digits = digits)
  }
  
  return(list(par = theta0, table = t17, gradient = g1, vcov = h1, ll = ll0, gg = gg, gHg = gHg, theta_rel_ch = theta_rel, converged = 1))
 } else {
  return(list(par = theta0, gradient = g1, vcov = h1, ll = ll0, gg = gg, gHg = gHg, theta_rel_ch = theta_rel, converged = 1))
 }
 
}

# Print the estimation results
.printpanel2nd = function(x, digits, kb, kvi, ku0, kdeli, model, eff.time.invariant, mean.u.0i.zero, na.print = "NA", max.name.length, mycutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), mysymbols = c("***", "**", "*", ".", " ")){
 
 Cf = cbind(ifelse(x[,1, drop = F]> 999, formatC(x[,1, drop = F], digits = 1, format = "e",width = 10), formatC(x[,1, drop = F], digits = digits, format = "f", width = 10)),
            ifelse(x[,2, drop = F]>999, formatC(x[,2, drop = F], digits = 1, format = "e", width = 10), formatC(x[,2, drop = F], digits = digits, format = "f", width = 10)),
            ifelse(x[,3, drop = F]>999, formatC(x[,3, drop = F], digits = 1, format = "e", width = 7), formatC(x[,3, drop = F], digits = 2, format = "f", width = 7)),
            ifelse(x[,4, drop = F]>999, formatC(x[,4, drop = F], digits = 1, format = "e", width = 10), formatC(x[,4, drop = FALSE], digits = digits, format = "f", width = 10)),
            formatC(mysymbols[findInterval(x = x[,4], vec = mycutpoints)], flag = "-"))
 
 row.names(Cf) <- formatC(row.names(Cf), width = max(nchar(row.names(Cf))), flag = "-")
 cat("",rep(" ", max.name.length+6),"Coef.        SE       z       P>|z|\n", sep = "")
 dimnames(Cf)[[2]] <- rep("", dim(Cf)[[2]])
 cat("",rep("_", max.name.length+42-1),"", "\n", "Frontier", "\n", sep = "")
 print.default(Cf[1:kb,,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 
 
 cat("",rep("-", max.name.length+42-1),"", "\n", "Random noise component: log(sigma_vit^2)", "\n", sep = "")
 print.default(Cf[(kb+1):(kb+kvi),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 
 cat("",rep("-", max.name.length+42-1),"", "\n", "Inefficiency component: log(sigma_u0i^2)", "\n", sep = "")
 
 cat("(id-specific time-invariant)", "\n", sep = "")
 print.default(Cf[(kb+kvi+1):(kb+kvi+ku0),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 
 if(!mean.u.0i.zero){
  cat("",rep("-", max.name.length+42-1),"", "\n", "Inefficiency component: conditional mean of the ", "\n", sep = "")
  cat("truncated-normal distribution", "\n", sep = "")
  cat("(id-specific time-invariant)", "\n", sep = "")
  print.default(Cf[(kb+kvi+ku0+1):(kb+kvi+ku0+kdeli),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 } else {
  kdeli <- 0
 }
 

 if(!eff.time.invariant){
  if(model == "BC1992"){
   Keta <- 1
   components <- "Component"
  } else {
   Keta <- 2
   components <- "Components"
  }
  cat(rep("-", max.name.length+42-1),"\n Function of inefficiency time change \n", sep = "")
  # cat("(function of inefficiency time evolution)", "\n", sep = "")
  print.default(Cf[(kb+kvi+ku0+kdeli+1):(kb+kvi+ku0+kdeli+Keta),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 } else {
  Keta <- 0
 }
 
 if(nrow(x) > kb+kvi+ku0+kdeli+Keta){
  cat("",rep("-", max.name.length+42-1),"", "\n", "Auxiliary parameters", "\n", sep = "")
  print.default(Cf[(kb+kvi+ku0+kdeli+Keta+1):nrow(x),,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
  
 }
 
 cat("",rep("_", max.name.length+42-1),"", "\n", sep = "")
 invisible(x)
}

# Technical efficiencies and prediction intervals
.u2efftnm.panel <- function( eit, sigv2it, sigu2i, mui, Git, alpha = alpha, prod = prod, cost.eff.less.one, ids, idvar, timevar, t0, it.names) {
 # if(prod){sn = -1} else {sn = 1}
 sn <- sn1 <- ifelse(prod, -1, 1)
 if(!prod & cost.eff.less.one) sn1 <- -1
 # sn <- -1
 
 u_mean <- te_jlms_mean <- te_jlms_mode <- te_bc <- te_l <- te_u <- c()
 
 # define indices
 m.end <- cumsum(t0)
 m.begin <- c(1,(m.end+1)[-length(t0)])
 
 for(i in 1:length(ids)){
  
  
 # for(i in ids){
  # sample.i <- idvar == i
  sample.i <- (m.begin[i]):(m.end[i])
  ei = eit[sample.i]; 
  sigv2i = sigv2it[sample.i]
  Gi = Git[sample.i];
  
  s2i <- (1/sigu2i[i] + sum(Gi^2/sigv2i))^(-1)
  mi <- (mui[i]/sigu2i[i] + sn*sum(ei*Gi/sigv2i))*s2i
  zi <- mi / sqrt(s2i)
  
  point.est.mean <- mi + sqrt(s2i) * dnorm(zi) / pnorm(zi)
  point.est.mode <- ifelse( mi >= 0, mi, 0 )
  
  # if(i < 2){
  #  cat("\n current i = ",i ,"\n")
  #  cat.print(sample.i)
  #  cat.print(ei)
  #  cat.print(mui[i])
  #  cat.print(sigv2i)
  #  cat.print(Gi)
  #  cat.print(s2i)
  #  cat.print(mi)
  #  cat.print(zi)
  #  cat("\n")
  # }
  
  u_mean <- append(u_mean, Gi*point.est.mean)
  te_jlms_mean <- append(te_jlms_mean, exp(sn1*Gi*point.est.mean))
  te_jlms_mode <- append(te_jlms_mode, exp(sn1*Gi*point.est.mode))
  te_bc <- append(te_bc, exp(sn1*mi*Gi + 0.5*s2i*Gi^2) * pnorm(zi + sn1 * sqrt(s2i)*Gi)/pnorm(zi))
  zl    <- qnorm( 1 - alpha / 2 * pnorm(zi) )
  zu    <- qnorm( 1 - ( 1 - alpha/2 ) * pnorm(zi) )
  te_l  <- append(te_l, exp(Gi*(sn1*mi - zl*sqrt(s2i))))
  te_u  <- append(te_u, exp(Gi*(sn1*mi - zu*sqrt(s2i))))
 }
 tymch <- data.frame(idvar, timevar, te_l, te_jlms_mean, te_jlms_mode, te_bc, te_u, u_mean)
 # tymch <- data.frame(idvar, timevar, te_l, te_jlms_mean, te_jlms_mode, te_bc, te_u, u_mean)
 colnames(tymch) <- c(it.names, "Lower bound","JLMS", "Mode", "BC","Upper bound", "u")
 row.names(tymch) <- NULL
 return(tymch)
}

# truncreg ---------------------------------------------------------------------

# a and b are lower and upper truncations
# .h.trunc <- function(al, be, a1 = 1, b1 = 1, a = -Inf, b = Inf, print = FALSE)
 # (ifelse( b == Inf, 0, b1*dnorm(be) ) - ifelse( a == -Inf, 0, a1*dnorm(al) )  ) /
 # ( pnorm(be) - pnorm(al) )

# log likelihood
# .ll.trunc <- function(y, x, beta, sigma, a, b) - 0.5*length(y)*log(2*pi) - length(y)*log(sigma) - 0.5*sigma^(-2)*sum((y - x %*% beta)^2) - sum( log( pnorm((b - x %*% beta)/sigma, log.p = FALSE) - pnorm((a - x %*% beta)/sigma, log.p = FALSE) ) )

# function with formula or matrices ---------------------------------------
.prepareYXllul <- function(formula, ll = -Inf, ul = Inf, data, subset, sysnframe = 1, ...) {
 # needed.frame <- sys.nframe() - 1
 needed.frame <- sys.nframe() - sysnframe
 # cat("sys.nframe() = ",sys.nframe(),"\n")
 # cat("needed.frame = ",needed.frame,"\n")
 mf0 <- match.call(expand.dots = FALSE, call = sys.call(sys.parent(n = needed.frame)))

 ll.clas <- class(ll)
 if (ll.clas == "formula"){
  if (length(all.vars(ll)) == 1){
   ll.form <- TRUE
  } else {
   stop("invalid formula for lower limit for left-truncation; 'll' must be '~ variableName'")
  }
 } else if (ll.clas == "numeric") {
  if (length(ll) == 1){
   ll.form <- FALSE
  } else {
   stop("invalid lower limit for left-truncation; 'll' must be a scalar")
  }
 } else {
  stop("invalid lower limit for left-truncation; 'll' must be a formula or a scalar")
 }

 ul.clas <- class(ul)
 if (ul.clas == "formula"){
  if (length(all.vars(ul)) == 1){
   ul.form <- TRUE
  } else {
   stop("invalid formula for upper limit for right-truncation; 'ul' must be '~ variableName'")
  }
 } else if (ul.clas == "numeric") {
  if (length(ul) == 1){
   ul.form <- FALSE
  } else {
   stop("invalid upper limit for right-truncation; 'ul' must be a scalar")
  }
 } else {
  stop("invalid upper limit for right-truncation; 'ul' must be a formula or a scalar")
 }

 form1 <- Formula::Formula(formula)

 # cat(" print(form1)\n", sep = "")
 # print(form1)

 if (ll.form) {
  form1 <- Formula(as.formula(paste("",deparse(form1, width.cutoff = 500L)," | ",ll[2]," ", sep = "")))
 }
 if (ul.form) {
  form1 <- Formula(as.formula(paste("",deparse(form1, width.cutoff = 500L)," | ",ul[2]," ", sep = "")))
 }

 # cat(" print(form1)\n", sep = "")
 # print(form1)

 # check if it is a matrix
 datasupplied <- !(match("data", names(mf0), 0) == 0)
 subssupplied <- !(match("subset", names(mf0), 0) == 0)

 if(subssupplied & !datasupplied){
  stop("Cannot specify 'subset' without specifying 'data'\n")
 }

 if(datasupplied){
  # begin get a logical vector equal TRUE if !missing
  # first using data and subset to get x without NA
  mf <- mf0
  mf$formula <- form1 #formula( form )
  m <- match(c("formula", "data", "subset"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  X <- as.matrix(model.matrix(mt, mf))
  # now get the names in the entire data
  # esample <- seq_len( nrow(data) ) %in% as.numeric(rownames(X))
  esample <- rownames(data) %in% rownames(X)
  if(sum(esample)==0){
   stop("no valid data points")
  }
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing

  # get the data
  # print(form1)
  dataesample <- model.frame(form1, data = data[esample,])
  # cat(" dim(dataesample)\n", sep = "")
  # print(dim(dataesample))
  # print(head(dataesample))

  # this is my full LHS
  Y <- model.part(form1, data = dataesample, lhs = 1, drop = TRUE)
  # cat(" print(Y)\n", sep = "")
  # print(Y)
  n.full <- length(Y)
  # Y <- as.matrix( model.matrix(formula(form1, lhs = 1, rhs = 0), data = data[esample,]))
  # print(Y)
  X <- as.matrix( model.matrix(formula(form1, lhs = 0, rhs = 1), data = dataesample))
  # print(X)

  if (ll.form) {
   LL <- model.matrix(formula(form1, lhs = 0, rhs = 2), data = dataesample)[,2]
  } else {
   LL <- rep(ll, n.full)
  }

  if (ul.form) {
   if (ll.form) {
    UL <- model.matrix(formula(form1, lhs = 0, rhs = 3), data = dataesample, drop = TRUE)[,2]
   } else {
    UL <- model.matrix(formula(form1, lhs = 0, rhs = 2), data = dataesample, drop = TRUE)[,2]
   }
  } else {
   UL <- rep(ul, n.full)
  }

  # cat(" LL\n", sep = "")
  # print(LL)
  # cat(" UL\n", sep = "")
  # print(UL)


  # flag those that are required for regression
  flag <- (Y > LL & Y < UL)
  # cat(" dim(X)\n", sep = "")
  # print(dim(X))
  # cat(" flag\n", sep = "")
  # print(flag)


  # get subsets
  # cat(" Y\n", sep = "")
  Y  <- Y[flag]
  n <- length(Y)
  # cat(" X\n", sep = "")
  X  <- X[flag, , drop = FALSE]
  # print(dim(X))
  # cat(" LL\n", sep = "")
  LL <- LL[flag]
  # print(LL)
  # cat(" UL\n", sep = "")
  UL <- UL[flag]
  # print(UL)
 }
 # if data are not supplied
 else {
  # begin get a logical vector equal TRUE if !missing

  # first using data and subset to get XZ without NA
  mf <- mf0
  mf$formula <- formula( form1 )
  subsetsupplied <- !(match("subset", names(mf), 0) == 0)
  if(subsetsupplied) stop("Subset with matrices is not allowed")
  m <- match(c("formula"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  # mf <- eval(mf, parent.frame())
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  Y <- as.matrix(model.response(mf))
  n <- nrow(Y)
  XZ <- as.matrix(model.matrix(mt, mf))
  # print(head(XZ))
  # get a logical vector equal TRUE if !missing
  with.na <- model.frame(mt, na.action = na.pass)
  esample <- rownames(with.na) %in% rownames(XZ)
  if(sum(esample)==0){
   stop("no valid data points")
  }
  # print(table(esample))
  # end get a logical vector equal TRUE if !missing

  # get the data

  # get number of vars in X
  mf <- mf0
  m <- match(c("formula"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf, sys.frame(sys.parent(n = needed.frame)))
  mt <- attr(mf, "terms")
  X <- as.matrix(model.matrix(mt, mf))
  # print(head(X))
  k <- ncol(X)

  Y <- Y[esample,]
  n.full <- length(Y)
  X <- XZ[esample,]

  if (ll.form) {
   LL <- XZ[esample,(k+1)]
  } else {
   LL <- rep(ll, n.full)
  }

  if (ul.form) {
   if (ll.form) {
    UL <- XZ[esample,(k+2)]
   } else {
    UL <- XZ[esample,(k+1)]
   }
  } else {
   UL <- rep(ul, n.full)
  }

  # cat(" LL\n", sep = "")
  # print(LL)
  # cat(" UL\n", sep = "")
  # print(UL)


  # flag those that are required for regression
  flag <- (Y > LL & Y < UL)
  # cat(" dim(X)\n", sep = "")
  # print(dim(X))
  # cat(" flag\n", sep = "")
  # print(flag)

  # get subsets
  # cat(" Y\n", sep = "")
  Y  <- Y[flag]
  n <- length(Y)
  # cat(" X\n", sep = "")
  X  <- X[flag, , drop = FALSE]
  # print(dim(X))
  # cat(" LL\n", sep = "")
  LL <- LL[flag]
  # print(LL)
  # cat(" UL\n", sep = "")
  UL <- UL[flag]
  # print(UL)

 }

 tymch <- list(Y = Y, X = X, LL = LL, UL = UL, n = n, n.full = n.full, ll.form = ll.form, ul.form = ul.form, esample = esample, nontruncsample = flag)
 class(tymch) <- "npsf"
 return(tymch)
}


# 4comp -------------------------------------------------------------------

.fourC_gtrelnls <- function(theta, prod = TRUE,
                       V0,U0,yit,xit,
                       ids,idvar,nt,my.n,R,idlenmax,
                       Ktheta,n_threads=1){
 do.prod <- ifelse( prod, 1, -1 )
 k <- Ktheta - 4
 log.lmd <- log.sig <- log.sigv0 <- log.sigu0 <- TRUE
 if(log.lmd)   theta[k+1] <- exp( theta[k+1] )
 if(log.sig)   theta[k+2] <- exp( theta[k+2] )
 if(log.sigv0) theta[k+3] <- exp( theta[k+3] )
 if(log.sigu0) theta[k+4] <- exp( theta[k+4] )
 t51 <- .C("gtre_ll", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(my.n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(Ktheta), lnls = as.double(1) )
 # t51 <- tryCatch( .C("gtre_ll", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(Ktheta), lnls = as.double(1) ), error = function(e) e )
 if(inherits(t51, "error")){
  return(NA)
 } else {
  return(t51$lnls)
 }
}


.fourC_gtregrad <- function(theta, prod = TRUE,
                       V0,U0,yit,xit,
                       ids,idvar,nt,my.n,R,idlenmax,
                       Ktheta,n_threads=1){
 do.prod <- ifelse( prod, 1, -1 )
 k <- Ktheta - 4
 log.lmd <- log.sig <- log.sigv0 <- log.sigu0 <- TRUE
 if(log.lmd)   theta[k+1] <- exp( theta[k+1] )
 if(log.sig)   theta[k+2] <- exp( theta[k+2] )
 if(log.sigv0) theta[k+3] <- exp( theta[k+3] )
 if(log.sigu0) theta[k+4] <- exp( theta[k+4] )
 # tt <- length(theta)
 t51 <- tryCatch( .C("gtre_grad", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(my.n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(Ktheta), lnls = as.double(1), grad = as.double(rep(17, Ktheta)) ), error = function(e) e )
 # t51 <- .C("gtre_grad", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(my.n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(Ktheta), lnls = as.double(1), grad = as.double(rep(17, Ktheta)) )
 if(inherits(t51, "error")){
  return(rep(NA,Ktheta))
 } else {
  mygrad <- t51$grad
  if(log.lmd) mygrad[k+1] <- mygrad[k+1] * theta[k+1]
  if(log.sig) mygrad[k+2] <- mygrad[k+2] * theta[k+2]
  if(log.sigv0) mygrad[k+3] <- mygrad[k+3] * theta[k+3]
  if(log.sigu0) mygrad[k+4] <- mygrad[k+4] * theta[k+4]
  # names(mygrad) <- mynames
  return(mygrad)
 }
}

.lower2full <- function(lower.tri.vec, names = NULL){
 kk <- sqrt(0.5^2 + 2 * length(lower.tri.vec) ) - 0.5
 full.m <- matrix(0, nrow = kk, ncol = kk)
 full.m[lower.tri(full.m, diag = TRUE)] <- lower.tri.vec
 full.m <- full.m + t(full.m)
 for(i in seq_len(kk)) full.m[i,i] <- full.m[i,i]/2
 if(!is.null(names)) colnames(full.m) <- rownames(full.m) <- names
 full.m
}

.fourC_gtrehess <- function(theta, BHHH = FALSE, prod = TRUE,
                       V0,U0,yit,xit,
                       ids,idvar,nt,my.n,R,idlenmax,
                       Ktheta,n_threads=1, 
                       grad.log = rep(0, length(theta)), ... ){
 # if(log.lmd)   theta[k+1] <- exp( theta[k+1] )
 # if(log.sig)   theta[k+2] <- exp( theta[k+2] )
 # if(log.sigv0) theta[k+3] <- exp( theta[k+3] )
 # if(log.sigu0) theta[k+4] <- exp( theta[k+4] )
 do.prod <- ifelse( prod, 1, -1 )	
 do.BHHH <- ifelse( BHHH, 1, 0)
 k <- Ktheta - 4
 log.lmd <- log.sig <- log.sigv0 <- log.sigu0 <- TRUE
 if(log.lmd)   theta[k+1] <- exp( theta[k+1] )
 if(log.sig)   theta[k+2] <- exp( theta[k+2] )
 if(log.sigv0) theta[k+3] <- exp( theta[k+3] )
 if(log.sigu0) theta[k+4] <- exp( theta[k+4] )
 tt <- length(theta)
 len.tri.tt <- ( tt^2 + tt ) / 2
 t51 <- tryCatch( .C("gtre", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(my.n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(k+4), BHHH = as.integer(do.BHHH),lnls = as.double(1), grad = as.double(rep(17, tt)), hesstri = as.double(rep(17, len.tri.tt)) ), error = function(e) e )
 # t51 <- .C("gtre", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(my.n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(k+4), BHHH = as.integer(do.BHHH), lnls = as.double(1), grad = as.double(rep(17, tt)), hesstri = as.double(rep(17, len.tri.tt)) )
 if(inherits(t51, "error")){
  return( matrix(NA,tt,tt) )
 } else {
  # grad
  mygrad <- t51$grad
  if(log.lmd) mygrad[k+1] <- mygrad[k+1] * theta[k+1]
  if(log.sig) mygrad[k+2] <- mygrad[k+2] * theta[k+2]
  if(log.sigv0) mygrad[k+3] <- mygrad[k+3] * theta[k+3]
  if(log.sigu0) mygrad[k+4] <- mygrad[k+4] * theta[k+4]
  # cat.print(mygrad)
  # hess
  myhess <- .lower2full( t51$hesstri, names = NULL )
  grad.log <- mygrad
  # cat.print(grad.log)
  if(log.lmd){
   myhess[k+1,] <- myhess[,k+1] <- myhess[k+1,] * theta[k+1]
   myhess[k+1,k+1] <- 2*myhess[k+1,k+1] * theta[k+1]^2 + grad.log[k+1]
  } 
  if(log.sig){
   myhess[k+2,] <- myhess[,k+2] <- myhess[k+2,] * theta[k+2]
   myhess[k+2,k+2] <- 2*myhess[k+2,k+2] * theta[k+2]^2 + grad.log[k+2]
  }
  if(log.sigv0){
   myhess[k+3,] <- myhess[,k+3] <- myhess[k+3,] * theta[k+3]
   myhess[k+3,k+3] <- 2*myhess[k+3,k+3] * theta[k+3]^2 + grad.log[k+3]
  }
  if(log.sigu0){
   myhess[k+4,] <- myhess[,k+4] <- myhess[k+4,] * theta[k+4]
   myhess[k+4,k+4] <- 2*myhess[k+4,k+4] * theta[k+4]^2 + grad.log[k+4]
  }
  return(myhess)
 }
}

.fourC_gtregr_gtrehess <- function(theta, BHHH = FALSE, prod = TRUE,
                            V0,U0,yit,xit,
                            ids,idvar,nt,my.n,R,idlenmax,
                            Ktheta,n_threads=1, 
                            grad.log = rep(0, length(theta)), ... ){
 # if(log.lmd)   theta[k+1] <- exp( theta[k+1] )
 # if(log.sig)   theta[k+2] <- exp( theta[k+2] )
 # if(log.sigv0) theta[k+3] <- exp( theta[k+3] )
 # if(log.sigu0) theta[k+4] <- exp( theta[k+4] )
 do.prod <- ifelse( prod, 1, -1 )	
 do.BHHH <- ifelse( BHHH, 1, 0)
 k <- Ktheta - 4
 log.lmd <- log.sig <- log.sigv0 <- log.sigu0 <- TRUE
 if(log.lmd)   theta[k+1] <- exp( theta[k+1] )
 if(log.sig)   theta[k+2] <- exp( theta[k+2] )
 if(log.sigv0) theta[k+3] <- exp( theta[k+3] )
 if(log.sigu0) theta[k+4] <- exp( theta[k+4] )
 tt <- length(theta)
 len.tri.tt <- ( tt^2 + tt ) / 2
 t51 <- tryCatch( .C("gtre", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(my.n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(k+4), BHHH = as.integer(do.BHHH), lnls = as.double(1), grad = as.double(rep(17, tt)), hesstri = as.double(rep(17, len.tri.tt)) ), error = function(e) e )
 # t51 <- .C("gtre", prod = as.integer(do.prod), V = as.double(V0), U = as.double(U0), N = as.integer(my.n), R = as.integer(R), Z = as.double(xit), NT = as.integer(nt), C = as.double(yit), ids = as.double(ids), idvar = as.double(idvar), idlenmax = as.integer(idlenmax), theta = as.double(theta), K4  = as.integer(k+4), BHHH = as.integer(do.BHHH), lnls = as.double(1), grad = as.double(rep(17, tt)), hesstri = as.double(rep(17, len.tri.tt)) )
 if(inherits(t51, "error")){
  return( matrix(NA,tt,tt) )
 } else {
  # grad
  mygrad <- t51$grad
  if(log.lmd) mygrad[k+1] <- mygrad[k+1] * theta[k+1]
  if(log.sig) mygrad[k+2] <- mygrad[k+2] * theta[k+2]
  if(log.sigv0) mygrad[k+3] <- mygrad[k+3] * theta[k+3]
  if(log.sigu0) mygrad[k+4] <- mygrad[k+4] * theta[k+4]
  # cat.print(mygrad)
  # hess
  myhess <- .lower2full( t51$hesstri, names = NULL )
  grad.log <- mygrad
  # cat.print(grad.log)
  if(log.lmd){
   myhess[k+1,] <- myhess[,k+1] <- myhess[k+1,] * theta[k+1]
   myhess[k+1,k+1] <- 2*myhess[k+1,k+1] * theta[k+1]^2 + grad.log[k+1]
  } 
  if(log.sig){
   myhess[k+2,] <- myhess[,k+2] <- myhess[k+2,] * theta[k+2]
   myhess[k+2,k+2] <- 2*myhess[k+2,k+2] * theta[k+2]^2 + grad.log[k+2]
  }
  if(log.sigv0){
   myhess[k+3,] <- myhess[,k+3] <- myhess[k+3,] * theta[k+3]
   myhess[k+3,k+3] <- 2*myhess[k+3,k+3] * theta[k+3]^2 + grad.log[k+3]
  }
  if(log.sigu0){
   myhess[k+4,] <- myhess[,k+4] <- myhess[k+4,] * theta[k+4]
   myhess[k+4,k+4] <- 2*myhess[k+4,k+4] * theta[k+4]^2 + grad.log[k+4]
  }
  # return(myhess)
  return(list(grad = mygrad, hessian1 = myhess))
 }
}

.ghk <- function(C, chol.vcov = TRUE,
                 from = rep(-Inf, ncol(C)), to = rep(Inf, ncol(C)),
                 draws = NULL,
                 simtype = c("halton", "random"), R = 10000, base.halton = 2){
 # if(R %% 2 == 1) R <- R+1
 # if(use.halton){
 # halton(17)
 # }
 if(!chol.vcov){
  C <- t(chol(C))
 }
 K <- ncol(C)
 z <- matrix(0, nrow = R, ncol = K-1)
 a1 <- rep(from[1], R)
 b1 <- rep(to[1], R)
 pi <- rep(1, R)
 # print(c(length(draws), R*(K-1), K))
 # cat.print(dim(draws))
 if(is.numeric(draws) & length(draws) >= R*(K-1)){
  myrand <- draws
  # cat("draws works\n")
 } else {
  if(simtype == "halton"){
   # myrand <- sHalton(floor(R*(K-1)), base = base.halton)
   myrand <- npsf::halton(R, n.bases = K-1, scale.coverage = TRUE, shuffle = TRUE)
   # cat("draws does not works\n")
  } else {
   # myrand <- runif(floor(R*(K-1)))
   myrand <- matrix( runif(floor(R*(K-1))), ncol = K-1, nrow = R)
  }
 }
 # cat("length in ghk = ",length(myrand),"\n")
 # cat.print(dim(myrand))
 for(k in seq_len(K)){
  if(k > 1){
   mysum <- z[, seq_len(k-1), drop = FALSE] %*% t(C[k, seq_len(k-1), drop = FALSE])
   a1 <- rep(from[k], R) - mysum
   b1 <- rep(to[k], R) - mysum
  }
  a1 <- pnorm( a1/C[k,k] )
  b1 <- pnorm( b1/C[k,k] )
  pi <- pi * (b1 - a1)
  if( k < K ){
   # myrand1 <- myrand[((k-1)*R+1) : (k*R)] #sHalton(floor(R/2), base = base.halton)
   myrand1 <- myrand[,k]
   # myrand <- c(myrand, 1 - myrand)
   z[, k] <- qnorm( a1 +  myrand1 * (b1 - a1) )
  }
 }
 mean(pi, na.rm = TRUE)
}

# svi2 : Ti x 1
# sv02 :  1 x 1
# sui2 : Ti x 1
# su02 :  1 x 1

.e_exp_tu <- function(ei, id, sui2, svi2, sv02, su02, prod = TRUE,
                      simtype = c("halton", "random"), seed = 17345168,
                      R = 500, halton.base = NULL, random.primes = FALSE,
                      inv.tol = .Machine$double.eps, print.level = 4){
 Ti <- length (ei)
 do.prod <- ifelse( prod, 1, -1)
 A <- -1 * do.prod * cbind (1, diag(Ti) )
 if(Ti == 1){
  Sig <- matrix(svi2 + sv02, 1, 1)
 }
 else {
  Sig <- diag(svi2) +  matrix(sv02, nrow = Ti, ncol = Ti)
 }
 # cat.print(c( su02, sui2 ))
 su02.zero <- FALSE
 if(su02 == 0 | su02 < sqrt(.Machine$double.eps)){
  su02 <- .Machine$double.eps
  su02.zero <- TRUE
 }
 sui2.zero <- sui2 == 0
 sui2 <- ifelse(sui2 == 0, .Machine$double.eps, sui2)
 V_1 <- diag ( 1/c( su02, sui2 ) )
 Sig_1 <- solve(Sig)
 Lmd <- tryCatch( solve( V_1 + t(A) %*% Sig_1 %*% A, tol = inv.tol ), error = function(e) e )
 # cat.print(Lmd)
 some_error <- inherits(Lmd, "error")
 if(!some_error){
  R_ <- Lmd %*% t(A) %*% Sig_1
  # cat("\n There was no error in the first type of LMD calculations\n")
  # # begin check
  # V <- diag ( c( su02, sui2 ) )
  # SAVA_1 <- tryCatch( solve(Sig + A %*% V %*% t(A), tol = inv.tol ), error = function(e) e )
  # if(inherits(Lmd, "error")){
  #   stop("Cannot invert some matrix, so no efficiencies for this ID")
  # }
  # R. <- V %*% t(A) %*% SAVA_1
  # Lmd. <- V - R. %*% A %*% V
  # cat("\nLmd's are equal",all.equal(Lmd,Lmd.),"; R's are equal",all.equal(R_,R.),"\n")
  # # end check
 } else {
  # cat("\n There was an error in the first type of LMD calculations\n")
  V <- diag ( c( su02, sui2 ) )
  SAVA_1 <- tryCatch( solve(Sig + A %*% V %*% t(A), tol = inv.tol ), error = function(e) e )
  if(inherits(Lmd, "error")){
   stop("Cannot invert some matrix, so no efficiencies for ID ",id," ", call. = FALSE)
  }
  R_ <- V %*% t(A) %*% SAVA_1
  Lmd <- V - R_ %*% A %*% V
  # Re <- R  %*% ei
 }
 Re <- R_ %*% ei
 # print(Lmd)
 # print(summary(Re))
 CholLmd <- t(chol(Lmd))
 # print(summary(CholLmd))
 if(Ti == 1){
  mydraws <- matrix(runif(R), ncol = 1)
 } else {
   if(simtype == "halton"){
     set.seed(seed)
     if(is.null(halton.base)){
       # deviates for each ID come from own base, but then scaled
       mydraws <- npsf::halton(R, n.bases = Ti, random.primes = random.primes, seed = seed, start = 1, scale.coverage = TRUE, shuffle = TRUE)
     } else {
       if(halton.base > 100008){
         stop("halton.base should be smaller than 100,008")
       } else {
         # deviates for all ID come from a sequence with said base
         mydraws <- matrix(  npsf::halton(R*Ti, bases = halton.base, random.primes = random.primes, seed = seed, start = 1, scale.coverage = TRUE, shuffle = TRUE), nrow = R, ncol = Ti)
         # cat("some halton.base",halton.base,"  ")
       }
     }
   } else {
     mydraws <- matrix( runif(floor(R*Ti)), nrow = R, ncol = Ti)
   }
 }
 # cat.print(head(mydraws))
 ghkReLmd <- .ghk(CholLmd, to = Re, simtype = simtype, draws = mydraws, R = R)
 te_it <- NULL
 for(tt in 1:(Ti+1) ){
  t00 <- -Re[tt] + 0.5 * Lmd[tt,tt]
  t11 <- exp( t00 )
  t22 <- .ghk( CholLmd, to = Re - Lmd[ , tt], simtype = simtype, draws = mydraws, R = R )
  temp <- t11 * t22 / ghkReLmd
  te_it <- c(te_it, temp )
 }
 # check if approximation was OK, correct if it wasn't
 te_it_finite <- is.finite(te_it)
 # if(any(!te_it_finite)) cat.print(te_it)
 te_it [!te_it_finite] <- NA
 range_star <- te_max <- NA
 if(sum(te_it_finite) > 0){
  te_max <- max(te_it, na.rm = TRUE)
  if(sum(te_it_finite) > 1) range_star <- te_max - sort(te_it [!te_it_finite])[2]
 }
 if(is.finite(te_max)){
  if(!is.na(te_max) && te_max > 1){
   # cat.print(id)
   # cat.print(id1)
   # cat.print(te_it)
   if(print.level > 5) {
    # cat.print(summary(te_it))
    .su(list(te_it, sui2, svi2, sv02, su02, ei), print = TRUE, width = 5, format = "fg", drop0trailing = FALSE, names = c("te_it", "te_it_alt", "sui2", "svi2", "sv02", "su02", "rez"))
    # cat.print(te_it)
    cat("Range is ",max(te_it, na.rm = TRUE) - min(te_it, na.rm = TRUE),": (",max(te_it, na.rm = TRUE)," ",min(te_it, na.rm = TRUE),"), ghkReLmd = ",ghkReLmd,"\n", sep = "")
    cat("Oops...\n")
    # print(te_it)
   }
  }
 }
 te_pers <- te_it[1]
 te_resi <- te_it[-1]
 if(su02.zero) te_pers <- 1
 te_resi <- ifelse(sui2.zero, 1, te_resi)
 tymch <- list(te_pers = te_pers, te_resi = te_resi, ghkReLmd = ghkReLmd, range_star = range_star)
 return(tymch)
}

# Print the estimation results
.printgtresfhet = function(x, digits, kb, kv0, ku0, kvi, kui, na.print = "NA", max.name.length, mycutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), mysymbols = c("***", "**", "*", ".", " ")){
 Cf = cbind(ifelse(x[,1, drop = F]> 999, formatC(x[,1, drop = F], digits = 1, format = "e",width = 10), formatC(x[,1, drop = F], digits = digits, format = "f", width = 10)),
            ifelse(x[,2, drop = F]>999, formatC(x[,2, drop = F], digits = 1, format = "e", width = 10), formatC(x[,2, drop = F], digits = digits, format = "f", width = 10)),
            ifelse(x[,3, drop = F]>999, formatC(x[,3, drop = F], digits = 1, format = "e", width = 7), formatC(x[,3, drop = F], digits = 2, format = "f", width = 7)),
            ifelse(x[,4, drop = F]>999, formatC(x[,4, drop = F], digits = 1, format = "e", width = 10), formatC(x[,4, drop = FALSE], digits = digits, format = "f", width = 10)),
            formatC(mysymbols[findInterval(x = x[,4], vec = mycutpoints)], flag = "-"))
 
 row.names(Cf) <- formatC(row.names(Cf), width = max(nchar(row.names(Cf))), flag = "-")
 cat("",rep(" ", max.name.length+6),"Coef.        SE       z       P>|z|\n", sep = "")
 dimnames(Cf)[[2]] <- rep("", dim(Cf)[[2]])
 cat("",rep("_", max.name.length+42-1),"", "\n", "Frontier", "\n", sep = "")
 print.default(Cf[1:kb,,drop=FALSE], quote = FALSE, right = TRUE, na.print = na.print)
 
 # cat("____________________________________________________", "\n", "Frontier", "\n")
 # cat("____________________________________________________", "\n", "Frontier")
 # print.default(Cf[1:kb,,drop=F], quote = FALSE, right = TRUE, na.print = na.print)
 if(kv0 > 0){
  # cat("----------------------------------------------------", "\n", "Random effects component: log(sigma_v0^2)", "\n")
  cat("",rep("-", max.name.length+42-1),"", "\n", "log(lambda)", "\n", sep = "")
  # dimnames(Cf)[[2]] <- rep("", dim(Cf)[[2]])
  print.default(Cf[(kb+1):(kb+kv0),,drop=F], quote = FALSE, right = TRUE, na.print = na.print)
 }
 if(ku0 > 0){
  # cat("----------------------------------------------------", "\n", "Persistent inefficiency component: log(sigma_u0^2)", "\n")
  cat("",rep("-", max.name.length+42-1),"", "\n", "log(sigma)", "\n", sep = "")
  print.default(Cf[(kb+kv0+1):(kb+kv0+ku0),,drop=F], quote = FALSE, right = TRUE, na.print = na.print)
 }
 # cat("----------------------------------------------------", "\n", "Random noise component: log(sigma_v^2)", "\n")
 cat("",rep("-", max.name.length+42-1),"", "\n", "Random noise component: log(sigma_v^2)", "\n", sep = "")
 print.default(Cf[(kb+kv0+ku0+1):(kb+kv0+ku0+kvi),,drop=F], quote = FALSE, right = TRUE, na.print = na.print)
 if(kui > 0){
  # cat("----------------------------------------------------", "\n", "Transient inefficiency component: log(sigma_u^2)", "\n")
  cat("",rep("-", max.name.length+42-1),"", "\n", "Persistent inefficiency component: log(sigma_u0^2)", "\n", sep = "")
  print.default(Cf[(kb+kv0+ku0+kvi+1):(kb+kv0+ku0+kvi+kui),,drop=F], quote = FALSE, right = TRUE, na.print = na.print)
 }
 cat("",rep("_", max.name.length+42-1),"", "\n", sep = "")
 invisible(x)
}

rescale <- function(x, lb = min(x), ub = max(x)){
  lb + ((x-min(x))/(max(x)-min(x)))*(ub-lb)
}