!pip install optuna

# Load datasets
career_batting = pd.read_csv('/content/career batting.csv')
contracts = pd.read_csv('/content/contracts.csv')
hitting = pd.read_csv('/content/hitting.csv')

# Step 1: Standardize and Align the 'Player' Column
# Fix Career Batting dataset
career_batting.rename(columns={'Name': 'Player'}, inplace=True)

# Fix Hitting dataset: Split 'last_name, first_name' into 'Player'
hitting['Player'] = hitting['last_name, first_name'].str.split(', ').str[::-1].str.join(' ')
hitting['Player'] = hitting['Player'].str.strip().str.title()

# Standardize 'Player' column in all datasets
contracts['Player'] = contracts['Player'].str.strip().str.title()
career_batting['Player'] = career_batting['Player'].str.strip().str.title()

# Step 2: Merge Datasets
data = pd.merge(hitting, contracts, on='Player', how='inner')
data = pd.merge(data, career_batting, on='Player', how='inner')

# Step 3: Display Merge Results
print(f"Merged dataset contains {data.shape[0]} rows and {data.shape[1]} columns.")
print(data.head())

# Rename columns for clarity and usability
data.rename(columns={
    'Team\n                        Currently With': 'Current Team',
    'Age\n                        At Signing': 'Age at Signing'
}, inplace=True)

# Display the updated column names to confirm changes
print(data.columns)

# Remove columns starting with 'x' (e.g., xWOBA, xwOBA, etc.)
columns_to_drop = [col for col in data.columns if col.lower().startswith('x')]
data = data.drop(columns=columns_to_drop, errors='ignore')

# Display updated dataset columns to confirm
print(f"Remaining columns in the dataset: {data.columns.tolist()}")

# Convert AAV to numeric by removing dollar signs and commas
data['AAV'] = data['AAV'].replace({'\$': '', ',': ''}, regex=True).astype(float)

# Plot the distribution of AAV
plt.figure(figsize=(10, 6))
plt.hist(data['AAV'], bins=30, edgecolor='k', alpha=0.7)
plt.title('Distribution of Average Annual Value (AAV)', fontsize=14)
plt.xlabel('AAV ($)', fontsize=12)
plt.ylabel('Frequency', fontsize=12)
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.show()

plt.figure(figsize=(10, 6))
sns.scatterplot(x=data['woba'], y=data['AAV'])
plt.title('Scatter Plot of AAV vs wOBA', fontsize=14)
plt.xlabel('wOBA', fontsize=12)
plt.ylabel('AAV ($)', fontsize=12)
plt.grid(alpha=0.3)
plt.show()

# Filter numeric columns for correlation calculation, excluding irrelevant features
numeric_data = data.select_dtypes(include=['float64', 'int64']).drop(columns=['player_id'], errors='ignore')

# Identify the top 10 features most correlated with AAV
correlation_with_aav = numeric_data.corr()['AAV'].abs().sort_values(ascending=False)
top_10_features = correlation_with_aav.index[1:11]  # Exclude 'AAV' itself

# Generate a heatmap for the top 10 features
plt.figure(figsize=(12, 8))
sns.heatmap(numeric_data[top_10_features].corr(), annot=True, cmap='coolwarm', fmt=".2f", linewidths=0.5)
plt.title('Correlation Matrix of Top 10 Features Most Correlated with AAV (Excluding Less Relevent Columns)', fontsize=14)
plt.show()

# Select the top features from the correlation heatmap for pairplot
top_features_correlation = ['barrel', 'home_run', 'HR', 'Off', 'RBI', 'R', 'WAR', 'PA', 'out_zone_swing', 'woba']

# Generate pairplot for the selected features
sns.pairplot(data[top_features_correlation], diag_kind='kde', plot_kws={'alpha': 0.7})
plt.suptitle('Pairwise Scatter Plots of Features from Correlation Heatmap', y=1.02, fontsize=14)
plt.show()

# Step 1: Drop irrelevant columns that won't affect clustering
columns_to_drop = ['Current Team', 'Value', 'Pos', 'player_id', 'Player']
numeric_data = data.drop(columns=columns_to_drop, errors='ignore')

# Step 2: Retain only numeric columns
numeric_data = numeric_data.select_dtypes(include=['float64', 'int64'])

# Step 3: Define X (features) and y (target - AAV for salary prediction)
X = numeric_data.drop(columns=['AAV'], errors='ignore')  # Drop 'AAV' from features
y = data['AAV']  # Target variable is the player's salary (AAV)

# Step 4: Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Step 5: Apply Ridge Regression with cross-validation to select the optimal alpha
ridge = RidgeCV(alphas=np.logspace(-4, 4, 100), cv=5).fit(X_train, y_train)

# Step 6: Output model performance
print(f"Optimal Alpha: {ridge.alpha_}")
print(f"Training Score: {ridge.score(X_train, y_train)}")
print(f"Testing Score: {ridge.score(X_test, y_test)}")

# Step 1: Standardize the features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("Features have been standardized.")

# Step 2: Determine the optimal number of clusters using the Elbow Method
wcss = []  # Within-cluster sum of squares
for i in range(1, 11):
    kmeans = KMeans(n_clusters=i, random_state=42)
    kmeans.fit(X_scaled)
    wcss.append(kmeans.inertia_)

# Plot the Elbow Method
plt.figure(figsize=(8, 5))
plt.plot(range(1, 11), wcss, marker='o')
plt.title('Elbow Method for Optimal Clusters')
plt.xlabel('Number of Clusters')
plt.ylabel('WCSS')
plt.show()

# Calculate cluster centers
cluster_centers = pd.DataFrame(kmeans.cluster_centers_, columns=X.columns)

# Calculate the range of each feature across clusters
feature_ranges = cluster_centers.max() - cluster_centers.min()

# Rank features by their range across clusters
cluster_feature_importance = feature_ranges.sort_values(ascending=False)
print("Feature Importance for Clusters (based on range):")
print(cluster_feature_importance)

# Select top 10 features by importance
top_features = cluster_feature_importance.head(10).index.tolist()

# Examine cluster centers for top features
top_cluster_centers = cluster_centers[top_features]
print("Cluster Centers for Top Features:")
print(top_cluster_centers)

# Visualize cluster centers to identify patterns
plt.figure(figsize=(10, 6))
top_cluster_centers.T.plot(kind='bar', figsize=(14, 8), colormap='viridis')
plt.title("Top Feature Averages Across Clusters", fontsize=16)
plt.xlabel("Features", fontsize=14)
plt.ylabel("Cluster Center Values", fontsize=14)
plt.legend(title="Cluster", bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()

# Add the numeric cluster labels from KMeans to the dataset
data['Cluster'] = cluster_labels  # This adds the numeric cluster labels (0, 1, 2) to the 'Cluster' column in the DataFrame

# Define a dictionary for cluster labels
cluster_labels_map = {
    0: "Power Hitters",      # Cluster 0 corresponds to Power Hitters
    1: "Balanced Hitters",   # Cluster 1 corresponds to Balanced Hitters
    2: "Utility Players"     # Cluster 2 corresponds to Utility Players
}

# Map the descriptive cluster labels to the numeric clusters
data['Cluster_Label'] = data['Cluster'].map(cluster_labels_map)

# Check if the mapping is successful by displaying a sample of the dataset
print(data[['Cluster', 'Cluster_Label']].drop_duplicates())

# Step 1: Predict salaries using the trained Ridge model
y_test_pred = ridge.predict(X_test)

# Step 2: Define a threshold for fairness (10% of average salary as an example)
threshold = 0.1 * y.mean()

# Step 3: Define salary fairness categories
def classify_salary(actual, predicted, threshold):
    if actual < predicted - threshold:
        return "Underpaid"
    elif actual > predicted + threshold:
        return "Overpaid"
    else:
        return "Fairly Paid"

# Step 4: Apply the classification logic to the test set
fairness_labels = [
    classify_salary(actual, predicted, threshold)
    for actual, predicted in zip(y_test, y_test_pred)
]

# Step 5: Combine the results into a new DataFrame
fairness_results = X_test.copy()
fairness_results['Actual Salary'] = y_test
fairness_results['Predicted Salary'] = y_test_pred
fairness_results['Fairness'] = fairness_labels

# Display a summary of the fairness classification
print(fairness_results[['Actual Salary', 'Predicted Salary', 'Fairness']].head(10))

# Plot the fairness distribution
plt.figure(figsize=(8, 5))
sns.countplot(data=fairness_results, x='Fairness', order=["Underpaid", "Fairly Paid", "Overpaid"])
plt.title("Salary Fairness Distribution", fontsize=16)
plt.xlabel("Fairness Category", fontsize=14)
plt.ylabel("Count", fontsize=14)
plt.tight_layout()
plt.show()

# Evaluate Ridge regression performance
mse = mean_squared_error(y_test, y_test_pred)
mae = mean_absolute_error(y_test, y_test_pred)
r2 = r2_score(y_test, y_test_pred)

print("Model Performance Metrics:")
print(f"Mean Squared Error (MSE): {mse:.2f}")
print(f"Mean Absolute Error (MAE): {mae:.2f}")
print(f"R-squared (R2): {r2:.2f}")

# Evaluate classification distribution
fairness_distribution = fairness_results['Fairness'].value_counts(normalize=True) * 100
print("\nFairness Classification Distribution:")
print(fairness_distribution)

# Visualize residuals
residuals = y_test - y_test_pred
plt.figure(figsize=(8, 5))
plt.scatter(y_test, residuals, alpha=0.7)
plt.axhline(y=0, color='r', linestyle='--')
plt.title("Residuals of Actual vs. Predicted Salary")
plt.xlabel("Actual Salary")
plt.ylabel("Residuals")
plt.show()

# Use RFE with Linear Regression
selector = RFE(estimator=LinearRegression(), n_features_to_select=10)  # Adjust the number of features
selector = selector.fit(X_train, y_train)

# Get selected features
selected_features = X.columns[selector.support_]
print("Selected Features by RFE:")
print(selected_features)

# Train the model using only selected features
X_train_selected = X_train[selected_features]
X_test_selected = X_test[selected_features]

linear_model_rfe = LinearRegression()
linear_model_rfe.fit(X_train_selected, y_train)
y_test_pred_rfe = linear_model_rfe.predict(X_test_selected)

# Evaluate the RFE model
mse_rfe = mean_squared_error(y_test, y_test_pred_rfe)
mae_rfe = mean_absolute_error(y_test, y_test_pred_rfe)
r2_rfe = r2_score(y_test, y_test_pred_rfe)

print("RFE Model Performance:")
print(f"Mean Squared Error (MSE): {mse_rfe:.2f}")
print(f"Mean Absolute Error (MAE): {mae_rfe:.2f}")
print(f"R-squared (R2): {r2_rfe:.2f}")

# Generate polynomial features
poly = PolynomialFeatures(degree=2, include_bias=False)
X_train_poly = poly.fit_transform(X_train)
X_test_poly = poly.transform(X_test)

# Train a linear model on polynomial features
poly_model = LinearRegression()
poly_model.fit(X_train_poly, y_train)
y_test_pred_poly = poly_model.predict(X_test_poly)

# Evaluate the polynomial model
mse_poly = mean_squared_error(y_test, y_test_pred_poly)
mae_poly = mean_absolute_error(y_test, y_test_pred_poly)
r2_poly = r2_score(y_test, y_test_pred_poly)

print("Polynomial Regression Model Performance:")
print(f"Mean Squared Error (MSE): {mse_poly:.2f}")
print(f"Mean Absolute Error (MAE): {mae_poly:.2f}")
print(f"R-squared (R2): {r2_poly:.2f}")

# Step 1: Apply Random Forest Regressor
rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
rf_model.fit(X_train, y_train)
y_pred_rf = rf_model.predict(X_test)

# Evaluate Random Forest Model
mse_rf = mean_squared_error(y_test, y_pred_rf)
mae_rf = mean_absolute_error(y_test, y_pred_rf)
r2_rf = r2_score(y_test, y_pred_rf)

print("Random Forest Model Performance:")
print(f"Mean Squared Error (MSE): {mse_rf:.2f}")
print(f"Mean Absolute Error (MAE): {mae_rf:.2f}")
print(f"R-squared (R2): {r2_rf:.2f}")

# Step 2: Apply Gradient Boosting Regressor
gb_model = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1, random_state=42)
gb_model.fit(X_train, y_train)
y_pred_gb = gb_model.predict(X_test)

# Evaluate Gradient Boosting Model
mse_gb = mean_squared_error(y_test, y_pred_gb)
mae_gb = mean_absolute_error(y_test, y_pred_gb)
r2_gb = r2_score(y_test, y_pred_gb)

print("\nGradient Boosting Model Performance:")
print(f"Mean Squared Error (MSE): {mse_gb:.2f}")
print(f"Mean Absolute Error (MAE): {mae_gb:.2f}")
print(f"R-squared (R2): {r2_gb:.2f}")

# Ensure only numeric columns are used for correlation
numeric_data = data.select_dtypes(include=['float64', 'int64'])

# Check correlations with AAV
if 'AAV' in numeric_data.columns:
    correlation_with_aav = numeric_data.corr()['AAV'].sort_values(ascending=False)
    print("Correlation of features with AAV:")
    print(correlation_with_aav.head(20))  # Display top 20 correlations
else:
    print("Error: 'AAV' column not found in the numeric data.")

# Step 2: Remove features with low correlation (e.g., less than 0.3)
low_corr_threshold = 0.3  # Adjust threshold as needed
low_corr_features = correlation_with_aav[correlation_with_aav.abs() < low_corr_threshold].index.tolist()
data = data.drop(columns=low_corr_features, errors='ignore')
print(f"Removed low correlation features: {low_corr_features}")

# Step 3: Add interaction terms for high-correlation features
data['home_run_RBI'] = data['home_run'] * data['RBI']  # Interaction of power metrics
data['barrel_WAR'] = data['barrel'] * data['WAR']      # Interaction of power and value metrics
data['PA_R'] = data['PA'] * data['R']                 # Interaction of plate appearances and runs

# Step 4: Prepare dataset again for modeling
numeric_data = data.select_dtypes(include=['float64', 'int64'])
X = numeric_data.drop(columns=['AAV'], errors='ignore')  # Features
y = data['AAV']  # Target variable

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Step 5: Train and evaluate Random Forest and Gradient Boosting again

# Random Forest
rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
rf_model.fit(X_train, y_train)
y_pred_rf = rf_model.predict(X_test)
print("\nRandom Forest Updated Metrics:")
print(f"R-squared: {r2_score(y_test, y_pred_rf):.2f}")
print(f"MAE: {mean_absolute_error(y_test, y_pred_rf):.2f}")
print(f"MSE: {mean_squared_error(y_test, y_pred_rf):.2f}")

# Gradient Boosting
gb_model = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1, random_state=42)
gb_model.fit(X_train, y_train)
y_pred_gb = gb_model.predict(X_test)
print("\nGradient Boosting Updated Metrics:")
print(f"R-squared: {r2_score(y_test, y_pred_gb):.2f}")
print(f"MAE: {mean_absolute_error(y_test, y_pred_gb):.2f}")
print(f"MSE: {mean_squared_error(y_test, y_pred_gb):.2f}")

# Random Forest Parameter Grid
rf_param_grid = {
    'n_estimators': [100, 200, 300],
    'max_depth': [None, 10, 20, 30],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4]
}

# Gradient Boosting Parameter Grid
gb_param_grid = {
    'n_estimators': [100, 200, 300],
    'learning_rate': [0.01, 0.1, 0.2],
    'max_depth': [3, 5, 7],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4]
}

# Random Forest Randomized Search
rf_random_search = RandomizedSearchCV(
    estimator=RandomForestRegressor(random_state=42),
    param_distributions=rf_param_grid,
    n_iter=20,  # Number of random combinations to try
    scoring='neg_mean_squared_error',
    cv=3,
    n_jobs=-1,
    random_state=42,
    verbose=2
)
rf_random_search.fit(X_train, y_train)
print("Best Random Forest Parameters:", rf_random_search.best_params_)
print("Best Random Forest Score (MSE):", -rf_random_search.best_score_)

# Gradient Boosting Randomized Search
gb_random_search = RandomizedSearchCV(
    estimator=GradientBoostingRegressor(random_state=42),
    param_distributions=gb_param_grid,
    n_iter=20,  # Number of random combinations to try
    scoring='neg_mean_squared_error',
    cv=3,
    n_jobs=-1,
    random_state=42,
    verbose=2
)
gb_random_search.fit(X_train, y_train)
print("\nBest Gradient Boosting Parameters:", gb_random_search.best_params_)
print("Best Gradient Boosting Score (MSE):", -gb_random_search.best_score_)

# Evaluate Random Forest
y_pred_rf = rf_random_search.best_estimator_.predict(X_test)

mse_rf = mean_squared_error(y_test, y_pred_rf)
rmse_rf = np.sqrt(mse_rf)
mae_rf = mean_absolute_error(y_test, y_pred_rf)
mape_rf = mean_absolute_percentage_error(y_test, y_pred_rf)
r2_rf = r2_score(y_test, y_pred_rf)
explained_variance_rf = explained_variance_score(y_test, y_pred_rf)

print("\nRandom Forest Final Metrics:")
print(f"MSE: {mse_rf:.2f}")
print(f"RMSE: {rmse_rf:.2f}")
print(f"MAE: {mae_rf:.2f}")
print(f"MAPE: {mape_rf:.2%}")  # Show as percentage
print(f"R-squared: {r2_rf:.2f}")
print(f"Explained Variance Score: {explained_variance_rf:.2f}")

# Evaluate Gradient Boosting
y_pred_gb = gb_random_search.best_estimator_.predict(X_test)

mse_gb = mean_squared_error(y_test, y_pred_gb)
rmse_gb = np.sqrt(mse_gb)
mae_gb = mean_absolute_error(y_test, y_pred_gb)
mape_gb = mean_absolute_percentage_error(y_test, y_pred_gb)
r2_gb = r2_score(y_test, y_pred_gb)
explained_variance_gb = explained_variance_score(y_test, y_pred_gb)

print("\nGradient Boosting Final Metrics:")
print(f"MSE: {mse_gb:.2f}")
print(f"RMSE: {rmse_gb:.2f}")
print(f"MAE: {mae_gb:.2f}")
print(f"MAPE: {mape_gb:.2%}")  # Show as percentage
print(f"R-squared: {r2_gb:.2f}")
print(f"Explained Variance Score: {explained_variance_gb:.2f}")

# Step 1: Ensure Clusters and Fairness Labels Exist
# Add the clusters from the K-Means model
data['Cluster'] = kmeans.labels_  # Use the labels from the K-Means model
cluster_labels_map = {
    0: "Power Hitters",      # High in RBI, HR, SLG, and related metrics
    1: "Balanced Hitters",   # Moderate across all metrics
    2: "Utility Players"     # Lower contributions across all metrics
}
data['Cluster_Label'] = data['Cluster'].map(cluster_labels_map)

# Add the 'Fairness' column if it doesn't exist
if 'Predicted Salary' not in data.columns:
    data['Predicted Salary'] = rf_model.predict(data[top_features + ['Cluster']])

threshold = 0.1 * data['AAV'].mean()

def classify_salary(actual, predicted, threshold):
    if actual < predicted - threshold:
        return "Underpaid"
    elif actual > predicted + threshold:
        return "Overpaid"
    else:
        return "Fairly Paid"

data['Fairness'] = [
    classify_salary(actual, predicted, threshold)
    for actual, predicted in zip(data['AAV'], data['Predicted Salary'])
]

# Map Fairness to numeric labels for SVM
data['Fairness_Label'] = data['Fairness'].map({'Fairly Paid': 0, 'Underpaid': 1, 'Overpaid': 2})

# Step 2: Prepare the SVM Data
X = data[top_features + ['Cluster']]  # Include cluster information as a feature
y = data['Fairness_Label']

# Split the dataset
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)

# Step 3: Simplified Random Search with Reduced Trials
reduced_svm_param_grid = {
    'C': [0.1, 1, 10],        # Limited options for regularization strength
    'kernel': ['linear', 'rbf'],  # Focus on common kernels
    'gamma': ['scale']         # Use default gamma for simplicity
}

svm_model = SVC(random_state=42)
svm_random_search = RandomizedSearchCV(
    estimator=svm_model,
    param_distributions=reduced_svm_param_grid,
    n_iter=5,  # Fewer trials for faster results
    scoring='accuracy',
    cv=2,  # Reduce cross-validation folds
    n_jobs=-1,
    random_state=42,
    verbose=1  # Reduce verbosity
)
svm_random_search.fit(X_train, y_train)

# Step 4: Evaluate the SVM Model
best_svm = svm_random_search.best_estimator_
y_pred = best_svm.predict(X_test)

print("Simplified Random Search Best Parameters:", svm_random_search.best_params_)
print("\nClassification Report:")
print(classification_report(y_test, y_pred, target_names=['Fairly Paid', 'Underpaid', 'Overpaid']))

# Step 5: Confusion Matrix
conf_matrix = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(8, 6))
sns.heatmap(conf_matrix, annot=True, fmt="d", cmap="Blues",
            xticklabels=['Fairly Paid', 'Underpaid', 'Overpaid'],
            yticklabels=['Fairly Paid', 'Underpaid', 'Overpaid'])
plt.title("Confusion Matrix")
plt.xlabel("Predicted Label")
plt.ylabel("True Label")
plt.show()

# Step 6: Visualize Fairness by Cluster
fairness_cluster_summary = data.groupby(['Cluster', 'Fairness']).size().unstack(fill_value=0)
fairness_cluster_summary.plot(kind='bar', stacked=True, figsize=(12, 6), colormap='viridis')
plt.title("Fairness Categories by Clusters")
plt.xlabel("Player Cluster")
plt.ylabel("Number of Players")
plt.legend(title="Fairness Category")
plt.tight_layout()
plt.show()

# Step 1: Ensure Clusters and Fairness Labels Exist
data['Fairness_Label'] = data['Fairness'].map({'Fairly Paid': 0, 'Underpaid': 1, 'Overpaid': 2})
X = data[top_features + ['Cluster']]  # Use top features and cluster as predictors
y = data['Fairness_Label']

# Step 2: Split the Data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)

# Step 3: Random Forest with Random Search
rf_param_grid = {
    'n_estimators': [100, 200, 300],
    'max_depth': [None, 10, 20],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4]
}

rf_random_search = RandomizedSearchCV(
    estimator=RandomForestClassifier(random_state=42),
    param_distributions=rf_param_grid,
    n_iter=10,  # Fewer iterations for faster results
    scoring='accuracy',
    cv=3,
    n_jobs=-1,
    random_state=42,
    verbose=2
)
rf_random_search.fit(X_train, y_train)
rf_best = rf_random_search.best_estimator_

# Random Forest Evaluation
rf_y_pred = rf_best.predict(X_test)
print("\nRandom Forest Classification Report:")
print(classification_report(y_test, rf_y_pred, target_names=['Fairly Paid', 'Underpaid', 'Overpaid']))

# Confusion Matrix for Random Forest
conf_matrix_rf = confusion_matrix(y_test, rf_y_pred)
plt.figure(figsize=(8, 6))
sns.heatmap(conf_matrix_rf, annot=True, fmt="d", cmap="Blues",
            xticklabels=['Fairly Paid', 'Underpaid', 'Overpaid'],
            yticklabels=['Fairly Paid', 'Underpaid', 'Overpaid'])
plt.title("Confusion Matrix - Random Forest")
plt.xlabel("Predicted Label")
plt.ylabel("True Label")
plt.show()

# Step 4: Gradient Boosting with Random Search
gb_param_grid = {
    'n_estimators': [100, 200, 300],
    'learning_rate': [0.01, 0.1],
    'max_depth': [3, 5, 7],
    'min_samples_split': [2, 5],
    'min_samples_leaf': [1, 2]
}

gb_random_search = RandomizedSearchCV(
    estimator=GradientBoostingClassifier(random_state=42),
    param_distributions=gb_param_grid,
    n_iter=10,  # Fewer iterations for faster results
    scoring='accuracy',
    cv=3,
    n_jobs=-1,
    random_state=42,
    verbose=2
)
gb_random_search.fit(X_train, y_train)
gb_best = gb_random_search.best_estimator_

# Gradient Boosting Evaluation
gb_y_pred = gb_best.predict(X_test)
print("\nGradient Boosting Classification Report:")
print(classification_report(y_test, gb_y_pred, target_names=['Fairly Paid', 'Underpaid', 'Overpaid']))

# Confusion Matrix for Gradient Boosting
conf_matrix_gb = confusion_matrix(y_test, gb_y_pred)
plt.figure(figsize=(8, 6))
sns.heatmap(conf_matrix_gb, annot=True, fmt="d", cmap="Blues",
            xticklabels=['Fairly Paid', 'Underpaid', 'Overpaid'],
            yticklabels=['Fairly Paid', 'Underpaid', 'Overpaid'])
plt.title("Confusion Matrix - Gradient Boosting")
plt.xlabel("Predicted Label")
plt.ylabel("True Label")
plt.show()

# Step 5: Fairness Distribution by Clusters
fairness_cluster_summary = data.groupby(['Cluster', 'Fairness']).size().unstack(fill_value=0)
fairness_cluster_summary.plot(kind='bar', stacked=True, figsize=(12, 6), colormap='viridis')
plt.title("Fairness Categories by Clusters")
plt.xlabel("Player Cluster")
plt.ylabel("Number of Players")
plt.legend(title="Fairness Category")
plt.tight_layout()
plt.show()

# Define the mapping for clusters to player categories
cluster_labels_map = {
    0: "Power Hitters",
    1: "Balanced Hitters",
    2: "Utility Players"
}

# Replace the cluster numbers in the table with their corresponding labels
fairness_cluster_summary.index = fairness_cluster_summary.index.map(cluster_labels_map)

# Display the updated table
print(fairness_cluster_summary)

fairness_proportions = fairness_cluster_summary.div(fairness_cluster_summary.sum(axis=1), axis=0)
print(fairness_proportions)

fairness_cluster_summary.plot(kind='bar', stacked=True, figsize=(12, 6), colormap='viridis')
plt.title("Fairness Categories by Clusters")
plt.xlabel("Player Cluster")
plt.ylabel("Number of Players")
plt.legend(title="Fairness Category")
plt.tight_layout()
plt.show()

sns.boxplot(x='Cluster', y='AAV', hue='Fairness', data=data)
plt.title("Actual Salary Distribution by Cluster and Fairness")
plt.show()

sns.boxplot(x='Cluster', y='Predicted Salary', hue='Fairness', data=data)
plt.title("Predicted Salary Distribution by Cluster and Fairness")
plt.show()

# Step 1: Normalize the feature set
scaler = StandardScaler()
numeric_features = ['barrel', 'home_run', 'HR', 'RBI', 'SLG', 'Predicted Salary', 'Cluster']  # Include relevant features
X = data[numeric_features]
X_scaled = scaler.fit_transform(X)

# Step 2: Encode the target variable (Fairness Labels)
label_encoder = LabelEncoder()
y_encoded = label_encoder.fit_transform(data['Fairness_Label'])  # Ensure 'Fairness_Label' is numeric (0, 1, 2)
y_categorical = to_categorical(y_encoded)  # Convert to one-hot encoding for multi-class classification

# Step 3: Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_categorical, test_size=0.2, random_state=42, stratify=y_encoded)

print("Data prepared:")
print(f"Training set: {X_train.shape}, {y_train.shape}")
print(f"Testing set: {X_test.shape}, {y_test.shape}")

# Step 1: Correct column alignment for X_train and X_test
X_train_df = pd.DataFrame(X_train, columns=data[top_features].columns[:X_train.shape[1]])
X_test_df = pd.DataFrame(X_test, columns=data[top_features].columns[:X_test.shape[1]])

# Ensure 'Cluster' is included in the datasets
X_train_df['Cluster'] = data.loc[X_train_df.index, 'Cluster'].values
X_test_df['Cluster'] = data.loc[X_test_df.index, 'Cluster'].values

# Convert back to numpy arrays
X_train_with_clusters = X_train_df.values
X_test_with_clusters = X_test_df.values

# Step 2: Build the Neural Network Model
model = Sequential([
    Dense(64, input_dim=X_train_with_clusters.shape[1], activation='relu'),  # Input layer with 64 neurons
    Dropout(0.3),  # Dropout to prevent overfitting
    Dense(32, activation='relu'),  # Hidden layer with 32 neurons
    Dropout(0.3),
    Dense(16, activation='relu'),  # Another hidden layer with 16 neurons
    Dense(y_train.shape[1], activation='softmax')  # Output layer with softmax for multi-class classification
])

# Compile the model
model.compile(optimizer='adam',  # Adam optimizer
              loss='categorical_crossentropy',  # Loss for multi-class classification
              metrics=['accuracy'])

print(model.summary())

# Step 3: Train the Model
history = model.fit(X_train_with_clusters, y_train,
                    validation_split=0.2,
                    epochs=50,  # Adjust epochs as needed
                    batch_size=16,  # Mini-batch size
                    verbose=1)

# Step 4: Evaluate the Model
loss, accuracy = model.evaluate(X_test_with_clusters, y_test)
print(f"\nTest Loss: {loss}")
print(f"Test Accuracy: {accuracy}")

# Step 5: Visualize Training Progress
plt.figure(figsize=(12, 6))
plt.plot(history.history['accuracy'], label='Training Accuracy')
plt.plot(history.history['val_accuracy'], label='Validation Accuracy')
plt.title('Model Accuracy')
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.legend()
plt.show()

plt.figure(figsize=(12, 6))
plt.plot(history.history['loss'], label='Training Loss')
plt.plot(history.history['val_loss'], label='Validation Loss')
plt.title('Model Loss')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()
plt.show()

# Step 6: Analyze Results by Cluster
# Predict fairness categories for the test set
y_pred = np.argmax(model.predict(X_test_with_clusters), axis=1)

# Map predictions back to fairness categories
fairness_pred = pd.Series(y_pred).map({0: 'Fairly Paid', 1: 'Underpaid', 2: 'Overpaid'})

# Add predictions and clusters to a summary dataframe
results_df = pd.DataFrame({
    'Cluster': X_test_df['Cluster'],
    'True Fairness': pd.Series(np.argmax(y_test, axis=1)).map({0: 'Fairly Paid', 1: 'Underpaid', 2: 'Overpaid'}),
    'Predicted Fairness': fairness_pred
})

# Analyze fairness distribution by cluster
fairness_by_cluster = results_df.groupby('Cluster')['Predicted Fairness'].value_counts(normalize=True).unstack()
print(fairness_by_cluster)

# Visualize fairness distribution by cluster
fairness_by_cluster.plot(kind='bar', stacked=True, figsize=(12, 6), colormap='viridis')
plt.title('Fairness Prediction Distribution by Cluster')
plt.xlabel('Player Cluster')
plt.ylabel('Proportion')
plt.legend(title='Fairness Category')
plt.tight_layout()
plt.show()

!cp "/content/drive/MyDrive/Colab Notebooks/final_project.ipynb" ./
!jupyter nbconvert --to html "final_project.ipynb"