CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: NumPy 配列の操作でメモリのコピーが生じているか調べる

Python NumPy

パフォーマンスの観点からいえば、データをコピーする機会は少ないほど望ましい。 コンピュータのバスの帯域幅は有限なので、データをコピーするには時間がかかる。

NumPy の配列 (ndarray) には、メモリを実際に確保している配列と、それをただ参照しているだけのビュー (view) がある。 そして、配列への操作によって、メモリが確保されて新しい配列が作られるか、それとも単なるビューになるかは異なる。 今回は NumPy の配列を操作するときにメモリのコピーが生じているか調べる方法について。

使った環境は次のとおり。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.15.5
BuildVersion:   19F101
$ python -V                             
Python 3.7.7
$ pip list | grep -i numpy
numpy              1.19.0

下準備

あらかじめ NumPy をインストールしておく。

$ pip install numpy

Python のインタプリタを起動する。

$ python

サンプルとなる配列を用意する。

>>> import numpy as np
>>> a = np.arange(10)

flags を使った調べ方

はじめに ndarray#flags を使った調べ方から。 NumPy の配列には flags というアトリビュートがあって、ここから配列の情報がいくらか得られる。 この中に owndata という情報があって、これはオブジェクトのメモリが自身のものか、それとも別のオブジェクトを参照しているかを表す。

numpy.org

最初に作った配列については、このフラグが True にセットされている。

>>> a.flags.owndata
True

では、上記に対してスライスを使って配列の一部を取り出した場合にはどうだろうか。

>>> b = a[1:]

スライスで取り出した配列の場合、フラグは False にセットされている。

>>> b.flags.owndata
False

つまり、自身でメモリを確保しているのではなく、別のオブジェクトを参照しているだけ。

ndarray#base を使った調べ方

同じように ndarray#base を使って調べることもできそうだ。 このアトリビュートは、オブジェクトが別のオブジェクトのメモリに由来している場合に、そのオブジェクトへの参照が入る。

numpy.org

先ほどの例では、最初に作った配列は None になっている。

>>> a.base is None
True

一方で、スライスを使って取り出した配列は、元になった配列への参照が入っている。

>>> b.base is a
True

インプレース演算

ところで、インプレース演算の場合は ndarray#flags や ndarray#base を使った判定ができないのかな、と思った。

たとえば配列を使った通常の加算 (__add__()) では、新しく配列が作られてメモリのコピーが生じる。

>>> c = a + 1
>>> c.flags.owndata
True
>>> c.base is None
True

一方で、インプレースの加算 (__iadd__()) を使ったときも、これまで紹介してきたアトリビュートは同じ見え方になる。

>>> a += 1
>>> a.flags.owndata
True
>>> a.base is None
True

では、メモリのコピーは生じているかというと生じていない。

NumPy の配列には __array_interface__ というアトリビュートがある。 その中にある data というキーからは、配列の最初の要素が格納されているメモリのアドレス情報が得られる。

>>> a.__array_interface__['data']
(140489713031584, False)

以下のように、インプレース演算をしてもアドレス情報に変化はない。

>>> a += 1
>>> a.__array_interface__['data']
(140489713031584, False)

つまり、新たにメモリは確保されていない。

numpy.org

なお、それ以外にも大きな配列を用意してベンチマークしたり、ソースコードを読んで調べることも考えられる。

NumPyによるデータ分析入門 ―配列操作、線形代数、機械学習のためのPythonプログラミング

NumPyによるデータ分析入門 ―配列操作、線形代数、機械学習のためのPythonプログラミング

ltrace(1) で共有ライブラリの呼び出しを追いかける

Linux Ubuntu18.04LTS

Linux システムでは ltrace(1) を使うことで共有ライブラリの呼び出しを調べることができる。 今回は、いくつかの例を用いて使い方についてざっくりと見ていく。

使った環境は次のとおり。

$ cat /etc/lsb-release 
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=18.04
DISTRIB_CODENAME=bionic
DISTRIB_DESCRIPTION="Ubuntu 18.04.4 LTS"
$ uname -r
4.15.0-111-generic
$ ltrace -V
ltrace version 0.7.3.
Copyright (C) 1997-2009 Juan Cespedes <cespedes@debian.org>.
This is free software; see the GNU General Public Licence
version 2 or later for copying conditions.  There is NO warranty.

もくじ

下準備

あらかじめ、利用するパッケージをインストールしておく。 尚、ltrace 以外は使い方のサンプルとして取り上げているだけ。

$ sudo apt-get -y install ltrace libc-bin gcc wget python3-requests python3-numpy

自作のプログラムを使った例

まずは最も単純な例として、ハローワールドするだけのプログラムを試してみよう。

次のような C 言語のソースコードを用意する。

$ cat << 'EOF' > greet.c
#include <stdio.h>
#include <stdlib.h>


int main(void) {
    printf("Hello, World!\n");
    return EXIT_SUCCESS;
}
EOF

上記をビルドする。

$ gcc -Wall -o greet greet.c

これで ELF フォーマットの実行可能ファイルがえきる。

$ file greet
greet: ELF 64-bit LSB shared object, x86-64, version 1 (SYSV), dynamically linked, interpreter /lib64/ld-linux-x86-64.so.2, for GNU/Linux 3.2.0, BuildID[sha1]=1771a2c7c6ded094f15254590870089080689968, not stripped

次のとおり、実行するとただメッセージを出力して終了する。

$ ./greet 
Hello, World!

このように、ただハローワールドするだけのプログラムでも標準 C ライブラリ (libc) はダイナミックリンクされている。

$ ldd greet
    linux-vdso.so.1 (0x00007ffd31d9e000)
    libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f7081218000)
    /lib64/ld-linux-x86-64.so.2 (0x00007f708180b000)

このプログラムを ltrace(1) 経由で実行してみよう。 すると、共有ライブラリの呼び出しに関する情報が出力される。

$ ltrace ./greet > /dev/null
puts("Hello, World!")                                                                                          = 14
+++ exited (status 0) +++

puts(3) は標準 C ライブラリの提供している API のひとつで、詳細は man を参照のこと。

$ man 3 puts

たとえば -l オプションを指定すると特定の共有ライブラリの呼び出しだけを追跡できる。 ただし、上記のサンプルでは libc の API しか呼び出していないので関係ない。

$ ltrace -l libc.so.6 ./greet > /dev/null
greet->puts("Hello, World!")                                                                                   = 14
+++ exited (status 0) +++

echo(1) を使った例

次は echo(1) を使ってみよう。 こちらも標準 C ライブラリがダイナミックリンクされている。

$ echo "Hello, World"
Hello, World
$ ldd $(which echo)
    linux-vdso.so.1 (0x00007ffe8dddf000)
    libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f87d2478000)
    /lib64/ld-linux-x86-64.so.2 (0x00007f87d2a72000)

ltrace(1) 経由で実行すると、先ほどよりも色々な API が呼び出されていることがわかる。

$ ltrace echo "Hello, World"
getenv("POSIXLY_CORRECT")                                                                                      = nil
strrchr("echo", '/')                                                                                           = nil
setlocale(LC_ALL, "")                                                                                          = "C.UTF-8"

...

__freading(0x7fc11d6c3680, 0, 4, 2880)                                                                         = 0
fflush(0x7fc11d6c3680)                                                                                         = 0
fclose(0x7fc11d6c3680)                                                                                         = 0
+++ exited (status 0) +++

wget(1) を使った例 (LibSSL)

もうちょっと複雑な例として wget(1) を試してみよう。

ただし、次は libc ではなく libssl の呼び出しを追跡する。 次のように wget(1) は libssl をダイナミックリンクしている。

$ ldd $(which wget) | grep ssl
    libssl.so.1.1 => /usr/lib/x86_64-linux-gnu/libssl.so.1.1 (0x00007f85b9051000)

libssl の呼び出しに限定して wget(1) を ltrace(1) 経由で実行してみよう。 たとえば Google のウェブサイトを取得させてみる。

$ ltrace -l libssl.so.1.1 wget -qO /dev/null https://google.co.jp
wget->OPENSSL_init_ssl(0, 0, 0x55af441e1ec0, 0)                                                                = 1
wget->OPENSSL_init_ssl(0x200002, 0, 0x7fffffff, 0)                                                             = 1
wget->TLS_client_method(0x7f9b550c891c, 129, 0, 0x55af42583904)                                                = 0x7f9b550c28e0

...

wget->SSL_pending(0x55af441efec0, 0x55af44201a30, 127, 0x55af441f5270)                                         = 2
wget->SSL_peek(0x55af441efec0, 0x55af44201a30, 127, 1)                                                         = 2
wget->SSL_read(0x55af441efec0, 0x55af44201a30, 2, 0x7f9b54447eb7)                                              = 2
+++ exited (status 0) +++

ちゃんと呼び出しがトレースできていることがわかる。

Python を使った例 (LibSSL)

次は Python のインタプリタについて、同じように libssl の呼び出しを追いかけてみる。 注意すべき点として、このような例では Python のバイナリは直接は LibSSL をリンクしていない。

$ ldd $(which python3)
    linux-vdso.so.1 (0x00007ffca4dd4000)
    libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f78ac216000)
    libpthread.so.0 => /lib/x86_64-linux-gnu/libpthread.so.0 (0x00007f78abff7000)
    libdl.so.2 => /lib/x86_64-linux-gnu/libdl.so.2 (0x00007f78abdf3000)
    libutil.so.1 => /lib/x86_64-linux-gnu/libutil.so.1 (0x00007f78abbf0000)
    libexpat.so.1 => /lib/x86_64-linux-gnu/libexpat.so.1 (0x00007f78ab9be000)
    libz.so.1 => /lib/x86_64-linux-gnu/libz.so.1 (0x00007f78ab7a1000)
    libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007f78ab403000)
    /lib64/ld-linux-x86-64.so.2 (0x00007f78ac607000)

代わりに、同梱されている共有ライブラリが間接的にリンクしている。

$ python3 -c "import _ssl; print(_ssl.__file__)"
/usr/lib/python3.6/lib-dynload/_ssl.cpython-36m-x86_64-linux-gnu.so
$ ldd /usr/lib/python3.6/lib-dynload/_ssl.cpython-36m-x86_64-linux-gnu.so | grep ssl
    libssl.so.1.1 => /usr/lib/x86_64-linux-gnu/libssl.so.1.1 (0x00007f181f27f000)

このような状況では ltrace(1) がブレークポイントを仕掛けるべき共有ライブラリを自動では認識できない。 そのため -l オプションで明示的に対象の共有ライブラリを指定しないと呼び出しが表示されないようだ。

$ ltrace -l libssl.so.1.1 python3 -c "import requests; requests.get('https://google.co.jp')"
--- SIGCHLD (Child exited) ---
--- SIGCHLD (Child exited) ---
--- SIGCHLD (Child exited) ---
_ssl.cpython-36m-x86_64-linux-gnu.so->TLS_method(0xa74560, 0, 0, 0)                                            = 0x7f2b498dd1a0
_ssl.cpython-36m-x86_64-linux-gnu.so->SSL_CTX_new(0x7f2b498dd1a0, 0, 0, 0)                                     = 0x25de450
_ssl.cpython-36m-x86_64-linux-gnu.so->SSL_CTX_get_verify_callback(0x25de450, 0x7f2b475ad5bc, 0, 0x7f2b498f8718) = 0

...

_ssl.cpython-36m-x86_64-linux-gnu.so->SSL_CTX_free(0x25de450, 0x9d2800, -3, 0x27f0070)                         = 0
_ssl.cpython-36m-x86_64-linux-gnu.so->SSL_free(0x27dbaf0, 0x7f2b475d15a8, -2, 2)                               = 0
_ssl.cpython-36m-x86_64-linux-gnu.so->SSL_CTX_free(0x22e8060, 0x9d2800, -3, 0x27f0070)                         = 0
+++ exited (status 0) +++

Python / NumPy を使った例 (LibBLAS)

最後に、おまけとして Python の NumPy が依存している libblas の呼び出しを追いかけてみる。 libblas の実装は Reference BLAS や ATLAS、OpenBLAS、Intel MKL など色々とある。 今回は apt-get(8) を使って NumPy をインストールして Reference BLAS が実装の例になる。

先ほどと同じように、libblas は NumPy に同梱されている共有ライブラリが間接的にリンクしている。 NumPy には Python/C API で書かれた、次のような共有ライブラリがある。

$ dpkg -L python3-numpy | grep so$ | grep -v test
/usr/lib/python3/dist-packages/numpy/core/_dummy.cpython-36m-x86_64-linux-gnu.so
/usr/lib/python3/dist-packages/numpy/core/multiarray.cpython-36m-x86_64-linux-gnu.so
/usr/lib/python3/dist-packages/numpy/core/umath.cpython-36m-x86_64-linux-gnu.so
/usr/lib/python3/dist-packages/numpy/fft/fftpack_lite.cpython-36m-x86_64-linux-gnu.so
/usr/lib/python3/dist-packages/numpy/linalg/_umath_linalg.cpython-36m-x86_64-linux-gnu.so
/usr/lib/python3/dist-packages/numpy/linalg/lapack_lite.cpython-36m-x86_64-linux-gnu.so
/usr/lib/python3/dist-packages/numpy/random/mtrand.cpython-36m-x86_64-linux-gnu.so

この中で ndarray を実装しているのが multiarray で、次のように libblas をリンクしている。

$ ldd /usr/lib/python3/dist-packages/numpy/core/multiarray.cpython-36m-x86_64-linux-gnu.so | grep blas
    libblas.so.3 => /usr/lib/x86_64-linux-gnu/libblas.so.3 (0x00007faf3e37a000)

試しに NumPy を使って配列のドット積を計算してみよう。 すると、内部的に LibBLAS の API を呼び出していることが確認できる。

$ ltrace -l libblas.so.3 python3 -c "import numpy as np; x = np.random.randn(100, 100); np.dot(x, x)"
multiarray.cpython-36m-x86_64-linux-gnu.so->cblas_dgemm(101, 111, 111, 100 <unfinished ...>
libblas.so.3->dgemm_(0x7ffc75477dd6, 0x7ffc75477dd7, 0x7ffc75477dc8, 0x7ffc75477dcc <unfinished ...>
libblas.so.3->lsame_(0x7ffc75477dd6, 0x7fd337a21140, 1, 1)                                                                  = 1
libblas.so.3->lsame_(0x7ffc75477dd7, 0x7fd337a21140, 1, 1)                                                                  = 1
<... dgemm_ resumed> )                                                                                                      = 0
<... cblas_dgemm resumed> )                                                                                                 = 0
+++ exited (status 0) +++

いじょう。

Linuxで動かしながら学ぶTCP/IPネットワーク入門

Linuxで動かしながら学ぶTCP/IPネットワーク入門

Python: LightGBM 開発環境メモ

LightGBM Python

最近 LightGBM にコントリビューションする機会を得たので、その際に調べたことの備忘録を残しておく。 現時点では、Python 周りの開発環境についてドキュメントは特に見当たらないようだった。 以下は CI 環境のスクリプトやエラーメッセージを読みながら雰囲気で作ったもの。 そのため、あまり信用せずオフィシャルなところは本家のリポジトリで確認してほしい。

使った環境は次のとおり。

$ sw_vers                            
ProductName:    Mac OS X
ProductVersion: 10.15.5
BuildVersion:   19F101
$ python -V                      
Python 3.7.7

下準備

拡張モジュールやドキュメントのビルドに必要なパッケージをあらかじめインストールしておく。

$ brew install cmake doxygen

そして、ソースコードのリポジトリをチェックアウトする。

$ git clone https://github.com/microsoft/LightGBM.git
$ cd LightGBM

Python のテストを通す

今回は手を入れるのが Python binding 部分だったので、ひとまず Python のユニットテストを通す必要がある。

まずはテストに必要なパッケージをインストールする。

$ pip install pandas psutil scipy

続いて、Python binding のあるディレクトリに移動する。

$ cd LightGBM/python-package

LightGBM の本体といえる拡張モジュール (lib_lightgbm.so) をビルドした上で、開発モードで LightGBM をインストールする。

$ python setup.py bdist develop

ユニットテストのあるディレクトリに移動する。

$ cd ../tests/python_package_test

テストランナーを実行して、ユニットテストがパスすることを確認する。

$ python -m unittest discover -v
...
test_stacking_regressor (test_sklearn.TestSklearn) ... ok
test_xendcg (test_sklearn.TestSklearn) ... ok

----------------------------------------------------------------------
Ran 107 tests in 146.355s

OK (skipped=5)

ドキュメントをビルドする

続いて、API に修正があった場合にはドキュメントにも手を入れる必要がある。 LightGBM では Sphinx を使ってドキュメントを書いている。

まずはドキュメントのディレクトリに移動する。

$ cd ../../docs/

ドキュメントをビルドするのに必要なパッケージをインストールする。

$ pip install -r requirements.txt

あとは通常の Sphinx の手順でターゲットを指定してドキュメントをビルドする。

$ make html

うまくいけば _build ディレクトリ以下に成果物ができるので修正した内容を確認する。

$ ls _build/html   
Advanced-Topics.html        Python-Intro.html
C-API.html          Quick-Start.html
Development-Guide.html      README.html
Experiments.html        _images
FAQ.html            _modules
Features.html           _sources
GPU-Performance.html        _static
GPU-Targets.html        gcc-Tips.html
GPU-Tutorial.html       genindex.html
GPU-Windows.html        index.html
Installation-Guide.html     objects.inv
Parallel-Learning-Guide.html    pythonapi
Parameters-Tuning.html      search.html
Parameters.html         searchindex.js
Python-API.html

いじょう。

Kaggleで勝つデータ分析の技術

Kaggleで勝つデータ分析の技術

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＰｙｔｈｏｎではじめるＫａｇｇｌｅスタートブック (ＫＳ情報科学専門書)

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＰｙｔｈｏｎではじめるＫａｇｇｌｅスタートブック (ＫＳ情報科学専門書)

xargs(1) でシェル関数を使いたい

ShellScript

コマンドラインの処理を並列実行したいときなどに使う xargs(1) だけど、引数にシェル関数を使おうとすると少し工夫する必要がある。 工夫しない場合に失敗する理由から説明しているので、うまくいくやり方だけ知りたいときは下までスクロールしてもらえると。

使った環境は次のとおり。 シェルとしては bash を想定する。 なお、xargs(1) は GNU 版か BSD 版かは問わない。

$ sw_vers       
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G5033
$ bash --version 
GNU bash, version 3.2.57(1)-release (x86_64-apple-darwin18)
Copyright (C) 2007 Free Software Foundation, Inc.

とりあえず失敗させてみる

たとえば、次のように greet() という名前でシェル関数を定義しておく。

$ greet() { echo "Hello, $1"; }
$ greet "World"
Hello, World

そして、特に何も考えず上記で定義した関数を xargs(1) 経由で実行しようとすると、次のようにエラーとなる。

$ echo "World" | xargs -I {} greet {}
xargs: greet {}: No such file or directory

コマンドが失敗する理由について

上記でコマンドが失敗する理由は複数ある。 まず、xargs(1) の基本的な動作原理は、プロセスを fork(2) して exec*(2) することに相当する 1 。 この動作原理から、少なくとも xargs(1) の引数は実行可能ファイルが先頭に来る必要があることがわかる。 シェル関数はあくまでシェルの中で利用できるものなので、まずはここが問題となる。

では、次のように bash をインラインで実行すればどうだろうか。 だいぶ確信には迫っているものの、まだ足りないのでエラーになる。

$ echo "World" | xargs -I {} bash -c "greet {}"
bash: greet: command not found

原因の 2 つ目は、xargs(1) 経由で実行したシェルの中でシェル関数が有効ではない点。 そのため、次のように export することでサブプロセスのシェルでもシェル関数が有効であるようにしなければいけない。

$ export -f greet

ようするに、xargs(1) は関係なくインラインで bash を実行したときに、ちゃんとシェル関数が使えるようになっている必要がある。

$ bash -c "greet 'World'"
Hello, World

うまくいくやり方

ということで、うまくいくやり方は次のとおり。

$ greet() { echo "Hello, $1"; }
$ export -f greet
$ echo "World" | xargs -I {} bash -c "greet {}"
Hello, World

いじょう。

Linuxで動かしながら学ぶTCP/IPネットワーク入門

Linuxで動かしながら学ぶTCP/IPネットワーク入門

  1. 少なくとも GNU 版の実装はそうなっているようだった↩

Python: 画像データをフーリエ変換して周波数領域で扱ってみる

Digital Signal Processing Python NumPy Matplotlib macOS

フーリエ変換は音声データに対して用いられることが多い手法だけど、画像データにも応用が効く。 音声データの場合、フーリエ変換を使うことで時間領域の情報を周波数領域の情報に直せる。 それに対し、画像データでは空間領域の情報を周波数領域の情報に直すことになる。 つまり、画像データの濃淡が複数の波形の合成によって作られていると見なす。 今回は、画像データをフーリエ変換して、周波数領域の情報にフィルタをかけたり元に戻したりして遊んでみる。

使った環境は次のとおり。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G5033
$ python -V        
Python 3.7.7

もくじ

下準備

あらかじめ、必要なパッケージをインストールしておく。

$ pip install pillow matplotlib

周波数領域の情報を可視化してみる

最初に、画像データをフーリエ変換して周波数領域で可視化してみる。 画像データは、二次元の NumPy 配列として読み込んだ上で np.fft.fft2() に渡すことでフーリエ変換できる。 また、周波数領域のデータは np.fft.ifft2() を使うことで空間領域のデータに戻せる。

以下のサンプルコードでは、元の画像と、周波数領域でのパワースペクトル、そして逆変換することで元に戻した画像を可視化している。 なお、読み込む画像は適当に用意しよう。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from PIL import Image
import numpy as np
from matplotlib import pyplot as plt


def main():
    # 画像を読み込む
    img = Image.open('lena.png')
    # グレイスケールに変換する
    gray_img = img.convert('L')
    # NumPy 配列にする
    f_xy = np.asarray(gray_img)

    # 2 次元高速フーリエ変換で周波数領域の情報を取り出す
    f_uv = np.fft.fft2(f_xy)
    # 画像の中心に低周波数の成分がくるように並べかえる
    shifted_f_uv = np.fft.fftshift(f_uv)

    # パワースペクトルに変換する
    magnitude_spectrum2d = 20 * np.log(np.absolute(shifted_f_uv))

    # 元の並びに直す
    unshifted_f_uv = np.fft.fftshift(shifted_f_uv)
    # 2 次元逆高速フーリエ変換で空間領域の情報に戻す
    i_f_xy = np.fft.ifft2(unshifted_f_uv).real  # 実数部だけ使う

    # 上記を画像として可視化する
    fig, axes = plt.subplots(1, 3, figsize=(8, 4))
    # 枠線と目盛りを消す
    for ax in axes:
        for spine in ax.spines.values():
            spine.set_visible(False)
        ax.set_xticks([])
        ax.set_yticks([])
    # 元画像
    axes[0].imshow(f_xy, cmap='gray')
    axes[0].set_title('Input Image')
    # 周波数領域のパワースペクトル
    axes[1].imshow(magnitude_spectrum2d, cmap='gray')
    axes[1].set_title('Magnitude Spectrum')
    # FFT -> IFFT した画像
    axes[2].imshow(i_f_xy, cmap='gray')
    axes[2].set_title('Reversed Image')
    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記をファイルに保存して実行する。

$ python imgfft2d.py

すると、次のようなグラフが得られる。

FFT -> IFFT

真ん中のプロットが、元となった画像データの周波数領域での表現になる。 この表現では、中心に近いほど低い周波数・遠いほど高い周波数の成分を含む。 白い部分ほど成分が多いことを示しているため、この画像は低周波の成分が比較的多いように見える。 また、右の画像を見ると、ちゃんと周波数領域のデータから空間領域のデータに復元できていることがわかる。

ローパスフィルタ (Low Pass Filter) をかけて復元してみる

先ほど述べた通り、周波数領域の表現では中心に近いほど低い周波数・遠いほど高い周波数の成分を含んでいる。 つまり、中心部分だけを取り出すような演算をすると、画像を構成する波形にローパスフィルタをかけることができる。

以下のサンプルコードでは、周波数領域のデータに対して中心部分を取り出すフィルタをかけることで低周波成分だけを抽出した。 そして、フィルタしたデータを逆フーリエ変換することで画像に戻している。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from PIL import Image
from PIL import ImageDraw
import numpy as np
from matplotlib import pyplot as plt


def main():
    # 画像を読み込む
    img = Image.open('lena.png')
    # グレイスケールに変換する
    gray_img = img.convert('L')
    # NumPy 配列にする
    f_xy = np.asarray(gray_img)

    # 2 次元高速フーリエ変換で周波数領域の情報を取り出す
    f_uv = np.fft.fft2(f_xy)
    # 画像の中心に低周波数の成分がくるように並べかえる
    shifted_f_uv = np.fft.fftshift(f_uv)

    # フィルタ (ローパス) を用意する
    x_pass_filter = Image.new(mode='L',  # 8-bit pixels, black and white
                              size=(shifted_f_uv.shape[0],
                                    shifted_f_uv.shape[1]),
                              color=0,  # default black
                              )
    # 中心に円を描く
    draw = ImageDraw.Draw(x_pass_filter)
    # 円の半径
    ellipse_r = 50
    # 画像の中心
    center = (shifted_f_uv.shape[0] // 2,
              shifted_f_uv.shape[1] // 2)
    # 円の座標
    ellipse_pos = (center[0] - ellipse_r,
                   center[1] - ellipse_r,
                   center[0] + ellipse_r,
                   center[1] + ellipse_r)
    draw.ellipse(ellipse_pos, fill=255)
    # フィルタ
    filter_array = np.asarray(x_pass_filter)

    # フィルタを適用する
    filtered_f_uv = np.multiply(shifted_f_uv, filter_array)

    # パワースペクトルに変換する
    magnitude_spectrum2d = 20 * np.log(np.absolute(filtered_f_uv))

    # 元の並びに直す
    unshifted_f_uv = np.fft.fftshift(filtered_f_uv)
    # 2 次元逆高速フーリエ変換で空間領域の情報に戻す
    i_f_xy = np.fft.ifft2(unshifted_f_uv).real  # 実数部だけ使う

    # 上記を画像として可視化する
    fig, axes = plt.subplots(1, 4, figsize=(12, 4))
    # 枠線と目盛りを消す
    for ax in axes:
        for spine in ax.spines.values():
            spine.set_visible(False)
        ax.set_xticks([])
        ax.set_yticks([])
    # 元画像
    axes[0].imshow(f_xy, cmap='gray')
    axes[0].set_title('Input Image')
    # フィルタ画像
    axes[1].imshow(filter_array, cmap='gray')
    axes[1].set_title('Filter Image')
    # フィルタされた周波数領域のパワースペクトル
    axes[2].imshow(magnitude_spectrum2d, cmap='gray')
    axes[2].set_title('Filtered Magnitude Spectrum')
    # FFT -> Band-pass Filter -> IFFT した画像
    axes[3].imshow(i_f_xy, cmap='gray')
    axes[3].set_title('Reversed Image')
    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python imglpf.py

すると、次のようなグラフが得られる。 分かりにくいかもしれないけど、復元した画像はちょっとぼやけている。 また、ぼやけ方も縞模様になっていて波を感じるものとなった。 これは、画像を構成するすべての周波数の中から、高周波の成分がフィルタによって取り除かれたことで生じている。

FFT -> LPF -> IFFT

中心のごく一部だけを取り出している (つまり、多くの情報が失われている) のに、それっぽい画像になるのはなんとも面白い。 なお、フィルタに使う円の半径を小さくすれば、それだけ高周波の成分が少なくなって、より分かりやすくぼやける。

ハイパスフィルタ (High Pass Filter) をかけて復元してみる

次は、同様に高周波成分だけを取り出すハイパスフィルタをかけてみよう。 つまり、中心だけを取り除くことになる。 先ほどとのコードの違いは、適用するフィルタの周波数特性が違うだけ。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from PIL import Image
from PIL import ImageDraw
import numpy as np
from matplotlib import pyplot as plt


def main():
    # 画像を読み込む
    img = Image.open('lena.png')
    # グレイスケールに変換する
    gray_img = img.convert('L')
    # NumPy 配列にする
    f_xy = np.asarray(gray_img)

    # 2 次元高速フーリエ変換で周波数領域の情報を取り出す
    f_uv = np.fft.fft2(f_xy)
    # 画像の中心に低周波数の成分がくるように並べかえる
    shifted_f_uv = np.fft.fftshift(f_uv)

    # フィルタ (ハイパス) を用意する
    x_pass_filter = Image.new(mode='L',  # 8-bit pixels, black and white
                              size=(shifted_f_uv.shape[0],
                                    shifted_f_uv.shape[1]),
                              color=255,  # default white
                              )
    # 中心に円を描く
    draw = ImageDraw.Draw(x_pass_filter)
    # 円の半径
    ellipse_r = 50
    # 画像の中心
    center = (shifted_f_uv.shape[0] // 2,
              shifted_f_uv.shape[1] // 2)
    # 円の座標
    ellipse_pos = (center[0] - ellipse_r,
                   center[1] - ellipse_r,
                   center[0] + ellipse_r,
                   center[1] + ellipse_r)
    draw.ellipse(ellipse_pos, fill=0)
    # フィルタ
    filter_array = np.asarray(x_pass_filter)

    # フィルタを適用する
    filtered_f_uv = np.multiply(shifted_f_uv, filter_array)

    # パワースペクトルに変換する
    magnitude_spectrum2d = 20 * np.log(np.absolute(filtered_f_uv))

    # 元の並びに直す
    unshifted_f_uv = np.fft.fftshift(filtered_f_uv)
    # 2 次元逆高速フーリエ変換で空間領域の情報に戻す
    i_f_xy = np.fft.ifft2(unshifted_f_uv).real  # 実数部だけ使う

    # 上記を画像として可視化する
    fig, axes = plt.subplots(1, 4, figsize=(12, 4))
    # 枠線と目盛りを消す
    for ax in axes:
        for spine in ax.spines.values():
            spine.set_visible(False)
        ax.set_xticks([])
        ax.set_yticks([])
    # 元画像
    axes[0].imshow(f_xy, cmap='gray')
    axes[0].set_title('Input Image')
    # フィルタ画像
    axes[1].imshow(filter_array, cmap='gray')
    axes[1].set_title('Filter Image')
    # フィルタされた周波数領域のパワースペクトル
    axes[2].imshow(magnitude_spectrum2d, cmap='gray')
    axes[2].set_title('Filtered Magnitude Spectrum')
    # FFT -> Band-pass Filter -> IFFT した画像
    axes[3].imshow(i_f_xy, cmap='gray')
    axes[3].set_title('Reversed Image')
    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を保存して実行する。

$ python imghpf.py

すると、今度は次のようなグラフが得られる。 またもや分かりにくいけど、復元した画像は元の画像の輪郭だけがぼんやりと浮かび上がったものとなっている。 これが画像を構成する高周波の成分で、先ほどの復元画像から取り除かれたものと解釈できる。

FFT -> HPF -> IFFT

カラー画像を扱う場合

なお、カラー画像を扱う場合、それぞれのチャネルごとに処理する必要がある。 以下のサンプルコードでは、最初のサンプルコードを RGB のチャネルでそれぞれ適用している。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from itertools import chain

from PIL import Image
import numpy as np
from matplotlib import pyplot as plt


def main():
    # 画像を読み込む
    img = Image.open('lena.png')
    # NumPy 配列にする
    f_xy_rgb = np.asarray(img)

    # 画像として可視化する
    fig, axes = plt.subplots(3, 5, figsize=(12, 6))
    # 枠線と目盛りを消す
    for ax in chain.from_iterable(axes):
        for spine in ax.spines.values():
            spine.set_visible(False)
        ax.set_xticks([])
        ax.set_yticks([])

    # 元画像
    axes[1][0].imshow(f_xy_rgb)
    axes[1][0].set_title('Input Image (RGB)')

    # 逆変換した画像を入れる領域を用意しておく
    i_f_xy_rgb = np.empty_like(f_xy_rgb)
    # RGB ごとに処理する
    channel_names = ('R', 'G', 'B')
    channel_indices = range(f_xy_rgb.shape[2])
    zipped_channel_info = zip(channel_names, channel_indices)
    for channel_name, channel_index in zipped_channel_info:
        # 各チャネルごとに配列として取り出す
        f_xy = f_xy_rgb[:, :, channel_index]

        # 各チャネルの元画像
        axes[channel_index][1].imshow(f_xy, cmap='gray')
        axes[channel_index][1].set_title(f'Input Image ({channel_name})')

        # 2 次元高速フーリエ変換で周波数領域の情報を取り出す
        f_uv = np.fft.fft2(f_xy)
        # 画像の中心に低周波数の成分がくるように並べかえる
        shifted_f_uv = np.fft.fftshift(f_uv)
        # パワースペクトルに変換する
        magnitude_spectrum2d = 20 * np.log(np.absolute(shifted_f_uv))
        # 元の並びに直す
        unshifted_f_uv = np.fft.fftshift(shifted_f_uv)
        # 2 次元逆高速フーリエ変換で空間領域の情報に戻す
        i_f_xy = np.fft.ifft2(unshifted_f_uv).real  # 実数部だけ使う

        # 周波数領域のパワースペクトル
        axes[channel_index][2].imshow(magnitude_spectrum2d, cmap='gray')
        axes[channel_index][2].set_title(f'Magnitude Spectrum ({channel_name})')

        # FFT -> Band-pass Filter -> IFFT した画像
        axes[channel_index][3].imshow(i_f_xy, cmap='gray')
        axes[channel_index][3].set_title(f'Reversed Image ({channel_name})')

        # 逆変換したチャネルを保存しておく
        i_f_xy_rgb[:, :, channel_index] = i_f_xy

    # 逆変換した RGB 画像
    axes[1][4].imshow(i_f_xy_rgb)
    axes[1][4].set_title('Reversed Image (RGB)')

    # グラフを表示する
    plt.show()


if __name__ == '__main__':
    main()

上記を保存して実行する。

$ python imgfftc.py

すると、次のようなグラフが得られる。

FFT (RGB) -> IFFT (RGB)

ちゃんと RGB それぞれのチャネルの情報から、カラー画像が復元できていることがわかる。

めでたしめでたし。

ディジタル画像処理[改訂第二版]

ディジタル画像処理[改訂第二版]

Amazon

Python: UMAP を使ってみる

Python scikit-learn 機械学習 macOS

UMAP (Uniform Manifold Approximation and Projection) は次元削減手法のひとつ。 似た手法としては t-SNE (t-distributed Stochastic Neighbor Embedding) があるけど、それよりも高速らしい。 公式のベンチマークが以下で紹介されていて、t-SNE に比べるとスケーラビリティに優れていることが伺える。

umap-learn.readthedocs.io

使った環境は次のとおり。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G5033
$ python -V                          
Python 3.7.7
$ pip list | grep -i umap            
umap-learn                0.4.4

もくじ

下準備

$ pip install umap-learn scikit-learn

Digit データセットを次元削減してみる

scikit-learn に付属している数字の画像データセットを使って動作を確認してみよう。 UMAP は scikit-learn に準拠したインターフェースを提供している。 以下のサンプルコードでは、元が 8 x 8 ピクセルから成る 64 次元の画像データを、UMAP で 2 次元にまで削減している。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
from matplotlib import pyplot as plt
from sklearn import datasets
import umap


def main():
    # データセットを読み込む
    dataset = datasets.load_digits()
    X, y = dataset.data, dataset.target

    # 次元削減する
    mapper = umap.UMAP(random_state=0)
    embedding = mapper.fit_transform(X)

    # 結果を二次元でプロットする
    embedding_x = embedding[:, 0]
    embedding_y = embedding[:, 1]
    for n in np.unique(y):
        plt.scatter(embedding_x[y == n],
                    embedding_y[y == n],
                    label=n)

    # グラフを表示する
    plt.grid()
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記を保存して実行する。

$ python digitsumap.py

すると、次のようなグラフが得られる。

t-SNE をはじめて使ったときも思ったけど、教師なしでここまで同じラベルのデータがまとまるのは不思議。

t-SNE も試してみる

一応、比較対象として t-SNE も試しておく。 以下のサンプルコードは、基本的に先ほどのコードで umap.UMAP() を sklearn.manifold.TSNE に変えただけ。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn import manifold


def main():
    # データセットを読み込む
    dataset = datasets.load_digits()
    X, y = dataset.data, dataset.target

    # 次元削減する
    mapper = manifold.TSNE(random_state=0)
    embedding = mapper.fit_transform(X)

    # 結果を二次元でプロットする
    embedding_x = embedding[:, 0]
    embedding_y = embedding[:, 1]
    for n in np.unique(y):
        plt.scatter(embedding_x[y == n],
                    embedding_y[y == n],
                    label=n)

    # グラフを表示する
    plt.grid()
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

上記を保存して実行する。

$ digitstsne.py

すると、次のようなグラフが得られる。

おそらく学習に使うパラメータにもよるだろうけど、UMAP の方がぎゅっとまとまってる感じかな。

いじょう。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

Amazon

Python: XGBoost の cv() 関数から学習済みモデルを取り出す

Python XGBoost 機械学習 scikit-learn macOS

今回は、以下のエントリを XGBoost で焼き直したもの。 つまり、XGBoost でも cv() 関数から学習済みモデルを取り出して Fold Averaging してみようという話。

blog.amedama.jp

使った環境は次のとおり。

$ sw_vers                   
ProductName:    Mac OS X
ProductVersion: 10.14.6
BuildVersion:   18G5033
$ python -V              
Python 3.7.7
$ pip list | grep xgboost               
xgboost                    1.1.1

　下準備

必要なパッケージをインストールしておく。

$ pip install xgboost scikit-learn numpy

学習済みモデルを取り出して Fold Averaging してみる

早速、以下にサンプルコードを示す。 乳がんデータセットをホールドアウトして、一方のデータで学習して、他方のデータを Fold Averaging している。

実装方法としては、LightGBM と同じようにコールバックを使って学習済みモデルへの参照を引っ張り出した。 cv() 関数のコールバックには cvfolds というパラメータ名で xgboost.training.CVPack のリストが渡される。 あとは CVPack#bst という名前で Booster オブジェクトにアクセスするだけ。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import inspect

import numpy as np
import xgboost as xgb
from sklearn import datasets
from sklearn.metrics import accuracy_score
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import train_test_split


class ModelExtractionCallback:
    """XGBoost の cv() 関数から学習済みモデルを取り出すためのコールバック"""

    def __init__(self):
        self.cvfolds = None

    def __call__(self, env):
        # コールバックの呼び出しで xgboost.training.CVPack のリストが得られる
        if self.cvfolds is None:
            self.cvfolds = env.cvfolds


class BoostersProxy:
    """コールバックから得られた CVPack のリストを使って Fold Averaging をやりやすくするクラス

    実用上は ModelExtractionCallback とニコイチしちゃっても良いけど一応分けた"""

    def __init__(self, cvfolds):
        self._cvfolds = cvfolds

    def __getattr__(self, name):
        def _wrap(*args, **kwargs):
            ret = []
            for cvpack in self._cvfolds:
                # それぞれの Booster から指定されたアトリビュートを取り出す
                attr = getattr(cvpack.bst, name)
                if inspect.ismethod(attr):
                    # オブジェクトがメソッドなら呼び出した結果を入れる
                    res = attr(*args, **kwargs)
                    ret.append(res)
                else:
                    # それ以外ならそのまま入れる
                    ret.append(attr)
            return ret
        return _wrap


def main():
    # データセットを読み込む
    dataset = datasets.load_breast_cancer()
    X, y = dataset.data, dataset.target

    # ホールドアウト検証のためにデータセットを分割する
    X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                        test_size=0.33,
                                                        shuffle=True,
                                                        random_state=42,
                                                        stratify=y)

    # XGBoost のデータセット表現に直す
    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test, label=y_test)

    # データの分割に使うコンテキスト
    folds = StratifiedKFold(n_splits=5,
                            shuffle=True,
                            random_state=42)
    # 学習に使うパラメータ
    xgb_params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
    }

    # モデルを取り出すのに使うコールバック
    extraction_cb = ModelExtractionCallback()
    callbacks = [
        extraction_cb,
    ]

    # 交差検証する
    xgb.cv(xgb_params,
           dtrain,
           num_boost_round=1000,
           early_stopping_rounds=100,
           folds=folds,
           verbose_eval=True,
           # コールバックを渡す
           callbacks=callbacks,
           )

    # コールバックから学習済みモデルを取り出してプロキシにくべる
    proxy = BoostersProxy(cvfolds=extraction_cb.cvfolds)

    # ホールドアウトしておいた検証データを Fold Averaging で予測する
    y_pred_proba_list = proxy.predict(dtest)
    y_pred_proba_avg = np.array(y_pred_proba_list).mean(axis=0)
    y_pred = np.where(y_pred_proba_avg > 0.5, 1, 0)
    accuracy = accuracy_score(y_test, y_pred)
    print('Fold averaging accuracy:', accuracy)


if __name__ == '__main__':
    main()

上記を実行してみよう。

$ python xgbcv.py 
[0]    train-logloss:0.46506+0.00234  test-logloss:0.48322+0.01040
[1]    train-logloss:0.33602+0.00369  test-logloss:0.36709+0.01240
[2]    train-logloss:0.25204+0.00529  test-logloss:0.29438+0.01717

...

[130]  train-logloss:0.00727+0.00008  test-logloss:0.11461+0.07024
[131]  train-logloss:0.00727+0.00008  test-logloss:0.11468+0.07019
[132]  train-logloss:0.00727+0.00008  test-logloss:0.11473+0.07017
Fold averaging accuracy: 0.9627659574468085

ホールドアウトしたデータを Fold Averaging で予測して約 0.962 という Accuracy スコアが得られた。

いじょう。

Kaggleで勝つデータ分析の技術

Kaggleで勝つデータ分析の技術

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＰｙｔｈｏｎではじめるＫａｇｇｌｅスタートブック (ＫＳ情報科学専門書)

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＰｙｔｈｏｎではじめるＫａｇｇｌｅスタートブック (ＫＳ情報科学専門書)