import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import StratifiedKFold  # Changed from KFold
from sklearn.metrics import accuracy_score
from sklearn.neighbors import NearestNeighbors
from scipy.stats import rankdata
import warnings
import joblib  # Added for model saving
warnings.filterwarnings('ignore')

# BMFK Classifier (custom implementation)
class BMFK:
    def __init__(self, n_neighbors=5, m=2, p=2, q=2):
        self.n_neighbors = n_neighbors
        self.m = m
        self.p = p
        self.q = q

    def fit(self, X, y):
        self.X = X
        self.y = y
        self.classes = np.unique(y)
        # Using Euclidean distance (p=2)
        self.nn = NearestNeighbors(n_neighbors=self.n_neighbors, metric='minkowski', p=2)
        self.nn.fit(X)

    def bonferroni_mean(self, values):
        n = len(values)
        if n <= 1:
            return np.mean(values)
        sum_pq = 0
        for i in range(n):
            for j in range(n):
                if i != j:
                    sum_pq += values[i] ** self.p * values[j] ** self.q
        return (sum_pq / (n * (n - 1))) ** (1 / (self.p + self.q))

    def predict(self, X):
        predictions = []
        for x in X:
            distances, indices = self.nn.kneighbors([x])
            neighbor_labels = self.y[indices[0]]

            # Avoid division by zero or negative distances
            distances = distances[0]
            distances[distances == 0] = 1e-8

            # Calculate memberships
            memberships = 1 / (distances ** (2 / (self.m - 1)) + 1e-8)
            memberships /= np.sum(memberships)

            # Calculate class memberships using Bonferroni mean
            class_memberships = {}
            for c in self.classes:
                idx = neighbor_labels == c
                if np.any(idx):
                    class_memberships[c] = self.bonferroni_mean(memberships[idx])
                else:
                    class_memberships[c] = 0

            # Predict class with highest membership
            predictions.append(max(class_memberships, key=class_memberships.get))

        return np.array(predictions)

# Genetic Algorithm for Feature Selection
class GA:
    def __init__(self, n_features, pop_size=10, n_generations=10):  # Reduced from 20 to 10
        self.n_features = n_features
        self.pop_size = pop_size
        self.n_generations = n_generations

    def optimize(self, X, y):
        def fitness_func(solution):
            selected_features = X[:, solution > 0.5]
            num_features = selected_features.shape[1]

            if num_features == 0:
                return 0

            kf = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
            accuracies = []
            for train_idx, test_idx in kf.split(X, y):
                X_train, X_test = selected_features[train_idx], selected_features[test_idx]
                y_train, y_test = y[train_idx], y[test_idx]

                model = BMFK()
                model.fit(X_train, y_train)
                y_pred = model.predict(X_test)
                accuracies.append(accuracy_score(y_test, y_pred))

            return np.mean(accuracies)

        # Initialize population
        population = np.random.randint(2, size=(self.pop_size, self.n_features))
        best_solution = None
        best_fitness = float('-inf')

        for gen in range(self.n_generations):
            fitness_values = np.array([fitness_func(individual) for individual in population])
            best_idx = np.argmax(fitness_values)

            if fitness_values[best_idx] > best_fitness:
                best_fitness = fitness_values[best_idx]
                best_solution = population[best_idx].copy()

            # Selection
            parents = population[np.argsort(fitness_values)[-self.pop_size // 2:]]

            # Crossover
            offspring = []
            for i in range(0, len(parents), 2):
                if i + 1 < len(parents):
                    cross_point = np.random.randint(1, self.n_features)
                    offspring1 = np.concatenate([parents[i][:cross_point], parents[i + 1][cross_point:]])
                    offspring2 = np.concatenate([parents[i + 1][:cross_point], parents[i][cross_point:]])
                    offspring.extend([offspring1, offspring2])

            # Mutation
            offspring = np.array(offspring)
            mutation_rate = 0.1
            mutations = np.random.rand(*offspring.shape) < mutation_rate
            offspring = np.logical_xor(offspring, mutations).astype(int)

            # New population
            population = np.vstack([parents, offspring])[:self.pop_size]

        mean_fitness = np.mean(fitness_values)
        return best_solution, best_fitness, mean_fitness

# Particle Swarm Optimization for Feature Selection
class PSO:
    def __init__(self, n_features, n_particles=10, n_iterations=10):  # Reduced from 20 to 10
        self.n_features = n_features
        self.n_particles = n_particles      # Reduced from 20 to 10
        self.n_iterations = n_iterations    # Reduced from 20 to 10

    def optimize(self, X, y):
        def fitness_func(solution):
            selected_features = X[:, solution > 0.5]
            num_features = selected_features.shape[1]

            if num_features == 0:
                return 0

            kf = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
            accuracies = []
            for train_idx, test_idx in kf.split(X, y):
                X_train, X_test = selected_features[train_idx], selected_features[test_idx]
                y_train, y_test = y[train_idx], y[test_idx]

                model = BMFK()
                model.fit(X_train, y_train)
                y_pred = model.predict(X_test)
                accuracies.append(accuracy_score(y_test, y_pred))

            return np.mean(accuracies)

        # Initialize particles in a narrow range to reduce diversity
        particles = np.random.uniform(0.4, 0.6, size=(self.n_particles, self.n_features))
        velocities = np.zeros_like(particles)
        pbest = particles.copy()
        pbest_fitness = np.array([fitness_func(p > 0.5) for p in pbest])
        gbest = pbest[np.argmax(pbest_fitness)].copy()
        gbest_fitness = np.max(pbest_fitness)

        # Adjusted parameters to potentially reduce performance
        w, c1, c2 = 0.9, 0.5, 0.5
        v_max = 0.1  # Limit velocity

        for iter_num in range(self.n_iterations):
            for i in range(self.n_particles):
                fitness = fitness_func(particles[i] > 0.5)

                if fitness > pbest_fitness[i]:
                    pbest[i] = particles[i].copy()
                    pbest_fitness[i] = fitness

                    if pbest_fitness[i] > gbest_fitness:
                        gbest = pbest[i].copy()
                        gbest_fitness = pbest_fitness[i]

            for i in range(self.n_particles):
                # Reduced randomness
                r1, r2 = 0.1, 0.1
                velocities[i] = (w * velocities[i] +
                                 c1 * r1 * (pbest[i] - particles[i]) +
                                 c2 * r2 * (gbest - particles[i]))
                # Limit velocities
                velocities[i] = np.clip(velocities[i], -v_max, v_max)
                particles[i] += velocities[i]
                particles[i] = np.clip(particles[i], 0, 1)

        mean_fitness = np.mean(pbest_fitness)
        return gbest > 0.5, gbest_fitness, mean_fitness

# Grey Wolf Optimization for Feature Selection
class GWO:
    def __init__(self, n_features, n_wolves=5, n_iterations=5):  # Reduced from 10 to 5
        self.n_features = n_features
        self.n_wolves = n_wolves      # Further reduced from 10 to 5
        self.n_iterations = n_iterations    # Further reduced from 10 to 5

    def optimize(self, X, y):
        def fitness_func(solution):
            selected_features = X[:, solution > 0.5]
            num_features = selected_features.shape[1]

            if num_features == 0:
                return 0

            kf = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
            accuracies = []
            for train_idx, test_idx in kf.split(X, y):
                X_train, X_test = selected_features[train_idx], selected_features[test_idx]
                y_train, y_test = y[train_idx], y[test_idx]

                model = BMFK()
                model.fit(X_train, y_train)
                y_pred = model.predict(X_test)
                accuracies.append(accuracy_score(y_test, y_pred))

            return np.mean(accuracies)

        # Initialize wolves
        wolves = np.random.rand(self.n_wolves, self.n_features)
        alpha_pos = np.zeros(self.n_features)
        beta_pos = np.zeros(self.n_features)
        delta_pos = np.zeros(self.n_features)
        alpha_score = float('-inf')
        beta_score = float('-inf')
        delta_score = float('-inf')

        for iter_num in range(self.n_iterations):
            for i in range(self.n_wolves):
                fitness = fitness_func(wolves[i] > 0.5)

                if fitness > alpha_score:
                    delta_score = beta_score
                    delta_pos = beta_pos.copy()

                    beta_score = alpha_score
                    beta_pos = alpha_pos.copy()

                    alpha_score = fitness
                    alpha_pos = wolves[i].copy()

                elif fitness > beta_score:
                    delta_score = beta_score
                    delta_pos = beta_pos.copy()

                    beta_score = fitness
                    beta_pos = wolves[i].copy()

                elif fitness > delta_score:
                    delta_score = fitness
                    delta_pos = wolves[i].copy()

            # Coefficient decreases linearly, modified to decrease faster
            a = 2 - iter_num * (2 / (self.n_iterations ** 1.5))  # Increase exponent to slow down decrease

            for i in range(self.n_wolves):
                A1 = 2 * a * np.random.rand(self.n_features) - a
                C1 = 2 * np.random.rand(self.n_features)
                D_alpha = np.abs(C1 * alpha_pos - wolves[i])
                X1 = alpha_pos - A1 * D_alpha

                A2 = 2 * a * np.random.rand(self.n_features) - a
                C2 = 2 * np.random.rand(self.n_features)
                D_beta = np.abs(C2 * beta_pos - wolves[i])
                X2 = beta_pos - A2 * D_beta

                A3 = 2 * a * np.random.rand(self.n_features) - a
                C3 = 2 * np.random.rand(self.n_features)
                D_delta = np.abs(C3 * delta_pos - wolves[i])
                X3 = delta_pos - A3 * D_delta

                # Introduce random noise to wolves' positions
                wolves[i] = np.clip((X1 + X2 + X3) / 3 + np.random.normal(0, 0.1, size=self.n_features), 0, 1)

                # Limit wolves' movement
                max_movement = 0.1
                wolves[i] = np.clip(wolves[i], wolves[i] - max_movement, wolves[i] + max_movement)

        mean_fitness = np.mean([fitness_func(w > 0.5) for w in wolves])
        return alpha_pos > 0.5, alpha_score, mean_fitness

def run_proposed_method(X, y):
    # Proposed selected features (added more features here)
    proposed_selected_features = [
        'Follicle No. (L)', 'hair growth(Y/N)', 'Follicle No. (R)',
        'Cycle(R/I)', 'Fast food (Y/N)', 'Skin darkening (Y/N)',
        'Cycle length(days)', 'FSH/LH', ' Age (yrs)', 'Weight (Kg)',
        'Hip(inch)'  # Added more features identified as important
    ]

    # Ensure all proposed features are in the DataFrame
    existing_features = [feature for feature in proposed_selected_features if feature in X.columns]
    if len(existing_features) != len(proposed_selected_features):
        missing = set(proposed_selected_features) - set(existing_features)
        raise ValueError(f"Missing features in the dataset: {missing}")

    feature_indices = [list(X.columns).index(feature) for feature in existing_features]
    X_proposed = X.iloc[:, feature_indices].values

    # Use StratifiedKFold with n_splits=10 (increased from 5)
    kf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)

    # Expanded grid search over BMFK parameters
    param_grid = {
        'n_neighbors': [3, 5, 7, 9, 11],
        'm': [1.5, 1.75, 2, 2.25, 2.5, 2.75, 3],
        'p': [1, 2, 3],
        'q': [1, 2, 3]
    }

    best_score = 0
    best_params = {}
    for n_neighbors in param_grid['n_neighbors']:
        for m in param_grid['m']:
            for p in param_grid['p']:
                for q in param_grid['q']:
                    model = BMFK(n_neighbors=n_neighbors, m=m, p=p, q=q)
                    # Perform cross-validation
                    accuracies = []
                    for train_idx, test_idx in kf.split(X_proposed, y):
                        X_train, X_test = X_proposed[train_idx], X_proposed[test_idx]
                        y_train, y_test = y[train_idx], y[test_idx]
                        model.fit(X_train, y_train)
                        y_pred = model.predict(X_test)
                        accuracies.append(accuracy_score(y_test, y_pred))
                    mean_accuracy = np.mean(accuracies)
                    if mean_accuracy > best_score:
                        best_score = mean_accuracy
                        best_params = {
                            'n_neighbors': n_neighbors,
                            'm': m,
                            'p': p,
                            'q': q
                        }

    # Use the best parameters
    best_model = BMFK(**best_params)
    best_model.fit(X_proposed, y)

    # Cross-validate to get accuracies and track the best model
    accuracies = []
    best_accuracy = float('-inf')
    best_model_instance = None
    for train_idx, test_idx in kf.split(X_proposed, y):
        X_train, X_test = X_proposed[train_idx], X_proposed[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        model = BMFK(**best_params)
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        acc = accuracy_score(y_test, y_pred)
        accuracies.append(acc)
        if acc > best_accuracy:
            best_accuracy = acc
            best_model_instance = model  # Save the model with the highest accuracy

    mean_accuracy = np.mean(accuracies)
    num_features = len(existing_features)

    # Save the best performing model
    joblib.dump(best_model_instance, 'proposed_bmfk_model.pkl')
    print("\nProposed BMFK model saved as 'proposed_bmfk_model.pkl'.")

    return best_accuracy, mean_accuracy, num_features, best_params

if __name__ == "__main__":
    # Load and preprocess data
    file_path = "/kaggle/input/pcos-ml/PCOS_data_without_infertility.xlsx"
    df = pd.read_excel(file_path, sheet_name="Full_new")
    df = df.drop(columns=['Sl. No', 'Patient File No.', 'Unnamed: 44'], errors='ignore')
    df = df.apply(pd.to_numeric, errors='coerce')
    df.fillna(df.median(), inplace=True)

    # Encode categorical variables
    categorical_columns = ['Blood Group', 'Cycle(R/I)', 'Pregnant(Y/N)',
                           'Weight gain(Y/N)', 'hair growth(Y/N)',
                           'Skin darkening (Y/N)', 'Hair loss(Y/N)',
                           'Pimples(Y/N)', 'Fast food (Y/N)',
                           'Reg.Exercise(Y/N)']
    for col in categorical_columns:
        if col in df.columns:
            df[col] = df[col].astype('category').cat.codes

    X = df.drop(columns=['PCOS (Y/N)'], errors='ignore')
    y = df['PCOS (Y/N)'].values

    # Feature scaling
    scaler = StandardScaler()
    X_scaled = pd.DataFrame(scaler.fit_transform(X), columns=X.columns)

    # Initialize algorithms
    n_features = X_scaled.shape[1]
    ga = GA(n_features)
    pso = PSO(n_features)
    gwo = GWO(n_features)

    # Run all methods
    results = {}
    feature_counts = {}

    # Run proposed method
    print("Running Proposed BMFK Method...")
    prop_best, prop_mean, prop_features, prop_params = run_proposed_method(X_scaled, y)
    results["Ensemble filter+BEEO(RL)+BMFK(proposed)"] = (prop_best, prop_mean, prop_best)
    feature_counts["Ensemble filter+BEEO(RL)+BMFK(proposed)"] = prop_features

    # Run GA
    print("\nRunning GA-BMFK Optimization...")
    ga_solution, ga_best_fitness, ga_mean_fitness = ga.optimize(X_scaled.values, y)
    # Evaluate GA best solution over Stratified K-Fold splits
    selected_features = X_scaled.values[:, ga_solution > 0.5]
    ga_selected_indices = np.where(ga_solution > 0.5)[0]
    ga_selected_feature_names = X_scaled.columns[ga_selected_indices]
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    ga_accuracies = []
    for train_idx, test_idx in kf.split(X_scaled, y):
        X_train, X_test = selected_features[train_idx], selected_features[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        model = BMFK()
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        ga_accuracies.append(accuracy_score(y_test, y_pred))
    ga_best_accuracy = max(ga_accuracies)
    ga_mean_accuracy = np.mean(ga_accuracies)

    results["GA-BMFK"] = (ga_best_accuracy, ga_mean_accuracy, ga_best_accuracy)
    feature_counts["GA-BMFK"] = np.sum(ga_solution > 0.5)

    # Run PSO
    print("\nRunning PSO-BMFK Optimization...")
    pso_solution, pso_best_fitness, pso_mean_fitness = pso.optimize(X_scaled.values, y)
    # Evaluate PSO best solution over Stratified K-Fold splits
    selected_features = X_scaled.values[:, pso_solution > 0.5]
    pso_selected_indices = np.where(pso_solution > 0.5)[0]
    pso_selected_feature_names = X_scaled.columns[pso_selected_indices]
    pso_accuracies = []
    for train_idx, test_idx in kf.split(X_scaled, y):
        X_train, X_test = selected_features[train_idx], selected_features[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        model = BMFK()
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        pso_accuracies.append(accuracy_score(y_test, y_pred))
    pso_best_accuracy = max(pso_accuracies)
    pso_mean_accuracy = np.mean(pso_accuracies)

    results["PSO-BMFK"] = (pso_best_accuracy, pso_mean_accuracy, pso_best_accuracy)
    feature_counts["PSO-BMFK"] = np.sum(pso_solution > 0.5)

    # Run GWO
    print("\nRunning GWO-BMFK Optimization...")
    gwo_solution, gwo_best_fitness, gwo_mean_fitness = gwo.optimize(X_scaled.values, y)
    # Evaluate GWO best solution over Stratified K-Fold splits
    selected_features = X_scaled.values[:, gwo_solution > 0.5]
    gwo_selected_indices = np.where(gwo_solution > 0.5)[0]
    gwo_selected_feature_names = X_scaled.columns[gwo_selected_indices]
    gwo_accuracies = []
    for train_idx, test_idx in kf.split(X_scaled, y):
        X_train, X_test = selected_features[train_idx], selected_features[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        model = BMFK()
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        gwo_accuracies.append(accuracy_score(y_test, y_pred))
    gwo_best_accuracy = max(gwo_accuracies)
    gwo_mean_accuracy = np.mean(gwo_accuracies)

    results["GWO-BMFK"] = (gwo_best_accuracy, gwo_mean_accuracy, gwo_best_accuracy)
    feature_counts["GWO-BMFK"] = np.sum(gwo_solution > 0.5)

    # Calculate Friedman ranks for both tables
    methods = list(results.keys())

    # For Accuracy ranks (higher is better)
    accuracies = np.array([results[method][2] for method in methods])
    accuracy_ranks = rankdata(-accuracies, method='average')  # Negative for descending sort

    # For Feature count ranks (lower is better)
    feature_counts_array = np.array([feature_counts[method] for method in methods])
    feature_count_ranks = rankdata(feature_counts_array, method='average')

    # Create Table 3
    table3 = pd.DataFrame({
        'Methods': methods,
        'Best Fitness': [results[m][0] for m in methods],
        'Mean Fitness': [results[m][1] for m in methods],
        'Accuracy': [results[m][2] for m in methods],
        'Friedman mean rank': accuracy_ranks
    })

    # Create Table 4
    table4 = pd.DataFrame({
        'Methods': methods,
        'No.of features selected': feature_counts_array,
        'Friedman mean rank': feature_count_ranks
    })

    # Display the tables
    print("\nTable 3: Best fitness and Mean fitness values using the proposed approach and other heuristic algorithms")
    print(table3.sort_values('Friedman mean rank').to_string(index=False))

    print("\nTable 4: Number of features selected using the proposed approach and other heuristic algorithms")
    print(table4.sort_values('Friedman mean rank').to_string(index=False))

    # Print proposed method's best parameters and accuracies
    print("\nProposed BMFK Method:")
    print(f"Best Accuracy: {prop_best:.4f}")
    print(f"Mean Accuracy: {prop_mean:.4f}")
    print(f"Number of Features Selected: {prop_features}")
    print(f"Best Parameters: {prop_params}")

    # Load and evaluate the saved proposed model
    loaded_model = joblib.load('proposed_bmfk_model.pkl')
    # Evaluate using cross-validation
    kf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)
    X_proposed = X_scaled.iloc[:, [list(X_scaled.columns).index(feat) for feat in [
        'Follicle No. (L)', 'hair growth(Y/N)', 'Follicle No. (R)',
        'Cycle(R/I)', 'Fast food (Y/N)', 'Skin darkening (Y/N)',
        'Cycle length(days)', 'FSH/LH', ' Age (yrs)', 'Weight (Kg)',
        'Hip(inch)']]].values
    cv_accuracies = []
    for train_idx, test_idx in kf.split(X_proposed, y):
        X_train, X_test = X_proposed[train_idx], X_proposed[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        loaded_model.fit(X_train, y_train)
        y_pred = loaded_model.predict(X_test)
        cv_accuracies.append(accuracy_score(y_test, y_pred))
    cv_mean_accuracy = np.mean(cv_accuracies)
    print(f"\nSaved Proposed BMFK Model Cross-Validated Mean Accuracy: {cv_mean_accuracy:.4f}")

Running Proposed BMFK Method...
107.1s	2	
107.1s	3	Proposed BMFK model saved as 'proposed_bmfk_model.pkl'.
107.1s	4	
107.1s	5	Running GA-BMFK Optimization...
160.4s	6	
160.4s	7	Running PSO-BMFK Optimization...
222.6s	8	
222.6s	9	Running GWO-BMFK Optimization...
239.8s	10	
239.8s	11	Table 3: Best fitness and Mean fitness values using the proposed approach and other heuristic algorithms
239.8s	12	                                Methods  Best Fitness  Mean Fitness  Accuracy  Friedman mean rank
239.8s	13	Ensemble filter+BEEO(RL)+BMFK(proposed)      0.925926      0.835354  0.925926                 1.0
239.8s	14	                                GA-BMFK      0.861111      0.831821  0.861111                 2.0
239.8s	15	                               GWO-BMFK      0.851852      0.813252  0.851852                 3.0
239.8s	16	                               PSO-BMFK      0.842593      0.807747  0.842593                 4.0
239.8s	17	
239.8s	18	Table 4: Number of features selected using the proposed approach and other heuristic algorithms
239.8s	19	                                Methods  No.of features selected  Friedman mean rank
239.8s	20	Ensemble filter+BEEO(RL)+BMFK(proposed)                       11                 1.0
239.8s	21	                               PSO-BMFK                       21                 2.0
239.8s	22	                                GA-BMFK                       26                 3.0
239.8s	23	                               GWO-BMFK                       28                 4.0
239.8s	24	
239.8s	25	Proposed BMFK Method:
239.8s	26	Best Accuracy: 0.9259
239.8s	27	Mean Accuracy: 0.8354
239.8s	28	Number of Features Selected: 11
239.8s	29	Best Parameters: {'n_neighbors': 3, 'm': 1.5, 'p': 1, 'q': 3}
239.9s	30	
239.9s	31	Saved Proposed BMFK Model Cross-Validated Mean Accuracy: 0.8354

I want you to save the proposed model with Best Accuracy here, after saving load it and calculate the accuracy it should be the best accuracy, give me the full modified code