NumPy221本ノックの本書いた

1～2ヶ月ぐらいかけてNumPyの基礎～実装力をつけるための本を書きました。問題演習形式で221問収録してあります。

着想

以前この本を書いたのですが、ディープラーニングに行く前のNumPyでつまずいている方が結構多かったのです。例えば「np.meanのaxisってどうするの？」や「shapeが何だかわからない」というケース。そこをどうにかケアしようというのが今回の目的です。

NumPyの問題演習は昔からありまして、100 numpy exercisesが有名です。昔初心者の頃自分もやったことがあります。最初の方は良いのですが、あとの方が「これイマイチ何に使うんだろう」という疑問があったのです。今振り返ってみても「NumPyの問題としてはよくても、そこまで実践的ではないな」と多少不満点はありました。「それならもう自分で作ってしまおう」ということで本を書きました。

大事にしたところ

本書を書く上で大事にしたところがいくつかあります。

関数暗記で終わらせない

NumPyの関数を英単語のように覚えても、ほとんど実装できるようにはならないと個人的には思います。そこで、覚える関数は最小限にして、「できることがこんなにあるんだよ」ということを演習問題の中で体感してもらうという点を重視しました。必要な関数は適宜ヒントとして入れる形式にしました。

ディープラーニングで使う関数を多く織り交ぜる

NumPyをはじめ数値計算ライブラリには関数が多くありますが、実際ディープラーニングでよく使う関数は限られています。ディープラーニングではNumPyではなく、tf...やtorch...のようなフレームワークの関数になりますが、NumPy関数と同じとなるようにフレームワーク側で設計されています。フレームワークとNumPyは密接な関係にあります。読者がNumPyからディープラーニングにステップアップしたときに、「NumPyではこう書くよな」というが頭に思い浮かぶのを理想形としました。

機械学習の要素を消す

「ディープラーニングで使う関数を多く織り交ぜる」と書いた一方で、本書では一切機械学習の要素は出てきません。これはNumPyがわからない方に追加で機械学習を教えると負荷が大きすぎるからです。ただし、機械学習やディープラーニングの実装で頻出するエッセンスだけ、NumPy関数で織り込みました。

その代わりに、データ処理やシミュレーション、Photoshop・GIMPに代表される画像編集ソフトの処理を配列計算で実装しています。配列計算という意味では、ディープラーニングと画像編集ソフトに共通点も多いです。最終章の動画編集は、最先端のディープラーニングのモデルを実装すると出てくるような黒魔術を出しています。NumPyを知っている方でも満足できると思います。

演習問題の自身の解答を自由にブログ等に貼ってよい

詳細はこちらの規約を読んでいただきたいのですが、「演習問題を解いたぞ」という勉強結果や、成果物をご自身のブログに貼り付けて構いません。ご自身の成果物に自分で解説をつけるのも構いません。リポジトリには模範解答や本文用のColabノートブックを貼り付けてあるので、リンクをたどりやすくなっています。ただし、問題文を丸写しての公開はご遠慮ください。

章構成

1. Pythonの基礎
2. NumPyの導入
3. 行列の導入
4. テンソルと画像処理
5. ブロードキャストの応用

1章がNumPy導入前のPythonのリスト。2章がNumPy導入をし、ベクトルの処理。3章で行列を導入します。4章がランク3の配列（テンソル）を導入し、画像処理をひたすら行います。5章がランク4になり動画を扱います。5章ではこれに加えてブロードキャストを多用するような問題を扱います。

5章の$(F, H, W, C)$というランク4のテンソルです（$F$はフレーム）。ディープラーニングの静止画のテンソル$(B, H, W, C)$であるため、5章ができればディープラーニングで画像は扱えるようになるという目論見です。

問題の一例

問題の一例を示します。模範解答つき。

初見の要素には、ヒントも入れています。1章はBoothで試し読みできます。

第1章

最初は簡単すぎるぐらいの問題です。後半になると少し実践的な問題になります。

Q1：はじめましてリスト

「こしあん、ごまあん、栗あん、かぼちゃあん」からなるリストをanko_a_nameで定義し、print関数で変数全体を表⽰せよ。

解答

anko_a_name = ["こしあん", "ごまあん", "栗あん", "かぼちゃあん"]
print(anko_a_name)

Q2：リストとforループ

Q1で定義したanko_a_nameをforループに⼊れ、print関数で要素ごとに表⽰せよ。

解答

## Q2
for anko in anko_a_name:
    print(anko)

Q3：⽂字列の結合

Q1で定義したanko_a_nameをforループに⼊れ、各要素の末尾に「好評発売中！」の⽂字を追加して表⽰せよ。

解答

## Q3
for anko in anko_a_name:
    print(anko+"好評発売中！")

Q24：北海道の市町村～データ読み込み～

ここから実際のデータで、北海道の市町村のデータを用いる。データ出典は北海道のホームページから　http://www.pref.hokkaido.lg.jp/ss/tuk/900brr/index2.htm。データは使いやすいように整形してGitHubに上がっているので、インポートのために以下のコードを実行しよう。Colab上での実行を想定している。

!git clone https://github.com/koshian2/numpy_book # 先頭に「!」をつけるのを忘れないように。Linuxコマンドを示す

「numpy_book/data/hokkaido_raw.txt」が今回使うデータだ。以下のようなカンマ区切りのテキストになっている。左から振興局（北海道内の地方）、市町村名、令和2年の人口を示している。これを以下のコードで読み込み、textという変数に代入して表示せよ。結果は行単位のリストになる。

with open("numpy_book/data/hokkaido_raw.txt") as fp:
    text = fp.readlines() # 行単位で区切ってリストにしなさいよという意味

解答

## Q24
!git clone https://github.com/koshian2/numpy_book
with open("numpy_book/data/hokkaido_raw.txt") as fp:
    text = fp.readlines()
print(text)

Q25：北海道の市町村～データ整形～

Q24のtextは行単位のデータだが、これを分析しやすいデータにしたい。最終的に作りたいのは、次のような二重リストの構造だ。

[
    ['空知', '夕張市', '市', 7769],
    ['空知', '岩見沢市', '市', 80410],
    : : : : :
    ['根室', '羅臼町', '町', 4955]
]

このようなデータをhokkaidoとして定義して表示せよ。ただし本書でまだ触れていない関数も使うためヒントを与える。ヒント

1. textの各要素は末尾の「\n」が不要だ。これを消すには、文字列.replace("\n", "")と文字列操作で置き換えるか、インデックスのスライスを使って文字列[:-1]と最後の1文字を取り除いたスライスを取る。
2. カンマ区切りの文字列をリストにするには、文字列.split(",")で分割する。「\n」を取り除いた各行のテキストを分割すると、3つの要素からなるリストができるはずだ。
3. 市名のあとに、「市、町、村」かどうかのフラグを追加したい。これは市町村名の最後の1文字を取ることで実現できる（なぜならどの自治体名も市町村のどれかの文字で終わるから）。
4. 人口はsplitの直後だと文字列として与えられている。しかし、これは後々数値に変換したほうが便利だ。文字列として与えられている数値を整数の型（詳細は次章以降）に変換するには、int(文字列の変数)とする。「整数にすると大小比較がやりやすい」ぐらいに今は理解しておけばOKだ。

解答

## Q25
hokkaido = []
for line in text:
    item = line[:-1].split(",") # itemは['空知', '夕張市', '7769']のようになっている
    item.insert(2, item[1][-1]) # 挿入箇所は人口の前、挿入する要素は市町村名の最終文字
    item[-1] = int(item[-1]) # 人口を数値に
    hokkaido.append(item) # 全体のリストへの追加を忘れずに
print(hokkaido)

Q26：北海道の市町村〜条件判定〜

（Q25で定義した）hokkaidoから、振興局名（各要素の1番⽬）が「⽯狩」のものだけ表⽰せよ。⾃治体単位のリストをprintで表⽰するだけで良い（ヒント参照）。

ヒント：「⽯狩のものだけ」の判定は次のようにする。等しいかどうかの判定はイコールが2個（==）であることに注意せよ。代⼊（=）とは異なる。

解答

## Q26
for city in hokkaido:
    if city[0] == "石狩":
        print(city)

第2章

2章はNumPy配列を使った初歩の問題が多いです。ただし最後の問題は簡単なゲームの戦闘のシミュレーションなのでかなり難しいです。難しい問題は（応用）とつけて、初心者が飛ばしても良いようにしています。

Q1：全て0と1

NumPyをインポートし、「要素数が3で全て1」「要素数が7で全て0」のNumPy配列をそれぞれ定義して表⽰せよ。

解答

## Q1
import numpy as np
print(np.ones(3))
print(np.zeros(7))

Q2：整数のNumPy配列

Q1と同じ配列をint32型で定義し、変数全体とそのdtypeを表⽰せよ。

解答

## Q2
x = np.ones(3, dtype=np.int32)
print(x, x.dtype)
x = np.zeros(7, dtype=np.int32)
print(x, x.dtype)

Q3：中間地点に代⼊

「要素数が5で全て0」のint32型のNumPy配列を定義し、中間地点にある1要素だけ1を代⼊して表⽰せよ。

解答

## Q3
x = np.zeros(5, dtype=np.int32)
x[2] = 1
print(x)

Q4：奇数と偶数のインデックス（1）

「奇数番⽬のみ1、偶数番⽬は0、要素数10」のint32型のNumPy配列を定義し表⽰せよ（本章で説明した内容でもいくつか書き⽅がある）。

解答

## Q4 Ans 1
x = np.zeros(10, dtype=np.int32)
x[np.arange(10) % 2 == 0] = 1
print(x)

## Q4 Ans 2
x = np.zeros(10, dtype=np.int32)
x[::2] = 1
print(x)

よくありがちな問題ですよね。

Q19（応用）：魔王と勇者のデータサイエンス

Q19からは応用問題の中でもかなり難しいので、難しいと思ったらパスして構わない。これまでのNumPyの内容をフル活用したシミュレーションだ。

あるところに魔王がいた、その魔王を倒すべく大量の勇者が戦いを挑む。その戦いの経過を見てみよう。戦いの詳細は以下のとおりだ。

• 魔王は1人、勇者は任意の数（num_heros人）派遣できる。魔王を倒せば勝利、勇者が全滅すれば敗北だ。戦いはRPG方式のコマンドバトルである。
• 勇者か魔王が倒れるまで戦うので引き分けは存在しない
• 勇者について：
  • どの勇者も必ず初期のHPは100である。HPが0になったら死亡する。一度死亡したら蘇生はできない
  • 勇者の行動は3パターンある。通常攻撃、捨て身攻撃、回復である
  • 勇者1人につき、通常攻撃は魔王に常に100のダメージを与える。捨て身攻撃は常に150ダメージを与えられるが、自身に20ダメージを負う
  • 勇者が1人回復行動を取ると、生存している勇者全員のHPが1回復する。3人が回復行動を取ると全員のHPが3回復する。死人のHPは回復しない
  • 回復があったとしても、HPが100を超えることはない
• 魔王について：
  • HPは100万であり、常に1人である
  • 1ターンにつき複数（satan_attack_freq）回行動できる
  • 魔王は通常攻撃か、範囲攻撃の2種類の行動パターンがある
  • 通常攻撃は、勇者1人に対してsatan_base_powerダメージを与える。範囲攻撃は勇者3人に対してsatan_base_power // 2ダメージ与える。例えば、satan_base_powerが100なら、通常攻撃のダメージは100（当たれば即死）、範囲攻撃のダメージは50×3人となる
• 常に「勇者から魔王への攻撃」「魔王から勇者への攻撃」の順で行動する。魔王の攻撃が終わった時点で1ターン経過とする。勇者の攻撃で魔王のHPが0になったら、そのターンの経過はカウントしない
• 各勇者の行動パターンは次の通りである。全ての勇者は共通の行動ルーチンを持つ：
  • ターン開始時にHPが60以上なら、常に捨て身攻撃をする（捨て身攻撃による自爆は発生しない）
  • ターン開始時にHPが0より大きく20以下なら、常に回復行動を取る
  • それ以外なら通常攻撃
• 魔王の行動パターンは次のとおりである
  • 通常攻撃→範囲攻撃→通常攻撃→範囲攻撃　と交互に行動する
  • 1ターンに複数回行動できるので、通常攻撃が先にくることもあれば、範囲攻撃が先にくることもある。魔王の行動のカウンターは、ターンを問わず戦闘を通じて持ち越して計算する
  • どの勇者を狙うかは、生存している勇者か先頭（配列的な意味で）を攻撃ターゲットとする。もし「1、 2、4、7、9」の勇者が生きているとすれば、通常攻撃なら「1」範囲攻撃なら「1、2、4」を狙う
  • 範囲攻撃で生存している勇者の数が3人未満なら、ターゲット数もそれに応じて減少する（範囲攻撃で1人に複数回攻撃がいかない）
• 魔王の行動カウンター、ターン数のカウントいずれも初期値0とする

def battle(num_heros, satan_base_power, satan_attack_freq):
    ## 戦闘処理
    ## : : :
    ##
    return 生存している勇者の数, ターン数

という関数を定義せよ、というのがQ19の問題である。生存している勇者の数が1人以上いれば勝利、0なら敗北を意味する。（以下詳細なヒントがあるが長いので省略。解答も省略）

Q20（応用）：魔王と勇者のデータサイエンス～最小討伐人数～

あなたは魔王である。satan_attack_freq=1, satan_base_power=80で、勇者の人数を1人ずつ増やしていったとき、何人以上いるとはじめて勇者側が勝利するのか計算せよ。以降、初めて勇者側が勝てる勇者の人数を「最小討伐人数」と呼ぶことにする。

解答

## Q20
n_heros = 1
while(True):
    result = battle(n_heros, 80, 1)
    if result[0] != 0:
        print("最小討伐人数", n_heros, "結果", result)
        break
    else:
        n_heros += 1

Q23（応用）：魔王と勇者のデータサイエンス～パワーかスピードか～

魔王は手持ちの宝石を使って自分の能力を強化できる。初期の能力はsatan_attack_freq=1, satan_base_power=10である。satan_attack_freqは「+1」ごとに、satan_base_powerは「+10」ごとにそれぞれ宝石を1個消費する。攻撃力+5のような端数の強化はできない。宝石は手元に20個あり、これを全て使うものとする。どのように能力強化を配分すれば、勇者側の勝利を最も遠ざけることができるか考えよ。ただし、ここで最適化すべきは「最小討伐人数の最大化」である。

解答

## Q23
min_winable_heros_counts = []
# 先に攻撃回数を固定する
freqs = np.arange(1, 22)

for freq in freqs:
    # 攻撃回数に振ったら残りの宝石は全て攻撃力に振る
    power = (21 - freq) * 10 + 10 
    n_heros = 1
    while(True):
        result = battle(n_heros, power, freq)
        if result[0] != 0:
            print("攻撃回数", freq, "基礎攻撃力", power, "最小討伐人数", n_heros, "結果", result)
            min_winable_heros_counts.append(n_heros)
            break
        else:
            n_heros += 1

plt.plot(freqs, min_winable_heros_counts)
plt.show()

第3章

行列を導入するための基本的な問題にはじまり、集約関数のaxisや行列計算の応用が中心です。コードで書くと簡単なのでここでは解答は省略します。

Q1：np.concatenate

次の各コードにある変数xのshapeを考えよ。エラーになる場合もあるので、そのときはエラーと答えよ。どういう結果になるか想像してからコードを実行せよ。

## Q1-1
a = np.arange(5)
b = np.arange(3)
x = np.concatenate([a, b])

## Q1-2
a = np.arange(5)
b = np.arange(1)
x = np.concatenate([a, b])

## Q1-3
a = np.arange(5)
b = 0
x = np.concatenate([a, b])

## Q1-4
a = np.arange(5)
b = np.arange(7)
x = np.concatenate([a, b], axis=0)

#（以下略）

Q2：np.stack

次の各コードにある変数xのshapeを考えよ。エラーになる場合もあるので、そのときはエラーと答えよ。どういう結果になるか想像してからコードを実行せよ。

## Q2-1
x = np.stack([np.arange(5) for i in range(2)], axis=0)

## Q2-2
x = np.stack([np.arange(4) for i in range(3)], axis=1)

## Q2-3
x = np.stack([np.arange(7) for i in range(5)])

## Q2-4
x = np.stack([np.arange(i+2) for i in range(4)], axis=0)

#（以下略）

Q7：集約関数のaxis

次のデータ$x$に対して、適切に集約関数を適用しなさい。全てのデータに欠損値はないものとする。

1. ある製品の品質チェックを行っている。チェック項目は10項目あり、この数値をもとに正常品かどうかを判定したい。製品3000個のテスト数値があり、x.shape=(3000, 10)で得られた。全製品間の各チェック項目の平均値を出すにはどうすればいいか。
2. あるテストの採点をしている。テストは5教科あり、40人の採点データがx.shape=(40, 5)で得られた。1人の全教科の合計点を出すにはどうすればいいか。
3. ある日の全国の気温データを扱っている。観測点は100点あり、1時間毎の気温データがx.shape=(100, 24)で得られた。各観測点の最高気温、最低気温を出すにはどうすればいいか。
4. 3において、全観測点の中での最高気温を出すにはどうすればいいか。
5. 3において、全観測点の中での最低気温を出した観測点のインデックスを出すにはどうすればいいか（インデックスは0、1、2という数値でよい）。

（以下省略）

Q22：あんこ商店

「numpy_book/data/anko_shop.npz」にはある個人商店の売上データが入っている。この商店には4つのものが売っている。

• アンパン（こしあん）→税抜：120円、税込：130円
• 大福（栗餡）→税抜：180円、税込：194円
• 缶ビール→税抜：270円、税込：297円
• 殺虫剤→税抜：1200円、税込：1320円

消費税は2020年基準とする。アンパンと大福は軽減税率対象（消費税8％）だが、缶ビールと殺虫剤は軽減税率対象ではない（消費税10％）。商品リストを行列として定義し、売上データとの行列計算1回で、1日単位の「税抜売上」「税込売上」を計算しなさい。ただし、売上データは以下のコードで読める。

anko_sell = np.load("numpy_book/data/anko_shop.npz")["sell"]

ヒント：全期間の売上を合計すると、税抜287970円、税込313214円となる。

第4章

章末問題だけで100本ノックです。画像処理が中心です。

Q1：ベタ塗りグレースケール

解像度128×128, 256×256, 512×512のグレースケール画像を作って表示せよ（表示の方法は好みで良い）。背景はそれぞれ黒（0）、50%灰色（128）、75%灰色（192）で塗りつぶそう。

解答

## Q1
fig = plt.figure(figsize=(16, 8))
ax = fig.add_subplot(1, 3, 1)
img = np.full((128, 128), 0, dtype=np.uint8) 
ax.imshow(img, cmap="gray", vmin=0, vmax=255)

ax = fig.add_subplot(1, 3, 2)
img = np.full((256, 256), 128, dtype=np.uint8) 
ax.imshow(img, cmap="gray", vmin=0, vmax=255)

ax = fig.add_subplot(1, 3, 3)
img = np.full((512, 512), 192, dtype=np.uint8) 
ax.imshow(img, cmap="gray", vmin=0, vmax=255)

Q2：四角形の塗りつぶし

横640×縦480でグレースケール、白背景の画像を作れ。以下の領域を塗り潰して表示せよ（右下の点を長方形の幅や高さに含めるかは好みで良い）。

• 左上が$(x, y) = (100, 100)$、横幅200、縦幅100の長方形。塗りつぶし色は黒（0）
• 左上が$(x, y) = (400, 200)$、横幅150、縦幅250の長方形。塗りつぶし色は50%灰色（128）

解答

## Q2
img = np.full((480, 640), 255, dtype=np.uint8)
img[100:200, 100:300] = 0
img[200:450, 400:550] = 128
plt.imshow(img, cmap="gray", vmin=0, vmax=255)
plt.show()

Q5：平方数の直線

横640×縦480でグレースケール、白背景の画像を作れ。以下の数列を満たす$x, y$の座標に線を引け。

$$1, 4, 9, 16, 25, 36, 49, 64, 81, 100, \cdots$$

つまり、$x=1, y=1, x=4, y=4, ...$という直線を引く。直線の太さは1、色は黒とする。

解答

## Q5
img = np.full((480, 640), 255, dtype=np.uint8)
for i in range(int(np.sqrt(img.shape[0]))):
    img[i**2] = 0
for i in range(int(np.sqrt(img.shape[1]))):
    img[:, i**2] = 0

plt.imshow(img, cmap="gray", vmin=0, vmax=255)
plt.show()

Q20: カラー画像入門

縦400、横600のカラー画像のキャンバスを作れ。横0～200を赤、横200～400を緑、横400～600を青で塗りつぶして表示しなさい。結果はR,G,Bの三色の国旗のようになる。

解答

## Q20
img = np.zeros((400, 600, 3), dtype=np.uint8)
img[:, 0:200, 0] = 255
img[:, 200:400, 1] = 255
img[:, 400:600, 2] = 255
plt.imshow(img)
plt.show()

Q30（応用）：セーシェルの国旗

以下の旗はインド洋に浮かぶ、アフリカの国・セーシェル共和国の国旗である。縦横比は1：2である。それ以外の正確な仕様は見つけられなかったが、以下のような特徴があった。この情報をもとに、横幅が1200ピクセルの「セーシェルの国旗」をNumPy配列で作成し、表示せよ。これを「seychelles.png」として保存せよ。

解答

## Q30
img = np.zeros((600, 1200, 3), dtype=np.uint8)
y_ind, x_ind = np.indices(img.shape[:-1])
aspect_ratio = -img.shape[1] / img.shape[0]
lines = [y_ind - (aspect_ratio / 3 * x_ind * 2 + img.shape[0]),
         y_ind - (aspect_ratio / 3 * x_ind  + img.shape[0]),
         y_ind - (1 / aspect_ratio / 3 * x_ind * 2 + img.shape[0]),
         y_ind - (1 / aspect_ratio / 3 * x_ind + img.shape[0])]

img[lines[0]<=0] = [0, 47, 108] # 青
img[lines[-1]>=0] = [0, 122, 51] # 緑
colors = [[254, 209, 65], [210, 39, 48], [255, 255, 255]] # 中間色
for i, c in zip(range(len(lines)-1), colors):
    img[np.logical_and(lines[i]>=0, lines[i+1]<=0)] = c

with Image.fromarray(img) as fp:
    fp.save("seychelles.png")

plt.imshow(img)
plt.show()

Q82（応用）：禍々しいハート

Q81の円形グラデーションを使いまわし、alphaの値を0.65つまり、$v=clip(1-0.65r)$とする。ここでclipは0-1のクリップである。今度の下地は、今までの問題で上乗せ画像として使っていた虹色のハートref_imgである。ref_imgとvをリニアライトで合成しプロットせよ。かなり禍々しいハートが出てくるはずだ。

解答

## Q82
y_ind, x_ind = np.indices(base_img.shape[:-1])
r = np.sqrt((y_ind-base_img.shape[0]//2)**2 + (x_ind-base_img.shape[1]//2)**2) / base_img.shape[0] * 2
v = np.expand_dims(np.clip(1-0.65*r, 0, 1), axis=-1)

a = ref_img + 2*v - 1
b = ref_img + 2*(v-0.5)
out = (v<0.5)*a + (v>=0.5)*b
out = np.clip(out, 0, 1)
out = (out*255.0).astype(np.uint8)
plt.imshow(out)
plt.show()

第5章

5章はブロードキャストがメインです。動画の問題もあります。最後の方はマニア向けなので、かなり黒魔術です。

Q1：ブロードキャスト後のshape

次のコードのxのshapeを予想しなさい。必要ならコードを実行して確認しなさい。

## Q1-1
x = np.random.rand(1, 5) * np.random.rand(3, 5)
## Q1-2
x = np.random.rand(1, 7) * np.random.rand(4, 1)
## Q1-3
x = np.random.rand(3, 5) * np.random.rand(5)
## Q1-4
x = np.random.rand(2, 1, 3) * np.random.rand(2, 4, 1)
## Q1-5
x = np.random.rand(1, 1, 2) * np.random.rand(4, 3, 2)
## Q1-6
x = np.random.rand(1, 5, 2) * np.random.rand(5, 1, 2)
## Q1-7
x = np.random.rand(3, 7, 1, 8) * np.random.rand(1, 7, 6, 1)
## Q1-8
x = np.random.rand(2, 8, 5) * np.random.rand(8, 5)

実行すればわかるのでここでは解答省略

Q3：掛け算九九の表

本章の内容を参考に「掛け算九九の表」の1～9の段までを、ブロードキャストを用いて求めなさい。余力があればワンライナーで実行すること。

解答

## Q3
print(np.arange(1, 10).reshape(-1, 1)*np.arange(1, 10).reshape(1, -1))

Q5：駅間距離

5駅からなる鉄道があり、始発駅からの累計距離（km）は「0, 1.1, 2.8, 4.7, 7.0」で表される。このとき、駅同士の距離を5×5の行列で示せ。

ヒント：ブロードキャストして差の絶対値を取る。

解答

## Q5
dists = np.array([0, 1.1, 2.8, 4.7, 7.0])
print(np.abs(dists[:, None]-dists[None, :]))

Q12（応用）：動くゲーミングあんこちゃん

この問題は4章の演習問題「Q77：ゲーミングあんこちゃん」を動画に拡張した問題である。Q11のグラデーションにオーバーレイ合成を入れる。背景$B$, 参照画像の$R$という2つの画像が出てくる。

背景の$B$から。「numpy_book/imgs/girl01.png」（RGBA画像）を読み込み、RGB色空間に変換し、0-1スケールにしたものを使う。この背景画像は全フレームで使いまわしする（したがってブロードキャストの要素が出てくる）。

参照画像の$R$について。Q11と似たように作るが、$B$と解像度は同じにする。縦方向に色相を$4i\to 4i+256$と変化させたグラデーションとする。これを300フレーム作る。また、オーバーレイの計算式は次の通りである。

$$\begin{cases}2BR & if B<0.5 \\ 1-2(1-B)(1-R) & otherwise \end{cases}$$

オーバーレイは4階テンソル同士のオーバーレイ計算を想定している。この計算式でブロードキャストを使えば、複数の画像のオーバーレイも同時に可能だ。手順としては、$B$を読み込む→$R$のフレームを作る→オーバーレイ計算をする→動画化するになる。

ただし、この問題、メモリ消費量に気をつける必要がある。なぜなら、解像度700×700近いfloatの画像を300枚近く作る必要があるからである。メモリ消費量を落とすためには次のように気をつけると良い。

• uint8の画像の配列を、255.0で割ると小数になるが、割る直前でnp.float32にキャストする。キャストしないとfloat64に勝手に変換される。キャストすることでメモリ消費量を半減できる（小数の精度低下は気にしなくて良い）。
• メモリ消費量が大きい配列（特にグラデーションのフレームを作った以降の画像）は、変数名を使い回すようにする。Jupyter Notebookでは変数に再代入されないとメモリが開放されないため、メモリ消費量が大きい配列の変数名を使い回すことで、メモリ内にゴミが散らからないようにする。
• どうしてもメモリが一杯一杯になったら「ランタイムの再起動」（ランタイム→ランタイムを再起動）をする。
• それでもメモリが溢れてしまうのなら、画像サイズを小さくしたり、フレーム単位でオーバーレイを計算をする。

ネタ pic.twitter.com/LlRHU4SRME

— 🔞こあん＠日高本、モザイク本好評通販中 (@koshian2) November 10, 2020

解答

## Q12
# 背景画像
with Image.open("numpy_book/imgs/girl01.png") as img:
    base_img = np.array(img.convert("RGB")).astype(np.float32) / 255.0
    base_img = base_img[None, :, :, :] # (1, H, W, 3)

# グラデーション
frames = []
for i in range(300):
    hsv = np.full(base_img.shape[1:], 255, dtype=np.uint8)
    hsv[:,:,0] = np.linspace(i*4, 256+i*4, num=hsv.shape[0], dtype=np.uint8)[:, None]
    with Image.fromarray(hsv, "HSV") as img:
        frames.append(np.array(img.convert("RGB")))
frames = np.stack(frames, axis=0).astype(np.float32) / 255.0 # (512, H, W, 3)

# オーバーレイ
frames = (base_img<0.5)*(2*base_img*frames) + (base_img>=0.5)*(1-2*(1-base_img)*(1-frames))
frames = (frames*255.0).astype(np.uint8)

# 動画にして保存
video_utils.save_video("videos/q12.mp4", frames)
moviepy.editor.VideoFileClip("videos/q12.mp4").ipython_display()

続きはこちらから

221問もあるので全部解こうとはせず、楽そうなところや、面白そうだなと思うところから解くと負担にならないはずです。自信のある方は4章あたりから解くとちょうどいいはずです。

ぜひお正月の暇なときにテレビやYouTubeでも見ながら解いてみてください。2020/12/26～開催される技術書典10で頒布します。詳しくはこちらから↓（オンラインマーケットは12/26～。Boothや委託書店はもうあります）

• Booth：https://koshian2.booth.pm/items/2462894
• 技術書典オンラインマーケット：https://techbookfest.org/product/5547509835366400
• その他まとめ：https://github.com/koshian2/numpy_book
• ISBN：978-4-910088-26-6

※Amazonでの物理書籍の取扱も調整中です