import pandas as pd
pd.options.display.max_columns = 100
pd.options.display.max_rows = 999
import numpy as np
import math
import warnings
warnings.filterwarnings('ignore')

import lightgbm as lgb
from sklearn.metrics import mean_squared_error

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import seaborn as sns

import datetime
import time

import gc

from google.colab import drive
drive.mount('/content/drive')

path = "/content/drive/MyDrive/tmp/7_probspace/3_taxi/"

train = pd.read_csv(path+'data/train_data.csv', parse_dates=["tpep_pickup_datetime"]).rename(columns={'tpep_pickup_datetime':'ds'})

# 前処理

# 異常値の補間
## 補間せずに、データから除外した方がいいのかもしれないというメモ
train.loc[((train["ds"]>="2017-03-12 02:00:00")&(train["ds"]<="2017-03-12 02:30:00")),list(map(str, range(79)))] = float('nan')
train.loc[((train["ds"]>="2018-03-11 02:00:00")&(train["ds"]<="2018-03-11 02:30:00")),list(map(str, range(79)))] = float('nan')
train.loc[((train["ds"]>="2019-03-10 02:00:00")&(train["ds"]<="2019-03-10 02:30:00")),list(map(str, range(79)))] = float('nan')

train = train.set_index("ds")
train = train.interpolate(method='linear', axis=0)
train = train.round().astype(int)
train = train.reset_index()

# trainに予測期間の時刻を追加
data = pd.concat([
    train,
    pd.DataFrame({
      "ds":pd.date_range(start='2019-12-01 00:00:00', end='2019-12-07 23:30:00', freq='30T')
    })
], axis=0).reset_index(drop=True)

# 縦持ち
data_df = data.copy()
data_df = pd.melt(
    data_df, id_vars=["ds"], value_vars=list(map(str, range(79))),
    var_name="id", value_name="y"
)
data_df.loc[:,"id"] = data_df["id"].astype("int64")

data_df

## 移動平均などの特徴量の追加するための関数

### 単純なlag
def create_lags(df, group_col, val_col, lags):
    lag_df = pd.DataFrame()
    for lag in lags:
        lag_df[f'lag_{lag}_{val_col}'] = df[group_col+[val_col]].copy().groupby(
            group_col
        )[val_col].transform(lambda x: x.shift(lag))
    return lag_df

# 以下の関数は現在未使用
### 移動平均
def create_rolls_ma(df, group_col, val_col, lags, rolls):
    roll_df = pd.DataFrame()
    for lag in lags:
        for roll in rolls:
            roll_df[f'ma_{lag}_{roll}_{val_col}'] = df[group_col+[val_col]].copy().groupby(
                group_col
            )[val_col].transform(lambda x: x.shift(lag).rolling(roll).mean())
    return roll_df

### 標準偏差
def create_rolls_std(df, group_col, val_col, lags, rolls):
    roll_df = pd.DataFrame()
    for lag in lags:
        for roll in rolls:
            roll_df[f'std_{lag}_{roll}_{val_col}'] = df[group_col+[val_col]].copy().groupby(
                group_col
            )[val_col].transform(lambda x: x.shift(lag).rolling(roll).std())
    return roll_df

### 移動最大
def create_rolls_max(df, group_col, val_col, lags, rolls):
    roll_df = pd.DataFrame()
    for lag in lags:
        for roll in rolls:
            roll_df[f'ma_{lag}_{roll}_{val_col}'] = df[group_col+[val_col]].copy().groupby(
                group_col
            )[val_col].transform(lambda x: x.shift(lag).rolling(roll).max())
    return roll_df

### 移動最小
def create_rolls_min(df, group_col, val_col, lags, rolls):
    roll_df = pd.DataFrame()
    for lag in lags:
        for roll in rolls:
            roll_df[f'ma_{lag}_{roll}_{val_col}'] = df[group_col+[val_col]].copy().groupby(
                group_col
            )[val_col].transform(lambda x: x.shift(lag).rolling(roll).min())
    return roll_df

# 目的変数の設定
## 一週間前（7*48行前）の同時刻との差分
data_df.loc[:,"y_diff"] =  data_df["y"] - data_df[["id", "y"]].groupby("id")["y"].transform(lambda x: x.shift(periods=48*7))

# 特徴量追加

## 時間関連
data_df.loc[:,"hour_minute"] = data_df["ds"].dt.strftime('%H%M')
data_df.loc[:,"ds_h"] = pd.to_datetime(data_df["ds"].dt.strftime('%Y-%m-%d %H')+":00:00")
data_df.loc[:,"ds_ymd"] = pd.to_datetime(data_df["ds"].dt.strftime('%Y-%m-%d ')+"00:00:00")
data_df.loc[:,"year"] = data_df["ds"].dt.year
data_df.loc[:,"month"] = data_df["ds"].dt.month
data_df.loc[:,"day"] = data_df["ds"].dt.day
data_df.loc[:,"hour"] = data_df["ds"].dt.hour
data_df.loc[:,"ampm"] = [ 1 if ampm == "AM" else 0 for ampm in data_df["ds"].dt.strftime('%p') ]
data_df.loc[:,"peak"] = 1 # 深夜4〜5時はとても少ないので、そこをフラグ化
data_df.loc[(4 <= data_df["hour"])&(data_df["hour"] < 6),"peak"] = 0

data_df.loc[:,"minute"] = data_df["ds"].dt.minute
data_df.loc[:,"quarter"] = data_df["ds"].dt.quarter
data_df.loc[:,"yobi"] = data_df["ds"].dt.weekday
data_df.loc[:,"week_number_year"] = data_df["ds"].dt.strftime('%U').astype(int) # 年における第何週か
data_df.loc[:,"week_number_month"] = [ (day - 1) // 7 + 1 for day in data_df["day"] ] # 月における第何週か

## 過去の情報を特徴量に追加

### 1,2,3,4…週間前の階差
data_df.loc[:,[
    "lag_7_y_diff","lag_8_y_diff","lag_14_y_diff","lag_15_y_diff","lag_21_y_diff","lag_28_y_diff","lag_35_y_diff","lag_42_y_diff","lag_49_y_diff","lag_364_y_diff","lag_182_y_diff",
]] = create_lags(data_df, ["id","hour_minute"], "y_diff", [7,8,14,15,21,28,35,42,49,7*52,7*26,])

### 1,2,3,4…週間前の原系列
data_df.loc[:,[
    "lag_7_y","lag_8_y","lag_14_y","lag_15_y","lag_21_y","lag_28_y","lag_35_y","lag_42_y","lag_49_y","lag_364_y","lag_182_y",
]] = create_lags(data_df, ["id","hour_minute"], "y", [7,8,14,15,21,28,35,42,49,7*52,7*26,])

### ratio
data_df.loc[:,"lag_7_y_ratio"] = data_df.loc[:,"lag_7_y_diff"] / data_df.loc[:,"lag_7_y"]
data_df.loc[:,"lag_8_y_ratio"] = data_df.loc[:,"lag_8_y_diff"] / data_df.loc[:,"lag_8_y"]
data_df.loc[:,"lag_14_y_ratio"] = data_df.loc[:,"lag_14_y_diff"] / data_df.loc[:,"lag_14_y"]
data_df.loc[:,"lag_15_y_ratio"] = data_df.loc[:,"lag_15_y_diff"] / data_df.loc[:,"lag_15_y"]
data_df.loc[:,"lag_21_y_ratio"] = data_df.loc[:,"lag_21_y_diff"] / data_df.loc[:,"lag_21_y"]
data_df.loc[:,"lag_28_y_ratio"] = data_df.loc[:,"lag_28_y_diff"] / data_df.loc[:,"lag_28_y"]
data_df.loc[:,"lag_35_y_ratio"] = data_df.loc[:,"lag_35_y_diff"] / data_df.loc[:,"lag_35_y"]
data_df.loc[:,"lag_42_y_ratio"] = data_df.loc[:,"lag_35_y_diff"] / data_df.loc[:,"lag_35_y"]
data_df.loc[:,"lag_49_y_ratio"] = data_df.loc[:,"lag_35_y_diff"] / data_df.loc[:,"lag_35_y"]
data_df.loc[:,"lag_364_y_ratio"] = data_df.loc[:,"lag_364_y_diff"] / data_df.loc[:,"lag_364_y"]
data_df.loc[:,"lag_182_y_ratio"] = data_df.loc[:,"lag_182_y_diff"] / data_df.loc[:,"lag_182_y"]

### 階差列のcumsum
data_df.loc[:,"lag_7_14_diff_sum"] = data_df["lag_7_y_diff"] + data_df["lag_14_y_diff"]
data_df.loc[:,"lag_7_14_21_diff_sum"] = data_df["lag_7_14_diff_sum"] + data_df["lag_21_y_diff"]
data_df.loc[:,"lag_7_14_21_28_diff_sum"] = data_df["lag_7_14_21_diff_sum"] + data_df["lag_28_y_diff"]
data_df.loc[:,"lag_7_14_21_28_35_diff_sum"] = data_df["lag_7_14_21_28_diff_sum"] + data_df["lag_35_y_diff"]
data_df.loc[:,"lag_7_14_21_28_35_42_diff_sum"] = data_df["lag_7_14_21_28_35_diff_sum"] + data_df["lag_42_y_diff"]
data_df.loc[:,"lag_7_14_21_28_35_42_49_diff_sum"] = data_df["lag_7_14_21_28_35_42_diff_sum"] + data_df["lag_49_y_diff"]

### 原系列系列のcumsum
data_df.loc[:,"lag_7_14_sum"] = data_df["lag_7_y"] + data_df["lag_14_y"]
data_df.loc[:,"lag_7_14_21_sum"] = data_df["lag_7_14_sum"] + data_df["lag_21_y"]
data_df.loc[:,"lag_7_14_21_28_sum"] = data_df["lag_7_14_21_sum"]  + data_df["lag_28_y"]
data_df.loc[:,"lag_7_14_21_28_35_sum"] = data_df["lag_7_14_21_28_sum"] + data_df["lag_35_y"]
data_df.loc[:,"lag_7_14_21_28_35_42_sum"] = data_df["lag_7_14_21_28_35_sum"] + data_df["lag_42_y"]
data_df.loc[:,"lag_7_14_21_28_35_42_49_sum"] = data_df["lag_7_14_21_28_35_42_sum"] + data_df["lag_49_y"]


# 確認のため、id==0の00:00、00:30だけ表示
data_df[((data_df["id"]==0)&((data_df["hour_minute"]=="0000")|(data_df["hour_minute"]=="0030")))].head(100)

# 手動でクロスバリデーション期間の設定
## valid期間は一週間で、それをずらしていく
t_start = datetime.datetime(2017, 7, 1)
t_end =   datetime.datetime(2019, 12, 1)

num_valid = 20

list_term = []
for i in range(num_valid):
  list_term += [
      [
          t_start + datetime.timedelta(days=7*i), # trainの最初の日
          t_end + datetime.timedelta(days=-7*(num_valid-i)),   # trainの最後の次の日 = validの最初の日
          t_end + datetime.timedelta(days=-7*(num_valid-i-1)), # validの最後の次の日
      ]
]

# testデータ作成
target = "y_diff"

feature_c = [
    "id","month","day","hour","yobi","week_number_year","peak","ampm","minute","quarter","week_number_month",
]
feature_r = [
    ## diff
    "lag_7_y_diff","lag_8_y_diff","lag_14_y_diff","lag_15_y_diff","lag_21_y_diff","lag_28_y_diff","lag_35_y_diff","lag_42_y_diff","lag_49_y_diff","lag_364_y_diff","lag_182_y_diff",
    ## y
    "lag_7_y","lag_8_y","lag_14_y","lag_15_y","lag_21_y","lag_28_y","lag_35_y","lag_42_y","lag_49_y","lag_364_y","lag_182_y",
    ## ratio
    "lag_7_y_ratio","lag_8_y_ratio","lag_14_y_ratio","lag_15_y_ratio","lag_21_y_ratio","lag_28_y_ratio","lag_35_y_ratio","lag_42_y_ratio","lag_49_y_ratio","lag_364_y_ratio","lag_182_y_ratio",
    ## sum
    "lag_7_14_diff_sum","lag_7_14_21_diff_sum","lag_7_14_21_28_diff_sum","lag_7_14_21_28_35_diff_sum","lag_7_14_21_28_35_42_diff_sum","lag_7_14_21_28_35_42_49_diff_sum",
    "lag_7_14_sum","lag_7_14_21_sum","lag_7_14_21_28_sum","lag_7_14_21_28_35_sum","lag_7_14_21_28_35_42_sum","lag_7_14_21_28_35_42_49_sum",
]
feature = feature_r + feature_c

test_X = data_df[data_df["ds"] >= t_end][["ds"]+feature].copy().reset_index(drop=True)
test_X.head()

start = time.time()
print(datetime.datetime.now())

params = {
    "objective" : "regression_l2",
    "metric" : "rmse",
    'verbosity': -1,
    'boosting_type': 'gbdt',
    "max_depth": 10,
    "seed" : 810, # 鳩
}

valid_scores = []
models = []
preds = []
for i,term in enumerate(list_term):
    print(i)
    print(datetime.datetime.now())
    X_train = data_df[(term[0] <= data_df["ds"])&(data_df["ds"] < term[1])][["ds"]+feature].copy().reset_index(drop=True)
    X_valid = data_df[(term[1] <= data_df["ds"])&(data_df["ds"] < term[2])][["ds"]+feature].copy().reset_index(drop=True)
    y_train = data_df[(term[0] <= data_df["ds"])&(data_df["ds"] < term[1])][[target]].copy().reset_index(drop=True)
    y_valid = data_df[(term[1] <= data_df["ds"])&(data_df["ds"] < term[2])][[target]].copy().reset_index(drop=True)

    train_set = lgb.Dataset(X_train.drop(columns=["ds"]), y_train, categorical_feature=feature_c)
    valid_set = lgb.Dataset(X_valid.drop(columns=["ds"]), y_valid, categorical_feature=feature_c)

    # optunaやろうかやるまいか。。あとでやる
    best_booster = lgb.train(
        params=params,
        num_boost_round=100,
        train_set=train_set,
        valid_sets=[valid_set],
        callbacks=[
            lgb.early_stopping(stopping_rounds=10,verbose=True),
            lgb.log_evaluation(0)
        ],
    )
    y_valid_pred = best_booster.predict(X_valid.drop(columns=["ds"]))
    score = np.sqrt(mean_squared_error(y_valid, y_valid_pred))
    print(score)

    valid_scores.append(score)
    models.append(best_booster)

    X_valid.loc[:,target] = y_valid
    X_valid.loc[:,target+"_pred"] = y_valid_pred
    preds += [X_valid[["ds","id",target,target+"_pred"]]]

df_pred = pd.concat(preds, axis=0)
print(datetime.datetime.now())
print(time.time() - start)

print(valid_scores)
print(f'cv: {np.mean(valid_scores):.4f} ± {np.std(valid_scores):.4f}')

importance = pd.DataFrame({
    "feature":feature# + feature_p,
})
for i in range(num_valid):
  importance.loc[:,"imp_"+str(i)] = models[i].feature_importance(importance_type='gain')

importance.loc[:,"imp"] = importance[["imp_"+str(i) for i in range(num_valid)]].mean(axis=1)
importance = importance.sort_values( ["imp"], ascending=False ).reset_index(drop=True)

# 可視化用に縦持ち
importance_p = pd.melt(
    importance, id_vars=["feature"], value_vars=["imp_"+str(i) for i in range(num_valid-2)],
    var_name="cv", value_name="imp"
)

fig, ax = plt.subplots(1, 1, figsize=(20,5))
sns.boxplot(x='feature', y='imp', data=importance_p, showfliers=True, ax=ax)
sns.stripplot(x='feature', y='imp', data=importance_p, jitter=True, color='black', ax=ax)
plt.xticks(rotation=90)
plt.show()
plt.close()
plt.clf()

df_p_test = pd.DataFrame({
    "ds":test_X["ds"],
    "id":test_X["id"],
})
for i in range(num_valid):
  df_p_test.loc[:,"p_"+str(i)] = models[i].predict(test_X.drop(columns=["ds"]))
df_p_test.loc[:,target] = df_p_test[["p_"+str(i) for i in range(num_valid)]].mean(axis=1)

df_p_test

# 原系列に戻す関数
## 超絶汚い希ガス・・・
def reverse_diff2org(data_df, d, colname):
  d = d[["ds","id",colname]].copy()

  # cumsumするための初期値
  init_s = d["ds"].min() + datetime.timedelta(days=-7)
  init_e = d["ds"].min()
  d_init = data_df.loc[((init_s<=data_df["ds"])&(data_df["ds"]<init_e)), ["ds","id","y"]].rename(columns={'y':colname})
  # 原系列に戻す
  d_org = pd.concat([d_init, d], axis=0).reset_index(drop=True)
  d_org.loc[:,"hour_minute"] = d_org["ds"].dt.strftime('%H%M')
  d_org.loc[:,"yobi"] = d_org["ds"].dt.weekday

  d_org.loc[:,"y_pred"] = d_org[["id","yobi","hour_minute",colname]].groupby(["id","yobi","hour_minute"]).cumsum().agg(np.ceil) # 切り上げ
  d_org.loc[d_org["y_pred"]<0,"y_pred"] = 0
  d_org.loc[:,"y_pred"] = d_org["y_pred"].astype(int)
  # 初期値を削除
  d_org = d_org[d_org["ds"]>=init_e]
  # 正解の値をマージ（validの評価用）
  d_org = pd.merge(d_org, data_df[["ds", "id", "y"]], how="left", on=["ds", "id"])

  return d_org

df_p_test_org = reverse_diff2org(data_df, d=df_p_test[["ds","id","y_diff"]], colname="y_diff")
df_p_test_org

# 予測結果の確認
plot_d_y = data_df[data_df["ds"]>="2019-11-01"]
plot_d_p = df_p_test_org

fig, ax = plt.subplots(79, 1, figsize=(20, 79*4))

for i in range(79):
    ax[i].plot(plot_d_y[plot_d_y["id"]==i]['ds'], plot_d_y[plot_d_y["id"]==i]["y"])
    ax[i].plot(plot_d_p[plot_d_p["id"]==i]['ds'], plot_d_p[plot_d_p["id"]==i]["y_pred"])
    ax[i].set_title(str(i))

plt.show()
plt.close()
plt.clf()

# 画像は添付をご確認ください

# sub用データ作成（たて→よこ）
df_sub = df_p_test_org.copy()
df_sub.loc[:,"id"] = df_sub["id"].astype(str)
df_sub = df_sub.pivot(
    index="ds",
    columns="id",
    values="y_pred" # , y_pred, y_maze
).reset_index()[["ds"]+list(map(str, range(79)))]
df_sub.columns = [["tpep_pickup_datetime"]+list(map(str, range(79)))]
df_sub.head()