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で割れば平均になります。

100 numpy exercisesの解説 1~50

NumPyでよく使うテクニックが網羅されたnumpy-100/100 Numpy exercises.ipynb at master · rougier/numpy-100 · GitHubを一通りやってみています。以下はその途中で調べたメモ書きです。自分が知らなかった機能や、一見して理解が難しい問題に絞ってメモしています。断らない限り64bit Python 3.6 + numpy 1.11で試しています。

問題数が多いので、この記事では50問目までを対象とします。51問目~75問目は 100 numpy exercisesの解説 51~75 - minus9d’s diary、76問目~100問目は 100 numpy exercisesの解説 76~100 - minus9d’s diary をご参照ください。

20. Consider a (6,7,8) shape array, what is the index (x,y,z) of the 100th element?

unravelは、「〈もつれた糸・編み物などを〉解く」という意味の動詞です(〈もつれた糸・編み物などを〉解く)。

26. What is the output of the following script?

python組み込みのsum()関数は、sum(iterable[, start])という引数を持ちます。これは以下のように計算されます。(参考:python sum function - `start` parameter explanation required - Stack Overflow)

def sum(values, start = 0):
    total = start
    for value in values:
        total = total + value
    return total

なので、print(sum(range(5),-1))を計算すると -1 + 0 + 1 + 2 + 3 + 4 = 9となります。

一方、from numpy as npすると、python組み込みのsum()関数の代わりに、np.sum()が呼ばれます。np.sum()の引数をnumpy.sum — NumPy v1.12 Manual で調べると、この問題の最終行は

print(np.sum(range(5), axis=-1))

と等価であることがわかります。

axisは何番目の軸を指定するかを意味します。負の値をいれた場合は、軸を逆順で数えます。今回の場合、range(5)は1次元なので、-1の指定は0の指定と同じです。結局、0 + 1 + 2 + 3 + 4 = 10となります。

27. Consider an integer vector Z, which of these expressions are legal?

Z <- Zが最初理解不能でした。<-という演算子があるわけではなく、Z < -Zをわざと変なスタイルで書いているだけだと思います。

>>> a = np.array([0,-3,4])
>>> a <- a
array([False,  True, False], dtype=bool)

28. What are the result of the following expressions?

まずnp.array(0)が理解できませんでした。np.array([0])であれば長さ1のベクトルですが、np.array(0)とは何なのでしょうか?

実はnp.array(0)は「0という値を持つ0次元のarray」を表すようです。コンソールで値をいじくってみると性質をつかめると思います。(参考:numpy - error extracting element from an array. python - Stack Overflow

>>> a = np.array(7)
>>> type(a)
<class 'numpy.ndarray'>

# 普通のndarrayが持つ関数はすべて使える
>>> a.ndim
0
>>> a.dtype
dtype('int64')
>>> a.min(), a.max()
(7, 7)

# 0次元なのでインデックスでアクセスはできない
>>> a[0]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
IndexError: too many indices for array

# item()関数を使うと値を取り出せる
>>> a.item()
7
>>> type(a.item())
<class 'int'>

0次元のarrayというのは数学的にはスカラーだと思うのですが、NumPyにおける0次元のndarrayは、NumPyにおけるscalarとは別物のようです。

>>> np.isscalar(np.array(0))
False
>>> np.isscalar(np.int_(4))
True

この問題に関しては python - Why are 0d arrays in Numpy not considered scalar? - Stack Overflow に詳しいですが、自分は深入りしていません。0-dimension array問題 - Qiita にも0次元のndarrayに関する詳細があります。

以上を踏まえて、28. を試してみます。実はPython 3とPython 2で挙動が変わります。まずはPython 3.6の場合。

>>> print(np.array(0) / np.array(0))
__main__:1: RuntimeWarning: invalid value encountered in true_divide
nan
>>> print(np.array(0) // np.array(0))
__main__:1: RuntimeWarning: divide by zero encountered in floor_divide
0

このRuntimeWarningは初回実行時のみ表示されるもので、同じ演算を2回目以降行っても表示されません。

次にPython 2.7 + numpy 1.21で試した場合。from __future__ import divisionは設定していません。

>>> print(np.array(0) / np.array(0))
__main__:1: RuntimeWarning: divide by zero encountered in divide
0
>>> print(np.array(0) // np.array(0))
__main__:1: RuntimeWarning: divide by zero encountered in floor_divide
0

ちょっと自分には理解不能な結果です。

3つ目はPython2, Python3ともに以下の結果になりました。

>>> print(np.array([np.nan]).astype(int).astype(float))
[ -9.22337204e+18]

np.float64型であるNaNをnp.int64型として解釈し、再びnp.float64型として解釈するとNaNに戻らないことを示した例なのですが、残念ながら何が起こっているのかここに書けるほどのことは分かりませんでした。

29. How to round away from zero a float array ?

“round away from zero”, つまり0から遠い方向に小数点を丸める方法について問うています。例えば、2.2は3に丸め、-7.4は-8に丸めることになります。 いったんndarrayの要素をすべて正の数にして、np.ceil()で0から遠い方向に丸めたあと、copysign()で元の符号を復活させています。

32. Is the following expressions true?

np.sqrt(-1)nanになり、np.emath.sqrt(-1)1jになります。ここで、jは複素数を表す記号です。np.emath複素数を扱うときに使うようですが、emathが何の略かはわかりませんでした。

38. Consider a generator function that generates 10 integers and use it to build an array

np.fromiter()は、イテレータからndarrayを生成する関数です。引数のcountは、要素数の上限です。

41. How to sum a small array faster than np.sum?

実際に試してみました。確かにnp.add.reduce()のほうがnp.sum()より高速です。

# setupで指定した部分は最初の1回のみ実行される。最初の引数は100万回実行される。
>>> timeit.timeit("np.add.reduce(a)", setup="import numpy as np; a = np.arange(10);")
0.7812885560124414
>>> timeit.timeit("np.sum(a)", setup="import numpy as np; a = np.arange(10);")
1.9236456239887048

python - What is the difference between np.sum and np.add.reduce? - Stack Overflow によると、``np.sum()は内部でnp.add.reduce()を呼んでいるので、オーバーヘッド分だけnp.sum()のほうが遅くなるようです。

43. Make an array immutable (read-only)

writableではなくwriteableなのに注意。

44. Consider a random 10x2 matrix representing cartesian coordinates, convert them to polar coordinates

R = np.sqrt(X**2+Y**2)の部分はR = np.hypot(X, Y)としたほうが簡潔です。

47. Given two arrays, X and Y, construct the Cauchy matrix C (Cij =1/(xi - yj))

C = 1.0 / np.subtract.outer(X, Y)とありますが、実際は添字の順序を考えるとC = 1.0 / np.subtract.outer(Y, X)が正しいような気がします(違っていたらすみません)。

このouterという関数は、ベクトルAに含まれる要素aと、ベクトルBに含まれる要素bとの全ペアに対して演算を施すときに使う関数のようです。(参照:numpy.ufunc.outer — NumPy v1.12 Manual

Pythonのitertoolsでできる全列挙のまとめ

Pythonitertoolsモジュールには、イテレーションに関する便利関数が多数用意されています。この記事では、その中でも競技プログラミングで全列挙に使える関数についてまとめます。Python 2.x, 3.xのどちらでも使えます。

Google Code JamのSmallでは全列挙するだけで解ける問題が出ることが多いです。一通りどんな関数があるか知っておきましょう。

import文

スクリプトの冒頭に import itertoolsと書いてあるものとします。

1, 2, 3と書かれた玉が入った袋から、2つの玉を取り出す (3P2の順列)

nums = [1, 2, 3]
for balls in itertools.permutations(nums, 2):
    print(balls)

結果は以下です。

(1, 2)
(1, 3)
(2, 1)
(2, 3)
(3, 1)
(3, 2)

ここで、第2引数を省略すると3P3 (= 3!)が列挙されます。C++std::next_permutation()でできることと同じです。

nums = [1, 2, 3]
for balls in itertools.permutations(nums):
    print(balls)

結果は以下です。

(1, 2, 3)
(1, 3, 2)
(2, 1, 3)
(2, 3, 1)
(3, 1, 2)
(3, 2, 1)

1, 2, 3と書かれた玉が入った袋から、2つの玉を区別なく同時に取り出す (3C2の組み合わせ)

nums = [1, 2, 3]
for balls in itertools.combinations(nums, 2):
    print(balls)

結果は以下です。

(1, 2)
(1, 3)
(2, 3)

1, 2, 3と書かれた玉が入った袋から、2つの玉を取り出す。ただし取り出した玉はまた袋に戻す (重複順列)

nums = [1, 2, 3]
for balls in itertools.product(nums, repeat=2):
    print(balls)

結果は以下です。

(1, 1)
(1, 2)
(1, 3)
(2, 1)
(2, 2)
(2, 3)
(3, 1)
(3, 2)
(3, 3)

1, 2, 3と書かれた玉が2セット入った袋から、2つの玉を区別なく同時に取り出す(重複組合せ)

nums = [1, 2, 3]
for balls in itertools.combinations_with_replacement(nums, 2):
    print(balls)

結果は以下です。

(1, 1)
(1, 2)
(1, 3)
(2, 2)
(2, 3)
(3, 3)

2つの集合の直積をとる

signs = ['a', 'b', 'c']
nums = [1, 2, 3]
for pairs in itertools.product(signs, nums):
    print(pairs)

結果は以下です。

('a', 1)
('a', 2)
('a', 3)
('b', 1)
('b', 2)
('b', 3)
('c', 1)
('c', 2)
('c', 3)

kaggleの練習コンペDogs vs. Catsに出場

Kaggleで開催されていた練習用コンテスト Dogs vs. Cats Redux: Kernels Edition | Kaggle に参加しました。このコンテストでは、画像の中に犬または猫のどちらが映っているかを識別する識別器を学習し、その性能を競います。通常では犬か猫かの2値で答える問題設定が多いかと思いますが、このコンペでは犬か猫かの確率値を0.0から1.0の間の実数値で答えます。

結果は1314人中74位で、上位6%に入ったといえば聞こえはいいのですが、実態は手も足も出なかったという印象です。以下、簡単なメモです。

実装1: Chainer

まずはChainerを使ってDeepNetの自力実装を試みました。ディープラーニングの様々なモデルを使ってCIFAR-10画像データセットの分類を行う - Qiita のVGG likeなモデルを参考に、与えられた学習データのみを使ってスクラッチからのモデルの学習を試みました。

入力画像をGTX 1070のメモリが許す限り大きくしたり(144x144の画像からランダムな位置で128x128の画像を切り出し)、on-the-flyで画像に摂動を加えて水増ししたり、出力層をsigmoid cross entropyにして確率値を出せるようにしたりなど、思いつきそうな工夫をあれこれ入れました。しかしエラー値0.13(正解率にして95%程度、順位にして約50%のあたり)程度で足踏みしてしまいました。感覚的にはもう少しよい値が出そうだと感じましたが、これがネットワークの限界なのか、学習方法に何らかの不備があるのかは不明です。

実装2: Keras

Kaggleではコンペ中も掲示板形式で参加者同士の議論が公開されています。その中に、Kerasの70行ほどのコードでエラー値0.06を出すサンプルが掲載されていました(Dogs vs. Cats Redux: Kernels Edition | Kaggle)。このコードでは、あらかじめ学習済のネットワークを特徴量抽出器として流用し、犬猫判別の判別器部分のみを学習する、いわゆる転移学習が行われています。

このコードを使うと、あっさりとエラー値0.059が出ました。転移学習恐るべしです。さらに小手先の閾値調整(値の範囲をクリッピングして、誤分類時のペナルティを下げるなど)を実施して、エラー値が0.05164(正解率にして99%程度)まで下がりました。

感想

学習データが数万枚オーダーで存在したためスクラッチで学習してもそれなりの精度が出ることを期待しましたが、実際は転移学習によるサンプルコードの性能に遥か及びませんでした。

今回はほぼサンプルコードを動かしただけで終わってしまいましたが、次回は最低一つくらいは何らかの工夫を加えてサンプルコードよりよい成績を狙いたいと思います。また、Discussionをくまなくチェックし、素早く巨人の肩に乗る練習も行いたいと思います。