R/Utility.r

Defines functions .sdea_cook_mpi_internal .sdea_mpi_internal .sdea_cook_internal .sdea_internal .Mal_Ben_Output .Mal_Ben_Input .checkoption .checkWeightRestriction .checkDataGood .checkData

#**************************************************************************
#
# Copyright Aurobindh Kalathil Puthanpura, ETA, PSU, 2015
# Use granted under BSD license terms
#
# R DEA Multiplier Model Package
#
# Author: Aurobindh Kalathil Puthanpura
# Revision: 1
# Date: 2015-20-9
#
# U t i l i t y f i l e
# internal use only
#
#***************************************************************************

#
# Package Global variables
#
#multiplierdea.version = "0.1.0"
#multiplierdea.id = paste0("multiplierdea version ", multiplierdea.version)

# Load ipsolverAPI library

# .loadLibrary <- function() {
#   if (requireNamespace("lpSolveAPI",quietly = TRUE))
#   {
#     library("lpSolveAPI")
#   } else
#   {
#     install.packages("lpSolveAPI",repos = "http://cran.us.r-project.org")
#     library("lpSolveAPI")
#   }
#
#   if (requireNamespace("lpSolveAPI" ,quietly = TRUE))
#   {
#     return(TRUE)
#   } else
#   {
#     return(FALSE)
#   }
#
# }

#
# Options Strings
#

options.orientation.l <- c ("input", "output")
options.rts.l <- c ("vrs" , "crs")
options.phase.l <- c("mal", "ben")

#
# Solver Status Codes
#
dict.solveStatus<-cbind(key= c(0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13),
            desc = c("optimal solution found","the model is sub-optimal", "the model is infeasible",
                     "the model is unbounded", "the model is degenerate", "numerical failure encountered",
                     "process aborted", "timeout", "the model was solved by presolve","the branch and bound routine failed",
                     "the branch and bound was stopped because of a break-at-first or break-at-value",
                     "a feasible branch and bound solution was found","no feasible branch and bound solution was found"))

# Function: .checkDate
# Check for valid types
# Check add names to unnamed variables

.checkData <- function(x, name) {
  # Is value a dataframe, array or vector?
  if (!is.data.frame(x) && !is.array (x) && (length(x) < 2))
    stop(name,
         " is not a matrix, array (1 or 2 dimensions) or data.frame",
         call. =
           FALSE)
  if (length(dim(x)) > 2)
    stop(name, " is greater then 2 dimensions", call. = FALSE)

  # If data.frame - convert to array
  if (is.data.frame(x))
    x <- data.matrix(x)

  # If vector convert to a 1 x N array
  if (is.vector (x))
    X <- array(x, c (length(x), 1))

  # Check that a l l numeric
  if (!is.numeric(x))
    stop(name, " must be numeric", call. = FALSE)

  # Add DMU names i f empty
  if (length(dim(x)) >= 1 && is.null(rownames(x)))
    rownames(x) <- paste0("DMU", 1:nrow(x))

  # Add Col names i f empty
  if (length(dim(x)) >= 2 && is.null(colnames(x)))

    colnames(x) <- paste0(name, 1:ncol(x))
  return(x)
}


# Function .CheckDataGood
# Check input & output data and warn about problems for dea
.checkDataGood <- function(x, y) {
  status <- TRUE

  # Check for any zero's
  # check for na' s *
  # Check for no positive values on input & output
  if (any(x == 0, na.rm = TRUE)) {
    cat(
      "Warning, data has DMU's with inputs that are zero, th is may cause numerical
      problems\n"
    )
    status <- FALSE
  }

  if (any(y == 0, na.rm = TRUE)) {
    cat(
      "Warning, data has DMU's with outputs that are zero, this may cause
      numerical problems\n"
    )
    status <- FALSE
  }

  for (k in (1:nrow (x))) {
    if (all(x[k, ] <= 0, na.rm = TRUE) & all(y[k, ] <= 0, na.rm = TRUE)) {
      cat(
        "warning, dmu # " ,
        k ,
        "has no positive inputs or outputs, this may cause
        numerical problems\n"
      )
      status <- FALSE
    }
  }
  return(status)
  }

# Function .checkWeightRestriction
# Check if the weight restriction dataframe has values
.checkWeightRestriction <- function(x, y, weightRestriction) {
  #Get the coloum names from the IP and OP data set
  colValues <- c(colnames(x), colnames(y))

  #Check lower bound values for NA
  if (any(is.na(weightRestriction[, 1]), na.rm = FALSE))
  {
    stop("Lower bound for the weight restriction table has NA values",
         call. = FALSE)
  }

  #Check upper bound values for NA
  # if(any(is.na(weightRestriction[,4]),na.rm = FALSE))
  # {
  #   stop("Upper bound for the weight restriction table has NA values", call. = FALSE)
  # }


  #Check numerator column names
  for (i in 1:nrow(weightRestriction)) {
    if (length(which(colValues == weightRestriction[i, 2])) == 0)
    {
      stop(
        "Numerator values does not match with the column names from input and output values",
        call. = FALSE
      )
    }
  }

  #Check denominator column names
  for (i in 1:nrow(weightRestriction)) {
    if (length(which(colValues == weightRestriction[i, 3])) == 0)
    {
      stop(
        "Denominator values does not match with the column names from input and output values",
        call. = FALSE
      )
    }
  }

}

# Function .checkoption()
# Check that value in legal l i s t options
# If on legal list, return lowercase short form option
#


.checkoption <- function(value, name, options.l = c(TRUE, FALSE)) {

  if (length(value) != 1) {
    stop("illegal vector for ",
         name,
         " option must be single value",
         call. = FALSE)
  }
  # If options are character
  if (is.character(options.l[1])) {
    if (is.character(value)) {
      tmp.value <- tolower(value)
      # if (tmp.value %in% options.1)
      i <- charmatch(tmp.value, options.l, nomatch = -1)
      if (i > 0)
        return(options.l[i])
    }

  } else if (is.logical(options.l[1])) {

    # Logical options
    if (is.logical(value)) {
      return(value)
    }
  } else {
    # Numeric options
    if (is.numeric(value) && is.finite(value) && (value >= 0)) {
      return(value)
    }
    options.l <- c("numeric >= 0")
  }
  stop(
    "Illegal value=",
    value,
    " for ",
    name,
    "; legal values are: ",
    paste(options.l, collapse = ", ") ,
    call. = FALSE
  )
}

.Mult_lP <- function (x_t, y_t, rts = 'crs', weightRestriction) {
  Result.rts <- 'crs'
  if (rts == 'vrs') {
    Result.rts <- 'vrs'
  } else
  {
    Result.rts <- 'crs'
  }

  Result.orientation <- "Input"
  Result.Inputvalues <- x_t
  Result.Outputvalues <- y_t
  Result.eff <- matrix(nrow = nrow(x_t))
  Result.Lambda <- matrix(nrow = nrow(x_t), ncol = nrow(x_t))
  Result.vx <- matrix(nrow = nrow(x_t), ncol = ncol(x_t))
  Result.uy <- matrix(nrow = nrow(y_t), ncol = ncol(y_t))
  Result.Free_weights <- matrix(nrow = nrow(x_t), ncol = 2)
  Result.StatusData <- vector(mode="logical", nrow(x_t))



  #free sign variable
  fs <- 0

  switch(rts,
         'crs' = fs <- 0,
         'vrs' = fs <- 2,
         fs <- 2)

  dmu <- row.names(x_t)

  # Assign names to Matrix for efficiency
  rownames(Result.eff) <- dmu
  colnames(Result.eff)[1] <- c("Eff")

  # Assign names to Matrix for efficiency
  rownames(Result.Lambda) <- dmu
  colnames(Result.Lambda) <- dmu

  # Assign names to Matrix for IP Weights
  rownames(Result.vx) <- dmu
  colnames(Result.vx) <- colnames(x_t)

  # Assign names to Matrix for OP Weights
  rownames(Result.uy) <- dmu
  colnames(Result.uy) <- colnames(y_t)

  # Assign names to Matrix for Free Weights
  rownames(Result.Free_weights) <- dmu
  colnames(Result.Free_weights) <- c("F1", "F2")



  for (d in 1:length(dmu)) {
    ip_v <- ncol(x_t)
    op_v <- ncol(y_t)

    # Set objective function for DMU under consideration
    c <- 0
    for (i in 1:ncol(y_t)) {
      c[i] <- c(y_t[d, i])
    }


    # if vrs add the free variables else omit
    if (rts == 'vrs') {
      # Adding the 1 and -1 for free variable
      if (length(c) != 0) {
        c[length(c) + 1] <- 1
        c[length(c) + 1] <- -1
      }
    }

    # indices for objective function
    ind_obj <- 0
    for (i in 1:ncol(y_t)) {
      ind_obj[i] <- c(ip_v + i)
    }

    # if vrs add the free variables indices else omit
    if (rts == 'vrs') {
      # indices for objective function for free variable
      if (length(ind_obj) != 0) {
        ind_obj[length(ind_obj) + 1] <- ind_obj[length(ind_obj)] + 1
        ind_obj[length(ind_obj) + 1] <- ind_obj[length(ind_obj)] + 1
      }
    }

    lps.model <- make.lp(0, ip_v + op_v + fs)
    set.objfn(lps.model, c, indices = ind_obj)

    # Input constraint for the DMU under consideration
    ip_con_dmu <- 0
    for (i in 1:ncol(x_t)) {
      ip_con_dmu[i] <- x_t[d, i]
    }

    #indices for DMU input constraints
    ind_ip_con <- 0
    for (i in 1:ncol(x_t)) {
      ind_ip_con[i] <- i
    }

    add.constraint(lps.model, ip_con_dmu, "=", rhs = 1, ind_ip_con)

    # Multiply input values by -1 for building the constraints for all OMU's
    x_temp <- -1 * x_t

    # Construct constraints for all dmus
    for (i in 1:nrow(x_t)) {
      dmu_con <- 0
      for (j in 1:ncol(x_temp)) {
        dmu_con[j] <- x_temp[i, j]
      }
      for (k in 1:ncol(y_t)) {
        dmu_con[ncol(x_temp) + k] <- y_t[i, k]
      }

      # If rts is vrs then add the free variable else omit
      if (rts == 'vrs') {
        #Adding free variable in the constraint
        if (length(dmu_con) != 0) {
          dmu_con[length(dmu_con) + 1] <- 1
          dmu_con[length(dmu_con) + 1] <- -1
        }

      }
      add.constraint(lps.model, dmu_con, "<=", rhs = 0)
    }
    #Weight Restrictions
    if (!missing(weightRestriction)) {
      colvalues <- c(colnames(x_t), colnames(y_t))

      for (i in 1:nrow(weightRestriction))
      {
        #First constraint
        add.constraint(lps.model,
                       c(weightRestriction[i, 1], -1),
                       ">=",
                       0,
                       indices = c(
                         which(colvalues == weightRestriction[i, 2]),
                         which(colvalues == weightRestriction[i, 3])
                       ))
        #Second constraint
        # Check if upper bound exists (Length ==4)
        if (length(weightRestriction) == 4) {
          #Check if upper bound is not Inf, NA, or NaN

          if (!is.nan(weightRestriction[i, 4]) &&
              !is.na(weightRestriction[i, 4]) &&
              !is.infinite(weightRestriction[i, 4])) {
            add.constraint(
              lps.model,
              c(1, -1 * weightRestriction[i, 4]),
              "<=",
              0,
              indices = c(
                which(colvalues == weightRestriction[i, 2]),
                which(colvalues == weightRestriction[i, 3])
              )
            )
          }
        }
      }
    }

    lp.control(lps.model, sense = 'max')

    #print(lps.model)
    Result.StatusData[d]<-solve(lps.model)

    Result.eff[d] <- get.objective(lps.model)
    temp_lamdba <- get.dual.solution(lps.model)
    Result.Lambda[d, ] <- temp_lamdba[3:(length(dmu) + 2)]


    temp_weight <- get.variables(lps.model)

    #Store the IP & OP weights in respective Matrices
    Result.vx[d, 1:ncol(x_t)] <- temp_weight[1:ncol(x_t)]
    Result.uy[d, 1:ncol(y_t)] <-
      temp_weight[(ncol(x_t) + 1):(ncol(x_t) + ncol(y_t))]
    if (rts == 'vrs') {
      Result.Free_weights[d, 1:2] <-
        temp_weight[(ncol(x_t) + ncol(y_t) + 1):length(temp_weight)]
    }
  }

  # Calculating the transpose of Lambda values for HCU calculation
  lambda_Transpose <- t(Result.Lambda)

  # Calculating HCU data for Input
  Hcu_input <- matrix(nrow = ncol(lambda_Transpose), ncol = ncol(x_t))
  colnames(Hcu_input) <- colnames(x_t)
  rownames(Hcu_input) <- rownames(lambda_Transpose)

  for (i in 1:ncol(lambda_Transpose)) {
    e <- x_t * lambda_Transpose[, i]
    for (j in 1:ncol(e)) {
      Hcu_input[i, j] <- sum(e[, j])
    }
  }

  # Calculating HCU data for Output
  Hcu_output <- matrix(nrow = ncol(lambda_Transpose), ncol = ncol(y_t))
  colnames(Hcu_output) <- colnames(y_t)
  rownames(Hcu_output) <- rownames(lambda_Transpose)

  for (i in 1:ncol(lambda_Transpose)) {
    e <- y_t * lambda_Transpose[, i]
    for (j in 1:ncol(e)) {
      Hcu_output[i, j] <- sum(e[, j])
    }
  }


  #Get the Status code and description from the LP model
  Result.StatusDesc<-sapply(Result.StatusData,function(x) dict.solveStatus[dict.solveStatus[,1] == x ,2])
  Result.status<-data.frame(code = Result.StatusData, desc = Result.StatusDesc)
  rownames(Result.status) <- dmu




  return(
    list(
      "rts" = Result.rts,
      "Orientation" = Result.orientation,
      "InputValues" =
        Result.Inputvalues,
      "OutputValues" = Result.Outputvalues,
      "Efficiency" = Result.eff,
      "Lambda" = Result.Lambda,
      "HCU_Input" = Hcu_input,
      "HCU_Output" = Hcu_output,
      "vx" = Result.vx,
      "uy" = Result.uy,
      "Free_Weights" = Result.Free_weights,
      "Model_Status" = Result.status
    )
  )
}

.Mult_OP <- function (x_t, y_t, rts = 'crs', weightRestriction) {
  Result.rts <- 'crs'

  if (rts == 'vrs') {
    Result.rts <- 'vrs'
  } else
  {
    Result.rts <- 'crs'
  }
  Result.orientation <- "Output"
  Result.Inputvalues <- x_t
  Result.Outputvalues <- y_t
  Result.eff_rec <- matrix(nrow = nrow(x_t))
  Result.Lambda <- matrix(nrow = nrow(x_t), ncol = nrow(x_t))
  Result.vx <- matrix(nrow = nrow(x_t), ncol = ncol(x_t))
  Result.uy <- matrix(nrow = nrow(y_t), ncol = ncol(y_t))
  Result.Free_weights <- matrix(nrow = nrow(x_t), ncol = 2)
  Result.StatusData <- vector(mode="logical", nrow(x_t))

  #free sign variable
  fs <- 0

  switch(rts,
         'crs' = fs <- 0,
         'vrs' = fs <- 2,
         fs <- 2)

  dmu <- row.names(x_t)

  # Assign names to Matrix for efficiency
  rownames(Result.eff_rec) <- dmu
  colnames(Result.eff_rec)[1] <- c("Eff")

  # Assign names to Matrix for efficiency
  rownames(Result.Lambda) <- dmu
  colnames(Result.Lambda) <- dmu

  # Assign names to Matrix for IP Weights
  rownames(Result.vx) <- dmu
  colnames(Result.vx) <- colnames(x_t)

  # Assign names to Matrix for OP Weights
  rownames(Result.uy) <- dmu
  colnames(Result.uy) <- colnames(y_t)

  # Assign names to Matrix for Free Weights
  rownames(Result.Free_weights) <- dmu
  colnames(Result.Free_weights) <- c("F1", "F2")


  for (d in 1:length(dmu)) {
    ip_v <- ncol(x_t)
    op_v <- ncol(y_t)

    # set objective function for DMU under consideration
    c <- 0
    for (i in 1:ncol(x_t)) {
      c[i] <- c(x_t[d, i])
    }

    #  for(i in 1:ncol(y_t)){
    # c[length(c) + 1] <- 0
    # }

    # If rts is vrs then add the free variable else omit
    if (rts == 'vrs') {
      # Adding the 1 and -1 for free variable
      if (length(c) != 0) {
        c[length(c) + 1] <- 1
        c[length(c) + 1] <- -1
      }
    }

    # indices for objective function
    ind_obj <- 0
    for (i in 1:ncol(x_t)) {
      ind_obj[i] <- i
    }

    # for(i in 1:ncol(y_t)){
    # ind_obj[ncol(x_t) + i ] <-nco1(x_t) + i
    # }

    # If rts is vrs then add the free variable indices else omit
    if (rts == 'vrs') {
      # indices for objective function fo r free variable
      if (length(ind_obj) != 0) {
        ind_obj[length(ind_obj) + 1] <- ip_v + op_v + 1
        ind_obj[length(ind_obj) + 1] <- ip_v + op_v + 2
      }
    }

    lps.model <- make.lp(0, ip_v + op_v + fs)
    set.objfn(lps.model, c, indices = ind_obj)


    # Output constraint for the DMU under consideration
    op_con_dmu <- 0
    for (i in 1:ncol(y_t)) {
      op_con_dmu[i] <- y_t[d, i]
    }

    #indices for DMU input constraints
    ind_op_con <- 0
    for (i in 1:ncol(y_t)) {
      ind_op_con[i] <- ip_v + i
    }
    add.constraint(lps.model, op_con_dmu, "=", rhs = 1, ind_op_con)

    # Multiply input values by -1 for building the constraints for all DMU's
    y_temp <- -1 * y_t

    # Construct constraints for all dmus
    for (i in 1:nrow(x_t)) {
      dmu_con <- 0
      for (j in 1:ncol(x_t)) {
        dmu_con[j] <- x_t[i, j]
      }
      for (k in 1:ncol(y_temp)) {
        dmu_con[ncol(x_t) + k] <- y_temp[i, k]
      }

      # if the rts is vrs then add free variable else omit
      if (rts == 'vrs') {
        if (length(dmu_con) != 0) {
          dmu_con[length(dmu_con) + 1] <- 1
          dmu_con[length(dmu_con) + 1] <- -1
        }
      }
      add.constraint(lps.model, dmu_con, ">=", rhs = 0)
    }

    #Weight Restrictions
    if (!missing(weightRestriction)) {
      colvalues <- c(colnames(x_t), colnames(y_t))

      for (i in 1:nrow(weightRestriction))
      {
        #First constraint
        add.constraint(lps.model,
                       c(weightRestriction[i, 1], -1),
                       ">=",
                       0,
                       indices = c(
                         which(colvalues == weightRestriction[i, 2]),
                         which(colvalues == weightRestriction[i, 3])
                       ))

        #Second constraint
        # Check if upper bound exists (Length ==4)
        if (length(weightRestriction) == 4) {
          #Check if upper bound is not Inf, NA, or NaN

          if (!is.nan(weightRestriction[i, 4]) &&
              !is.na(weightRestriction[i, 4]) &&
              !is.infinite(weightRestriction[i, 4])) {
            add.constraint(
              lps.model,
              c(1, -1 * weightRestriction[i, 4]),
              "<=",
              0,
              indices = c(
                which(colvalues == weightRestriction[i, 2]),
                which(colvalues == weightRestriction[i, 3])
              )
            )
          }
        }

      }
    }

    lp.control(lps.model, sense = 'min')


    Result.StatusData[d] <- solve(lps.model)

    Result.eff_rec[d] <- get.objective(lps.model)
    temp_lamdba <- get.dual.solution(lps.model)
    Result.Lambda[d, ] <- temp_lamdba[3:(length(dmu) + 2)]

    temp_weight <- get.variables(lps.model)
    #Store the IP & OP weights in respective Matrices
    Result.vx[d, 1:ncol(x_t)] <- temp_weight[1:ncol(x_t)]
    Result.uy[d, 1:ncol(y_t)] <-
      temp_weight[(ncol(x_t) + 1):(ncol(x_t) + ncol(y_t))]
    if (rts == 'vrs') {
      Result.Free_weights[d, 1:2] <-
        temp_weight[(ncol(x_t) + ncol(y_t) + 1):length(temp_weight)]
    }
  }
  Result.eff <- 1 / Result.eff_rec


  # Calculating the transpose of Lambda values for HCU calculation
  lambda_Transpose <- t(Result.Lambda)

  # Calculating HCU data for Input
  Hcu_input <- matrix(nrow = ncol(lambda_Transpose), ncol = ncol(x_t))
  colnames(Hcu_input) <- colnames(x_t)
  rownames(Hcu_input) <- rownames(lambda_Transpose)

  for (i in 1:ncol(lambda_Transpose)) {
    e <- x_t * lambda_Transpose[, i]
    for (j in 1:ncol(e)) {
      Hcu_input[i, j] <- sum(e[, j])
    }
  }

  # Calculating HCU data for Output
  Hcu_output <- matrix(nrow = ncol(lambda_Transpose), ncol = ncol(y_t))
  colnames(Hcu_output) <- colnames(y_t)
  rownames(Hcu_output) <- rownames(lambda_Transpose)

  for (i in 1:ncol(lambda_Transpose)) {
    e <- y_t * lambda_Transpose[, i]
    for (j in 1:ncol(e)) {
      Hcu_output[i, j] <- sum(e[, j])
    }
  }

  #Get the Status code and description from the LP model
  Result.StatusDesc<-sapply(Result.StatusData,function(x) dict.solveStatus[dict.solveStatus[,1] == x ,2])
  Result.status<-data.frame(code = Result.StatusData, desc = Result.StatusDesc)
  rownames(Result.status) <- dmu

  return(
    list(
      "rts" = Result.rts,
      "Orientation" = Result.orientation,
      "InputValues" =
        Result.Inputvalues,
      "OutputValues" = Result.Outputvalues,
      "Efficiency" = Result.eff,
      "Lambda" = Result.Lambda,
      "HCU_Input" = Hcu_input,
      "HCU_Output" = Hcu_output,
      "vx" = Result.vx,
      "uy" = Result.uy,
      "Free_Weights" = Result.Free_weights,
      "Model_Status" = Result.status
    )
  )
}

.Mal_Ben_Input<- function(x = x, y = y, rts ="crs",phase1_efficiency, phase = "mal", weightRestriction, include = TRUE){

  dmu <- row.names(x)

  Result.CEM <- matrix(nrow = nrow(x))
  rownames(Result.CEM) <- dmu
  colnames(Result.CEM)[1] <- c("Eff")

  Result.Lambda <- matrix(nrow = nrow(x), ncol = nrow(x))
  # Assign names to Matrix for efficiency
  rownames(Result.Lambda) <- dmu
  colnames(Result.Lambda) <- dmu

  Result.vx <-matrix(nrow = nrow(x), ncol = ncol(x))
  # Assign names to Matrix for IP Weights
  rownames(Result.vx) <- dmu
  colnames(Result.vx) <- colnames(x)

  Result.uy <-matrix(nrow = nrow(y), ncol = ncol(y))
  # Assign names to Matrix for OP Weights
  rownames(Result.uy) <- dmu
  colnames(Result.uy) <- colnames(y)

  Result.Free_weights <-matrix(nrow = nrow(x), ncol = 2)
  # Assign names to Matrix for Free Weights
  rownames(Result.Free_weights) <- dmu
  colnames(Result.Free_weights) <- c("F1","F2")

  Result.StatusData <- vector(mode="logical", nrow(x))

  sense <- ''
  equality <- ''

  # Set the Equality and equality based on Phase
  if(phase == 'mal'){
    sense <- 'min'
    equality <- '<='
  }else if (phase == 'ben'){
    sense <- 'max'
    equality <- '<='
  }

  #free sign variable
  fs<-0

  switch(rts,
         'crs' = fs <-0,
         'vrs' = fs <-2,
         fs<-2)

  #Second phase iteration over DMU efficiency scores
  #Setting input and output for each dmu
  for(d in 1 : length(dmu)){

    ip_v <- ncol(x)
    op_v <- ncol(y)

    # Set objective function for VIrtual DMU which is sum of the output variables for all DMUs
    c<-0
    for(i in 1:ncol(y)){

      if(include == TRUE){
        c[i]<-c(sum(y[,i]))
      }else if (include == FALSE) {
        c[i] <-c(sum(y[,i]) - y[d,i])
      }
    }


    # if vrs add the free variables else omit
    if(rts == 'vrs'){
      # Adding the 1 and -1 for free variable
      if(length(c) != 0){
        c[length(c) + 1] <- 1
        c[length(c) + 1] <- -1
      }
    }

    # indices for objective function
    ind_obj<-0
    for(i in 1:ncol(y)){
      ind_obj[i] <-c(ip_v + i)
    }

    # if vrs add the free variables indices else omit
    if(rts == 'vrs') {
      # indices for objective function for free variable
      if(length(ind_obj) !=0){
        ind_obj[length(ind_obj) + 1] <- ind_obj[length(ind_obj)] +1
        ind_obj[length(ind_obj) + 1] <- ind_obj[length(ind_obj)] +1
      }
    }

    lps.model <- make.lp(0, ip_v + op_v + fs)
    set.objfn(lps.model, c, indices = ind_obj)

    # Input constraint for the Virtual DMU which is sum of the Input variables for all DMUs
    ip_con_dmu<-0
    for(i in 1:ncol(x)){
      if(include == TRUE){
        ip_con_dmu[i] <- c(sum(x[,i]))
      }else if (include == FALSE) {
        ip_con_dmu[i] <- c(sum(x[,i]) - x[d,i])
      }
    }

    #indices for DMU input constraints
    ind_ip_con<-0
    for(i in 1:ncol(x)){
      ind_ip_con[i]<-i
    }

    add.constraint(lps.model, ip_con_dmu, "=", rhs = 1,ind_ip_con)

    # Multiply input values by -1 for building the constraints for all OMU's
    x_temp <- -1 *x

    # Construct constraints for all dmus excluding Rating DMU
    for(i in 1:nrow(x)){
      dmu_con<-0

      if(i != d){

        for(j in 1:ncol(x_temp)){
          dmu_con[j]<-x_temp[i,j]
        }
        for(k in 1:ncol(y)){
          dmu_con[ncol(x_temp) + k] <- y[i,k]
        }

        # If rts is vrs then add the free variable else omit
        if(rts == 'vrs') {
          #Adding free variable in the constraint
          if(length(dmu_con) !=0 ) {
            dmu_con[length(dmu_con) + 1] <- 1
            dmu_con[length(dmu_con) + 1] <- -1
          }
        }
        add.constraint(lps.model, dmu_con, equality, rhs = 0)
      }


      # Add constrains for the rating DMU
      if(i == d){
        dmu_con_Rating<-0
        for(j in 1:ncol(x_temp)){
          dmu_con_Rating[j]<- phase1_efficiency[d] * x_temp[d,j]
        }
        for(k in 1:ncol(y)){
          dmu_con_Rating[ncol(x_temp) + k] <- y[d,k]
        }

        # If rts is vrs then add the free variable else omit
        if(rts == 'vrs') {
          #Adding free variable in the constraint
          if(length(dmu_con_Rating) !=0 ) {
            dmu_con_Rating[length(dmu_con_Rating) + 1] <- 1
            dmu_con_Rating[length(dmu_con_Rating) + 1] <- -1
          }
        }
        add.constraint(lps.model, dmu_con_Rating, "=", rhs = 0)

      }
    }

    #Weight Restrictions
    if (!missing(weightRestriction)) {
      colvalues <- c(colnames(x), colnames(y))

      for (i in 1:nrow(weightRestriction))
      {
        #First constraint
        add.constraint(lps.model,
                       c(weightRestriction[i, 1], -1),
                       ">=",
                       0,
                       indices = c(
                         which(colvalues == weightRestriction[i, 2]),
                         which(colvalues == weightRestriction[i, 3])
                       ))
        #Second constraint
        # Check if upper bound exists (Length ==4)
        if (length(weightRestriction) == 4) {
          #Check if upper bound is not Inf, NA, or NaN

          if (!is.nan(weightRestriction[i, 4]) &&
              !is.na(weightRestriction[i, 4]) &&
              !is.infinite(weightRestriction[i, 4])) {
            add.constraint(
              lps.model,
              c(1, -1 * weightRestriction[i, 4]),
              "<=",
              0,
              indices = c(
                which(colvalues == weightRestriction[i, 2]),
                which(colvalues == weightRestriction[i, 3])
              )
            )
          }
        }
      }
    }

    lp.control(lps.model,sense=sense)

    Result.StatusData[d]<-solve(lps.model)


    Result.CEM[d]<-get.objective(lps.model)
    temp_lamdba <- get.dual.solution(lps.model)
    Result.Lambda[d,] <- temp_lamdba[3:(length(dmu)+2)]


    temp_weight <- get.variables(lps.model)

    ##Store the IP & OP weights in respective Matrices
    Result.vx[d,1:ncol(x)] <- temp_weight[1:ncol(x)]
    Result.uy[d,1:ncol(y)] <- temp_weight[(ncol(x)+1):(ncol(x) + ncol(y))]
    if(rts == 'vrs'){
      Result.Free_weights[d,1:2] <- temp_weight[(ncol(x) + ncol(y) + 1) : length(temp_weight)]
    }
  }

##Second Phase Cross Efficiency Calculation
  #Create the Cross Efficiency Matrix
  CrossEvaluation_Matrix <-
    matrix(nrow = length(dmu), ncol = length(dmu))
  rownames(CrossEvaluation_Matrix) <- dmu
  colnames(CrossEvaluation_Matrix) <- dmu

  #create Cross efficiency Average
  ce_ave <- matrix(nrow = 1, ncol = length(dmu))
  colnames(ce_ave) <- dmu
  rownames(ce_ave) <- "Average"

  #create Cross efficiency Max
  ceva_max <- matrix(nrow = 1, ncol = length(dmu))
  colnames(ceva_max) <- dmu
  rownames(ceva_max) <- "Max"

  #create Cross efficiency Min
  ceva_min <- matrix(nrow = 1, ncol = length(dmu))
  colnames(ceva_min) <- dmu
  rownames(ceva_min) <- "Min"


  # Input oriented CR efficiency
  for (i in 1:length(dmu)) {
    for (j in 1:length(dmu)) {
      #Calculate for IP values of DMU i * Ip weight of other DMUs J
      cr_ip <- x[i, ] * Result.vx[j, ]

      #Calculate for OP values of DMU i * OP weight of other DMUs J
      cr_op <- y[i, ] * Result.uy[j, ]

      #Calculate and store the values in Cross efficiency Matrix
      CrossEvaluation_Matrix[j, i] <- sum(cr_op) / sum(cr_ip)

    }
  }

  #Calculate cross efficiency Average
  for (i in 1:length(dmu)) {
    ce_ave[, i] <- mean(CrossEvaluation_Matrix[, i])
  }

  #Calculate cross efficiency Max
  for (i in 1:length(dmu)) {
    ceva_max[, i] <- max(CrossEvaluation_Matrix[, i])
  }

  #Calculate cross efficiency Max
  for (i in 1:length(dmu)) {
    ceva_min[, i] <- min(CrossEvaluation_Matrix[, i])
  }


  #Get the Status code and description from the LP model
  Result.StatusDesc<-sapply(Result.StatusData,function(x) dict.solveStatus[dict.solveStatus[,1] == x ,2])
  Result.status<-data.frame(code = Result.StatusData, desc = Result.StatusDesc)
  rownames(Result.status) <- dmu

  return(list("Phase2_Efficiency" = Result.CEM, "Phase2_Lambda" = Result.Lambda, "Phase2_vx" = Result.vx, "Phase2_uy" = Result.uy, "Phase2_Free_weights" = Result.Free_weights,"ceva_matrix" = CrossEvaluation_Matrix,
              "ce_ave" = ce_ave, "ceva_max" = ceva_max, "ceva_min" = ceva_min, "Phase2_Model_Status" = Result.status))
}

.Mal_Ben_Output<- function(x = x, y = y, rts ="crs",phase1_efficiency, phase = "mal", weightRestriction, include = TRUE){

  dmu <- row.names(x)

  Result.CEM <- matrix(nrow = nrow(x))
  rownames(Result.CEM) <- dmu
  colnames(Result.CEM)[1] <- c("Eff")

  Result.Lambda <- matrix(nrow = nrow(x), ncol = nrow(x))
  # Assign names to Matrix for efficiency
  rownames(Result.Lambda) <- dmu
  colnames(Result.Lambda) <- dmu

  Result.vx <-matrix(nrow = nrow(x), ncol = ncol(x))
  # Assign names to Matrix for IP Weights
  rownames(Result.vx) <- dmu
  colnames(Result.vx) <- colnames(x)

  Result.uy <-matrix(nrow = nrow(y), ncol = ncol(y))
  # Assign names to Matrix for OP Weights
  rownames(Result.uy) <- dmu
  colnames(Result.uy) <- colnames(y)

  Result.Free_weights <-matrix(nrow = nrow(x), ncol = 2)
  # Assign names to Matrix for Free Weights
  rownames(Result.Free_weights) <- dmu
  colnames(Result.Free_weights) <- c("F1","F2")

  Result.StatusData <- vector(mode="logical", nrow(x))

  sense <- ''
  equality <- ''

  # Set the Equality and equality based on Phase
  if(phase == 'mal'){
    sense <- 'min'
    equality <- '>='
  }else if (phase == 'ben'){
    sense <- 'max'
    equality <- '>='
  }

  #free sign variable
  fs<-0

  switch(rts,
         'crs' = fs <-0,
         'vrs' = fs <-2,
         fs<-2)

  #Second phase iteration over DMU efficiency scores
  #Setting input and output for each dmu
  for(d in 1:length(dmu)){

    ip_v <- ncol(x)
    op_v <- ncol(y)

    # set objective function for Virtual DMU which is sum of the output variables for all DMUs
    c<-0
    for(i in 1:ncol(x)){
      if(include == TRUE){
        c[i]<- c(sum(x[,i]))
      }else if (include == FALSE){
        c[i]<- c(sum(x[,i]) - x[d,i])
      }
    }


    # If rts is vrs then add the free variable else omit
    if(rts == 'vrs'){
      # Adding the 1 and -1 for free variable
      if(length(c) !=0){
        c[length(c) + 1] <- 1
        c[length(c) + 1] <- -1
      }
    }

    # indices for objective function
    ind_obj<-0
    for(i in 1:ncol(x)){
      ind_obj[i] <- i
    }


    # If rts is vrs then add the free variable indices else omit
    if(rts == 'vrs'){
      # indices for objective function fo r free variable
      if(length(ind_obj) != 0){
        ind_obj[length(ind_obj) + 1] <- ip_v + op_v + 1
        ind_obj[length(ind_obj) + 1] <- ip_v + op_v + 2
      }
    }

    lps.model <- make.lp(0, ip_v + op_v + fs)
    set.objfn(lps.model, c, indices = ind_obj)


    # Output constraint for the Virtual DMU which is sum of the Input variables for all DMUs
    op_con_dmu<-0
    for(i in 1:ncol(y)){
      if(include == TRUE){
        op_con_dmu[i] <- c(sum(y[,i]))
      }else if (include == FALSE) {
        op_con_dmu[i] <- c(sum(y[,i]) - y[d,i])
      }
    }

    #indices for DMU output constraints
    ind_op_con<-0
    for(i in 1:ncol(y)){
      ind_op_con[i]<- ip_v + i
    }
    add.constraint(lps.model, op_con_dmu, "=", rhs = 1,ind_op_con)

    # Multiply output values by -1 for building the constraints for all DMU's
    y_temp <- -1 * y

    # Construct constraints for all dmus excluding Rating DMU
    for(i in 1:nrow(x)){
      dmu_con<-0

      if(i != d){

        for(j in 1:ncol(x)){
          dmu_con[j] <-x[i,j]
        }
        for(k in 1:ncol(y_temp)){
          dmu_con[ncol(x) + k] <- y_temp[i,k]
        }

        # if the rts is vrs then add free variable else omit
        if(rts == 'vrs'){
          if(length(dmu_con) != 0 ) {
            dmu_con[length(dmu_con) + 1 ] <- 1
            dmu_con[length(dmu_con) + 1] <- -1
          }
        }
        add.constraint(lps.model, dmu_con, equality, rhs = 0)
      }

      # Add constrains for the rating DMU

      if(i == d){
        dmu_con_Rating<-0

        for(j in 1:ncol(x)){
          dmu_con_Rating[j] <-x[d,j]
        }

        for(k in 1:ncol(y_temp)){
          dmu_con_Rating[ncol(x) + k] <- phase1_efficiency[d] * y_temp[d,k]
        }

        # if the rts is vrs then add free variable else omit
        if(rts == 'vrs'){
          #Adding free variable in the constraint
          if(length(dmu_con_Rating) != 0 ) {
            dmu_con_Rating[length(dmu_con_Rating) + 1 ] <- 1
            dmu_con_Rating[length(dmu_con_Rating) + 1] <- -1
          }
        }
        add.constraint(lps.model, dmu_con_Rating, "=", rhs = 0)

      }
    }

    #Weight Restrictions
    if (!missing(weightRestriction)) {
      colvalues <- c(colnames(x), colnames(y))

      for (i in 1:nrow(weightRestriction))
      {
        #First constraint
        add.constraint(lps.model,
                       c(weightRestriction[i, 1], -1),
                       ">=",
                       0,
                       indices = c(
                         which(colvalues == weightRestriction[i, 2]),
                         which(colvalues == weightRestriction[i, 3])
                       ))

        #Second constraint
        # Check if upper bound exists (Length ==4)
        if (length(weightRestriction) == 4) {
          #Check if upper bound is not Inf, NA, or NaN

          if (!is.nan(weightRestriction[i, 4]) &&
              !is.na(weightRestriction[i, 4]) &&
              !is.infinite(weightRestriction[i, 4])) {
            add.constraint(
              lps.model,
              c(1, -1 * weightRestriction[i, 4]),
              "<=",
              0,
              indices = c(
                which(colvalues == weightRestriction[i, 2]),
                which(colvalues == weightRestriction[i, 3])
              )
            )
          }
        }

      }
    }

    lp.control(lps.model,sense=sense)
    Result.StatusData[d]<-solve(lps.model)


    Result.CEM[d]<- get.objective(lps.model)
    temp_lamdba <- get.dual.solution(lps.model)
    Result.Lambda[d,] <- temp_lamdba[3:(length(dmu)+2)]

    temp_weight <- get.variables(lps.model)

    #Store the IP & OP weights in respective Matrices
    Result.vx[d,1:ncol(x)] <- temp_weight[1:ncol(x)]
    Result.uy[d,1:ncol(y)] <- temp_weight[(ncol(x)+1):(ncol(x) + ncol(y))]
    if(rts == 'vrs'){
      Result.Free_weights[d,1:2] <- temp_weight[(ncol(x) + ncol(y) + 1) : length(temp_weight)]
    }
  }

  ##Second Phase Cross Efficiency Calculation
  #Create the Cross Efficiency Matrix
  CrossEvaluation_Matrix <-
    matrix(nrow = length(dmu), ncol = length(dmu))
  rownames(CrossEvaluation_Matrix) <- dmu
  colnames(CrossEvaluation_Matrix) <- dmu

  #create Cross efficiency Average
  ce_ave <- matrix(nrow = 1, ncol = length(dmu))
  colnames(ce_ave) <- dmu
  rownames(ce_ave) <- "Average"

  #create Cross efficiency Max
  ceva_max <- matrix(nrow = 1, ncol = length(dmu))
  colnames(ceva_max) <- dmu
  rownames(ceva_max) <- "Max"

  #create Cross efficiency Min
  ceva_min <- matrix(nrow = 1, ncol = length(dmu))
  colnames(ceva_min) <- dmu
  rownames(ceva_min) <- "Min"

  # Output oriented CR efficiency
  for (i in 1:length(dmu)) {
    for (j in 1:length(dmu)) {
      #Calculate for IP values of DMU i * Ip weight of other DMUs J
      cr_ip <- x[i, ] * Result.vx[j, ]

      #Calculate for OP values of DMU i * OP weight of other DMUs J
      cr_op <- y[i, ] * Result.uy[j, ]

      #Calculate and store the values in Cross efficiency Matrix
      CrossEvaluation_Matrix[j, i] <- 1 / (sum(cr_ip) / sum(cr_op))

    }
  }


  #Calculate cross efficiency Average
  for (i in 1:length(dmu)) {
    ce_ave[, i] <- mean(CrossEvaluation_Matrix[, i])
  }

  #Calculate cross efficiency Max
  for (i in 1:length(dmu)) {
    ceva_max[, i] <- max(CrossEvaluation_Matrix[, i])
  }

  #Calculate cross efficiency Max
  for (i in 1:length(dmu)) {
    ceva_min[, i] <- min(CrossEvaluation_Matrix[, i])
  }

  #Get the Status code and description from the LP model
  Result.StatusDesc<-sapply(Result.StatusData,function(x) dict.solveStatus[dict.solveStatus[,1] == x ,2])
  Result.status<-data.frame(code = Result.StatusData, desc = Result.StatusDesc)
  rownames(Result.status) <- dmu

  return(list("Phase2_Efficiency" = 1/Result.CEM, "Phase2_Lambda" = Result.Lambda, "Phase2_vx" = Result.vx, "Phase2_uy" = Result.uy, "Phase2_Free_weights" = Result.Free_weights, "ceva_matrix" = CrossEvaluation_Matrix,
              "ce_ave" = ce_ave, "ceva_max" = ceva_max, "ceva_min" = ceva_min, "Phase2_Model_Status" = Result.status))
}


.sdea_internal<-function(x, y, orientation="input", rts ="crs"){
  #Results setup
  Results <- list()

  Results$Input <- x
  Results$Output <- y
  Results$Orientation <- orientation
  Results$RTS <- rts

  Results$Efficiency <- matrix(nrow = nrow(x))
  Results$Lambda <- matrix(nrow = nrow(x), ncol = nrow(x))
  Results$StatusData <- vector(mode="logical", nrow(x))

  dmu <- row.names(x)

  # Assign names to Matrix for efficiency
  rownames(Results$Efficiency) <- dmu
  colnames(Results$Efficiency)[1] <- c("Eff")

  # Assign names to Matrix for efficiency
  rownames(Results$Lambda) <- dmu
  colnames(Results$Lambda) <- dmu


  nd<- nrow(x)
  nx <- ncol(x)
  ny <- ncol(y)

  #Model setup
  mobj <- "min"
  if(orientation=="input")
  {
    mobj <- "min"
  }else if(orientation=="output")
  {
    mobj <- "max"
  }

  #Create model
  for(k in 1:nd)
  {
    vlambda <- vbeta <- vtheta <- j <- NULL

    result <- ompr::MIPModel()
    result <- ompr::add_variable(result, vlambda[j], j=1:nd, type = "continuous", lb = 0)
    result <- ompr::add_variable(result, vtheta, type = "continuous", lb=0)
    result <- ompr::set_objective(result, vtheta, mobj)


    # input constraints
    if(orientation=="input")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * x[j,i], j = 1:nd) <=  vtheta * x[k,i])
      }
    }else if(orientation=="output")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * x[j,i], j = 1:nd) <= x[k,i])
      }
    }


    # output constraints
    if(orientation=="input")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * y[j,o], j = 1:nd) >=  y[k,o])
      }
    }else if(orientation=="output")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * y[j,o], j = 1:nd) >= vtheta * y[k,o])
      }
    }

    # rts constraints
    if(rts == "vrs")
    {
      result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j], j=1:nd) == 1)
    }

    # add lambda constraints for dmu under sonsideration as 0 for SDEA
    result<- ompr::add_constraint(result, vlambda[k] == 0)

    result <- ompr::solve_model(result, ompr.roi::with_ROI(solver = "glpk"))

    rtheta <- ompr::get_solution(result,vtheta)
    rlambda <-  ompr::get_solution(result,vlambda[j])
    rstatus <- ompr::solver_status(result)

    if(rstatus == "success"){
      Results$Efficiency[k] <- rtheta
      Results$Lambda[k, ] <- rlambda$value
      Results$StatusData[k] <- rstatus
    }else
    {
      Results$Efficiency[k] <- NaN
      Results$Lambda[k, ] <- 0
      Results$StatusData[k] <- rstatus
    }

  }

  if(orientation=="output")
  {
    Results$Efficiency <- 1/Results$Efficiency
  }


  return(Results)
}


.sdea_cook_internal<-function(x, y, orientation="input", rts ="crs", cook=TRUE){

  Results <- list()

  Results$Input <- x
  Results$Output <- y
  Results$Orientation <- orientation
  Results$RTS <- rts

  Results$Efficiency <- matrix(nrow = nrow(x))
  Results$Theta <- matrix(nrow = nrow(x))
  Results$Beta <- matrix(nrow = nrow(x))
  Results$Lambda <- matrix(nrow = nrow(x), ncol = nrow(x))
  Results$StatusData <- vector(mode="logical", nrow(x))

  dmu <- row.names(x)

  # Assign names to Matrix for efficiency
  rownames(Results$Efficiency) <- dmu
  colnames(Results$Efficiency)[1] <- c("Eff")

  # Assign names to Matrix for theta
  rownames(Results$Theta) <- dmu
  colnames(Results$Theta)[1] <- c("Theta")

  # Assign names to Matrix for Beta
  rownames(Results$Beta) <- dmu
  colnames(Results$Beta)[1] <- c("Beta")

  # Assign names to Matrix for efficiency
  rownames(Results$Lambda) <- dmu
  colnames(Results$Lambda) <- dmu


  nd<- nrow(x)
  nx <- ncol(x)
  ny <- ncol(y)


  #Model setup
  mobj <- "min"

  # Cook user defined large number
  vm <- 10^5

  #Create model
  for(k in 1:nd)
  {
    vlambda <- vbeta <- vtheta <- j <- NULL

    result <- ompr::MIPModel()
    result <- ompr::add_variable(result, vlambda[j], j=1:nd, type = "continuous", lb = 0)
    result <- ompr::add_variable(result, vbeta, type = "continuous", lb=0)
    result <- ompr::add_variable(result, vtheta, type = "continuous")
    result <- ompr::set_objective(result, (vtheta + vm * vbeta), mobj)

    # input constraints
    if(orientation=="input")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * x[j,i], j = 1:nd) <=  (1 + vtheta) * x[k,i])
      }
    }else if(orientation=="output")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * x[j,i], j = 1:nd) <=  (1 + vbeta) * x[k,i])
      }
    }


    # output constraints
    if(orientation=="input")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * y[j,o], j = 1:nd) >= (1 - vbeta) * y[k,o])
      }
    }else if(orientation=="output")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * y[j,o], j = 1:nd) >= (1 - vtheta) * y[k,o])
      }
    }

    # rts constraints
    if(rts == "vrs")
    {
      result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j], j=1:nd) == 1)
    }

    # add lambda constraints for dmu under sonsideration as 0 for SDEA
    result<- ompr::add_constraint(result, vlambda[k] == 0)

    result <- ompr::solve_model(result, ompr.roi::with_ROI(solver = "glpk"))

    rtheta <- ompr::get_solution(result,vtheta)
    rbeta <- ompr::get_solution(result,vbeta)
    rlambda <-  ompr::get_solution(result,vlambda[j])
    rstatus <- ompr::solver_status(result)


    if(orientation == "input"){
      if(rbeta == 0){
        Results$Efficiency[k] <- 1 + rtheta
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }else
      {
        Results$Efficiency[k] <- 1 + rtheta + 1/(1-rbeta)
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }
    }else if(orientation == "output"){
      if(rbeta == 0){
        Results$Efficiency[k] <- 1/(1-rtheta)
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }else
      {
        Results$Efficiency[k] <- 1 + rbeta + 1/(1-rtheta)
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }
    }

  }


  return(Results)

}



.sdea_mpi_internal<-function(x, y, orientation="input", rts ="crs", Cook=FALSE, ref.x=NULL, ref.y=NULL){

  if(ncol(x) != ncol(ref.x) && ncol(y) != ncol(ref.y)){
    stop("input or output columns for x,y ref.x amd ref.y does not match")
  }
  #Results setup
  Results <- list()

  Results$Input <- x
  Results$Output <- y
  Results$Orientation <- orientation
  Results$RTS <- rts
  Results$ref.x <- ref.x
  Results$ref.y <- ref.y

  Results$Efficiency <- matrix(nrow = nrow(x))
  Results$Lambda <- matrix(nrow = nrow(x), ncol = nrow(ref.x))
  Results$StatusData <- vector(mode="logical", nrow(x))

  dmu <- row.names(x)
  dmu_ref <- row.names(ref.x)

  # Assign names to Matrix for efficiency
  rownames(Results$Efficiency) <- dmu
  colnames(Results$Efficiency)[1] <- c("Eff")

  # Assign names to Matrix for efficiency
  rownames(Results$Lambda) <- dmu
  colnames(Results$Lambda) <- dmu_ref

  #DMUs under consideration
  nd<- nrow(x)
  nx <- ncol(x)
  ny <- ncol(y)

  #Reference DMUS
  ndr <- nrow(ref.x)
  nxr <- ncol(ref.x)
  nyr <- ncol(ref.y)

  #Model setup
  mobj <- "min"
  if(orientation=="input")
  {
    mobj <- "min"
  }else if(orientation=="output")
  {
    mobj <- "max"
  }

  #Create model
  for(k in 1:nd)
  {
    vlambda <- vbeta <- vtheta <- j <- NULL

    result <- ompr::MIPModel()
    result <- ompr::add_variable(result, vlambda[j], j=1:ndr, type = "continuous", lb = 0)
    result <- ompr::add_variable(result, vtheta, type = "continuous", lb=0)
    result <- ompr::set_objective(result, vtheta, mobj)


    # input constraints
    if(orientation=="input")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.x[j,i], j = 1:ndr) <=  vtheta * x[k,i])
      }
    }else if(orientation=="output")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.x[j,i], j = 1:ndr) <= x[k,i])
      }
    }


    # output constraints
    if(orientation=="input")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.y[j,o], j = 1:ndr) >=  y[k,o])
      }
    }else if(orientation=="output")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.y[j,o], j = 1:ndr) >= vtheta * y[k,o])
      }
    }

    # rts constraints
    if(rts == "vrs")
    {
      result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j], j=1:ndr) == 1)
    }



    result <- ompr::solve_model(result, ompr.roi::with_ROI(solver = "glpk"))

    rtheta <- ompr::get_solution(result,vtheta)
    rlambda <-  ompr::get_solution(result,vlambda[j])
    rstatus <- ompr::solver_status(result)

    if(rstatus == "success"){
      Results$Efficiency[k] <- rtheta
      Results$Lambda[k, ] <- rlambda$value
      Results$StatusData[k] <- rstatus
    }else
    {
      Results$Efficiency[k] <- NaN
      Results$Lambda[k, ] <- 0
      Results$StatusData[k] <- rstatus
    }

  }

  if(orientation=="output")
  {
    Results$Efficiency <- 1/Results$Efficiency
  }


  return(Results)
}


.sdea_cook_mpi_internal<-function(x, y, orientation="input", rts ="crs", Cook=TRUE, ref.x=NULL, ref.y=NULL){
  if(ncol(x) != ncol(ref.x) && ncol(y) != ncol(ref.y)){
    stop("input or output columns for x,y ref.x amd ref.y does not match")
  }
  Results <- list()

  Results$Input <- x
  Results$Output <- y
  Results$Orientation <- orientation
  Results$RTS <- rts
  Results$ref.x <- ref.x
  Results$ref.y <- ref.y

  Results$Efficiency <- matrix(nrow = nrow(x))
  Results$Theta <- matrix(nrow = nrow(x))
  Results$Beta <- matrix(nrow = nrow(x))
  Results$Lambda <- matrix(nrow = nrow(x), ncol = nrow(ref.x))
  Results$StatusData <- vector(mode="logical", nrow(x))

  dmu <- row.names(x)
  dmu_ref <- row.names(ref.x)

  # Assign names to Matrix for efficiency
  rownames(Results$Efficiency) <- dmu
  colnames(Results$Efficiency)[1] <- c("Eff")

  # Assign names to Matrix for theta
  rownames(Results$Theta) <- dmu
  colnames(Results$Theta)[1] <- c("Theta")

  # Assign names to Matrix for Beta
  rownames(Results$Beta) <- dmu
  colnames(Results$Beta)[1] <- c("Beta")

  # Assign names to Matrix for efficiency
  rownames(Results$Lambda) <- dmu
  colnames(Results$Lambda) <- dmu_ref

  #DMUs under consideration
  nd<- nrow(x)
  nx <- ncol(x)
  ny <- ncol(y)

  #Reference DMUS
  ndr <- nrow(ref.x)
  nxr <- ncol(ref.x)
  nyr <- ncol(ref.y)


  #Model setup
  mobj <- "min"

  # Cook user defined large number
  vm <- 10^5

  #Create model
  for(k in 1:nd)
  {
    vlambda <- vbeta <- vtheta <- j <- NULL

    result <- ompr::MIPModel()
    result <- ompr::add_variable(result, vlambda[j], j=1:ndr, type = "continuous", lb = 0)
    result <- ompr::add_variable(result, vbeta, type = "continuous", lb=0)
    result <- ompr::add_variable(result, vtheta, type = "continuous")
    result <- ompr::set_objective(result, (vtheta + vm * vbeta), mobj)

    # input constraints
    if(orientation=="input")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.x[j,i], j = 1:ndr) <=  (1 + vtheta) * x[k,i])
      }
    }else if(orientation=="output")
    {
      for (i in 1:nx) {
        result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.x[j,i], j = 1:ndr) <=  (1 + vbeta) * x[k,i])
      }
    }


    # output constraints
    if(orientation=="input")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.y[j,o], j = 1:ndr) >= (1 - vbeta) * y[k,o])
      }
    }else if(orientation=="output")
    {
      for(o in 1:ny){
        result<- ompr::add_constraint(result, ompr::sum_expr(vlambda[j] * ref.y[j,o], j = 1:ndr) >= (1 - vtheta) * y[k,o])
      }
    }

    # rts constraints
    if(rts == "vrs")
    {
      result <- ompr::add_constraint(result, ompr::sum_expr(vlambda[j], j=1:ndr) == 1)
    }

    result <- ompr::solve_model(result, ompr.roi::with_ROI(solver = "glpk"))

    rtheta <- ompr::get_solution(result,vtheta)
    rbeta <- ompr::get_solution(result,vbeta)
    rlambda <-  ompr::get_solution(result,vlambda[j])
    rstatus <- ompr::solver_status(result)


    if(orientation == "input"){
      if(rbeta == 0){
        Results$Efficiency[k] <- 1 + rtheta
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }else
      {
        Results$Efficiency[k] <- 1 + rtheta + 1/(1-rbeta)
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }
    }else if(orientation == "output"){
      if(rbeta == 0){
        Results$Efficiency[k] <- 1/(1-rtheta)
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }else
      {
        Results$Efficiency[k] <- 1 + rbeta + 1/(1-rtheta)
        Results$Theta[k] <- rtheta
        Results$Beta[k] <- rbeta
        Results$Lambda[k, ] <- rlambda$value
        Results$StatusData[k] <- rstatus
      }
    }

  }


  return(Results)
}

Try the MultiplierDEA package in your browser

Any scripts or data that you put into this service are public.

MultiplierDEA documentation built on Sept. 2, 2022, 1:06 a.m.