僕の逆引きnumpy・scipyメモ

# 逆引きnumpyメモ
numpy? ムズいよね。使い慣れるまで好きになれないと思うよ。だけど、俺は負けないよ。
え〜、numだっ、numpyが躍動する俺のpythonを、皆さんに見せたいね。
**※ この記事は少しずつ加筆修正すると思います。**

## 俺的マイナーイディオム
下記は解説記事とかであまり見ない書き方だけど、
僕の中でほとんどイディオムなやつなのでこの記事では使うね。

あと、xとyは全部np.ndarrayの行列とします。

```python
# グラム行列ベースに何かする
# この書き方なら、例えばカーネル関数を明確に分けて定義できる。
# *は関数にぶち込む時に展開して貰うための修飾
fuga = hoge(*np.meshgrid(x, y))

# Aという行列のx軸に対してhogeを行ったものを取り出す
# map的な何か。何故か皆for文を使うけど、numpyは関数型っぽく書いたほうがわかり易くない？
fuga = np.apply_along_axis(hoge, x, A)

# イディオムというか、普通にdot積なんだけど、
# 何故皆これ使わずにnp.dotばかりなんですか？サクッと書く時は楽でいいんだけど…
hoge @ fuga
```

# 初期化系
## Iを作る
この辺は流石に要らないかと思ったけど、意外とググってる自分が居るんだな。

```
np.eye(5)
```

## ただの等差数列をつくる
0から100までの大きさの、互いの差が1の等差数列。

```python
np.arange(0, 100, 1)
```

といっても、最後の1とかは黙ってれば1にしてくれるので

```python
np.arange(100)
```
とか。ま、基本やね。

## 等比数列を作る
これは、1から10の3乗までを1000分割するよという意味。
適当に読み替えてくださいな。

```
np.logspace(1, 3, 1000, base=10.0)
```


# 基本的なところ
## 一番大きい数を取ってくる
np.nanmaxが良い。
この関数はNanを無視してくれる。
まぁ、Nan引っ掛けて欲しい場合もあるやろうけどね。

## 逆順にする

np.flipでおｋ

## 連結
並列というか、縦につなげるやつはnp.vstack
直列というか、横につなげるやつnp.hstack
なんだけど、この関数はnp.vstack(hoge, fuga)ではなくて
np.vstack((hoge, fuga))みたいに使うのです。
あとは、np.concatenateとかもあります。

## グラム行列を作る
普通にグラム行列を作る。こんな時はmeshgridがいい。
…なんだけど、meshgrid自体はnumpy芸とも言いたくなるような
技巧たっぷりな書き方になってしまうのは僕だけ？

```python
x = np.arange(10)
y = np.arange(5)
xg, yg = np.meshgrid(x, y)
gram = xg * yg  # ここがやりたい演算
```

だけど、こう書くことも出来る。この場合、下記のカーネル法によく似ているという意味で
わかりやすいので、この方が僕は好き。(カーネル法は一般化内積ですしね)

```python
from operator import mul
x = np.arange(10)
y = np.arange(5)
gram = mul(*np.meshgrid(x, y))  # ここではまだ大差ない
```
# 回帰とか
## カーネル法
取り敢えずガウシアンカーネルで書いてみた。
上記グラム行列のmulを適当にカーネルで置き換えてあげればいいだけ。
ちなみに、この書き方なら色々なカーネル関数…
例えば、ガウシアンや線形等を簡単に入れ替えられるので便利。

```python
x = np.arange(10)
y = np.arange(5)

def kernel(a: np.ndarray, b: np.ndarray, v: float) -> np.ndarray:
    return np.exp(- v * np.square(a - b))

gram = kernel(*np.meshgrid(x, y), 3)
```
二次正則化付きの場合は

```python
lmd = 0.1  # 適当
gram = kernel(*np.meshgrid(x, x)) + lmd * np.eye(x.shape[0])
```

というのを応用して、以前、ふざけて書いたカーネル回帰ゴルフ。

```python
def solve_kernel(x: np.ndarray, y: np.ndarray,
                 kernel: Callable, lmd: float) -> np.ndarray:
    K = kernel(*np.meshgrid(x, x)) + lmd * np.eye(x.shape[0])
    alpha = np.linalg.solve(K, y)
    return np.vectorize(lambda xs: np.sum(alpha * kernel(x, xs)))
```
カーネル回帰はPEP8守っても実質たった4行で書ける。numpyはえらいなぁ…。

## 三次元で山っぽいのを作ってみる

100行100列で、ガウシアンな山を作ってみる。
これにもmeshgridが使える。若干黒魔術っぽいのが悲しいが、大したことはない。

```python
def mount(x: int, y: int) -> np.ndarray:
    return np.exp(-(np.square(x) + np.square(y)) / 1000)

x, y = 100, 100
g = mount(*np.meshgrid(*(np.arange(-n/2, n/2, 1) for n in [x, y])))
```

何故こんなの作ったかって？クラスタリングの実験用に書いたやつです。
meshgridと*の組み合わせは黄金かも知れないし、コーディング規約違反と言われて殴られるのかも知れない。

# 関数型的な何か
## for文が嫌な時
複数やり方がある。まぁ、mapとかは今更なのでいいでしょう。
ちなみに、mapに限らず、内包表記もpythonなので糞遅いことを自覚すること。
そんなものに頼っているようではnumpy使いとして二流である。
(ごく稀に、そのほうが速いこともありはする)

### ユニバーサル関数を使う
ユニバーサル関数という、numpyっぽい挙動の関数を作れるので、それを使う方法が一つ。
ただし、速度が速くなるわけではない。

```python
A = [[2, 4, 8], [1, 1, 1]]
def test(x):
    return x * 2

np.vectorize(test)(A)
> array([[4, 8, 16], [2, 2, 2]])
```
ちなみに、デコレータとしても使えるけど、コーディング規約的な意味で許されるかは分からん。
ちゃんと関数になるためには引数がただの数じゃないと駄目とか、色々あるっぽい。

### np.apply_along_axisを使う
この方法は軸を指定できるのでかなり応用が効く。
しかも、公式曰く速いらしい。vectorizeとどっちのほうが速いかは知らない。

```python
# 上記と同様にして
np.apply_along_axis(test, 0, A)
```

ちゃんと動くためには次元の制約があったりもするらしい。
他にも関数はいくつかある。https://numpy.org/devdocs/reference/routines.functional.html

### np.tileの応用
例えば、あるarrayの列一つ一つに、違う事をしていきたいときがあるとする。
こんな感じで。

```python
matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)

result = []
for r, num in zip(matrix, nums):
    result.append(r * num)
```

このコードは愚の骨頂である。非常に遅いからだ。だけど、apply_along hogehogeもどうも思い浮かばなかった。
そこで、こうする。

```python
matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)

num_arr = tile(i, (matrix.shape[0], 1))
result = matrix * num_arr
```
勿論、圧倒的に速い。しかし、どうもnumpy芸である…
そこで、関数にしてしまうのはどうか。
ここでは、先に書いたnp.meshgridグラム行列風のインターフェイスを目指してみる。

```python
def for_each_row(x, y):
    return np.tile(x, (y.shape[0], 1)), y

matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)
def func(x, y): return x * y

result = func(*for_each_row(nums, matrix))
```

ちなみに、実はこのケースではこれで良い。

```python
matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)
matrix * nums
```
上記のやり方は、より複雑なものを実装しないといけなくなった時に役立つ…かもしれない。

# 逆問題
## 連立方程式を解く(解ける場合)
普通の連立方程式。例えばこう。
$5x + 3 = y$
$x + 1 = y$
行列表現ではこう。
$$AX=Y$$
これを解く場合は下記のようにする。

```python
X = np.linalg.solve(A, Y)
```
この場合、未知の数は変数Xだった。
しかし、実用的には測定値Xから定数Aを知りたい事がほとんどと思う。
だから、最小二乗法なんかではここのAとXが逆の立場になる。
逆行列は明示的に出さずにこの関数を使うほうが良いらしい。

## 最小二乗法を解く(解がない場合)
解が存在しないパターンの場合は最小二乗法を使うことになる。
「解が存在しない」とは例えば下記の様なやつ。

$8x + 4 = y$
$5x + 3 = y$
$x + 1 = y$

この連立方程式は互いに矛盾する解となるため、解がない。
最小二乗法の解は以下の式なので

$$X^TXA = X^TY$$


```python
A = np.linalg.solve(X.T @ X, X.T @ Y)
```

疑似逆行列を使うことも出来る。

```python
A = pinv(x) @ y
```
なんだけど、既にちゃんと関数があるから上記ではなく下記が正解。

```python
A = np.linalg.lstsq(X, Y)[0]
```

## 正則化付き最小二乗法

正則化付きの場合はどうするのかよく分からない。
ググったけど、上記のような正解を見つけられなかったから、適当に自分で書くことになるのかな？
$$(X^TX + \lambda I)A = X^TY$$

```python
lmd = 1
seisoku = lmd * np.eye(min(X.shape))
A = np.linalg.solve(X.T @ X + seisoku , X.T @ Y)
```

## 最小ノルム解を解く(解が無限にある場合)
このような時は解が無限にある。
$8x_1 + x_2 + 6x_3 = 4$
$5x_1 + x_2 + 2X_3 = 3$
この方程式の最適解を求めてみる。
$$A = X^T{(XX^T)}^{-1}Y$$

```python
A = X.T @ np.linalg.inv(X @ X.T) @ Y
```
ってなりそうじゃん？
だけど、ムーアペンローズの逆行列は凄いので、普通にいける。

```python
A = pinv(X) @ Y
```

でも、これも正解は下記。

```python
A = np.linalg.lstsq(X, Y)[0]
```
この関数はすごいねー。

## 正則化付き最小二乗法(解が無限にある場合)
これも正攻法が分からん。だけど、式は簡単に作れる。
$$A = X^T{(XX^T + \lambda I)}^{-1}Y$$

```python
lmd = 1
seisoku = lmd * np.eye(min(X.shape))
A = X.T @ np.linalg.inv(X @ X.T + seisoku) @ Y
```

うーん…(･ัω･ั)
正則化付き探せなかったのは僕の勉強不足か？
大して難しくないから良いんだけどさ…

# その他、いろんな最適化問題
https://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize


# 速くする
numpyでは```**```とかの演算子は比較的遅い。```e**hoge```とかやるべきではなくて、np.expとかを使うほうが速い。
僕の書いたWavelet変換のパッケージではそうだった。(なので可読性が悪い…)

```python
hoge = np.arange(1000)
e ** hoge  # 遅い
np.exp(hoge)  # 速い

hoge ** 2  # 遅い
np.square(hoge)  # 速い

hoge ** (0.5)  # 遅い
np.sqart(hoge)  # 速い

hoge ** (3)  # 遅い
np.power(hoge, 3)  # 速い

hoge ** (3.4)  # 遅い
np.float_power(hoge, 3.4)  # 速い
```

この他、np.powerとか色々あるけど、基本、最適化された関数を使うほうが速い。

# 楽をするために
## 行列の積
何故皆@演算子使わないのかよくわからない。np.dotとかより読みやすいと思うんだよな…。
遅かったりするんですか？

```python
np.arange(10).T@np.arange(10)
```

# 複素数関連
## 複素共役
np.conjでおｋ

## 実軸虚軸関連
hoge.real, hoge.imag
という感じで書けますから、こいつを代入する時は

```python
fuga.real, fuga.imag = hoge.real, hoge.imag
```

何故これが欲しいかと言うと、虚軸だけ、実軸だけに何らかの処理をすることが意外とあるから。

## フーリエ変換
フーリエ変換はどっちかと言うとscipyを使うことが多いと思う。numpyにもあるんだけどね。

```python
scipy.fftpack.fft(wave)
```

逆変換

```python
comp  # フーリエ変換後の複素数の行列とする
scipy.fftpack.ifft(comp)
```

# cupy関連
cupyへ変換

```python
cp.asarray(numpy_array)
```

numpyへ変換

```python
cp.asnumpy(cupy_array)
```

cupy自体、numpyと違う所もそれなりにあるので、色々気をつけようね。

# 場合分け
## 場合分け
np.whereを使う。
もしもxが条件式に合致すれば第二引数を、合致しなければ第三引数を出す。

```python
x = np.arange(10)
np.where(x > 5, x * 10, x)
> array([0, 1, 2, 3, 4, 5, 60, 70, 80, 90])
```

## 条件文
条件式をぶっこむとRっぽい挙動になる。

```python
x = np.arange(10)
x > 5
> array([False, False, False, False, False, False, True, True, True, True])
```

# 数学の基礎っぽいやつ
## 行列に関して
```python
np.det(hoge)  # 行列式
np.norm(hoge)  # ノルム
```

## 積分
これもscipy
これ、マジカヨすげーってなった。

```python
from scipy.integrate import quad
def test(x):
    return x ** 2 - x
print(quad, test, -5, 5)
```
他に、2回積分するのとかもある。それはdblquad


## 特殊関数
さては君は数学好きだね？
このリンクを見ると良い。
https://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special

# ゼロ埋め
hogeは配列
次のは3番目以前と5番目移行を埋めるということ。
constantは定数で埋めるという意味。
最後のは左右で埋める定数。
他にも色々モードがある。

```python
np.pad(hoge, [3, 5], 'constant' ,[0, 0])
```

# 応用
## 画像のくり抜き
値の大きな所だけくり抜いてみる。
さっきのようにTrue, Falseの配列を作って、条件式でくり抜く。
例ではnumpyで山みたいな物を作った上で、山のてっぺんを切り取っている。
この時に、np.nanが必要になったりする。matplotlibのimshowではnanが無視されるからだ。

```python
def mount(x: int, y: int) -> np.ndarray:
    return np.exp(-(np.square(x) + np.square(y)) / 1000)


# くり抜く関数。
# valuesが普通の行列で、clusterがboolの行列
def make_single_cluster(values, cluster):
    img = np.nan * np.zeros_like(values)
    img[cluster] = tvalues[cluster]
    return img


x_y = 100, 100
# 今回は一寸ノイズを入れてみた
X = np.random.rand(*x_y) / 100
# 山を作る
X += mount(*np.meshgrid(*(np.arange(-n/2, n/2, 1) for n in x_y)))

plt.imshow(make_single_cluster(X, X > 0.5))
plt.show()
```
![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/139153/86370117-9f62-c82a-98c9-d75c568b5e71.png)

# 互換性
## matlab読み書き
scipyにmatlab読み書き機能がある。

```python
import scipy.io as sio
sio.loadmat('hoge.mat')

data = {'fuga', np.arange(100)}
sio.savemat(fuga.mat, data)
```

残念ながら僕はmatlab使ったことないからよく分からん。