菜

周波数ワーピング

インターネットに解説があんまり無くて永遠に理解できてなかったけどやっとわかってきた (本を探せばあったのかな……)　z 変換をちゃんと理解してなかった

参考にしたページ
https://www.rle.mit.edu/dspg/documents/Proc_Frequency_1999.pdf

z 変換

インパルス応答とかの離散信号をフーリエ変換すると正規化角周波数 \omega (単位が rad / sample のやつ) の式 F(\omega) が得られる。 \omega は -\infty から \infty までの数直線上の値 (2\pi ごとにループするけど) である。 数直線上から自由に角周波数 \omega を選んで F(\omega) を計算すると、元の信号がインパルス応答ならその角周波数でのフィルタの振幅特性や位相特性がわかる。

これを、何を思ったか突然 z=\cos\omega + j\sin\omega で変数変換してみたものが z 変換らしい。このとき \omega と z の対応は下のようになる。


\begin{array}{c|ccccccc}
\omega & \cdots & 0 & \frac{\pi}{2} & \pi & \frac{3\pi}{2} & 2\pi & \cdots \\ \hline
z & \cdots & 1 & j & -1 & -j & 1 & \cdots
\end{array}

つまり、 \omega が数直線上を動くのは、 z が単位円上をぐるぐるするのと同じ。 \omega は数直線上しか動かないので z は必ず単位円上の値だと考えていい。(本当はだめかも？とりあえず今は大丈夫) ただの変数変換なので、z 領域の式 X(z) で、 z に j を代入することで正規化角周波数 \omega=\frac{\pi}{2}のときの周波数特性が求められたりする。

周波数ワーピング

\tilde z^{-1}=\frac{z^{-1}-\alpha}{1-\alpha z^{-1}}(|a|\lt1) によって、またしても変数変換をする。(この式はオールパスフィルタそのもの。) z^{-1} について解けば、 z^{-1}=\frac{\tilde z^{-1}+\alpha}{1+\alpha \tilde z^{-1}} となる。

変数変換なので、対応を確認することにする。


\begin{array}{c|ccccccc}
z &  \cdots & 1 & j & -1 & -j & 1 & \cdots \\ \hline
\tilde z & \cdots & 1 & \frac{-2\alpha + j(1-\alpha^2)}{1+\alpha^{2}} & -1 & \frac{-2\alpha - j(1-\alpha^2)}{1+\alpha^{2}} & 1 & \cdots
\end{array}

\frac{-2\alpha + j(1-\alpha^{2})}{1+\alpha^{2}} ってなんだよという感じだが、この値は絶対値を取ると \alpha に関わらず 1 になる。実は、 z が単位円上にあれば、 \tilde z も必ず単位円上にある。これは以下のように確認できる。(この式はオールパスフィルタの振幅特性が平坦なことを示すものと同じ。)


\begin{align}
|\tilde z^{-1}|^{2} &= \frac{|z^{-1}-\alpha |^2}{| 1-\alpha z^{-1} |^2} \\
&= \frac{\left(z^{-1}-\alpha\right)\left(\overline{z^{-1}}-\alpha\right)}{\left( 1-\alpha z^{-1} \right)\left( 1-\alpha \overline{z^{-1}} \right)} \\
&= \frac{ z^{-1} \overline{z^{-1}} + \alpha \left( z^{-1} + \overline{z^{-1}} \right) + \alpha^{2} }{ 1 + \alpha \left( z^{-1} + \overline{z^{-1}} \right) + \alpha^{2} z^{-1} \overline{z^{-1}} } \\
&= \frac{ |z^{-1}|^2 + \alpha \left( z^{-1} + \overline{z^{-1}} \right) + \alpha^{2} }{ 1 + \alpha \left( z^{-1} + \overline{z^{-1}} \right) + \alpha^{2} |z^{-1}|^2 }
\end{align}

式でも確認できるけど、 \tilde z^{-1}=\frac{z^{-1}+\alpha}{1+\alpha z^{-1}} の分母の絶対値と分子の絶対値が等しいことを言えばいいだけので、図を書いたほうがわかりやすい。 図は \alpha \lt 0 なので注意。

f:id:nagiss:20220201171836p:plain

z=\cos \omega + j\sin\omega のとき、 \tilde z の絶対値は1であることがわかったので、新しい変数 \tilde\omega を導入して \tilde z = \cos \tilde\omega + j\sin\tilde\omega と表すことができる。 z=\cos \omega + j\sin\omega のとき \tilde\omega はどんな値になるかも確認したい。(この計算はオールパスフィルタの位相特性の計算と同じ。) さっきと同じように図で考えると、 \omega と \tilde\omega は下図の角度になる。

f:id:nagiss:20220201175536p:plain

図から、高校数学を駆使して \omega と \tilde\omega の関係を計算する。筋のいい方法がなかなか分からなくて大変だった。ネタバレをすると、図の二等辺三角形の各辺の長さを求めて余弦定理を使うと良い。結果は、
\displaystyle \cos\tilde\omega = \frac{-2\alpha + (1+\alpha^{2})\cos\omega}{(1+\alpha^{2})-2\alpha\cos\omega}
となる。-\pi\lt\omega\lt\pi の範囲では、\omega と \tilde\omega の符号が等しいことを使えば、\omega から \tilde\omega が一意に定まる。左辺を \tan\tilde\omega にした場合は、
\displaystyle \tan\tilde\omega = \frac{(1-\alpha^{2})\sin\omega}{(1+\alpha^{2})\cos\omega-2\alpha}
となる。

\omega と \tilde\omega の対応を描画してみるとこんな感じ。

f:id:nagiss:20220201221057p:plain

描画コード

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(4, 4), dpi=200)
for alpha in [0.8, 0.4, 0.0, -0.4, -0.8]:
    omega = np.linspace(0, np.pi, 1000)
    omega_tilde = np.arccos(
        (-2 * alpha + (1 + alpha * alpha) * np.cos(omega)) / ((1 + alpha * alpha) - 2 * alpha * np.cos(omega))
    )
    plt.plot(omega, omega_tilde, label=fr"$\alpha={alpha:.1f}$")
plt.legend()
plt.xlabel(r"$\omega$")
plt.ylabel(r"$\widetilde\omega$")
plt.xlim(0, np.pi)
plt.ylim(0, np.pi)
plt.xticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.yticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.grid()
plt.show()

結局、ここまで長い道のりを辿ってやったのは、\omega で表されたフーリエ変換の式 F(\omega) を、変数変換によって \tilde\omega で表したことである。 \omega と \tilde\omega の関係はさっき計算した通りなので、 \tilde\omega の値から周波数特性を求められるようになった。

オールパスフィルタ

適当なフィルタ H(z) を考える。イメージしやすいようにとりあえず H(z)=1-0.2z^{-1} のローパスフィルタとしておく。このフィルタの周波数特性 F_{H}(\omega) は、
F_{H}(\omega)=H(\cos\omega+j\sin\omega)=1-\frac{0.2}{\cos\omega+j\sin\omega}
と求められる。

ここで、唐突にフィルタ内の遅延要素 z^{-1} をオールパスフィルタ \frac{z^{-1}-\alpha}{1-\alpha z^{-1}} で置き換えたフィルタ G(z) を考えると、
G(z) = 1-0.2\frac{z^{-1}-\alpha}{1-\alpha z^{-1}}
となる。

この式にさっきの変数変換を適用すれば、
G(z)=1-0.2\tilde z^{-1} = H(\tilde z)
と表せる。

H(\tilde z) の周波数特性を考える。置き換える前のフィルタ H(z) の周波数特性は F_{H}(\omega)=1-\frac{0.2}{\cos\omega+j\sin\omega} なので、置き換えたフィルタの周波数特性は同様の計算で F_{G}(\omega)=1-\frac{0.2}{\cos\tilde\omega+j\sin\tilde\omega} となる。

これにより、遅延要素をオールパスフィルタで置き換えたフィルタの特性は、元のフィルタの周波数特性の、周波数の軸を伸縮させたものになっているとわかる。

メル目盛の近似

サンプリング周波数が 16kHz のとき、\alpha=0.42 くらいにするとメル目盛をよく近似する周波数目盛になる……とよく言われるけど 0.42 の根拠をあんまり見たことが無かった気がするので、描画して確認してみる。

f:id:nagiss:20220201223757p:plain

描画コード

plt.figure(figsize=(4, 4), dpi=200)

alpha = 0.42
sample_rate = 16000

omega = np.linspace(0, np.pi, 1000)
omega_tilde = np.arccos(
    (-2 * alpha + (1 + alpha * alpha) * np.cos(omega)) / ((1 + alpha * alpha) - 2 * alpha * np.cos(omega))
)
plt.plot(omega, omega_tilde, label=fr"$\alpha={alpha:.2f}$")

freq = omega * sample_rate / (2.0 * np.pi)
mel = 1000 * np.log2(freq / 1000.0 + 1.0)
max_mel = mel[-1]
target_omega = mel * (np.pi / max_mel)
plt.plot(omega, target_omega, color="black", label="mel")

plt.legend()
plt.xlabel(r"$\omega$")
plt.ylabel(r"$\widetilde\omega$")
plt.xlim(0, np.pi)
plt.ylim(0, np.pi)
plt.xticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.yticks(np.linspace(0, np.pi, 5), [r"$0$", r"$\pi/4$", r"$\pi/2$", r"$3\pi/4$", r"$\pi$"])
plt.grid()
plt.show()

よく近似できてるっぽい。

おわりに

SPTK の freqt がどういう仕組みになっているのかわかりません　だれかおしえてください

[競プロ] 割と真面目に Python から C++ の set を使えるようにしたのでみんな使ってみてほしい

C/C++ 拡張は Python の機能のひとつなので当然 C++ は Python

この記事は何

この間の記事のすごくはやい set を真面目に使いやすくした

真面目ではない使い方説明

これを見ると大体わかる

colab.research.google.com

(20/9/14追記) C 側でクラスも作ってちょっと速くなったけど絶対バグあるだろバージョン

colab.research.google.com

(21/09/05追記) 最新版はこのあたりを見て snippet/wrap_cpp_set.py at master · Lgeu/snippet · GitHub

割と真面目な使い方説明

↑のひとつめのセルの中身をコピペすれば、普通に Python のクラスとして使えるようになる

長過ぎると思ったときは、下の方のセルに圧縮する方法を書いたのでそれを使うと良いかもしれない

CppSet オブジェクトひとつにつき C++ 側では set オブジェクトひとつとイテレータをひとつ保持する
イテレータは基本的に最後に操作を行った位置を指し、これによって一部の操作を高速に行えることがある

g++ 拡張で高速になるメソッドを使わない場合、C++ のコードの 3 行目くらいのコメントアウトを外すと定数倍高速になる

以下各メソッドの説明（計算量の説明で、 n は要素数とする）

CppSet (コンストラクタ)

  • 引数無しで呼ぶと、空の CppSet オブジェクトを作る
  • リストまたはタプルを引数に入れると、その要素が入った CppSet オブジェクトを作る
  • 計算量は \mathcal{O}(n \log n)
    • リストがソート済みの場合は \mathcal{O}(n)
  • イテレータは最小の要素を指す

add

  • 要素を追加して True を返す
  • ただし、既に同じ要素が入っていた場合、要素を追加せず False を返す
  • 計算量は \mathcal{O}(\log n)
  • イテレータは追加した要素（または既に入っていた要素）を指す

remove

  • 指定した要素を削除して、削除した要素の次の要素を返す
    • 最も大きい要素を削除した場合、None を返す
  • 計算量は \mathcal{O}(\log n)
  • イテレータは返した要素を指す
    • None を返した場合、最大の要素の次を指す
  • 指定した要素が入っていなければ KeyError を出す

search_higher_equal

  • 指定した要素と同じかそれより大きい要素を返す
    • そのような要素がなければ None を返す
  • 計算量は \mathcal{O}(\log n)
  • イテレータは返した要素を指す
    • None を返した場合、最大の要素の次を指す

min

  • 最小の要素を返す
  • イテレータは返した要素を指す
  • 計算量は \mathcal{O}(1)
  • 要素数が 0 の場合、 IndexError を出す

max

  • 最大の要素を返す
  • イテレータは返した要素を指す
  • 計算量は \mathcal{O}(1)
  • 要素数が 0 の場合、 IndexError を出す

pop_min

  • 最小の要素を削除し、その値を返す
  • イテレータは削除した要素の次の要素を指す
    • すなわち、削除後の最小の要素を指す
  • 計算量は \mathcal{O}(1)
  • 要素数が 0 の場合、 IndexError を出す

pop_max

  • 最大の要素を削除し、その値を返す
  • イテレータは削除した要素の次の要素を指す
    • すなわち、削除後の最大の要素の次を指す
  • 計算量は \mathcal{O}(1)
  • 要素数が 0 の場合、 IndexError を出す

__len__

  • 要素数を返す
  • 計算量は \mathcal{O}(1)

next

prev

to_list

  • 要素をリストにして返す
  • 計算量は \mathcal{O}(n)

get

erase

  • イテレータの指す要素を削除し、その次の要素を返す
    • そのような要素がなければ None を返す
  • イテレータは返した要素を指す
    • None を返した場合、最大の要素の次を指す
  • 計算量は \mathcal{O}(1)
  • イテレータが最大の要素の次を指していた場合、 KeyError を出す

__getitem__

  • k 番目の要素を返す
    • 負の値もいける
  • 計算量は \mathcal{O}(\log n)
  • イテレータは返した要素を指す
    • g++ 環境でない場合、計算量は \mathcal{O}(n)
  • k 番目の要素が存在しない場合、 IndexError を出す

pop

  • k 番目の要素を削除し、その値を返す
    • 負の値もいける
  • 計算量は \mathcal{O}(\log n)
  • イテレータは返した要素の次の要素を指す
    • g++ 環境でない場合、計算量は \mathcal{O}(n)
  • k 番目の要素が存在しない場合、 IndexError を出す

index

  • 要は bisect_left
  • 計算量は \mathcal{O}(\log n)
    • g++ 環境でない場合、計算量は \mathcal{O}(n)
  • イテレータは変化しない

copy

  • 自身のコピーを返す
  • 計算量は \mathcal{O}(n)
  • 新しいオブジェクトのイテレータは最初の要素を指す

割と真面目な使用例

AtCoder Regular Contest 033 C - データ構造

k 番目の要素

atcoder.jp

圧縮版

atcoder.jp

CPSCO2019 Session1 E - Exclusive OR Queries

イテレータ操作で高速化、g++ 拡張を使わず定数倍高速化

atcoder.jp

割と真面目な注意点

実行する .py があるのと同じ階層に C++ のファイルと setup.py を作るので、同じ名前のファイルを置かないように気をつける

set は C++ だけど乗せてるものとか比較関数とかは Python のものなので速さに限度がある

PyPy では使えない

バグの検証がしきれないのでみんな使って教えてくれるとうれしい

(20/9/23追記) AtCoder 以外のオンラインジャッジで使えるか？

Codeforces

ルールで禁止されているので使ってはいけない

AOJ

利用規約で禁止されているので使ってはいけない

yukicoder

使えるっぽい

ただし、setup部分は↓のように書き換えないといけない

それと、どこかのテストケースにはコンパイル時間が乗るはず

import os
import sys
os.chdir("__pycache__")
with open("set_wrapper.cpp", "w") as f:
    f.write(code_cppset)
with open("setup.py", "w") as f:
    f.write(code_setup)
os.system(f"CC=g++7 {sys.executable} setup.py build_ext --inplace > /dev/null")
os.chdir("../")

from __pycache__.cppset import CppSet

使用例

yukicoder.me

（↑の使用例は pb_ds の tree の比較関数をいじって multiset として使うテクを使っている（具体的には 13 行目の最後を Py_LT から Py_LE に変えている）ので注意（実際のところ完全に multiset にはならなくて、lower_bound とかの以上・超過があべこべになったりする））

まとめ

C++ は Python なので当然 Python 何もわからない

Numba のコンパイルが通らなかった時の対処

Numba はいいぞ

この記事は何

ふつうの Python なら動くけど Numba では動かないようなコードを列挙して、対処法を書いたもの
主に AtCoder 目的だけどそれ以外でも役に立つはず

Numba のバージョン 0.48.0 くらいの情報なので将来的にいろいろ変わってくると思うので注意（2020 年 8 月現在で AtCoder に入ってるのも 0.48.0）

先に読んでおくといいかもしれない記事

qiita.com

ikatakos.com

Numba で使えないもの

2 次元以上の ndarray のイテレーション

できない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    array = np.random.rand(5, 2)  # 5x2 array
    for a in array:  # コンパイルエラー
        ...
エラーメッセージ
Direct iteration is not supported for arrays with dimension > 1. Try using indexing instead.
[1] During: typing of intrinsic-call at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)

File "untitled.py", line 7:
def solve():
    <source elided>
    array = np.random.rand(5, 2)  # 5x2 array
    for a in array:  # コンパイルエラー
    ^
対処

range で回す

@numba.njit("void()", cache=True)
def solve():
    array = np.random.rand(5, 2)  # 5x2 array
    for i in range(len(array)):
        a = array[i]
        ...

変な(?)方法での空のリスト作成

Numba が型を推測できないとエラーが出る
list 以外に dict と set でも起こる
結構ありがちでエラーメッセージもわかりにくかったりするので、とりあえず型を明示しておくのがいいかもしれない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    lst = [[] for _ in range(10)]  # コンパイルエラー
    lst[0].append(0)
エラーメッセージ
Undecided type $26load_method.11 := <undecided>
[1] During: resolving caller type: $26load_method.11
[2] During: typing of call at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)


File "untitled.py", line 7:
def solve():
    <source elided>
    lst = [[] for _ in range(10)]  # コンパイルエラー
    lst[0].append(0)
    ^
対処

どうにかして型を Numba に教える

@numba.njit("void()", cache=True)
def solve():
    lst = [[0] * 0 for _ in range(10)]  # 型を明示する
    lst[0].append(0)

dict の場合はいらない要素を入れておくとか
set の場合は {0}-{0} とかすれば Numba くんはわかってくれる

辞書の内包表記

dict と set は内包表記で生成できない
あと dict はリストからの生成とかもできない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    fib = [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]
    inv_fib = {v: i for i, v in enumerate(fib)}  # コンパイルエラー
エラーメッセージ
Use of unsupported opcode (MAP_ADD) found

File "untitled.py", line 7:
def solve():
    <source elided>
    fib = [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]
    inv_fib = {v: i for i, v in enumerate(fib)}  # コンパイルエラー
    ^
対処

ひとつずつ入れる

@numba.njit("void()", cache=True)
def solve():
    fib = [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]
    inv_fib = {}
    for i, v in enumerate(fib):
        inv_fib[v] = i

リストを値に取る辞書

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    dictionary = {3023: [0, 1, 2], 4006: [3, 4, 5]}  # コンパイルエラー
    d = dictionary[3023]
エラーメッセージ
list(int64) as value is forbidden
[1] During: typing of dict at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)

File "untitled.py", line 7:
def solve():
    dictionary = {3023: [0, 1, 2], 4006: [3, 4, 5]}  # コンパイルエラー
    ^
対処

numpy.ndarray でもいいならそれを使う
そうでないなら、値は別にリストか何かで持っておいて、辞書にはそのインデックスを入れる

@numba.njit("void()", cache=True)
def solve():
    dictionary = {3023: 0, 4006: 1}
    container = [[0, 1, 2], [3, 4, 5]]
    d = container[dictionary[3023]]

pow の第 3 引数

使えない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    mod = 10 ** 9 + 7
    inv10 = pow(10, mod-2, mod)  # コンパイルエラー
エラーメッセージ
Invalid use of Function(<built-in function pow>) with argument(s) of type(s): (Literal[int](10), int64, Literal[int](1000000007))
Known signatures:
 * (int64, int64) -> int64
 * (int64, uint64) -> int64
 * (uint64, int64) -> int64
 * (uint64, uint64) -> uint64
 * (float32, int32) -> float32
 * (float32, int64) -> float32
 * (float32, uint64) -> float32
 * (float64, int32) -> float64
 * (float64, int64) -> float64
 * (float64, uint64) -> float64
 * (float32, float32) -> float32
 * (float64, float64) -> float64
 * (complex64, complex64) -> complex64
 * (complex128, complex128) -> complex128
In definition 0:
    All templates rejected with literals.
In definition 1:
    All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: resolving callee type: Function(<built-in function pow>)
[2] During: typing of call at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)


File "untitled.py", line 7:
def solve():
    <source elided>
    mod = 10 ** 9 + 7
    inv10 = pow(10, mod-2, mod)  # コンパイルエラー
    ^
対処

代わりのものを作っておく

@numba.njit("i8(i8,i8,i8)", cache=True)
def pow_mod(base, exp, mod):
    exp %= mod - 1
    res = 1
    while exp:
        if exp & 1:
            res = res * base % mod
        base = base * base % mod
        exp >>= 1
    return res

@numba.njit("void()", cache=True)
def solve():
    mod = 10 ** 9 + 7
    inv10 = pow_mod(10, mod-2, mod)

built-in の sum 関数

max は使えるのに sum は何故か使えない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    a = np.random.rand(5)
    s = sum(a)  # コンパイルエラー
エラーメッセージ
Untyped global name 'sum': cannot determine Numba type of <class 'builtin_function_or_method'>

File "untitled.py", line 7:
def solve():
    <source elided>
    a = np.random.rand(5)
    s = sum(a)  # コンパイルエラー
    ^
対処

numpy.sum か numpy.ndarray.sum を使う
リストの場合は numpy.sum でもエラーになるのでそのときは numpy.ndarray に変換するとかひとつずつ足すとかする

@numba.njit("void()", cache=True)
def solve():
    a = np.random.rand(5)
    s = np.sum(a)  # s = a.sum() でもいい

numpy.max とか numpy.ndarray.max とかの axis

numpy.sum は axis が使えるのに numpy.max は axis が使えない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    array = np.random.rand(4, 5)
    m = array.max(1)
エラーメッセージ
[1] During: resolving callee type: BoundFunction(array.max for array(float64, 2d, C))
[2] During: typing of call at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)

Enable logging at debug level for details.

File "untitled.py", line 7:
def solve():
    <source elided>
    array = np.random.rand(4, 5)
    m = array.max(1)
    ^
対処

for を回す

@numba.njit("void()", cache=True)
def solve():
    array = np.random.rand(4, 5)
    m = np.empty(4, dtype=array.dtype)
    for i in range(4):
        m[i] = array[i].max()

2 次元以上の ndarray の boolean indexing

できない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    array = np.random.rand(4, 5)
    array[array < 0.5] = 0  # コンパイルエラー
エラーメッセージ
Invalid use of Function(<built-in function setitem>) with argument(s) of type(s): (array(float64, 2d, C), array(bool, 2d, C), Literal[int](0))
 * parameterized
In definition 0:
    All templates rejected with literals.
In definition 1:
    All templates rejected without literals.
In definition 2:
    All templates rejected with literals.
In definition 3:
    All templates rejected without literals.
In definition 4:
    All templates rejected with literals.
In definition 5:
    All templates rejected without literals.
In definition 6:
    All templates rejected with literals.
In definition 7:
    All templates rejected without literals.
In definition 8:
    TypeError: unsupported array index type array(bool, 2d, C) in [array(bool, 2d, C)]
    raised from C:\Users\nagiss\Anaconda3\lib\site-packages\numba\typing\arraydecl.py:71
In definition 9:
    TypeError: unsupported array index type array(bool, 2d, C) in [array(bool, 2d, C)]
    raised from C:\Users\nagiss\Anaconda3\lib\site-packages\numba\typing\arraydecl.py:71
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of setitem at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)

File "untitled.py", line 7:
def solve():
    <source elided>
    array = np.random.rand(4, 5)
    array[array < 0.5] = 0
    ^
対処

numpy.where を使う

@numba.njit("void()", cache=True)
def solve():
    array = np.random.rand(4, 5)
    array = np.where(array < 0.5, 0, array)

ndarray の None による次元の追加

できない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    a = np.random.rand(4, 5)
    a = a[:, None, :]  # コンパイルエラー
    assert a.shape == (4, 1, 5)
エラーメッセージ
Invalid use of Function(<built-in function getitem>) with argument(s) of type(s): (array(float64, 2d, C), Tuple(slice<a:b>, none, slice<a:b>))
 * parameterized
In definition 0:
    All templates rejected with literals.
In definition 1:
    All templates rejected without literals.
In definition 2:
    All templates rejected with literals.
In definition 3:
    All templates rejected without literals.
In definition 4:
    All templates rejected with literals.
In definition 5:
    All templates rejected without literals.
In definition 6:
    All templates rejected with literals.
In definition 7:
    All templates rejected without literals.
In definition 8:
    All templates rejected with literals.
In definition 9:
    All templates rejected without literals.
In definition 10:
    All templates rejected with literals.
In definition 11:
    All templates rejected without literals.
In definition 12:
    TypeError: unsupported array index type none in Tuple(slice<a:b>, none, slice<a:b>)
    raised from C:\Users\nagiss\Anaconda3\lib\site-packages\numba\typing\arraydecl.py:71
In definition 13:
    TypeError: unsupported array index type none in Tuple(slice<a:b>, none, slice<a:b>)
    raised from C:\Users\nagiss\Anaconda3\lib\site-packages\numba\typing\arraydecl.py:71
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)
[2] During: typing of static-get-item at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (7)

File "untitled.py", line 7:
def solve():
    <source elided>
    a = np.random.rand(4, 5)
    a = a[:, None, :]  # コンパイルエラー
    ^
対処

reshape か expand_dims を使う
ただし expand_dims の第 2 引数に tuple は使えない

@numba.njit("void()", cache=True)
def solve():
    a = np.random.rand(4, 5)
    a = np.expand_dims(a, 1)
    assert a.shape == (4, 1, 5)

int.bit_length

使えない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    b = (998244353).bit_length()
エラーメッセージ
Unknown attribute 'bit_length' of type Literal[int](998244353)

File "untitled.py", line 6:
def solve():
    b = (998244353).bit_length()
    ^

[1] During: typing of get attribute at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (6)

File "untitled.py", line 6:
def solve():
    b = (998244353).bit_length()
    ^
対処

np.log2 とかでなんとかする（249-1 以上は誤差に注意、あと 0 も正しく動かない）

@numba.njit("void()", cache=True)
def solve():
    b = int(np.log2(998244353)) + 1
    print(b)

collections

ほぼ使えない

対処

defaultdict -> 値があるか自力で確認する
deque -> リングバッファみたいなのを適当に実装する
Counter -> 自力で数える

itertools

使えない

対処

itertools を使う部分はコンパイルしないように切り分けるか、あらかじめ代わりになりそうなものを用意しておく？

@numba.jit("i8[:,:](i8[:],i8)", cache=True)
def combinaions(arr, r):
    n = len(arr)
    assert 0 <= r <= n
    res_length = 1
    for i in range(r):
        res_length = res_length * (n-i) // (1+i)
    res = np.empty((res_length, r), dtype=arr.dtype)
    idxs_arr = np.arange(r)
    for idx_res in range(res_length):
        res[idx_res] = arr[idxs_arr]
        i = 1
        while idxs_arr[r-i] == n-i:
            i += 1
        idxs_arr[r-i] += 1
        for j in range(r-i+1, r):
            idxs_arr[j] = idxs_arr[j-1] + 1
    return res

↑の実装はジェネレータじゃないので少し探索してやめるような場合には効率が悪くなってしまう
これを嫌うなら C++ の next_permutation みたいなのを用意しておくと汎用性も高くて良さそう

string が返る関数

str とか bin とか format とか "%d" % 42 とかは使えない

エラーになるコード
@numba.njit("void()", cache=True)
def solve():
    popcnt = bin(4047).count("1")  # コンパイルエラー
エラーメッセージ
Untyped global name 'bin': cannot determine Numba type of <class 'builtin_function_or_method'>

File "untitled.py", line 6:
def solve():
    popcnt = bin(4047).count("1")  # コンパイルエラー
    ^
対処

Numba で文字列を扱おうとしない

popcount についてはあらかじめ用意しておく（参考: Python 3でpopcountを計算する - にせねこメモ）

@numba.njit("u8(u8)", cache=True)
def popcount(n):
    n = (n & 0x5555555555555555) + (n>>1 & 0x5555555555555555)
    n = (n & 0x3333333333333333) + (n>>2 & 0x3333333333333333)
    n = (n & 0x0f0f0f0f0f0f0f0f) + (n>>4 & 0x0f0f0f0f0f0f0f0f)
    n = (n & 0x00ff00ff00ff00ff) + (n>>8 & 0x00ff00ff00ff00ff)
    n = (n & 0x0000ffff0000ffff) + (n>>16 & 0x0000ffff0000ffff)
    n = (n & 0x00000000ffffffff) + (n>>32 & 0x00000000ffffffff)
    return n

@numba.njit("void()", cache=True)
def solve():
    popcnt = popcount(4047)

関数外の変数の書き換え

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

@numba.njit("void()", cache=True)
def solve():
    array[0] = 4  # コンパイルエラー
エラーメッセージ
Invalid use of Function(<built-in function setitem>) with argument(s) of type(s): (readonly array(int32, 1d, C), Literal[int](0), Literal[int](4))
 * parameterized
In definition 0:
    All templates rejected with literals.
In definition 1:
    All templates rejected without literals.
In definition 2:
    All templates rejected with literals.
In definition 3:
    All templates rejected without literals.
In definition 4:
    All templates rejected with literals.
In definition 5:
    All templates rejected without literals.
In definition 6:
    All templates rejected with literals.
In definition 7:
    All templates rejected without literals.
In definition 8:
    TypeError: Cannot modify value of type readonly array(int32, 1d, C)
    raised from C:\Users\nagiss\Anaconda3\lib\site-packages\numba\typing\arraydecl.py:179
In definition 9:
    TypeError: Cannot modify value of type readonly array(int32, 1d, C)
    raised from C:\Users\nagiss\Anaconda3\lib\site-packages\numba\typing\arraydecl.py:179
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of staticsetitem at C:/Users/nagiss/PycharmProjects/untitled/untitled.py (8)

File "untitled.py", line 8:
def solve():
    array[0] = 4  # コンパイルエラー
    ^
対処

引数で渡す

@numba.njit("void(i4[:])", cache=True)
def solve(array):
    array[0] = 4

標準入力

できないので諦める

関数を返す関数

できないので諦める

関数内の関数の再帰

できないので諦める

他色々

できないと思って諦める

まとめ

Numba はいいぞ

[Python]平衡二分木が必要な時に代わりに何とかする変なテク[競プロ]

はじめに

読んで

qiita.com

変なテクが足りない

と思ったので勝手に付け足す

解決法5: 変なデータ構造を作る

平衡二分木みたいな複雑なデータ構造は動作が重いが、Python の標準ライブラリの関数でほとんど処理が終わるようなデータ構造なら十分実用できる速さになる

実装が簡単なので機能を足すのも簡単

例

新ジャッジでの実行時間なので比較するときは注意

CPSCO2019 Session1 E - Exclusive OR Queries (1344 ms)

自分的ベンチマーク問題

atcoder.jp

ARC033 C - データ構造 (500 ms)

k 番目の値を取り出す

atcoder.jp

すぬけのお誕生日コンテスト F - Lake (685 ms)

tuple を乗せる

atcoder.jp

解決法6: C++ の set に Python のオブジェクトを乗せる

コンパイル時間にコンパイルする　何もおかしいことはしていない

爆速だけど機能を足すのは面倒、あと PyPy では使えない

例

CPSCO2019 Session1 E - Exclusive OR Queries (916 ms)

int と float を同時に入れるとバグるので注意

atcoder.jp

ARC033 C - データ構造 (317 ms)

k 番目の要素は g++ 拡張で

atcoder.jp

すぬけのお誕生日コンテスト F - Lake (630 ms)

tuple も乗る

atcoder.jp

20/08/20 追記: これ前に誰かがやってた気がしていて記事書いたときは見つけられなかったんですが、今先駆者を見つけました（ごめんなさい）
自分のやり方と結構違ってるっぽい？

atcoder.jp

(20/07/14追記) 解決法7: SQLite を使う

そういえばデータベースは B 木だか B+ 木だかで実装されているので利用できる

例

CPSCO2019 Session1 E - Exclusive OR Queries (1785 ms)

上手く書けばもっと早くなるかも

atcoder.jp

k 番目の要素は取り出せなさそう

おわりに

おわりです

Kaggleとかのユーティリティのメモ

毎回探すのアホだろと思ったので自分用まとめ
適宜更新する

  • タイマー、seed固定、混同行列、feature importance 等
    qiita.com

  • ヒートマップ高速描画
    spcx8.hatenablog.com

  • trainとtestのdataframeを結合したり切り離したりするやつ

def merge_train_test(df_train, df_test):
    if "target" not in df_test.columns.values:
        df_test["target"] = -1  # df_train["target"] に存在しない値にする
    res = pd.concat([df_train, df_test])
    res.reset_index(inplace=True, drop=True)
    return res

def split_train_test(df):
    df_train = df[df["target"] != -1]
    df_test = df[df["target"] == -1]
    df_train.reset_index(inplace=True, drop=True)
    df_test.reset_index(inplace=True, drop=True)
    return df_train, df_test
  • count encoding
class CountEncoder:
    def fit(self, series):
        self.counts = series.groupby(series).count()
    
    def transform(self, series):
        return series.map(self.counts).fillna(0)
  • pytorchで異なるシーケンス長のデータをひとつのミニバッチにまとめる
from torch.nn.utils.rnn import pad_sequence, pack_padded_sequence, pad_packed_sequence

def pad_collate(batch):
    feature, target = zip(*batch)
    lengths = torch.tensor([len(f) for f in feature])
    feature = pad_sequence(feature, batch_first=True)
    target = pad_sequence(target, batch_first=True)
    return feature, target, lengths

dataloader = torch.utils.data.DataLoader(dataset, batch_size=BATCH_SIZE, shuffle=True, collate_fn=pad_collate)

RNNとかに通すときは以下のようにpackする

lstm = nn.LSTM(n_features, 32, batch_first=True, bidirectional=True)
for feature, target, lengths in dataloader:
    h0, c0 = torch.randn(2, BATCH_SIZE, 32, device=device), torch.randn(2, BATCH_SIZE, 32, device=device)

    feature = pack_padded_sequence(feature, lengths, enforce_sorted=False, batch_first=True)
    x, _ = lstm(feature, (h0, c0))
    x, l = pad_packed_sequence(x, batch_first=True)
    assert (l == lengths).all()

    break

巡回セールスマン問題(TSP)のメモ（未完）

1ヶ月くらい前に途中まで書いたけど飽きて放置されていたものを適当に処理して投下する

前置き

対称巡回セールスマン問題(STSP)のみ扱う
（実は対称性が無い(ATSP)のほうが枝刈りしやすく簡単らしい）

巡回セールスマン問題: グラフの最も短いハミルトン閉路を見つける問題
ハミルトン閉路: グラフのすべての頂点を一度だけ通る閉路

厳密解関連

Held-Karpアルゴリズム (1962)

いわゆるbitDP
頂点数nに対して \mathcal{O}(n ^ 2 2 ^ n)で厳密解を求められる

Held-Karp下界 (1970)

The Traveling Salesman Problem and Minimum Spanning Trees (Held and Karp, 1970)

一般に、枝刈り探索は良い下界（最適解はこの値より小さくならないよという値）があると高速になる

まず、グラフの最小1-木の大きさはTSPの解の下界になっている

1-木とは、頂点1, 2, 3, ..., nで構成されるグラフについて、頂点1を除いた2, 3, ..., nのグラフの全域木と、頂点1からそれ以外の頂点への2本の辺で構成されるグラフのことである（1-木という名前だが木ではない）

ここで、グラフのすべてのハミルトン閉路は1-木であるから、最小1-木は最小のハミルトン閉路よりも小さい（巡回セールスマン問題の下界となる）

最小1-木が運良くハミルトン閉路であった場合、その1-木はTSPの解に等しい

グラフの最小1-木を求めるには、まず頂点1を除いたグラフの最小全域木をクラスカル法などで求め、そこに頂点1に繋がる辺のうち最も短い2つを足せば良く、辺の数 m に対して \mathcal{O}(m \log m) で求められる

Held-Karp下界では、グラフの辺の重みをうまく変えて最小1-木を求めることによってより良い下界を求める

具体的には、グラフの各頂点についてペナルティ \pi _ 1 , \pi _ 2 , ..., \pi _ n を定め、
頂点 i, j を結ぶ辺の重み c _ {ij} を、すべて c _ {ij} + \pi _ i + \pi _ j に変えたグラフで最小1-木を求める

この最小1-木の大きさから 2 \sum _ {i=1} ^ n \pi _ i を引いた値は、TSPの下界となっている

なぜなら、上記の方法で辺の重みを変えると任意のハミルトン閉路の大きさは 2 \sum _ {i=1} ^ n \pi _ i だけ等しく変わり、
最短路自体は変化しないためである

式で書くと以下の通り

\ (重み修正グラフのTSPの解) \geq (重み修正グラフの最小1 \unicode{x2d} 木の大きさ) \\ \displaystyle \Leftrightarrow (元のグラフのTSPの解) + 2 \sum _ {i=1} ^ n \pi _ i \geq (重み修正グラフの最小1 \unicode{x2d} 木の大きさ) \\ \displaystyle \Leftrightarrow (元のグラフのTSPの解) \geq (重み修正グラフの最小1 \unicode{x2d} 木の大きさ) - 2 \sum _ {i=1} ^ n \pi _ i \\

良い下界（厳密解に近い値をとる下界）が得られるように \pi _ i を定めるには、まずペナルティ無しで最小1-木を作り、次数が2より大きい頂点についてはペナルティを増やし、次数が2より小さい頂点についてはペナルティを減らして、最小1-木を再計算する…といったことを何度か繰り返す

参考: 巡回セールスマン問題 - 日本オペレーションズ・リサーチ学会 (1979)

整数計画問題を使った解法

特殊な形の整数計画問題に帰着して解く　つよい　でもよくわからない

組合せ最適化の黄色い本に載ってた気がする

85,900都市の問題を解いたTSPソルバのConcordeもこれを元にしている

Certification of an Optimal TSP Tour Through 85,900 Cities (Applegate et al., 2009)

近似解

構築法と改善法があるけど改善法のが重要なはず

2-opt, 3-opt, Or-opt

いつもの

Lin-Kernighan (LK) (1973)

k-optの一般化と呼ばれるやつ

現時点での使う辺と使わない辺を交互に辿る閉路を作って、使うか使わないかを反転させる

効率を高めるヒューリスティックな工夫が色々ある

Iterated Lin-Kernighan (ILK) (1991)

Lin-Kernighanで実現できない遷移(double bridge move)によって局所最適から脱する

Chained Lin-Kernighan (CLK) (2003)

ILKの進化したやつらしい

concordeに実装があるけどよくわからない

github.com

Lin-Kernighan-Helsgaun (LKH) (2000)

Lin-Kernighanを効率的にやる　つよい

論文が教科書みたいに丁寧
引用数がすごい

An Effective Implementation of the Lin-Kernighan Traveling Salesman Heuristic

LKH (Keld Helsgaun)

巡回セールスマンみたいな問題のなるべく良い近似解が欲しくなったらとりあえずこれを調べると良さそう

最近の

EAXとかACOとかめっちゃ強いのがあるらしい

まとめ

マラソンマッチはたいへん

ICPC 2019 アジア地区横浜大会　参加記

チームPICOPICOPONで参加してきた

とまと視点

25toma.hatenablog.com

PICOPICOPON

メンバー

yamad

AOJ-ICPCを17万点埋めてる実装のプロ
見てきた問題数がすごいので、典型っぽい問題を相談すると大体なんとかしてくれる
去年はWAsedACとしてWFに出場していた
DPが得意
早起きのプロ

toma25

とまと
去年はMohuhu Universityの一員だった
実装が得意
DPが得意
本番前日に青から黄色になるエンターテイナー
チーム練ではよくデバッグとか高速化にまわってた

nagiss

ICPC初参加
普段はPythonを使ってるけどC++を書くと他の人がデバッグしてくれることに気付いたのでICPCではC++を書くことにした
実装力がない
DPができない

PICOPICOPON is one of Jiro Ramen.

tabelog.com

おいしい

1日目

昼食

早稲田の人たちで中華街に行って中華を食べる

おいしい

会場入り

Eigo is murimuri of muriなのでyamad-sanに全てを任せる

practice contest

環境構築担当のyamadは色々と確認に忙しそうだったが、後ろに座っていた自分は無だった

夜

ABCに出ると、EとFがDPなのでつらい気持ちになる
とまとがperf 2400を叩き出してチーム全員黄色になった　めでたい

こどふぉdiv2に出ると、EがDPなのでつらい気持ちになる
やめてくれ？

2日目

朝

なんか日が差してきたので自然に目覚めた

昨日数秒間に合わなかったこどふぉのEを投げたら普通にWAが出た　DPはつらい

コンテスト

問題文

Aizu Online Judge Arena

0:00 [0/11]

作戦通り、yamadが環境構築、とまとがAを読み、自分がBを読む

環境構築が終わったあたりでAもBも読み終わっており、どうやらAが少しややこしいらしくBがやるだけだったので、yamadにBの概要を伝えて実装してもらう

Q. なんで自分で実装しないんですか？
A. 3重ループなんて自分で書いたらバグりそうなので

Bの概要を伝えたあと、とまとにA問題の様子を聞いて実装しやすい方法を一緒に考える

0:21 [1/11]

B問題が通ったので、A問題をとまとが実装する

自分は横から覗き込んでバグを見つける役割をした

0:26 [2/11]

A問題が無事通る

次に読む問題を決めようと持ち込んだサイコロを振ると5が出たので、Cから数えて5番目のF（？？？？？）を読む

オーとまとンっぽいから後でとまとに聞いてみるか…と思いながら読み終わると、とまとは前から、yamadは後ろから読んでいるようだったので、次にEを読むことにする

0:45

Eの問題文が長くて苦労しているとHを解くチームが現れ始めている

yamadがHを考えつつ実装を始めているようだった

0:55

Eを解くチームが現れ始める

EはDPっぽい雰囲気なので「多分DPだと思うんですけど～」と言いながらyamadさんに問題を共有すると、Hより早そうらしいのでHを中断してEの実装を始めるようだった

代わりにHの問題を伝えられて、だいたい解法わかってるけどうまく実装する方法考えてみたいなことを言われた

yamadはEの実装に忙しそうだったのでだいたいわかってる解法を自分で考え直すしか無く、この時間は意味があるのか…？と思いながらHを考察するのはしんどかった

1:10

yamadがEに思い違いがったらしく、Hの調子を聞きに来たのでその時点の考察を答えると納得したようでHの実装を始めたようだった

引き続いてHを詰めていく

1:38 [3/11]

Hが通ってようやく3完

とまとが他の問題を色々読んでいたのでCとDあたりを共有してもらう

Gが数チームに解かれていたが文字列なのでyamadの担当らしい

（このあたりの時間帯で自分が何をしてたかあまり覚えていない）

PCが空く時間帯が続く　もったいないけど仕方ない

しばらくするとyamadがGの実装を始めた

2:20

E、DP得意な2人がなかなか詰まってるしDPじゃないんじゃないの～？と思ってグラフで考えるとなんかうまくいきそうだったので2人に考察を共有した

yamadがその方針でEの実装を始めた

2:31 [4/11]

Eが通った　早い

調子が良いのでお菓子を取りに行くとパンとおにぎりがあったのでパンをもらう
ホイップクリームパンのホイップクリームの量が多くて良かった

Gの実装をしているyamadに何か取ってきますか？と聞いたがいらないらしい
お茶くみ係への道は長い

自分たちが通しておらず上位の他のチームが通している問題はあとCGIなので自分はCとIを考えることにする

とまとはIは入力が木ならHL分解でできると主張していて、自分はHL分解を知らないのではえ～となっていた
少し考えると入力にサイクルがあればそれは逐次的に潰せるので実質入力は木だろと思ったのでそれを伝えようとするが伝達ミスをしてしまい時間をロスする

3:10

yamadがGを実装したがWAが出る

とまとがデバッグに回る

3:28 [5/11]

Gが通った

そのまま考察はできているが実装激重のIをyamadに伝える
後半部分トポロジカルソートでいけませんか？と口を挟んだけど即座に反例を出されたのでごめんなさいになった

yamadが実装をしてとまとと自分はCを考える

おいおいCこれ座標圧縮してDPでは？とわかってくる　今日もDPかよ勘弁してくれ～～

Iの実装が煮詰まっていくと、とまとがIのテストケース作成に引き抜かれる　きつい

4:33 [6/11]

Iが一発で通る（天才）

Cの考察が大体できていたので急いでyamadに伝える

考察はできていたけど実装レベルまで落とせていなかったので、yamadがセグ木を写経したりする横で自分はセグ木に乗せる値の整理をしたりする

4:57

とても間に合いません　諦めです

yamad「指がつりそう」

結果

https://icpcsec.firebaseapp.com/standings/

f:id:nagiss:20191121200610p:plain

6完797分で11位だった

Cは解くべき問題だった気がしていたので惜しかったなあと思っていたけど、蓋を開けてみると思ったより良い位置にいた

因縁

模擬国内 f:id:nagiss:20191121202324p:plain

夏合宿Day1 f:id:nagiss:20191121202557p:plain

夏合宿Day2 f:id:nagiss:20191121202915p:plain

ああ、やっぱり今回もだめだったよ

懇親会

11位ということでいい生活さんとレトリバさんから企業賞をいただけた　うれしい

おいしい

懇親会では企業ブースが並んでいて、DeNAのブースではパネルにシールを貼るタイプのアンケートをやっていた
自分はアンケートに参加しなかったが、「どんなコンテストに参加していますか？」という質問のKaggleの項目にシールが2枚貼ってあるのを見かけた　誰だろう

その後はGoogleのブースで配布されていた問題に全ての時間を溶かした　は？

3日目

企業見学の日

NEC

グループC（東京近郊組らしい）はNECでご飯だったが、これがめちゃくちゃおいしかった
Twitterを見ても任意の人がNECに胃袋を掴まれていたようだった

説明をしていただいた方に機械学習をやっていて音楽を作っている方がいて親近感が湧いたので昼食のときにお話を聞きに行った
NECの中でもKaggleをやっている人たちがいて集まってチームを組んでメダルを獲ったりしているらしい
競プロの集まりはその方の部署にはまだ無いらしかった

しながわ水族館

ビズリーチの見学ができなくなって代わりに行くことになったらしい

たのしい

Huawei

まとめ

NECのご飯(来客用)はおいしい……覚えましたし