2018/2/28にPydata Osakaで'SciPyの概要と各モジュールの紹介'というタイトルで発表しました

from IPython.display import IFrame
IFrame("https://www.hiromasa.info/slide/4.slides.html", "100%", "450px")

# azure notebooksデフォルトのSciPyをpipで更新
!pip install -U scipy

# 本発表では以下を仮定
import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt

# version
print("numpy version", np.version.full_version)
print("scipy version", sp.version.full_version)
print("matplotlib version", mpl.__version__)

# 本発表では以下を仮定
import scipy.io as spio
import scipy.special as special
import scipy.linalg as la
import scipy.fftpack as fft
import scipy.optimize as opt
import scipy.signal as sg
import scipy.interpolate as interp
import scipy.stats as stats

a = np.ones((3, 3))
spio.savemat('file.mat', {'a': a}) # savematには辞書を渡す
data = spio.loadmat('file.mat')
data['a']

a = np.ones(3)
spio.savemat('file.mat', {'a': a})

print(a)
print(spio.loadmat('file.mat')['a']) # matlabには1次元配列がない
print(spio.loadmat('file.mat', squeeze_me=True)["a"])

special?

%load_ext cython

%%cython -a

cimport scipy.special.cython_special as csc

cdef:
    double x = 1
    double complex z = 1 + 1j
    double si, ci, rgam
    double complex cgam

rgam = csc.gamma(x)
print(rgam)
cgam = csc.gamma(z)
print(cgam)
csc.sici(x, &si, &ci)
print(si, ci)

# 行列式の計算
A = np.array([[1, 2],
              [3, 4]])
print(la.det(A))  # 1*4-2*3=-2

A = np.array([[3, 2],
              [6, 4]])
print(la.det(A) )  # 3*4-2*6=0

print(la.det(np.ones((3, 4))))  # 正方行列ではない

# 逆行列の計算
A = np.array([[1, 2],
              [3, 4]])
A_inv = la.inv(A)
print(A_inv)
print(A @ A_inv)  # AA^-1 ~= I
print(np.allclose(A @ A_inv, np.eye(2)))

# 特異行列(行列式が0)の場合は逆行列が求められないのでLinAlgError例外を投げる
A = np.array([[3, 2],
              [6, 4]])
la.inv(A)  

import numpy as np
from scipy import linalg
A = np.array([[1, 2], [3, 4]])
b = np.array([[5], [6]])

print("逆行列による解法")
print(linalg.inv(A) @ b)  # 遅い
print((A @ linalg.inv(A) @ b) - b)  # 確認

print("solveによる解法")
print(np.linalg.solve(A, b))  # 高速
print((A @ np.linalg.solve(A, b)) - b)  # 確認

A = np.array([[1,2],[3,4]])
print(A)
print(la.norm(A))
print(la.norm(A, 'fro'))  # デフォルト
print(np.sqrt(np.trace(A.T @ A)))
print(la.norm(A,1))       # L1ノルム(各列の絶対値和のうち最大のもの)
print(la.norm(A,-1))      # L1ノルム(各列の絶対値和のうち最小のもの)
print(la.norm(A,np.inf))  # L∞ノルム(各行の絶対値和のうち最大のもの)

c1, c2 = 5.0, 2.0
i = np.r_[1:11]
xi = 0.1*i  # 説明変数
yi = c1*np.exp(-xi) + c2*xi
zi = yi + 0.05 * np.max(yi) * np.random.randn(len(yi))  # 観測値

A = np.c_[np.exp(-xi)[:, np.newaxis],
          xi[:, np.newaxis]]  # 計画行列
print(A)

# 計画行列と観測値を与えてモデルパラメータを推定
c, resid, rank, sigma = linalg.lstsq(A, zi)

# パラメータ、残差の総和、Aのrank, Aの特異値
print(c, resid, rank, sigma)

xi2 = np.r_[0.1:1.0:100j]
yi2 = c[0]*np.exp(-xi2) + c[1]*xi2
yi_true = c1*np.exp(-xi2) + c2*xi2

plt.plot(xi,zi,'x',xi2,yi2, xi2, yi_true)
plt.axis([0,1.1,3.0,5.5])
plt.xlabel('$x_i$')
plt.title('Data fitting with linalg.lstsq')
plt.show()

A = np.array([[1, 2], [3, 4]])
(l1, l2), v = la.eig(A)

print(l1, l2)   # 固有値
print(v)        # 各列ベクトルが固有ベクトル
print(la.norm(v, axis=0))  # 固有ベクトルは正規化されている

# Ax = λxとなることを確認
# |1 2|  =  |v1 v2 ||l1  0| = |l1*v1 l2*v2|
# |3 4|            | 0 l2| 
print(A@v - v@np.diag((l1,l2)))

# 正方行列に対する特異値分解
A = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1])
U, s, Vh = la.svd(A)
Sig = np.diag(s)
print(A)
print(U.shape, s.shape, Sig.shape, Vh.shape)
print(U @ Sig @ Vh)

# 非正方行列に対する特異値分解
A = np.array([[1,2,3],[4,5,6]])
M,N = A.shape

U,s,Vh = linalg.svd(A)
Sig = la.diagsvd(s, M, N)
print("")
print(A)
print(U.shape, s.shape, Sig.shape, Vh.shape)
print(U@Sig@Vh)

import IPython.display

fs = 48000.0; nfft = 2048
xs = np.arange(1.0 * fs)  # 1[s]
ts = xs / fs
hz = 2000.0
ys = np.sin(2.0 * np.pi * hz * ts) + 0.1 * np.random.randn(int(fs))

fig = plt.figure(1, figsize=(12,4))
ax = fig.add_subplot(1, 1, 1)
ax.plot(ts[:nfft], ys[:nfft])
ax.set_xlabel("time [s]"); ax.set_ylabel("amplitude")

IPython.display.Audio(data=ys, rate=fs)

# 窓関数を取得
win = sg.get_window("hanning", nfft)

# FFTを実行
Ys = fft.fft(win * ys[:nfft], nfft)[:nfft//2+1]

# 周波数軸を取得
f = fft.fftfreq(nfft, d=1.0/fs)[:nfft//2+1]
f[-1] = fs//2  # デフォルトだと-24kHz

fig = plt.figure(1, figsize=(12,4));
ax = fig.add_subplot(1, 2, 1)
ax.plot(ts[:nfft], win * ys[:nfft])
ax.set_xlabel("time [s]"); ax.set_ylabel("amplitude")
ax = fig.add_subplot(1, 2, 2)
ax.semilogx(f, np.abs(Ys))
ax.set_xlabel("frequency [Hz]"); ax.set_ylabel("amplitude");

# 観測値
xs = np.linspace(0, 1, 10)
ys = np.sin(2.0 * np.pi * xs) + 0.01 * np.random.randn(len(xs))

# 異なる2種類の補間関数の生成
linear_interp = interp.interp1d(xs, ys)
cubic_interp = interp.interp1d(xs, ys, kind='cubic')

# 補間値を算出
xs_new = np.linspace(0, 1, 50)
ys_lin = linear_interp(xs_new)
ys_cubic = cubic_interp(xs_new)

fig = plt.figure(figsize=(12,6)); ax = fig.add_subplot(111)
ax.plot(xs, ys, "x-", label="obserbed")
ax.plot(xs_new, ys_lin, label="linear")
ax.plot(xs_new, ys_cubic, label="cubic")
ax.grid(); ax.legend()

# ホワイトノイズの生成
time_s = 5.0
ys = 0.1 * np.random.randn(int(time_s * fs))
IPython.display.Audio(ys, rate=fs)

nlpf = 96
cutoff_hz = 4000.0
lpf = sg.firwin(nlpf, cutoff=cutoff_hz, fs=fs)
ys_lpf = sg.lfilter(lpf, [1.0], ys)
IPython.display.Audio(ys_lpf, rate=fs)

fig = plt.figure(1, figsize=(12,4));
ax = fig.add_subplot(1, 2, 1)
ax.plot(f, 20.0*np.log10(np.abs(fft.fft(ys, nfft)[:nfft//2+1])))
ax.set_xlabel("time [s]"); ax.set_ylabel("power [dB]")

ax = fig.add_subplot(1, 2, 2)
ax.plot(f, 20.0*np.log10(np.abs(fft.fft(ys_lpf, nfft))[:nfft//2+1]))
ax.set_xlabel("frequency [Hz]"); ax.set_ylabel("power [dB]")

# ノイズの乗った正弦波を観測値とする
x_data = np.linspace(-5, 5, num=50)
y_data = 2.9 * np.sin(1.5 * x_data) + np.random.normal(size=50)

def test_func(x, a, b):
    return a * np.sin(b * x)

# curve_fitを用いてa, bを推定
params, params_covariance = opt.curve_fit(test_func, x_data, y_data, p0=[2, 2])
print(params, params_covariance)

fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(111)
ax.plot(x_data, y_data, "x")
ax.plot(x_data, test_func(x_data, params[0], params[1]))

# 最小値を見つけたい対象
def f(x):
    return x**2 + 10*np.sin(x)

# 大域解が-1.3, 局所解が3.8
x = np.arange(-10.0, 10.0, 0.1)
fig = plt.figure(); ax = fig.add_subplot(111)
ax.plot(x, f(x))

result = opt.minimize(f, x0=0.0)  # x0で探索開始点を与える
print(result)
print(result.x)

# nfev(評価回数)が12に減る
opt.minimize(f, x0=0, method="L-BFGS-B")

# 3.0から始めると局所解に陥る
result = opt.minimize(f, x0=3.0, method="L-BFGS-B")
print(result)
print(result.x)

# `scipy.optimize.basinhopping`を用いて大域解を探す
opt.basinhopping(f, 1.0)

# 探索範囲に制約を与えることも可能
res = opt.minimize(f, x0=1,
                   bounds=((0, 10), ))
res.x    

# f=0となる点をoptimize.rootで探す
root = opt.root(f, x0=1)
print(root)
print(root.x)

# もうひとつの根を探す
root = opt.root(f, x0=-3.0)
print(root)
print(root.x)

# サンプリングした点のヒストグラム(頻度分布)と確率分布を同時に描画
samples = np.random.normal(size=1000) # sampling

bins = np.arange(-4, 5)
histogram = np.histogram(samples, bins=bins, normed=True)[0] # plot histogram

bins = 0.5*(bins[1:] + bins[:-1])
pdf = stats.norm.pdf(bins)  # normは正規分布の意味

fig = plt.figure(); ax = fig.add_subplot(111)
ax.plot(bins, histogram, label="histgram")
ax.plot(bins, pdf, label="probability density")
ax.grid(); ax.legend()

# 観測値を確率分布にフィッティング
loc, std = stats.norm.fit(samples)
print(loc, std)
xs = np.linspace(-3.0, 3.0, 1000)
f = lambda x: (1.0/np.sqrt(2.0*np.pi*std))*np.exp(-(x - loc)**2/ (2.0 * std**2))
plt.plot(xs, f(xs))

# 平均、中央値、四分位数
print(np.mean(samples))
print(np.median(samples))
print(stats.scoreatpercentile(samples, 50))
print(stats.scoreatpercentile(samples, 90))

from scipy.spatial import distance
coords = [(35.0456, -85.2672),
          (35.1174, -89.9711),
          (35.9728, -83.9422),
          (36.1667, -86.7833)]

print(np.sqrt((coords[0][0] - coords[1][0])**2
             +(coords[0][1] - coords[1][1])**2))
distance.cdist(coords, coords, 'euclidean')

# pdistは対角を除いた上三角成分のみ返す
presult = distance.pdist(coords)
print(presult)

# squareformで距離行列に変換
print(distance.squareform(presult))

import scipy.misc
import matplotlib.cm as cm
face = scipy.misc.face()

print(face.shape)
print(face.max())
print(face.dtype)

fig = plt.figure(1)
ax = fig.add_subplot(1, 1, 1)
ax.imshow(face)

plt.show()

face = scipy.misc.face(gray=True)

fig = plt.figure(1)
ax = fig.add_subplot(1, 1, 1)
ax.imshow(face, cmap=cm.gray)

plt.show()

# グレースケールのデモ用画像
ascent = scipy.misc.ascent()
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.imshow(ascent, cmap=cm.gray)
plt.show()

from scipy.misc import derivative

# x0の点におけるn次の中心差分を返す
# n==1の場合はf(x+h)-f(x-h)/2h
# 内部的にcentral_diff_weightを使う

# derivative(func, x0, dx=1.0, n=1, args=(), order=3)
f = lambda x: x**3 + x**2
derivative(f, 1.0, dx=1e-6)

print((f(1.0 + 1e-6) - f(1.0 - 1e-6)) / (2.0 * 1.e-6))

!pip install librosa

import librosa
import librosa.display

y, sr = librosa.load(librosa.util.example_audio_file())
IPython.display.Audio(data=y, rate=sr)

np.pad(np.arange(5), 2, "reflect")

# strideトリックを使った信号を短時間波形フレームの配列に変換する関数
from numpy.lib.stride_tricks import as_strided

def split_frames(xs, frame_width=1024, noverlap=2):

    nsample = np.alen(xs)
    nshift = (frame_width // noverlap)
    nframe = int(nsample / nshift) - noverlap + 1
    
    nbyte = xs.dtype.itemsize

    # as_strided(xs, ys_shape, new_stride)
    xs_split = as_strided(xs, (nframe, frame_width), (nbyte * nshift, nbyte))
    
    return xs_split

# 信号の短時間波形への切り出し
import matplotlib.gridspec as gridspec

fig = plt.figure(1, figsize=(12,8))
gs = gridspec.GridSpec(2, 21)

ax = fig.add_subplot(gs[0, :-1])
xs = np.arange(len(y)) / sr
ax.plot(xs[:sr], y[:sr])
ymax = np.max(y[:sr])
ax.set_xlim([0.0, 1.0]); ax.set_ylim([-ymax, ymax]);
ax.set_xlabel("time [s]"); ax.set_ylabel("amplitude"); ax.grid()

frame_width = sr * 0.1
xs = np.arange(frame_width) / sr
for i, y_cut in enumerate(split_frames(y, int(frame_width), 1)[:10]):
    ax = fig.add_subplot(gs[1, 2*i:2*(i+1)])
    ax.plot(xs, y_cut)
    ax.set_ylim([-ymax, ymax]); ax.grid()

# 切り出し間隔を1/2オーバーラップした場合
import matplotlib.gridspec as gridspec

fig = plt.figure(1, figsize=(12,8))
gs = gridspec.GridSpec(3, 21)

ax = fig.add_subplot(gs[0, :-1])
ax.plot(y[:sr])
ymax = np.max(y[:sr])
ax.set_ylim([-ymax, ymax]); ax.grid()

frame_width = sr * 0.1
half_width = int(frame_width // 2)
xs = np.arange(frame_width)
for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
    if i % 2 == 0:
        ax = fig.add_subplot(gs[1, i:i+2])
        ax.plot(xs[:half_width], y_cut[:half_width])
        ax.plot(xs[half_width:], y_cut[half_width:])
    else:
        ax = fig.add_subplot(gs[2, i:i+2])
        ax.plot(xs[half_width:], y_cut[half_width:])
        ax.plot(xs[:half_width], y_cut[:half_width])

    ax.set_ylim([-ymax, ymax]); ax.grid()

frame_width = int(sr * 0.1)
win = sg.get_window("hanning", frame_width)
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.plot(win)
ax.grid()

# 更に各短時間波形に窓掛けを行った場合
import matplotlib.gridspec as gridspec

fig = plt.figure(1, figsize=(12,8))
gs = gridspec.GridSpec(3, 21)

ax = fig.add_subplot(gs[0, :-1])
ax.plot(y[:sr])
ymax = np.max(y[:sr])
ax.set_ylim([-ymax, ymax]); ax.grid()

frame_width = sr * 0.1
half_width = int(frame_width // 2)
win = sg.get_window("hanning", int(frame_width))
xs = np.arange(frame_width)
for i, y_cut in enumerate(split_frames(y, int(frame_width), 2)[:20]):
    y_win = win * y_cut
    if i % 2 == 0:
        ax = fig.add_subplot(gs[1, i:i+2])
        ax.plot(xs[:half_width], y_win[:half_width])
        ax.plot(xs[half_width:], y_win[half_width:])
    else:
        ax = fig.add_subplot(gs[2, i:i+2])
        ax.plot(xs[half_width:], y_win[half_width:])
        ax.plot(xs[:half_width], y_win[:half_width])

    ax.set_ylim([-ymax, ymax]); ax.grid()

# 1/4オーバーラップして窓掛けした各波形をそれぞれFFT
nfft = 2048
win = sg.get_window("hanning", nfft)
y_center = np.pad(y, int(nfft//2), "reflect")
fft_list = []
phase_list = []
for y_cut in split_frames(y_center, nfft, 4):
    y_end = min(len(y_cut), nfft)
    y_win = np.zeros(nfft)
    y_win[:y_end] = win[:y_end] * y_cut[:y_end]
    Y = fft.fft(y_win, nfft)
    fft_list.append(Y[:nfft//2+1])
    phase_list.append(np.exp(1.j * np.angle(Y[:nfft//2+1])))

my_D = np.array(fft_list, dtype=np.complex64).T
my_phase = np.array(phase_list, dtype=np.complex64).T

fig = plt.figure(1, figsize=(12,4))
ax = fig.add_subplot(2,1,1)
librosa.display.specshow(librosa.amplitude_to_db(np.abs(my_D), ref=np.max), y_axis="log")

# 得られたスペクトログラムnp.abs(my_D.T)に対してNMFを実行
from sklearn.decomposition import NMF
model = NMF(n_components=8)
W = model.fit_transform(np.abs(my_D.T))
H = model.components_

components = H.T
activations = W.T

fig = plt.figure(2, figsize=(12,4))

ax = fig.add_subplot(1,2,1)
librosa.display.specshow(activations, x_axis='time')
ax.set_xlabel('Time'); ax.set_ylabel('Component'); ax.set_title('Activations')

ax = fig.add_subplot(1,2,2)
librosa.display.specshow(librosa.amplitude_to_db(components, ref=np.max), y_axis='log')
ax.set_xlabel('Component'); ax.set_ylabel('Frequency'); ax.set_title('Components')

# 1.e-12以下となっている要素の数をコンポーネントごとに算出
for i in range(8):
    print(sum(np.where(components[:, i] < 1.e-12, 1, 0)))

# Play back the reconstruction
# Reconstruct a spectrogram by the outer product of component k and its activation
D_k = components @ activations

# invert the stft after putting the phase back in
#y_k = librosa.istft(D_k * phase)
y_k = librosa.istft(D_k * my_phase)

# And playback
print('Full reconstruction')

IPython.display.Audio(data=y_k, rate=sr)

components_ = np.zeros(components.shape)
components_[:, 0] = components[:, 0]
components_[:, 2] = components[:, 2]
components_[:, 6] = components[:, 6]

D_k = components_ @ activations

# invert the stft after putting the phase back in
y_k = librosa.istft(D_k * my_phase)

# And playback
print('Partly reconstruction')

IPython.display.Audio(data=y_k, rate=sr)

from scipy.ndimage.interpolation import shift
xs = np.array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.]) 
print(shift(xs, -3, cval=0.0))
print(np.roll(xs, -3))