#!/usr/bin/env python
# coding: utf-8

# # Import Libraries

# In[1]:


import pandas as pd
import numpy as np

from sklearn.preprocessing import StandardScaler
import statsmodels.api as sm

import random
import time
import os

import datetime

import xgboost
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import explained_variance_score
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV

from statsmodels.stats.outliers_influence import variance_inflation_factor


import matplotlib.pyplot as plt
import seaborn as sns

get_ipython().run_line_magic('matplotlib', 'inline')

import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.family'] = 'NanumGothic'
plt.rcParams['font.size'] = 13


# # Fixed Random Seed

# In[2]:


def seed_everything(seed):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)

seed_everything(42) # Seed 고정


# # Load Dataset

# In[3]:


train = pd.read_csv('data/train.csv')
test = pd.read_csv('data/test.csv')


# In[4]:


train.head(5)


# In[5]:


test.head(5)


# # At a glance

# Dataset Info.
# 
# #### train.csv
# 57920개의 데이터  
# ID : 샘플 별 고유 id  
# 월: 데이터가 기록된 달을 나타냅니다.  
# 일: 데이터가 기록된 날짜를 나타냅니다.  
# 측정 시간대: 데이터가 측정된 시간대를 나타냅니다. 오전, 오후, 저녁, 새벽으로 구분되어 있습니다.  
# 섭씨 온도 (° ⁣C)  
# 절대 온도 (K)  
# 이슬점 온도 (° ⁣C)  
# 상대 습도 (%)  
# 대기압 (mbar)  
# 포화 증기압 (mbar)  
# 실제 증기압 (mbar)  
# 증기압 부족량 (mbar)  
# 수증기 함량 (g/kg): 공기 1 kg당 수증기의 질량을 그램(g) 단위로 나타냅니다.  
# 공기 밀도 (g/m**3): 1 m³의 부피에 들어있는 공기의 질량을 그램(g) 단위로 나타냅니다.  
# 풍향 (deg): 바람의 향하는 방향을 각도(degree)로 나타냅니다.  
# 풍속 (m/s)

# train 데이터프레임은 총 26129개의 행(RangeIndex)과 16개의 열(Data columns)로 구성되어 있습니다.
# 
# 수치형 데이터가 12개(float64)와 2개(int64)의 컬럼으로 구성되어 있으며, object 타입의 컬럼도 2개 있습니다.
# 따라서 'ID'와 '측정 시간대' 컬럼이 문자열 데이터라는 것을 알 수 있습니다.

# In[6]:


train.info()


# In[7]:


train.describe()


# # Background Knowledge

# #### 풍속과 관련 있는 변수
# 
# 1. 대기압 (mbar): 기압은 대기의 무게로 인해 발생하므로, 대기압은 기압에 직접적인 영향을 미칩니다.
# 
# 2. 이슬점 온도 (° ⁣C): 이슬점 온도는 공기가 포화 상태로 수증기를 포함하는데 필요한 온도를 의미하며, 공기가 포화되면 대기압에 영향을 줄 수 있습니다.
# 
# 3. 상대 습도 (%): 상대 습도는 공기가 포화 상태에 얼마나 가까운지를 나타내는데, 수증기가 포화되면 기압에 영향을 미칠 수 있습니다.
# 
# 4. 수증기 함량 (g/kg): 공기 내의 수증기의 양은 대기압과 관련하여 영향을 미칠 수 있습니다.
# 
# 5. 섭씨 온도 (° ⁣C) 또는 절대 온도 (K): 온도의 변화는 대기의 열적 팽창과 함께 기압에 영향을 미칠 수 있습니다.
# 
# 6. 증기압 부족량 (mbar): 공기 중의 수증기가 포화 상태에 미치지 못하는 압력 차이로, 기압에 영향을 줄 수 있습니다.
# 
# 7. 공기 밀도 (g/m**3): 공기의 밀도는 공기의 질량과 부피에 영향을 미치며, 이는 일부 경우에 기압에 영향을 미칠 수 있습니다.
# 
# 8. 측정 시간대: 측정 시간대는 기압에 직접적인 영향은 없지만, 기상 조건과 관련하여 다른 변수들과 함께 영향을 미칠 수 있습니다.
# 
# 풍향 (deg)은 바람의 방향을 나타내기 때문에 기압에 직접적인 영향을 미치지 않지만, 바람은 기압의 차이로 인해 발생하고 이동하므로 간접적으로 기압과 연관이 있습니다. 따라서 바람은 기압 변화와 관련하여 중요한 요소입니다. 하지만 바람의 방향 자체가 기압에 직접적인 영향을 주지는 않습니다.

# #### 풍속과 관련 적은 변수
# 1. 섭씨 온도 (° ⁣C) 또는 절대 온도 (K): 풍속은 주로 기압 차이에 의해 발생하고, 온도는 풍속에 직접적인 영향을 주지 않습니다. 따라서 풍속을 결정하는데 온도는 크게 영향을 미치지 않습니다.
# 
# 2. 증기압 부족량 (mbar): 증기압 부족량은 공기 중의 수증기가 포화 상태에 미치지 못하는 압력 차이를 나타냅니다. 풍속을 결정하는데 이 변수는 직접적인 영향을 주지 않습니다.
# 
# 3. 공기 밀도 (g/m**3): 풍속은 대기의 움직임과 기압 차이에 의해 발생하고, 공기의 밀도는 풍속에 큰 영향을 미치지 않습니다.
# 
# 4. 측정 시간대: 풍속은 대기의 움직임과 기압 차이에 따라 형성되므로, 측정 시간대는 풍속에 직접적인 영향을 주지 않습니다.
# 
# 위 변수들은 풍속에 직접적인 영향을 미치지 않는다는 것을 기억하는 것이 중요합니다. 풍속은 주로 기압 차이에 의해 발생하며, 바람은 기압이 낮은 곳에서 기압이 높은 곳으로 흐르는 경향을 가지기 때문에 기압과 관련된 변수들이 풍속에 큰 영향을 미칩니다. 따라서 풍속을 예측하거나 이해하기 위해서는 주로 기압과 관련된 변수들이 중요하며, 기온, 증기압, 공기 밀도, 측정 시간대는 풍속과의 직접적인 관련성이 낮습니다.

# # Data Split

# 모델을 학습하기 위해서는 독립변수(X)와 종속변수(y)로 나누어야 합니다.  
# 독립변수는 Feature, 종속변수는 Target이라고 불리기도 합니다.  
# 
# Feature : 예측을 위해 활용되는 데이터  
# Target : 예측하고자 하는 대상  
# (이 때, 분석에 활용하지 않는 데이터인 'id'는 제거하겠습니다.)

# In[8]:


# train_x는 독립변수이므로 종속변수(풍속 (m/s))를 제거합니다. 
# 또한 target 이외의 분석에 활용하지 않는 데이터(id)를 제거합니다.
train_x = train.drop(columns=['ID', '풍속 (m/s)'], axis = 1)

# train_y는 종속변수로 값을 설정합니다.
train_y = train['풍속 (m/s)']


# In[9]:


# train에서와 마찬가지로 분석에 활용하지 않는 데이터(id)를 제거합니다.
test_x = test.drop(columns=['ID'])


# # Data Pre-processing

# ### LabelEncoder

# In[10]:


from sklearn.preprocessing import LabelEncoder

le = LabelEncoder()
le = le.fit(train_x['측정 시간대'])
train_x['측정 시간대'] = le.transform(train_x['측정 시간대'])
test_x['측정 시간대'] = le.transform(test_x['측정 시간대'])


# # 단계별 선택법

# In[11]:


df = train_x.copy()


# In[12]:


## 전진 단계별 선택법
variables = df.columns.tolist() ## 설명 변수 리스트
 
y = train_y ## 반응 변수
selected_variables = [] ## 선택된 변수들
sl_enter = 0.05
sl_remove = 0.05

sv_per_step = [] ## 각 스텝별로 선택된 변수들
adjusted_r_squared = [] ## 각 스텝별 수정된 결정계수
steps = [] ## 스텝
step = 0
while len(variables) > 0:
    remainder = list(set(variables) - set(selected_variables))
    pval = pd.Series(index=remainder) ## 변수의 p-value
    ## 기존에 포함된 변수와 새로운 변수 하나씩 돌아가면서 
    ## 선형 모형을 적합한다.
    for col in remainder: 
        X = df[selected_variables+[col]]
        X = sm.add_constant(X)
        model = sm.OLS(y,X).fit()
        pval[col] = model.pvalues[col]
 
    min_pval = pval.min()
    if min_pval < sl_enter: ## 최소 p-value 값이 기준 값보다 작으면 포함
        selected_variables.append(pval.idxmin())
        ## 선택된 변수들에대해서
        ## 어떤 변수를 제거할지 고른다.
        while len(selected_variables) > 0:
            selected_X = df[selected_variables]
            selected_X = sm.add_constant(selected_X)
            selected_pval = sm.OLS(y,selected_X).fit().pvalues[1:] ## 절편항의 p-value는 뺀다
            max_pval = selected_pval.max()
            if max_pval >= sl_remove: ## 최대 p-value값이 기준값보다 크거나 같으면 제외
                remove_variable = selected_pval.idxmax()
                selected_variables.remove(remove_variable)
            else:
                break
        
        step += 1
        steps.append(step)
        adj_r_squared = sm.OLS(y,sm.add_constant(df[selected_variables])).fit().rsquared_adj
        adjusted_r_squared.append(adj_r_squared)
        sv_per_step.append(selected_variables.copy())
    else:
        break


# In[13]:


# 원래 컬럼
variables


# In[14]:


# 단계별 선택법으로 선택된 컬럼
selected_variables


# In[15]:


# 선택되지 않은 변수
list(set(variables)-set(selected_variables))


# In[16]:


fig = plt.figure(figsize=(10,10))
fig.set_facecolor('white')
 
font_size = 15
plt.xticks(steps,[f'step {s}\n'+'\n'.join(sv_per_step[i]) for i,s in enumerate(steps)], fontsize=12)
plt.plot(steps,adjusted_r_squared, marker='o')
    
plt.ylabel('Adjusted R Squared',fontsize=font_size)
plt.grid(True)
plt.show()


# ## Column Drop

# 중요도순
# 
# 풍향, 대기압, 일, 섭씨 온도, 이슬점 온도, 상대 습도, 월, 공기 밀도, 증기압 부족량, 측정 시간대, 절대 온도, 수증기 함량, 실제 증기압, 포화 증기압

# In[17]:


# train_x = train_x.drop(['실제 증기압(mbar)', '포화 증기압(mbar)', '이슬점 온도(°C)', '절대 온도(K)'], axis=1)
train_x = train_x.drop(['증기압 부족량(mbar)', '측정 시간대'], axis=1)


# In[18]:


# test_x = test_x.drop(['실제 증기압(mbar)', '포화 증기압(mbar)', '이슬점 온도(°C)', '절대 온도(K)'], axis=1)
test_x = test_x.drop(['증기압 부족량(mbar)', '측정 시간대'], axis=1)


# In[19]:


len(train_x.columns), len(test_x.columns)


# # Explore Datset

# ## 결측치 분석

# In[20]:


train_x.isnull().sum()


# ## 이상치 분석

# In[21]:


train_x


# 모두 측정값의 숫자 분포라 특별히 이상치라고 할 만한 게 없다. 말이 안되게 튀는 숫자는 없는 것 같다.

# ## 날짜 컬럼

# 풍속이 월에 따라 바뀔 수는 있어도 일에 따라 바뀌긴 어려울 것 같다. 계절을 주기로 돌아갈 순 있어도 한 달 안에서 날짜 주기로 돌진 않을테니까.

# ** 일 컬럼 drop 시 성능이 낮아지는 것을 확인함

# ** 계절로 1, 2, 3, 4 인코딩 시 성능이 낮아지는 것을 확인함

# ## 풍향

# 풍향은 degree, 각도이기 때문에 값의 크고 작음이 의미가 없다. 사인 함수를 활용해서 변환해본다.

# In[22]:


train_x


# In[23]:


# 각도를 라디안으로 변환하여 사인 함수 적용 후 새로운 컬럼에 저장
train_x['풍향 (sin)'] = np.sin(np.deg2rad(train_x['풍향 (deg)']))
train_x


# In[24]:


test_x['풍향 (sin)'] = np.sin(np.deg2rad(test_x['풍향 (deg)']))
test_x


# ## Column Drop

# In[25]:


train_x = train_x.drop('풍향 (deg)', axis=1)
test_x = test_x.drop('풍향 (deg)', axis=1)


# ## 상관관계 분석 및 다중공선성 분석

# In[26]:


# 상관관계 분석
correlation_matrix = train_x.corr()

plt.figure(figsize=(8, 6))

sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f', linewidths=0.5)
plt.title('Correlation Heatmap', fontsize=16)

plt.show()


# In[27]:


# VIF 계산
vif = pd.DataFrame()
vif["Variables"] = train_x.columns
vif["VIF"] = [variance_inflation_factor(train_x.values, i) for i in range(train_x.shape[1])]

vif

# 일반적으로 VIF 값이 5 이상인 변수는 다중공선성의 문제가 있을 수 있으며, VIF 값이 10 이상이면 심각한 다중공선성이 발생하는 것으로 간주


# In[28]:


# 상관관계가 1인 변수들의 분산 확인
# 두 개 중 분산이 작은 것을 삭제
train_x.describe()


# ## Column Drop

# 절대 온도와 섭씨 온도는 상관관계가 1이므로 절대 온도 Drop
# 
# 수증기 함량과 실제 증기압 중 분산이 작은 수증기 함량 삭제
# 
# 다중공선성이 높은 이슬점 온도 Drop

# In[29]:


train_x = train_x.drop(['절대 온도(K)', '이슬점 온도(°C)', '공기 밀도 (g/m**3)', '포화 증기압(mbar)', '수증기 함량 (g/kg)'], axis=1)
test_x = test_x.drop(['절대 온도(K)', '이슬점 온도(°C)', '공기 밀도 (g/m**3)', '포화 증기압(mbar)', '수증기 함량 (g/kg)'], axis=1)


# In[30]:


# 컬럼 삭제 후 상관관계 재확인
correlation_matrix = train_x.corr()

plt.figure(figsize=(8, 6))

sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f', linewidths=0.5)
plt.title('Correlation Heatmap', fontsize=16)

plt.show()


# # 분포 시각화

# In[31]:


sns.pairplot(train_x)
plt.title("Pair Plot")
plt.show()


# # Data Split

# In[32]:


# 학습 데이터와 테스트 데이터로 분리
X_train, X_test, y_train, y_test = train_test_split(train_x, train_y, test_size=0.2, shuffle=True, random_state=42)


# # Scaling

# In[33]:


# 학습 데이터로 StandardScaler 객체 생성 및 적용
scaler = StandardScaler()
scaler.fit(X_train)

# 학습 데이터와 테스트 데이터를 각각 표준화
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)


# In[34]:


# # 'col1', 'col2', 'col3' 컬럼만 선택하여 스케일링할 경우
# selected_cols = ['col1', 'col2', 'col3']

# # StandardScaler 객체 생성
# scaler = StandardScaler()

# # 'col1', 'col2', 'col3' 컬럼만 스케일링하기 위해 해당 컬럼들을 추출하여 스케일링
# X_train_scaled_selected = X_train[selected_cols].copy()
# X_train_scaled_selected = scaler.fit_transform(X_train_scaled_selected)

# # 스케일링된 컬럼들을 다시 원래 X_train에 삽입
# X_train[selected_cols] = X_train_scaled_selected

# # 나머지 컬럼들은 스케일링하지 않으므로 그대로 사용
# # ...

# # 위와 동일한 방식으로 X_test에 대해서도 스케일링 수행
# X_test_scaled_selected = X_test[selected_cols].copy()
# X_test_scaled_selected = scaler.transform(X_test_scaled_selected)
# X_test[selected_cols] = X_test_scaled_selected

# # 나머지 컬럼들은 스케일링하지 않으므로 그대로 사용
# # ...


# # Hyperopt

# In[35]:


from hyperopt import fmin, tpe, hp, Trials

# 데이터를 준비하고 X_train, y_train, X_test, y_test로 나누는 등의 전처리 작업을 수행합니다.


# In[36]:


# Objective 함수 정의 (하이퍼파라미터 튜닝 대상)
def objective(params):
    model = XGBRegressor(
        n_estimators=int(params['n_estimators']),
        max_depth=int(params['max_depth']),
        learning_rate=params['learning_rate'],
        reg_lambda=params['reg_lambda'],
        colsample_bytree=params['colsample_bytree'],
        max_leaves=int(params['max_leaves']),
        subsample=params['subsample'],
        gamma=params['gamma'],
        
        random_state=42,
        n_jobs=-1,
        # tree_method='gpu_hist', gpu_id=0
    )
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    mae = mean_absolute_error(y_test, y_pred)
    return mae


# In[37]:


# 하이퍼파라미터 공간 정의
space = {
    'n_estimators': hp.quniform('n_estimators', 150, 250, 1),
    'max_depth': hp.quniform('max_depth', 10, 20, 1),
    'learning_rate': hp.loguniform('learning_rate', np.log(0.001), np.log(0.1)),
    'reg_lambda': hp.loguniform('reg_lambda', 0.01, 5.0),
    'colsample_bytree': hp.uniform('colsample_bytree', 0.1, 1.0),
    'max_leaves': hp.quniform('max_leaves', 20, 40, 1),
    'subsample': hp.uniform('subsample', 0.1, 1.0),
    'gamma': hp.uniform('gamma', 0.0, 0.5)
}


# In[38]:


# Hyperopt를 사용하여 하이퍼파라미터 튜닝
trials = Trials()
best = fmin(fn=objective,
            space=space,
            algo=tpe.suggest,
            max_evals=200,  # 시도할 횟수
            trials=trials,
            verbose=1)


# In[39]:


# 최적의 하이퍼파라미터 출력
print("Best Hyperparameters:")
print(best)


# # Model Definition

# In[40]:


# 최적의 하이퍼파라미터로 모델 생성
best_model = XGBRegressor(
    n_estimators=int(best['n_estimators']),
    max_depth=int(best['max_depth']),
    learning_rate=best['learning_rate'],
    reg_lambda=best['reg_lambda'],
    colsample_bytree=best['colsample_bytree'],
    max_leaves=int(best['max_leaves']),
    subsample=best['subsample'],
    gamma=best['gamma'],
    
    random_state=42,
    n_jobs=-1,
    # tree_method='gpu_hist', gpu_id=0
)


# # Model Fit

# In[41]:


# 최적의 하이퍼파라미터로 모델 훈련
best_model.fit(X_train, y_train)


# In[42]:


# 예측 결과 출력
y_pred = best_model.predict(X_test)
mae = mean_absolute_error(y_test, y_pred)
print("MAE: {:.4f}".format(mae))
# 0.4252 {'colsample_bytree': 0.9, 'learning_rate': 0.03, 'max_depth': 19, 'max_leaves': 29, 'n_estimators': 195, 'subsample': 0.7}


# In[43]:


xgboost.plot_importance(best_model)


# # Prediction

# In[44]:


# dt_pred = dt_clf.predict(test_x) 


# In[45]:


# 테스트 데이터로 예측
xgb_pred = best_model.predict(test_x)


# # Submit

# In[46]:


submission = pd.read_csv('./data/sample_submission.csv')


# In[47]:


submission['풍속 (m/s)'] = xgb_pred
submission.head()


# In[48]:


submission.to_csv('submission/submission_'+datetime.datetime.now().strftime("%Y%m%d%H%M")+'.csv', index=  False)


# //끝//