# Import Necessary Libraries
import pandas as pd
import numpy as np
import os
import warnings
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score
from sklearn.tree import DecisionTreeRegressor, plot_tree, export_graphviz
from sklearn.ensemble import RandomForestRegressor, StackingRegressor
from sklearn.linear_model import Ridge
from sklearn.neighbors import KNeighborsRegressor
from sklearn.metrics import mean_squared_error, r2_score
import graphviz

# Suppress Warnings for Cleaner Output
warnings.filterwarnings('ignore')

# Function to Check Graphviz Installation
def check_graphviz_installed():
    try:
        if os.name == 'nt':  # Windows
            from shutil import which
            if which("dot") is not None:
                print("Graphviz is installed correctly.")
            else:
                print("Please install Graphviz from https://graphviz.org/download/ and add it to your PATH.")
        else:
            # Unix/Linux/Mac
            exit_code = os.system("dot -V")
            if exit_code == 0:
                print("Graphviz is installed correctly.")
            else:
                print("Please install Graphviz from https://graphviz.org/download/")
    except Exception as e:
        print(f"Error checking Graphviz installation: {e}")

# Check Graphviz Installation
check_graphviz_installed()

# 1. Data Loading and Inspection
def load_data(filepath):
    try:
        df = pd.read_csv(filepath)
        print("\nData successfully loaded.")
        print("\nFirst 5 rows of the dataset:")
        print(df.head())
        return df
    except FileNotFoundError:
        print(f"\nFile not found at {filepath}. Please check the path and try again.")
        exit()
    except Exception as e:
        print(f"\nAn error occurred while loading the data: {e}")
        exit()

# 2. Data Preprocessing
def preprocess_data(df):
    print("\n--- Data Preprocessing ---\n")
    
    # Drop Irrelevant Columns
    columns_to_drop = ['LGA Name', 'LGA']
    df = df.drop(columns=columns_to_drop, axis=1, errors='ignore')
    print(f"Dropped columns: {columns_to_drop}")
    
    # Check for Duplicates
    duplicate_count = df.duplicated().sum()
    if duplicate_count > 0:
        print(f"Found {duplicate_count} duplicate rows. Removing duplicates.")
        df = df.drop_duplicates()
    else:
        print("No duplicate rows found.")
    
    # Identify and Remove Features Identical to Target (Data Leakage)
    target = 'EGM expenditure per capita'
    identical_features = [col for col in df.columns if col != target and df[col].equals(df[target])]
    if identical_features:
        print(f"\nWarning: The following feature(s) are identical to the target and will be removed to prevent data leakage: {identical_features}")
        df = df.drop(columns=identical_features)
    else:
        print("\nNo features identical to the target found.")
    
    # Encode Categorical Variables (if any)
    categorical_columns = df.select_dtypes(include=['object', 'category']).columns.tolist()
    if 'Region' in categorical_columns:
        df = pd.get_dummies(df, columns=['Region'], drop_first=True)
        print("\nOne-hot encoded 'Region' column.")
    else:
        print("\nNo categorical variables found to encode.")
    
    return df

# 3. Exploratory Data Analysis (EDA)
def perform_eda(df, features, target):
    print("\n--- Exploratory Data Analysis (EDA) ---\n")
    
    # Descriptive Statistics
    print("Descriptive Statistics:")
    print(df[features + [target]].describe())
    
    # Check Missing Values
    print("\nMissing Values:")
    print(df[features + [target]].isnull().sum())
    
    # Correlation Matrix
    plt.figure(figsize=(20, 15))
    correlation_matrix = df[features + [target]].corr()
    sns.heatmap(correlation_matrix, annot=True, fmt=".2f", cmap='coolwarm')
    plt.title('Correlation Matrix')
    plt.show()
    
    # Identify Highly Correlated Features (Correlation > 0.8)
    upper_tri = correlation_matrix.where(np.triu(np.ones(correlation_matrix.shape), k=1).astype(bool))
    high_corr_features = [column for column in upper_tri.columns if any(upper_tri[column] > 0.8)]
    if high_corr_features:
        print(f"\nHighly correlated features identified (Correlation > 0.8): {high_corr_features}")
    else:
        print("\nNo highly correlated features identified.")
    
    return high_corr_features

# 4. Feature Engineering and Selection
def feature_engineering(df, features, high_corr_features):
    print("\n--- Feature Engineering and Selection ---\n")
    
    # Remove Highly Correlated Features
    if high_corr_features:
        print(f"Removing highly correlated features: {high_corr_features}")
        features = [feature for feature in features if feature not in high_corr_features]
    else:
        print("No highly correlated features to remove.")
    
    # Additional Feature Engineering (Optional)
    # Example: Creating Interaction Terms
    # df['Income_IRSD'] = df['Personal income <$400/week, persons'] * df['IRSD (avg)']
    # features.append('Income_IRSD')
    
    return features

# 5. Model Training and Evaluation with Stacking Ensemble
def train_and_evaluate_models(X_train, X_test, y_train, y_test, features):
    print("\n--- Training and Evaluating Models ---\n")
    
    # Define Base Models for Stacking
    base_models = [
        ('rf', RandomForestRegressor(n_estimators=100, random_state=42)),
        ('dt', DecisionTreeRegressor(random_state=42)),
        ('knn', KNeighborsRegressor())
    ]
    
    # Define Meta Learner
    meta_learner = Ridge()
    
    # Initialize Stacking Regressor
    stacking_regressor = StackingRegressor(
        estimators=base_models,
        final_estimator=meta_learner,
        passthrough=False,
        cv=5,
        n_jobs=-1
    )
    
    # Train Stacking Regressor
    stacking_regressor.fit(X_train, y_train)
    print("Stacking Regressor trained successfully.")
    
    # Make Predictions
    y_pred_stacking = stacking_regressor.predict(X_test)
    
    # Evaluate Stacking Regressor
    mse_stacking = mean_squared_error(y_test, y_pred_stacking)
    r2_stacking = r2_score(y_test, y_pred_stacking)
    print(f"\nStacking Regressor - MSE: {mse_stacking:.2f}, R² Score: {r2_stacking:.2f}")
    
    # Feature Importances from Random Forest (One of the Base Models)
    rf_model = stacking_regressor.named_estimators_['rf']
    importances_rf = rf_model.feature_importances_
    feature_importances_rf = pd.DataFrame({
        'Feature': features,
        'Importance': importances_rf
    }).sort_values(by='Importance', ascending=False)
    print("\nRandom Forest Feature Importances within Stacking Regressor:")
    print(feature_importances_rf)
    
    # Visualize Feature Importances
    plt.figure(figsize=(12, 10))
    sns.barplot(x='Importance', y='Feature', data=feature_importances_rf)
    plt.title('Random Forest Feature Importances (Stacking Regressor)')
    plt.xlabel('Importance')
    plt.ylabel('Feature')
    plt.show()
    
    # Residual Analysis for Stacking Regressor
    residuals = y_test - y_pred_stacking
    residuals_df = pd.DataFrame({
        'Actual': y_test,
        'Predicted': y_pred_stacking,
        'Residual': residuals
    }).reset_index(drop=True)
    
    print("\nStacking Regressor Residuals (First 5 Instances):")
    print(residuals_df.head())
    
    # Scatter Plot of Residuals
    plt.figure(figsize=(12, 6))
    sns.scatterplot(x=y_pred_stacking, y=residuals, alpha=0.7)
    plt.axhline(y=0, color='r', linestyle='--')
    plt.xlabel('Predicted Values (Stacking Regressor)')
    plt.ylabel('Residuals')
    plt.title('Stacking Regressor Residuals')
    plt.show()
    
    return stacking_regressor

# 6. Hyperparameter Tuning for Stacking Regressor
def hyperparameter_tuning_stacking(X_train, y_train):
    print("\n--- Hyperparameter Tuning for Stacking Regressor ---\n")
    
    # Define Parameter Grid
    param_grid = {
        'rf__n_estimators': [100, 200],
        'rf__max_depth': [None, 10, 20],
        'dt__max_depth': [None, 10, 20],
        'knn__n_neighbors': [5, 10, 15],
        'final_estimator__alpha': [0.1, 1.0, 10.0]
    }
    
    # Initialize Stacking Regressor
    base_models = [
        ('rf', RandomForestRegressor(random_state=42)),
        ('dt', DecisionTreeRegressor(random_state=42)),
        ('knn', KNeighborsRegressor())
    ]
    meta_learner = Ridge()
    stacking_regressor = StackingRegressor(
        estimators=base_models,
        final_estimator=meta_learner,
        passthrough=False,
        cv=5,
        n_jobs=-1
    )
    
    # Initialize GridSearchCV
    grid_search = GridSearchCV(
        estimator=stacking_regressor,
        param_grid=param_grid,
        cv=3,
        scoring='r2',
        n_jobs=-1,
        verbose=2
    )
    
    # Perform Grid Search
    grid_search.fit(X_train, y_train)
    print(f"\nBest Parameters from Grid Search: {grid_search.best_params_}")
    
    # Best Estimator
    best_stacking = grid_search.best_estimator_
    print("\nBest Stacking Regressor Model trained successfully.")
    
    return best_stacking

# 7. Cross-Validation for Best Stacking Regressor
def perform_cross_validation(model, X_train, y_train):
    print("\n--- Cross-Validation for Best Stacking Regressor ---\n")
    cv_scores = cross_val_score(model, X_train, y_train, cv=5, scoring='r2', n_jobs=-1)
    print(f"Cross-Validation R² Scores: {cv_scores}")
    print(f"Average CV R² Score: {cv_scores.mean():.2f}")

# 8. Main Function to Execute All Steps
def main():
    # Define Data Path
    data_path = '/Users/rahulkulkarni/Downloads/Data/Proportional_Data_CSV.csv'
    
    # 1. Load Data
    df = load_data(data_path)
    
    # 2. Preprocess Data
    df = preprocess_data(df)
    
    # 3. Define Feature Columns
    feature_columns = [
        'IRSD (avg)',
        'Personal income <$400/week, persons',
        'Number of families',
        'Female-headed lone parent families',
        'Male-headed lone parent families',
        'Holds degree or higher, persons',
        'Did not complete year 12, persons',
        'Unemployed, persons',
        'Aboriginal or Torres Strait Islander, persons',
        'Born in non-English speaking country, persons',
        'population',
        '2012 ERP age 0-4, persons',
        '2012 ERP age 5-9, persons',
        '2012 ERP age 10-14, persons',
        '2012 ERP age 15-19, persons',
        '2012 ERP age 20-24, persons',
        '2012 ERP age 25-44, persons',
        '2012 ERP age 45-64, persons',
        '2012 ERP age 65-69, persons',
        '2012 ERP age 70-74, persons',
        '2012 ERP age 75-79, persons',
        '2012 ERP age 80-84, persons',
        '2012 ERP age 85+, persons'
    ]
    
    # Include One-Hot Encoded 'Region_' Columns
    categorical_features = [col for col in df.columns if 'Region_' in col]
    feature_columns.extend(categorical_features)
    
    # Check if Target Variable Exists and is Correct
    target = 'EGM expenditure per capita'
    if target not in df.columns:
        print(f"Target variable '{target}' not found in the dataset.")
        exit()
    
    # **Important:** Ensure that you are NOT overwriting the target variable with random values.
    # If present, remove or comment out the following line:
    # df[target] = np.random.rand(len(df)) * 1000  # Remove this line!
    
    # Feature Selection (Remove any missing features and duplicates were already handled)
    available_features = [col for col in feature_columns if col in df.columns]
    missing_features = set(feature_columns) - set(available_features)
    if missing_features:
        print(f"\nWarning: The following features are missing and will be ignored: {missing_features}")
    
    # Set Features and Target
    X = df[available_features]
    y = df[target]
    
    print(f"\nSelected Features ({len(available_features)}):")
    print(available_features)
    
    # 4. Perform EDA
    high_corr_features = perform_eda(df, available_features, target)
    
    # 5. Feature Engineering and Selection
    selected_features = feature_engineering(df, available_features, high_corr_features)
    
    # 6. Split Data
    X_selected = X[selected_features]
    X_train, X_test, y_train, y_test = train_test_split(
        X_selected,
        y,
        test_size=0.2,
        random_state=42
    )
    
    print(f"\nTraining Set: {X_train.shape}, Testing Set: {X_test.shape}")
    
    # 7. Train and Evaluate Initial Stacking Regressor
    stacking_regressor = train_and_evaluate_models(
        X_train,
        X_test,
        y_train,
        y_test,
        selected_features
    )
    
    # 8. Hyperparameter Tuning for Stacking Regressor
    best_stacking = hyperparameter_tuning_stacking(X_train, y_train)
    
    # 9. Evaluate Tuned Stacking Regressor
    print("\n--- Evaluating Tuned Stacking Regressor ---\n")
    y_pred_best_stacking = best_stacking.predict(X_test)
    mse_best_stacking = mean_squared_error(y_test, y_pred_best_stacking)
    r2_best_stacking = r2_score(y_test, y_pred_best_stacking)
    print(f"Tuned Stacking Regressor - MSE: {mse_best_stacking:.2f}, R² Score: {r2_best_stacking:.2f}")
    
    # Residual Analysis for Tuned Stacking Regressor
    residuals = y_test - y_pred_best_stacking
    residuals_df = pd.DataFrame({
        'Actual': y_test,
        'Predicted': y_pred_best_stacking,
        'Residual': residuals
    }).reset_index(drop=True)
    
    print("\nTuned Stacking Regressor Residuals (First 5 Instances):")
    print(residuals_df.head())
    
    # Scatter Plot of Residuals
    plt.figure(figsize=(12, 6))
    sns.scatterplot(x=y_pred_best_stacking, y=residuals, alpha=0.7)
    plt.axhline(y=0, color='r', linestyle='--')
    plt.xlabel('Predicted Values (Tuned Stacking Regressor)')
    plt.ylabel('Residuals')
    plt.title('Tuned Stacking Regressor Residuals')
    plt.show()
    
    # 10. Cross-Validation for Tuned Stacking Regressor
    perform_cross_validation(best_stacking, X_train, y_train)
    
    print("\n--- Model Training and Evaluation Complete ---\n")
    
    # Optional: Save the Best Stacking Regressor Model
    # import joblib
    # joblib.dump(best_stacking, 'best_stacking_regressor_model.pkl')
    # print("Best Stacking Regressor model saved as 'best_stacking_regressor_model.pkl'.")

# Execute the Main Function
if __name__ == '__main__':
    main()
How many datapoints are requred for a decision tree to output good results. I have 50 sets of data and when I do the decision tree or anything else no discernible trend is observed