import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from scipy.stats import norm, skew

# Set visualization style
sns.set_style('whitegrid')
plt.rcParams['figure.figsize'] = (12, 8)

from sklearn.linear_model import (
    ElasticNet, Lasso, BayesianRidge, LassoLarsIC, Ridge
)
from sklearn.ensemble import (
    RandomForestRegressor, GradientBoostingRegressor, StackingRegressor
)
from sklearn.kernel_ridge import KernelRidge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import RobustScaler
from sklearn.model_selection import (
    KFold, cross_val_score, train_test_split, GridSearchCV
)
from sklearn.metrics import mean_squared_error
import xgboost as xgb
import lightgbm as lgb

# Load training and test datasets
train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')

# Display basic information
print(f'Training set shape: {train.shape}')
print(f'Test set shape: {test.shape}')
print(f'\nFirst few rows of training data:')
print(train.head())

# Descriptive statistics
print("Sale Price Statistics:")
print(train['SalePrice'].describe())

# Visualize distribution
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Histogram with KDE
sns.histplot(train['SalePrice'], kde=True, ax=axes[0])
axes[0].set_title('Sale Price Distribution')
axes[0].set_xlabel('Sale Price ($)')

# Box plot to identify outliers
sns.boxplot(y=train['SalePrice'], ax=axes[1])
axes[1].set_title('Sale Price Box Plot')
axes[1].set_ylabel('Sale Price ($)')

plt.tight_layout()
plt.show()

# Separate numerical and categorical features
numerical_features = train.select_dtypes(
    include=['int64', 'float64']
).columns.tolist()

categorical_features = train.select_dtypes(
    include=['object']
).columns.tolist()

# Remove Id and SalePrice from numerical
numerical_features.remove('Id')
if 'SalePrice' in numerical_features:
    numerical_features.remove('SalePrice')

print(f'Numerical features: {len(numerical_features)}')
print(f'Categorical features: {len(categorical_features)}')

# Calculate correlations with SalePrice
correlations = train[numerical_features + ['SalePrice']].corr()
saleprice_corr = correlations['SalePrice'].sort_values(
    ascending=False
)

# Display top correlations
print("Top 10 features correlated with SalePrice:")
print(saleprice_corr[1:11])

# Visualize correlation heatmap
top_features = saleprice_corr[1:11].index.tolist() + ['SalePrice']
plt.figure(figsize=(12, 10))
sns.heatmap(
    train[top_features].corr(),
    annot=True,
    fmt='.2f',
    cmap='coolwarm',
    center=0
)
plt.title('Correlation Matrix of Top Features')
plt.show()

# Store IDs for later use in submissions
train_id = train['Id']
test_id = test['Id']

# Drop Id from datasets
train = train.drop('Id', axis=1)
test = test.drop('Id', axis=1)

print(f'Training set shape after removing Id: {train.shape}')
print(f'Test set shape after removing Id: {test.shape}')

# Get numerical columns and count unique values
numeric_cols = train.select_dtypes(include=[np.number]).columns
unique_counts = {}

for col in numeric_cols:
    if col != 'SalePrice':
        unique_counts[col] = train[col].nunique()

# Sort by number of unique values (descending)
sorted_features = sorted(
    unique_counts.items(), 
    key=lambda x: x[1], 
    reverse=True
)

# Select top features with high variance
high_variance_features = [f[0] for f in sorted_features[:19]]

# Check correlation with SalePrice
correlations = train[high_variance_features + ['SalePrice']].corr()['SalePrice']
print(correlations.sort_values(ascending=False))

def check_outliers(feature_name):
    """
    Visualize potential outliers for a feature:
    1. Test set scatter (are similar values present?)
    2. Train/Test distribution comparison
    3. Train set scatter against SalePrice
    """
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    
    # Plot 1: Test set scatter
    axes[0].scatter(
        range(len(test)), 
        test[feature_name],
        alpha=0.5
    )
    axes[0].set_title(f'{feature_name} in Test Set')
    axes[0].set_xlabel('Sample Index')
    axes[0].set_ylabel(feature_name)
    
    # Plot 2: Distribution comparison
    axes[1].hist(train[feature_name], bins=50, 
                alpha=0.5, label='Train')
    axes[1].hist(test[feature_name], bins=50, 
                alpha=0.5, label='Test')
    axes[1].set_title('Distribution Comparison')
    axes[1].legend()
    
    # Plot 3: Train set scatter vs SalePrice
    axes[2].scatter(
        train[feature_name], 
        train['SalePrice'],
        alpha=0.5
    )
    axes[2].set_title(f'{feature_name} vs SalePrice (Train)')
    axes[2].set_xlabel(feature_name)
    axes[2].set_ylabel('SalePrice')
    
    plt.tight_layout()
    plt.show()

# Visualize GrLivArea
check_outliers('GrLivArea')

# Identify outliers: large GrLivArea but low SalePrice
outliers_grlivarea = train[
    (train['GrLivArea'] > 4000) & 
    (train['SalePrice'] < 300000)
]
print(f'Found {len(outliers_grlivarea)} outliers in GrLivArea')
print(outliers_grlivarea[['GrLivArea', 'SalePrice']])

# Remove outliers
train = train.drop(outliers_grlivarea.index)
print(f'\nNew training set size: {train.shape[0]}')

import numpy as np
from scipy.stats import norm, skew

# Calculate skewness before transformation
original_skew = skew(train['SalePrice'])
print(f'Original skewness: {original_skew:.4f}')

# Apply log transformation: log(1+x) to handle any zeros
train['SalePrice'] = np.log1p(train['SalePrice'])

# Calculate skewness after transformation
transformed_skew = skew(train['SalePrice'])
print(f'Transformed skewness: {transformed_skew:.4f}')

# Visualize the transformation
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Original distribution (we'll recreate from backup)
sns.histplot(np.expm1(train['SalePrice']), kde=True, ax=axes[0])
axes[0].set_title('Original Sale Price Distribution')
axes[0].set_xlabel('Sale Price ($)')

# Transformed distribution
sns.histplot(train['SalePrice'], kde=True, ax=axes[1])
axes[1].set_title('Log-Transformed Sale Price')
axes[1].set_xlabel('log(Sale Price)')

plt.tight_layout()
plt.show()

# Store the target variable separately
y_train = train['SalePrice'].copy()
train = train.drop('SalePrice', axis=1)

# Concatenate train and test
n_train = train.shape[0]
n_test = test.shape[0]
all_data = pd.concat([train, test], axis=0, sort=False)

print(f'Training samples: {n_train}')
print(f'Test samples: {n_test}')
print(f'Combined dataset shape: {all_data.shape}')
print(f'\nNumber of features: {all_data.shape[1]}')

# Calculate missing value percentages
missing = all_data.isnull().sum()
missing_percent = (100 * missing / len(all_data)).sort_values(
    ascending=False
)

# Filter to show only features with missing values
missing_data = pd.DataFrame({
    'Missing Count': missing,
    'Percentage': missing_percent
})
missing_data = missing_data[missing_data['Missing Count'] > 0]

print(f'Features with missing values: {len(missing_data)}')
print('\nTop 10 features with most missing data:')
print(missing_data.head(10))

# Separate low percentage missing value columns (LPMVC)
low_missing_cols = missing_data[
    missing_data['Percentage'] < 5
].index.tolist()

print(f'Features with <5% missing: {len(low_missing_cols)}')

# Distinguish numeric vs categorical
numeric_lpmvc = all_data[low_missing_cols].select_dtypes(
    include=[np.number]
).columns.tolist()

categorical_lpmvc = all_data[low_missing_cols].select_dtypes(
    include=['object']
).columns.tolist()

# Fill numeric features with 0 (structural zeros)
for col in numeric_lpmvc:
    all_data[col].fillna(0, inplace=True)
    
# Fill categorical features with mode (most frequent value)
for col in categorical_lpmvc:
    if col != 'Utilities':  # Handle Utilities separately
        mode_value = all_data[col].mode()[0]
        all_data[col].fillna(mode_value, inplace=True)
        print(f'{col}: Filled with "{mode_value}"')

# High percentage missing value columns (HPMVC)
high_missing_cols = missing_data[
    missing_data['Percentage'] >= 5
].index.tolist()

# These are mostly "absence of feature" type missing values
# Fill numeric with 0, categorical with 'None'

for col in high_missing_cols:
    if all_data[col].dtype == 'object':
        all_data[col].fillna('None', inplace=True)
    else:
        all_data[col].fillna(0, inplace=True)

# Fill LotFrontage based on median of neighborhood
all_data['LotFrontage'] = all_data.groupby('Neighborhood')[
    'LotFrontage'
].transform(
    lambda x: x.fillna(x.median())
)

print('LotFrontage missing values after imputation:', 
      all_data['LotFrontage'].isnull().sum())

# Check for any remaining missing values
remaining_missing = all_data.isnull().sum().sum()
print(f'Total missing values remaining: {remaining_missing}')

if remaining_missing == 0:
    print('✅ Successfully handled all missing values!')
else:
    print(f'⚠️ Still have {remaining_missing} missing values to address')

# Total square footage = basement + first floor + second floor
all_data['TotalSF'] = (
    all_data['TotalBsmtSF'] + 
    all_data['1stFlrSF'] + 
    all_data['2ndFlrSF']
)

# Calculate correlation with target (on training data)
correlation = train_with_target[
    ['TotalSF', 'SalePrice']
].corr().iloc[0,1]

print(f'TotalSF correlation with SalePrice: {correlation:.4f}')

# Total bathrooms = full baths + (0.5 × half baths)
all_data['TotalBathrooms'] = (
    all_data['FullBath'] + 
    (0.5 * all_data['HalfBath']) +
    all_data['BsmtFullBath'] + 
    (0.5 * all_data['BsmtHalfBath'])
)

print('TotalBathrooms statistics:')
print(all_data['TotalBathrooms'].describe())

# Basement Quality - Ordinal Encoding
# Higher number = better quality
basement_quality_map = {
    'None': 0,
    'Po': 1,    # Poor
    'Fa': 2,    # Fair
    'TA': 3,    # Average/Typical
    'Gd': 4,    # Good
    'Ex': 5     # Excellent
}

all_data['BsmtQual'] = all_data['BsmtQual'].map(basement_quality_map)
all_data['BsmtCond'] = all_data['BsmtCond'].map(basement_quality_map)

# Basement Exposure - Ordinal with different meaning
basement_exposure_map = {
    'None': 0,
    'No': 1,   # No exposure
    'Mn': 2,   # Minimum exposure
    'Av': 3,   # Average exposure
    'Gd': 4    # Good exposure
}

all_data['BsmtExposure'] = all_data['BsmtExposure'].map(
    basement_exposure_map
)

# Standard quality mapping used across multiple features
quality_map = {
    'None': 0, 'Po': 1, 'Fa': 2, 
    'TA': 3, 'Gd': 4, 'Ex': 5
}

# Apply to multiple quality-related features
quality_features = [
    'KitchenQual', 'FireplaceQu', 'ExterQual', 
    'ExterCond', 'HeatingQC', 'GarageQual', 
    'GarageCond', 'PoolQC'
]

for feature in quality_features:
    all_data[feature] = all_data[feature].map(quality_map)
    print(f'Encoded {feature}')

# Garage Finish - Ordinal encoding
garage_finish_map = {
    'None': 0,
    'Unf': 1,    # Unfinished
    'RFn': 2,    # Rough Finished
    'Fin': 3     # Finished
}

all_data['GarageFinish'] = all_data['GarageFinish'].map(
    garage_finish_map
)

# GarageType - Keep as categorical for one-hot encoding later
# Types: Attchd, Detchd, BuiltIn, Basment, CarPort, 2Types, None

# Create interaction: Garage Score = Cars × Quality
all_data['GarageScore'] = (
    all_data['GarageCars'] * all_data['GarageQual']
)

print('Created GarageScore interaction feature')

# House Age at time of sale
all_data['HouseAge'] = all_data['YrSold'] - all_data['YearBuilt']

# Years since remodeling
all_data['YearsSinceRemodel'] = (
    all_data['YrSold'] - all_data['YearRemodAdd']
)

# Was the house remodeled? (Binary feature)
all_data['WasRemodeled'] = (
    all_data['YearRemodAdd'] != all_data['YearBuilt']
).astype(int)

# Garage Age
all_data['GarageAge'] = all_data['YrSold'] - all_data['GarageYrBlt']
all_data['GarageAge'].fillna(0, inplace=True)  # No garage = 0

print('Created age-related features')
print(f'Average house age: {all_data["HouseAge"].mean():.1f} years')

# Analyze neighborhood impact on price
neighborhood_stats = train.groupby('Neighborhood')['SalePrice'].agg(
    ['mean', 'median', 'count']
).sort_values('mean', ascending=False)

print('Top 5 most expensive neighborhoods:')
print(neighborhood_stats.head())

print('\nBottom 5 least expensive neighborhoods:')
print(neighborhood_stats.tail())

# Neighborhood will be one-hot encoded later
# But we can create a "NeighborhoodScore" based on average price
neighborhood_scores = train.groupby('Neighborhood')[
    'SalePrice'
].mean().to_dict()

all_data['NeighborhoodScore'] = all_data['Neighborhood'].map(
    neighborhood_scores
)

from sklearn.preprocessing import PolynomialFeatures

# Select key features for polynomial expansion
poly_features = ['OverallQual', 'GrLivArea', 'TotalSF']

# Create squared terms (degree=2)
for feature in poly_features:
    all_data[f'{feature}_Squared'] = all_data[feature] ** 2
    print(f'Created {feature}_Squared')

# Create interaction term: Quality × Size
all_data['QualitySize'] = (
    all_data['OverallQual'] * all_data['GrLivArea']
)

print('\nCreated polynomial and interaction features')

# Has swimming pool?
all_data['HasPool'] = (all_data['PoolArea'] > 0).astype(int)

# Has second floor?
all_data['Has2ndFloor'] = (all_data['2ndFlrSF'] > 0).astype(int)

# Has basement?
all_data['HasBasement'] = (all_data['TotalBsmtSF'] > 0).astype(int)

# Has garage?
all_data['HasGarage'] = (all_data['GarageArea'] > 0).astype(int)

# Has fireplace?
all_data['HasFireplace'] = (all_data['Fireplaces'] > 0).astype(int)

# Has porch?
all_data['HasPorch'] = (
    all_data['OpenPorchSF'] + 
    all_data['EnclosedPorch'] + 
    all_data['3SsnPorch'] + 
    all_data['ScreenPorch']
) > 0

all_data['HasPorch'] = all_data['HasPorch'].astype(int)

print('Created binary indicator features')

from scipy.stats import skew

# Get all numerical features
numeric_features = all_data.dtypes[
    all_data.dtypes != 'object'
].index

# Calculate skewness for each feature
skewness = all_data[numeric_features].apply(
    lambda x: skew(x.dropna())
)

# Filter features with high skewness (|skew| > 0.75)
skewed_features = skewness[abs(skewness) > 0.75].sort_values(
    ascending=False
)

print(f'Found {len(skewed_features)} skewed features')
print('\nTop 10 most skewed features:')
print(skewed_features.head(10))

from scipy.special import boxcox1p
from scipy.stats import boxcox_normmax

# Box-Cox parameter lambda=0.15 works well for most features
lam = 0.15

for feature in skewed_features.index:
    # boxcox1p applies Box-Cox transformation to (1 + x)
    # This handles zero values gracefully
    all_data[feature] = boxcox1p(all_data[feature], lam)

# Verify skewness reduction
new_skewness = all_data[skewed_features.index].apply(
    lambda x: skew(x.dropna())
)

print('\nSkewness before and after transformation:')
comparison = pd.DataFrame({
    'Before': skewed_features,
    'After': new_skewness
})
print(comparison.head(10))

# One-hot encode categorical variables
all_data = pd.get_dummies(all_data)

# Separate back into train and test
train = all_data[:len(y_train)]
test = all_data[len(y_train):]

print(f'Final training set shape: {train.shape}')
print(f'Final test set shape: {test.shape}')
print(f'Total features after encoding: {train.shape[1]}')

import tensorflow as tf

def build_dnn_model(input_dim, layers=[256, 128, 64], 
                       dropout_rate=0.3, l2_reg=0.01):
    """Build and compile DNN model"""
    
    model = tf.keras.Sequential()
    
    # Input layer
    model.add(tf.keras.layers.Dense(
        layers[0], 
        activation='relu',
        input_dim=input_dim,
        kernel_regularizer=tf.keras.regularizers.l2(l2_reg)
    ))
    model.add(tf.keras.layers.Dropout(dropout_rate))
    
    # Hidden layers
    for neurons in layers[1:]:
        model.add(tf.keras.layers.Dense(
            neurons,
            activation='relu',
            kernel_regularizer=tf.keras.regularizers.l2(l2_reg)
        ))
        model.add(tf.keras.layers.Dropout(dropout_rate))
    
    # Output layer
    model.add(tf.keras.layers.Dense(1, activation='linear'))
    
    # Compile model
    model.compile(
        optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),
        loss='mse',
        metrics=['mae']
    )
    
    return model

# Build and train
dnn_model = build_dnn_model(input_dim=train.shape[1])

# Early stopping to prevent overfitting
early_stop = tf.keras.callbacks.EarlyStopping(
    monitor='val_loss',
    patience=50,
    restore_best_weights=True
)

# Train model
history = dnn_model.fit(
    train, y_train,
    epochs=500,
    batch_size=32,
    validation_split=0.2,
    callbacks=[early_stop],
    verbose=0
)

# Make predictions
dnn_predictions = dnn_model.predict(test)
dnn_predictions = np.expm1(dnn_predictions.flatten())

print(f'DNN Final Validation Loss: {history.history["val_loss"][-1]:.5f}')