function main()
    % 清除環境
    clear; close all; clc;
    
    fprintf('=== 基於粒子羣優化算法的特徵選擇SVM分類 ===\n');
    
    % 加載數據
    [data, labels] = loadSampleData();
    
    % 數據預處理
    [normalized_data, label_encoded] = preprocessData(data, labels);
    
    % 設置PSO參數
    pso_params = setPSOParameters(size(normalized_data, 2));
    
    % 執行PSO特徵選擇
    fprintf('開始PSO特徵選擇優化...\n');
    [best_features, best_fitness, convergence_curve] = ...
        PSO_FeatureSelection(normalized_data, label_encoded, pso_params);
    
    % 顯示結果
    displayResults(best_features, best_fitness, convergence_curve, ...
                  normalized_data, label_encoded, pso_params);
end

function [data, labels] = loadSampleData()
    % 加載示例數據 - 這裏使用UCI葡萄酒數據集作為示例
    % 實際使用時可以替換為自己的數據
    try
        data = readtable('wine.data', 'FileType', 'text');
        labels = data{:, 1};
        data = data{:, 2:end};
        fprintf('成功加載葡萄酒數據集: %d個樣本, %d個特徵\n', size(data, 1), size(data, 2));
    catch
        % 如果文件不存在，生成模擬數據
        fprintf('使用模擬數據...\n');
        [data, labels] = generateSimulatedData();
    end
end

function [data, labels] = generateSimulatedData()
    % 生成模擬數據用於演示
    rng(42); % 設置隨機種子保證可重複性
    
    n_samples = 200;
    n_features = 30;
    n_informative = 10; % 只有10個特徵是有用的
    
    % 生成特徵數據
    data = randn(n_samples, n_features);
    
    % 只有前n_informative個特徵與標籤相關
    informative_weights = randn(n_informative, 1);
    linear_combination = data(:, 1:n_informative) * informative_weights;
    
    % 生成標籤 (3類)
    probabilities = 1 ./ (1 + exp(-linear_combination));
    labels = discretize(probabilities, [0, 0.33, 0.66, 1]);
    
    % 添加噪聲特徵
    data(:, n_informative+1:end) = data(:, n_informative+1:end) + 0.5 * randn(n_samples, n_features - n_informative);
    
    fprintf('生成模擬數據: %d個樣本, %d個特徵, %d個類別\n', ...
            n_samples, n_features, length(unique(labels)));
end

function [normalized_data, label_encoded] = preprocessData(data, labels)
    % 數據預處理
    
    % 標準化特徵
    normalized_data = zscore(data);
    
    % 確保標籤從1開始連續編號
    unique_labels = unique(labels);
    label_encoded = zeros(size(labels));
    for i = 1:length(unique_labels)
        label_encoded(labels == unique_labels(i)) = i;
    end
    
    fprintf('數據預處理完成: 樣本數=%d, 特徵數=%d, 類別數=%d\n', ...
            size(normalized_data, 1), size(normalized_data, 2), length(unique_labels));
end

function params = setPSOParameters(n_features)
    % 設置PSO參數
    params.n_particles = 30;           % 粒子數量
    params.max_iter = 50;              % 最大迭代次數
    params.w = 0.9;                    % 慣性權重
    params.w_damp = 0.99;              % 慣性權重衰減係數
    params.c1 = 2.0;                   % 個體學習因子
    params.c2 = 2.0;                   % 社會學習因子
    params.n_features = n_features;    % 總特徵數
    params.velocity_max = 0.5;         % 最大速度
    params.cost_function = @fitnessFunction; % 適應度函數
    
    fprintf('PSO參數設置: %d個粒子, %d次迭代, %d個特徵\n', ...
            params.n_particles, params.max_iter, params.n_features);
end

function [best_features, best_fitness, convergence_curve] = ...
    PSO_FeatureSelection(data, labels, params)
    % PSO特徵選擇主函數
    
    % 初始化粒子羣
    particles = initializeParticles(params);
    
    % 初始化個體最優和全局最優
    personal_best.position = particles.position;
    personal_best.cost = inf(params.n_particles, 1);
    personal_best.feature_count = zeros(params.n_particles, 1);
    
    global_best.cost = inf;
    global_best.position = [];
    global_best.feature_count = 0;
    
    convergence_curve = zeros(params.max_iter, 1);
    
    % PSO主循環
    for iter = 1:params.max_iter
        for i = 1:params.n_particles
            % 計算當前粒子的適應度
            current_position = particles.position(i, :);
            [current_cost, feature_count] = params.cost_function(...
                current_position, data, labels);
            
            % 更新個體最優
            if current_cost < personal_best.cost(i)
                personal_best.position(i, :) = current_position;
                personal_best.cost(i) = current_cost;
                personal_best.feature_count(i) = feature_count;
            end
            
            % 更新全局最優
            if current_cost < global_best.cost
                global_best.position = current_position;
                global_best.cost = current_cost;
                global_best.feature_count = feature_count;
            end
        end
        
        % 記錄收斂曲線
        convergence_curve(iter) = global_best.cost;
        
        % 更新粒子速度和位置
        particles = updateParticles(particles, personal_best, global_best, params, iter);
        
        % 顯示進度
        if mod(iter, 10) == 0
            fprintf('迭代 %d/%d: 最佳適應度 = %.4f, 選擇特徵數 = %d\n', ...
                    iter, params.max_iter, global_best.cost, global_best.feature_count);
        end
    end
    
    % 獲取最終結果
    best_features = global_best.position > 0.5;
    best_fitness = global_best.cost;
end

function particles = initializeParticles(params)
    % 初始化粒子羣
    particles.position = rand(params.n_particles, params.n_features);
    particles.velocity = zeros(params.n_particles, params.n_features);
    
    % 隨機初始化部分特徵被選中
    for i = 1:params.n_particles
        % 每個粒子隨機選擇約50%的特徵
        threshold = rand(1) * 0.3 + 0.2; % 閾值在0.2-0.5之間
        particles.position(i, :) = particles.position(i, :) > threshold;
    end
end

function particles = updateParticles(particles, personal_best, global_best, params, iter)
    % 更新粒子位置和速度
    
    % 更新慣性權重
    w = params.w * (params.w_damp ^ (iter-1));
    
    for i = 1:params.n_particles
        % 更新速度
        r1 = rand(1, params.n_features);
        r2 = rand(1, params.n_features);
        
        cognitive_component = params.c1 * r1 .* (personal_best.position(i, :) - particles.position(i, :));
        social_component = params.c2 * r2 .* (global_best.position - particles.position(i, :));
        
        particles.velocity(i, :) = w * particles.velocity(i, :) + cognitive_component + social_component;
        
        % 限制速度範圍
        particles.velocity(i, :) = max(particles.velocity(i, :), -params.velocity_max);
        particles.velocity(i, :) = min(particles.velocity(i, :), params.velocity_max);
        
        % 更新位置
        particles.position(i, :) = particles.position(i, :) + particles.velocity(i, :);
        
        % 使用sigmoid函數將位置映射到[0,1]範圍
        sigmoid_pos = 1 ./ (1 + exp(-particles.position(i, :)));
        
        % 轉換為二進制位置（特徵選擇）
        particles.position(i, :) = sigmoid_pos > 0.5;
    end
end

function [fitness, feature_count] = fitnessFunction(feature_mask, data, labels)
    % 適應度函數：基於SVM分類性能
    
    % 獲取選擇的特徵索引
    selected_features = find(feature_mask);
    feature_count = length(selected_features);
    
    % 如果沒有選擇任何特徵，給予懲罰
    if feature_count == 0
        fitness = inf;
        return;
    end
    
    % 如果選擇特徵太多，給予懲罰（鼓勵稀疏性）
    total_features = size(data, 2);
    if feature_count > total_features * 0.8
        penalty = (feature_count / total_features) * 2;
    else
        penalty = 1;
    end
    
    % 提取選擇的特徵
    X_selected = data(:, selected_features);
    
    % 使用5折交叉驗證評估SVM性能
    cv = cvpartition(labels, 'KFold', 5);
    cv_accuracy = zeros(cv.NumTestSets, 1);
    
    for fold = 1:cv.NumTestSets
        train_idx = cv.training(fold);
        test_idx = cv.test(fold);
        
        % 訓練SVM模型
        svm_model = fitcsvm(X_selected(train_idx, :), labels(train_idx), ...
                           'KernelFunction', 'rbf', 'BoxConstraint', 1, ...
                           'KernelScale', 'auto', 'Standardize', true);
        
        % 預測
        predictions = predict(svm_model, X_selected(test_idx, :));
        
        % 計算準確率
        cv_accuracy(fold) = sum(predictions == labels(test_idx)) / length(predictions);
    end
    
    % 平均準確率
    mean_accuracy = mean(cv_accuracy);
    
    % 適應度 = 1/準確率 + 特徵數量懲罰
    % 目標是最小化適應度（即最大化準確率同時最小化特徵數量）
    alpha = 0.01; % 特徵數量懲罰係數
    fitness = (1 - mean_accuracy) + alpha * (feature_count / total_features);
    
    % 應用懲罰項
    fitness = fitness * penalty;
end

function displayResults(best_features, best_fitness, convergence_curve, ...
                       data, labels, params)
    % 顯示和分析結果
    
    fprintf('\n=== 特徵選擇結果 ===\n');
    
    % 選擇的特徵信息
    selected_indices = find(best_features);
    n_selected = length(selected_indices);
    n_total = params.n_features;
    
    fprintf('總特徵數: %d\n', n_total);
    fprintf('選擇的特徵數: %d (%.1f%%)\n', n_selected, n_selected/n_total*100);
    fprintf('最佳適應度: %.4f\n', best_fitness);
    fprintf('選擇的特徵索引: %s\n', mat2str(selected_indices));
    
    % 性能比較
    comparePerformance(data, labels, best_features);
    
    % 繪製收斂曲線
    figure('Position', [100, 100, 1200, 800]);
    
    subplot(2, 3, 1);
    plot(convergence_curve, 'b-', 'LineWidth', 2);
    xlabel('迭代次數');
    ylabel('適應度值');
    title('PSO收斂曲線');
    grid on;
    
    % 繪製特徵選擇情況
    subplot(2, 3, 2);
    feature_importance = calculateFeatureImportance(data, labels);
    bar(feature_importance);
    xlabel('特徵索引');
    ylabel('重要性得分');
    title('特徵重要性排序');
    grid on;
    
    % 標記選擇的特徵
    hold on;
    plot(selected_indices, feature_importance(selected_indices), 'ro', ...
         'MarkerSize', 8, 'LineWidth', 2);
    legend('所有特徵', '選擇的特徵', 'Location', 'best');
    
    % 繪製特徵數量變化
    subplot(2, 3, 3);
    hist_data = randi([1, n_total], 1000, 1);
    histogram(hist_data, 20, 'FaceColor', [0.8, 0.8, 0.8]);
    hold on;
    xline(n_selected, 'r-', 'LineWidth', 3);
    xlabel('特徵數量');
    ylabel('頻次');
    title('選擇的特徵數量');
    legend('隨機選擇', 'PSO選擇', 'Location', 'best');
    grid on;
    
    % 繪製SVM分類邊界（使用前兩個主要特徵）
    if n_selected >= 2
        subplot(2, 3, [4, 5, 6]);
        plotSVMBoundary(data(:, selected_indices(1:min(2, n_selected))), labels);
        title('SVM分類邊界 (使用選擇的特徵)');
    end
    
    % 顯示詳細的性能指標
    showDetailedMetrics(data, labels, best_features);
end

function comparePerformance(data, labels, best_features)
    % 比較使用所有特徵和選擇特徵的性能
    
    fprintf('\n=== 性能比較 ===\n');
    
    % 使用所有特徵
    [all_acc, all_time] = evaluateSVM(data, labels);
    
    % 使用選擇的特徵
    selected_data = data(:, best_features);
    [selected_acc, selected_time] = evaluateSVM(selected_data, labels);
    
    fprintf('所有特徵 - 準確率: %.2f%%, 訓練時間: %.4f秒\n', all_acc*100, all_time);
    fprintf('選擇特徵 - 準確率: %.2f%%, 訓練時間: %.4f秒\n', selected_acc*100, selected_time);
    fprintf('特徵減少: %.1f%%, 時間減少: %.1f%%\n', ...
            (1 - nnz(best_features)/size(data,2))*100, ...
            (1 - selected_time/all_time)*100);
    
    if selected_acc > all_acc
        fprintf('✅ 特徵選擇提高了分類性能!\n');
    else
        fprintf('⚠️  特徵選擇略微降低了分類性能，但提高了效率\n');
    end
end

function [accuracy, training_time] = evaluateSVM(data, labels)
    % 評估SVM性能
    cv = cvpartition(labels, 'KFold', 5);
    accuracies = zeros(cv.NumTestSets, 1);
    times = zeros(cv.NumTestSets, 1);
    
    for fold = 1:cv.NumTestSets
        train_idx = cv.training(fold);
        test_idx = cv.test(fold);
        
        tic;
        svm_model = fitcsvm(data(train_idx, :), labels(train_idx), ...
                           'KernelFunction', 'rbf', 'Standardize', true);
        times(fold) = toc;
        
        predictions = predict(svm_model, data(test_idx, :));
        accuracies(fold) = sum(predictions == labels(test_idx)) / length(predictions);
    end
    
    accuracy = mean(accuracies);
    training_time = mean(times);
end

function importance = calculateFeatureImportance(data, labels)
    % 計算特徵重要性（基於隨機森林）
    try
        tree_bagger = TreeBagger(50, data, labels, 'Method', 'classification', ...
                                'OOBPredictorImportance', 'on');
        importance = tree_bagger.OOBPermutedPredictorDeltaError;
    catch
        % 如果TreeBagger不可用，使用簡單的方差方法
        importance = std(data)';
    end
end

function plotSVMBoundary(data, labels)
    % 繪製SVM分類邊界（僅適用於2D）
    if size(data, 2) < 2
        return;
    end
    
    % 訓練SVM模型
    svm_model = fitcsvm(data(:, 1:2), labels, 'KernelFunction', 'rbf');
    
    % 創建網格
    h = 0.02;
    x1 = min(data(:,1)):h:max(data(:,1));
    x2 = min(data(:,2)):h:max(data(:,2));
    [X1, X2] = meshgrid(x1, x2);
    
    % 預測網格點
    predictions = predict(svm_model, [X1(:), X2(:)]);
    Z = reshape(predictions, size(X1));
    
    % 繪製決策邊界
    contourf(X1, X2, Z, 'Alpha', 0.3);
    hold on;
    
    % 繪製數據點
    unique_labels = unique(labels);
    colors = ['r', 'g', 'b', 'c', 'm', 'y'];
    for i = 1:length(unique_labels)
        idx = labels == unique_labels(i);
        scatter(data(idx, 1), data(idx, 2), 50, colors(i), 'filled');
    end
    
    xlabel('特徵 1');
    ylabel('特徵 2');
    legend('決策區域', '類別1', '類別2', '類別3', 'Location', 'best');
    grid on;
end

function showDetailedMetrics(data, labels, best_features)
    % 顯示詳細的性能指標
    fprintf('\n=== 詳細性能指標 ===\n');
    
    selected_data = data(:, best_features);
    cv = cvpartition(labels, 'KFold', 5);
    
    accuracies = zeros(cv.NumTestSets, 1);
    precisions = zeros(cv.NumTestSets, 1);
    recalls = zeros(cv.NumTestSets, 1);
    f1_scores = zeros(cv.NumTestSets, 1);
    
    for fold = 1:cv.NumTestSets
        train_idx = cv.training(fold);
        test_idx = cv.test(fold);
        
        svm_model = fitcsvm(selected_data(train_idx, :), labels(train_idx));
        predictions = predict(svm_model, selected_data(test_idx, :));
        
        % 計算多分類指標（使用宏平均）
        cm = confusionmat(labels(test_idx), predictions);
        
        % 計算每個類別的精度和召回率
        n_classes = size(cm, 1);
        class_precision = zeros(n_classes, 1);
        class_recall = zeros(n_classes, 1);
        
        for i = 1:n_classes
            class_precision(i) = cm(i,i) / sum(cm(:, i));
            class_recall(i) = cm(i,i) / sum(cm(i, :));
        end
        
        % 宏平均
        precisions(fold) = mean(class_precision, 'omitnan');
        recalls(fold) = mean(class_recall, 'omitnan');
        accuracies(fold) = sum(diag(cm)) / sum(cm(:));
        f1_scores(fold) = 2 * (precisions(fold) * recalls(fold)) / ...
                          (precisions(fold) + recalls(fold));
    end
    
    fprintf('平均準確率: %.2f%% ± %.2f\n', mean(accuracies)*100, std(accuracies)*100);
    fprintf('平均精確率: %.2f%% ± %.2f\n', mean(precisions)*100, std(precisions)*100);
    fprintf('平均召回率: %.2f%% ± %.2f\n', mean(recalls)*100, std(recalls)*100);
    fprintf('平均F1分數: %.2f%% ± %.2f\n', mean(f1_scores)*100, std(f1_scores)*100);
end

function quick_test()
    % 快速測試函數
    clear; close all; clc;
    
    fprintf('快速測試PSO-SVM特徵選擇...\n');
    
    % 生成小型測試數據
    rng(123);
    [data, labels] = generateSimulatedData();
    [normalized_data, label_encoded] = preprocessData(data, labels);
    
    % 簡化參數用於快速測試
    params.n_particles = 20;
    params.max_iter = 20;
    params.w = 0.9;
    params.w_damp = 0.95;
    params.c1 = 2.0;
    params.c2 = 2.0;
    params.n_features = size(normalized_data, 2);
    params.velocity_max = 0.5;
    params.cost_function = @fitnessFunction;
    
    % 執行PSO
    [best_features, best_fitness, convergence_curve] = ...
        PSO_FeatureSelection(normalized_data, label_encoded, params);
    
    % 顯示簡要結果
    selected_count = nnz(best_features);
    fprintf('\n快速測試完成!\n');
    fprintf('選擇特徵: %d/%d\n', selected_count, params.n_features);
    fprintf('最佳適應度: %.4f\n', best_fitness);
    
    % 繪製收斂曲線
    figure;
    plot(convergence_curve, 'b-', 'LineWidth', 2);
    title('PSO收斂曲線 (快速測試)');
    xlabel('迭代次數');
    ylabel('適應度值');
    grid on;
end

% 運行完整版本
main();

% 運行快速測試
quick_test();

% 替換數據加載部分
function [data, labels] = loadYourData()
    % 加載您的數據
    % data: n_samples × n_features 矩陣
    % labels: n_samples × 1 向量
    data = your_feature_matrix;
    labels = your_labels;
end