# Import required libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
from sklearn.model_selection import (
    train_test_split, cross_val_score, GridSearchCV, LeaveOneOut, cross_val_predict
)
from sklearn.preprocessing import (
    StandardScaler, PolynomialFeatures, RobustScaler, PowerTransformer
)
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression, Lasso, Ridge, ElasticNet
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import (
    RandomForestRegressor, GradientBoostingRegressor, VotingRegressor
)
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from statsmodels.stats.outliers_influence import variance_inflation_factor

# Suppress warnings for cleaner output
warnings.filterwarnings('ignore')

# ---------------------------
# 1. Data Loading and Preparation
# ---------------------------

# Define a function to clean 'LGA NAME'
def clean_lga_name(name):
    if isinstance(name, str):
        return name.strip().upper().replace('-', ' ').replace('&', 'AND')
    return name

# Load the initial data files
lga_population_data = pd.read_excel('/Users/rahulkulkarni/Downloads/Data/LGA_population_data.xlsx', sheet_name=None)
egm_data = pd.read_csv('/Users/rahulkulkarni/Downloads/Data/EGM.csv')

# Load the correct population data from the 'LGA Population Data' sheet
lga_population_actual_data = lga_population_data['LGA Population Data']

# Standardize the column names for consistency
lga_population_actual_data.rename(columns=lambda x: x.strip().upper(), inplace=True)

# Rename LGA column for clarity
lga_population_actual_data['LGA NAME'] = lga_population_actual_data['LGA FROM WIKIPEDIA'].str.upper()

# Extract relevant columns for population data
lga_population_actual_data = lga_population_actual_data[['LGA NAME', 'TOTAL POPULATION']]

# Extract only 2012 data from EGM data and standardize column names
egm_data.rename(columns=lambda x: x.strip().upper(), inplace=True)
egm_data_2012 = egm_data[['LGA NAME', '2012']]

# Clean 'LGA NAME' in EGM data
egm_data_2012['LGA NAME'] = egm_data_2012['LGA NAME'].apply(clean_lga_name)

# Merge EGM data with the actual population data using inner join
merged_data_2012 = pd.merge(egm_data_2012, lga_population_actual_data, on='LGA NAME', how='inner')

# Calculate per capita spending
merged_data_2012['PER CAPITA EGM SPENDING'] = merged_data_2012['2012'] / merged_data_2012['TOTAL POPULATION']

# Now load the latest uploaded file containing LGA percentage data
latest_lga_population_data = pd.read_excel('/Users/rahulkulkarni/Downloads/Data/LGA_population_data.xlsx', sheet_name=None)

# Extract the relevant "LGA Percentage Data" sheet
latest_lga_percentage_data = latest_lga_population_data['LGA Percentage Data']

# Standardize the column names to upper case
latest_lga_percentage_data.rename(columns=lambda x: x.strip().upper(), inplace=True)

# Clean 'LGA NAME' in percentage data
latest_lga_percentage_data['LGA NAME'] = latest_lga_percentage_data['LGA FROM WIKIPEDIA'].str.upper()
latest_lga_percentage_data['LGA NAME'] = latest_lga_percentage_data['LGA NAME'].apply(clean_lga_name)

# -----------------------------------
# *** Start of Column Removal Section ***
# -----------------------------------

# Define the columns to remove (ensure they are in uppercase to match the standardized column names)
columns_to_remove = [
    '2007 ERP AGE 0-4, PERSONS',
    '2007 ERP AGE 0-4, %',
    '2007 ERP AGE 5-9, PERSONS',
    '2007 ERP AGE 5-9, %',
    '2007 ERP AGE 10-14, PERSONS',
    '2007 ERP AGE 10-14, %',
    '2007 ERP AGE 15-19, PERSONS',
    '2007 ERP AGE 15-19, %',
    '2007 ERP AGE 20-24, PERSONS',
    '2007 ERP AGE 20-24, %',
    '2007 ERP AGE 25-44, PERSONS',
    '2007 ERP AGE 25-44, %',
    '2007 ERP AGE 45-64, PERSONS',
    '2007 ERP AGE 45-64, %',
    '2007 ERP AGE 65-69, PERSONS',
    '2007 ERP AGE 65-69, %',
    '2007 ERP AGE 70-74, PERSONS',
    '2007 ERP AGE 70-74, %',
    '2007 ERP AGE 75-79, PERSONS',
    '2007 ERP AGE 75-79, %',
    '2007 ERP AGE 80-84, PERSONS',
    '2007 ERP AGE 80-84, %',
    '2007 ERP AGE 85+, PERSONS',
    '2007 ERP AGE 85+, %',
    '2007 ERP, TOTAL',
    '% CHANGE, 2007-2012, AGE 0-4',
    '% CHANGE, 2007-2012, AGE 5-9',
    '% CHANGE, 2007-2012, AGE 10-14',
    '% CHANGE, 2007-2012, AGE 15-19',
    '% CHANGE, 2007-2012, AGE 20-24',
    '% CHANGE, 2007-2012, AGE 25-44',
    '% CHANGE, 2007-2012, AGE 45-64',
    '% CHANGE, 2007-2012, AGE 65-69',
    '% CHANGE, 2007-2012, AGE 70-74',
    '% CHANGE, 2007-2012, AGE 75-79',
    '% CHANGE, 2007-2012, AGE 80-84',
    '% CHANGE, 2007-2012, AGE 85+',
    '% CHANGE, 2007-2012, TOTAL'
]

# Drop the specified columns from the percentage data
latest_lga_percentage_data.drop(columns=columns_to_remove, inplace=True, errors='ignore')

# -----------------------------------
# *** End of Column Removal Section ***
# -----------------------------------

# Identify all columns that represent percentages (typically columns containing "%")
percentage_columns = [col for col in latest_lga_percentage_data.columns if '%' in col]

# Select only the percentage columns and the 'LGA NAME'
latest_lga_percentage_data = latest_lga_percentage_data[['LGA NAME'] + percentage_columns]

# Merge the comprehensive percentage data with the EGM per capita data using a 'outer' join to include all LGAs
final_comprehensive_data = pd.merge(merged_data_2012, latest_lga_percentage_data, on='LGA NAME', how='outer')

# Check the number of unique LGAs after merging
print(f"\nTotal Unique LGAs after Merge: {final_comprehensive_data['LGA NAME'].nunique()}")

# Check missing values after merge
print("\nMissing Values After Merge:")
print(final_comprehensive_data.isnull().sum())

# Fill missing values in percentage columns with mean (imputation)
X = final_comprehensive_data[percentage_columns]
X.fillna(X.mean(), inplace=True)

# Fill missing 'PER CAPITA EGM SPENDING' if any (unlikely with 'outer' merge)
y = final_comprehensive_data['PER CAPITA EGM SPENDING'].fillna(final_comprehensive_data['PER CAPITA EGM SPENDING'].mean())

# Log transformation of the target variable to handle skewness
y_log = np.log1p(y)

# ---------------------------
# 2. Data Cleaning and Preprocessing
# ---------------------------

# (Already handled missing values and log transformation above)

# ---------------------------
# 3. Multicollinearity Analysis
# ---------------------------

# Function to calculate VIF
def calculate_vif(X):
    vif_data = pd.DataFrame()
    vif_data["feature"] = X.columns
    vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(len(X.columns))]
    return vif_data

# Initial VIF calculation
vif_data = calculate_vif(X)
print("\nInitial Variance Inflation Factors (VIF):")
print(vif_data.sort_values(by="VIF", ascending=False).head(10))  # Show top 10 for brevity

# Plotting pairwise correlation heatmap
plt.figure(figsize=(20, 15))
correlation_matrix = X.corr()
sns.heatmap(correlation_matrix, annot=False, cmap='coolwarm')
plt.title('Pairwise Correlation Heatmap')
plt.show()

# Iteratively remove features with high VIF
threshold = 15.0  # Adjusted threshold to prevent excessive feature removal
high_vif_features = vif_data[vif_data['VIF'] > threshold]['feature'].tolist()

while len(high_vif_features) > 0:
    feature_to_remove = vif_data.sort_values('VIF', ascending=False).iloc[0]['feature']
    print(f"Removing feature '{feature_to_remove}' with VIF: {vif_data['VIF'].max()}")
    X = X.drop(columns=[feature_to_remove])
    
    # Recalculate VIF
    vif_data = calculate_vif(X)
    high_vif_features = vif_data[vif_data['VIF'] > threshold]['feature'].tolist()

print("\nFinal VIF after removing multicollinear features:")
print(vif_data.sort_values(by="VIF", ascending=False))

# Optionally, visualize the reduced correlation matrix
plt.figure(figsize=(20, 15))
correlation_matrix = X.corr()
sns.heatmap(correlation_matrix, annot=False, cmap='coolwarm')
plt.title('Reduced Pairwise Correlation Heatmap')
plt.show()

# Define feature_names based on the remaining columns in X
feature_names = X.columns.tolist()

# ---------------------------
# 4. Dimensionality Reduction (PCA)
# ---------------------------

X.fillna(X.mean(), inplace=True)

# Fill missing values in y if any
y = final_comprehensive_data['PER CAPITA EGM SPENDING'].fillna(final_comprehensive_data['PER CAPITA EGM SPENDING'].mean())

# ---------------------------
# Address Outliers and Perform Scaling
# ---------------------------

# Apply Robust Scaling to X to handle outliers
scaler = RobustScaler()
X_scaled = scaler.fit_transform(X)

# Apply Power Transformation (Yeo-Johnson) to y to stabilize variance
pt = PowerTransformer(method='yeo-johnson')
y_transformed = pt.fit_transform(y.values.reshape(-1, 1)).flatten()

# ---------------------------
# PCA for Dimensionality Reduction
# ---------------------------

# Fit PCA to capture the maximum variance
pca = PCA()
X_pca_full = pca.fit_transform(X_scaled)

# Plot cumulative explained variance to determine the optimal number of components
cumulative_variance = np.cumsum(pca.explained_variance_ratio_)
plt.figure(figsize=(10, 6))
plt.plot(range(1, len(cumulative_variance) + 1), cumulative_variance, marker='o', linestyle='--')
plt.xlabel('Number of Components')
plt.ylabel('Cumulative Explained Variance')
plt.title('Cumulative Explained Variance by PCA Components')
plt.axhline(y=0.95, color='r', linestyle='-')  # 95% variance threshold
plt.text(0.5, 0.90, '95% Variance Threshold', color='red', fontsize=12)
plt.grid(True)
plt.show()

# Determine number of components to retain 95% variance
n_components_95 = np.argmax(cumulative_variance >= 0.95) + 1
print(f"Number of components to retain 95% variance: {n_components_95}")

# Transform data using optimal number of components
pca_optimal = PCA(n_components=n_components_95)
X_pca = pca_optimal.fit_transform(X_scaled)

# ---------------------------
# Train-Test Split
# ---------------------------

# Perform Train-Test Split (e.g., using 80% for training and 20% for testing)
X_train_pca, X_test_pca, y_train, y_test = train_test_split(
    X_pca, y_transformed, test_size=0.2, random_state=14, shuffle=True
)

print(f"Training set shape: {X_train_pca.shape}")
print(f"Testing set shape: {X_test_pca.shape}")

# ---------------------------
# Model Training and Cross-Validation (Using LOOCV)
# ---------------------------

# Initialize Ridge Regression model with regularization
ridge = Ridge(alpha=1.0)
loo = LeaveOneOut()

# Fit model using Leave-One-Out Cross-Validation
ridge.fit(X_train_pca, y_train)

# Get cross-validated predictions using LOOCV
y_pred_cv = cross_val_predict(ridge, X_pca, y_transformed, cv=loo)

# Inverse transform predictions and actual values to original scale
y_pred_original_cv = pt.inverse_transform(y_pred_cv.reshape(-1, 1)).flatten()
y_original = pt.inverse_transform(y_transformed.reshape(-1, 1)).flatten()

# ---------------------------
# Evaluate Performance of LOOCV Predictions
# ---------------------------

# ### NEW CODE STARTS HERE ###
import pandas as pd
import numpy as np

# Convert predictions to a pandas Series for easier manipulation
y_pred_series = pd.Series(y_pred_original_cv)

# Calculate Q1 (25th percentile) and Q3 (75th percentile)
Q1 = y_pred_series.quantile(0.25)
Q3 = y_pred_series.quantile(0.75)
IQR = Q3 - Q1

# Define lower and upper bounds for outliers
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR

print(f"\nOutlier Detection Using IQR:")
print(f"Q1: {Q1}")
print(f"Q3: {Q3}")
print(f"IQR: {IQR}")
print(f"Lower Bound: {lower_bound}")
print(f"Upper Bound: {upper_bound}")

# Identify outliers
outliers = y_pred_series[(y_pred_series < lower_bound) | (y_pred_series > upper_bound)]
print(f"\nNumber of Outliers Detected: {outliers.count()}")

# Replace outliers with NaN
y_pred_cleaned = y_pred_series.mask((y_pred_series < lower_bound) | (y_pred_series > upper_bound))

# Convert back to NumPy array
y_pred_cleaned = y_pred_cleaned.values

# Calculate residuals, which will have NaN where predictions were outliers
residuals = y_original - y_pred_cleaned

# Calculate performance metrics excluding NaN values
mask = ~np.isnan(y_pred_cleaned)

mse_cv = mean_squared_error(y_original[mask], y_pred_cleaned[mask])
mae_cv = mean_absolute_error(y_original[mask], y_pred_cleaned[mask])
r2_cv = r2_score(y_original[mask], y_pred_cleaned[mask])

print("\nCross-Validated Model Performance (Ridge Regression) After Removing Outliers:")
print("Mean Squared Error (MSE):", mse_cv)
print("Mean Absolute Error (MAE):", mae_cv)
print("R^2 Score:", r2_cv)
### NEW CODE ENDS HERE ###

# ---------------------------
# Identify Outliers and Cap Predictions
# ---------------------------

# Cap predictions to avoid extreme values beyond observed EGM spending
max_spending = y_original.max()
y_pred_capped = np.minimum(y_pred_original_cv, max_spending)

# Re-evaluate metrics after capping
mse_capped = mean_squared_error(y_original, y_pred_capped)
mae_capped = mean_absolute_error(y_original, y_pred_capped)
r2_capped = r2_score(y_original, y_pred_capped)

print("\nModel Performance after Capping Predictions:")
print("Capped Mean Squared Error (MSE):", mse_capped)
print("Capped Mean Absolute Error (MAE):", mae_capped)
print("Capped R^2 Score:", r2_capped)

# ---------------------------
# Residual Analysis
# ---------------------------

# Residuals for analysis
residuals = y_original - y_pred_capped

# Plot residuals vs. predicted values
plt.figure(figsize=(10, 6))
sns.scatterplot(x=y_pred_capped, y=residuals)
plt.axhline(0, color='red', linestyle='--')
plt.title('Residuals vs Predicted Values (Ridge Regression with Capped Predictions)')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
plt.show()

# ---------------------------
# Feature Importance Analysis using PCA Components
# ---------------------------

# Note: Since we are using PCA, feature importance corresponds to the principal components
feature_importance = ridge.coef_
feature_names_pca = [f'PC{i+1}' for i in range(n_components_95)]

# Create a DataFrame for feature importances
importance_df = pd.DataFrame({
    'Principal Component': feature_names_pca,
    'Importance': feature_importance
})

# Sort by absolute importance
importance_df['Absolute Importance'] = np.abs(importance_df['Importance'])
importance_df = importance_df.sort_values(by='Absolute Importance', ascending=False)

# Display feature importance
plt.figure(figsize=(12, 8))
sns.barplot(x='Importance', y='Principal Component', data=importance_df)
plt.title('Feature Importance from Ridge Regression (PCA Components)')
plt.xlabel('Importance')
plt.ylabel('Principal Component')
plt.show()

# ---------------------------
# Train and Evaluate an Ensemble Model (Voting Regressor Example)
# ---------------------------

# Initialize an ensemble model with Voting Regressor
ensemble = VotingRegressor([
    ('ridge', ridge),
    ('random_forest', RandomForestRegressor(n_estimators=100, random_state=42)),
    ('gradient_boosting', GradientBoostingRegressor(random_state=42))
])

# Fit the ensemble model using the training data
ensemble.fit(X_train_pca, y_train)

# Predict on the test set using the ensemble model
y_pred_ensemble = ensemble.predict(X_test_pca)
y_pred_ensemble_original = pt.inverse_transform(y_pred_ensemble.reshape(-1, 1)).flatten()

# Cap the ensemble predictions
y_pred_ensemble_capped = np.minimum(y_pred_ensemble_original, max_spending)

# Evaluate the ensemble model performance
ensemble_mse = mean_squared_error(y_test, y_pred_ensemble_capped)
ensemble_mae = mean_absolute_error(y_test, y_pred_ensemble_capped)
ensemble_r2 = r2_score(y_test, y_pred_ensemble_capped)

print("\nEnsemble Model Performance on Test Set (Capped Predictions):")
print("Mean Squared Error (MSE):", ensemble_mse)
print("Mean Absolute Error (MAE):", ensemble_mae)
print("R^2 Score:", ensemble_r2)

# ---------------------------
# Optional: Analyze PCA Loadings for Feature Interpretation
# ---------------------------

def analyze_pca_loadings(pca, feature_names, top_n=5):
    loadings = pd.DataFrame(pca.components_.T, index=feature_names, columns=[f'PC{i+1}' for i in range(pca.n_components_)])
    for i in range(pca.n_components_):
        pc = loadings.columns[i]
        sorted_features = loadings[pc].abs().sort_values(ascending=False)
        top_features = sorted_features.head(top_n).index.tolist()
        print(f"\nTop {top_n} features contributing to {pc}:")
        for feature in top_features:
            print(f"- {feature} (Loading: {loadings.loc[feature, pc]:.4f})")

# Analyze PCA Loadings using the defined feature_names
analyze_pca_loadings(pca_optimal, feature_names, top_n=5)
def remove_outliers_iqr(predictions, multiplier=1.5):
    """
    Removes outliers from predictions using the IQR method.
    
    Parameters:
        predictions (array-like): The predicted values.
        multiplier (float): The multiplier for IQR to define bounds.
        
    Returns:
        np.ndarray: Predictions with outliers replaced by NaN.
        int: Number of outliers detected.
    """
    pred_series = pd.Series(predictions)
    Q1 = pred_series.quantile(0.25)
    Q3 = pred_series.quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - multiplier * IQR
    upper_bound = Q3 + multiplier * IQR
    
    outliers = pred_series[(pred_series < lower_bound) | (pred_series > upper_bound)]
    num_outliers = outliers.count()
    
    pred_cleaned = pred_series.mask((pred_series < lower_bound) | (pred_series > upper_bound)).values
    
    return pred_cleaned, num_outliers

# Usage
y_pred_cleaned, num_outliers = remove_outliers_iqr(y_pred_original_cv, multiplier=1.5)
print(f"Number of Outliers Detected: {num_outliers}")

Now plot the residuals with the removed outliers