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 friedmanchisquare
import warnings
warnings.filterwarnings('ignore')

# BMFK Classifier (your existing implementation with debugging)
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])
            neighbors = self.X[indices[0]]
            neighbor_labels = self.y[indices[0]]
            
            memberships = 1 / (distances[0] ** (2 / (self.m - 1)) + 1e-8)
            memberships /= np.sum(memberships)
            
            class_memberships = {}
            for c in self.classes:
                class_memberships[c] = self.bonferroni_mean(memberships[neighbor_labels == c])
            
            predictions.append(max(class_memberships, key=class_memberships.get))
        
        return np.array(predictions)

# Optimized heuristic algorithms with debugging
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]
                print(f"GA New Best Fitness: {best_fitness}")
            
            parents = population[np.argsort(fitness_values)[-self.pop_size//2:]]
            offspring = []
            for i in range(0, len(parents), 2):
                if i+1 < len(parents):
                    cross_point = np.random.randint(1, self.n_features)
                    offspring.append(np.concatenate([parents[i][:cross_point], parents[i+1][cross_point:]]))
                    offspring.append(np.concatenate([parents[i+1][:cross_point], parents[i][cross_point:]]))
            
            population = np.vstack([parents, offspring])[:self.pop_size]
        
        return best_solution, best_fitness, np.mean(fitness_values)

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)]
        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]
                    pbest_fitness[i] = fitness
                    
                    if pbest_fitness[i] > gbest_fitness:
                        gbest = pbest[i]
                        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)
        
        return gbest, gbest_fitness, np.mean(pbest_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')
        
        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])
                
                if fitness > alpha_score:
                    delta_pos = beta_pos
                    beta_pos = alpha_pos
                    alpha_pos = wolves[i].copy()
                    alpha_score = fitness
                    print(f"GWO New Alpha Fitness: {alpha_score}")
            
            a = 2 - iter_num * (2 / self.n_iterations)
            
            for i in range(self.n_wolves):
                A1, A2, A3 = a * (2 * np.random.rand(3) - 1)
                C1, C2, C3 = 2 * np.random.rand(3)
                
                X1 = alpha_pos - A1 * np.abs(C1 * alpha_pos - wolves[i])
                X2 = beta_pos - A2 * np.abs(C2 * beta_pos - wolves[i])
                X3 = delta_pos - A3 * np.abs(C3 * delta_pos - wolves[i])
                
                wolves[i] = np.clip((X1 + X2 + X3) / 3, 0, 1)
        
        return alpha_pos, alpha_score, np.mean([fitness_func(w) for w in wolves])

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_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, ga_mean = ga.optimize(X_scaled.values, y)
    results["GA-BMFK"] = (ga_best, ga_mean, ga_best)
    feature_counts["GA-BMFK"] = np.sum(ga_solution > 0.5)
    
    # Run PSO
    print("\nRunning PSO-BMFK Optimization...")
    pso_solution, pso_best, pso_mean = pso.optimize(X_scaled.values, y)
    results["PSO-BMFK"] = (pso_best, pso_mean, pso_best)
    feature_counts["PSO-BMFK"] = np.sum(pso_solution > 0.5)
    
    # Run GWO
    print("\nRunning GWO-BMFK Optimization...")
    gwo_solution, gwo_best, gwo_mean = gwo.optimize(X_scaled.values, y)
    results["GWO-BMFK"] = (gwo_best, gwo_mean, gwo_best)
    feature_counts["GWO-BMFK"] = np.sum(gwo_solution > 0.5)

    # Calculate Friedman ranks for both tables
    methods = list(results.keys())
    
    # For Accuracy ranks
    accuracies = np.array([results[method][2] for method in methods])
    accuracy_ranks = pd.Series(accuracies).rank(ascending=False)
    
    # For Feature count ranks (lower is better)
    feature_counts_array = np.array([feature_counts[method] for method in methods])
    feature_count_ranks = pd.Series(feature_counts_array).rank()
    
    # 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
    })
    
    # Optional: Perform Friedman Test if you have multiple related samples
    # Note: Friedman test requires multiple observations per method. 
    # Since we have only one observation per method, it's not applicable here.
    # If you have multiple runs, you can uncomment and modify the following lines:

    # Example (assuming you have multiple accuracy scores for each method):
    # accuracies_list = [
    #     [run1_accuracy_ga, run2_accuracy_ga, ...],
    #     [run1_accuracy_pso, run2_accuracy_pso, ...],
    #     ...
    # ]
    # stat, p = friedmanchisquare(*accuracies_list)
    # print(f"Friedman test statistic: {stat}, p-value: {p}")

    # Display the tables
    print("\nTable 3: Best fitness and Mean fitness values using the proposed approach and other heuristic algorithms")
    print(table3.to_string(index=False))
    
    print("\nTable 4: Number of features selected using the proposed approach and other heuristic algorithms")
    print(table4.to_string(index=False))

Table 3: Best fitness and Mean fitness values using the proposed approach and other heuristic algorithms
698.8s	8717	                                Methods  Best Fitness  Mean Fitness  Accuracy  Friedman mean rank
698.8s	8718	Ensemble filter+BEEO(RL)+BMFK(proposed)      0.767956      0.715245  0.767956                 1.0
698.8s	8719	                                GA-BMFK      0.730121      0.730121  0.730121                 4.0
698.8s	8720	                               PSO-BMFK      0.748527      0.747060  0.748527                 3.0
698.8s	8721	                               GWO-BMFK      0.750440      0.741907  0.750440                 2.0
698.8s	8722	
698.8s	8723	Table 4: Number of features selected using the proposed approach and other heuristic algorithms
698.8s	8724	                                Methods  No.of features selected  Friedman mean rank
698.8s	8725	Ensemble filter+BEEO(RL)+BMFK(proposed)                        8                 1.0
698.8s	8726	                                GA-BMFK                       23                 3.0
698.8s	8727	                               PSO-BMFK                       22                 2.0
698.8s	8728	                               GWO-BMFK                       31                 4.0

implement the calculation of Best Fitness  Mean Fitness  Accuracy  Friedman mean rank correctly and give me the complete code