Boost.Pythonを使ってC++からSparse Codingする．

人生のうちで，１度や２度，C++のプログラムでSparseCodingしたいときがあると思います．しかし，昨今機械学習の分野で優遇されているのはPythonばかりです．

Sparse CodingについてはPythonではscikit learnを使えば簡単に使えます．
sklearn.decomposition.DictionaryLearning — scikit-learn 0.15.2 documentation

しかし，C++だとMLPACK等のライブラリを別途入れる必要があります．
MLPACKは動作環境がwindowsだったりすると，導入は非常に困難です（私は諦めました）．
また，人生のうち１度や２度しか使わないのでいちいち新しいライブラリを導入するのは正直面倒です．

こんなときは，Boost.Pythonを介してscikit learnを使うという荒業を使えば，
手軽にC++からSparseCodingをすることが可能です．
これはscikit learnやscipyに実装されているものならなんでもC++から使えてしまう方法なので，非常に便利だと思います．

今回は，Boost.Pythonを使ったSparse Codingを使って，手書き数字認識に挑戦してみたいと思います．

概要

手書き数字認識

今回は，データセットとしてMNISTを利用します．
訓練画像・テスト画像とも28x28の手書き数字（0〜9）です．
噂のDeep Learningを使えば99.7%以上の精度が出るみたいです．
すごいですね（小並感）

Sparse Coding

理論面は私もよくわかってないので省略します．
参考：scikit-learnでSparse CodingとDictionary Learning -理論編- - takminの書きっぱなし備忘録

Sparse Codingによる手書き文字認識

文献[1]のインスパイアで，訓練データをそのまま辞書として使ってスパースコーディングします．
テストデータと同じクラスの訓練データの基底の係数が非ゼロになることを期待しています．
もしかすると同じような考えのもとに誰かが手書き文字認識されているかもしれませんが，あまり気にしません．

Boost.Python

本来はPythonからC/C++のコードを呼び出すものですが，逆も可能です．
Python自体にC/C++用のインタフェースがありますが，Boost.Pythonを使うとより簡単に使えます．

実装

sparse_coding.py

# -*- coding: utf-8 -*-
import numpy as np
from sklearn.decomposition import SparseCoder

def sparse_coding(dict,data,transform_n_nonzero_coefs=10,transform_algorithm='lars',transform_alpha=1.0):
    coder = SparseCoder(dictionary=dict,transform_algorithm=transform_algorithm,transform_n_nonzero_coefs=transform_n_nonzero_coefs,transform_alpha=transform_alpha)
    return coder.transform(data)

こちらの関数をC++から呼び出します．

mnist.cpp

OpenCVとBoostを利用しています．

#define BOOST_PYTHON_STATIC_LIB
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

#include <python.h>
#include <numpy/arrayobject.h>
#include <boost/python.hpp>
#include <boost/python/numeric.hpp>
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>

#include <boost/progress.hpp>

#include <iostream>
#include <fstream>
#include <vector>
#include <opencv2/opencv.hpp>

typedef unsigned char BYTE;
const unsigned int NUM_TRAIN = 5000;
const unsigned int NUM_TEST = 500;

inline void endian_swap(unsigned int& x) {
    x = (x>>24) | ((x<<8) & 0x00FF0000) | ((x>>8) & 0x0000FF00) | (x<<24);
}

void loadImages(const std::string& path,cv::OutputArray image,unsigned int n_max=INT_MAX) {
  std::ifstream ifs(path,std::ios::binary);
  if(!ifs){
    std::cout << path <<" is not found." << std::endl;
    exit(1);
  }

  unsigned int magic_number,num_of_items,height,width;
  ifs.read((char*)&magic_number,sizeof(unsigned int));
  ifs.read((char*)&num_of_items,sizeof(unsigned int));
  ifs.read((char*)&height,sizeof(unsigned int));
  ifs.read((char*)&width,sizeof(unsigned int));

  //big-endianなのでlittle-endianに変更
  endian_swap(magic_number);
  endian_swap(num_of_items);
  endian_swap(height);
  endian_swap(width);

  num_of_items = std::min(num_of_items,n_max);
  unsigned int n_bytes=num_of_items*height*width;
  std::vector<BYTE> buf(n_bytes);
  ifs.read((char*)&buf[0],n_bytes);

  std::cout << "=== load images ===" << std::endl;
  std::cout << "path : "<< path  << std::endl;
  std::cout << "magic_number: "<< magic_number  << std::endl;
  std::cout << "num_of_items: "<< num_of_items  << std::endl;
  std::cout << "height: "<< height  << std::endl;
  std::cout << "width: "<< width  << std::endl;

  cv::Mat m(buf);

  m = m.reshape(1,num_of_items);
  m.convertTo(m,CV_64FC1);
  m.copyTo(image);
}

void loadLabels(const std::string& path,cv::OutputArray labels,unsigned int n_max=INT_MAX){
  std::ifstream ifs(path,std::ios::binary);
  if(!ifs){
    std::cout << path <<" is not found." << std::endl;
    exit(1);
  }

  unsigned int magic_number,num_of_items;
  ifs.read((char*)&magic_number,sizeof(unsigned int));
  ifs.read((char*)&num_of_items,sizeof(unsigned int));
  endian_swap(magic_number);
  endian_swap(num_of_items);
  num_of_items = std::min(num_of_items,n_max);

  unsigned int n_bytes=num_of_items;
  std::vector<BYTE> buf(n_bytes);
  ifs.read((char*)&buf[0],n_bytes);

  std::cout << "=== load labels ===" << std::endl;
  std::cout << "path : "<< path  << std::endl;
  std::cout << "magic_number: "<< magic_number  << std::endl;
  std::cout << "num_of_items: "<< num_of_items  << std::endl;

  cv::Mat m(buf);
  m = m.reshape(1,1);
  m.convertTo(m,CV_32SC1);
  m.copyTo(labels);
}


int main(int argc,char* argv[]){
  if(argc!=2){
    std::cout << "usage : " << argv[0] << " MNIST_path" <<std::endl;
    return -1;
  }

  std::string MINIST_path(argv[1]);
  std::string train_images_path = MINIST_path+"/train-images-idx3-ubyte";
  std::string train_labels_path = MINIST_path+"/train-labels-idx1-ubyte";
  std::string test_images_path  = MINIST_path+"/t10k-images-idx3-ubyte";
  std::string test_labels_path  = MINIST_path+"/t10k-labels-idx1-ubyte";

  //dataset load
  cv::Mat train_images, train_labels;
  loadImages(train_images_path,train_images,NUM_TRAIN);
  loadLabels(train_labels_path,train_labels,NUM_TRAIN);

  cv::Mat test_images, test_labels;
  loadImages(test_images_path,test_images,NUM_TEST);
  loadLabels(test_labels_path,test_labels,NUM_TEST);

  //predictのために工夫 : [1,2,3]を[1,0,0; 0,1,0; 0,0,1]に変換
  cv::Mat dot_mat;
  for(int i = 0;i < 10;i++){
    cv::Mat mask = (train_labels == i);
    mask.convertTo(mask,CV_64FC1);
    dot_mat.push_back(mask);
  }
  dot_mat/=255;

  std::cout << "=== test ===" << std::endl;
  {
    //sparce coding
    cv::Mat X;

    try{
      boost::progress_timer t;
      namespace bp = boost::python;

      Py_Initialize();

      // おまじない
      bp::numeric::array::set_module_and_type("numpy", "ndarray");
      bp::class_<std::vector<double>>("std::vector<double>")
        .def(bp::vector_indexing_suite<std::vector<double>>());

      // 利用するpyファイルの指定
      bp::object global_ns = bp::import("__main__").attr("__dict__");
      bp::exec_file("sparse_coding.py", global_ns, global_ns);

      // cv::Matからboost::python::numeric::arrayに変換
      std::vector<double>train_images_vec = train_images.reshape(1,1);
      bp::numeric::array dict(train_images_vec);
      dict.resize(train_images.rows,train_images.cols);

      std::vector<double>test_images_vec = test_images.reshape(1,1);
      bp::numeric::array data(test_images_vec);
      data.resize(test_images.rows,test_images.cols);

      // 呼び出す関数を指定
      bp::object py_func = global_ns["sparse_coding"];

      // pythonの関数を呼び出し
      bp::object py_val = py_func(dict,data,9); 
      
      // 戻り値をcv::Matに変換
      bp::numeric::array py_arr_val = bp::extract <bp::numeric::array>(py_val);

      PyArrayObject* py_arr_obj = (PyArrayObject*)(py_arr_val.ptr());

      double* data_ptr = (double*)(PyArray_DATA(py_arr_obj));

      npy_intp* dims = PyArray_DIMS(py_arr_obj);

      std::vector<double> vec_val(data_ptr, data_ptr + (dims[0]*dims[1]));

      std::cout << vec_val.size() << std::endl;
      cv::Mat m(vec_val);
      m = m.reshape(1,dims[0]);
      m.copyTo(X);// vec_valが開放されても大丈夫なようにdeep copy

      std::cout << "calc time : ";//boost::progres_timer
    }
    catch (...){
      PyErr_Print();
    }

    //predict
    int nT = 0, nF = 0;

    int* gtl=test_labels.ptr<int>(0);
    for(int i = 0;i < test_images.rows;i++){
      std::cout << "[";
      double max_p = 0.0;
      int max_d = -1;
      for(int d = 0;d < 10;d++){
        double p = X.row(i).dot(dot_mat.row(d));
        if(p > max_p){
          max_p = p;
          max_d = d;
        }
        std::cout << p << ",";
      }
      std::cout << "] predict("<< max_d <<")";
      std::cout << " GT(" << gtl[i] << ") " << std::endl;

      (gtl[i] == max_d)? nT++ : nF++;
    }

    std::cout << "Accuracy : " << 1.0*nT/(nT+nF) << std::endl;
  }
}

全データを使うと学習に時間がかかるので，学習データ5000，テストデータ500でやっています．
また，predict部分では係数の合計を認識スコアとして使っています．

ビルド＆使い方

clang++ -std=c++11 -lboost_python $(shell python-config --cflags) $(shell pkg-config --libs opencv) mnist.cpp -o mnist.out
./mnist.out data

↑でdataはデータセットを入れたディレクトリへのパスです．

python-config
||
**出力
>||
（前略）
calc time : 86.88 s
[0,0,0,0,0,0,0,0.658375,0,0,] predict(7) GT(7)
[0,0,0.539294,0.0822132,0,0,0,0,0,0,] predict(2) GT(2)
[0,0.167486,0,0,0,0,0,0,0,0,] predict(1) GT(1)
[0.426151,0,0,0,0,0,0,0,0,0,] predict(0) GT(0)
[-0.0978972,0,0,0,0.235536,0,0,0,0,0.0391382,] predict(4) GT(4)
（中略）
[0,0,0,0,0,0,0,0,0,0.426833,] predict(9) GT(9)
[0,0,0,0,0.147261,0,0,0.0295174,0.275997,0,] predict(8) GT(4)
[0.394697,0,0,0,0,0,0,0,0.0479334,0,] predict(0) GT(0)
[0,0,0,0,0,0,0.19342,0,0.0134353,0,] predict(6) GT(6)
Accuracy : 0.902

精度は90.2%になりました．
うーん・・・学習データを全部使っていないとはいえイマイチですね．
やはりオーソドックスに辞書学習と組み合わせるほうが良さそうですね．

また，sparse_codingのtransform_n_nonzero_coefsの値によって精度が大きく変わりました．
大きくても小さくてもダメなようで，今回は一番精度の良かった9を使いました（これはひどい）

今回のまとめ

Boost.Pythonを使うことでPythonで実装されている関数をC++から使うことができました．
今回はMacで動かしましたが，WindowsでもPyhtonをインストールすれば動かす事ができます．

しかし，体感ですがやはりC++⇔Pythonのデータのやり取りにはそれなりに時間がかかるようです（大きい行列だと1分とか・・・）．
何かをちょこっと試したいという場合に，ライブラリのインストールと格闘とする前に試してみるくらいには使えると思います．

でも，機械学習をしたい人は最初からPythonでプログラミングしましょう．

参考にしたWebサイト

参考文献

[1] J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry, and Y. Ma: "Robust face recognition via sparse representation." (2009)

のどあめ

ググってもなかなか出てこなかった情報を垂れ流すブログ。