docs: def PCA(Y, components):
  	"""
  	run PCA, retrieving the first (components) principle components
  	return [s0, eig, w0]
  	s0: factors
  	w0: weights
  	"""
  
  	N,D = Y.shape
  	sv = linalg.svd(Y, full_matrices=0);
  	[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
  	v = s0.std(axis=0)
  	s0 /= v;
  	w0 *= v;
  	return [s0, w0]
  
  	if N>D:
  		sv = linalg.svd(Y, full_matrices=0);
  		[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
  		v = s0.std(axis=0)
  		s0 /= v;
  		w0 *= v;
  		return [s0, w0]
  	else:
  		K=np.cov(Y)
  		sv = linalg.eigh(K)
  		std_var = np.sqrt(sv[0])
  		pc = sv[1]*std_var[np.newaxis(),0]
  		#import ipdb
  		#ipdb.set_trace()
  		return [pc,std_var]