Illustration2VecをONNX経由で使う
趣味プロジェクトでIllustration2Vecを使いたくなったのですが、これは2015年の論文なのでモデルをCaffeかChainerで使うことになっています。 github.com
残念ながらCaffeもChainerも既に開発が終了しているため、Illustration2VecのモデルをONNXという共通フォーマットに変換して今後も使えるようにしました。 利用方法だけ知りたい人は「モデルの変換」を飛ばして「使い方」を見てください。
モデルの変換
まずはオリジナルのIllustration2Vecのモデルをダウンロードします。以下を実行するとCaffeのモデルがダウンロードできます。
git clone https://github.com/rezoo/illustration2vec.git cd illustration2vec ./get_models.sh
このスクリプトでは特徴抽出モデルのprototxtがダウンロードできなかったので、Illustration2VecのInternet ArchiveからNetwork configuration file (feature vectors) illust2vec.prototxt
を追加でダウンロードしました。
必要なライブラリをインストールします。
pip install onnx coremltools onnxmltools
以下のPythonスクリプトを実行すると、タグ予測モデルのillust2vec_tag_ver200.onnx
と特徴ベクトル抽出モデルのillust2vec_ver200.onnx
が作成されます。
import os import onnx import coremltools import onnxmltools # CaffeモデルをONNX形式で読み込む関数 # https://github.com/onnx/onnx-docker/blob/master/onnx-ecosystem/converter_scripts/caffe_coreml_onnx.ipynb def caffe_to_onnx(proto_file, input_caffe_path): output_coreml_model = 'model.mlmodel' # 中間ファイル名 # 中間ファイルが既に存在したら例外を送出 if os.path.exists(output_coreml_model): raise FileExistsError('model.mlmodel already exists') # CaffeのモデルをCore MLに変換 coreml_model = coremltools.converters.caffe.convert( (input_caffe_path, proto_file)) # Core MLモデルを保存 coreml_model.save(output_coreml_model) # Core MLモデルを読み込む coreml_model = coremltools.utils.load_spec(output_coreml_model) # Core MLモデルをONNXに変換 onnx_model = onnxmltools.convert_coreml(coreml_model) # Core MLモデルを削除 os.remove(output_coreml_model) return onnx_model # タグ予測モデルの変換・保存 onnx_tag_model = caffe_to_onnx( 'illust2vec_tag.prototxt', 'illust2vec_tag_ver200.caffemodel') onnxmltools.utils.save_model(onnx_tag_model, 'illust2vec_tag_ver200.onnx') # 特徴ベクトル抽出モデルの変換・保存 onnx_model = caffe_to_onnx('illust2vec.prototxt', 'illust2vec_ver200.caffemodel') # encode1レイヤーをONNXから利用できるようにする # https://github.com/microsoft/onnxruntime/issues/2119 intermediate_tensor_name = 'encode1' intermediate_layer_value_info = onnx.helper.ValueInfoProto() intermediate_layer_value_info.name = intermediate_tensor_name onnx_model.graph.output.extend([intermediate_layer_value_info]) onnx.save(onnx_model, 'illust2vec_ver200.onnx')
使い方
上のようにしてONNX形式に変換したモデルと、それを利用するためのコードを用意しました。 github.com
ONNX形式を他のフレームワークで読み込んで実行してもいいのですが、ONNX RuntimeというMicrosoft製のパフォーマンスを重視した推論専用のライブラリがあったのでこれを使うことにしました。 UbuntuでONNX RuntimeをCPU向けにインストールするコマンドは以下の通りです。
sudo apt install libgomp1 pip install onnxruntime
例を実行する前にコードとpre-trainedモデルのダウンロードを行ってください。
git clone https://github.com/kivantium/illustration2vec.git cd illustration2vec ./get_onnx_models.sh
タグ予測
コード
import i2v from PIL import Image from pprint import pprint illust2vec = i2v.make_i2v_with_onnx( "illust2vec_tag_ver200.onnx", "tag_list.json") img = Image.open("images/miku.jpg") pprint(illust2vec.estimate_plausible_tags([img], threshold=0.5))
入力
Hatsune Miku (初音ミク), © Crypton Future Media, INC., http://piapro.net/en_for_creators.html. This image is licensed under the Creative Commons - Attribution-NonCommercial, 3.0 Unported (CC BY-NC).
出力
[{'character': [('hatsune miku', 0.9999994039535522)], 'copyright': [('vocaloid', 0.9999999403953552)], 'general': [('thighhighs', 0.9956372976303101), ('1girl', 0.9873461723327637), ('twintails', 0.9812833666801453), ('solo', 0.9632900953292847), ('aqua hair', 0.9167952537536621), ('long hair', 0.8817101716995239), ('very long hair', 0.8326565027236938), ('detached sleeves', 0.7448851466178894), ('skirt', 0.6780778169631958), ('necktie', 0.560835063457489), ('aqua eyes', 0.5527758598327637)], 'rating': [('safe', 0.9785730242729187), ('questionable', 0.02053523063659668), ('explicit', 0.0006299614906311035)]}]
Chainer版とほとんど同じ結果が出力されました。Chainerではこの処理に6秒かかっていましたが、onnx-runtimeだと2秒で実行できたのでたしかにパフォーマンスにも優れているようです(ChainerではCaffeのモデルを変換する手間が掛かっているので1枚を処理する時間で比較するのは公平ではないですが)。
特徴ベクトルの抽出
コード
import i2v from PIL import Image illust2vec = i2v.make_i2v_with_onnx("illust2vec_ver200.onnx") img = Image.open("images/miku.jpg") # extract a 4,096-dimensional feature vector result_real = illust2vec.extract_feature([img]) print("shape: {}, dtype: {}".format(result_real.shape, result_real.dtype)) print(result_real) # i2v also supports a 4,096-bit binary feature vector result_binary = illust2vec.extract_binary_feature([img]) print("shape: {}, dtype: {}".format(result_binary.shape, result_binary.dtype)) print(result_binary)
先ほどと同じ入力に対する出力
shape: (1, 4096), dtype: float32 [[ 7.474596 3.6860986 0.537967 ... -0.14563629 2.7182112 7.3140917 ]] shape: (1, 512), dtype: uint8 [[246 215 87 107 249 190 101 32 187 18 124 90 57 233 245 243 245 54 229 47 188 147 161 149 149 232 59 217 117 112 243 78 78 39 71 45 235 53 49 77 49 211 93 136 235 22 150 195 131 172 141 253 220 104 163 220 110 30 59 182 252 253 70 178 148 152 119 239 167 226 202 58 179 198 67 117 226 13 204 246 215 163 45 150 158 21 244 214 245 251 124 155 86 250 183 96 182 90 199 56 31 111 123 123 190 79 247 99 89 233 61 105 58 13 215 159 198 92 121 39 170 223 79 245 83 143 175 229 119 127 194 217 207 242 27 251 226 38 204 217 125 175 215 165 251 197 234 94 221 188 147 247 143 247 124 230 239 34 47 195 36 39 111 244 43 166 118 15 81 177 7 56 132 50 239 134 78 207 232 188 194 122 169 215 124 152 187 150 14 45 245 27 198 120 146 108 120 250 199 178 22 86 175 102 6 237 111 254 214 107 219 37 102 104 255 226 206 172 75 109 239 189 211 48 105 62 199 238 211 254 255 228 178 189 116 86 135 224 6 253 98 54 252 168 62 23 163 177 255 58 84 173 156 84 95 205 140 33 176 150 210 231 221 32 43 201 73 126 4 127 190 123 115 154 223 79 229 123 241 154 94 250 8 236 76 175 253 247 240 191 120 174 116 229 37 117 222 214 232 175 255 176 154 207 135 183 158 136 189 84 155 20 64 76 201 28 109 79 141 188 21 222 71 197 228 155 94 47 137 250 91 195 201 235 249 255 176 245 112 228 207 229 111 232 157 6 216 228 55 153 202 249 164 76 65 184 191 188 175 83 231 174 158 45 128 61 246 191 210 189 120 110 198 126 98 227 94 127 104 214 77 237 91 235 249 11 246 247 30 152 19 118 142 223 9 245 196 249 255 0 113 2 115 149 196 59 157 117 252 190 120 93 213 77 222 215 43 223 222 106 138 251 68 213 163 57 54 252 177 250 172 27 92 115 104 231 54 240 231 74 60 247 23 242 238 176 136 188 23 165 118 10 197 183 89 199 220 95 231 61 214 49 19 85 93 41 199 21 254 28 205 181 118 153 170 155 187 60 90 148 189 218 187 172 95 182 250 255 147 137 157 225 127 127 42 55 191 114 45 238 228 222 53 94 42 181 38 254 177 232 150 99]]
Chainer版と同じbinary vectorが出力されていました。
次回はこれを使ってイラストの機械学習をします。
閉殻Hartree-Fock法によるHeH+のエネルギー計算
量子化学計算を勉強するために新しい量子化学―電子構造の理論入門〈上〉を読んでいます。
サンプルコードがFortranだったので勉強がてらC言語に移植しました。
- 作者:Attila Szabo,Neil S. Ostlund
- 出版社/メーカー: 東京大学出版会
- 発売日: 1987/07/01
- メディア: 単行本
理論
本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。
ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。
難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。
ある本で30ページで説明されてて何も分からなかった内容を100ページで説明してる別の本読んだら完全に理解して何が分からなかったのか分からなくなった
— きばん (@kivantium) 2019年3月21日
実装
本のFortran実装を素直にC言語で書き直しました。
本のコードを打ち込むのが面倒だったので、検索して見つけたページにあったコードと見比べながら書きました。
最終的に得られたエネルギーはFortran実装の結果と小数第5位まで一致しました。
この差がどこで生じているのか謎ですが、実装は正しいんじゃないでしょうか……
#include <math.h> #include <stdio.h> // プロトタイプ宣言 void HFCALC(int, int, double, double, double, double, double); void INTGRL(int, int, double, double, double, double, double); void COLECT(int, int, double, double, double, double, double); void SCF(int, int, double, double, double, double, double); double calcS(double, double, double); double calcT(double, double, double); double F0(double); double calcV(double, double, double, double, double); double TWOE(double, double, double, double, double, double, double); void FORMG(void); void DIAG(double [][2], double [][2], double [][2]); void MULT(double [][2], double [][2], double [][2], int, int); void MATOUT(double [][2], int, int, int, int, char*); double S12, T11, T12, T22, V11A, V12A, V22A, V11B, V12B, V22B, V1111, V2111, V2121, V2211, V2221, V2222; // Slater型関数に対する最小2乗近似の短縮係数と軌道指数 double COEF[3][3] = {{1.0, 0.678914, 0.444635}, {0.0, 0.430129, 0.535328}, {0.0, 0.0, 0.154329}}; double EXPON[3][3] = {{0.27095, 0.151623, 0.109818}, {0.0, 0.851819, 0.405771}, {0.0, 0.0, 2.22766}}; double D1[3], A1[3], D2[3], A2[3]; double PI = 3.1415926535898; double S[2][2], X[2][2], XT[2][2], H[2][2], F[2][2], G[2][2], C[2][2], FPRIME[2][2], CPRIME[2][2], P[2][2], OLDP[2][2], TT[2][2][2][2], E[2][2]; int main(void) { int IOP = 2; // 基底関数STO-NG int N = 3; // 基底状態の平衡結合長(a.u.) double R = 1.4632; // Heに対するSTO-3G軌道指数値 // He原子に対する最良の軌道指数は1.6875 // HeH+では電子雲が縮んでいることを反映して1.24倍 double ZETA1 = 2.0925; // Hに対するSTO-3G軌道指数値 // H原子に対する最良の軌道指数は1.0 // HeH+では電子雲が縮んでいることを反映して1.24倍 double ZETA2 = 1.24; double ZA = 2.0; double ZB = 1.0; HFCALC(IOP, N, R, ZETA1, ZETA2, ZA, ZB); } /* * 1s基底関数に対応するSTO-NG最小基底を用いて2電子系2原子分子に対する * Hartree-Fock計算を行う関数 * * IOP=0 何も表示しない * IOP=1 収束の結果のみ表示 * IOP=2 全繰り返しの結果を表示 * N STO-NG を使う (N=1, 2, 3) * R 結合長 (a.u.) * ZETA1 原子Aの基底関数のSlater軌道指数値 * ZETA2 原子Bの基底関数のSlater軌道指数値 * ZA 原子Aの原子番号 * ZB 原子Bの原子番号 */ void HFCALC(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA, double ZB) { if (IOP != 0) { printf(" STO-%dG for atomic numbers %5.2lf and %5.2lf\n\n\n", N, ZA, ZB); } INTGRL(IOP, N, R, ZETA1, ZETA2, ZA, ZB); COLECT(IOP, N, R, ZETA1, ZETA2, ZA, ZB); SCF(IOP, N, R, ZETA1, ZETA2, ZA, ZB); } /* SCF計算に必要な積分を計算する関数*/ void INTGRL(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA, double ZB) { double R2 = R * R; // 軌道指数のスケーリングと短縮係数の規格化 for (int i = 0; i < N; ++i) { A1[i] = EXPON[i][N-1] * pow(ZETA1, 2.0); D1[i] = COEF[i][N-1] * pow(2.0 * A1[i] / PI, 0.75); A2[i] = EXPON[i][N-1] * pow(ZETA2, 2.0); D2[i] = COEF[i][N-1] * pow(2.0 * A2[i] / PI, 0.75); } S12 = 0.0; T11 = 0.0; T12 = 0.0; T22 = 0.0; V11A = 0.0; V12A = 0.0; V22A = 0.0; V11B = 0.0; V12B = 0.0; V22B = 0.0; V1111 = 0.0; V2111 = 0.0; V2121 = 0.0; V2211 = 0.0; V2221 = 0.0; V2222 = 0.0; // 1電子積分の計算 // V12A = 原子核Aへの各引力の非対角項 // RAP2 = 中心Aと中心Pの距離の2乗 for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) { double RAP = A2[j] * R / (A1[i] + A2[j]); double RAP2 = pow(RAP, 2.0); double RBP2 = pow(R - RAP, 2.0); S12 += calcS(A1[i], A2[j], R2) * D1[i] * D2[j]; T11 += calcT(A1[i], A1[j], 0.0) * D1[i] * D1[j]; T12 += calcT(A1[i], A2[j], R2) * D1[i] * D2[j]; T22 += calcT(A2[i], A2[j], 0.0) * D2[i] * D2[j]; V11A += calcV(A1[i], A1[j], 0.0, 0.0, ZA) * D1[i] * D1[j]; V12A += calcV(A1[i], A2[j], R2, RAP2, ZA) * D1[i] * D2[j]; V22A += calcV(A2[i], A2[j], 0.0, R2, ZA) * D2[i] * D2[j]; V11B += calcV(A1[i], A1[j], 0.0, R2, ZB) * D1[i] * D1[j]; V12B += calcV(A1[i], A2[j], R2, RBP2, ZB) * D1[i] * D2[j]; V22B += calcV(A2[i], A2[j], 0.0, 0.0, ZB) * D2[i] * D2[j]; } } // 2電子積分の計算 for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) { for (int k = 0; k < N; ++k) { for (int l = 0; l < N; ++l) { double RAP = A2[i] * R / (A2[i] + A1[j]); //double RBP = R - RAP; double RAQ = A2[k] * R / (A2[k] + A1[l]); double RBQ = R - RAQ; double RPQ = RAP - RAQ; double RAP2 = RAP * RAP; //double RBP2 = RBP * RBP; //double RAQ2 = RAQ * RAQ; double RBQ2 = RBQ * RBQ; double RPQ2 = RPQ * RPQ; V1111 += TWOE(A1[i], A1[j], A1[k], A1[l], 0.0, 0.0, 0.0) * D1[i] * D1[j] * D1[k] * D1[l]; V2111 += TWOE(A2[i], A1[j], A1[k], A1[l], R2, 0.0, RAP2) * D2[i] * D1[j] * D1[k] * D1[l]; V2121 += TWOE(A2[i], A1[j], A2[k], A1[l], R2, R2, RPQ2) * D2[i] * D1[j] * D2[k] * D1[l]; V2211 += TWOE(A2[i], A2[j], A1[k], A1[l], 0.0, 0.0, R2) * D2[i] * D2[j] * D1[k] * D1[l]; V2221 += TWOE(A2[i], A2[j], A2[k], A1[l], 0.0, R2, RBQ2) * D2[i] * D2[j] * D2[k] * D1[l]; V2222 += TWOE(A2[i], A2[j], A2[k], A2[l], 0.0, 0.0, 0.0) * D2[i] * D2[j] * D2[k] * D2[l]; } } } } if (IOP != 0) { printf(" R ZETA1 ZETA2 S12 T11\n\n"); printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", R, ZETA1, ZETA2, S12, T11); printf(" T12 T22 V11A V12A V22A\n\n"); printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", T12, T22, V11A, V12A, V22A); printf(" V11B V12B V22B V1111 V2111\n\n"); printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", V11B, V12B, V22B, V1111, V2111); printf(" V2121 V2211 V2221 V2222\n\n"); printf("%11.6lf%11.6lf%11.6lf%11.6lf\n", V2121, V2211, V2221, V2222); } } // 重なり積分 double calcS(double A, double B, double RAB2) { return pow(PI / (A + B), 1.5) * exp(-A * B * RAB2 / (A + B)); } // 運動エネルギー積分 double calcT(double A, double B, double RAB2) { return A * B / (A + B) * (3.0 - 2.0 * A * B * RAB2 / (A + B)) * pow(PI / (A + B), 1.5) * exp(-A*B*RAB2/(A+B)); } // F0関数 double F0(double ARG) { double ret; if (ARG < 1.0e-6) { ret = 1.0 - ARG / 3.0; } else { ret = sqrt(PI / ARG) * erf(sqrt(ARG)) / 2.0; } return ret; } // 核引力積分 double calcV(double A, double B, double RAB2, double RCP2, double ZC) { double v = 2.0 * PI / (A + B) * F0((A + B) * RCP2) * exp(-A * B * RAB2 / (A + B)); return -v * ZC; } // 2電子積分 double TWOE(double A, double B, double C, double D, double RAB2, double RCD2, double RPQ2) { return 2.0 * pow(PI, 2.5) / ((A + B) * (C + D) * sqrt(A + B + C + D)) * F0((A + B) * (C + D) * RPQ2 / (A + B + C + D)) * exp(-A * B * RAB2 / (A + B) - C * D * RCD2 / (C + D)); } // 積分の結果から行列を求める void COLECT(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA, double ZB) { // 核ハミルトニアン cf. 式(3.153) H[0][0] = T11 + V11A + V11B; H[0][1] = T12 + V12A + V12B; H[1][0] = H[0][1]; H[1][1] = T22 + V22A + V22B; // 重なり行列 cf. 式(3.136) S[0][0] = 1.0; S[0][1] = S12; S[1][0] = S[0][1]; S[1][1] = 1.0; // 正準直交化 cf. 式(3.169), p.191 X[0][0] = 1.0 / sqrt(2.0 * (1.0 + S12)); X[1][0] = X[0][0]; X[0][1] = 1.0 / sqrt(2.0 * (1.0 - S12)); X[1][1] = -X[0][1]; // 変換行列の転置 XT[0][0] = X[0][0]; XT[0][1] = X[1][0]; XT[1][0] = X[0][1]; XT[1][1] = X[1][1]; // 2電子積分の行列 TT[0][0][0][0] = V1111; TT[1][0][0][0] = V2111; TT[0][1][0][0] = V2111; TT[0][0][1][0] = V2111; TT[0][0][0][1] = V2111; TT[1][0][1][0] = V2121; TT[0][1][1][0] = V2121; TT[1][0][0][1] = V2121; TT[0][1][0][1] = V2121; TT[1][1][0][0] = V2211; TT[0][0][1][1] = V2211; TT[1][1][1][0] = V2221; TT[1][1][0][1] = V2221; TT[1][0][1][1] = V2221; TT[0][1][1][1] = V2221; TT[1][1][1][1] = V2222; if (IOP != 0) { MATOUT(S, 2, 2, 2, 2, "S "); MATOUT(X, 2, 2, 2, 2, "X "); MATOUT(H, 2, 2, 2, 2, "H "); printf("\n"); for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { for (int k = 0; k < 2; ++k) { for (int l = 0; l < 2; ++l) { printf(" (%2d%2d%2d%2d )%10.6lf\n", i+1, j+1, k+1, l+1, TT[i][j][k][l]); } } } } } } // SCFを実行する void SCF(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA, double ZB) { // 密度行列の収束規準 double CRIT = 1.0e-4; // 繰り返しの最大回数 int MAXIT = 25; // 核ハミルトニアンを最初の推定に用いる(つまりPを零行列で初期化) for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { P[i][j] = 0.0; } } if (IOP == 2) { MATOUT(P, 2, 2, 2, 2, "P "); } int ITER = 0; double EN; for (;;) { ITER++; if (IOP == 2) { printf("\n START OF ITERATION NUMBER = %2d\n", ITER); } // Fock行列の2電子部分を計算 FORMG(); if (IOP == 2) { MATOUT(G, 2, 2, 2, 2, "G "); } // 核ハミルトニアンを足してFock行列を得る for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { F[i][j] = H[i][j] + G[i][j]; } } // 電子エネルギーを計算する EN = 0.0; for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { EN += 0.5 * P[i][j] * (H[i][j] + F[i][j]); } } if (IOP == 2) { MATOUT(F, 2, 2, 2, 2, "F "); printf("\n\n\n ELECTRONIC ENERGY = %20.12lf\n", EN); } // Gを用いてFock行列を変換 MULT(F, X, G, 2, 2); MULT(XT, G, FPRIME, 2, 2); // Fock行列の対角化 DIAG(FPRIME, CPRIME, E); // 固有ベクトルから行列Cを得る MULT(X, CPRIME, C, 2, 2); // 新しい密度行列を計算する for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { OLDP[i][j] = P[i][j]; P[i][j] = 0.0; for (int k = 0; k < 1; ++k) { // ここはk<1で本当にいいのか不明 P[i][j] += 2.0 * C[i][k] * C[j][k]; } } } if (IOP == 2) { MATOUT(FPRIME, 2, 2, 2, 2, "F' "); MATOUT(CPRIME, 2, 2, 2, 2, "C' "); MATOUT(E, 2, 2, 2, 2, "E "); MATOUT(C, 2, 2, 2, 2, "C "); MATOUT(P, 2, 2, 2, 2, "P "); } // deltaを計算する double DELTA = 0.0; for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { DELTA += pow(P[i][j] - OLDP[i][j], 2.0); } } DELTA = sqrt(DELTA / 4.0); if (IOP != 0) { printf("\n DELTA(CONVERGENCE OF DENSITY MATRIX) = %10.6lf\n", DELTA); } // 収束判定 if (DELTA < CRIT) break; if (ITER < MAXIT) continue; // 収束しなかった場合 printf(" NO CONVERGENCE IN SCF"); return; } double ENT = EN + ZA * ZB / R; if (IOP != 0) { printf("\n\n CALCULATION CONVERGED" "\n\n ELECTORONIC ENERGY = %20.12lf" "\n\n TOTAL ENERGY = %20.12f", EN, ENT); } if (IOP == 1) { MATOUT(G, 2, 2, 2, 2, "G "); MATOUT(F, 2, 2, 2, 2, "F "); MATOUT(E, 2, 2, 2, 2, "E "); MATOUT(C, 2, 2, 2, 2, "C "); MATOUT(P, 2, 2, 2, 2, "P "); } MULT(P, S, OLDP, 2, 2); if (IOP != 0) { MATOUT(OLDP, 2, 2, 2, 2, "PS "); } } void FORMG(void) { for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { G[i][j] = 0.0; for (int k = 0; k < 2; ++k) { for (int l = 0; l < 2; ++l) { G[i][j] += P[k][l] * (TT[i][j][k][l] - 0.5 * TT[i][l][k][j]); } } } } } // Fを対角化し、Cに固有ベクトル、Eに固有値を入れる void DIAG(double F[][2], double C[][2], double E[][2]) { double theta; if (fabs(F[0][0] - F[1][1]) > 1.0e-20) { theta = 0.5 * atan(2.0 * F[0][1] / (F[0][0] - F[1][1])); } else { theta = PI / 4.0; } C[0][0] = cos(theta); C[1][0] = sin(theta); C[0][1] = sin(theta); C[1][1] = -cos(theta); E[0][0] = F[0][0] * pow(cos(theta), 2.0) + F[1][1] * pow(sin(theta), 2.0) + F[0][1] * sin(2.0 * theta); E[1][1] = F[1][1] * pow(cos(theta), 2.0) + F[0][0] * pow(sin(theta), 2.0) - F[0][1] * sin(2.0 * theta); E[1][0] = 0.0; E[0][1] = 0.0; // 固有値と固有ベクトルを整列する if (E[1][1] > E[0][0]) return; double temp = E[1][1]; E[1][1] = E[0][0]; E[0][0] = temp; temp = C[0][1]; C[0][1] = C[0][0]; C[0][0] = temp; temp = C[1][1]; C[1][1] = C[1][0]; C[1][0] = temp; } // AとBの積をCに入れる void MULT(double A[][2], double B[][2], double C[][2], int IM, int M) { for (int i = 0; i < M; ++i) { for (int j = 0; j < M; ++j) { C[i][j] = 0.0; for (int k = 0; k < M; ++k) { C[i][j] += A[i][k] * B[k][j]; } } } } void MATOUT(double A[][2], int IM, int IN, int M, int N, char *LABEL) { int IHIGH = 0; for (;;) { int LOW = IHIGH + 1; IHIGH = IHIGH + 5; IHIGH = (IHIGH < N) ? IHIGH : N; printf("\n\n\n THE %s ARRAY\n", LABEL); printf(" "); for (int i = LOW; i <= IHIGH; ++i) { printf(" %3d ", i); } printf("\n"); for (int i = 1; i <= M; ++i) { printf("%10d ", i); for (int j = LOW; j <= IHIGH; ++j) { printf(" %18.10lf", A[i - 1][j - 1]); } printf("\n"); } if (N <= IHIGH) break; } }
最終的な出力は
CALCULATION CONVERGED ELECTORONIC ENERGY = -4.227529304014 TOTAL ENERGY = -2.860662163500 THE PS ARRAY 1 2 1 1.5296370381 1.1199277981 2 0.6424383868 0.4703629619
でした。参照したFortranコードの出力は
CALCULATION CONVERGED ELECTRONIC ENERGY = -0.422752913203D+01 TOTAL ENERGY = -0.286066199152D+01 THE PS ARRAY 1 2 1 0.1529637121D+01 0.1119927796D+01 2 0.6424383099D+00 0.4703628793D+00
でした。
FortranからCに移植するときのtips
詰まったことを挙げておきます。
配列の添字
Fortranでは添字が1から始まるので配列にアクセスするときに添字をうまく調整する必要があります。
多次元配列の並び順
このページに詳しく書いてありますが、C言語の多次元配列はメモリ上でa[0][0], a[0][1], a[1][0], a[1][1], a[2][0], a[2][1]
のような順番で並びますが、Fortranではa(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2)
のような順番に並びます。多次元配列の初期化を書き直すときに気をつける必要があります。
整合配列
Fortranではサブルーチンに配列を渡した配列のサイズをサブルーチン側で定義することができます。(整合配列という)
dimension a(6) a(1) = 1 a(2) = 2 a(3) = 3 a(4) = 4 a(5) = 5 a(6) = 6 call sub(a,3,2) end subroutine sub(a,n,m) dimension a(n,m) <---- a(n,*) でも同じ結果を得る do i=1,n write(6,*) (a(i,j),j=1,m) end do return end
コードはFortran プログラミングの基礎知識から引用しました。
C99で導入された可変長配列の機能を使うと同じことができるそうです。cf. 3.3 整合配列、可変長配列
今回は配列のサイズが決め打ちできたのでvoid MULT(double A[][2], double B[][2], double C[][2], int IM, int M) {
のように関数を定義しました。
TwitterにMP4動画をアップロードするにはyuv420pを使う必要がある(らしい)
以前の記事でMP4をアニメーションGIFに変換する方法を紹介しました。
kivantium.hateblo.jp
その後Twitterの機能変更でMP4の動画がそのままアップロードできるようになったのですが、アップロードに失敗することが非常に多かったです。
Twitterの公式ドキュメント(How to share and watch videos on Twitter)には
モバイルアプリではMP4とMOVの動画形式をサポートしています。
ブラウザではMP4(H264形式、AACオーディオ)をサポートしています。アップロードできる動画のサイズは最大512MBです。ただし長さは2分20秒間以下にしてください。
最小解像度: 32 x 32
最大解像度: 1920 x 1200(および1200 x 1900)
縦横比: 1:2.39~2.39:1の範囲(両方の値を含む)
最大フレームレート: 40fps
最大ビットレート: 25Mbps
としか書いていないのですが、実際にはフォーマットをyuv420pにしないとアップロードできないらしいです。[要出典]
(そのように書いてあるTwitter公式のドキュメントは見当たらないのですが、多くのサイトでそう指摘されていました。)
今まで失敗していたのはffmpegのデフォルトだとyuv444になってしまうからのようです。
H.264で作ったのにアップロードに失敗した動画をyuv420pにするのは簡単で、
ffmpeg -i input.mp4 -pix_fmt yuv420p output.mp4
とするだけです。
コーデックの変換も含めて行うには
ffmpeg -i input.webm -vcodec libx264 -pix_fmt yuv420p -strict -2 -acodec aac output.mp4
のようにすればいいようです。
ボトムアップ型自動微分の実験
Deep Learningの関係で自動微分が脚光を浴びつつあるような気がしますが、自動微分を解説したページは少なくまだまだマイナーな分野だと思います。先日ようやく「アルゴリズムの自動微分と応用」を一通り眺めたのでいろいろ実験して遊んでいます。今日は自動微分のうち、ボトムアップ型自動微分をオペレータオーバーロードを用いて実現する方法について書きます。
自動微分とは
自動微分は、アルゴリズムによって定義された関数からその関数の偏導関数値を計算するアルゴリズムを導出するための技術です。一般的にはのような関数が先にあって、その関数を計算するアルゴリズムやプログラムがあるというように考えますが、自動微分ではどちらかというとアルゴリズムが先にあってアルゴリズムが表現する関数が生まれるというような考え方をします。
プログラムで微分を扱う上でよく知られている技術には数値微分と数式微分があります。
数値微分
数値微分は、小さな値を用いてやを計算することで微分の近似値を得る技術です。数値微分はあくまで微分の近似値を計算する技術ですが、自動微分ではアルゴリズムで定義された関数から解析的に得た偏導関数の値を計算します。また、数値微分では摂動幅を設定する必要がありますが、が大きすぎると近似精度が悪くなり、小さすぎると丸め誤差の影響で精度が悪くなるという問題があります。
この問題を理解するために簡単なプログラムを書きました。におけるの微分をを計算することで求めます。解析的にはが解なので、この2倍を出力するようにしてよく知っているの値にどれくらい近づくかを見てみます。
#include <iostream> #include <iomanip> #include <cmath> using namespace std; int main(void){ // 小数点以下15桁表示する cout << fixed << setprecision(15); // 刻み幅を小さくして数値微分 double h = 0.1; for(int i=0; i<15; ++i) { cout << 2 * (sin(M_PI/4 + h) - sin(M_PI/4))/h << " (h=" << h << ")" << endl; h /= 10; } }
結果は
1.341205945807979 (h=0.100000000000000) 1.407118983378419 (h=0.010000000000000) 1.413506219948735 (h=0.001000000000000) 1.414142849338607 (h=0.000100000000000) 1.414206491290315 (h=0.000010000000000) 1.414212855488372 (h=0.000001000000000) 1.414213490757987 (h=0.000000100000000) 1.414213568473599 (h=0.000000010000000) 1.414213635086980 (h=0.000000001000000) 1.414215411443819 (h=0.000000000100000) 1.414224293228016 (h=0.000000000010000) 1.414202088767524 (h=0.000000000001000) 1.414424133372449 (h=0.000000000000100) 1.421085471520200 (h=0.000000000000010) 1.554312234475219 (h=0.000000000000001)
のようになりました。大きすぎても小さすぎても良くなくて、h=0.000000010000000のときの結果が最も良いようです。
一般に、の計算誤差をなるべく小さくするためには関数の計算値に含まれる計算誤差をとしてとするのがよいことが知られています。しかし、を求めるのにの値が必要になってしまうのでこの式を使うのは現実的ではありません。
また、入力変数が多い関数に対しては入力変数の数だけ関数を計算する必要があり計算量が増えてしまいますが、自動微分を使うことで計算量を抑えることができます。
自動微分の種類
自動微分にはボトムアップ型(フォーワードモードとも)とトップダウン型(リバースモードとも)の2種類があります。どちらも微分の連鎖律を用いて関数を単純な演算に分解することで偏導関数値を求めますが、分解の方法が異なります。
ボトムアップ型自動微分
ボトムアップ型自動微分では一つの入力変数に対して、全ての中間変数の偏導関数値を計算していく手法です。
例としてWikipediaの自動微分のページにあるの例を考えます。この関数の計算をステップごとに中間変数に保存するとすれば、
のように5つの中間変数を用いて表すことができます。
ここで入力変数について各中間変数の微分を計算すると、
のように分解することができます。これは関数値を求める際に必要な各中間変数の計算と同時に行うことができます。
計算グラフを使って表すと以下のようになります。(Wikipediaから引用)
トップダウン型自動微分
トップダウン型自動微分では一つの出力変数について、全ての中間変数に対する偏導関数値を計算していく手法です。
先ほどと同じ例を考えます。出力変数について各中間変数に対する微分を計算すると、
のように分解することができます。
についての分解が分かりにくいですが、の変化は中間変数の変化を経由してしか起こらないことに注目すると納得できるかもしれません。計算グラフは以下の通りです。
トップダウン型自動微分では一度計算した後で各中間変数がどのように計算されたのかの履歴を逆向きに辿って偏導関数値を計算する必要があるため、計算グラフを保存するためのメモリが必要になります。
勾配を求める際には、入力変数の方が多いときはトップダウン型、出力変数の方が多い時はボトムアップ型を使うのが計算量の点で有利になります。ニューラルネットワークのバックプロパゲーションはトップダウン型の自動微分の一種です。
自動微分の実現方法
自動微分の実現方法には、ソースコードを解析して偏導関数値を求めるプログラムを作る方法とオペレータオーバーロードによる方法があります。ソースコードを解析する方法はコンパイラを作るのと同じような労力が必要になるので難しいですが、オペレータオーバーロードによる方法は比較的容易なので今回はこちらを採用します。
オペレータオーバーロードによるボトムアップ型自動微分
ボトムアップ型自動微分では関数値を計算するのと同時に偏導関数値を計算することができます。そのためdoubleを自動微分用に拡張したBUdoubleというクラスをつくって計算を行うと同時に偏導関数値を保存するようにします。それを実現したのが以下のコードです。
#include <iostream> #include <cmath> #include <iomanip> using namespace std; // ボトムアップ型微分積分用doubleクラス class BU_double { // 変数値 double val; // 偏導関数値 double d_val; public: // コンストラクタ BU_double(double v=0, double dv=0) { val = v; d_val = dv; } // 微分する入力変数として選択する関数 void select(void) { d_val = 1.0; } // 変数値を返す double get_value(void) { return val; } // 偏導関数値を返す double get_d_value(void) { return d_val; } // 各種演算子の定義 friend BU_double operator + (BU_double x, BU_double y) { return BU_double(x.val+y.val, x.d_val+y.d_val); } friend BU_double operator - (BU_double x, BU_double y) { return BU_double(x.val-y.val, x.d_val-y.d_val); } friend BU_double operator * (BU_double x, BU_double y) { return BU_double(x.val*y.val, x.d_val*y.val+x.val*y.d_val); } friend BU_double operator / (BU_double x, BU_double y) { double w = x.val/y.val; return BU_double(w, (x.d_val-w*y.d_val)/y.val); } friend BU_double operator + (BU_double x) { return BU_double(x.val, x.d_val); } friend BU_double operator - (BU_double x) { return BU_double(-x.val, -x.d_val); } friend bool operator < (BU_double x, BU_double y) { return x.val < y.val; } friend bool operator <= (BU_double x, BU_double y) { return x.val <= y.val; } friend bool operator > (BU_double x, BU_double y) { return x.val > y.val; } friend bool operator >= (BU_double x, BU_double y) { return x.val >= y.val; } // 基本関数の定義 friend BU_double sqrt(BU_double x) { return BU_double(sqrt(x.val), 0.5*x.d_val/sqrt(x.val)); } friend BU_double exp(BU_double x) { return BU_double(exp(x.val), x.d_val*exp(x.val)); } friend BU_double log(BU_double x) { return BU_double(log(x.val), x.d_val/x.val); } friend BU_double sin(BU_double x) { return BU_double(sin(x.val), cos(x.val)); } friend BU_double cos(BU_double x) { return BU_double(cos(x.val), -sin(x.val)); } // coutに出力するフォーマットの定義 friend ostream& operator<<(ostream &s, BU_double x) { return s << "BU_double("<< x.val << ", " << x.d_val << ") "; } }; int main(void){ cout << setprecision(15); BU_double x, y; // sin(x)の微分 x = M_PI / 4; x.select(); y = sin(x); cout << 2 * y.get_d_value() << endl; // 2016-07-13に修正 // y = 10*x^2 のfor文を使ったアホな定義 x = 5; y = 0; x.select(); for(int i=0; i<10; ++i) y = y + x*x; // x = 5.0のときのyとdy/dxを出力 cout << y << endl; // ニュートン法の実行 x = 2.0; // 適当な初期値 for(int i=0; i<10000; ++i) { x.select(); // y = (x-sqrt(2))*(x^3+1) y = (x-sqrt(2))*(x*x+x+1); // yの値がdoubleの丸め誤差と同程度に小さければ終了 if (y.get_value() < 1e-15) break; // ニュートン法の更新式で次のxを求める x = x - y.get_value()/y.get_d_value(); } // x(=sqrt(2))を出力 cout << x.get_value() << endl; }
出力は以下の通りです。
1.41421356237309 BU_double(250, 100) 1.4142135623731
出力の1行目はの微分ののときの値の2倍(すなわち)を表しています。最初に示した数値微分より高い精度で計算できていることが分かります。
2行目はという関数にを入れて計算したときの関数値と導関数値です。for文を使った繰り返し文として定義しているにも関わらずきちんと導関数値が求まっていることが分かります
3行目はニュートン法を用いて求めたの解を表しています。ニュートン法はの解をという更新式を用いて反復することで求めるというアルゴリズムです。導関数の計算が自動的に行われるため簡潔に記述することができます。
参考文献
- 作者: 久保田光一,伊理正夫
- 出版社/メーカー: コロナ社
- 発売日: 1998/06
- メディア: 単行本
- クリック: 1回
- この商品を含むブログを見る
- 作者: 伊理正夫,藤野和建
- 出版社/メーカー: 共立出版
- 発売日: 1985/06/03
- メディア: 単行本
- 購入: 9人 クリック: 41回
- この商品を含むブログ (10件) を見る