library(lpSolve)
library(lpSolveAPI)

##정해진 위치에 -1, 나머지는 0인 벡터를 생성해주는 함수 #lpSolveAPI를 편하게 쓰기 위해 짜는 것임

mvector <- function(n, x)
{
  m <- rep(0, n)
  m[x] = -1
  return(m)
}

##crs 투입 함수

crs_input <- function(data, i_start, i_end, o_start, o_end, ac)
#data     분석대상데이터
#i_start  투입변수가 시작되는 열 번호
#i_end    투입변수가 종료되는 열 번호
#o_start  산출변수가 시작되는 열 번호
#o_end    산출변수가 종료되는 열 번호
#ac       아르키메데스 상수

{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    #선형계획을 생성한다. 
    #목적함수 및 제약식에 사용되는 변수 숫자만큼의 열을 생성해야 한다. 
    #1은 효율성점수 theta에 할당되는 열의 숫자이다. 
    #nrow(data)는 lambda에 해당되는 dmu의 숫자만큼의 열을 생성해야 한다. 
    #i_end - i_start + 1는 투입변수의 갯수로서, 제약식에 사용되는 투입변수의 slack을 의미한다. 
    #o_end - o_start + 1는 산출변수의 갯수로서, 제약식에 사용되는 산출변수의 slack을 의미한다. 
    
    lp.control(dea1, sense = "min")
    #선형계획의 목적함수를 최대화 또는 최소화할것인지를 결정한다. 
    
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    #선형계획에서 사용되는 각 열의 변수명을 입력한다. 
    #효율성점수를 나타내는 첫 열은 theta
    #참조하는 dmu의 가중치인 lambda는 "l"과 dmu 번호의 조합으로 나타낸다(l1, l2, l3, ...)
    #투입변수의 slack은 "is"와 투입요소 번호의 조합으로 나타낸다(i1, i2, i3, ...)
    #산출변수의 slack은 "os"와 산출요소 번호의 조합으로 나타낸다(o1, o2, o3, ...)
    
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(-ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    #목적함수를 만든다. 
    #theta에 해당하는 첫 열에 1을 할당한다.
    #rep(0, nrow(data)) 사용되지 않는 열(변수)인 lambda에는 0을 할당한다 
    #아르키메데스상수와 투입 및 산출요소의 합과 slack의 곱을 theat에 빼는 것이 목적함수이다.  
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(data[i,j], -data[,j], mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", 0)
    }
    #투입요소 갯수만큼의 제약식을 추가한다. 
    #제약식은 좌변에 변수항을, 우변에 0을 나타내도록 고친다
    #제약식: theta*투입요소 - lambda*투입요소 - 투입요소 slack = 0
    #data[i,j]  제약식 중 theta와 투입요소의 곱을 의미한다. 
    #-data[,j]  제약식 중 lambda와 투입요소의 곱을 의미한다. 
    #mvector((i_end - i_start + 1), (j - i_start + 1))  제약식 중 투입요소 slack의 음수값을 의미한다. 
    #rep(0, o_end - o_start + 1)  산출요소 slack은 사용되지 않는다. 
    
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(0, data[,k], rep(0, i_end - i_start + 1), mvector((o_end - o_start + 1), (k - o_start + 1))), "=", data[i,k])
    }
    #산출요소 갯수만큼의 제약식을 추가한다.
    #제약식은 좌변에 변수항을, 우변에 산출요소를 나타내도록 고친다
    #제약식: lambda*산출요소 - 산출요소 slack = 산출요소
    #0        제약식에 theta는 사용되지 않는다. 
    #data[,k] 제약식 중 lambda와 산출요소의 곱을 의미한다. 
    #rep(0, i_end - istart + 1) 제약식 중 투입요소 slack은 사용되지 않는다. 
    #mvector((o_end - o_start + 1), (k - o_start + 1))) 제약식 중 산출요소 slack의 음수값을 의미한다. 
    #data[i, k] 제약식 중 산출요소를 의미한다. 
    
    solve(dea1)
    #최적화 작업을 수행한다(목적함수 최소화)
    
    output <- rbind(output, get.variables(dea1))
    #최적화 결과를 저장한다. 
  }
  colnames(output) <- colnames(dea1)
  #열 이름을 저장한다. 
  
  output <- as.data.frame(output)
  return(output)
  #전체 분석결과를 리턴한다. 
}

##vrs 투입 함수

vrs_input <- function(data, i_start, i_end, o_start, o_end, ac)
{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    lp.control(dea1, sense = "min")
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(-ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(data[i,j], -data[,j], mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", 0)
    }
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(0, data[,k], rep(0, i_end - i_start + 1), mvector((o_end - o_start + 1), (k - o_start + 1))), "=", data[i,k])
    }

    add.constraint(dea1, c(0, rep(1, nrow(data)), rep(0, o_end - i_start + 1)), "=", 1)
    
    solve(dea1)
    output <- rbind(output, get.variables(dea1))
    colnames(output) <- colnames(dea1)
    }
  output <- as.data.frame(output)
  return(output)
}

##crs 산출 함수

crs_output <- function(data, i_start, i_end, o_start, o_end, ac)

{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    
    lp.control(dea1, sense = "max")
    
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(0, data[,j], -mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", data[i,j])
    }
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(data[i,k], -data[,k], rep(0, i_end - i_start + 1), -mvector((o_end - o_start + 1), (k - o_start + 1))), "=", 0)
    }
    
    solve(dea1)
    output <- rbind(output, get.variables(dea1))
    colnames(output) <- colnames(dea1)
  }
  output <- as.data.frame(output)
  output$theta2 <- 1/output$theta
  output <- output[, c(ncol(output), 1:ncol(output)-1)]
  return(output)
}

##vrs 산출 함수

vrs_output <- function(data, i_start, i_end, o_start, o_end, ac)
  
{
  output <- c()
  
  for(i in 1:nrow(data))
  {
    dea1 <- make.lp(0, 1 + nrow(data) + (i_end - i_start + 1) + (o_end - o_start + 1)) 
    
    lp.control(dea1, sense = "max")
    
    colnames(dea1) <- c("theta", paste("l", 1:nrow(data), sep = ""), 
                        paste("is", 1:(i_end - i_start + 1), sep = ""), 
                        paste("os", 1:(o_end - o_start + 1), sep = ""))
    set.objfn(dea1, c(1, rep(0, nrow(data)), rep(ac, (i_end - i_start + 1) + (o_end - o_start + 1))))
    
    
    for (j in i_start:i_end)
    {
      add.constraint(dea1, c(0, data[,j], -mvector((i_end - i_start + 1), (j - i_start + 1)), rep(0, o_end - o_start + 1)), "=", data[i,j])
    }
    
    for (k in o_start:o_end)
    {
      add.constraint(dea1, c(data[i,k], -data[,k], rep(0, i_end - i_start + 1), -mvector((o_end - o_start + 1), (k - o_start + 1))), "=", 0)
    }
    
    add.constraint(dea1, c(0, rep(1, nrow(data)), rep(0, o_end - i_start + 1)), "=", 1)
    
    solve(dea1)
    output <- rbind(output, get.variables(dea1))
    colnames(output) <- colnames(dea1)
  }
  output <- as.data.frame(output)
  output$theta2 <- 1/output$theta
  output <- output[, c(ncol(output), 1:ncol(output)-1)]
  return(output)
}

##분석예시

i1 <- c(8.939   ,   8.625   ,   10.813  ,   10.638  ,   6.24    ,   4.719)
i2 <- c(64.3    ,   99  ,   99.6    ,   96  ,   96.2    ,   79.9)
o1 <- c(25.2    ,   28.2    ,   29.4    ,   26.4    ,   27.2    ,   25.5)
o2 <- c(223 ,   287 ,   317 ,   291 ,   295 ,   222)

school <- data.frame(i1, i2, o1, o2)

crs_input(school, 1, 2, 3, 4, 10^-20)
##       theta        l1 l2 l3 l4        l5 l6 is1 is2      os1 os2
## 1 1.0000000 1.0000000  0  0  0 0.0000000  0   0   0 0.000000   0
## 2 0.9096132 0.4203319  0  0  0 0.6551389  0   0   0 0.212143   0
## 3 0.9635345 0.8795257  0  0  0 0.4097145  0   0   0 3.908281   0
## 4 0.9143053 0.8458087  0  0  0 0.3470666  0   0   0 4.354592   0
## 5 1.0000000 0.0000000  0  0  0 1.0000000  0   0   0 0.000000   0
## 6 1.0000000 0.0000000  0  0  0 0.0000000  1   0   0 0.000000   0
vrs_input(school, 1, 2, 3, 4, 10^-20)
##       theta        l1 l2        l3 l4         l5         l6 is1 is2      os1
## 1 1.0000000 1.0000000  0 0.0000000  0 0.00000000 0.00000000   0   0 0.000000
## 2 0.9788329 0.0000000  0 0.5020750  0 0.43641614 0.06150884   0   0 0.000000
## 3 1.0000000 0.0000000  0 1.0000000  0 0.00000000 0.00000000   0   0 0.000000
## 4 0.9394746 0.2595833  0 0.6677271  0 0.07268964 0.00000000   0   0 1.749833
## 5 1.0000000 0.0000000  0 0.0000000  0 1.00000000 0.00000000   0   0 0.000000
## 6 1.0000000 0.0000000  0 0.0000000  0 0.00000000 1.00000000   0   0 0.000000
##       os2
## 1  0.0000
## 2 14.5555
## 3  0.0000
## 4  0.0000
## 5  0.0000
## 6  0.0000
crs_output(school, 1, 2, 3, 4, 10^-20)
##      theta2    theta        l1 l2 l3 l4        l5 l6 is1 is2       os1 os2
## 1 1.0000000 1.000000 1.0000000  0  0  0 0.0000000  0   0   0 0.0000000   0
## 2 0.9096132 1.099368 0.4620996  0  0  0 0.7202390  0   0   0 0.2332233   0
## 3 0.9635345 1.037846 0.9128118  0  0  0 0.4252204  0   0   0 4.0561922   0
## 4 0.9143053 1.093727 0.9250835  0  0  0 0.3795960  0   0   0 4.7627330   0
## 5 1.0000000 1.000000 0.0000000  0  0  0 1.0000000  0   0   0 0.0000000   0
## 6 1.0000000 1.000000 0.0000000  0  0  0 0.0000000  1   0   0 0.0000000   0
vrs_output(school, 1, 2, 3, 4, 10^-20)
##      theta2    theta       l1 l2        l3 l4        l5 l6        is1      is2
## 1 1.0000000 1.000000 1.000000  0 0.0000000  0 0.0000000  0 0.00000000 0.000000
## 2 0.9948007 1.005226 0.000000  0 0.5215395  0 0.4784605  0 0.00000000 1.026766
## 3 1.0000000 1.000000 0.000000  0 1.0000000  0 0.0000000  0 0.00000000 0.000000
## 4 0.9466074 1.056404 0.101983  0 0.8980170  0 0.0000000  0 0.01611615 0.000000
## 5 1.0000000 1.000000 0.000000  0 0.0000000  0 1.0000000  0 0.00000000 0.000000
## 6 1.0000000 1.000000 0.000000  0 0.0000000  0 0.0000000  1 0.00000000 0.000000
##        os1      os2
## 1 0.000000  0.00000
## 2 0.000000 17.97387
## 3 0.000000  0.00000
## 4 1.082603  0.00000
## 5 0.000000  0.00000
## 6 0.000000  0.00000