def C(w, Xs):
n = len(Xs)<EOL>P = projection_matrix(w)<EOL>Ys = [np.dot(P, X) for X in Xs]<EOL>A = calc_A(Ys)<EOL>A_hat = calc_A_hat(A, skew_matrix(w))<EOL>return np.dot(A_hat, sum(np.dot(Y, Y) * Y for Y in Ys)) / np.trace(np.dot(A_hat, A))<EOL>
Calculate the cylinder center given the cylinder direction and 
    a list of data points.