import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import KFold
from sklearn.metrics import accuracy_score
from sklearn.neighbors import NearestNeighbors
from scipy.stats import rankdata
import warnings
warnings.filterwarnings('ignore')

# BMFK Classifier (your existing 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)
        self.nn = NearestNeighbors(n_neighbors=self.n_neighbors, metric='minkowski', p=2)
        self.nn.fit(X)
        # print(f"BMFK: Fitted NearestNeighbors with {X.shape[1]} features.")

    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):
        # print(f"BMFK: Predicting with {X.shape[1]} features.")
        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

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

            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

            predictions.append(max(class_memberships, key=class_memberships.get))

        return np.array(predictions)

# Optimized heuristic algorithms
class GA:
    def __init__(self, n_features, pop_size=20, n_generations=20):
        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]
            # print(f"GA Fitness Function: Selected {num_features} features.")

            if num_features == 0:
                return 0

            kf = KFold(n_splits=3, shuffle=True, random_state=42)
            accuracies = []
            for train_idx, test_idx in kf.split(X):
                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)

        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):
            # print(f"GA Generation {gen + 1}/{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()
                # print(f"GA New Best Fitness: {best_fitness}")

            # 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)

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

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

class PSO:
    def __init__(self, n_features, n_particles=20, n_iterations=20):
        self.n_features = n_features
        self.n_particles = n_particles
        self.n_iterations = n_iterations

    def optimize(self, X, y):
        def fitness_func(solution):
            selected_features = X[:, solution > 0.5]
            num_features = selected_features.shape[1]
            # print(f"PSO Fitness Function: Selected {num_features} features.")

            if num_features == 0:
                return 0

            kf = KFold(n_splits=3, shuffle=True, random_state=42)
            accuracies = []
            for train_idx, test_idx in kf.split(X):
                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)

        particles = np.random.rand(self.n_particles, self.n_features)
        velocities = np.zeros_like(particles)
        pbest = particles.copy()
        pbest_fitness = np.array([fitness_func(p) for p in pbest])
        gbest = pbest[np.argmax(pbest_fitness)].copy()
        gbest_fitness = np.max(pbest_fitness)

        w, c1, c2 = 0.7, 1.5, 1.5

        for iter_num in range(self.n_iterations):
            # print(f"PSO Iteration {iter_num + 1}/{self.n_iterations}")
            for i in range(self.n_particles):
                fitness = fitness_func(particles[i])

                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]
                        # print(f"PSO New Best Fitness: {gbest_fitness}")

            for i in range(self.n_particles):
                r1, r2 = np.random.rand(2)
                velocities[i] = (w * velocities[i] +
                                 c1 * r1 * (pbest[i] - particles[i]) +
                                 c2 * r2 * (gbest - particles[i]))
                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

class GWO:
    def __init__(self, n_features, n_wolves=20, n_iterations=20):
        self.n_features = n_features
        self.n_wolves = n_wolves
        self.n_iterations = n_iterations

    def optimize(self, X, y):
        def fitness_func(solution):
            selected_features = X[:, solution > 0.5]
            num_features = selected_features.shape[1]
            # print(f"GWO Fitness Function: Selected {num_features} features.")

            if num_features == 0:
                return 0

            kf = KFold(n_splits=3, shuffle=True, random_state=42)
            accuracies = []
            for train_idx, test_idx in kf.split(X):
                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)

        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):
            # print(f"GWO Iteration {iter_num + 1}/{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()
                    # print(f"GWO New Alpha Fitness: {alpha_score}")

                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()

            a = 2 - iter_num * (2 / self.n_iterations)

            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

                wolves[i] = np.clip((X1 + X2 + X3) / 3, 0, 1)

        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 = ['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']

    feature_indices = [list(X.columns).index(feature) for feature in proposed_selected_features]
    X_proposed = X.iloc[:, feature_indices].values
    # print(f"Proposed Method: Selected {X_proposed.shape[1]} features.")

    kf = KFold(n_splits=3, shuffle=True, random_state=42)
    accuracies = []
    for train_idx, test_idx in kf.split(X_proposed):
        X_train, X_test = X_proposed[train_idx], X_proposed[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.max(accuracies), np.mean(accuracies), len(proposed_selected_features)

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'])
    df = df.apply(pd.to_numeric, errors='coerce')
    df.fillna(df.median(), inplace=True)

    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)'])
    y = df['PCOS (Y/N)'].values

    scaler = StandardScaler()
    X_scaled = pd.DataFrame(scaler.fit_transform(X), columns=X.columns)
    # print(f"X_scaled shape: {X_scaled.shape}")

    # 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 = run_proposed_method(X_scaled, y)
    results["Ensemble filter+BEEO(RL)+BMFK(proposed)"] = (prop_mean, 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 KFold splits
    selected_features = X_scaled.values[:, ga_solution > 0.5]
    kf = KFold(n_splits=3, shuffle=True, random_state=42)
    ga_accuracies = []
    for train_idx, test_idx in kf.split(X_scaled):
        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_fitness, ga_mean_fitness, 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 KFold splits
    selected_features = X_scaled.values[:, pso_solution > 0.5]
    kf = KFold(n_splits=3, shuffle=True, random_state=42)
    pso_accuracies = []
    for train_idx, test_idx in kf.split(X_scaled):
        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_fitness, pso_mean_fitness, 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 KFold splits
    selected_features = X_scaled.values[:, gwo_solution > 0.5]
    kf = KFold(n_splits=3, shuffle=True, random_state=42)
    gwo_accuracies = []
    for train_idx, test_idx in kf.split(X_scaled):
        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_fitness, gwo_mean_fitness, 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))

Adjust the parameter and try to increase our proposed or reduce the pso values