Principal Component Analysis

#include <Eigen/Dense>
#include <iostream>
#include <random>

using namespace std;
using namespace Eigen;

int main() {
    // Set up random number generators
    std::mt19937 gen(1);  // Seed for reproducibility
    std::normal_distribution<> normal(0.0, 1.0);
    std::uniform_real_distribution<> uniform(0.0, 1.0);

    // Generate a 3x3 random matrix for transformation
    MatrixXd rand_transform(3, 3);
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            rand_transform(i, j) = uniform(gen);

    // Generate 3x150 matrix of normal random numbers
    MatrixXd rand_norm(3, 150);
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 150; ++j)
            rand_norm(i, j) = normal(gen);

    // Create the dataset: shape (150, 3)
    MatrixXd X = (rand_transform * rand_norm).transpose();

    // (Optional) Print first few rows
    cout << "First 5 data points:\n" << X.topRows(5) << endl;

    // ... continue with next steps
}

// Calculate the mean for each column
RowVectorXd X_mean = X.colwise().mean();

// Center the data (subtract the mean)
MatrixXd X_centered = X.rowwise() - X_mean;

// Calculate the standard deviation for each column
RowVectorXd X_std(X.cols());
for (int j = 0; j < X.cols(); ++j) {
    double sum_sq = 0.0;
    for (int i = 0; i < X.rows(); ++i) {
        sum_sq += pow(X_centered(i, j), 2);
    }
    X_std(j) = sqrt(sum_sq / (X.rows() - 1));
}

// Standardize the data (divide by std deviation)
MatrixXd X_standardized = X_centered;
for (int j = 0; j < X.cols(); ++j) {
    for (int i = 0; i < X.rows(); ++i) {
        X_standardized(i, j) /= X_std(j);
    }
}

// (Optional) Print first few standardized rows
cout << "First 5 standardized data points:\n" << X_standardized.topRows(5) << endl;

// Compute the covariance matrix (3x3)
MatrixXd cov_matrix = (X_standardized.adjoint() * X_standardized) / double(X_standardized.rows() - 1);

cout << "Covariance matrix:\n" << cov_matrix << endl;

#include <Eigen/Eigenvalues>

// Compute eigenvalues and eigenvectors
SelfAdjointEigenSolver<MatrixXd> eig_solver(cov_matrix);
VectorXd eigenvalues = eig_solver.eigenvalues();
MatrixXd eigenvectors = eig_solver.eigenvectors(); // Columns are eigenvectors

cout << "Eigenvalues:\n" << eigenvalues << endl;
cout << "Eigenvectors (columns):\n" << eigenvectors << endl;

#include <vector>
#include <algorithm>

// Pair each eigenvalue with its eigenvector
vector<pair<double, VectorXd>> eigen_pairs;
for (int i = 0; i < eigenvalues.size(); ++i) {
    eigen_pairs.emplace_back(eigenvalues(i), eigenvectors.col(i));
}

// Sort pairs in descending order of eigenvalue
std::sort(eigen_pairs.begin(), eigen_pairs.end(),
          [](const auto& a, const auto& b) { return a.first > b.first; });

// (Optional) Print sorted eigenvalues
cout << "Sorted eigenvalues:\n";
for (const auto& p : eigen_pairs) {
    cout << p.first << " ";
}
cout << endl;

// Form the projection matrix using the top 2 eigenvectors
MatrixXd W(X.cols(), 2);
W.col(0) = eigen_pairs[0].second;
W.col(1) = eigen_pairs[1].second;

// Project the data onto the new 2D space
MatrixXd X_pca = X_standardized * W;

// (Optional) Print first few projected points
cout << "First 5 projected data points:\n" << X_pca.topRows(5) << endl;

#include <matplotlibcpp.h>
namespace plt = matplotlibcpp;

// Prepare data for plotting
vector<double> x_vals, y_vals;
for (int i = 0; i < X_pca.rows(); ++i) {
    x_vals.push_back(X_pca(i, 0));
    y_vals.push_back(X_pca(i, 1));
}

// Create scatter plot
plt::scatter(x_vals, y_vals, 20.0);
plt::xlabel("Principal Component 1");
plt::ylabel("Principal Component 2");
plt::title("Scatter Plot of Transformed Dataset Using PCA");
plt::grid(true);
plt::save("static/images/plot.png"); // Save the plot as an image