Cifar

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");
PARAM_STRING_IN("dataset", "path/to/dataset.", "i", "mnist.arm");
PARAM_STRING_IN("label", "path/to/label.", "l", "mnist.arm");
PARAM_INT_IN("dataset-max-cols", "dataset max cols.", "m",  -1);
PARAM_INT_IN("poolSize", "poolSize", "p",  1);
PARAM_DOUBLE_IN("slabPenalty", "slabPenalty.", "a",  10.5);
PARAM_INT_IN("num-patches", "num patches", "np", 80);
PARAM_INT_IN("hidden-layer-size", "hidden layer size.", "G",  128);
PARAM_INT_IN("batch-size", "batch-size.", "b", 100);
PARAM_DOUBLE_IN("step-size", "learning rate.", "r", 0.01);
PARAM_INT_IN("num-epoches", "num-epoches.", "e", 20);
PARAM_DOUBLE_IN("tolerance", "tolerance.", "t", 1e-5);
PARAM_DOUBLE_IN("lambdaBias", "lambdaBias", "q",  10);
PARAM_FLAG("shuffle", "shuffle or not", "s");


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

  arma::mat dataset, labels;
  size_t datasetMaxCols = 500;
  size_t hiddenLayerSize = 256;
  size_t numPatches = 49 * 500;
  size_t batchSize = 10;
  double stepSize = 0.0001;
  size_t numEpoches = 10;
  double tolerance = 1e-5;
  bool shuffle = true;
  size_t poolSize = 4;
  double slabPenalty = 10.5;
  double lambdaBias = 10;

  if (!CLI::HasParam("dataset"))
    Log::Fatal << "Input dataset is undefined!" << std::endl;
  if (!CLI::HasParam("label"))
    Log::Fatal << "Input lable  is undefined!" << std::endl;

  std::string datafile = CLI::GetParam<std::string>("dataset");
  std::string labelfile = CLI::GetParam<std::string>("label");


  if (CLI::HasParam("dataset-max-cols"))
    datasetMaxCols = CLI::GetParam<int>("dataset-max-cols");

  if (CLI::HasParam("hidden-layer-size"))
    hiddenLayerSize = CLI::GetParam<int>("hidden-layer-size");

  if (CLI::HasParam("num-patches"))
    numPatches = CLI::GetParam<int>("num-patches");

  if (CLI::HasParam("batch-size"))
    batchSize = CLI::GetParam<int>("batch-size");

  if (CLI::HasParam("step-size"))
    stepSize = CLI::GetParam<double>("step-size");

  if (CLI::HasParam("num-epoches"))
    numEpoches = CLI::GetParam<int>("num-epoches");

  if (CLI::HasParam("tolerance"))
    tolerance = CLI::GetParam<double>("tolerance");
  if (CLI::HasParam("lambdaBias"))
    lambdaBias = CLI::GetParam<double>("lambdaBias");
  if (CLI::HasParam("slabPenalty"))
    slabPenalty = CLI::GetParam<double>("slabPenalty");
  if (CLI::HasParam("poolSize"))
    poolSize = CLI::GetParam<double>("poolSize");


  shuffle = CLI::HasParam("shuffle");

  Log::Info << "dataset = '" << dataset << "'" << std::endl;
  
  Log::Info << std::boolalpha
      << "HiddenLayerSize = " << hiddenLayerSize
      << " batchSize = " << batchSize
      << " stepSize = " << stepSize
      << " numEpoches = " << numEpoches
      << " tolerance = " << tolerance
      << " shuffle = " << shuffle << std::endl;

  arma::mat tempData;
  double radius = 0;
  double tempRadius = 0;
  // number of patches
  size_t numImages = numPatches / 49;
  std::cout << "numImages" << numImages << std::endl;
  // 10%
  size_t trainSize = numImages - numImages / 10;
  // 10%
  size_t labelSize = numImages - numImages / 10;
  size_t testSize = numPatches - trainSize;

  
  dataset.load(datafile);
  labels.load(labelfile);

  std::cout << "dataset size "<< arma::size(dataset) << std::endl;
  std::cout << dataset.n_cols << std::endl;
  std::cout << "trainSize" << trainSize << std::endl;
  std::cout << "label size " << arma::size(labels) << std::endl;

  assert(labels.n_rows >= numImages);
  std::cout << "numPatches" << numPatches << std::endl;
  std::cout << "dataset.n_cols" << dataset.n_cols << std::endl;
  assert(dataset.n_cols >= numPatches);

  arma::mat tempTrainData(192, numPatches);

  arma::mat trainData = arma::mat(dataset.memptr(), 192, trainSize, false, false);
  arma::mat testData = arma::mat(dataset.memptr() + trainData.n_elem, 192, testSize,
      false, false);

  arma::mat trainLabelsTemp = arma::mat(labels.memptr(), labelSize , 1, false, false);
  arma::mat testLabelsTemp = arma::mat(labels.memptr() + trainLabelsTemp.n_elem,
      numImages - labelSize, 1, false, false);

  std::cout << trainLabelsTemp(trainLabelsTemp.n_elem - 1) << std::endl;
  std::cout << testLabelsTemp(testLabelsTemp.n_elem - 1) << std::endl;
  
  GaussianInitialization gaussian(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)
  {
    testLabelsTemp(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 *=3;

  
  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(10);
  modelssRBM.Policy().SpikeBias().fill(-5);
  modelssRBM.Train(trainData, msgd);

  // hiddenLayer size * number of patches in one image
  arma::cube hiddenMeanTest(hiddenLayerSize , 49, testData.n_cols);
  arma::cube hiddenMeanTrain(hiddenLayerSize , 49, trainData.n_cols);
  arma::mat ssRBMFeaturesTest(hiddenLayerSize * 49, testData.n_cols);
  arma::mat ssRBMFeaturesTrain(hiddenLayerSize * 49, trainData.n_cols);

  

  ;
  for (size_t i = 0, j = 0; i < testData.n_cols; i++, j++)
  {
    j = j % 49;
    modelssRBM.SampleHidden(std::move(testData.col(i)), 
        std::move(hiddenMeanTest.slice(i).col(j)));
    j++;
  }

  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, j = 0; i < trainData.n_cols; i++, j++)
  {
    j = j % 49;
    modelssRBM.SampleHidden(std::move(trainData.col(i)), 
        std::move(hiddenMeanTrain.slice(i).col(j)));
    j++;
  }
  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;

  arma::mat normalTrainFeat(trainData.n_rows * 49, trainData.n_cols);
  arma::mat normalTestFeat(trainData.n_rows * 49, testData.n_cols);
  
  trainData.resize(trainData.n_rows * 49, trainData.n_cols);
  testData.resize(trainData.n_rows * 49, testData.n_cols);

  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 regressor2(normalTrainFeat, trainLabels, numClasses,
        0.001, false, optimizer1);

  SoftmaxRegression regressor1(ssRBMFeaturesTrain, trainLabels, numClasses,
        0.001, false, optimizer1);
  
  double classificationAccurayssRBM = regressor1.ComputeAccuracy(ssRBMFeaturesTest, testLabels);
  double classificationAccurayNormal = regressor2.ComputeAccuracy(normalTestFeat, testLabels);

  std::cout << "RBM Accuracy" << classificationAccurayssRBM << std::endl;
  std::cout << "noraml Accuracy" << classificationAccurayNormal << std::endl;

}

generate_patches.py

import numpy as np
import cPickle

def wrapper(file):
  with open(file, 'rb') as fo:
      dict = cPickle.load(fo)
  return getPatches(np.array(dict['data']))
  

def getPatches(image):
  print image.shape[0]
  patches = np.zeros((8, 8, 3, 49, image.shape[0]))
  for img in range(0, image.shape[0]):
    img_data = image[img].reshape(32, 32, 3)
    for i in range(0, 7):
      for j in range(0, 7):
        for channel in range(0, 3):
          patches[:,:,channel, i * 7 + j, img] = img_data[i*4 : i*4 + 8, j*4 : j*4+8, channel]

  return preprocess(patches)


def preprocess(data):
  for img in range(0, data.shape[4]):
    for i in range(0, data.shape[3]):
      for channel in range(0, data.shape[2]):
        data[:, :, channel, i, img] -= data[:, :, i, channel, img].mean()
        data[:, :, channel, i, img] /= np.var(data[:, :, i, channel, img])
        data[:, :, channel, i, img] = np.dot(zca_whitening_matrix(data[:, :, i, channel, img]), data[:, :, i, channel, img])
        return data

def zca_whitening_matrix(X):
    """
    Function to compute ZCA whitening matrix (aka Mahalanobis whitening).
    INPUT:  X: [M x N] matrix.
        Rows: Variables
        Columns: Observations
    OUTPUT: ZCAMatrix: [M x M] matrix
    """
    # Covariance matrix [column-wise variables]: Sigma = (X-mu)' * (X-mu) / N
    sigma = np.cov(X, rowvar=True) # [M x M]
    # Singular Value Decomposition. X = U * np.diag(S) * V
    U,S,V = np.linalg.svd(sigma)
        # U: [M x M] eigenvectors of sigma.
        # S: [M x 1] eigenvalues of sigma.
        # V: [M x M] transpose of U
    # Whitening constant: prevents division by zero
    epsilon = 0.1
    # ZCA Whitening matrix: U * Lambda * U'
    ZCAMatrix = np.dot(U, np.dot(np.diag(1.0/np.sqrt(S + epsilon)), U.T)) # [M x M]
    return ZCAMatrix


def main():
  data = wrapper("/Users/kris/Downloads/cifar-10-batches-py/data_batch_1")
  
  data = data[:, :, :, :, 0:100]
  data = data.reshape(192, 49 * 100)
  print data.shape
  np.savetxt('batch1patches.txt', data)

if __name__ == '__main__':
  main()