StatRed ass2: Implemented Reconstruct() function.

statred/ass2/pca.py

@@ -7,24 +7,21 @@ def sortedeig(M):
     return (d[si], U[:,si])
 
 def calc_PCA(**kwargs):
-    data = kwargs['data']
-
-    if data == 'natural':
+    if kwargs['data'] == 'natural':
         X = loadtxt('natural400_700_5.asc').T
         N = 219
-    elif data == 'munsell':
+    elif kwargs['data'] == 'munsell':
         X = loadtxt('munsell380_800_1.asc').reshape(1269, 421).T
         N = 1269
     else:
         raise ValueError('Undefined data set.')
 
-    Xzm = X - tile(mean(X, 1)[:, newaxis], N)
+    Xzm = X - tile(mean(X, 1)[:,newaxis], N)
     S = dot(Xzm, Xzm.T) / (N - 1)
     return sortedeig(S)
 
 def PCA(**kwargs):
     d, U = calc_PCA(**kwargs)
-
     figure(1)
     plot(d)
     show()
@@ -35,19 +32,41 @@ def EigenImages(k, **kwargs):
         min, max, step = 400, 701, 5
     elif kwargs['data'] == 'munsell':
         min, max, step = 380, 801, 1
-    else:
-        raise ValueError('Undefined data set.')
-
     figure(2)
     a = arange(min, max, step)
     for i in range(k):
-        plot(a, U[:,i], label=str(i+1))
+        plot(a, U[:,i], label=str(i + 1))
     legend()
     show()
 
+def Reconstruct(k, sample, **kwargs):
+    d, U = calc_PCA(**kwargs)
+    if kwargs['data'] == 'natural':
+        X = loadtxt('natural400_700_5.asc').T
+        min, max, step = 400, 701, 5
+    elif kwargs['data'] == 'munsell':
+        X = loadtxt('munsell380_800_1.asc').reshape(1269, 421).T
+        min, max, step = 380, 801, 1
+
+    # Select the specified vector, subtract the mean from it and multiply with
+    # the transposed eigenvector basis to get the coordinates with respect to U
+    # Then, take the first k components and try to reconstruct the original
+    # spectrum
+    x = X[:,sample]
+    xbar = mean(X, 1)
+    yzm = dot(U.T, x - xbar)[:k]
+    x_k = dot(U[:,:k], yzm[:k]) + xbar
+
+    a = arange(min, max, step)
+    figure(3)
+    plot(a, x, label='original')
+    plot(a, x_k, label='reconstructed')
+    legend()
+    show()
 
 #PCA(data='natural')
 #PCA(data='munsell')
 
-EigenImages(5, data='natural')
+#EigenImages(5, data='natural')
 
+Reconstruct(5, 23, data='natural')