アイデア

プロット用のデータを毎回計算すると時間がかかりすぎるため、事前計算しておきたい。
そこで、バイナリファイルに書き出すようにしてみた。

複素数関数の、引数の実部と虚部をそれぞれ等差的に変化させたときの格子点に対しての、関数の値の実部と虚部をdouble型で格納。
複数の関数を組み合わせて追加計算したい場合を想定して、複数の関数を同じファイルに格納。

データ構造のメモ書き

ascii4 "PLOT"
uint2  format version info
uint2  number of functions
uint4  stepcount of x(real part) of numpy.linspace(,,*)
uint4  stepcount of y(imag part)
ascii128 description of func[0]
ascii128 description of func[1]
:
ascii128 description of func[n-1]
double start x of numpy.linspace(*,,)
double end x of numpy.linspace(,*,)
double start y
double end y
double func[0].Re[xindex=0,yindex=0]
double func[0].Im[xindex=0,yindex=0]
double func[1].Re[xindex=0,yindex=0]
double func[1].Im[xindex=0,yindex=0]
:
double func[n-1].Re[xindex=0,yindex=0]
double func[n-1].Im[xindex=0,yindex=0]
double func[0].Re[xindex=1,yindex=0]
double func[0].Im[xindex=1,yindex=0]
double func[1].Re[xindex=1,yindex=0]
double func[1].Im[xindex=1,yindex=0]
:
double func[0].Re[xindex=0,yindex=1]
double func[0].Im[xindex=0,yindex=1]
double func[1].Re[xindex=0,yindex=1]
double func[1].Im[xindex=0,yindex=1]
:

ソースコード

プログラム上の要点としては、struct.pack を使ってdouble等の値をバイナリに変換している。

今回使っている複素関数とプロットについては下記参照。

なお以下のコードは、グローバル変数使いまくりだったり、なぐり書きのソースなので作法的にはマネしない方がよいかも。

計算してファイルに出力

おおよそ5秒経過毎に進捗率を表示するようにしている。
※for k in range(len(x[0])):の部分が長くなると5秒ごとに出力できなくなる可能性があるが、処理速度優先で外側のforのほうに入れた。

calc.py

import time
import struct
import numpy as np
import mpmath
mpmath.mp.dps = 15
# mpmath.mp.pretty = True

# -----------------------

def dumpFuncValues(f, x, y):
    global previous_time

    # z = np.zeros((len(x), len(x[0])))
    for i in range(len(x)):
        current_time = time.time()
        if previous_time+5 < current_time :
            previous_time = current_time
            print(str(round(100*i/len(x),3))+'%')
        
        for k in range(len(x[0])):
            t = x[i][k] + y[i][k]*1j
            ze    = mpmath.zeta(t)
            ze_d1 = mpmath.zeta(t,derivative=1)
            ga    = mpmath.gamma(t*0.5)
            diga  = mpmath.digamma(t*0.5)
            p = mpmath.power(mpmath.pi,-0.5*t)
            ect = ga*p*ze

            re = struct.pack('@d', ze.real)
            im = struct.pack('@d', ze.imag)
            f.write(re)
            f.write(im)
            
            re = struct.pack('@d', ze_d1.real)
            im = struct.pack('@d', ze_d1.imag)
            f.write(re)
            f.write(im)
            
            re = struct.pack('@d', ga.real)
            im = struct.pack('@d', ga.imag)
            f.write(re)
            f.write(im)
            
            re = struct.pack('@d', diga.real)
            im = struct.pack('@d', diga.imag)
            f.write(re)
            f.write(im)

            re = struct.pack('@d', p.real)
            im = struct.pack('@d', p.imag)
            f.write(re)
            f.write(im)

            re = struct.pack('@d', ect.real)
            im = struct.pack('@d', ect.imag)
            f.write(re)
            f.write(im)
    return
# -----------------------
def dump():
    x = np.linspace(X_START, X_END, X_COUNT)
    y = np.linspace(Y_START, Y_END, Y_COUNT)

    X, Y = np.meshgrid(x, y)

    func_comments = [
        'RiemannZeta(z)',
        'RiemannZeta(z) derivative',
        'Gamma(z/2)',
        'Digamma(z/2)',
        'pow(pi,-z/2)',
        'prod of (#1,#3,#5)'
    ]

    f = open('myfile.dat', 'wb')
    f.write('PLOT'.encode('utf-8'))
    f.write(struct.pack('@H', 0)) # version info
    f.write(struct.pack('@H', len(func_comments))) # number of functions
    f.write(struct.pack('@L', X_COUNT))
    f.write(struct.pack('@L', Y_COUNT))


    for fc in func_comments:
        f.write(fc.ljust(FUNC_COMMENT_LEN,' ').encode('utf-8'))
    f.write(struct.pack('@d', X_START))
    f.write(struct.pack('@d', X_END))
    f.write(struct.pack('@d', Y_START))
    f.write(struct.pack('@d', Y_END))
    dumpFuncValues(f, X, Y)
    f.close()
    return
# -----------------------

FUNC_COMMENT_LEN = 128

X_STEP = 1.0/32.0
Y_STEP = 1.0/2.0

X_START = 3/8
X_END   = 5/8

Y_START = 10
Y_END   = 200

X_COUNT = round(((X_END-X_START)/X_STEP)+1)
Y_COUNT = round(((Y_END-Y_START)/Y_STEP)+1)

previous_time = time.time()

if X_COUNT * Y_COUNT < 1000000:
    dump()
else:
    f = open('mesh_is_toobig.txt', 'wb')
    f.write('X:')
    f.write(str(X_COUNT))
    f.write(',Y:')
    f.write(str(Y_COUNT))
    f.close()

読み出しとプロット

load_and_plot.py

import sys
import os
import struct
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt

def load_header_info(f):
    global x_count
    global y_count
    global x_start
    global x_end
    global y_start
    global y_end
    global seekpos
    global count_of_funcs
    f.read(4) # skip 'PLOT'
    f.read(2) # skip version info
    count_of_funcs = struct.unpack('@H', f.read(2))[0]
    print(count_of_funcs)
    x_count = struct.unpack('@L', f.read(4))[0]
    y_count = struct.unpack('@L', f.read(4))[0]
    # func_comments = []
    for i in range(count_of_funcs):
        fc = f.read(128).decode()
        # func_comments.append(fc)
    x_start = struct.unpack('@d', f.read(8))[0]
    x_end   = struct.unpack('@d', f.read(8))[0]
    y_start = struct.unpack('@d', f.read(8))[0]
    y_end   = struct.unpack('@d', f.read(8))[0]
    seekpos = f.tell() # get current seek pos

    tmpx = np.linspace(x_start, x_end, x_count)
    tmpy = np.linspace(y_start, y_end, y_count)
    return tmpx,tmpy

def load_nth_funcs_real(f,index):
    global x
    global y
    global count_of_funcs
    z = np.zeros((len(y), len(x)))
    for i in range(len(y)):
        for k in range(len(x)):
            for funci in range(count_of_funcs):
                if funci == index:
                    tmp_re = struct.unpack('@d', f.read(8))[0]
                    f.read(8)
                    z[i][k] = tmp_re
                else:
                    f.read(8)
                    f.read(8)
    return z

def load_nth_funcs_imag(f,index):
    global x
    global y
    global count_of_funcs
    z = np.zeros((len(y), len(x)))
    for i in range(len(y)):
        for k in range(len(x)):
            for funci in range(count_of_funcs):
                if funci == index:
                    f.read(8)
                    tmp_im = struct.unpack('@d', f.read(8))[0]
                    z[i][k] = tmp_im
                else:
                    f.read(8)
                    f.read(8)
    return z

x_count=0
y_count=0
x_start=0
x_end=0
y_start=0
y_end=0
seekpos=0
count_of_funcs=0


args = sys.argv
path = None
if len(args)>1:
    tmp_path = args[1]
    dummy,ext = os.path.splitext(tmp_path)
    if ext == ".dat":
        path = tmp_path
else:
    path = 'myfile.dat'

f = open(path, 'rb')
x,y = load_header_info(f)
#Z = load_nth_funcs_real(f,0)
# f.seek(seekpos)
Z = load_nth_funcs_real(f,0)
f.close()
X, Y = np.meshgrid(x, y)

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(elev=20, azim=225)
# ax.set_zlim3d(top=Z_LIMIT)
# ax.set_ylim3d(top=50,bottom=10)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis', edgecolor='none')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()