Visual C++にて、CPUにAVX命令があるかどうかを実行中に知る方法

表題のことは、__cpuid, __cpuidexにあるサンプルコードで実現できます。Visual Studio 2017で確認しました。

自分の初代Core i7では、下記のようにAVXに対応していないことがわかります。

GenuineIntel
Intel(R) Core(TM) i7 CPU         870  @ 2.93GHz
3DNOW not supported
3DNOWEXT not supported
ABM not supported
ADX not supported
AES not supported
AVX not supported
AVX2 not supported
AVX512CD not supported
AVX512ER not supported
AVX512F not supported
AVX512PF not supported
BMI1 not supported
BMI2 not supported
CLFSH supported
CMPXCHG16B supported
CX8 supported
ERMS not supported
F16C not supported
FMA not supported
FSGSBASE not supported
FXSR supported
HLE not supported
INVPCID not supported
LAHF supported
LZCNT not supported
MMX supported
MMXEXT not supported
MONITOR supported
MOVBE not supported
MSR supported
OSXSAVE not supported
PCLMULQDQ not supported
POPCNT supported
PREFETCHWT1 not supported
RDRAND not supported
RDSEED not supported
RDTSCP supported
RTM not supported
SEP supported
SHA not supported
SSE supported
SSE2 supported
SSE3 supported
SSE4.1 supported
SSE4.2 supported
SSE4a not supported
SSSE3 supported
SYSCALL not supported
TBM not supported
XOP not supported
XSAVE not supported

g++によるhello world再訪

g++でhello.cppから実行ファイルを生成するときに何が起こっているかを調べてまとめました。全体的に gcc Compilation Process and Steps of C Program in Linux を参考にしています。元記事はCファイルを対象にしていますが、本記事ではC++ファイルを対象としています。

ソースコード

以下のコード hello.cpp をステップバイステップでビルドしていきます。このコードはC++11相当のコードです。

#include <iostream>
#include <vector>
#define MESSAGE "hello  "

int main(void)
{
    std::vector<int> vec{1,2,3,4,5};
    for(auto e: vec) {
        std::cout << MESSAGE << __FILE__ << " " << e << std::endl;
    }
    
    return 0;
}

1. プリプロセッシング

プリプロセッシングの結果を見る

まず最初に、プリプロセッシングと呼ばれる前処理が行われます。具体的には、

  • #includeでインクルードしたヘッダファイルの展開
  • #defineで定義した文字列の置換処理
  • __FILE____LINE__などの記号定数の置換処理

が行われます。(参考:マクロ展開の例

$ cpp -x c++ -std=c++11 hello.cpp > hello.i

または

$ g++ -E -std=c++11 hello.cpp > hello.i

で、プロプロセッシング後のファイルをhello.iという名前で保存できます。

hello.iの抜粋を以下に示します。

# 1 "hello.cpp"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/usr/include/stdc-predef.h" 1 3 4
# 1 "<command-line>" 2
# 1 "hello.cpp"
# 1 "/usr/include/c++/5/iostream" 1 3
# 36 "/usr/include/c++/5/iostream" 3
...(中略)...
# 5 "hello.cpp"
int main(void)
{
    std::vector<int> vec{1,2,3,4,5};
    for(auto e: vec) {
        std::cout << "hello  " << "hello.cpp" << " " << e << std::endl;
    }

    return 0;
}

ヘッダファイルが検索されるディレクトリの場所を調べる

ヘッダファイルが捜索されるディレクトリの場所は、Cの場合は

$ gcc -xc -E -v -

C++の場合は

$ gcc -xc++ -E -v -

で調べることができます(参考:c++ - What are the GCC default include directories? - Stack Overflow)。私の場合は以下のように表示されました。

()
#include "..." search starts here:
#include <...> search starts here:
 /usr/include/c++/5
 /usr/include/x86_64-linux-gnu/c++/5
 /usr/include/c++/5/backward
 /usr/lib/gcc/x86_64-linux-gnu/5/include
 /usr/local/include
 /usr/lib/gcc/x86_64-linux-gnu/5/include-fixed
 /usr/include/x86_64-linux-gnu
 /usr/include
End of search list.

includeされるヘッダファイルのフルパスを調べる

ヘッダファイルは芋づる式にincludeされていくわけですが、具体的にどのパスにあるどのヘッダファイルがincludeされたかは、

$ g++ -std=c++11 -H hello.cpp

で調べることができます(参考:c++ - How to tell where a header file is included from? - Stack Overflow)。私の環境だと以下のようにツリーが表示されました。

$ g++ -H hello.cpp       
. /usr/include/c++/5/iostream
.. /usr/include/x86_64-linux-gnu/c++/5/bits/c++config.h
... /usr/include/x86_64-linux-gnu/c++/5/bits/os_defines.h
.... /usr/include/features.h
..... /usr/include/x86_64-linux-gnu/sys/cdefs.h
...... /usr/include/x86_64-linux-gnu/bits/wordsize.h
..... /usr/include/x86_64-linux-gnu/gnu/stubs.h
...... /usr/include/x86_64-linux-gnu/gnu/stubs-64.h
...(略)...

2. アセンブリコードへの変換

次に、コンパイラは、hello.iアセンブリコードに変換します。

$ g++ -std=c++11 -S hello.i

とすると、アセンブリコードhello.sが生成されます。

        .file   "hello.cpp"
        .section        .rodata
        .type   _ZStL19piecewise_construct, @object
        .size   _ZStL19piecewise_construct, 1
_ZStL19piecewise_construct:
        .zero   1
        .local  _ZStL8__ioinit
        .comm   _ZStL8__ioinit,1,1
...(略)...

3. オブジェクトファイルへの変換

次に、asコマンドを使って、アセンブリコードhello.sからマシン語に変換し、オブジェクトファイルhello.oを得ます。

$ as hello.s -o hello.o

hello.oの中身は

$ objdump -d hello.o

で調べることができます。以下は実行例です。

0000000000000000 <main>:
   0:   55                      push   %rbp
   1:   48 89 e5                mov    %rsp,%rbp
   4:   41 55                   push   %r13
   6:   41 54                   push   %r12
   8:   53                      push   %rbx
   9:   48 83 ec 58             sub    $0x58,%rsp
   d:   64 48 8b 04 25 28 00    mov    %fs:0x28,%rax
...(略)...

4. リンク

実行ファイルを生成する

最後に、さっき生成したオブジェクトファイルhello.oおよび、他に必要なオブジェクトファイルを集めてリンクし、実行ファイルを生成します。

$ g++ hello.o -o hello

とすると、実行ファイルhelloが生成されます。

ライブラリが検索されるディレクトリの場所、リンクされるライブラリのフルパスを調べる

$ g++ -v hello.o -o hello

とすると、以下のように詳細情報が表示されます。

LIBRARY_PATH=/usr/lib/gcc/x86_64-linux-gnu/5/:/usr/lib/gcc/x86_64-linux-gnu/5/../../../x86_64-linux-gnu/:/usr/lib/gcc/x86_64-linux-gnu/5/../../../../lib/:/lib/x86_64-linux-gnu/:/lib/../lib/:/usr/lib/x86_64-linux-gnu/:/usr/lib/../lib/:/usr/lib/gcc/x86_64-linux-gnu/5/../../../:/lib/:/usr/lib/

は、ライブラリが検索されるディレクトリの場所を表していると思います。

また、

COLLECT_GCC_OPTIONS='-v' '-o' 'hello' '-shared-libgcc' '-mtune=generic' '-march=x86-64'
 /usr/lib/gcc/x86_64-linux-gnu/5/collect2 -plugin /usr/lib/gcc/x86_64-linux-gnu/5/liblto_plugin.so -plugin-opt=/usr/lib/gcc/x86_64-linux-gnu/5/lto-wrapper -plugin-opt=-fresolution=/tmp/ccvdOxv1.res -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lc -plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lgcc --sysroot=/ --build-id --eh-frame-hdr -m elf_x86_64 --hash-style=gnu --as-needed -dynamic-linker /lib64/ld-linux-x86-64.so.2 -z relro -o hello /usr/lib/gcc/x86_64-linux-gnu/5/../../../x86_64-linux-gnu/crt1.o /usr/lib/gcc/x86_64-linux-gnu/5/../../../x86_64-linux-gnu/crti.o /usr/lib/gcc/x86_64-linux-gnu/5/crtbegin.o -L/usr/lib/gcc/x86_64-linux-gnu/5 -L/usr/lib/gcc/x86_64-linux-gnu/5/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/5/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/5/../../.. hello.o -lstdc++ -lm -lgcc_s -lgcc -lc -lgcc_s -lgcc /usr/lib/gcc/x86_64-linux-gnu/5/crtend.o /usr/lib/gcc/x86_64-linux-gnu/5/../../../x86_64-linux-gnu/crtn.o

は、リンク時に使われたコマンドの詳細を表していると思います。collect2とはgcc/g++が内部で呼び出すリンクのためのコマンドのようです。ldとの違いはよくわかりません…。

その他の情報

$ g++ -save-temps -std=c++11 hello.cpp

とすると、上記でステップバイステップで生成した中間ファイル hello.ii, hello.s, hello.o が一気に生成されます。

Ubuntuの「壊れた変更禁止パッケージがあります」エラーにaptitudeで対処する方法

Ubuntu 16.04にて、少し前からapt-getで新しいパッケージを入れようとすると「壊れた変更禁止パッケージがあります」エラーが出てインストールできない事象に悩まされていました。

$ sudo apt-get install libgfortran3  
パッケージリストを読み込んでいます... 完了
依存関係ツリーを作成しています                
状態情報を読み取っています... 完了
インストールすることができないパッケージがありました。おそらく、あり得
ない状況を要求したか、(不安定版ディストリビューションを使用しているの
であれば) 必要なパッケージがまだ作成されていなかったり Incoming から移
動されていないことが考えられます。
以下の情報がこの問題を解決するために役立つかもしれません:

以下のパッケージには満たせない依存関係があります:
 libgfortran3 : 依存: gcc-5-base (= 5.3.1-14ubuntu2) しかし、5.4.0-6ubuntu1~16.04.4 はインストールされようとしています
E: 問題を解決することができません。壊れた変更禁止パッケージがあります。

aptitudeを使ってこの問題への対応を試みました(参考:apt - E: Unable to correct problems, you have held broken packages - Ask Ubuntu)。以下にその方法について記します。自分でもこれが正しい方法なのかよくわかっていません。試す場合は自己責任でお願いします。

$ sudo aptitude install libgfortran3

とすると、以下のように、壊れたパッケージを修正するための改善案が表示されます。改善案を受け入れるならY、他の改善案を探すならnです。

以下の新規パッケージがインストールされます:
  libgfortran3{b} 
0 個のパッケージを更新、 1 個を新たにインストール、 0 個を削除予定、0 個が更新されていない。
260 k バイトのアーカイブを取得する必要があります。 展開後に 1,290 k バイトのディスク領域が新たに消費されます。
以下のパッケージには満たされていない依存関係があります:
 libgfortran3 : 依存: gcc-5-base (= 5.3.1-14ubuntu2) [5.4.0-6ubuntu1~16.04.4 が既にインストール済みです]
以下のアクションでこれらの依存関係の問題は解決されます:

     以下のパッケージを現在のバージョンに一時固定する:
1)     libgfortran3 [インストールされていません]      



この解決方法を受け入れますか? [Y/n/q/?]

私の場合、Yとしてもパッケージをインストールできなかったので、nで次の改善案を探しました。その結果が以下です。

以下のアクションでこれらの依存関係の問題は解決されます:

      以下のパッケージを削除する:                                              
1)      build-essential                                                        
2)      g++                                                                    
3)      g++-5                                                                  
4)      gcc                                                                    
5)      gcc-5                                                                  
6)      libasan2                                                               
7)      libatomic1                                                             
8)      libcilkrts5                                                            
9)      libgcc-5-dev                                                           
10)     libitm1                                                                
11)     liblsan0                                                               
12)     libmpx0                                                                
13)     libstdc++-5-dev                                                        
14)     libtsan0                                                               
15)     libubsan0                                                              

      以下のパッケージをインストールする:                                      
16)     tcc [0.9.27~git20151227.933c223-1 (xenial)]                            

      以下のパッケージをダウングレードする:                                    
17)     cpp-5 [5.4.0-6ubuntu1~16.04.4 (now) -> 5.3.1-14ubuntu2 (xenial)]       
18)     gcc-5-base [5.4.0-6ubuntu1~16.04.4 (now) -> 5.3.1-14ubuntu2 (xenial)]  
19)     libcc1-0 [5.4.0-6ubuntu1~16.04.4 (now) -> 5.3.1-14ubuntu2 (xenial)]    
20)     libgomp1 [5.4.0-6ubuntu1~16.04.4 (now) -> 5.3.1-14ubuntu2 (xenial)]    
21)     libquadmath0 [5.4.0-6ubuntu1~16.04.4 (now) -> 5.3.1-14ubuntu2 (xenial)]
22)     libstdc++6 [5.4.0-6ubuntu1~16.04.4 (now) -> 5.3.1-14ubuntu2 (xenial)]  

      以下の依存関係を未解決のままにする:                                      
23)     codeblocks が gcc | g++ を推奨                                         
24)     dpkg-dev が build-essential を推奨                                     
25)     cmake が gcc を推奨                                                    


この解決方法を受け入れますか? [Y/n/q/?] 

ここでYとすると、無事所望のパッケージをインストールできました。しかし、今度はgccとg++がアンインストールされ、利用できなくなってしまいました。

sudo aptitude install g++

以下の新規パッケージがインストールされます:
  g++ g++-5{a} gcc{a} gcc-5{a} libasan2{a} libatomic1{a} libc-dev-bin{a} 
  libc6-dev{ab} libcilkrts5{a} libgcc-5-dev{a} libitm1{a} liblsan0{a} 
  libmpx0{a} libstdc++-5-dev{a} libtsan0{a} libubsan0{a} linux-libc-dev{a} 
  manpages-dev{a} 
0 個のパッケージを更新、 18 個を新たにインストール、 0 個を削除予定、0 個が更新されていない。
26.5 M バイトのアーカイブを取得する必要があります。 展開後に 101 M バイトのディスク領域が新たに消費されます。
以下のパッケージには満たされていない依存関係があります:
 libc6-dev : 依存: libc6 (= 2.23-0ubuntu3) [2.23-0ubuntu7 が既にインストール済みです]
以下のアクションでこれらの依存関係の問題は解決されます:

     以下のパッケージを現在のバージョンに一時固定する:    
1)     g++ [インストールされていません]                   
2)     g++-5 [インストールされていません]                 
3)     libc6-dev [インストールされていません]             
4)     libstdc++-5-dev [インストールされていません]       

     以下の依存関係を未解決のままにする:                  
5)     gcc が libc6-dev | libc-dev を推奨                 
6)     gcc-5 が libc6-dev (>= 2.13-0ubuntu6) を推奨       
7)     libgcc-5-dev が libc6-dev (>= 2.13-0ubuntu6) を推奨

Yで答えましたが、まだgccもg++も利用できません。再びsudo aptitude install g++で、いちどnしたあとの

以下のアクションでこれらの依存関係の問題は解決されます:

     以下のパッケージをダウングレードする:                      
1)     libc6 [2.23-0ubuntu7 (now) -> 2.23-0ubuntu3 (xenial)]    
2)     libc6-dbg [2.23-0ubuntu7 (now) -> 2.23-0ubuntu3 (xenial)]

Yで受け入れると、gccとg++が使えるようになりました。

apt-get updateで 「ターゲット ○○ は 複数回設定されています」というエラーに対応する

Ubuntu 16.04でsudo apt-get updateすると、以下のようなエラーが出ました。

W: ターゲット Packages (apps/binary-amd64/Packages) は /etc/apt/sources.list.d/getdeb.list:1 と /etc/apt/sources.list.d/getdeb.list:2 で複数回設定されています
W: ターゲット Packages (apps/binary-i386/Packages) は /etc/apt/sources.list.d/getdeb.list:1 と /etc/apt/sources.list.d/getdeb.list:2 で複数回設定されています
…
W: ターゲット DEP-11 (main/dep11/Components-amd64.yml) は /etc/apt/sources.list.d/google-chrome.list:3 と /etc/apt/sources.list.d/google.list:1 で複数回設定されています
W: ターゲット DEP-11-icons (main/dep11/icons-64x64.tar) は /etc/apt/sources.list.d/google-chrome.list:3 と /etc/apt/sources.list.d/google.list:1 で複数回設定されています

apt - How can I automatically fix W: Target Packages ... is configured multiple times? - Ask Ubuntu にあるとおりPythonスクリプトをダウンロードして実行すると、

$ sudo ./apt-remove-duplicate-source-entries.py 
Overlapping source entries:
  1. /etc/apt/sources.list.d/getdeb.list: deb http://archive.getdeb.net/ubuntu yakkety-getdeb apps
  2. /etc/apt/sources.list.d/getdeb.list: deb http://archive.getdeb.net/ubuntu yakkety-getdeb apps
I disabled the latter entry.

Overlapping source entries:
  1. /etc/apt/sources.list.d/google.list: deb [arch=amd64] http://dl.google.com/linux/chrome/deb/ stable main
  2. /etc/apt/sources.list.d/google-chrome.list: deb [arch=amd64] http://dl.google.com/linux/chrome/deb/ stable main
I disabled the latter entry.


2 source entries were disabled:
  # deb http://archive.getdeb.net/ubuntu yakkety-getdeb apps
  # deb [arch=amd64] http://dl.google.com/linux/chrome/deb/ stable main

Do you want to save these changes? (y/N)

と出たので、yとするとエラーが出なくなりました。実のところ、何が起こっているのかよくわかっていません。試す場合は自己責任でお願いします。

Fluent PythonはPython中級者にお勧めの名著

少し前のことですが、一年以上かかってのろのろ読んでいたFluent Python - O'Reilly Mediaをようやく読み終わりました。792ページもの厚さにもかかわらず密度が非常に濃く、どのページにも発見がある名著だと思います。オライリー公式での高評価も頷けます(今はなぜかレビューページがなくなっているようですが)。

本の特徴

本が執筆された当時のPython 3.4を対象に執筆されていますが、asyncioなど一部Python 3.5相当の機能についての言及もあります。また、Python 2.7との比較も要所でなされています。

この本の特徴は、Pythonicな書き方とは何か、なぜPythonの言語仕様がこう決められたのか、という疑問を、数多くの資料を援用して解き明かしていることだと思います。資料には、書籍のみならず、PEPsや、Stack Overflow、PyConのチュートリアル資料などのWeb上のリソースが数多く含まれます。いかにPythonが多くの議論のもと進化してきたかがよく分かるとともに、著者の博識さにも驚かされます。リンクが多いので電子書籍のほうが紙より楽かもしれません。

また、著者であるLuciano RamalhoはPythonの講師の経験も豊富とのことで、説明もわかりやすいと感じました。非ネイティブには辛い難解な言い回しやユーモアなどもなく、素直で読みやすい英語だと思います。

誰向けの本か

初心者向けの本ではなく、仕事や趣味ですでにある程度Pythonを使っている人向けです。Pythonicなコードとは何かを知りたい人におすすめできます。

勉強になった箇所

どの章も濃密ですが、以下は、特にこの本でしか読めない情報があると私が感じた箇所です。

ユニコード関係

第4章 “Text versus Bytes"は、ユニコードを中心としたテキスト処理に関する話題がよくまとまっていて、Pythonを使う人でなくても読んでは損はない内容だと思います。

自然言語処理系の人なら、テキストの前処理に使える情報を見つけられるかもしれません。例えばunicode.normalize('NFKC', ½)と呼ぶことで½1/2に置換できることなどが紹介されています。

並行処理関係

第17章 “Concurrency with Futures”, 第18章 “Concurrency with asyncio"と、かなりの部分が並行処理の説明に割かれています。特にasyncioはあまり解説がWebにないような気がするので、貴重な情報源だと思います。ただ、告白すると、このあたりはあまり消化できていない箇所です。

100 numpy exercisesの解説 76~100

100 numpy exercisesの解説 51~75 - minus9d’s diary の続きです。引き続き、 numpy-100/100 Numpy exercises.ipynb at master · rougier/numpy-100 · GitHub を片手にご覧ください。

76. Consider a one-dimensional array Z, build a two-dimensional array whose first row is (Z[0],Z[1],Z[2]) and each subsequent row is shifted by 1 (last row should be (Z[-3],Z[-2],Z[-1])

例えばZ = np.array([0, 1, …, 9])のとき、以下のような行列を作るのが目標です。

0 1 2
1 2 3
2 3 4
...
7 8 9

この回答例を理解するためには、strideという概念を理解する必要があります。

stridesとは、「Pythonによるデータ分析入門」を読む NumPy編 - 車輪を再発明 によると、"「『隣』の要素にアクセスするために、ポインタをいくら動かせばいいか」を表す数値"です。例えば、ある3x4行列のstrideをnp.arange(12).reshape(3,4).stridesというコードで調べると、(32,8)とわかります。8の意味は、「ポインタを一個右に動かすためには8バイト移動する」という意味で、32の意味は、「ポインタを一個下に動かすためには32バイト移動する」という意味です。np.arange(12).reshape(3,4)で生成した行列のdtypeはint64なので、計算が合っています。

回答例で使われているstride_tricks.as_stridedは、stridesの値を任意に指定することで、ポインタの進め方を偽装する関数です。回答例を以下に示します。

def rolling(a, window):
    shape = (a.size - window + 1, window)  # この場合 (8, 3)
    strides = (a.itemsize, a.itemsize)     # この場合 (8, 8)
    return stride_tricks.as_strided(a, shape=shape, strides=strides)  # (a)
Z = rolling(np.arange(10), 3)

式(a)は、「shapeが(10,)であるndarray aを、stridesを(8, 8)と偽装して、shapeが(8, 3)になるように舐めていく」という意味になります。stridesが(8, 8)ということは、「ポインタを一個右に動かすにも、ポインタを一個下に動かすにも、8バイト移動する」ということになります。

79. Consider 2 sets of points P0,P1 describing lines (2d) and a point p, how to compute distance from p to each line i (P0[i],P1[i])?

P0とP1とを結ぶ直線と、点pとの距離を求める問題です。想定解では10組分のP0, P1を考えているのでややこしいコードになっていますが、解法自体は高校数学の範囲です。

はてなブログで数式を綴るのが苦痛なので、ここからは画像で。

80. Consider an arbitrary array, write a function that extract a subpart with a fixed shape and centered on a given element (pad with a fill value when necessary)

理解できていません…。

81. Consider an array Z = [1,2,3,4,5,6,7,8,9,10,11,12,13,14], how to generate an array R = [[1,2,3,4], [2,3,4,5], [3,4,5,6], …, [11,12,13,14]]?

なぜかQ76と実質問題です。

82. Compute a matrix rank

matrix_rank()という関数が提供されている (numpy.linalg.matrix_rank — NumPy v1.13.dev0 Manual ) のでこれを使うのが良さそうです。

84. Extract all the contiguous 3x3 blocks from a random 10x10 matrix

Q76, Q81の応用編です。

86. Consider a set of p matrices wich shape (n,n) and a set of p vectors with shape (n,1). How to compute the sum of of the p matrix products at once? (result has shape (n,1))

テンソル積を行う関数tensordot()を用いて、pnnのテンソルと、pn1のテンソルの積をとり、n*1の結果をえる問題です。tensordot()についてはpython - Understanding tensordot - Stack Overflowを見るとよくイメージがわくと思います。

88. How to implement the Game of Life using numpy arrays?

ライフゲームの実装。各セルについて隣り合う8セルの状況を調べるのに、8通りにずらしたarrayを足し合わせる技に感心しました。

89. How to get the n largest values of an array

np.argsort(arr)は、配列arrを昇順にソートしたときの、配列のインデックスを返します。

Z = [9, 12, 16, 10, 14, 1, 4, 18, 13, 5]
print(np.argsort(Z))  # "[5 6 9 0 3 1 8 4 2 7]" とprintされる

np.argpartition(arr, idx)は、配列arrのうち、値が小さいものidx + 1個分のインデックスを返します。

Z = [9, 12, 16, 10, 14, 1, 4, 18, 13, 5]
print(np.argpartition(Z, 3))  # "[5 9 6 0 1 3 2 7 8 4] とprintされる

この場合、[5 9 6 0]までが、もっとも小さい要素ベスト4に対応するインデックスです。必ずしも小さい順になっている保証はないことに注意です。また、以降の[1 3 2 7 8 4]は何の意味もない値のはずです。

90. Given an arbitrary number of vectors, build the cartesian product (every combinations of every item)

デカルト積(=直積)を実装する問題です。pythonの標準ライブラリだと

list(itertools.product([1, 2, 3], [4, 5], [6, 7]))

というコードで直積を実現できますが、これをnumpyで実装しようという問題です。

回答例は理解できていません…。

92. Consider a large vector Z, compute Z to the power of 3 using 3 different methods

自分の手元では以下の結果でした。Q69に続き、np.einsum()は高速です。

1 loop, best of 3: 2.68 s per loop
1 loop, best of 3: 452 ms per loop
1 loop, best of 3: 299 ms per loop

93. Consider two arrays A and B of shape (8,3) and (2,2). How to find rows of A that contain elements of each row of B regardless of the order of the elements in B?

Aを構成する各行のうち、「Bを構成するどの行とも、共通する要素が存在する」ような行を探すという問題です。 一見してわかりにくい問題なので、例を見ます。

np.random.seed(0)
A = np.random.randint(0,5,(8,3))
B = np.random.randint(0,5,(2,2))

で生成されるのは

A
[[4 0 3]
 [3 3 1]
 [3 2 4]
 [0 0 4]
 [2 1 0]
 [1 1 0]
 [1 4 3]
 [0 3 0]]
B
[[2 3]
 [0 1]]

です。Aの最初の行[4 0 3]は、3がBの0行目に、1がBの1行目に存在するので、題意を満たします。一方、Aの3行目[3 2 4]は、Bの1行目の要素を一つも持たないため、題意を満たしません。

回答例をうまく説明することが難しいので、中間値を示すことでお茶を濁します。

C = A[..., np.newaxis, np.newaxis] == B

により、8x3x3x3のTrue, False表を作成します。さらに

C.any((3,1))

として

[[ True  True]
 [ True  True]
 [ True False]
 [False  True]
 [ True  True]
 [False  True]
 [ True  True]
 [ True  True]]
C.any((3,1)).all()

として

[ True  True False False  True False  True  True]

となります。

96. Given a two dimensional array, how to extract unique rows?

理解できていません…。

98. Considering a path described by two vectors (X,Y), how to sample it using equidistant samples (★★★)?

ある軌跡の座標列がベクトルX, Yで与えられるとき、この軌跡をリサンプルして補間する?という問題のようです。要再読

99. Given an integer n and a 2D array X, select from X the rows which can be interpreted as draws from a multinomial distribution with n degrees, i.e., the rows which only contain integers and which sum to n.

行列Xの各行のうち、整数のみを含み、合計がnになる行を探す問題です。M = np.logical_and.reduce(np.mod(X, 1) == 0, axis=-1)が全要素が整数であるという条件を、(X.sum(axis=-1) == n)が合計nであるという条件を表しています。

100. Compute bootstrapped 95% confidence intervals for the mean of a 1D array X (i.e., resample the elements of an array with replacement N times, compute the mean of each sample, and then compute percentiles over the means).

bootstrapとは、与えられた一つのサンプルから、シミュレーションにより統計値を求める方法です。

  • X = np.random.randn(100)
    • 長さ100のサンプルが一つで与えられます。
  • N = 1000; idx = np.random.randint(0, X.size, (N, X.size))
    • 1000回、復元抽出ありで、100個の要素を取り出すシミュレーションを実施します。
  • means = X[idx].mean(axis=1)
    • 1000回のシミュレーションのそれぞれで取り出した要素の平均値をとります
  • confint = np.percentile(means, [2.5, 97.5])
    • 1000個分の平均値をソートして、中央の95%の範囲をとったものが、95%の信頼区間になります。

100 numpy exercisesの解説 51~75

100 numpy exercisesの解説 1~50 - minus9d’s diary の続きです。コードの難易度が高く、一問一問を理解するための調査量が多すぎて、この記事で完了させることができませんでした。 numpy-100/100 Numpy exercises.ipynb at master · rougier/numpy-100 · GitHub を片手にご覧ください。

76問目~100問目は 100 numpy exercisesの解説 76~100 - minus9d’s diary に書きました。

52. Consider a random vector with shape (100,2) representing coordinates, find point by point distances

100点間の距離をすべて求めるという問題です。ここでは簡単化のために3点で考えます。

>>> Z = np.array([[1,1],[3,1],[2,5]])

np.at_least_2d()は、入力されたndarrayの次元数が最低でも2になるように次元を拡張する関数のようです。

3点のXとYを分離するのに、普通だと以下のようになるところが、

>>> X, Y = Z[:,0], Z[:,1]
>>> X.shape, Y.shape
((3,), (3,))

np.at_least_2d()を使うと、以下のようになります。

>>> X, Y = np.atleast_2d(Z[:,0], Z[:,1])
>>> X.shape, Y.shape
((1, 3), (1, 3))

次に、(X-X.T)**2の部分です。Xはarray([[1, 3, 2]])で1x3の形、X.Tはその転置なので

array([[1],
       [3],
       [2]])

で3x1の形です。

X-X.Tでは1x3の行列から3x1の行列を引くことになります。shapeが一致していないのに引き算なんてできるのかと思いますが、問題なくできます。これはNumPyのbroadcastingと呼ばれる性質に基づくものです。1.3.2. Numerical operations on arrays — Scipy lecture notes の “1.3.2.3. Broadcasting” 冒頭にあるイラストを見ると演算の中身が理解しやすいと思います。

結果、X-X.Tでは以下の結果が得られます。

array([[ 0,  2,  1],
       [-2,  0, -1],
       [-1,  1,  0]])

ここまでくればD = np.sqrt((X-X.T)**2 + (Y-Y.T)**2)の演算はできたも同然で、結果は

array([[ 0.        ,  2.        ,  4.12310563],
       [ 2.        ,  0.        ,  4.12310563],
       [ 4.12310563,  4.12310563,  0.        ]])

となります。

この問題、自分だったら以下のように書きたくなりそうです。

>>> Z = np.array([[1,1],[3,1],[2,5]])
>>> X, Y = Z[:,0], Z[:,1]
>>> X.shape, Y.shape
((3,), (3,))
>>> X_dist = np.subtract.outer(X, X)
>>> Y_dist = np.subtract.outer(Y, Y)
>>> X_dist.shape, Y_dist.shape
((3, 3), (3, 3))
>>> np.hypot(X_dist, Y_dist)
array([[ 0.        ,  2.        ,  4.12310563],
       [ 2.        ,  0.        ,  4.12310563],
       [ 4.12310563,  4.12310563,  0.        ]])

61. Find the nearest value from a given value in an array

以下の回答ではarrayとして長さ10の1次元arrayを例にとっています。flat()は任意次元のarrayを1次元とみなしてアクセスできるようにする関数ですが、この例の場合もとから1次元なので、Z.flatを使う必要性がなくなっています。

Z = np.random.uniform(0,1,10)
z = 0.5
m = Z.flat[np.abs(Z - z).argmin()]
print(m)

そこで、Zを

Z = np.random.uniform(0,1,(10,10))

とすると、Z.flatの有り難みがわかるようになります。

62. Considering two arrays with shape (1,3) and (3,1), how to compute their sum using an iterator?

np.nditer()は多次元の配列に対して要素をなめるときに使う関数で、numpy.nditer — NumPy v1.12 Manual を見てわかる通り、非常に多機能です。

まず基本形として、2x2の行列2つをイテレートする例を見てみます。

A = np.array([[1,2],[3,4]])
B = np.array([[5,6],[7,8]])
it = np.nditer([A,B])
for x, y in it:
    print("x y:", x, y)

を実行すると、

x y: 1 5
x y: 2 6
x y: 3 7
x y: 4 8

となります。

次に、上のコードで、Aが3x1の行列, Bが1x3の行列と変えることを考えます。このときも、Q52と同様、Broadcastingのルールが働きます(参考:NumPy Iterating Over Arrayの"Broadcasting Iteration")。

A = np.array([[1],[2],[3]])
B = np.array([[5,6,7]])
print(A.shape)
print(B.shape)
it = np.nditer([A,B])
for x, y in it:
    print("x y:", x, y)

実行結果は以下です。

x y: 1 5
x y: 1 6
x y: 1 7
x y: 2 5
x y: 2 6
x y: 2 7
x y: 3 5
x y: 3 6
x y: 3 7

ここまで理解した上で、元のコードを見てみます。

A = np.arange(3).reshape(3,1)
B = np.arange(3).reshape(1,3)
it = np.nditer([A,B,None])
for x,y,z in it: z[...] = x + y
print(it.operands[2])

A, Bに加えてNoneをも同時にイテレートして、値の代入までしてしまっているのがなんとも不思議なコードです。

63. Create an array class that has a name attribute

ndarrayをサブクラス化して自作クラスを作成する例です。Subclassing ndarray — NumPy v1.12 Manualに公式の説明がありますが、読んでいません。既存クラスのサブクラス化はハマりどころが多そうで、自分でやる気にはなれません。

64. Consider a given vector, how to add 1 to each element indexed by a second vector (be careful with repeated indices)?

np.bincount()は、値の生起回数をカウントする関数です。 まず基本形。np.bincount([3,5,3,8])を実行すると、array([0, 0, 0, 2, 0, 1, 0, 0, 1])が返ります。3が2回登場するので、3番目のindexが2になっていることがわかります。また、返る配列の長さは、登場する値の最大値で決まることもわかります。

これをふまえて問題文のコードを見てみます。I = np.random.randint(0,len(Z),20)array([8, 3, 8, 2, 0, 2, 4, 4, 8, 3, 8, 6, 7, 3, 5, 4, 0, 7, 8, 6])が生成されたものとします。0から9までの値をランダムに20個並べたものですが、今回はたまたま9が出なかったことに注意してください。

引数無しで単にnp.bincount(I)とすると、さきほど見たとおり、各数値をカウントした結果であるarray([2, 0, 2, 3, 3, 1, 2, 2, 5])が返ります。

ここでもしZ += np.bincount(I)とすると、Zの長さは10なのに対してnp.bincount(I)の長さは9のため、ValueErrorが起こります。

解答例のように、np.bincount(I, minlength=len(Z))minlength=10と指定することで、np.bincount()が返す配列の長さが最低でも10確保されるようになり、ValueErrorを防げます。

次に、別解であるnp.add.at(Z, I, 1)を読み解きます。numpy.ufunc.at — NumPy v1.12 Manualによると、この関数呼び出しはZ[I] += 1と等価です。ここで、Iはインデックス列です。

65. How to accumulate elements of a vector (X) to an array (F) based on an index list (I)?

  1. に引き続きbincount()に関する出題です。
X = [1,2,3,4,5,6]
I = [1,3,9,3,4,1]
F = np.bincount(I,X)
print(F)

np.bincount(I,X)は、np.bincount(I, weights=X)と等価です。つまり、値の生起回数に重みを与えます。実行結果は[ 0. 7. 0. 6. 5. 0. 0. 0. 0. 3.]です。

66. Considering a (w,h,3) image of (dtype=ubyte), compute the number of unique colors

F = I[...,0]*(256*256) + I[...,1]*256 +I[...,2]により(R, G, B)の値をR * 256 * 256 + G * 256 + Bという単一の値に変換したあと、np.unique()で同一の色を持つ要素を削除して、残った要素数を数えています。

C++unique()は隣り合う同一の数値しか削除しませんが、np.unique()の場合は隣り合わない同一の数値も削除するようです。

69. How to get the diagonal of a dot product?

第1の方法 np.diag(np.dot(A, B)) は愚直にAとBの積をとり、その対角成分を取得しています。対角成分以外の要素も計算しているので無駄が多いと言えます。

第2の方法 np.sum(A * B.T, axis=1) では、要素ごとの積をとる*を使用しています。紙と鉛筆で確かめるとわかりますが、AとBの積をとったときの対角成分に相当する部分のみを演算しているので高速化が期待できます。

第3の方法 np.einsum('ij,ji->i', A, B) は初めて見ました。アインシュタインの縮約記号にのっとって演算を行うための関数だそうです。詳細は【einsum】アインシュタインの縮約記法のように使えるnumpyの関数。性能と使い方を解説。 - プロクラシスト などをご覧ください。こちらも

実際に時間を計測した結果は以下の通りです。

>>> import timeit
>>> test_string = "import numpy as np; size = 5; A = np.random.uniform(0,1,(size,size)); B = np.random.uniform(0,1,(size,size));"
>>> timeit.timeit("np.diag(np.dot(A, B))", setup=test_string)
2.1652195840001696
>>> timeit.timeit("np.sum(A * B.T, axis=1)", setup=test_string)
3.46571562500003
>>> timeit.timeit("np.einsum('ij,ji->i', A, B)", setup=test_string)
0.8595781789999819

第3の方法がもっとも高速というのはよいとして、第1の方法より第2の方法が遅いのは予想外でした。そこで、行列のサイズを100x100に変えて再測定してみます。

>>> import timeit
>>> test_string = "import numpy as np; size = 100; A = np.random.uniform(0,1,(size,size)); B = np.random.uniform(0,1,(size,size));"
>>> timeit.timeit("np.diag(np.dot(A, B))", setup=test_string)
18.509561333999955
>>> timeit.timeit("np.sum(A * B.T, axis=1)", setup=test_string)
15.581413918999715
>>> timeit.timeit("np.einsum('ij,ji->i', A, B)", setup=test_string)
9.336061797000184

第2の方法が第1の方法より早くなりましたが、計算量から予想されるほどの差は出ませんでした。私の環境だと、第1の方法の実行中は4つの物理コアがすべてぶん回るのに対し、第2の方法と第3の方法では1つの物理コアしか回らないことがtopコマンドで確かめられたので、第1の方法には何らかの最適化が効いているのかもしれません。あくまで推測です。

71. Consider an array of dimension (5,5,3), how to mulitply it by an array with dimensions (5,5)?

(5,5)のshapeを持つBに対してB[:,:,None]とすると、(5,5,1)のshapeになるようです。

72. How to swap two rows of an array?

>>> A = np.arange(25).reshape(5,5)
>>> A
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])

上記のような(5,5)のshapeを持つarrayAに対して、A0,1というアクセスを行うと

>>> A[[0,1]]
array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]])

が返ります。このアクセス方法は"Advanced Indexing"と呼ばれる方法です(参考:Indexing — NumPy v1.12 Manual)。ふつうのindexingであるA[(0,1)]や、それの簡略系であるA[0,1]との違いに注意してください。対比のために書いておくと

>>> A[0,1]
1

です。Advanced Indexingについてまだ勉強不足ですが、公式リファレンスよりはとっつきやすい解説が NumPy Advanced Indexing にあります。slicingだと規則正しい

A0,1と一見ほぼ同等の結果を、以下のようにslicingを使って得ることもできます。

>>> A[0:2,:]
array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]])

ただし、advanced Indexingでは常にコピーが返るのに対し、slicingではビューが返るという違いがあるので注意です。

0行目と1行目を入れ替える、slicingを使った別解を以下に示します。

>>> A = np.arange(25).reshape(5,5)
>>> A[0:2,:] = A[1::-1,:]
>>> A
array([[ 5,  6,  7,  8,  9],
       [ 0,  1,  2,  3,  4],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])

時間を計測すると、advanced indexingよりslicingのほうが高速でした。これがコピーとビューの差でしょうか?

>>> timeit.timeit("A[0:2,:] = A[1::-1,:]", setup="import numpy as np; A = np.arange(25).reshape(5,5)")
0.7406019989994093
>>> timeit.timeit("A[[0,1]] = A[[1,0]]", setup="import numpy as np; A = np.arange(25).reshape(5,5)")
3.9081818149988976

73. Consider a set of 10 triplets describing 10 triangles (with shared vertices), find the set of unique line segments composing all the triangles

告白すると問題の意味がよく分かりませんでした。問題文には「10個の三角形を構成する10組のtriplet」を考えよ、と述べていると思いますが、想定解

faces = np.random.randint(0,100,(10,3))

で生成している10組の数値はランダムで、当然ながら三角形を構成しないtripletも生成し得ます。

その後の"the set of unique line segments composing all the triangle"の意味もよく分かりませんでした。想定解から察するに、三角形の各角を構成する2辺を、重複を排除してprintしているようです。

75. How to compute averages using a sliding window over an array?

あるarrayが与えられたとき、arrayの上で固定ウィンドウを動かし、各位置でのウィンドウ内の数値の平均値を求める問題です。

公式回答は以下になりますが、順を追わないとなかなか理解が難しいです。

def moving_average(a, n=3) :
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n

仮に、a = [A_0, A_1, A_2, ..., A_(K-1)]とします。 np.cumsum(a)は累積和をとる関数で、ret = [S_0, S_1, ..., S_(k-1)]です。ここで、S_i = A_0 + A_1 + ... + A_iです。

ret[n:] - ret[:-n]と、ちょうど累積和の数列をn個ずらして差をとることで、数列Aの隣り合うn項の和を順に並べたものを生成できています。 より具体的に書くと、

  • ret[n:][S_n, S_(n+1), ..., S_(K-1)]
  • ret[:-n][S_0, S_1, ..., S_(K-n)]

なので、

  • ret[n:] - ret[:-n][S_n - S_0, S_(n+1) - S_1, ..., S_(K-1) - S_(K-n)]

です。

ただ、最後の計算で得られた配列には、その最初に来るべきA_0 + A_1 + ... + A_(n-1)が抜けています。A_0 + A_1 + ... + A_(n-1)ret[n-1]に等しいので、ret[n-1]に、最後の計算で得られた配列ret[n:] - ret[:-n]を連結すると、無事数列aの上で固定ウィンドウを動かしてウィンドウ内の和をとった数列が求まります。あとはnで割れば平均になります。