import numpy as np
from sklearn import datasets
from sklearn.decomposition import PCA
def PCA_eig(X,k, center=True, scale=False):
  n,p = X.shape
  ones = np.ones([n,1])
  h = ((1/n) * np.matmul(ones, ones.transpose())) if center  else np.zeros([n,n])
  H = np.eye(n) - h
  X_center =  np.matmul(H, X)
  covariance = 1/(n-1) * np.matmul(X_center.transpose(), X_center)
  scaling =  np.sqrt(1/np.diag(covariance)) if scale else np.ones(p)
  scaled_covariance = np.matmul(np.diag(scaling), covariance)
  w,v = np.linalg.eig(scaled_covariance)
  components = v[:, :k]
  explained_variance = w[:k]
  return { 'X':X, 'k':k, 'components':components,
    'explained_variance':explained_variance }