############################################################
# Problem 2

# All images will be of this size (imageDim[1] rows, imageDim[2] cols)
# Returns the data matrix, where each column is an image
loadFaces = function(directory) {
  files = dir(directory, full.names = TRUE)

  # Read in data
  faces = c()
  n = length(files)
  for (i in 1:n) {
    print(files[i])
    pic = getChannels(read.pnm(files[i])) # Get the matrix of pixels
    imageDim <<- dim(pic)
    pic = c(pic); d = length(pic) # Linearize the pixels
    if (length(faces) == 0) faces = matrix(0, d, n)
    faces[,i] = pic
  }
  faces
}

errorRate = function(model, dataset) {
  confusion = table(predict(model, dataset$x), dataset$y)
  (confusion[1,2]+confusion[2,1])/length(dataset$y)
}

# augmentFeatures: takes a data matrix and transforms it
evaluateFeatures = function(augmentFeatures, trainPos, trainNeg, testPos, testNeg) {
  createSet = function(pos, neg) {
    list(
      x = t(cbind(pos, neg)),
      y = t(cbind(t(rep(1, dim(pos)[2])), t(rep(0, dim(neg)[2]))))
    )
  }
  
  # Augment feature vectors
  trainPos = augmentFeatures(trainPos)
  trainNeg = augmentFeatures(trainNeg)
  testPos = augmentFeatures(testPos)
  testNeg = augmentFeatures(testNeg)

  # Combine positive and negative examples together
  train = createSet(trainPos, trainNeg)
  test = createSet(testPos, testNeg)

  # Run SVM classifier
  c = 1/mean(sqrt(rowSums(train$x * train$x)))
  model = svm(train$x, factor(train$y), type="C", kernel="linear", cost=c)
  list(
    trainError = errorRate(model, train),
    testError = errorRate(model, test)
  )
}

# Given a vector of pixels, plot it.
plotImg = function(pixels) {
  plot(pixmapGrey(matrix((pixels-min(pixels))/(max(pixels)-min(pixels)), imageDim[1], imageDim[2])))
}

loadAllFaces = function() {
  trainPos <<- loadFaces('data/faces/train/positive')
  trainNeg <<- loadFaces('data/faces/train/negative')
  testPos  <<- loadFaces('data/faces/test/positive')
  testNeg  <<- loadFaces('data/faces/test/negative')
}

# Assue loadAllFaces() has already been called.
p2 = function() {
  library(pixmap)
  library(e1071)

  computePCA = function(X) {
    svd(X)
  }

  # Comput PCA features on positive training examples
  pca = computePCA(trainPos)

  augmentFeatures = function(X) {
    # INSERT CODE HERE (use the pca variable computed above)
    X # Use original pixel features
  }
  print(evaluateFeatures(augmentFeatures, trainPos, trainNeg, testPos, testNeg))
}

############################################################
# Problem 4

p4 = function() {
  library(tuneR)
  library(fastICA)

  d = 4 # Number of sources
  n = 88202 # Number of time examples

  # Read the wav files into the data matrix
  print("Reading wave files...")
  X = matrix(0, n, d)
  for (i in 1:d) {
    X[,i] = readWave(paste("data/sounds/x", i, ".wav", sep=""))@left
  }

  # Compute ICA
  print("Computing ICA...")
  # REPLACE CODE BELOW (should be one line of code)
  S = X # Just use original signals as sources

  # Output sources to current directory (s1.wav, s2.wav...)
  print("Writing sources...")
  for (i in 1:d) {
    writeWave(normalize(Wave(S[,i], bit=16, samp.rate=44100), unit=16), paste("s", i, ".wav", sep=""))
  }
}