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

            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=10, n_iterations=10):  # Reduced from 20 to 10
        self.n_features = n_features
        self.n_wolves = n_wolves      # 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 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
    ]

    feature_indices = [list(X.columns).index(feature) for feature in proposed_selected_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
    model = BMFK(**best_params)
    # Re-evaluate on 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))
    best_accuracy = max(accuracies)
    mean_accuracy = np.mean(accuracies)
    num_features = len(proposed_selected_features)
    return best_accuracy, mean_accuracy, num_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)

    # 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)'])
    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 = 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]
    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_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_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))

GA& GWO mean fitness is higher than the proposed one, please increase the proposed mean fitness, give me the complete code