# 計算每類使用者類別及使用者數
Type <- table(data_FB[, 3])

# 方法1：基礎繪圖
p <- barplot(Type, space = 0, ylim = c(0, 30), col = rainbow(7), xaxt = "n",
             ylab = "計數", main = "竊電使用者用電類別分佈分析")
df <- data.frame(Type)
axis(1, p, df$Var1, las = 2)
text(p, Type, labels = Type,pos = 3)  # 新增數值標籤
# xaxt="n"設定是否顯示x軸資訊,axes=F不顯示座標軸
# 顏色可以用heat.colors,terrain.colors,cm.colors等
 

# 方法2：餅圖
Type <- data.frame(Type)
pct <- round(Type$Freq / sum(Type$Freq) * 100, 1)
lbls <- paste0(Type$Var1, pct, "%")
pie(Type$Freq, labels = lbls)  # 普通餅圖
library(plotrix)  # 3D餅圖
pie3D(Type$Freq, labels = lbls, main = "竊電使用者用電類別分佈",
      labelrad = 1.4, start = 3)
pie3D(Type$Freq, labels = lbls, explode = 0.1, radius = 1) 
# radius半徑，explode分離度
 

# ----------------------------正常使用者用電量趨勢分析---------------------------

Regular <- read.csv("正常用電量資料.csv",header=T)

# 基礎繪圖
plot(1:59, Regular[, 2], type = "l", col = "blue", 
     main = "正常使用者的用電量趨勢")  # 主標題，x,y軸標題

# ---------------------------竊漏電使用者用電量趨勢分析--------------------------

Unusual <- read.csv("竊電用電量資料.csv", header = T)

# 方法1：基礎繪圖
plot(Unusual[, 2], type = "l", col = "blue")  # 設定主標題，x,y軸標題

 

# 方法2：在一張圖上對比
plot(Regular[,2], col = "blue", lty = 1, type = "l", 
     main = "正常使用者與竊電使用者的用電趨勢比較", ylim = c(0,8000), 
     ylab = "用電量", xlab = "", xaxt = "n")  # axes = F

axis(1, at = 1:59, Regular[, 1], las = 2)  # 設定x軸
lines(Unusual[,2], col = "red", lty = 2, type = "l")  # 新增竊電資料
legend("bottomleft", legend = c("正常使用者","竊電使用者"), 
       lty = 1:2, col = c("blue", "red"))  # 新增圖例

# ---------------------------------求斜率--------------------------------------
# 讀取資料
Power <- read.csv("使用者日用電量.csv")

# 設定一個向量k存放求出的每天的斜率值
k <- rep(0, nrow(Power))
Down <- Power$日電量

for (i in 1:nrow(Power)) {  # 迴圈所有天，求出每天的前後5天共計11天的平均斜率k
  if (i <= 5) {
    l <- 1:(i + 5)  # 前面不足5天時求平均斜率用到的天數
  }
  if (i >5 & i < (nrow(Power) - 5)) {
    l <- (i - 5):(i + 5)  # 前後均滿足5天時的求平均斜率用到的天數
  }
  if (i >= (nrow(Power) - 5)) {
    l <- (i-5):nrow(Power)  # 後面不足5天時求平均斜率用到的天數
  }
  k[i] <- cov(Down[l], l) / var(l)  # 計算第i天的斜率，公式類似求協方差除方差
}


# -----------------------------標記用電量趨勢----------------------------------

Decrease <- rep(0, nrow(Power))  # 設定變數D，用於存放電量趨勢標記1或0

for (i in 2:nrow(Power)) {  # 從第二天開始迴圈，求出所有天的電量趨勢標記
  if (k[i] < k[i - 1]) {
    Decrease[i] <- 1  # 當天比前一天斜率低，標記為1
  } 
  if (k[i] >= k[i - 1]) { 
    Decrease[i] <- 0  # 當天高於或等於前一天斜率，標記為0
  }
}

# --------------------------統計11天內的趨勢下降次數---------------------------

Total <- rep(0, nrow(Power)) #設定變數T，存放對第i天前4後5共計10天的電量趨勢彙總值

for (i in 1:nrow(Power)) {  # 迴圈所有天，求出每天前4後5共計10天的趨勢彙總
  if (i < 5) {
    m <- 1:(i + 5)  # 前面不足5天要彙總的天數
  } 
  if (i >= 5 & i <= nrow(Power) - 5) {
    m <- (i - 4):(i + 5)  # 滿足前後5天時研彙總的天數
  }  
  if (i > nrow(Power) - 5) {
    m <- (i - 4):nrow(Power)  # 後面不足5天時用到的天數
  }  
  Total[i] <- sum(Decrease[m])
}

# 讀取資料
data_alarm <- read.csv("告警.csv")
data <- read.csv("使用者.csv")

# 構造ID&date屬性
data_alarm$ID_date <- paste(data_alarm[, 1], data_alarm[, 2])
data$ID_date <- paste(data[, 1], data[, 2])

# 統計使用者每天的告警次數
D <- data.frame(matrix(0, nrow(data), nrow(data_alarm)))

for (i in (1:nrow(data))) {
  for (k in (1:nrow(data_alarm))) {
    if (data$ID_date[i] == data_alarm$ID_date[k]) {
      D[i, k] <- 1
    } else {
      D[i, k] <- 0}
  }
}

D$sum <- apply(D, 1, sum)  # 按行計算總和
data$alarm_ind <- D$sum
data <- data[, c(1, 2, 6)]  # 去除不需要的ID，日期和告警次數

library(XLConnect)
missing_data <- XLConnect::readWorksheetFromFile(file = "missing_data.xls",
                                                 sheet = 1, header = FALSE)

lagrange <- function(x, xi, yi) {
  n <- length(xi)
  lage <- 0
  for (i in 1:n) {
    li <- 1
    for (j in 1:n) {
      if (i != j) {
        li <- li * (x - xi[j]) / (xi[i] - xi[j])
      }
    }
    lage <- li * yi[i] + lage
  }
  return(lage)
}

missdata <- missing_data

for (k in 1:3) {
  x <- which(is.na(missing_data[, k]))
  x1 <- c(0, x)
  x2 <- c(x, nrow(missing_data))
  x12 <- x2 - x1 - 1
  xx1 <- x12[1:(length(x12) - 1)]  # 缺失值前面的行數
  xx2 <- x12[2:(length(x12))]  # 缺失值後面的行數
  
  j <- 1
  for (m in x) {
    if (xx1[j] >= 5) {  # 空值前的判斷
      xi <- (m - 5):(m - 1)
    } else {
      xi <- (m - xx1[j]):(m - 1)
    }
    if (xx2[j] >= 5) {  # 空值後的判斷
      xi <- c(xi, (m + 1):(m + 5))
    } else {
      xi <- c(xi, (m + 1):(m + xx2[j]))
    }
    yi <- missing_data[xi, k]
    missdata[m, k] <- lagrange(m, xi, yi) 
    print(c(m, missdata[m, k]))
    j <- j + 1
  }
}

# 讀取資料
data_loss <- read.csv("線損.csv")

# 構造線損
data_loss$日線損率 <- (data_loss[, 3] - data_loss[, 4]) / data_loss[, 3]


# 便於程式碼呼叫，將日線損率資料賦予變數v
V <- data_loss$日線損率

# Vb為當天與後5天共6天的線損率平均值
# Vf為當天與前5天共6天的線損率平均值
n <- nrow(data_loss)
Vb <- rep(0, n) 
Vf <- rep(0, n) 
E <- rep(0, n)   # 設定變數E，存放線損指標
for (i in 1:n) {  # 迴圈所有天，求出每天的線損指標
  if (i <= 5) {  # 前面不足5天的情況
    Vb[i] <- mean(V[i:(i + 5)])
    Vf[i] <- mean(V[1:i])
  }
  if (i > 5 & i < n - 5) {  # 前後均滿足5天的情況
    Vb[i] <- mean(V[i:(i + 5)])
    Vf[i] <- mean(V[(i - 5):i])
  }
  if (i >= n - 5) {  # 後面不足5天的情況
    Vb[i] <- mean(V[i:n])
    Vf[i] <- mean(V[(i - 5):i])
  }
  if ((Vb[i] - Vf[i]) / Vf[i] > 0.01) {
    E[i] <- 1  # Vb比Vf的增長率判斷，並標記
  }
  if ((Vb[i]-Vf[i]) / Vf[i] <= 0.01) {
    E[i] <- 0
  }
}