使用PSO-NN解决双螺旋分类TSNP问题

2025-03-14
Course Assignment
核心思想

粒子群算法PSO优化神经网络权重和偏置参数，取代传统的梯度下降法。
神经网络参数矩阵展平，作为PSO算法中粒子的速度参数进行优化。
Presentation PPT

报告

代码

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <random>
#include <limits>
#include <algorithm>
#include <fstream>

using namespace std;

// data parameter
const int num_spiral_rounds = 5;
// PSO-Parameter
const double w = 0.2;    
const double c1 = 1;
const double c2 = 2;
const double LR = 1.2;
const int num_particles = 1000;
const int max_iterations = 10000; 
// NN-Parameter
const int M_hidden_nodes = 32;


struct DataPoint {
    vector<double> input; // Position:(x, y)
    double label;         // Lable: 0 / 1
};

struct Particle {
    vector<double> position;  
    vector<double> velocity;  
    vector<double> best_position; 
    double best_fitness;     
};

struct GlobalBest {
    vector<double> position; 
    double fitness;         
};

double Tanh_activation_function(double x) {
    return tanh(x);
}
double Sigmiod_activation_function(double x) {
    return 1.0 / (1.0 + exp(-x));
}
double ReLU_activation_function(double x) {
    return max(0.0, x);
}

double fitness_function(const vector<double>& position, const vector<DataPoint>& data, int hidden_nodes) {
    //==================== 2-layer NN dim ====================
    int input_dim = 2; 
    int hidden_dim = hidden_nodes; 
    int output_dim = 1; 

    //==================== Extract (W,b) from particles' weight sequence ====================
    vector<vector<double>> input_to_hidden(input_dim, vector<double>(hidden_dim));
    vector<double> hidden_to_output(hidden_dim);
    vector<double> hidden_bias(hidden_dim);
    double output_bias;

    int idx = 0;
    // 2*M: W (Input->Hidden)
    for (int i = 0; i < input_dim; i++) {
        for (int j = 0; j < hidden_dim; j++) {
            input_to_hidden[i][j] = position[idx++];
        }
    }
    // 1*M: W (Hidden->Output)
    for (int j = 0; j < hidden_dim; j++) {
        hidden_to_output[j] = position[idx++];
    }
    // 1*M: b (Hidden)
    for (int j = 0; j < hidden_dim; j++) {
        hidden_bias[j] = position[idx++];
    }
    // 1: b(Output)
    output_bias = position[idx];

    //==================== Forward-Propagate: Build up the NN / then get error for PSO to optimize position(W,b) ====================
    double error = 0.0;
    int correct_predictions = 0;
    //for every data
    for (const auto& sample : data) {
        const vector<double>& input = sample.input;
        double true_output = sample.label;
        //Forward-Propagate: HiddenLayer
        vector<double> hidden_output(hidden_dim, 0.0);
        for (int j = 0; j < hidden_dim; j++) {
            for (int i = 0; i < input_dim; i++) {
                //y = w1*x1 + w2*x2 + b
                hidden_output[j] += input[i] * input_to_hidden[i][j];
                //y = w1*(x1^2+x2^2) + w2*(x1*x2) + b 
                /*
                if(i == 1)
                    hidden_output[j] += input_to_hidden[i-1][j] * (pow(input[i-1],2) + pow(input[i],2)) 
                                      + input_to_hidden[i][j] * (input[i-1] * input[i-1]);
                else
                    continue;
                */
                //y = w1*sinx1 + w2*cosx2 + b 
                /*
                if(i == 1)
                    hidden_output[j] += sin(input[i-1]) * input_to_hidden[i-1][j] + cos(input[i]) * input_to_hidden[i][j];
                else
                    continue;
                */   
            }
            hidden_output[j] += hidden_bias[j];
            hidden_output[j] = Tanh_activation_function(hidden_output[j]);
            //hidden_output[j] = ReLU_activation_function(hidden_output[j]);
            //hidden_output[j] = Sigmiod_activation_function(hidden_output[j]);
        }
        //Forward-Propagate: OutputLayer
        double output = 0.0;
        for (int j = 0; j < hidden_dim; j++) {
            output += hidden_output[j] * hidden_to_output[j];
        }
        output += output_bias;
        output = Sigmiod_activation_function(output); 
        //Error Calculation
        error += pow(output - true_output, 2) / data.size();
        // Accuracy calculation
        int predicted_label = (output >= 0.5) ? 1 : 0;
        if (predicted_label == true_output) {
            correct_predictions++;
        }
    }
    double accuracy = static_cast<double>(correct_predictions) / data.size();
    
    //double E = error + (1-accuracy);
    double E = error;
    return E;
}


GlobalBest pso(const vector<DataPoint>& data, int hidden_nodes, int num_particles, int max_iterations) {

    int dim = 2 * hidden_nodes + hidden_nodes + hidden_nodes + 1;  // Particle.Position dim: 2*M + 1*M + 1*M + 1;
    
    //==================== Initialization of Particles ====================
    random_device rd;
    mt19937 gen(rd());
    uniform_real_distribution<> dist(0.0, 1.0);

    vector<Particle> particles(num_particles);
    GlobalBest global_best;
    global_best.fitness = numeric_limits<double>::infinity();

    for (auto& particle : particles) {
        particle.position.resize(dim);
        particle.velocity.resize(dim);
        particle.best_position.resize(dim);
        for (int i = 0; i < dim; i++) {
            particle.position[i] = dist(gen);
            particle.velocity[i] = dist(gen) * 0.1;
        }
        particle.best_position = particle.position;
        particle.best_fitness = fitness_function(particle.position, data, hidden_nodes);

        if (particle.best_fitness < global_best.fitness) {
            global_best.position = particle.best_position;
            global_best.fitness = particle.best_fitness;
        }
    }

    // PSO迭代
    ofstream outfile("result/fitness_results_ARM5.txt");
    for (int iter = 0; iter < max_iterations; iter++) {
        for (auto& particle : particles) {

            double fitness = fitness_function(particle.position, data, hidden_nodes);

            if (fitness < particle.best_fitness) {
                particle.best_fitness = fitness;
                particle.best_position = particle.position;
            }

            if (fitness < global_best.fitness) {
                global_best.fitness = fitness;
                global_best.position = particle.position;
            }
            for (int i = 0; i < dim; i++) {
                double r1 = dist(gen), r2 = dist(gen);
                particle.velocity[i] = w * particle.velocity[i] +
                                       c1 * r1 * (particle.best_position[i] - particle.position[i]) +
                                       c2 * r2 * (global_best.position[i] - particle.position[i]);
                particle.position[i] += particle.velocity[i] * LR;
            }
        }
        cout << "Iteration " << iter + 1 << ": Best Fitness(ErrorRate) = " << global_best.fitness << endl;
        outfile << "Iteration " << iter + 1 << ": Best Fitness(ErrorRate) = " << global_best.fitness << endl;
    }

    outfile.close(); 
    return global_best;
}

// Two-Nested-Spirals
vector<DataPoint> generate_spiral_data(int num_points, int num_turns) {
    vector<DataPoint> data;
    double pi = 3.141592653589793;
    for (int i = 0; i < num_points; i++) {
        double t = i * pi * num_turns / num_points;
        data.push_back(\{\{t * cos(t), t * sin(t)\}, 0.0\}); 
        data.push_back(\{\{-t * cos(t), -t * sin(t)\}, 1.0\}); 
    }
    return data;
}

void normalize_data(vector<DataPoint>& data, double max_value) {
    for (auto& point : data) {
        point.input[0] /= max_value;
        point.input[1] /= max_value; 
    }
}

pair<vector<DataPoint>, vector<DataPoint>> split_data_for_validation(const vector<DataPoint>& data, double validation_ratio) {
    int validation_size = static_cast<int>(data.size() * validation_ratio);

    vector<DataPoint> shuffled_data = data;
    random_device rd;
    mt19937 gen(rd());
    shuffle(shuffled_data.begin(), shuffled_data.end(), gen);

    vector<DataPoint> validation_data(shuffled_data.begin(), shuffled_data.begin() + validation_size);
    vector<DataPoint> training_data(shuffled_data.begin() + validation_size, shuffled_data.end());

    return {training_data, validation_data};
}



double test_model(const vector<double>& best_position, const vector<DataPoint>& test_data, int hidden_nodes) {

    int input_dim = 2;
    int hidden_dim = hidden_nodes;

    //==================== Extract (W,b) from the trained weights ====================
    vector<vector<double>> input_to_hidden(input_dim, vector<double>(hidden_dim));
    vector<double> hidden_to_output(hidden_dim);
    vector<double> hidden_bias(hidden_dim);
    double output_bias;

    int idx = 0;
    // 2*M: W (Input->Hidden)
    for (int i = 0; i < input_dim; i++) {
        for (int j = 0; j < hidden_dim; j++) {
            input_to_hidden[i][j] = best_position[idx++];
        }
    }
    // 1*M: W (Hidden->Output)
    for (int j = 0; j < hidden_dim; j++) {
        hidden_to_output[j] = best_position[idx++];
    }
    // 1*M: b (Hidden)
    for (int j = 0; j < hidden_dim; j++) {
        hidden_bias[j] = best_position[idx++];
    }
    // 1: b(Output)
    output_bias = best_position[idx];

    //==================== Forward-Propagate for Testing ====================
    int correct_predictions = 0;

    for (const auto& sample : test_data) {
        const vector<double>& input = sample.input;
        double true_output = sample.label;

        // Forward-Propagate: Hidden Layer 
        vector<double> hidden_output(hidden_dim, 0.0);
        for (int j = 0; j < hidden_dim; j++) {
            for (int i = 0; i < input_dim; i++) {
                //y = w1*x1 + w2*x2 + b
                hidden_output[j] += input[i] * input_to_hidden[i][j];
                //y = w1*(x1^2+x2^2) + w2*(x1*x2) + b 
                /*
                if(i == 1)
                    hidden_output[j] += input_to_hidden[i-1][j] * (pow(input[i-1],2) + pow(input[i],2)) 
                                      + input_to_hidden[i][j] * (input[i-1] * input[i-1]);
                else
                    continue;
                */
                //y = w1*sinx1 + w2*cosx2 + b 
                /*
                if(i == 1)
                    hidden_output[j] += sin(input[i-1]) * input_to_hidden[i-1][j] + cos(input[i]) * input_to_hidden[i][j];
                else
                    continue;
                */
            }
            hidden_output[j] += hidden_bias[j];
            hidden_output[j] = Tanh_activation_function(hidden_output[j]);
            //hidden_output[j] = ReLU_activation_function(hidden_output[j]);
            //hidden_output[j] = Sigmiod_activation_function(hidden_output[j]);
        }

        // Forward-Propagate: Output Layer
        double output = 0.0;
        for (int j = 0; j < hidden_dim; j++) {
            output += hidden_output[j] * hidden_to_output[j];
        }
        output += output_bias;
        output = Sigmiod_activation_function(output); // Sigmoid activation

        // Classify and Count Correct Predictions
        int predicted_label = (output >= 0.5) ? 1 : 0; // Sigmoid threshold at 0.5
        if (predicted_label == true_output) {
            correct_predictions++;
        }
    }

    //==================== Calculate Accuracy ====================
    double accuracy = static_cast<double>(correct_predictions) / test_data.size();
    return accuracy;
}





int main() {
    
    int num_points = 100 * num_spiral_rounds;  
    vector<DataPoint> data = generate_spiral_data(num_points,num_spiral_rounds);
    normalize_data(data,30.0);

    double validation_ratio = 0.2;
    pair<vector<DataPoint>, vector<DataPoint>> split_result = split_data_for_validation(data, validation_ratio);
    vector<DataPoint> train_data = split_result.first;
    vector<DataPoint> validation_data = split_result.second;

    GlobalBest best_solution = pso(train_data, M_hidden_nodes, num_particles, max_iterations);
    cout << "Best Fitness(Minimum-ErrorRate): " << best_solution.fitness << endl;

    double validation_accuracy = test_model(best_solution.position, validation_data, M_hidden_nodes);
    cout << "Model Accuracy on Validation Data: " << validation_accuracy * 100 << "%" << endl;

    return 0;