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
    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]
        best_model.fit(X_train, y_train)
        y_pred = best_model.predict(X_test)
        accuracies.append(accuracy_score(y_test, y_pred))
    best_accuracy = max(accuracies)
    mean_accuracy = np.mean(accuracies)
    num_features = len(existing_features)

    # Save the trained best model
    joblib.dump(best_model, '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
    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}")

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

save the proposed one with high accuracy that is the one with the best accuracy achieved, here it is saving Saved Proposed BMFK Model Cross-Validated Mean Accuracy: 0.8354
but the best accuracy is Best Accuracy: 0.9259. I want that model to be saved, save that one
give me the full modified code