Go to Qiita Advent Calendar 2024 Top
LoginSignup
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@7of9

ADDA > chpointをPythonで読込む > v0.5: PEP8対応 / Cと同じ実装はできた ( print書式を除く) // v0.6: print時に少数桁数を変更

Last updated at Posted at 2016-12-23
動作環境
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 14.04 LTS desktop amd64
TensorFlow v0.11
cuDNN v5.1 for Linux
CUDA v8.0
Python 2.7.6
IPython 5.1.0 -- An enhanced Interactive Python.
gcc (Ubuntu 4.8.4-2ubuntu1~14.04.3) 4.8.4

数値光散乱シミュレーションで出力されるチェックポイントファイルの読込みをCで実装した。
http://qiita.com/7of9/items/0fbea38b48cdfe4588dd

これをmatplotlibで読むためにPython実装に変更しようとしている。

参考: C / Python > Cでfwrite()した値をPythonで読込む

v0.1-v0.3: http://qiita.com/7of9/items/68673ac3239532064c28

v0.5

  • PEP8対応
  • Cと同じ内容まで実装
    • read scalars[], xvec[], rvec[], pvec[], vectors[][]

C実装: http://qiita.com/7of9/items/0fbea38b48cdfe4588dd

read_chpoint.py
import numpy as np
import array
import sys

'''
v0.5 Dec. 23, 2016
	- read scalars[], xvec[], rvec[], pvec[], vectors[][]
    - refactor according to PEP8 (except for E501 line too long for long comment)
    - wrap_fromfile() takes [num] arg
v0.4 Dec. 21, 2016
    - read from "auxiliary" file [auxfp]
        + add read_auxiliary_info()
v0.3 Dec. 19, 2016
    - add wrap_fromfile()
    - read [inprodR],[prev_err],[resid_scale]
v0.2 Dec. 19, 2016
    - fix bug > read [local_nRows] for "size_t" type
v0.1 Dec. 18, 2016
    - read_chpoint_file() > read [ind_m],[local_nRows],[niter],[counter]
'''

def wrap_fromfile(rfp, typecode, num):
    res = array.array(typecode)
    res.fromfile(rfp, num)
    return res


def read_auxiliary_info(auxFilename):
    with open(auxFilename, "rb") as auxfp:
        sc_N = wrap_fromfile(auxfp, 'i', num=1)  # the number of scalars' sizes

        sc_sizes = array.array('i')  #
        sc_sizes.fromfile(auxfp, *sc_N)  # '*':unpack list

        vec_N = wrap_fromfile(auxfp, 'i', num=1)  # the number of vercotrs' sizes

        vec_sizes = array.array('i')  #
        vec_sizes.fromfile(auxfp, *vec_N)  # '*':unpack list

        return sc_N, sc_sizes, vec_N, vec_sizes
    # TODO: 0m > error handling when file is NULL


def read_chpoint_file(chpFilename, auxFilename):
    print(chpFilename)
    print(auxFilename)

    # TODO: 0m > error handling when file is NULL
    with open(chpFilename, "rb") as chpfp:
        ind_m = wrap_fromfile(chpfp, 'i', num=1)  # index of iterative method
        print('ind_m:', ind_m)
        # number of local rows of decomposition (only real dipoles)
        #    where L: unsigned long = size_t
        local_nRows = wrap_fromfile(chpfp, 'L', num=1)
        print('local_nRows:', local_nRows)
        niter = wrap_fromfile(chpfp, 'i', num=1)  # iteration count
        print('niter:', niter)
        # number of successive iterations without residual decrease
        counter = wrap_fromfile(chpfp, 'i', num=1)
        print('counter:', counter)
        # used as |r_0|^2 and best squared norm of residual up to some iteration
        inprodR = wrap_fromfile(chpfp, 'd', num=1)
        print('inprodR:', inprodR)
        # previous relative error; used in ProgressReport, initialized in IterativeSolver
        prev_err = wrap_fromfile(chpfp, 'd', num=1)
        print('prev_err:', prev_err)
        # scale to get square of relative error
        resid_scale = wrap_fromfile(chpfp, 'd', num=1)
        print('resid_scale:', resid_scale)  # scale to get square of relative error

        # sc_N: the number of scalars' sizes
        # sc_sizes:
        # vec_N: the number of vercotrs' sizes
        # vec_sizes:
        sc_N, sc_sizes, vec_N, vec_sizes = read_auxiliary_info(auxFilename)
        print('sc_N:', sc_N)
        print('sc_sizes', sc_sizes)
        print('vec_N:', vec_N)
        print('vec_sizes:', vec_sizes)
        # scalars[]
        scalars = []
        for idx in range(sc_N[0]):
            if sc_sizes[idx] is 8:  # double
                scl = wrap_fromfile(chpfp, 'd', num=1)
                scalars.append(scl)
            elif sc_sizes[idx] is 16:  # double complex
                scl = wrap_fromfile(chpfp, 'd', num=2)
                scalars.append(scl)
        for idx in range(sc_N[0]):
       		print('scalars:', scalars[idx])
       	# xvec[] : total electric field on the dipoles
       	# TODO: 0m > replace 'local_nRows[0]' with function?
       	# TODO: 0m > replace 'vec_N[0]' with function?
       	# TODO: 0m > define function for read [drel,dimg]??
       	xvec = [[] for idx in range(local_nRows[0])]
       	for idx in range(local_nRows[0]):
       		drel = wrap_fromfile(chpfp, 'd', 1)
       		dimg = wrap_fromfile(chpfp, 'd', 1)
       		xvec[idx] = [drel, dimg]
       	# rvec[] : current residual
       	rvec = [[] for idx in range(local_nRows[0])]
       	for idx in range(local_nRows[0]):
       		drel = wrap_fromfile(chpfp, 'd', 1)
       		dimg = wrap_fromfile(chpfp, 'd', 1)
       		rvec[idx] = [drel, dimg]       
       	# pvec[] : polarization of dipoles, also an auxiliary vector in iterative solvers
       	pvec = [[] for idx in range(local_nRows[0])]
       	for idx in range(local_nRows[0]):
       		drel = wrap_fromfile(chpfp, 'd', 1)
       		dimg = wrap_fromfile(chpfp, 'd', 1)
       		pvec[idx] = [drel, dimg]       
       	# vectors[] : ???
       	MAXNUM_VECTORS = 20
       	vectors = [[] for idx in range(MAXNUM_VECTORS)]
       	for idx_vecN in range(vec_N[0]):
       	    alist = []
            for idx_LnRows in range(local_nRows[0]):
                if vec_sizes[idx_vecN] is 16:
                    drel = wrap_fromfile(chpfp, 'd', 1)
                    dimg = wrap_fromfile(chpfp, 'd', 1)
                    alist.append([drel, dimg])
            vectors[idx_vecN] = alist
       	#
       	for idx in range(3):
       		print('xvec:',xvec[idx])
       		print('rvec:',rvec[idx])
       		print('pvec:',pvec[idx])
       		print('vectors:',vectors[idx][0])
       		print('vectors:',vectors[idx][1])
       		print('vectors:',vectors[idx][2])

argvs = sys.argv
argc = len(argvs)

print argvs
if (argc < 3):
    print("ERROR: chpoint file is not specified\r\n")
    print("   [cmd] [chpoint file] [auxiliary file]\r\n")
    sys.exit()

read_chpoint_file(chpFilename=argvs[1], auxFilename=argvs[2])

Cとの比較

C実行

ind_m:5
local_nRows:27984
niter:209
counter:139
inprodR:4.317
prev_err:0.320
resid_scale:0.020

Python実行

('ind_m:', array('i', [5]))
('local_nRows:', array('L', [27984L]))
('niter:', array('i', [209]))
('counter:', array('i', [139]))
('inprodR:', array('d', [4.317011265934948]))
('prev_err:', array('d', [0.3198525404317458]))
('resid_scale:', array('d', [0.020384616751203153]))

C実行

sc_N:8
sc_sizes[]:8
sc_sizes[]:8
sc_sizes[]:8
sc_sizes[]:8
sc_sizes[]:16
sc_sizes[]:16
sc_sizes[]:16
sc_sizes[]:16
vec_N:3
vec_sizes[]:16
vec_sizes[]:16
vec_sizes[]:16

Python実行

('sc_N:', array('i', [8]))
('sc_sizes', array('i', [8, 8, 8, 8, 16, 16, 16, 16]))
('vec_N:', array('i', [3]))
('vec_sizes:', array('i', [16, 16, 16]))

C実行

xvec:(0.000916,-0.001571), rvec:(-0.002100,0.002827), pvec:(-0.001287,-0.008525)
xvec:(0.032771,0.090957), rvec:(0.004384,-0.023016), pvec:(-0.005860,0.025872)
xvec:(-0.004613,-0.000545), rvec:(0.005385,-0.006732), pvec:(-0.000544,0.022825)
xvec:(-0.000916,0.001571), rvec:(0.002100,-0.002827), pvec:(0.001287,0.008525)
xvec:(0.032771,0.090957), rvec:(0.004384,-0.023016), pvec:(-0.005860,0.025872)
vectors 1:(-0.203491,0.033598), 2:(0.513358,-0.235685), 3:(0.548095,0.061967)
vectors 1:(0.100941,-0.212583), 2:(-0.220141,0.644602), 3:(-0.386403,0.429201)
vectors 1:(-0.000433,0.008886), 2:(0.007084,-0.024328), 3:(0.006453,-0.022867)

Python実行

('xvec:', [array('d', [0.0009161390559393775]), array('d', [-0.001570522583473328])])
('rvec:', [array('d', [-0.0020997633140212722]), array('d', [0.0028267192546817064])])
('pvec:', [array('d', [-0.0012872008298920184]), array('d', [-0.008525276560774427])])
('vectors:', [array('d', [-0.20349115145116184]), array('d', [0.03359754374072836])])
('vectors:', [array('d', [0.5133580879635192]), array('d', [-0.23568549427068541])])
('vectors:', [array('d', [0.548094579220016]), array('d', [0.061967135670311146])])
('xvec:', [array('d', [0.03277106724041005]), array('d', [0.09095666869781775])])
('rvec:', [array('d', [0.004384056712075103]), array('d', [-0.023015736282339546])])
('pvec:', [array('d', [-0.0058598173209674954]), array('d', [0.025872335594481318])])
('vectors:', [array('d', [0.1009411988714421]), array('d', [-0.21258258378933656])])
('vectors:', [array('d', [-0.22014111334002046]), array('d', [0.644601864206821])])
('vectors:', [array('d', [-0.3864030960612225]), array('d', [0.42920073139253534])])
('xvec:', [array('d', [-0.0046125547011494445]), array('d', [-0.00054452849251024])])
('rvec:', [array('d', [0.005384873980389105]), array('d', [-0.0067320912238461586])])
('pvec:', [array('d', [-0.0005437524384725202]), array('d', [0.02282484518076862])])
('vectors:', [array('d', [-0.00043270205801244216]), array('d', [0.008886034848580922])])
('vectors:', [array('d', [0.007084160219588482]), array('d', [-0.02432821023841541])])
('vectors:', [array('d', [0.00645272013550734]), array('d', [-0.022867155594089517])])

本当は違うんだ日記

  • read_chpoint_file()の中身が長すぎる問題
  • double complexの読込みはDon't repeat yourselfを無視
    • 関数にするのが良い方法かで躊躇している
  • print時にCと同じ形式で表示しようとして失敗
    • %.4fをxvec[idx]などに使おうとしたが失敗
    • １年後にはできるようになる

今日はmatplotlib表示に挑みたいがどうだろうか。

v0.6 print時に少数桁数を変更

@shiracamus さんにコメントいただいたnp.set_printoptions(precision=6)を使い、C実装の結果と同じ少数桁数にしました。

code

read_chpoint.py
import numpy as np
import array
import sys

'''
v0.6 Dec. 23, 2016
	- try to print fractional part of xvec[] in 4 digits
	   + xvec[idx] = np.array([drel, dimg]), also for rvec[], pvec[], vectors[][]
	   + add [np.set_printoptions(precision=4)]
v0.5 Dec. 23, 2016
	- read scalars[], xvec[], rvec[], pvec[], vectors[][]
    - refactor according to PEP8 (except for E501 line too long for long comment)
    - wrap_fromfile() takes [num] arg
v0.4 Dec. 21, 2016
    - read from "auxiliary" file [auxfp]
        + add read_auxiliary_info()
v0.3 Dec. 19, 2016
    - add wrap_fromfile()
    - read [inprodR],[prev_err],[resid_scale]
v0.2 Dec. 19, 2016
    - fix bug > read [local_nRows] for "size_t" type
v0.1 Dec. 18, 2016
    - read_chpoint_file() > read [ind_m],[local_nRows],[niter],[counter]
'''

np.set_printoptions(precision=6)

def wrap_fromfile(rfp, typecode, num):
    res = array.array(typecode)
    res.fromfile(rfp, num)
    return res


def read_auxiliary_info(auxFilename):
    with open(auxFilename, "rb") as auxfp:
        sc_N = wrap_fromfile(auxfp, 'i', num=1)  # the number of scalars' sizes

        sc_sizes = array.array('i')  #
        sc_sizes.fromfile(auxfp, *sc_N)  # '*':unpack list

        vec_N = wrap_fromfile(auxfp, 'i', num=1)  # the number of vercotrs' sizes

        vec_sizes = array.array('i')  #
        vec_sizes.fromfile(auxfp, *vec_N)  # '*':unpack list

        return sc_N, sc_sizes, vec_N, vec_sizes
    # TODO: 0m > error handling when file is NULL


def read_chpoint_file(chpFilename, auxFilename):
    print(chpFilename)
    print(auxFilename)

    # TODO: 0m > error handling when file is NULL
    with open(chpFilename, "rb") as chpfp:
        ind_m = wrap_fromfile(chpfp, 'i', num=1)  # index of iterative method
        print('ind_m:', ind_m)
        # number of local rows of decomposition (only real dipoles)
        #    where L: unsigned long = size_t
        local_nRows = wrap_fromfile(chpfp, 'L', num=1)
        print('local_nRows:', local_nRows)
        niter = wrap_fromfile(chpfp, 'i', num=1)  # iteration count
        print('niter:', niter)
        # number of successive iterations without residual decrease
        counter = wrap_fromfile(chpfp, 'i', num=1)
        print('counter:', counter)
        # used as |r_0|^2 and best squared norm of residual up to some iteration
        inprodR = wrap_fromfile(chpfp, 'd', num=1)
        print('inprodR:', inprodR)
        # previous relative error; used in ProgressReport, initialized in IterativeSolver
        prev_err = wrap_fromfile(chpfp, 'd', num=1)
        print('prev_err:', prev_err)
        # scale to get square of relative error
        resid_scale = wrap_fromfile(chpfp, 'd', num=1)
        print('resid_scale:', resid_scale)  # scale to get square of relative error

        # sc_N: the number of scalars' sizes
        # sc_sizes:
        # vec_N: the number of vercotrs' sizes
        # vec_sizes:
        sc_N, sc_sizes, vec_N, vec_sizes = read_auxiliary_info(auxFilename)
        print('sc_N:', sc_N)
        print('sc_sizes', sc_sizes)
        print('vec_N:', vec_N)
        print('vec_sizes:', vec_sizes)
        # scalars[]
        scalars = []
        for idx in range(sc_N[0]):
            if sc_sizes[idx] is 8:  # double
                scl = wrap_fromfile(chpfp, 'd', num=1)
                scalars.append(scl)
            elif sc_sizes[idx] is 16:  # double complex
                scl = wrap_fromfile(chpfp, 'd', num=2)
                scalars.append(scl)
        for idx in range(sc_N[0]):
       		print('scalars:', scalars[idx])
       	# xvec[] : total electric field on the dipoles
       	# TODO: 0m > replace 'local_nRows[0]' with function?
       	# TODO: 0m > replace 'vec_N[0]' with function?
       	# TODO: 0m > define function for read [drel,dimg]??
       	xvec = [[] for idx in range(local_nRows[0])]
       	for idx in range(local_nRows[0]):
       		drel = wrap_fromfile(chpfp, 'd', 1)
       		dimg = wrap_fromfile(chpfp, 'd', 1)
       		xvec[idx] = np.array([drel, dimg])
       	# rvec[] : current residual
       	rvec = [[] for idx in range(local_nRows[0])]
       	for idx in range(local_nRows[0]):
       		drel = wrap_fromfile(chpfp, 'd', 1)
       		dimg = wrap_fromfile(chpfp, 'd', 1)
       		rvec[idx] = np.array([drel, dimg] )      
       	# pvec[] : polarization of dipoles, also an auxiliary vector in iterative solvers
       	pvec = [[] for idx in range(local_nRows[0])]
       	for idx in range(local_nRows[0]):
       		drel = wrap_fromfile(chpfp, 'd', 1)
       		dimg = wrap_fromfile(chpfp, 'd', 1)
       		pvec[idx] = np.array([drel, dimg])
       	# vectors[] : ???
       	MAXNUM_VECTORS = 20
       	vectors = [[] for idx in range(MAXNUM_VECTORS)]
       	for idx_vecN in range(vec_N[0]):
       	    alist = []
            for idx_LnRows in range(local_nRows[0]):
                if vec_sizes[idx_vecN] is 16:
                    drel = wrap_fromfile(chpfp, 'd', 1)
                    dimg = wrap_fromfile(chpfp, 'd', 1)
                    alist.append(np.array([drel, dimg]))
            vectors[idx_vecN] = alist
       	#
       	for idx in range(3):
       		print('xvec:',np.array_repr(xvec[idx]).replace('\n',''))
       		print('rvec:',np.array_repr(rvec[idx]).replace('\n',''))
       		print('pvec:',np.array_repr(pvec[idx]).replace('\n',''))
       		print('vectors:',np.array_repr(vectors[idx][0]).replace('\n',''))
       		print('vectors:',np.array_repr(vectors[idx][1]).replace('\n',''))
       		print('vectors:',np.array_repr(vectors[idx][2]).replace('\n',''))

argvs = sys.argv
argc = len(argvs)

print argvs
if (argc < 3):
    print("ERROR: chpoint file is not specified\r\n")
    print("   [cmd] [chpoint file] [auxiliary file]\r\n")
    sys.exit()

read_chpoint_file(chpFilename=argvs[1], auxFilename=argvs[2])

結果

['read_chpoint.py', 'LN-CHP', 'LN-AUX']
LN-CHP
LN-AUX
('ind_m:', array('i', [5]))
('local_nRows:', array('L', [27984L]))
('niter:', array('i', [209]))
('counter:', array('i', [139]))
('inprodR:', array('d', [4.317011265934948]))
('prev_err:', array('d', [0.3198525404317458]))
('resid_scale:', array('d', [0.020384616751203153]))
('sc_N:', array('i', [8]))
('sc_sizes', array('i', [8, 8, 8, 8, 16, 16, 16, 16]))
('vec_N:', array('i', [3]))
('vec_sizes:', array('i', [16, 16, 16]))
('scalars:', array('d', [26.040533833425627]))
('scalars:', array('d', [21.695179646708983]))
('scalars:', array('d', [0.01879089326833417]))
('scalars:', array('d', [0.01622964090344006]))
('scalars:', array('d', [2.393088070369695, 5.705068471314085]))
('scalars:', array('d', [0.2463390807093872, 0.2808495718014472]))
('scalars:', array('d', [0.8838701376516526, -0.4673547711298418]))
('scalars:', array('d', [0.34430517188699605, -0.9387175013645006]))
('xvec:', 'array([[ 0.000916],       [-0.001571]])')
('rvec:', 'array([[-0.0021  ],       [ 0.002827]])')
('pvec:', 'array([[-0.001287],       [-0.008525]])')
('vectors:', 'array([[-0.203491],       [ 0.033598]])')
('vectors:', 'array([[ 0.513358],       [-0.235685]])')
('vectors:', 'array([[ 0.548095],       [ 0.061967]])')
('xvec:', 'array([[ 0.032771],       [ 0.090957]])')
('rvec:', 'array([[ 0.004384],       [-0.023016]])')
('pvec:', 'array([[-0.00586 ],       [ 0.025872]])')
('vectors:', 'array([[ 0.100941],       [-0.212583]])')
('vectors:', 'array([[-0.220141],       [ 0.644602]])')
('vectors:', 'array([[-0.386403],       [ 0.429201]])')
('xvec:', 'array([[-0.004613],       [-0.000545]])')
('rvec:', 'array([[ 0.005385],       [-0.006732]])')
('pvec:', 'array([[-0.000544],       [ 0.022825]])')
('vectors:', 'array([[-0.000433],       [ 0.008886]])')
('vectors:', 'array([[ 0.007084],       [-0.024328]])')
('vectors:', 'array([[ 0.006453],       [-0.022867]])')

C実装結果と比較しやすくなりました。
情報感謝です。

0
1
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?