TensorFlow | ADDA > 線形方程式を解き直してみた > 計算効率の改善は大きくない

```txt:動作環境
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.1.0
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
```

学習コードv0.1 http://qiita.com/7of9/items/5819f36e78cc4290614e

http://qiita.com/7of9/items/49d073f964e5689d9b96
の続き。

### 概要

This article is related to [ADDA (light scattering simulator based on the discrete dipole approximation)](https://github.com/adda-team/adda/).

- TFRecordsを読込んで学習する
- input: 5 nodes
- output: 6 nodes
- サンプル数: 223,872
- 学習データ: ADDAにより計算した値
    - #input
    - x,y,z: dipole position
    - refractive index: real and imaginary part
    - #output
    - initial values for linear equation solution for (x,y,z),(real,imaginary)

5次元補間を学習させようとしている。

上記のパラメータのうち、refractive index は以下の値を用いた。

```
real part [1.45, 1.5, 1.33, 1.4]
imaginary part [0.001, 0.050000001, 0.029999999, 0.1, 0.0099999998, 9.9999997e-05]
```

ネットワークはhidden layerが30x100x100の３層とした。

## 今回

前回: http://qiita.com/7of9/items/49d073f964e5689d9b96
にて学習したネットワークを用いて、実際に線形方程式の初期値を計算してみる。
その初期値を用いて、線形方程式を解き直す。

### codeの実行引数

refractive indexの実部、虚部、dipoleの分割通の指定。

```
./adda -m 1.45 0.0001 -grid 26 
```

上記に加えて、線形方程式の初期値のファイルの読込みを加えた例。

```
./adda -m 1.45 0.0001 -grid 26 -init_field read REPLACE_NN_170729/IntField-Y.out_170729
```

### コード

refractive indexは1.45 + 0.0001i固定としている。

LN-IntField-Y.inを読込んで、EXR, EXI, EYR, EYI, EZR, EZIを学習結果に置き換えたファイルIntField-Y.out_170729を出力する。

ネットワークファイルはlearning rate=0.0001, batch size=2, step=3,000,000のファイルを用いた (プリフィックス: 170728_t2311)。

```replace_by_NN_170729.py
import numpy as np
import tensorflow as tf
import tensorflow.contrib.slim as slim
import sys
import os
import math

"""
v0.4 Aug. 05, 2017
  - remove os.remove() and use open(OUT_FILE, 'wb+')
v0.3 Jul. 29, 2017
  - replace [exr, exi, eyr, eyi, ezr, ezi]
  - add convolution calculations
    + calc_sigmoid()
    + add calc_conv()
v0.2 Jul. 29, 2017
  -  add calc_e2()
v0.1 Jul. 29, 2017
  - read and output IntField-Y file without modification
"""

# on
#   Ubuntu 16.04 LTS
#   TensorFlow v1.1
#   Python 3.5.2
#   IPython 6.0.0 -- An enhanced Interactive Python.


def calc_sigmoid(x):
    return 1.0 / (1.0 + math.exp(-x))


def calc_conv(src, weight, bias, applyActFnc):
    wgt = weight.shape
    conv = [0.0] * bias.size

    # weight
    for idx2 in range(wgt[1]):
        tmp_vec = weight[:, idx2] * src[:]
        conv[idx2] = tmp_vec.sum()

    # bias
    for idx2 in range(wgt[1]):
        conv[idx2] = conv[idx2] + bias[idx2]
    # activation function
    if applyActFnc:
        for idx2 in range(wgt[1]):
            conv[idx2] = calc_sigmoid(conv[idx2])
    return conv  # return list


def calc_e2(xr, xi, yr, yi, zr, zi):
    xx = xr*xr + xi*xi
    yy = yr*yr + yi*yi
    zz = zr*zr + zi*zi
    return (xx + yy + zz)

INP_FILE = 'LN-IntField-Y.in'
OUT_FILE = 'IntField-Y.out_170729'
NETWORK_FILE = 'LN-model_variables_170722.npy'

HEADER_TXT = 'x y z |E|^2 Ex.r Ex.i Ey.r Ey.i Ez.r Ez.i'

dat = np.genfromtxt(INP_FILE, delimiter=' ',
                    skip_header=1)  # skip header line

model_var = np.load(NETWORK_FILE)

# TODO: 0m > runtime parameter for real(m) and imag(m)
amr = 1.5
ami = 0.1

with open(OUT_FILE, 'wb+') as fhndl:
    # line 1
    hdr = np.array(HEADER_TXT).reshape(1,)
    np.savetxt(fhndl, hdr, delimiter=' ', fmt='%s')
    # line 2 ..
    for aline in dat:
        ax, ay, az = aline[0:3]
        # e2 = aline[3]
        exr, exi = aline[4:6]
        eyr, eyi = aline[6:8]
        ezr, ezi = aline[8:10]

        # input layer (5 nodes)
        inlist = (ax, ay, az, amr, ami)
        # hidden layer 1
        cnv = calc_conv(inlist, model_var[0], model_var[1], applyActFnc=True)
        # hidden layer 2
        cnv = calc_conv(cnv, model_var[2], model_var[3], applyActFnc=True)
        # hidden layer 3
        cnv = calc_conv(cnv, model_var[4], model_var[5], applyActFnc=True)
        # output layer
        cnv = calc_conv(cnv, model_var[6], model_var[7], applyActFnc=False)

        # print(exr, exi, eyr, eyi, ezr, ezi)
        # print(cnv)
        # sys.exit()
        exr, exi = cnv[0:2]
        eyr, eyi = cnv[2:4]
        ezr, ezi = cnv[4:6]
        e2 = calc_e2(exr, exi, eyr, eyi, ezr, ezi)

        res = ax, ay, az, e2, exr, exi, eyr, eyi, ezr, ezi
        np.savetxt(fhndl, np.c_[res], delimiter=' ', fmt='%.10f')

```

### 初期値ファイル

学習元ファイル。

```
$ head LN-IntField-Y.in 
x y z |E|^2 Ex.r Ex.i Ey.r Ey.i Ez.r Ez.i
-0.2166615618 -1.516630933 -5.416539045 0.5516786267 -0.02931061187 0.02652841941 0.1694307002 0.7044162375 -0.09245370677 -0.1290700275
0.2166615618 -1.516630933 -5.416539045 0.5516786268 0.0293106112 -0.02652841834 0.1694306999 0.7044162377 -0.09245370673 -0.1290700275
-1.083307809 -1.083307809 -5.416539045 0.5009708155 -0.1210777782 0.1210667177 0.3201499555 0.5959123603 -0.06902233265 -0.09634427508
-0.6499846854 -1.083307809 -5.416539045 0.5264383492 -0.11184864 0.05898299433 0.1712903828 0.6818062435 -0.06862119958 -0.1074254128
-0.2166615618 -1.083307809 -5.416539045 0.7409394521 -0.04313270512 0.006539303629 0.1179657327 0.8426239137 -0.06420657934 -0.1047988574
0.2166615618 -1.083307809 -5.416539045 0.7409394524 0.04313270474 -0.00653930278 0.1179657325 0.8426239139 -0.06420657931 -0.1047988574
0.6499846854 -1.083307809 -5.416539045 0.5264383496 0.1118486397 -0.0589829935 0.1712903821 0.6818062441 -0.0686211995 -0.1074254129
1.083307809 -1.083307809 -5.416539045 0.5009708155 0.1210777777 -0.121066717 0.3201499543 0.5959123613 -0.06902233252 -0.09634427512
-1.083307809 -0.6499846854 -5.416539045 0.6583579786 -0.0961687486 0.04927804902 0.3522909973 0.7189310502 -0.04061949379 -0.0637218901
```

学習結果ファイル。

```
$ head IntField-Y.out_170729 
x y z |E|^2 Ex.r Ex.i Ey.r Ey.i Ez.r Ez.i
-0.2166615618 -1.5166309330 -5.4165390450 0.5675500801 -0.0193769136 0.0257250985 0.1257248657 0.7376280124 -0.0201897607 -0.0787615494
0.2166615618 -1.5166309330 -5.4165390450 0.5649162936 0.0557718589 -0.0353116427 0.1250714218 0.7332224673 -0.0245850118 -0.0818316035
-1.0833078090 -1.0833078090 -5.4165390450 0.5099061225 -0.1043974991 0.1243685041 0.2900170477 0.6306092703 -0.0021055274 -0.0419211991
-0.6499846854 -1.0833078090 -5.4165390450 0.6183017704 -0.0621416877 0.0667333275 0.1454713949 0.7648010827 -0.0139817942 -0.0608993483
-0.2166615618 -1.0833078090 -5.4165390450 0.7161159885 -0.0154258305 0.0209483584 0.0586382421 0.8403647786 -0.0219150020 -0.0728527733
0.2166615618 -1.0833078090 -5.4165390450 0.7160066138 0.0388290424 -0.0211947038 0.0539279894 0.8395532457 -0.0242294204 -0.0755298526
0.6499846854 -1.0833078090 -5.4165390450 0.6210297700 0.0995246920 -0.0723342541 0.1434001493 0.7617724520 -0.0185808021 -0.0684561446
1.0833078090 -1.0833078090 -5.4165390450 0.5302857067 0.1531068315 -0.1190136232 0.3017937814 0.6316543950 -0.0074256838 -0.0505752936
-1.0833078090 -0.6499846854 -5.4165390450 0.5640130338 -0.0614287574 0.0871006316 0.2572461696 0.6969037454 0.0020706913 -0.0282544088
```

プラスマイナスが変わっている点などもあり思ったよりもよくないかも。

## 結果

### オリジナルの計算

```
$./adda -m 1.45 0.0001 -grid 26
all data is saved in 'run389_sphere_g26_m1.45'
box dimensions: 26x26x26
lambda: 6.283185307   Dipoles/lambda: 14.5
Required relative residual norm: 1e-05
Total number of occupied dipoles: 9328
Memory usage for MatVec matrices: 5.3 MB
Calculating Green's function (Dmatrix)
Fourier transform of Dmatrix......
Initializing FFTW3
Total memory usage: 9.0 MB

here we go, calc Y

CoupleConstant:0.005318998739+1.990081942e-05i
x_0 = 0
RE_000 = 1.0000000000E+00
RE_001 = 9.9949888119E-01  + 
RE_002 = 6.6288438828E-01  + 
RE_003 = 5.7899435298E-01  + 
RE_004 = 3.8132777430E-01  + 
RE_005 = 3.7959771539E-01  + 
RE_006 = 3.7506505878E-01  + 
RE_007 = 3.1110177742E-01  + 
RE_008 = 2.7547296224E-01  + 
RE_009 = 2.7167173776E-01  + 
RE_010 = 2.5735390676E-01  + 
RE_011 = 2.4526366657E-01  + 
RE_012 = 2.4474777355E-01  + 
RE_013 = 1.6142200852E-01  + 
RE_014 = 1.2335888902E-01  + 
RE_015 = 1.2844200469E-01  - 
RE_016 = 9.7791700094E-02  + 
RE_017 = 6.1123861971E-02  + 
RE_018 = 4.9122562721E-02  + 
RE_019 = 5.3089721290E-02  - 
RE_020 = 5.0091507229E-02  -+
RE_021 = 1.5335196047E-02  + 
RE_022 = 1.2012757440E-02  + 
RE_023 = 9.6921655623E-03  + 
RE_024 = 1.0083335345E-02  - 
RE_025 = 4.6322051344E-03  + 
RE_026 = 1.9955243149E-03  + 
RE_027 = 2.1216783263E-03  - 
RE_028 = 1.5285732569E-03  + 
RE_029 = 1.2434615253E-03  + 
RE_030 = 1.4109246275E-03  - 
RE_031 = 4.9466095450E-04  + 
RE_032 = 5.9617103286E-04  - 
RE_033 = 7.0183505780E-04  - 
RE_034 = 6.7331242455E-04  -+
RE_035 = 5.2314565357E-04  -+
RE_036 = 2.7672873054E-04  + 
RE_037 = 1.9327532420E-04  + 
RE_038 = 1.4347435105E-04  + 
RE_039 = 1.1264635015E-04  + 
RE_040 = 6.6036490566E-05  + 
RE_041 = 4.3628388762E-05  + 
RE_042 = 3.0040318652E-05  + 
RE_043 = 2.3526363249E-05  + 
RE_044 = 1.6655541961E-05  + 
RE_045 = 1.2351988965E-05  + 
RE_046 = 1.2181461237E-05  + 
RE_047 = 4.4332533947E-06  + 
Cext	= 357.7373256
Qext	= 3.55623704
Cabs	= 0.2646844797
Qabs	= 0.002631206428
```

### 学習結果を用いた計算

```
$./adda -m 1.45 0.0001 -grid 26 -init_field read REPLACE_NN_170729/IntField-Y.out_170729 
all data is saved in 'run390_sphere_g26_m1.45'
box dimensions: 26x26x26
lambda: 6.283185307   Dipoles/lambda: 14.5
Required relative residual norm: 1e-05
Total number of occupied dipoles: 9328
Memory usage for MatVec matrices: 5.3 MB
Calculating Green's function (Dmatrix)
Fourier transform of Dmatrix......
Initializing FFTW3
Total memory usage: 9.0 MB

here we go, calc Y

CoupleConstant:0.005318998739+1.990081942e-05i
x_0 = from file REPLACE_NN_170729/IntField-Y.out_170729
RE_000 = 1.1679879124E-01
RE_001 = 6.4230412955E-02  + 
RE_002 = 7.1799992180E-02  - 
RE_003 = 4.7224492616E-02  + 
RE_004 = 2.6556103969E-02  + 
RE_005 = 1.7048702747E-02  + 
RE_006 = 1.3900485613E-02  + 
RE_007 = 1.0175642224E-02  + 
RE_008 = 9.7214647819E-03  + 
RE_009 = 8.5396162000E-03  + 
RE_010 = 7.5743086840E-03  + 
RE_011 = 6.9725526454E-03  + 
RE_012 = 5.8677519068E-03  + 
RE_013 = 4.8047038760E-03  + 
RE_014 = 4.0509568037E-03  + 
RE_015 = 4.1032080549E-03  - 
RE_016 = 4.1832924557E-03  - 
RE_017 = 2.5945649758E-03  + 
RE_018 = 1.6861527244E-03  + 
RE_019 = 1.2652969230E-03  + 
RE_020 = 7.2534589225E-04  + 
RE_021 = 6.9663758963E-04  + 
RE_022 = 6.8015398100E-04  + 
RE_023 = 6.7451996447E-04  + 
RE_024 = 7.3766749041E-04  - 
RE_025 = 7.4189053102E-04  - 
RE_026 = 6.1413205555E-04  + 
RE_027 = 4.0672190391E-04  + 
RE_028 = 1.9807999097E-04  + 
RE_029 = 1.3867484335E-04  + 
RE_030 = 1.2718760642E-04  + 
RE_031 = 7.2156403839E-05  + 
RE_032 = 2.5132396270E-05  + 
RE_033 = 2.2675149604E-05  + 
RE_034 = 2.2036500883E-05  + 
RE_035 = 1.9574129179E-05  + 
RE_036 = 1.7800932787E-05  + 
RE_037 = 7.0289329649E-06  + 
Cext	= 357.7374614
Qext	= 3.556238389
Cabs	= 0.2646839962
Qabs	= 0.002631201621
```

### 考察

オリジナルの計算は`RE_000 = 1.0000000000E+00`の誤差から始まる。　
一方で、学習結果を用いた計算は`RE_000 = 1.1679879124E-01`の誤差から始まる。

初期誤差は確かに小さくなった。

しかしながら、オリジナルの計算は47回のiterationで収束しているのに対して、学習結果を用いた計算は37回のiterationで収束している。

思ったよりも計算効率は改善していない。

線形方程式を解き直すため、両者の結果(Cext, Qext, Cabs, Qabs)の値に差異はない。この点は予定通り。