import pandas as pd
import matplotlib.pyplot as plt

import tensorflow as tf
from tensorflow.keras import layers
import seaborn as sns

from sklearn.datasets import fetch_california_housing

data_train = pd.read_csv('/kaggle/input/playground-series-s3e1/train.csv', index_col='id')
data_test = pd.read_csv('/kaggle/input/playground-series-s3e1/test.csv', index_col='id')
features = data_train.columns[:-1]

original_df, original_y = fetch_california_housing(return_X_y=True)
original_df = pd.DataFrame(original_df, columns=features)
original_df['MedHouseVal'] = original_y

def add_secondary_plot(
    df, column, target_column, ax, n_bins, color=3,
    show_yticks=False,
):
    secondary_ax = ax.twinx()
    bins = pd.cut(df[column], bins=n_bins)
    bins = pd.IntervalIndex(bins)
    bins = (bins.left + bins.right) / 2
    target = df.groupby(bins)[target_column].mean()
    target.plot(
        ax=secondary_ax, linestyle='',
        marker='.', color=colors[color], label=f'Mean {target_column}'
    )
    secondary_ax.grid(visible=False)
    
    if not show_yticks:
        secondary_ax.get_yaxis().set_ticks([])
        
    return secondary_ax

import matplotlib.colors as mpl_colors

def hex_to_rgb(h):
    h = h.lstrip('#')
    return tuple(int(h[i:i+2], 16)/255 for i in (0, 2, 4))

palette = ['#b4d2b1', '#568f8b', '#1d4a60', '#cd7e59', '#ddb247', '#d15252']
palette_rgb = [hex_to_rgb(x) for x in palette]
cmap = mpl_colors.ListedColormap(palette_rgb)
colors = cmap.colors
bg_color= '#fdfcf6'

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
from sklearn.model_selection import train_test_split
from sklearn.datasets import fetch_california_housing

train_df = pd.read_csv('/kaggle/input/playground-series-s3e1/train.csv')
test_df = pd.read_csv('/kaggle/input/playground-series-s3e1/test.csv')
submission = pd.read_csv('/kaggle/input/playground-series-s3e1/sample_submission.csv')
train_df = train_df.drop('id', axis=1)

train_df['r'] = np.sqrt(train_df['Latitude']**2 + train_df['Longitude']**2)
train_df['theta'] = np.arctan2(train_df['Latitude'], train_df['Longitude'])

test_df['r'] = np.sqrt(test_df['Latitude']**2 + test_df['Longitude']**2)
test_df['theta'] = np.arctan2(test_df['Latitude'], test_df['Longitude'])

df = pd.concat([train_df, test_df], axis=0, ignore_index=True)

emb_size = 20
precision = 1e6 

latlon = np.expand_dims(df[['Latitude', 'Longitude']].values, axis=-1) 

m = np.exp(np.log(precision) / emb_size) 
angle_freq = m ** np.arange(emb_size) 
angle_freq = angle_freq.reshape(1, 1, emb_size) 

latlon = latlon * angle_freq 
latlon[..., 0::2] = np.cos(latlon[..., 0::2]) 
latlon[..., 1::2] = np.sin(latlon[..., 1::2]) 
latlon = latlon.reshape(-1, 2 * emb_size) 

df['exp_latlon1'] = [lat[0] for lat in latlon]
df['exp_latlon2'] = [lat[1] for lat in latlon]

from sklearn.decomposition import PCA

def pca(data):
    '''
    input: dataframe containing Latitude(x) and Longitude(y)
    '''
    coordinates = data[['Latitude','Latitude']].values
    pca_obj = PCA().fit(coordinates)
    pca_x = pca_obj.transform(data[['Latitude', 'Longitude']].values)[:,0]
    pca_y = pca_obj.transform(data[['Latitude', 'Longitude']].values)[:,1]
    return pca_x, pca_y

# train_df['pca_x'], train_df['pca_y'] = pca(train_df)
# test_df['pca_x'], test_df['pca_y'] = pca(test_df)
df['pca_x'], df['pca_y'] = pca(df)

from umap import UMAP

# UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction
coordinates = df[['Latitude', 'Longitude']].values
umap = UMAP(n_components=2, n_neighbors=50, random_state=228).fit(coordinates)
df['umap_lat'] = umap.transform(coordinates)[:,0]
df['umap_lon'] = umap.transform(coordinates)[:,1]

def crt_crds(df):
    df['rot_15_x'] = (np.cos(np.radians(15)) * df['Longitude']) + \
                      (np.sin(np.radians(15)) * df['Latitude'])
    
    df['rot_15_y'] = (np.cos(np.radians(15)) * df['Latitude']) + \
                      (np.sin(np.radians(15)) * df['Longitude'])
    
    df['rot_30_x'] = (np.cos(np.radians(30)) * df['Longitude']) + \
                      (np.sin(np.radians(30)) * df['Latitude'])
    
    df['rot_30_y'] = (np.cos(np.radians(30)) * df['Latitude']) + \
                      (np.sin(np.radians(30)) * df['Longitude'])
    
    df['rot_45_x'] = (np.cos(np.radians(45)) * df['Longitude']) + \
                      (np.sin(np.radians(45)) * df['Latitude'])
    return df

# train_df = crt_crds(train_df)
# test_df = crt_crds(test_df)
df = crt_crds(df)

from haversine import haversine

Sac = (38.576931, -121.494949)
SF = (37.780080, -122.420160)
SJ = (37.334789, -121.888138)
LA = (34.052235, -118.243683)
SD = (32.715759, -117.163818)

df['dist_Sac'] = df.apply(lambda x: haversine((x['Latitude'], x['Longitude']), Sac, unit='ft'), axis=1)
df['dist_SF'] = df.apply(lambda x: haversine((x['Latitude'], x['Longitude']), SF, unit='ft'), axis=1)
df['dist_SJ'] = df.apply(lambda x: haversine((x['Latitude'], x['Longitude']), SJ, unit='ft'), axis=1)
df['dist_LA'] = df.apply(lambda x: haversine((x['Latitude'], x['Longitude']), LA, unit='ft'), axis=1)
df['dist_SD'] = df.apply(lambda x: haversine((x['Latitude'], x['Longitude']), SD, unit='ft'), axis=1)
df['dist_nearest_city'] = df[['dist_Sac', 'dist_SF', 'dist_SJ', 
                              'dist_LA', 'dist_SD']].min(axis=1)

from shapely.geometry import LineString, Point

coast_points = LineString([(32.6644, -117.1613), (33.2064, -117.3831),
                           (33.7772, -118.2024), (34.4634, -120.0144),
                           (35.4273, -120.8819), (35.9284, -121.4892),
                           (36.9827, -122.0289), (37.6114, -122.4916),
                           (38.3556, -123.0603), (39.7926, -123.8217),
                           (40.7997, -124.1881), (41.7558, -124.1976)])

df['dist_to_coast'] = df.apply(lambda x: Point(x['Latitude'], x['Longitude']).distance(coast_points), axis=1)

# combine latitude and longitude
# codes from 
# https://datascience.stackexchange.com/questions/49553/combining-latitude-longitude-position-into-single-feature
from math import radians, cos, sin, asin, sqrt

def single_pt_haversine(lat, lng, degrees=True):
    """
    'Single-point' Haversine: Calculates the great circle distance
    between a point on Earth and the (0, 0) lat-long coordinate
    """
    r = 6371 # Earth's radius (km). Have r = 3956 if you want miles

    # Convert decimal degrees to radians
    if degrees:
        lat, lng = map(radians, [lat, lng])

    # 'Single-point' Haversine formula
    a = sin(lat/2)**2 + cos(lat) * sin(lng/2)**2
    d = 2 * r * asin(sqrt(a)) 

    return d
# add more metric 
# referred to this discussion
# https://www.kaggle.com/competitions/playground-series-s3e1/discussion/376210

def manhattan(lat,lng):
    return np.abs(lat) + np.abs(lng)
def euclidean(lat,lng):
    return (lat**2 + lng**2) **0.5

def add_combine(df):      
    df['haversine'] = [single_pt_haversine(x, y) for x, y in zip(df.Latitude, df.Longitude)]
    df['manhattan'] = [manhattan(x,y) for x,y in zip(df.Latitude, df.Longitude)]
    df['euclidean'] = [euclidean(x,y) for x,y in zip(df.Latitude,df.Longitude)]
    return df

df = add_combine(df)

train_df = df.iloc[:-len(test_df),:]
test_df = df.iloc[-len(test_df):,:].drop('MedHouseVal', axis=1).reset_index(drop=True)

X = train_df.drop(['MedHouseVal', 'id'], axis=1)
y = train_df.MedHouseVal
X_test = test_df.drop('id', axis=1)

import catboost
from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error

n_folds = 15

MAX_ITER = 15000
PATIENCE = 1000
DISPLAY_FREQ = 100

eval_predsCB = []
predsCB = []

k_fold = KFold(n_splits=n_folds, random_state=42, shuffle=True)

MODEL_PARAMS = {
                'random_seed': 1234,    
#                 'learning_rate': 0.1,   # 0.15: 0.5678, 0.12: 0.5685, 0.1: 0.56757, 0.05: 0.57, 0.01, 0.57             
                'iterations': MAX_ITER,
                'early_stopping_rounds': PATIENCE,
#                 'metric_period': DISPLAY_FREQ,
                'use_best_model': True,
                'eval_metric': 'RMSE',
                'verbose': 1000,
#                 'task_type': 'GPU'
               }


for train_index, test_index in k_fold.split(X, y):
    X_train, X_valid = X.iloc[train_index], X.iloc[test_index]
    y_train, y_valid = y.iloc[train_index], y.iloc[test_index]
    
    model = catboost.CatBoostRegressor(**MODEL_PARAMS)
    
    model.fit(X=X_train, y=y_train,
          eval_set=[(X_valid, y_valid)],
          early_stopping_rounds = PATIENCE,
#           metric_period = DISPLAY_FREQ
         )
    predsCB.append(model.predict(X_test))
#     eval_predsCB.append(model.predict(X))
#     print("RMSE valid = {}".format(mean_squared_error(y_valid, model.predict(X_valid))))
#     print("RMSE full = {}".format(mean_squared_error(y, model.predict(X))))

from xgboost import XGBRegressor

# n_folds = 20
k_fold = KFold(n_splits=n_folds, random_state=42, shuffle=True)

eval_predsXB = []
predsXB = []

PATIENCE = 200

MODEL_PARAMS = {       'n_estimators': 1000, #1000, 5000
#                        'learning_rate': 0.05,
                       'max_depth': 4, # 3
                       'colsample_bytree': 0.9, # 0.95
                       'subsample': 1,
                       'reg_lambda': 20,
                       'early_stopping_rounds': PATIENCE,
#                        'tree_method': 'gpu_hist',
                       'seed': 1
}

for train_index, test_index in k_fold.split(X, y):
    X_train, X_valid = X.iloc[train_index], X.iloc[test_index]
    y_train, y_valid = y.iloc[train_index], y.iloc[test_index]
    
    model = XGBRegressor(**MODEL_PARAMS)
    
    model.fit(X=X_train, y=y_train,
          eval_set=[(X_valid, y_valid)],
#           early_stopping_rounds = PATIENCE,
          verbose = 100
         )
    predsXB.append(model.predict(X_test))
#     eval_predsXB.append(model.predict(X))

import lightgbm as lgbm
from lightgbm.sklearn import LGBMRegressor

# n_folds = 20
k_fold = KFold(n_splits=n_folds, random_state=42, shuffle=True)

eval_predsLB = []
predsLB = []

MODEL_PARAMS = {
                       'learning_rate': 0.01,
                       'max_depth': 9,
                       'num_leaves': 90,
                       'colsample_bytree': 0.8,
                       'subsample': 0.9,
                       'subsample_freq': 5,
                       'min_child_samples': 36,
                       'reg_lambda': 28,
                       'n_estimators': 20000,
                       'metric': 'rmse',
                       'random_state': 1
}

callbacks = [lgbm.early_stopping(30, verbose=1), lgbm.log_evaluation(period=0)]

for train_index, test_index in k_fold.split(X, y):
    X_train, X_valid = X.iloc[train_index], X.iloc[test_index]
    y_train, y_valid = y.iloc[train_index], y.iloc[test_index]
    
    model = lgbm.LGBMRegressor(**MODEL_PARAMS)
    
    model.fit(X=X_train, y=y_train,
          eval_set=[(X_valid, y_valid)],
#           early_stopping_rounds = PATIENCE,
          callbacks=callbacks
         )
    predsLB.append(model.predict(X_test))
#     eval_predsLB.append(model.predict(X))

a = 0.4
b = 0.2
c = 0.4

predCB = np.average(np.array(predsCB),axis=0)
predXB = np.average(np.array(predsXB),axis=0)
predLB = np.average(np.array(predsLB),axis=0)
pred = predCB * a + predXB * b + predLB * c

submission['MedHouseVal'] = pred
submission

vals = train_df['MedHouseVal'].unique().tolist()
submission['MedHouseVal'] = submission['MedHouseVal'].apply(lambda x: min(vals, key=lambda v: abs(v - x)))
submission.MedHouseVal.clip(0, 5, inplace=True)

submission.to_csv('submission_v3.1.4.csv', index=False)