normal_cifar.cpp

#include <mlpack/core.hpp>

#include <mlpack/core/optimizers/rmsprop/rmsprop.hpp>
#include <mlpack/methods/ann/init_rules/gaussian_init.hpp>
#include <mlpack/methods/ann/rbm/spike_slab_rbm_policy.hpp>
#include <mlpack/methods/ann/rbm.hpp>
#include <mlpack/methods/ann/rbm/binary_rbm_policy.hpp>
#include <mlpack/methods/softmax_regression/softmax_regression.hpp>
#include <mlpack/core/optimizers/minibatch_sgd/minibatch_sgd.hpp>
#include <mlpack/core/optimizers/sgd/sgd.hpp>
#include <mlpack/core/optimizers/lbfgs/lbfgs.hpp>

using namespace mlpack;
using namespace mlpack::ann;
using namespace mlpack::optimization;
using namespace mlpack::regression;

PROGRAM_INFO("ssRBM", "cifar");


int main(int argc, char** argv)
{
  CLI::ParseCommandLine(argc, argv);

  size_t hiddenLayerSize = 256;
  size_t numPatches = 100;
  size_t batchSize = 20;
  double stepSize = 0.001;
  size_t numEpoches = 10;
  double tolerance = 1e-5;
  bool shuffle = true;
  size_t poolSize = 4;
  double slabPenalty = 10.5;
  double lambdaBias = 10;
  double radius = 0;
  double tempRadius = 0;

  /*
   * Training Phase
   * We have 40000 patches for 400 training images and 192 feature vector.
   * 
   */

  std::string trainDataFile = "./random_train.txt";
  std::string testDataFile = "./random_test.txt";
  std::string LabelsFile = "./labels.txt";

  arma::mat trainData, testData, trainLabelsTemp, testLabelsTemp;
  
  trainData.load(trainDataFile);
  testData.load(testDataFile);
  trainLabelsTemp.load(LabelsFile);
  // since we are training on 500 images
  testLabelsTemp = trainLabelsTemp.rows(400, 499);
  trainLabelsTemp = trainLabelsTemp.rows(0, 399);

  std::cout << testLabelsTemp(0) << std::endl;

  trainData = trainData.cols(0, 10);
  testData = testData.cols(0, 10);


  trainData = trainData.t();
  testData = testData.t();
  
  GaussianInitialization gaussian(0, 0.1);

  arma::Row<size_t> trainLabels = arma::zeros<arma::Row<size_t>>(1,
      trainLabelsTemp.n_rows);
  arma::Row<size_t> testLabels = arma::zeros<arma::Row<size_t>>(1,
      testLabelsTemp.n_rows);

  for (size_t i = 0; i < trainLabelsTemp.n_rows; ++i)
  {
    trainLabels(i) = arma::as_scalar(trainLabelsTemp.row(i));
  }

  for (size_t i = 0; i < testLabelsTemp.n_rows; ++i)
  {
    testLabels(i) = arma::as_scalar(testLabelsTemp.row(i));
  }

  // Calculate radius
  for (size_t i = 0; i < trainData.n_cols; i++)
  {
    tempRadius = arma::norm(trainData.col(i));
    if (radius < tempRadius)
      radius = tempRadius;
  }
  radius *= 1.5;

  
  size_t iter =  (numEpoches * trainData.n_cols) / batchSize;
  std::cout << "iter#" << iter << std::endl;  

  SpikeSlabRBMPolicy<> rbmPolicy(trainData.n_rows, hiddenLayerSize, poolSize,
      slabPenalty, radius);

  RBM<GaussianInitialization, SpikeSlabRBMPolicy<>> modelssRBM(trainData,
      gaussian, rbmPolicy, 1, 1, true, false);

  MiniBatchSGD msgd(batchSize, stepSize, iter, tolerance, shuffle);
  modelssRBM.Reset();
  modelssRBM.Policy().VisiblePenalty().fill(lambdaBias);
  modelssRBM.Policy().SpikeBias().fill(-5);
  modelssRBM.Train(trainData, msgd);
  
  // hiddenLayer size * number of patches in one image
  // 256 * 100 * 1000
  arma::cube hiddenMeanTest(hiddenLayerSize , numPatches, testData.n_cols / numPatches);
  arma::cube hiddenMeanTrain(hiddenLayerSize , numPatches, trainData.n_cols / numPatches);
  arma::mat ssRBMFeaturesTest(hiddenLayerSize * numPatches, testData.n_cols / numPatches);
  arma::mat ssRBMFeaturesTrain(hiddenLayerSize * numPatches, trainData.n_cols / numPatches);
  // slice size 256 * 100 * 1000
  for (size_t i = 0; i < hiddenMeanTest.n_slices; i++)
    for (size_t k = 0; k < numPatches; k++)
    {
      arma::mat tmp;
      modelssRBM.SampleHidden(std::move(testData.col(i * numPatches + k)),
          std::move(tmp));
      hiddenMeanTest.slice(i).col(k) = tmp.rows(0, hiddenLayerSize - 1);
    }

  for (size_t i = 0; i < hiddenMeanTest.n_slices; ++i)
  {
    ssRBMFeaturesTest.col(i) = arma::vectorise(hiddenMeanTest.slice(i), 1).t();
  }

  for (size_t i = 0; i < hiddenMeanTrain.n_cols; i++)
    for (size_t k = 0; k < numPatches; k++)
    {
      
      arma::mat tmp;
      modelssRBM.SampleHidden(std::move(trainData.col(i * numPatches + k)), 
          std::move(tmp));
      hiddenMeanTrain.slice(i).col(k) = tmp.rows(0, hiddenLayerSize - 1);;
      
    }
  for (size_t i = 0; i < hiddenMeanTrain.n_slices; ++i)
  {
    ssRBMFeaturesTrain.col(i) = arma::vectorise(hiddenMeanTrain.slice(i), 1).t();
  }

  std::cout << "ssRBMFeaturesTrain = " << arma::size(ssRBMFeaturesTrain) << std::endl;
  std::cout << "ssRBMFeaturesTest = " << arma::size(ssRBMFeaturesTest) << std::endl;
  

  const size_t numClasses = 10; // Number of classes.
  const size_t numBasis = 5; // Parameter required for L-BFGS algorithm.
  const size_t numIterations = 1000; // Maximum number of iterations.


  arma::Row<size_t> predictions1, predictions2;
  L_BFGS optimizer1(numBasis, numIterations);
  SoftmaxRegression regressor1(ssRBMFeaturesTrain, trainLabels, numClasses,
  0.001, false, optimizer1);
  double classificationAccurayssRBM = regressor1.ComputeAccuracy(ssRBMFeaturesTest, testLabels);
  arma::Row<size_t> tmpLabel;
  regressor1.Predict(ssRBMFeaturesTest, tmpLabel);
  std::cout << "tmpLabel Size = " << arma::size(tmpLabel) << std::endl;
  tmpLabel.print();
  std::cout << "Real lables" << std::endl;
  testLabels.print();
  std::cout << "RBM Accuracy" << classificationAccurayssRBM << std::endl;
  
}

patches.py

import numpy as np
import cPickle
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

def wrapper(file):
  with open(file, 'rb') as fo:
      dict = cPickle.load(fo)
  print np.array(dict['data']).shape
  data = np.array(dict['data'])
  data = data.reshape((-1,3,32,32))
  # output is of shape (N, 32, 32, 3)
  data = data.transpose(0,2,3,1)
  plt.imshow(data[6])
  plt.show()
  # output is of shape (N, 32, 32, 3)
  # put data back into a design matrix (N, 3072)
  data = data.reshape(-1, 3072)
  data = preprocess(data[0:5000][:])
  return getPatches(data)
  

def getPatches(image):
  print "Inside get patches........"
  data = image.reshape((-1, 32, 32, 3))
  print data.shape
  patches = np.zeros((image.shape[0] * 49, 8*8*3))
  for img in range(0, image.shape[0]):
    for i in range(0, 7):
      for j in range(0, 7):
          patch1 = data[img, i*4 : i*4 + 8, j*4 : j*4+8, 0].reshape(64)
          patch2 = data[img, i*4 : i*4 + 8, j*4 : j*4+8, 1].reshape(64)
          patch3 = data[img, i*4 : i*4 + 8, j*4 : j*4+8, 2].reshape(64)
          patches[img*49 + i*7 + j] = np.hstack((patch1, patch2, patch3))
  return patches


def preprocess(data):
  data = data - data.mean(axis = 0)
  data = data / np.sqrt((data**2).sum(axis = 1))[:, None]
  pca = PCA(whiten = True)
  pca.fit_transform(data)
  return data



def main():
  data = wrapper("/Users/kris/Downloads/cifar-10-batches-py/data_batch_1")

  data = data.reshape((-1, 8 * 8 * 3 * 49))
  print data.shape
  np.savetxt("batch1pathces.txt", data)
  data = data.reshape(-1, 49, 3, 8, 8)
  data = data.transpose(0,1,4,3,2 )
  print data.shape
  
  m,M = data[6][0].min(), data[6][0].max()
  plt.imshow((data[6][0] - m) / (M - m))
  plt.show()
  
  
  # print data[:, 1].mean()
  # np.savetxt('batch1patches.txt', data)

if __name__ == '__main__':
  main()

random_patch.py

import numpy as np
import matplotlib.pyplot as plt
import cPickle
from sklearn.decomposition import PCA

"""Get Random samples"""
def wrapper(file, random_samples):
  with open(file, 'rb') as fo:
      dict = cPickle.load(fo)
  print np.array(dict['data']).shape
  data = np.array(dict['data'])
  data = data.reshape((-1,3,32,32))
  # output is of shape (N, 32, 32, 3)
  data = data.transpose(0,2,3,1)
  plt.imshow(data[6])
  plt.show()
  # put data back into a design matrix (N, 3072)
  data = data.reshape(-1, 3072)
  data = preprocess(data[:5000][:])
  return getPatches(data, random_samples)

def getPatches(data, random_samples):
  data = data.reshape((-1, 32, 32, 3))
  allpatch = np.zeros((data.shape[0] * random_samples, 64 * 3))
  for i in range(0, data.shape[0]):
    for j in range(0, random_samples):
      starting_point = np.random.randint(data.shape[1] - 8) # get starting point(row #)
      patch1 = data[i, starting_point : starting_point + 8, starting_point:starting_point + 8, 0].reshape(64)
      patch2 = data[i, starting_point : starting_point + 8, starting_point:starting_point + 8, 1].reshape(64)
      patch3 = data[i, starting_point : starting_point + 8, starting_point:starting_point + 8, 2].reshape(64)
      allpatch[i*random_samples + j] = np.hstack((patch1, patch2, patch3))
  return allpatch

def preprocess(data):
  data = data - data.mean(axis = 0)
  data = data / np.sqrt((data**2).sum(axis = 1))[:, None]
  pca = PCA(whiten = True)
  pca.fit_transform(data)
  return data


def main():
  random_samples = 10
  data = wrapper("/Users/kris/Downloads/cifar-10-batches-py/data_batch_1", random_samples)
  np.savetxt("randompathces.txt", data)
  data = data.reshape(-1, 3, 8, 8)
  data = data.transpose(0,2,3,1)
  m,M = data[6].min(), data[6].max()
  plt.imshow((data[6] - m) / (M - m))
  plt.show()
  

if __name__ == '__main__':
  main()