"""From-scratch principal component analysis via the singular value decomposition.

PCA finds the orthogonal directions of maximum variance in the (centered) data.
The clean way to compute it is the SVD of the centered matrix X = U S V^T:

    * the rows of V^T (columns of V) are the principal directions / loadings,
    * the scores (coordinates in the new basis) are X V = U S,
    * the variance explained by component k is S_k^2 / (n-1).

This is numerically preferable to eigen-decomposing the covariance matrix and is
exactly what scikit-learn's PCA does; the notebook validates the match.

Factor analysis (used via scikit-learn in the notebook) is the probabilistic
cousin: it models the data as a few COMMON latent factors plus per-variable
UNIQUE noise, X = L f + e, separating shared covariation from idiosyncratic
variance -- whereas PCA lumps all variance together into orthogonal directions.
"""
import numpy as np


def pca(X, standardize=False):
    """Principal component analysis by SVD.

    Parameters
    ----------
    X : (n_samples, n_features) data.
    standardize : if True, scale each feature to unit variance first
                  (correlation-based PCA); else covariance-based.

    Returns dict with scores, loadings (components as rows), explained-variance
    ratio, singular values, and the center/scale used.
    """
    X = np.asarray(X, float)
    mean = X.mean(0)
    Xc = X - mean
    scale = np.ones(X.shape[1])
    if standardize:
        scale = Xc.std(0, ddof=1); scale[scale == 0] = 1.0
        Xc = Xc / scale
    U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
    var = S ** 2 / (len(X) - 1)
    return {
        "scores": U * S,            # coordinates of each sample on the components
        "loadings": Vt,             # principal directions (one per row)
        "explained_variance": var,
        "explained_variance_ratio": var / var.sum(),
        "singular_values": S,
        "mean": mean, "scale": scale,
    }


def reconstruct(res, n_components):
    """Rank-`n_components` reconstruction of the (centered, scaled) data."""
    Z = res["scores"][:, :n_components] @ res["loadings"][:n_components]
    return Z * res["scale"] + res["mean"]