Pythonの最適化ライブラリーOptunaを使ったあとその結果を可視化してグラフに保存したいと思ったのでそのメモ．

必要なライブラリーのインストール

今回必要になるライブラリは Optuna，それと可視化のためのPlotlyそしてKaleidoの３つ．KaleidoはPlotlyのグラフを保存するのに必要みたいで，直接呼び出したりはしない．ではそれぞれpipでインストールしておく

$ pip install optuna
$ pip install plotly 
$ pip install kaleido

Plotlyの簡単な使い方は過去の記事にある． Plotlyで簡単なグラフを描いてみる（Python） - My Memo

Optunaで最適化

次はOptunaで最適化を行う．今回はz=x²+y²+xyが最小になるように x, yを最適化する（最小/最大化させたい関数zを目的関数という）．ちなみに，z=x²+y²+xyはx=0, y=0で最小となるので，０に近いx, yが得られれば最適化できたことになる．ちなみにこの式のグラフはこんな感じ．

この図の描画のPlotlyでの仕方は後述する．まずは，今回必要なライブラリをインポートする．

import optuna
import numpy as np
import plotly.graph_objects as go

次に最小化（あるいは最大化）する関数を書く．今回はz=x²+y²+xy．

def func(x, y):
    return x**2+y**2+x*y

次に，最適化したい変数の範囲を設定して目的関数のラッパーを定義する．最適化したい変数の範囲に関して，今回はx, yともに 0から1の値を取るようにする．

def objective(trial):
    x = trial.suggest_uniform('x', -1, 1)  # xの範囲指定
    y = trial.suggest_uniform('y', -1, 1)  # yの範囲指定
    return func(x, y)

この関数は必ずtrialを引数にとり，最適化したい目的関数を返すようにする．なお範囲を指定するとき，整数のみをとることや，リストの要素から選ぶようにすることも可能

メソッド
suggest_categorical(name, choices)	choicesというリストの中から選ぶ	x=trial.suggest_categorical('x', ['a', 'b', 'c'])
suggest_int(name, low, high[, step, log])	整数から選ぶ	x=trial.suggest_categorical('x', 1, 10, 2)　1から10まで2刻み

のように．

そして，目的関数と変数の範囲を指定したら最適化する．

if __name__ == "__main__":
    # search optimal x and y
    study = optuna.create_study(direction='minimize')
    study.optimize(objective, n_trials=50)

今回は目的関数の値を最小化させたいので，direction='minimize'と書いたが，目的関数を最大化させたい場合はdirection='maximize'とすると良い．なおdirectionを省略した場合自動的にdirection='minimize'になる．study.optimize()では目的関数と試行回数を引数として渡す．これで実行すれば最適化をしてくれる．実行すると標準出力に結果が逐次表示されるが，optunaの標準出力をファイルに保存したい場合は実行時に

$ python optimize.py &> output.txt

とすればファイルに保存できる．＆を忘れずに．

結果の可視化

ここまで来るとあとはかんたん．最適化したあとの数値的な結果はstudyが持っているので数値はそれを出せばよい．以下のコードで標準出力してくれる．

    # print result
    print(study.best_params)  # x, yの最適な値
    print(study.best_value)  #  zの最適値
    print(study.best_trial)  # x, y, zの最適値

次は結果のグラフを出力する．まずは， x，yの値に対するzの値を等高線図で表示する．コードは以下．

    # print result
    # contour plot
    fig = optuna.visualization.plot_contour(study)
    fig.write_html('image/contour.html')  # save as html file
    fig.write_image('image/contour.png')  # save as png file

optuna.visualization.plot_contour(study)メソッドは Plotlyのグラフインスタンスを返すので，保存するには上のようにする． Plotlyのグラフインスタンスの扱い方の詳細は Plotlyで簡単なグラフを描いてみる（Python） - My Memo

Plotlyで線グラフをたくさん描く - My Memo

するとこのようなグラフが保存される．

次はどの変数が最も最小化に寄与したか（今回はx，yは対称なので，理想的には同じ数値になる）．

    # Hyperparameter importances
    fig = optuna.visualization.plot_param_importances(study)
    fig.write_html('image/param_importance.html')  # save as html file
    fig.write_image('image/param_importance.png')  # save as png file

グラフは

となる．他にも色々かける

    # Empirical Distribution Function Plot
    fig = optuna.visualization.plot_edf(study)
    fig.write_html('image/edf.html')  # save as html file
    fig.write_image('image/edf.png')  # save as png file

    # Optimization History Plot
    fig = optuna.visualization.plot_optimization_history(study)
    fig.write_html('image/optimization_history.html')  # save as html file
    fig.write_image('image/optimization_history.png')  # save as png file

    # parallel coordinate plot
    fig = optuna.visualization.plot_parallel_coordinate(study)
    fig.write_html('image/parallel_coordinate.html')  # save as html file
    fig.write_image('image/parallel_coordinate.png')  # save as png file

    # slice
    fig = optuna.visualization.plot_slice(study)
    fig.write_html('image/slice.html')  # save as html file
    fig.write_image('image/slice.png')  # save as png file

それぞれグラフは

といろいろ描ける．

今回のコード全体

今回のコードの全体を貼っておく．なお，目的関数の三次元グラフ描画する部分も付け足しておいた．

import optuna
import numpy as np
import plotly.graph_objects as go


def func(x, y):
    return x**2+y**2+x*y


def objective(trial):
    x = trial.suggest_uniform('x', -1, 1)
    y = trial.suggest_uniform('y', -1, 1)
    return func(x, y)


if __name__ == "__main__":
    # plot func(x, y) by plotly
    x = np.linspace(-1, 1, 100)
    y = np.linspace(-1, 1, 100)
    xx, yy = np.meshgrid(x, y)
    z = func(xx, yy)
    fig = go.Figure()
    fig.add_trace(go.Surface(z=z, x=x, y=y))
    fig.write_html('objective.html')

    # search optimal x and y
    study = optuna.create_study(direction='minimize')
    study.optimize(objective, n_trials=50)

    # print result
    print(study.best_params)
    print(study.best_value)
    print(study.best_trial)

    # visualize
    # contour plot
    fig = optuna.visualization.plot_contour(study)
    fig.write_html('image/contour.html')  # save as html file
    fig.write_image('image/contour.png')  # save as png file

    # Empirical Distribution Function Plot
    fig = optuna.visualization.plot_edf(study)
    fig.write_html('image/edf.html')  # save as html file
    fig.write_image('image/edf.png')  # save as png file

    # Optimization History Plot
    fig = optuna.visualization.plot_optimization_history(study)
    fig.write_html('image/optimization_history.html')  # save as html file
    fig.write_image('image/optimization_history.png')  # save as png file

    # parallel coordinate plot
    fig = optuna.visualization.plot_parallel_coordinate(study)
    fig.write_html('image/parallel_coordinate.html')  # save as html file
    fig.write_image('image/parallel_coordinate.png')  # save as png file

    # Hyperparameter importances
    fig = optuna.visualization.plot_param_importances(study)
    fig.write_html('image/param_importance.html')  # save as html file
    fig.write_image('image/param_importance.png')  # save as png file

    # slice
    fig = optuna.visualization.plot_slice(study)
    fig.write_html('image/slice.html')  # save as html file
    fig.write_image('image/slice.png')  # save as png file

2021-09-28

線形計画問題を単体法（シンプレックス法）によって解くコードを Pythonで実装してみた

現在執筆中..

線形計画問題
単体法
単体法の実装
おまけ
- 線形計画問題の変形

線形計画問題

今回解く線形計画問題は不等式標準形とする．なお，任意の線形計画問題は不等式標準形に変換可能なので，このコードを適用させるには不等式標準形に変換（→おまけ）してから適用させる必要がある．ただし，任意の線形計画問題は不等式標準形に変形できるという事実があるので結構使えると思う．
不等式標準形の線形計画問題
$\displaystyle{\text 条件} A\bf{x}\ge\bf{b}, \bf{x}\ge \bf{0}{\text のもとで}\bf{c}^T\bf{x}{\text を最小にする}\bf{x}{\text を求める.}$
また，このような線形計画問題を以下のように表すことにする．
最小化: $\bf{c}^T\bf{x}$
条件: $A\bf{x}\ge\bf{b}, \bf{x}\ge\bf{0}$

具体的な例としては下の3次元空間での線形計画問題を取り扱う．
例1
最小化: $-x_1+2x_2-x_3$
条件: $-2x_1+2x_2\ge-6\\ -x_1-x_2+2x_3\ge-3\\ -x_1-x_2-x_3\ge-3\\x_1, x_2, x_3\ge0$
この場合
$\displaystyle{A=\left(\begin{array}{ccc} -2&2&0 \\-1&-1&2 \\-1&-1&-1\end{array}\right)}$ , $\displaystyle{\bf{b}=\left(\begin{array}{c} -6 \\-3 \\-3\end{array}\right)}$ , $\displaystyle{\bf{c}=\left(\begin{array}{c} -1 \\2 \\-1\end{array}\right)}$ , $\displaystyle{\bf{x}=\left(\begin{array}{c} x_1 \\x_2 \\x_3\end{array}\right)}$
とすれば，上の
最小化: $\bf{c}^T\bf{x}$
条件: $A\bf{x}\ge\bf{b}, \bf{x}\ge\bf{0}$
という形で表す事ができるとわかる．とにかく条件を満たす $x_1, x_2, x_3$ で目標関数を最小にするものを求めるということを行う．そのような線形計画問題のための解法の1つが単体法（シンプレックス法）である．というわけで最初に単体法の解説を簡単に書いておいて，その後pythonのコードに焼き直す．

単体法

単体法の

単体法の実装

とにかく上のアルゴリズムを Pythonでコーディングすると↓のようになる．

from copy import deepcopy
import numpy as np
import traceback


class Simplex():
    def __init__(self, A, b, c):
        self.A = A
        self.b = b
        self.c = c
        self.m, self.n = A.shape  # m: 方程式の数, n: 変数の次元

        # 辞書の定義
        self.dic = np.zeros((self.m + 1, self.n + 1), dtype='float64')
        for i in range(self.m):
            for j in range(self.n):
                self.dic[i + 1, j + 1] += A[i, j]
        for j in range(self.n):
            self.dic[0, j + 1] += c[j]
        for i in range(self.m):
            self.dic[i + 1, 0] += -b[i]
        self.x_n = np.zeros(self.n, dtype='float64')  # 非基底変数
        self.x_m = np.zeros(self.m, dtype='float64')  # 基底変数
        self.x_n_index = np.arange(0, self.n, dtype='int64')
        self.x_m_index = np.arange(self.n, self.m + self.n, dtype='int64')

    def run(self):
        # 原点が許容解であるか確認
        for j in range(self.n):
            self.x_n[j] = 0.0
        self.x_m = self.update_basis()  # 非基底変数が原点のときの基底変数を代入
        is_feasible = True
        for x in self.x_m:
            is_feasible = is_feasible and (x >= 0.0)
        if not is_feasible:
            try:  # 原点が非許容である場合エラー出力
                raise Exception('The origin is not feasible.')
            except BaseException:
                print(traceback.format_exc())

        z_sol = self.loop()

        x_sol = [0 for _ in range(self.n)]
        x_concat = np.concatenate([self.x_n, self.x_m])
        index_concat = np.concatenate([self.x_n_index, self.x_m_index])
        for j, x in zip(index_concat, x_concat):
            if j < self.n:
                x_sol[j] = x

        return z_sol, np.array(x_sol)

    def loop(self):
        while True:
            self.show_dic()
            # 非基底変数の選択
            j_pivot = -1
            for j, a in enumerate(self.dic[0, 1:]):
                if a < 0:
                    j_pivot = j + 1
                    break
            if j_pivot == -1:  # 負である係数が存在しない場合終了
                return self.get_cost()

            # 基底変数の選択
            base_candidate = {}
            for i in range(1, self.m + 1):
                if self.dic[i, j_pivot] < 0:
                    base_candidate[i] = self.dic[i, 0] \
                        / abs(self.dic[i, j_pivot])
            if len(base_candidate) == 0:
                # 基底の候補がない場合はこの線形計画問題は非有界
                raise Exception('This problem is unbounded.')
            d_min = min(base_candidate.values())
            i_min_list = [kv[0] for kv in base_candidate.items()
                          if kv[1] == d_min]
            # 基底の候補から添え字の最小となる基底を選択
            i_pivot = self.m + self.n
            index_pivot = self.m + self.n
            for i_min in i_min_list:
                if index_pivot >= self.x_m_index[i_min - 1]:
                    i_pivot = i_min
                    index_pivot = self.x_m_index[i_min - 1]

            # 選択された非基底変数を許容性を保存しつつ最大限まで増加させる
            self.x_n[j_pivot - 1] = d_min
            self.x_m[i_pivot - 1] = 0.0

            # ピボット演算でdic, x_m, x_n, x_m_index, x_n_indexを更新
            self.pivot(i_pivot, j_pivot)

    def pivot(self, i_pivot, j_pivot):
        self.x_m[i_pivot - 1], self.x_n[j_pivot - 1] = \
            self.x_n[j_pivot - 1], self.x_m[i_pivot - 1]
        self.x_m_index[i_pivot - 1], self.x_n_index[j_pivot - 1] = \
            self.x_n_index[j_pivot - 1], self.x_m_index[i_pivot - 1]

        dic_prev = deepcopy(self.dic)
        for i in range(self.m + 1):
            if i == i_pivot:
                c = -1.0 / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] * c
                self.dic[i_pivot, j_pivot] = -c
            else:
                c = dic_prev[i, j_pivot] / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] - dic_prev[i_pivot, :] * c
                self.dic[i, j_pivot] = c

    def update_basis(self):
        return np.dot(self.dic[1:, 1:], self.x_n) + self.dic[1:, 0]

    def get_cost(self):
        return np.dot(self.dic[0, 1:], self.x_n) + self.dic[0, 0]

    def show_dic(self):
        print('    |          |', end='')
        for x in self.x_n_index:
            print('%10.2f ' % (x), end='')
        print('')
        print('-' * (4 + 11 + 11 * self.n))

        print('cost|%10.2f|' % (self.dic[0, 0]), end='')
        for d in self.dic[0, 1:]:
            print('%10.2f ' % (d), end='')
        print('')

        for i in range(1, self.m + 1):
            print('%4d|%10.2f|' % (self.x_m_index[i - 1], self.dic[i, 0]), end='')
            for d in self.dic[i, 1:]:
                print('%10.2f ' % (d), end='')
            print('')
        print('')


# 実際に実行してみる
if __name__ == "__main__":
    A = np.array([[-2, -2, 1], [-2, 0, -4], [4, -3, 1]], dtype='float64')
    b = np.array([4, 4, 1], dtype='float64')
    c = np.array([-2, -1, -1], dtype='float64')
    print(Simplex(A, b, c).run())

まずはSimplexクラスのコンストラクタの説明から．

    def __init__(self, A, b, c):
        self.A = A
        self.b = b
        self.c = c
        self.m, self.n = A.shape  # m: 方程式の数, n: 変数の次元

とりあえず不等式標準形の各定数をメンバ変数として持っておく（結局使うことはなかったが...）．そして，行列Aの形状から方程式の数（不等式の数）： m，変数の次元：n，をメンバ変数に格納する．次に，辞書を定義する．

        # 辞書の定義
        self.dic = np.zeros((self.m + 1, self.n + 1), dtype='float64')
        for i in range(self.m):
            for j in range(self.n):
                self.dic[i + 1, j + 1] += A[i, j]
        for j in range(self.n):
            self.dic[0, j + 1] += c[j]
        for i in range(self.m):
            self.dic[i + 1, 0] += -b[i]
        self.x_n = np.zeros(self.n, dtype='float64')  # 非基底変数
        self.x_m = np.zeros(self.m, dtype='float64')  # 基底変数
        self.x_n_index = np.arange(0, self.n, dtype='int64')
        self.x_m_index = np.arange(self.n, self.m + self.n, dtype='int64')

まず，辞書の形状はm+1×n+1の行列であるので，そのような形状ですべての要素が0である行列self.dicを定義して，それにA, b, cの数値を代入して行っている．また，非基底変数と基底変数を宣言する．単体法はそのアルゴリズムを通して非基底変数と基底変数をスワップすることを何回も行う．したがって，現在の非基底変数または基底変数の何番目の変数が $x_1$ や $x_2$ に対応するのか記憶しておかなければならない．そこで，k番目（0オリジン）の非基底変数に対応する変数が $x_i$ であるとすると，self.x_n_index[k]=iというふうに，非基底変数に対してはself.x_n_index，基底変数に対してはself.x_m_indexに対応する変数の添字を記憶しておく．コンストラクタは以上で終了．
メソッドrun()にて単体法のアルゴリズムを実行する．

    def run(self):
        # 原点が許容解であるか確認
        for j in range(self.n):
            self.x_n[j] = 0.0
        self.x_m = self.update_basis()  # 非基底変数が原点のときの基底変数を代入
        is_feasible = True
        for x in self.x_m:
            is_feasible = is_feasible and (x >= 0.0)
        if not is_feasible:
            try:  # 原点が非許容である場合エラー出力
                raise Exception('The origin is not feasible.')
            except BaseException:
                print(traceback.format_exc())

        z_sol = self.loop()

ここでは最初に原点が許容解であるか確認する．（原点が許容解ではないならば二段階単体法を用いる必要がある．）まずは，全ての非基底変数に0を代入し，それと辞書を用いて基底変数の値を求める．そして，基底変数（今は1度もピボットしていないのでスラック変数）が負であるならば原点は非許容であるので，例外を吐いて終了する．許容解であるならば単体法のループ部分（loop()メソッド，戻り値：目的関数）に入る．なお，このコードでは辞書と非基底変数から基底変数，目的関数を求めるためにそれぞれ，update_basis(), get_cost()というメソッドを定義した．

    def update_basis(self):
        return np.dot(self.dic[1:, 1:], self.x_n) + self.dic[1:, 0]

    def get_cost(self):
        return np.dot(self.dic[0, 1:], self.x_n) + self.dic[0, 0]

さて，単体法の心臓部分のループ．なお，計算の過程がわかるように1ループごとに辞書の内容をメソッドshow_dic()にて出力する．

   def loop(self):
        while True:
            self.show_dic()
            # 非基底変数の選択
            j_pivot = -1
            for j, a in enumerate(self.dic[0, 1:]):
                if a < 0:
                    j_pivot = j + 1
                    break
            if j_pivot == -1:  # 負である係数が存在しない場合終了
                return self.get_cost()

            # 基底変数の選択
            base_candidate = {}
            for i in range(1, self.m + 1):
                if self.dic[i, j_pivot] < 0:
                    base_candidate[i] = self.dic[i, 0] \
                        / abs(self.dic[i, j_pivot])
            if len(base_candidate) == 0:
                # 基底の候補がない場合はこの線形計画問題は非有界
                raise Exception('This problem is unbounded.')
            d_min = min(base_candidate.values())
            i_min_list = [kv[0] for kv in base_candidate.items()
                          if kv[1] == d_min]
            # 基底の候補から添え字の最小となる基底を選択
            i_pivot = self.m + self.n
            index_pivot = self.m + self.n
            for i_min in i_min_list:
                if index_pivot >= self.x_m_index[i_min - 1]:
                    i_pivot = i_min
                    index_pivot = self.x_m_index[i_min - 1]

            # 選択された非基底変数を許容性を保存しつつ最大限まで増加させる
            self.x_n[j_pivot - 1] = d_min
            self.x_m[i_pivot - 1] = 0.0

            # ピボット演算でdic, x_m, x_n, x_m_index, x_n_indexを更新
            self.pivot(i_pivot, j_pivot)

最初に非基底変数の選択．選択できない場合，現在の値が最適解であるのでループ終了．次に基底変数の選択．ここで選択できない場合は線形計画問題が非有界であるため最適解はない．ループ終了．さて，無事に非基底変数・基底変数のペアが選択されたら次は基底変数が0になるまで非基底を増加させる．それ以上増加させたら基底がマイナスになって非許容となるためそこで止める．そしてピボットで諸々の値を更新する．
ピボットの部分はこんな感じ．

    def pivot(self, i_pivot, j_pivot):
        self.x_m[i_pivot - 1], self.x_n[j_pivot - 1] = \
            self.x_n[j_pivot - 1], self.x_m[i_pivot - 1]
        self.x_m_index[i_pivot - 1], self.x_n_index[j_pivot - 1] = \
            self.x_n_index[j_pivot - 1], self.x_m_index[i_pivot - 1]

        dic_prev = deepcopy(self.dic)
        for i in range(self.m + 1):
            if i == i_pivot:
                c = -1.0 / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] * c
                self.dic[i_pivot, j_pivot] = -c
            else:
                c = dic_prev[i, j_pivot] / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] - dic_prev[i_pivot, :] * c
                self.dic[i, j_pivot] = c

さて，ループが終わったらrun()メソッドに戻ろう．

        z_sol = self.loop()

        x_sol = [0 for _ in range(self.n)]
        x_concat = np.concatenate([self.x_n, self.x_m])
        index_concat = np.concatenate([self.x_n_index, self.x_m_index])
        for j, x in zip(index_concat, x_concat):
            if j < self.n:
                x_sol[j] = x

        return z_sol, np.array(x_sol)

ループが無事終了したら変数z_solに最適値が帰ってくるはずだ．そしてz_solを与える解をx_solに代入してそれらを返そう．x_solの求め方は基底・非基底に対応する添字を小さい方からn個取り出して．その値をx_solに格納して終了！（投げやり）
さて，上の例で示した線形計画問題を解いてみる．
例1
最小化: $-x_1+2x_2-x_3$
条件: $-2x_1+2x_2\ge-6\\ -x_1-x_2+2x_3\ge-3\\ -x_1-x_2-x_3\ge-3\\x_1, x_2, x_3\ge0$

if __name__ == "__main__":
    A = np.array([[-2, 2, 0], [-1, -1, 2], [-1, -1, -1]], dtype='float64')
    b = np.array([-6, -3, -3], dtype='float64')
    c = np.array([-1, 2, -1], dtype='float64')
    print(Simplex(A, b, c).run())

すると標準出力は

    |          |      0.00       1.00       2.00 
------------------------------------------------
cost|      0.00|     -1.00       2.00      -1.00 
   3|      6.00|     -2.00       2.00       0.00 
   4|      3.00|     -1.00      -1.00       2.00 
   5|      3.00|     -1.00      -1.00      -1.00 

    |          |      3.00       1.00       2.00 
------------------------------------------------
cost|     -3.00|      0.50       1.00      -1.00 
   0|      3.00|     -0.50       1.00       0.00 
   4|      0.00|      0.50      -2.00       2.00 
   5|      0.00|      0.50      -2.00      -1.00 

    |          |      3.00       1.00       5.00 
------------------------------------------------
cost|     -3.00|      0.00       3.00       1.00 
   0|      3.00|     -0.50       1.00      -0.00 
   4|      0.00|      1.50      -6.00      -2.00 
   2|      0.00|      0.50      -2.00      -1.00 

(-3.0, array([3., 0., 0.]))

こんな感じ．辞書の更新がされて最後に目的関数-3.0, 最適解 $(x_1, x_2, x_3)=(3, 0, 0)$ ということでいい感じに計算できてましたということでおしまい．

おまけ

線形計画問題の変形

等式標準形から，不等式標準形への変形を行う．
不等式標準形
最小化: $\bf{c}^T\bf{x}$
条件: $A\bf{x}\ge\bf{b}, \bf{x}\ge\bf{0}$
等式標準形
最小化: $\bf{\tilde{c}}^T\bf{\tilde{x}}$
条件: $A\bf{\tilde{x}}=\bf{\tilde{b}}, \bf{\tilde{x}}\ge\bf{0}$
等式標準系から不等式標準系への変形には次の式変形を用いる．

$\sum_{j=1}^na_jx_j=b \to \sum_{j=1}^na_jx_j\ge b,\ \sum_{j=1}^na_jx_j\le b$
$x （非負制限なし）\to x=x_1-x_2, x_1\ge0, x_2\ge0$

2021-05-03

サンプリング定理の証明

サンプリング定理の証明
- 証明
- エイリアシング
矩形波のフーリエ展開
低域通過フィルター

サンプリング定理の証明

　最近勉強しているとサンプリング定理（標本化定理）が出てきて，それについての証明をしてみたので未来の自分に向けてメモを残しておこうと思う．
　今回はこういうふうに考える．ある信号波 $S(t)$ に対して搬送波として矩形波 $p(t)$ を用いて変調を行い，PAM波 $f(t)=S(t)p(t)$ を得る（パルス振幅変調：PAM）．このPAM波を復調したときに信号波を完全に復元できるためには搬送波である矩形波の周波数は，
$\displaystyle{\omega_c\ge2\omega_{max}}$
（ $\omega_c$ ：矩形波の周波数． $\omega_{max}$ ：信号波に含まれる周波数の最大値）
を満たす必要がある．これを示してみる．

f:id:Ibetoge:20210503233513p:plain — 信号波，矩形波，PAM波

証明

　まず，搬送波として振幅 $A$ ，周期 $T=2\pi/\omega_c$ ，一周期内での非ゼロの時間 $\tau$ である矩形波 $p(t)$ を考える．

f:id:Ibetoge:20210503215342p:plain — 矩形波

この矩形波をフーリエ級数展開すると，以下のようになる（矩形波のフーリエ級数展開のやり方はこのページの後ろを参照）．
$\displaystyle{p(x)=\frac{A\omega_c\tau}{\pi}\left(1+\sum_{m=1}^{\infty}Sinc(m\pi\tau/2)\cos(m\omega_ct)\right)}$
　さらに，信号波をフーリエ展開すると以下のようになるとする．
$\displaystyle{S(t)=\sum_1^N\left\{\alpha_n\cos(n\omega_st)+\beta_n\sin(n\omega_st)\right\}}$
このとき信号波に含まれる周波数の最大値は， $\omega_{max}=N\omega_s$ となる．
そうするとPAM波が得られ，
$\displaystyle{ f(t)=S(t)p(t)=S(t)\frac{A\omega_c\tau}{\pi}\left(1+\sum_{m=1}^{\infty}Sinc(m\pi\tau/2)\cos(m\omega_ct)\right)}$

$\displaystyle{ =\frac{A\omega_c\tau}{\pi}S(t)+\frac{A\omega_c\tau}{\pi}S(t)\sum_{m=1}^{\infty}Sinc(m\pi\tau/2)\cos(m\omega_ct)}$
となる．ここで，式を見やすくするために $B=\frac{A\omega_c\tau}{\pi}, a_m=Sinc(m\pi\tau/2)$ と置いて，式の第二項を展開すると，

$\displaystyle{ f(t)=BS(t)+B\sum_{m=1}^{\infty}\sum_{n=1}^N\left\{a_m\alpha_n\cos(m\omega_ct)\cos(n\omega_st)+a_m\beta_n\cos(m\omega_ct)\sin(n\omega_st)\right\}}$
$\displaystyle{ =BS(t)+}$
$\displaystyle{ B\sum_{m=1}^{\infty}\sum_{n=1}^N\left\{\frac{a_m\alpha_n}{2}(\cos(m\omega_ct+n\omega_st)+\cos(m\omega_ct-n\omega_st))+\frac{a_m\beta_n}{2}(\sin(m\omega_ct+n\omega_st)-\sin(m\omega_ct-n\omega_st))\right\}}$
さて，ここで上式の第一項に関しては，Bの値が既知なので第二項さえ取り除くことができれば，PAM波をもとの信号波に完全に復元できたと言える．ではどのような手法で第二項を取り除くか？またその手法を用いることのできる条件は何か？
　ここで用いるのが低域通過フィルタである．これは，合成波のなかで，ある周波数以上の周波数成分を取り除く電子回路である（このページの後ろに簡単な説明を載せときます）．ではここではいくつ以上の周波数成分を取り除けばよいかというと，もちろん $\omega\gt\omega_{max}$ となる成分を取り除けば良い．なぜなら， $\omega_{max}$ 以下の周波数成分を取り除くと上式の第一項の成分も取り除かれて，PAM波から信号波を完全に取り除くことができなくなるからである．更にいうと， $\omega\gt\omega_{max}$ を取り除く低域通過フィルタを作用させたときに第二項が完全に除去される必要がある．つまり，第二項の成分のうち，周波成の最小値が $\omega\gt\omega_{max}$ を満たしていれば良い．では，上式第二項の最小周波数はいくつなのか．整数 $m$ の範囲は $1\le m \le \infty$ で， $n$ は $1\le n \le N$ なので，最小周波数は $m\omega_c-n\omega_s$ に $m=1, n=N$ を代入した $\omega_c-N\omega_s=\omega_c-\omega_{max}$ となる．したがって，低域通過フィルタを用いてPAM波から信号波を完全に復調できる条件は，低域通過フィルタ $\omega\gt\omega_{max}$ で上式第二項の最小周波数成分 $\omega_c-\omega_{max}$ を除去できることであり，数式にすると，
$\displaystyle{\omega_c-\omega_{max}\gt\omega_{max}\Longleftrightarrow\omega_c\gt2\omega_{max}}$
つまり，搬送波として用いる矩形周波数は，信号波に含まれる周波数の最大値の二倍より大きくなくてはならない．また， $\omega_c\gt2\omega_{max}$ を満たすような周波数 $\omega_c$ をナイキスト周波数とよんだりする．

エイリアシング

　写真などでよくわからない周期的な模様ができたりする．これをエイリアシングとか折返し雑音とか呼んだりする．これは上の式で矩形波の周波数がナイキスト周波数より小さかったとき，第二項が十分に取り除かれずに残ってしまったときに生じるものであるらしい．この場合上の式のt が写真内の距離xに対応して，周波数はピクセル数に対応するのだと思う．

矩形波のフーリエ展開

　矩形波は周期Tとして，区間 $-T/2\le t\le T/2$ において，
$\displaystyle{ p(t) = \left\{ \begin{array}{cc} A,& -\tau/2\le t\le\tau/2 \\ 0, & otherwise \end{array} \right. }$
で表される．周期関数をフーリエ展開する場合，
$\displaystyle{ f(t)=a_0+\sum_{n=1}^\infty(a_n\cos(2n\pi t/T)+b_n\sin(2n\pi t/T)), }$
$\displaystyle{ a_n=\frac{2}{T}\int_{-T/2}^{T/2}f(t)\cos(2n\pi t/T)dt, }$
$\displaystyle{ b_n=\frac{2}{T}\int_{-T/2}^{T/2}f(t)\sin(2n\pi t/T)dt, }$
である．ここで，上の式の $f(x)$ に矩形波 $p(x)$ を代入するが， $p(x)$ は偶関数なので $b_n=0$ となる．次に $a_n$ を求める．
$\displaystyle{a_n=\frac{2}{T}\int_{-T/2}^{T/2}p(t)\cos(2n\pi t/T)dt}$
$\displaystyle{=\frac{2}{T}\int_{-\tau/2}^{\tau/2}A\cos(2n\pi t/T)dt}$
$\displaystyle{=\frac{A}{n\pi}\left\{\sin(n\pi\tau/T)-sin(-n\pi\tau/T)\right\}}$
$\displaystyle{=\frac{2A}{n\pi}\sin(n\pi\tau/T)}$
$\displaystyle{=\frac{2A}{n\pi}\sin(n\pi\tau/T)\cdot\frac{\tau/T}{\tau/T}}$
$\displaystyle{=\frac{2A\tau}{T}\frac{\sin(n\pi\tau/T)}{n\pi\tau/T}}$
$\displaystyle{=\frac{2A\tau}{T}Sinc(n\pi\tau/T)}$

低域通過フィルター

低域通過フィルタはローパスフィルタ，LPSなどと呼ばれる電子回路である．WIkipediaによると，「ローパスフィルタ（英語: Low-pass filter: LPF、低域通過濾波器）とは、フィルタの一種で、なんらかの信号のうち、遮断周波数より低い周波数の成分はほとんど減衰させず、遮断周波数より高い周波数の成分を逓減させるフィルタである。ハイカットフィルタ等と呼ぶ場合もある。電気回路・電子回路では、フィルタ回路の一種である。」らしい．

2021-04-27

Plotlyで線グラフをたくさん描く

Python Programming Plotly

Plotlyとは？
線グラフ
Plotlyのグラフをファイルに保存する
ソースコード
- いろいろな線グラフ

Plotlyとは？

Pythonのグラフ描画ライブラリーで描画したグラフの拡大縮小，軸のスケールの変更などなどができる．
Plotlyのグラフ描画モジュールにはexpressとgraph_objectsがあるらしいが，基本はgraph_objectsでexpressモジュールはサクッとグラフを書くためのモジュールで内部ではgraph_objectsの関数を呼んでいるらしい．（公式ドキュメントPlotly Express | Python | Plotlyの雑和訳）

f:id:Ibetoge:20210427124756g:plain — Plotlyの動くグラフ

今回はgraph_objectsでいい感じのグラフを色々書いてみる

expressでの描写はこちら→
ibetoge.hatenablog.com

線グラフ

基本的なグラフを描写してみる．とりあえず使うモジュールをインポートする．

import numpy as np
import plotly.graph_objects as go

ひとまず $0\le t\le2\pi$ で $\sin{4t},e^t$ を描く．

# データを用意
t = np.linspace(0, 2*np.pi, 10)
sint = np.sin(4*t)
ex = np.exp(t/np.pi)

# グラフのインスタンス作成
fig = go.Figure(data=[
    go.Scatter(x=t, y=sint, name=r'$\sin t$'),
    go.Scatter(x=t, y=ex, name=r'$e^t$')
])

# 動かせるグラフをブラウザーで表示
fig.show()

基本的にはこれでグラフが描ける．

f:id:Ibetoge:20210427170338p:plain — 基本的な線グラフ

あとは，軸に名前をつけるのと，作成したグラフのインスタンスに新しいグラフを追加してみたい．

# グラフ追加
fig.add_trace(go.Scatter(x=t, y=t**2/4, name=r'$t^2/4$', mode='markers'))  # ドットのみ
fig.add_trace(go.Scatter(x=t, y=sint*ex, name=r'$e^t\sin t$', mode='lines'))  # 線のみ
fig.add_trace(go.Scatter(x=t, y=1/t, name=r'$1/t$', mode='lines+markers'))  # ドットと線

# 軸の名前を追加
fig.update_xaxes(title="t")
fig.update_yaxes(title=r"$f(t)$")

こんな感じで追加できる．このとき，modeで線グラフの種類を選べる．modeを指定しないとPlotlyがいい感じに決めてくれる．グラフを追加できるので，まず空のインスタンスを作ってから追加していってもOK．つまり，

# 空のインスタンス作成
fig = go.Figure()

# グラフ追加
fig.add_trace(go.Scatter(x=t, y=t**2/4, name=r'$t^2/4$')) 
fig.add_trace(go.Scatter(x=t, y=sint*ex, name=r'$e^t\sin t$'))
fig.add_trace(go.Scatter(x=t, y=1/t, name=r'$1/t$' )) 

fig.show()

みたいな感じで書いてもいい．

f:id:Ibetoge:20210427171237p:plain — いろいろな線グラフ

f:id:Ibetoge:20210427171830p:plain — いろいろな線グラフ（横軸のメッシュ細かめ）

Plotlyのグラフをファイルに保存する

fig.show()

を用いると，実行した際に動かせるグラフがブラウザーに表示されるが，グラフをhtmlファイルとして保存したい場合は，

#htmlで保存
fig.write_html("figure.html")

でOK．画像ファイルとして保存するには，kaleidoパッケージが必要なので，次のどちらかのコマンドを叩いてkaleidoをインストールする．(kaleidoの代わりにorcaパッケージでも良いらしいが非推奨)
pipの場合

 pip install -U kaleido

anacondaの場合

 conda install -c conda-forge python-kaleido

kaleidoをインストールしたら，

#pngで保存
fig.write_image("figure.png")
#jpegで保存
fig.write_image("figure.jpeg")

と書けばすれば良し．こうすれば保存したhtmlファイルをブラウザで開けばいつでも動かせるグラフが見られる．もちろん画像ファイルとして保存した場合は動かせない．

棒グラフとは円グラフとかも書こうと思ったけど疲れたからここまで

ソースコード

いろいろな線グラフ

import numpy as np
import plotly.graph_objects as go

# データを用意
t = np.linspace(0, 2*np.pi, 100)
sint = np.sin(4*t)
ex = np.exp(t/np.pi)

# グラフの空インスタンス作成
fig = go.Figure()

# グラフ追加
fig.add_trace(go.Scatter(
    x=t, y=t**2/4, name=r'$t^2/4$', mode='markers'))  # ドットのみ
fig.add_trace(go.Scatter(x=t, y=sint*ex,
                         name=r'$e^t\sin t$', mode='lines'))  # 線のみ
fig.add_trace(go.Scatter(x=t, y=1/t, name=r'$1/t$',
                         mode='lines+markers'))  # ドットと線
# 軸の名前を追加
fig.update_xaxes(title="t")
fig.update_yaxes(title=r"$f(t)$")

# 動かせるグラフをブラウザーで表示
fig.show()
#htmlで保存
fig.write_html("figure.html")
#pngで保存
fig.write_image("figure.png")

2021-04-27

Plotlyで簡単なグラフを描いてみる（Python）

Python Plotly Programming

Plotlyとは？
plotly.expressでとりあえず描く
Plotlyのグラフをファイルに保存する
最後にソースコード

Plotlyとは？

今回はexpressでサインカーブを書いてみる．
graph_objectsでの描写はこちら→
ibetoge.hatenablog.com

plotly.expressでとりあえず描く

環境　macOS，python 3.7.6，Plotly 4.14.3
Plotlyがない人はターミナルで下のコマンドを叩いてインストール．

$ pip install plotly

モジュールをインポートする．

import numpy as np
import plotly.express as px

まずはデータを用意する．今回は横軸 $t$ の範囲を $0\le t\le 4\pi$ としてサインカーブ $f(t)=\sin t$ を書くのでそれを用意する．

t = np.linspace(0, 4*np.pi, 100)
sint = np.sin(t)

次にplotly.expressモジュールで描画する．

#グラフのインスタンスを作成
fig = px.line(x=t, y=sint)
#グラフをブラウザーで表示
fig.show()

これで動かせるグラフが描画できる．
グラフの軸に名前を設定するには，インスタンス作成時に

fig = px.line(x=t, y=sint, labels={'x': r"$t$", 'y': r"$\sin t$"})

とすればOK．ここでr"$\sin t$"とかのr"$$"の部分ですが，Plotlyで文字列にLatex数式を使うときはr"$$"を使用します．
もっとしっかり目盛りなども設定したい人向け→

f:id:Ibetoge:20210427144425p:plain — Expressで書いたサインカーブ(これはPNGなので動きません)

Plotlyのグラフをファイルに保存する

fig.show()

を用いると，実行した際に動かせるグラフがブラウザーに表示されるが，グラフをhtmlファイルとして保存したい場合は，

#htmlで保存
fig.write_html("figure.html")

 $ pip install -U kaleido

anacondaの場合

 $ conda install -c conda-forge python-kaleido

kaleidoをインストールしたら，

#pngで保存
fig.write_image("figure.png")
#jpegで保存
fig.write_image("figure.jpeg")

最後にソースコード

import plotly.express as px
import numpy as np

# 横軸は0<t<4pi，グラフはsin(t)
t = np.linspace(0, 4*np.pi, 100)
sint = np.sin(t)

# グラフのインスタンス作成（軸ラベル付き）
fig = px.line(x=t, y=sint, labels={'x': r"$t$", 'y': r"$\sin t$"})
fig.update_xaxes(title="t")
fig.update_yaxes(title=r"$\sin t$")

# グラフをブラウザーに表示
fig.show()
# htmlで保存
fig.write_html("figure.html")

2020-12-19

極座標形式の2次元データをMatplotlibで描写する

Python Matplotlib Programming

極座標からデカルト座標に直して描写する．

f:id:Ibetoge:20201219144408p:plain — 右の図を作ろう

f:id:Ibetoge:20201219144020p:plain — 右の図を作ろう

数値解析を極座標で行った後，2次元の極座標形式データを描写したい時がありました．
例えば，
${\displaystyle \begin{eqnarray} \left\{ \begin{array}{l} r=x^2+y^2 \\ \theta=\arccos\left(\frac{x}{\sqrt{x^2+y^2}}\right) \end{array} \right. \end{eqnarray} }$
のような座標でのデータを描写したい時があります．つまり

double r[imax];
double theta[jmax];
double T[jmax][imax];

のような配列データを格納したcsvファイルやらが手元にあります．それらをpythonを用いてmatplotlibで描写しましょう．
( $r,~\theta$ 方向にそれぞれ ${\rm imax,~jmax}$ 分割してるような空間を配列r,thetaによって表現していて，
そのグリッド点で物理量T[jmax][imax]を定義しています．)
上の配列はそれぞれr.csv，theta.csv，T.csvに入っていることにしておきます．
今回は
${\displaystyle \begin{equation} T=\mathrm{e}^{-r^2} \end{equation} }$
として，Tを描写します．

流れとしては，

極座標配列をcsvから読み込む
x-y平面のデータに変換(スプライン補完)
描写

という感じです．
まずは必要なモジュールたちをインポートします．

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

scipyはスプライン補完に使い，mpl_toolkits.mplot3dと matplotlib.pyplotは描写，numpyは色々なことに使います．

次にnumpyでcsvファイルを読み込んで配列に格納しましょう．

r = np.loadtxt("r.csv", delimiter=',')
theta = np.loadtxt("theta.csv", delimiter=',')
T_polar = np.loadtxt("T.csv", delimiter=',')

描写したいx-y平面とx-y平面でのTを定義します．今回は配列r,thetaで表せる部分をすべて含めて，グリッド幅は配列rと同じにしておきます．

size_r = r.shape[0]
x = np.linspace(-r[-1], r[-1], 2*size_r)
y = np.linspace(-r[-1], r[-1], 2*size_r)
T_cartesian = np.empty((2*size_r, 2*size_r))

スプライン補完でx-y平面でのTを計算

spl_T = interpolate.interp2d(theta, r, T_polar, kind='cubic')
for ix in range(0, 2*size_r):
    for iy in range(0, 2*size_r):
        r_tmp = (x[ix]**2+y[iy]**2)**0.5
        costh = x[ix]/(x[ix]**2+y[iy]**2)**0.5
        th_tmp = np.arccos(costh)
        T_cartesian[ix, iy] = spl_T(th_tmp, r_tmp)

最後に描写

x, y = np.meshgrid(x, y, indexing='ij')
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.contourf(x, y, T_cartesian)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

これで完成．

割と一般化できるやり方なので，これで任意の座標系から任意の座標系へ変換して描写できるようになります．

ちなみに今回はカラーマップの等値線図を書きましたが

ax = fig.add_subplot(1, 1, 1)
ax.contourf(x, y, T_cartesian)

の部分を変えるといろいろな図にできます．

グラフの種類

等値線図

ax = fig.add_subplot(1, 1, 1)
ax.contour(x, y, T_cartesian)

f:id:Ibetoge:20201219150159p:plain — 等値線図

3次元ワイヤーフレーム

ax = fig.add_subplot(1, 1, 1, projection="3d")
ax.plot_wireframe(x, y, T_cartesian)

f:id:Ibetoge:20201219150408p:plain — ワイヤー

３次元曲面

ax = fig.add_subplot(1, 1, 1, projection="3d")
ax.plot_surface(x, y, T_cartesian)

f:id:Ibetoge:20201219150501p:plain — 三次元曲面

コード

最後にcsv作るように書いたコードとプロットするのに使ったコードを載せときます．

CSV作るよう

#include <cmath>
#include <fstream>
#include <iostream>

using namespace std;

int main() {
  int imax = 200;
  int jmax = 128;
  double r[imax], theta[jmax], T[jmax][imax];
  double pi;
  pi = acos(-1.0);

  for (int i = 0; i < imax; i++) {
    r[i] = (i + 0.50) / imax;
  }
  for (int j = 0; j < jmax; j++) {
    theta[j] = (j + 0.50) / jmax * 2.0 * pi;
  }

  for (int i = 0; i < imax; i++) {
    for (int j = 0; j < jmax; j++) {
      T[j][i] = exp(-1.0 * pow(r[i], 2));
    }
  }
  ofstream ofsr("r.csv");
  for (int i = 0; i < imax - 1; i++) {
    ofsr << r[i] << ",";
  }
  ofsr << r[imax - 1] << endl;

  ofstream ofst("theta.csv");
  for (int j = 0; j < jmax - 1; j++) {
    ofst << theta[j] << ",";
  }
  ofst << theta[jmax - 1] << endl;

  ofstream ofsv("T.csv");
  for (int i = 0; i < imax; i++) {
    for (int j = 0; j < jmax; j++) {
      if (j == jmax - 1) {
        ofsv << T[j][i] << endl;
      } else {
        ofsv << T[j][i] << ",";
      }
    }
  }
}

プロットする

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

# csvを読み込む
r = np.loadtxt("r.csv", delimiter=',')
theta = np.loadtxt("theta.csv", delimiter=',')
T_polar = np.loadtxt("T.csv", delimiter=',')

# x-y平面を定義してスプライン補完
size_r = r.shape[0]
x = np.linspace(-r[-1], r[-1], 2*size_r)
y = np.linspace(-r[-1], r[-1], 2*size_r)
T_cartesian = np.empty((2*size_r, 2*size_r))

spl_T = interpolate.interp2d(theta, r, T_polar, kind='cubic')
for ix in range(0, 2*size_r):
    for iy in range(0, 2*size_r):
        r_tmp = (x[ix]**2+y[iy]**2)**0.5
        costh = x[ix]/(x[ix]**2+y[iy]**2)**0.5
        th_tmp = np.arccos(costh)
        T_cartesian[ix, iy] = spl_T(th_tmp, r_tmp)

# matplotlibで描写する
x, y = np.meshgrid(x, y, indexing='ij')
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.contourf(x, y, T_cartesian)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('cartesian.png')

# 極座標での描写
theta, r = np.meshgrid(theta, r, indexing='ij')
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.contourf(theta, r, T_polar)
plt.xlabel(r'$\theta$')
plt.ylabel('r')
plt.savefig('polar.png')

My Memo

私のメモ

Optunaの最適化結果をグラフで可視化して保存する

必要なライブラリーのインストール

Optunaで最適化

結果の可視化

今回のコード全体

線形計画問題を単体法（シンプレックス法）によって解くコードを Pythonで実装してみた

線形計画問題

単体法

単体法の実装

おまけ

線形計画問題の変形

サンプリング定理の証明

サンプリング定理の証明

証明

エイリアシング

矩形波のフーリエ展開

低域通過フィルター

Plotlyで線グラフをたくさん描く

Plotlyとは？

線グラフ

Plotlyのグラフをファイルに保存する

ソースコード

いろいろな線グラフ

Plotlyで簡単なグラフを描いてみる（Python）

Plotlyとは？

plotly.expressでとりあえず描く

Plotlyのグラフをファイルに保存する

最後にソースコード

極座標形式の2次元データをMatplotlibで描写する

グラフの種類

等値線図

3次元ワイヤーフレーム

３次元曲面

コード

CSV作るよう

プロットする