main.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_openml

def compute_covariance_matrix(data):
    return np.cov(data, rowvar=False)

def calculate_eigenvalues_and_eigenvectors(cov_matrix):
    eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
    return eigenvalues, eigenvectors

def sort_eigenvectors(eigenvalues, eigenvectors):
    sorted_indices = np.argsort(eigenvalues)[::-1]
    sorted_eigenvalues = eigenvalues[sorted_indices]
    sorted_eigenvectors = eigenvectors[:, sorted_indices]
    return sorted_eigenvalues, sorted_eigenvectors

def select_top_k_eigenvectors(eigenvectors, k):
    return eigenvectors[:, :k]

def project_onto_principal_components(data, principal_components):
    return np.dot(data, principal_components)

class PCA:
    def __init__(self, n_components=None):
        self.n_components = n_components
        self.eigenvalues = None
        self.eigenvectors = None
        self.principal_components = None
    
    def fit(self, data):
        cov_matrix = compute_covariance_matrix(data)
        self.eigenvalues, self.eigenvectors = calculate_eigenvalues_and_eigenvectors(cov_matrix)
        sorted_eigenvalues, sorted_eigenvectors = sort_eigenvectors(self.eigenvalues, self.eigenvectors)
        self.eigenvalues = sorted_eigenvalues
        self.eigenvectors = sorted_eigenvectors
        if self.n_components:
            self.principal_components = select_top_k_eigenvectors(sorted_eigenvectors, self.n_components)
        else:
            self.principal_components = sorted_eigenvectors
    
    def fit_transform(self, data):
        self.fit(data)
        return project_onto_principal_components(data, self.principal_components)
    
    def inverse_transform(self, data):
        return np.dot(data, self.principal_components.T)
    
    def explained_variance_ratio_(self):
        return self.eigenvalues / np.sum(self.eigenvalues)
    
    def plot_variance_explained(self):
        plt.figure(figsize=(8, 4))
        plt.plot(np.cumsum(self.explained_variance_ratio_()))
        plt.xlabel('Number of Components')
        plt.ylabel('Variance')
        # plt.title('Variance Explained by Principal Components')
        plt.grid(True)
        plt.show()
    
    def transform(self, data):
        return project_onto_principal_components(data, self.principal_components)

def load_and_normalize_mnist():
    mnist = fetch_openml('mnist_784', version=1)
    images = mnist.data / 255.
    labels = mnist.target.astype(int)
    return images, labels


def plot_images(original, reconstructed, n_images=10):
    if isinstance(original, pd.DataFrame):
        original = original.values

    print(f"Type of original: {type(original)}, Shape of original: {original.shape}")
    print(f"Type of reconstructed: {type(reconstructed)}, Shape of reconstructed: {reconstructed.shape}")

    plt.figure(figsize=(20, 4))
    for i in range(n_images):
        ax = plt.subplot(2, n_images, i + 1)
        plt.imshow(original[i].reshape(28, 28), cmap='gray')
        ax.axis('off')

        ax = plt.subplot(2, n_images, i + 1 + n_images)
        plt.imshow(reconstructed[i].reshape(28, 28), cmap='gray')
        ax.axis('off')
    # plt.show()

def compute_mse(original, reconstructed):
    return np.mean((original - reconstructed) ** 2)

def main():
    images, labels = load_and_normalize_mnist()
    
    pca = PCA()
    pca.fit(images)
    pca.plot_variance_explained()

    # Determine the optimal number of components (e.g., 95% variance explained)
    cumulative_variance = np.cumsum(pca.explained_variance_ratio_())
    optimal_components = np.argmax(cumulative_variance >= 0.99) + 1

    print(f"Optimal number of components to explain 95% variance: {optimal_components}")

    pca = PCA(n_components=optimal_components)
    principal_components = pca.fit_transform(images)
    reconstructed_images = pca.inverse_transform(principal_components)
    
    print(f"Number of components: {optimal_components}")
    plot_images(images, reconstructed_images)
    # plt.suptitle(f"Original and Reconstructed Images with {optimal_components} Components")
    plt.show()


    # mse_values = []
    # components_list = [10, 50, 100, 200]
    # for n_components in components_list:
    #     pca = PCA(n_components=n_components)
    #     principal_components = pca.fit_transform(images)
    #     reconstructed_images = pca.inverse_transform(principal_components)
        
    #     mse = compute_mse(images, reconstructed_images)
    #     mse_values.append(mse)

    #     print(f"Number of components: {n_components}, MSE: {mse}")
    #     plot_images(images, reconstructed_images)
    #     plt.suptitle(f"Original and Reconstructed Images with {n_components} Components")
    #     plt.show()

    # print(f"Number of components: {optimal_components}, MSE: {mse}")
    # plot_images(images, reconstructed_images)
    # plt.suptitle(f"Original and Reconstructed Images with {optimal_components} Components")
    # plt.show()

    mse_values = []
    components_list = range(10, 300, 10)
    for n_components in components_list:
        pca = PCA(n_components=n_components)
        principal_components = pca.fit_transform(images)
        reconstructed_images = pca.inverse_transform(principal_components)
        
        mse = compute_mse(images, reconstructed_images)
        mse_values.append(mse)

        # print(f"Number of components: {n_components}, MSE: {mse}")
        # plot_images(images, reconstructed_images)
        # plt.suptitle(f"Original and Reconstructed Images with {n_components} Components")
        # plt.show()

    plt.figure(figsize=(8, 4))
    plt.plot(components_list, mse_values, marker='o')
    plt.xlabel('Number of Components')
    plt.ylabel('Mean Squared Error')
    # plt.title('MSE vs. Number of Components')
    plt.grid(True)
    plt.show()

if __name__ == "__main__":
    main()