<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Procrustes Analysis
Step1:  PCA
Step2:  Solving b vector
$$b = \Phi^T \left(x - \bar{x}\right)$$
<ASSISTANT_TASK:>
Python Code:
import pandas
df = pandas.read_csv('muct76-opencv.csv', header=0, usecols=np.arange(2,154), dtype=float)
df.head()
X = df.iloc[:, ::2].values
Y = df.iloc[:, 1::2].values
d = np.hstack((X,Y))
d.shape
import sys
threshold = 1.0e-8
def center(vec):
    pivot = int(vec.shape[0]/2)
    meanx = np.mean(vec[:pivot])
    meany = np.mean(vec[pivot:])
    return(meanx, meany)
def calnorm(vec):
    vsqsum = np.sum(np.square(vec))
    return(vsqsum)
def scale(vec):
    vcopy = vec.copy()
    vmax = np.max(vec)
    if vmax > 2.0:
        vcopy = vcopy / vmax
    vnorm = calnorm(vcopy)
    return (vcopy / np.sqrt(vnorm))
def caldiff(pref, pcmp):
    return np.mean(np.sum(np.square(pref - pcmp), axis=1))
def simTransform(pref, pcmp, showerror = False):
    err_before = np.mean(np.sum(np.square(pref - pcmp), axis=1))
    ref_mean = np.mean(pref, axis=0)
    prefcentered = np.asmatrix(pref) - np.asmatrix(ref_mean)
    
    cmp_mean = np.mean(pcmp, axis=0)
    pcmpcentered = np.asmatrix(pcmp) - np.asmatrix(cmp_mean)   
    
    Sxx = np.sum(np.square(pcmpcentered[:,0]))
    Syy = np.sum(np.square(pcmpcentered[:,1]))
    Sxxr = prefcentered[:,0].T * pcmpcentered[:,0] #(ref_x, x)
    Syyr = prefcentered[:,1].T * pcmpcentered[:,1] #(ref_y, y)
    Sxyr = prefcentered[:,1].T * pcmpcentered[:,0] #(ref_y, x)
    Syxr = prefcentered[:,0].T * pcmpcentered[:,1] #(ref_x, y)
    a = (Sxxr + Syyr)/(Sxx + Syy) #(Sxxr + Syyr) / (Sxx + Syy)
    b = (Sxyr - Syxr) / (Sxx + Syy)
    a = np.asscalar(a)
    b = np.asscalar(b)
    Rot = np.matrix([[a, -b],[b, a]])
    translation = -Rot * np.asmatrix(cmp_mean).T + np.asmatrix(ref_mean).T
    outx, outy = [], []
    res = Rot * np.asmatrix(pcmp).T + translation
    err_after = np.mean(np.sum(np.square(pref - res.T), axis=1))
    if showerror:
        print("Error before: %.4f    after: %.4f\n"%(err_before, err_after))
    return (res.T, err_after)
    
    
def align2mean(data):
    d = data.copy()
    pivot = int(d.shape[1]/2)
    for i in range(d.shape[0]):
        cx, cy = center(d[i,:])
        d[i,:pivot] = d[i,:pivot] - cx
        d[i,pivot:] = d[i,pivot:] - cy
        #print(cx, cy, center(d[i,:]))
        d[i,:] = scale(d[i,:])
        norm = calnorm(d[i,:])
    
    d_aligned = d.copy()
    pref = np.vstack((d[0,:pivot], d[0,pivot:])).T
    print(pref.shape)
    mean = pref.copy()
    mean_diff = 1
    while mean_diff > threshold:
        err_sum = 0.0
        for i in range(1, d.shape[0]):
            p = np.vstack((d[i,:pivot], d[i,pivot:])).T
            p_aligned, err = simTransform(mean, p)
            d_aligned[i,:] = scale(p_aligned.flatten(order='F'))
            err_sum += err
        
        oldmean = mean.copy()
        mean = np.mean(d_aligned, axis=0)
        mean = scale(mean)
        mean = np.reshape(mean, newshape=pref.shape, order='F')
        d = d_aligned.copy()
        mean_diff = caldiff(oldmean, mean)
        sys.stdout.write("SumError: %.4f MeanDiff: %.6f\n"%(err_sum, mean_diff))
    return (d_aligned, mean)
        
    
d_aligned, mean = align2mean(d)
plt.figure(figsize=(7,7))
plt.gca().set_aspect('equal')
plotFaceShape(mean)
Explanation: Procrustes Analysis:
Translate each example such that its  centroid is at origin, and scale them so that $|\mathbf{x}|=1$
Choose the first example as the initial estimate of the mean shape ($\bar{\mathbf{x}}$), and save a copy as a reference ($\mathbf{\bar{x}_0}$)
Align all the examples with the current estimate of the mean ($\bar{\mathbf{x}}$)
Re-estimate mean from the aligned shapes
Apply constraints on the current estiate of the mean by aligning it with reference $\mathbf{\bar{x}_0}$ and scaling so that $|\bar{\mathbf{x}}|=1$
Repeat steps 3, 4, 5 until convergence
End of explanation
d_aligned.shape
from sklearn.decomposition import PCA
pca = PCA(n_components=8)
pca.fit(d_aligned)
print(pca.explained_variance_ratio_) 
cov_mat = np.cov(d_aligned.T)
print(cov_mat.shape)
eig_values, eig_vectors = np.linalg.eig(cov_mat)
print(eig_values.shape, eig_vectors.shape)
num_eigs = 8
Phi_matrix = eig_vectors[:,:num_eigs]
Phi_matrix.shape
Explanation: PCA
End of explanation
# * ()
mean_matrix = np.reshape(mean, (152,1), 'F')
d_aligned_matrix = np.matrix(d_aligned)
delta = d_aligned_matrix.T - mean_matrix
b = (np.matrix(Phi_matrix).T * delta).T
b.shape
mean.dump('models/meanshape-ocvfmt.pkl')
eig_vectors.dump('models/eigenvectors-ocvfmt.pkl')
eig_values.dump('models/eigenvalues-ocvfmt.pkl')
Phi_matrix.dump('models/phimatrix.pkl')
b.dump('models/bvector.pkl')
d_aligned.dump('models/alignedfaces.pkl')
mean_matrix.dump('models/meanvector.pkl')
Explanation: Solving b vector
$$b = \Phi^T \left(x - \bar{x}\right)$$
End of explanation
<END_TASK>