###############################################################################
### 1
###############################################################################

## The image matrix for training sample. 256x1707
azip <- read.table("azip.dat")

## The true digits given in the training sample. lenght = 1707
dzip <- as.numeric(read.table("dzip.dat"))

## The testing image matrix. 256x2007
testzip <- read.table("testzip.dat")

## The true digits for the testing sample. length = 2007
dtest <- read.table("dtest.dat")

## Display the image
i <- 120
image(matrix(azip[ , i], ncol = 16)[ , 1:16], col = gray(255:0/255))

###----------------------------------------------------------------------------
### The naive method
###----------------------------------------------------------------------------

## The mean of training sample of a single digit
digits <- 0:9 # The possible digits in the US postal code
img.mean <- matrix(0, 256, length(digits))

for(i in digits)
  {
    idx <- (i == dzip) # the location indicator for the ith digit
    imgi <- azip[, idx, drop = FALSE]
    imgi.mean <- rowMeans(imgi)
    img.mean[,  i+1] <- imgi.mean
  }


## Plot the mean image
par(mfrow = c(2, 5))
for(i in 1:10)
  {
    image(matrix(img.mean[, i], 16)[, 16:1],
          col = gray(255:0/255))
  }


## Sketch a distance function to compute the Euclidean distance between two
## matrices in row wise.
rdist <- function(X, Y)
  {
    dim.X <- dim(X)
    dim.Y <- dim(Y)
    sum.X <- matrix(rowSums(X^2), dim.X[1], dim.Y[1])
    sum.Y <- matrix(rowSums(Y^2), dim.X[1], dim.Y[1], byrow = TRUE)
    dist0 <- sum.X + sum.Y - 2*tcrossprod(X, Y)
    out <- sqrt(dist0)
    return(out)
  }

## For an unknown testing digit image, compare the distance to the means
test.sample <- 1:5

## Let's first plot those testing image
par(mfcol = c(ceiling(length(test.sample)/5), 5)) # five columns

for(i in test.sample)
  {
    image(matrix(testzip[, i], 16)[, 16:1], col = gray(255:0/255))
  }

## Calculate the distance from testing sample to the mean in the training
## sample.
img.dist <- rdist(t(testzip[, test.sample]), t(img.mean))


## The classification results by the naive method
apply(img.dist, 1, which.min)-1

###----------------------------------------------------------------------------
### The SVD method
###----------------------------------------------------------------------------

## Compute the singular matrix of a single digit in the training sample
digit <- 9

## Subtract the matrix for that digit
img.mat <- azip[, digit == dzip, drop = FALSE]
img.matSVD <- svd(img.mat)


X = img.matSVD$u
y = testzip[, 1]


img.mat <- azip[, dzip == 7, drop = FALSE]
img.matSVD <- svd(img.mat)
X = img.matSVD$u


## Plot the singular matrix under different basis.
par(mfrow = c(2, 5))
for(i in 1:10)
  {
    image(matrix(img.matSVD$u[, i], 16)[, 16:1], col = gray(255:0/255), main =
          paste("singular image ", i, sep = ""))
  }

## Do the least square method with different basis and find the minimal
## residuals.

## The testing digit matrix
test.idx <- 1
image(matrix(testzip[, test.idx], 16)[, 16:1], col = gray(255:0/255), main =
      paste("Testing digit"))

resid.norm <- matrix(NA, 10, 1, dimnames = list(0:9,"resid"))
for(i in 0:9)
  {
    img.mat <- azip[, i == dzip, drop = FALSE]
    img.matSVD <- svd(img.mat)
    basis.max <- ncol(img.matSVD$u)
    resid.norm[i+1, ] <- norm(matrix(lm(testzip[, test.idx]~
                                        0+img.matSVD$u[, 1:basis.max])$resid), "F")
  }
resid.norm

## Find the optimum basis function.

test.idx <- 1:100
basis.max <- 1:20
nCorrectOut <- numeric(length(test.idx))
resid.norm <- matrix(0, 10, 1, dimnames = list(0:9,"resid"))
for(b in basis.max)
  {
    nCorrect <- 0
    for(k in test.idx)
      {
        for(i in 0:9)
          {
            img.mat <- azip[, i == dzip, drop = FALSE]
            img.matSVD <- svd(img.mat)
            resid.norm[i+1, ] <- norm(matrix(lm(testzip[, test.idx[k]]~
                                                0+img.matSVD$u[, 1:basis.max[b]])$resid), "F")
          }
        if((which.min(resid.norm)-1)  == dtest[k])
          {
            nCorrect <- nCorrect+1
          }
      }
    nCorrectOut[b] <- nCorrect
    print(b)
  }


###----------------------------------------------------------------------------
### Regression with QR,  SVD
###----------------------------------------------------------------------------

y <- matrix(1:5)
X <- cbind(1, rnorm(5), runif(5), 1+runif(5)/1000)

## Traditional method
betaHat1 <- solve(t(X)%*%X)%*%t(X)%*%y

## QR
XQR <- qr(X)
Q <- qr.Q(XQR)
R <- qr.R(XQR)

XPlus2 <- solve(R)%*%t(Q)

XPlus2%*%X

betaHat2 <- XPlus2%*%y


## SVD

XSVD <- svd(X)
U <- XSVD$u
V <- XSVD$v
Sigma <- diag(XSVD$d)
Sigma1 <- diag(1/XSVD$d)

XPlus3 <- V%*%Sigma1%*%t(U)

betaHat3 <- XPlus3%*%y


cbind(betaHat1, betaHat2, betaHat3)