if __name__ == "__main__":
    rng = np.random.default_rng(0)

    # Generate data with 2 true informative dimensions embedded
    # in 10 dimensions via a random linear map plus noise.
    n, true_dim, ambient_dim = 300, 2, 10
    Z_true = rng.normal(size=(n, true_dim))
    labels = (Z_true[:, 0] > 0).astype(int)
    random_map = rng.normal(size=(true_dim, ambient_dim))
    X = Z_true @ random_map + 0.1 * rng.normal(size=(n, ambient_dim))

    model = pca_fit(X, n_components=5)
    print("Explained variance ratio by component:")
    print(np.round(model["explained_variance_ratio"], 3))
    cumulative = np.cumsum(model["explained_variance_ratio"])
    print("Cumulative:", np.round(cumulative, 3))

    Z_2d = pca_transform(X, model)[:, :2]
    X_reconstructed = pca_inverse_transform(
        np.hstack([Z_2d, np.zeros((n, 3))]), model
    )
    reconstruction_error = np.mean((X - X_reconstructed) ** 2)
    print(f"Reconstruction error (2 components): {reconstruction_error:.4f}")

    class_separation = np.abs(
        Z_2d[labels == 1, 0].mean() - Z_2d[labels == 0, 0].mean()
    )
    print(f"Class separation along PC1: {class_separation:.3f}")
