Numo::NArrayをpocketfftでフーリエ変換するGemを作成した

はじめに

NumPyが、ver. 1.17から、フーリエ変換にFFTPACKではなくpocketfftを使う様になった。pocketfftはFFTPACKを修正したもので、速度と精度が改善されている。

「NumPy 1.17」リリース | OSDN Magazine

ENH: Add pocketfft sources to numpy for testing, benchmarks, etc. by mreineck · Pull Request #11888 · numpy/numpy · GitHub

https://gitlab.mpcdf.mpg.de/mtr/pocketfft

pocketfftは、単一のC言語のソースファイルからなるので「これはソースを同梱する形でNumo::NArrayな配列をフーリエ変換するなにかが作れるぞ」と思い、早速作ってみた。Numo::NArrayでは、フーリエ変換のAPIは提供されておらず、別でNumo::FFTWやNumo::FFTEをインストールする形になっている。そこで、Numo::Pocketfftとして公開した。

numo-pocketfft | RubyGems.org | your community gem host

使い方

pocketfftを同梱しているので、普通にgemコマンドでインストールできる。Mac、Ubuntu、Windows 10（RubyInstallerな環境）でインストールできることを確認した。

$ gem install numo-pocketfft

メソッド名や動きはnumpy.fftに従った。普通の高速フーリエ変換をfftメソッド、ifftメソッドで、実数な配列を想定した高速フーリエ変換をrfftメソッド、irfftメソッドで行う。

Module: Numo::Pocketfft — Documentation for numo-pocketfft (0.1.0)

まずは、実数の1次元配列のフーリエ変換・逆変換は以下のようになる。

irb(main):001:0> require 'numo/pocketfft'
=> true
irb(main):002:0> a = Numo::DFloat.new(4).rand
=> Numo::DFloat#shape=[4]
[0.0617545, 0.373067, 0.794815, 0.201042]
irb(main):003:0> x = Numo::Pocketfft.rfft(a)
=> Numo::DComplex#shape=[3]
[1.43068+0i, -0.733061-0.172025i, 0.282461+0i]
irb(main):004:0> Numo::Pocketfft.irfft(x)
=> Numo::DFloat#shape=[4]
[0.0617545, 0.373067, 0.794815, 0.201042]
irb(main):005:0>

逆変換の結果が複素数で返ってきても良いなら、全てfft系メソッドでもよい。例えば実数の2次元配列のフーリエ変換・逆変換は以下のようになる。fft2と2をつける（これはnumpy.fftにあわせた）。

irb(main):001:0> require 'numo/pocketfft'
=> true
irb(main):002:0> a = Numo::DFloat.new(4, 4).rand
=> Numo::DFloat#shape=[4,4]
[[0.0617545, 0.373067, 0.794815, 0.201042],
 [0.116041, 0.344032, 0.539948, 0.737815],
 [0.165089, 0.0508827, 0.108065, 0.0687079],
 [0.904121, 0.478644, 0.342969, 0.164541]]
irb(main):003:0> Numo::Pocketfft.ifft2(Numo::Pocketfft.fft2(a))
=> Numo::DComplex#shape=[4,4]
[[0.0617545+0i, 0.373067+0i, 0.794815+0i, 0.201042+0i],
 [0.116041+0i, 0.344032+0i, 0.539948+0i, 0.737815+0i],
 [0.165089+0i, 0.0508827+0i, 0.108065+0i, 0.0687079+0i],
 [0.904121+0i, 0.478644+0i, 0.342969+0i, 0.164541+0i]]

任意の次元数を受け付けるfftnとrfftnメソッドもある。

irb(main):001:0> require 'numo/pocketfft'
=> true
irb(main):002:0> a = Numo::DFloat.new(2,2,2).rand + Complex::I * Numo::DFloat.new(2,2,2).rand
=> Numo::DComplex#shape=[2,2,2]
[[[0.0617545+0.165089i, 0.373067+0.0508827i],
  [0.794815+0.108065i, 0.201042+0.0687079i]],
 [[0.116041+0.904121i, 0.344032+0.478644i],
  [0.539948+0.342969i, 0.737815+0.164541i]]]
irb(main):003:0> Numo::Pocketfft.ifftn(Numo::Pocketfft.fftn(a))
=> Numo::DComplex#shape=[2,2,2]
[[[0.0617545+0.165089i, 0.373067+0.0508827i],
  [0.794815+0.108065i, 0.201042+0.0687079i]],
 [[0.116041+0.904121i, 0.344032+0.478644i],
  [0.539948+0.342969i, 0.737815+0.164541i]]]

といった感じで、フーリエ変換・逆変換できる。

git submodule による外部ライブラリの同梱

pocketfftの同梱は、Gitlab上のpocketfftのリポジトリを、submoduleとして取り込むことで実現した。submoduleの追加は、特別な工夫はなくgitコマンド通りになる。

$ cd ext/numo/pocketfft
$ git submodule add https://gitlab.mpcdf.mpg.de/mtr/pocketfft pocketfft

Ruby extensionでは、mkmfというライブラリでMakefileを作成する。mkmfでは、基本的には、全てのコードが同一のディレクトリ下（Numo::Pocketfftを例にすればext/numo/pocketfft下）にあることを想定している。submoduleで追加したpocketfftは、サブディレクトリにある。この階層構造を、mkmfに教えるようなことが必要になる。

# extconf.rb
# ...省略

# ext/numo/pocketfft下のコードを$srcsに追加する
$srcs = Dir.glob("#{$srcdir}/*.c").map { |path| File.basename(path) }

# pocketfftのコードを$srcsに追加する
$srcs << 'pocketfft.c'

# ext/numo/pocketfft/pocketfftをインクルードパスやmakeの探索パスに追加する
# ※ インクルードパスはいらないかも
Dir.glob("#{$srcdir}/*/") do |path|
  dir = File.basename(path)
  $INCFLAGS << " -I$(srcdir)/#{dir}"
  $VPATH << "$(srcdir)/#{dir}"
end

# ...省略

これで、rake compileでpocketfftも含めてコンパイルできる。そして、gemファイルを作る際に、submoduleなpocketfftを含めるために、gemspecファイルで、Gem::Specificationのfilesにsubmoduleへのパスを追加する。

# numo-pocketfft.gemspec
# ...省略

Gem::Specification.new do |spec|
# ...省略

  # submoduleのディレクトリ下のファイルをfilesに追加する
  # 参考: https://gist.github.com/mattconnolly/5875987
  gem_dir = __dir__ + '/'
  `git submodule --quiet foreach pwd`.split($OUTPUT_RECORD_SEPARATOR).each do |submodule_path|
    Dir.chdir(submodule_path) do
      submodule_relative_path = submodule_path.sub gem_dir, ''
      `git ls-files`.split($OUTPUT_RECORD_SEPARATOR).each do |filename|
        spec.files << "#{submodule_relative_path}/#{filename}"
      end
    end
  end

# ...省略

これでsubmoduleにより外部ライブリをgemに同梱できた。後で発見したが、同様のことがruby-eigenでも行われていた。ruby-eigenでは、C++の線形代数ライブラリEigenが、submoduleにより同梱されている。

おわりに

画像にフィルタを掛ける場合、画像とフィルタとの畳み込み演算なる。畳み込み演算は、フーリエ変換により複素数の世界に持っていくと、シンプルな乗算になる。これを利用すると、畳み込み演算をフーリエ変換・逆変換で実装できる。画像処理ライブラリのMagroの開発で、このフーリエ変換による畳み込みを必要としていたので、良い機会だと思いpocketfftによるフーリエ変換のgemを作った。NumPyのフーリエ変換では、正規化のオプションがあったりするが、Numo::Pocketfftでは、まずは、シンプルなフーリエ変換・逆変換を実装した。今後bugfixとかも含めて、バージョンアップしていくつもりでいる。

github.com

画像をNumo::NArrayで表現するRubyの画像処理ライブラリを作り始めた

はじめに

PythonではOpenCVとScikit-imageという画像処理ライブラリがある。これらライブラリでは、画像データをNumPyのndarrayで表す。

>>> from skimage import io
>>> img = io.imread('lena.png')
>>> type(img)
<class 'numpy.ndarray'>

同じような感じで、画像データをNumo::NArrayで扱うものが欲しくなったので作り始めた。名前は、Image ProcessingからとってMagroとした（Gem名を考えているときに、ふと「そういえば鮪ペイントって昔あったよな」と思い出したのもあって）。実はまだ、pngとjpegを読み書きすることしかできないが...

magro | RubyGems.org | your community gem host

インストール

Magroでは、pngとjpegの読み書きにlibpngとlibjpegを使用するため、これらを先にインストールする。メジャーなライブラリなので、なにかしらパッケージがある（環境によってはすでにインストールされている場合もある）。

macOSであればbrewで:

$ brew install libpng libjpeg

Ubuntuなどlinux系であれば適当なパッケージマネージャで:

$ sudo apt-get install libpng-dev libjpeg-dev

Magroは、gemコマンドでインストールできる。

$ gem install magro

使い方

Magroでは、まだ、pngやjpeg画像の読み書きしかできない。画像の形式はファイル名の拡張子で判定している。例えば、有名なRGBなPNG画像ファイルをグレイスケールにするならば、以下のようになる。

[1] pry(main)> require 'magro'
=> true
[2] pry(main)> img = Magro::IO.imread('lena.png')
=> Numo::UInt8#shape=[512,512,3]
[[[226, 137, 125],
  [226, 137, 125],
  [223, 137, 133],
...
[3] pry(main)> gray = img.median(axis: 2)
=> Numo::UInt8#shape=[512,512]
[[137, 137, 137, 136, 138, 129, 138, 134, 140, 136, 135, 134, 130, 139, ...],
 [137, 137, 137, 136, 138, 129, 138, 134, 140, 136, 135, 134, 130, 139, ...],
 [137, 137, 137, 136, 138, 129, 138, 134, 140, 136, 135, 134, 130, 139, ...],
...
[4] pry(main)> Magro::IO.imsave('lena_gray.png', gray)
=> true

画像データはNumo::UInt8で表現されるので、Red ChainerのConvolution Netsに投げるとか、Numo::DFloatなベクトルに変換してRumaleのRandom Forestに投げるとかできる。

おわりに

Rubyの画像処理ライブラリとしては、RMagickをはじめImageMagickをbindingしたものや、OpenCVをbindingしたものがある。OpenCVでは、画像をCvMatというOpenCVの行列型で表現する。行列の表現がNumo::NArrayでなくてもOKで、画像を行列で欲しい場合は、RubyであればOpenCVのbindingが良い（ただし対応しているOpenCVのバージョンは2系の様子）。

Magroは、PythonのScikit-Imageや、C言語のVLFeatあたりの雰囲気を目指して、あまり大きくならない程度に開発する予定。次は、画像にフィルタをかけるには、畳み込み演算が必要で、そのためのフーリエ変換をなにかしら組み込むとかかなぁ。

github.com

Numo::NArrayでLIBLINEARを使うGemを作成した

はじめに

Numo::Libsvmと同様に、特徴ベクトルやラベルをNumo::NArray形式で扱えるLIBLINEARのbinding gemを作成した。LIBLINEARは、SVMやロジスティック回帰による線形の分類器や回帰が実装されている。

numo-liblinear | RubyGems.org | your community gem host

使い方

準備

Numo::Liblinearでは、LIBLINEAR自体を同梱していないので、別途LIBLINEARをインストールする必要がある（必要になるのはliblinear.soとlinear.hである）。

例えば、macOSであれば、

$ brew install liblinear

Ubuntuであれば、

$ sudo apt-get install liblinear-dev

となる。

インストール

numo-liblinearは、gemコマンドでインストールできる。

$ gem install numo-liblinear

L2正則化ロジスティック回帰による分類

LIBSVM DataにあるPendigitsデータセットを使って、ロジスティック回帰による分類を行う。データの取得にはred-datasetsを用いる。

$ gem install red-datasets-numo-narray 

まず、ロジスティック回帰による分類器を訓練する。

require 'numo/liblinear'
require 'datasets-numo-narray'

# Pendigitsデータの訓練データセットをダウンロードする.
puts 'Download dataset.'
pendigits = Datasets::LIBSVM.new('pendigits').to_narray
x = pendigits[true, 1..-1]
y = pendigits[true, 0]

# LIBLINEARのパラメータを定義する.
# solver_typeでL2正則化ロジスティック回帰を示す.
param = {
  solver_type: Numo::Liblinear::SolverType::L2R_LR_DUAL
}

# L2正則化ロジスティック回帰を訓練する.
puts 'Train logistic regression.'
model = Numo::Liblinear.train(x, y, param)

# パラメータと訓練したモデルをMarshalでファイルに保存する.
puts 'Save parameters and model with Marshal'
File.open('pendigits.dat', 'wb') { |f| f.write(Marshal.dump([param, model])) }

これを実行すると、以下のようになる。

$ ruby train.rb
Download dataset.
Train logistic regression.
Save parameters and model with Marshal

次に、訓練したモデルで、テストデータの分類を行う。

require 'numo/liblinear'
require 'datasets-numo-narray'

# Pendigitsデータのテストデータセットをダウンロードする.
puts 'Download dataset.'
pendigits_test = Datasets::LIBSVM.new('pendigits', note: 'testing').to_narray
x = pendigits_test[true, 1..-1]
y = pendigits_test[true, 0]

# パラメータと訓練したモデルをMarshalで読み込む.
puts 'Load parameter and model.'
param, model = Marshal.load(File.binread('pendigits.dat'))

# テストデータのラベルを推定する.
puts 'Predict labels.'
predicted = Numo::Liblinear.predict(x, param, model)

# 推定結果を正確度で評価する.
mean_accuracy = y.eq(predicted).count.fdiv(y.size)
puts "Accuracy: %.1f %%" % (100 * mean_accuracy)

これを実行すると、以下のようになる。正確度（Accuracy）が87.9%となり、線形分類としてはまずまずこんなモノかなといった感じ。

$ ruby test.rb
Download dataset.
Load parameter and model.
Predict labels.
Accuracy: 87.9 %

おわりに

Numo::Libsvmと同様、ひとまず、最低限の動きができるようになったので公開した。また、これもNumo::Libsvmと同様で、LIBLINEAR自体を同梱していないので、RubyInstallerによるWindows環境では動作しないかもしれない（liblinearライブラリとlinear.hがあればWindowsでも動くと思われる）。ひとまず、広く使われている機械学習関連ライブラリの薄いラッパーを用意して、リッチなインターフェースをRumaleかなにかで提供できれば、と考えている。次は、なにか近似最近傍探索ライブリかな〜。

github.com

データセットをNumo::NArrayであつかうLIBSVMのGemを作成した

はじめに

特徴ベクトルやラベルをNumo::NArray形式で扱えるLIBSVMのbinding gemが欲しかったので、Ruby拡張ライブラリの勉強もかねて作成した。開発に際しては、Numo::FFTWをとても参考にさせて頂いた。gem名は、勝手ながらnumo-libsvmとした。

numo-libsvm | RubyGems.org | your community gem host

使い方

準備

多くのLIBSVM関連ライブラリがLIBSVMのコード自体を同梱しているが、Numo::Libsvmでは同梱しないことにした。brewやaptなどのパッケージマネージャでインストールしたLIBSVMと、バージョンがズレるのはどうかなと思ったためである。というわけで、LIBSVMをインストールする必要がある（必要になるのはlibsvm.soとsvm.hである）。

例えば、macOSであれば、

$ brew install libsvm

Ubuntuであれば、

$ sudo apt-get install libsvm-dev

となる。

インストール

numo-libsvmは、gemコマンドでインストールできる。

$ gem install numo-libsvm

C-SVCによる分類

LIBSVM DataにあるPendigitsデータセットを使って、C-SVCによる分類を行う。データの取得にはred-datasetsを用いる。

$ gem install red-datasets-numo-narray 

では、まず、C-SVCによる分類器を訓練する。

require 'numo/narray'
require 'numo/libsvm'
require 'datasets-numo-narray'

# Pendigitsデータの訓練データセットをダウンロードする.
puts 'Download dataset.'
pendigits = Datasets::LIBSVM.new('pendigits').to_narray
x = pendigits[true, 1..-1] # 特徴ベクトル
y = pendigits[true, 0]     # ラベル

# RBFカーネルによるC-SVCを実現するパラメータを定義する.
param = {
  svm_type: Numo::Libsvm::SvmType::C_SVC,
  kernel_type: Numo::Libsvm::KernelType::RBF,
  gamma: 0.0001, # RBFカーネルのパラメータ
  C: 10, # C-SVCのパラメータ
  shrinking: true
}

# C-SVCを訓練する.
puts 'Train support vector machine.'
model = Numo::Libsvm.train(x, y, param)

# パラメータと訓練したモデルをMarshalでファイルに保存する。
puts 'Save parameters and model with Marshal.'
File.open('pendigits.dat', 'wb') { |f| f.write(Marshal.dump([param, model])) }

これを実行すると、以下のようになる。

$ ruby train.rb
Download dataset.
Train support vector machine.
Save paramters and model with Marshal.

次に、訓練したモデルで、テストデータの分類を行う。

require 'numo/narray'
require 'numo/libsvm'
require 'datasets-numo-narray'

# Pendigitsデータのテストデータセットをダウンロードする.
puts 'Download dataset.'
pendigits_test = Datasets::LIBSVM.new('pendigits', note: 'testing').to_narray
x = pendigits_test[true, 1..-1]
y = pendigits_test[true, 0]

# パラメータと訓練したモデルをMarshalで読み込む.
puts 'Load parameter and model.'
param, model = Marshal.load(File.binread('pendigits.dat'))

# テストデータのラベルを推定する.
puts 'Predict labels.'
predicted = Numo::Libsvm.predict(x, param, model)

# 推定結果を正確度で評価する.
mean_accuracy = y.eq(predicted).count.fdiv(y.size)
puts "Accuracy: %.1f %%" % (100 * mean_accuracy)

これを実行すると、以下のようになる。正確度（Accuracy）が98.3%となり、上手く分類できていることがわかる。

$ ruby test.rb
Download dataset.
Load parameter and model.
Predict labels.
Accuracy: 98.3 %

おわりに

ひとまず、最低限の動きができるようになったので公開した。これから、ドキュメントの整備も含めてアップデートしていきたい。また、LIBSVMのコードを同梱していないので、RubyInstallerによるWindows環境では動作しないかもしれない（libsvmのライブラリとヘッダーが指定できればWindowsでも動くと思われる）。

Rumaleなインターフェースも用意しようかと考えたが、それは今後、例えばrumale-libsvmのような形で、Rumale側でできればと思う。

github.com

Rubyのパイプラインを考えてみた

はじめに

Rubyにパイプライン演算子が追加されるかも？というのが、少し前に話題になった。他方、R言語にはmagrittrというパイプラインを実現するパッケージがあり、これが便利だ。

# h(g(f(x)))
x %>% f %>% g %>% h

データ分析では、カテゴリ変数を整数値化して欠損値を補ってアレしてコレして〜など、データを段階的に変換する処理がよく発生する。そういった処理が、magrittrでスッキリ書ける。これがRubyにも欲しい。

やってみた

パイプラインを実現するために、ひとまず「自らを第一引数としてメソッドを呼び出す」ことができれば良さそうだ（magrittrではドットで任意の引数に入れることができるが、今回は見送る）。 これは、以下のようにシンプルに書ける。受け取ったメソッドfnを、自らを第一引数としてsendで呼び出す。メソッド名は、絵文字でそれっぽいのにした。

module Pipeline
  def 👉(fn, *args)
    send(fn, self, *args)
  end
end

これをStringで使ってみる。scatメソッドは、スペースをはさんで文字列を連結する。

class String
  include Pipeline
end

def scat(a, b)
  "#{a} #{b}"
end

"Hello"
  .👉(:scat, "world")
  .👉(:scat, "!!")
  .👉(:puts)

これを実行すると、次のようになる。Helloにworldが連結され、それにさらに!!が連結され、putsで出力されている。

Hello world !!

数値もいける。

class Numeric
  include Pipeline
end

def add(a, b)
  a + b
end

1
  .👉(:add, 2)
  .👉(:add, 0.5)
  .👉(:puts)

# 1 + 2 + 0.5 で 3.5 が表示される

もうKernelにぶっこんじゃえばいいかも。

module Kernel
  def 👉(fn, *args)
    send(fn, self, *args)
  end
end

メソッド名の指定にシンボルを使っているのと、👉の手前にドットがあるのが、演算子っぽくない。TracePointとか黒魔術を試みたけど、なんか上手くいかなかった。やってるうちに、そもそも演算子じゃないからコレで良いかな、という心になってしまった。途中に、パイプラインと関係ないメソッドを挟んでも、自然に見えるので。

"Hello"
  .👉(:scat, "world")
  .👉(:scat, "!!")
  .upcase
  .👉(:puts)

# HELLO WORLD !! と表示される

おわりに

こだわれば、もっと複雑で高度なモノにできると思う。ひとまず、データの前処理がmagrittrぽく書けるからOK 👌