ページ

2013年12月1日日曜日

Mathmaticaを用いた数値対角化の計算時間

Mathematicaを用いた対角化の計算速度をplotしてみました.

対角化した行列のclassは
ランダムではなくスパーズではない一般のエルミート行列(赤ドット)と
ランダムではなくスパースではない一般のユニタリー行列(青ドット).
また,行列は200桁精度で近似評価されています.

厳密評価のまま計算するとキットモット遅い.

行列の構成方法は気にしないと言う事でご勘弁.

実線は計算時間を3次の多項式でfittingしたもの.
最後に描画とfittingのプログラムを置いときます.
(もちろん計算時間のデータファイルがないと動きませんが・・・
欲しいなら上げますが,欲しいか?)




import numpy as np
import matplotlib.pyplot as plt

def f(x,para):
    return para[0]*x**3 + para[1]*x**2 + para[2]*x + para[3]
    
data = np.loadtxt("time_prec200.dat").transpose()
fig = plt.figure(figsize=(8,4))
ax0 = fig.add_subplot(1,2,1)
ax1 = fig.add_subplot(1,2,2)
colors=["red", "blue"]
for i, ax in enumerate([ax0,ax1]):
    ax.plot(data[0],data[1]/3600,'o', label='Hermitian matrix',markersize=10, color=colors[0])
    ax.plot(data[0],data[2]/3600,'o', label='Unitary matrix',markersize=10,color=colors[1])
#    ax.plot(data[0],data[3]/3600,'og')
    for j in range(1,3):
        x = np.arange(100, 10000,1)    
        para = np.polyfit(data[0],data[j],3)
        y = f(x,para)
        ax.plot(x,y/3600,'-',linewidth=3,color=colors[j-1])
    if i==1:
        ax.loglog()
        ax.set_xlim(50,10000)
    else:
        ax.set_ylim(0,5)
        ax.set_xlim(0,1000)
        ax.set_ylabel("calculation time (hours)",fontsize=16)            
    le = ax.legend(loc=2,numpoints=1)
    le.draw_frame(False)
    ax.grid()
    ax.set_xlabel("matrix dimension",fontsize=16)
sup = plt.suptitle("diagonalization time with Mathematica (precision=200)",fontsize=15)
sup.set_position((0.5,1.0))
plt.tight_layout()
plt.savefig("calc_time.png")
plt.show()

2013年11月21日木曜日

Maximaによる任意精度計算

Mathematicaはいつ使えなくなるか(ライセンス的な関係で)分らないので,Maximaを任意精度数値計算に使えないか試してみる.


精度(precision)のコントロールは

fpprec:50

演算による多倍長浮動小数点の有理数への変換を制御するため

bftorat: true

にする必要がある.デフォルトではfalse.
falseの場合多倍長浮動小数点は小さい有理数で近似される.

また,出力する桁のコントロールは

fpprintprec:50;

で可能.
maxima 5.30.0 だとload("fft")が出来なかったので,古いバージョンを使う.
参考:http://maxima.sourceforge.jp/maxima_5.html
$ maxima
Maxima 5.21post http://maxima.sourceforge.net
using Lisp SBCL 1.0.39
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) fpprec;
(%o1) 16
(%i2) fpprec:50;
(%o2) 50
(%i5) bftorat:true;
(%o5) true
(%i6) x:bfloat([0,1/3,0,0]);
(%o6) [0.0b0,3.3333333333333333333333333333333333333333333333333b−1,0.0b0,0.0b0]
(%i7) load("fft");
(%o7) "/Applications/Maxima.app/share/maxima/5.21post/share/numeric/fft.lisp"
(%i10) res:inverse_fft(fft(x));
(%o10) [0.0,.3333333333333333,0.0,0.0]
(%i11) bfloatp(res[1]);
(%o11) false

FFTは任意精度計算やってくれないっぽい.一方で行列mat=((i,pi),(pi,i))の対角化は任意精度でやってくれるっぽい.上のmatの固有値は(±pi, i)である.

(%i12) mat:bfloat(matrix([%i,%pi],[%pi,%i]));
(%o12) matrix([%i,3.1415926535897932384626433832795028841971693993751b0],[3.1415926535897932384626433832795028841971693993751b0,%i])
(%i13) load("eigen");
(%o13) "/Applications/Maxima.app/share/maxima/5.21post/share/matrix/eigen.mac"
(%i14) res:eigenvalues(mat);
`rat' replaced -9.8696044010893586188344909998761511353136994072408B0 by -41494057121218467653886768/4204226981644629723547549 = -9.8696044010893586188344909998761511353136994072408B0
(%o14) [[(4204226981644629723547549*%i−4*sqrt(10903152157933135742303986703335184117163956870727))/4204226981644629723547549,(4204226981644629723547549*%i+4*sqrt(10903152157933135742303986703335184117163956870727))/4204226981644629723547549],[1,1]]
(%i24) bfloat(realpart(res[1][1])); imagpart(res[1][1]);
(%o24) −3.1415926535897932384626433832795028841971693993751b0
(%o25) 1
(%i22) -bfloat(%pi);
(%o22) −3.1415926535897932384626433832795028841971693993751b0
(%i26) bfloat(realpart(res[1][2])); imagpart(res[1][2]);
(%o26) +3.1415926535897932384626433832795028841971693993751b0
(%o27) 1
---------- Mathmatica res ----------
In[2]:= N[Pi, 50]
Out[2]= 3.1415926535897932384626433832795028841971693993751

maxima対角化は時間が掛かりすぎるから現実的ではない・・・か
(%i1) f[i,j]:=random(1.0);
(%o1) f[i,j]:=random(1.0)
(%i6) showtime:true;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o6) true
(%i12) mat:genmatrix(f,4,4);
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o12) matrix([.9138095996128959,.1951846177977887,.4823905248516196,.3554400571592862],[.4155627227214829,.6125738892538246,.8418126553359213,.5055894446438118],[0.0363933158999219,.6252975183418072,.9937314578694818,.04158924615444981],[.5931521815169765,.6546577661145074,.4707106740336644,.05192355366451751])
(%i13) eigenvalues(mat);
Evaluation took 9.2310 seconds (14.0640 elapsed) using 2764.567 MB.

2013年11月16日土曜日

Mathematicaの糖衣構文

Mathematicaの関数へのアクセスは糖衣構文によって色々な人に受け入れやすい事になっています.
多分,同じ数学概念でも棲んでいる世界によって解釈が異なったりするからだと思います.ただ,私には甘すぎては苦手です.

可読性って概念はMathematicaにはないのかね…
他人のプログラム読む時,絶望を覚えるんですが・・・

例えば以下の例は全てSqrt[3]を返します.

In[1]:= f[x_, y_: 2] := x^(1/y)

In[2]:= f[3]
f[3, 2]
3~f~2
3 // f
3 // f[#, 2] &
f @ 3
f @@ {3}
f @@ {3, 2}
f[#, 2] &  @ 3
Apply[f, {3}]
Apply[f, {3, 2}]
Apply[f, {##}] &  @@ {3, 2}


Out[2]= Sqrt[3]

Out[3]= Sqrt[3]

Out[4]= Sqrt[3]

Out[5]= Sqrt[3]

Out[6]= Sqrt[3]

Out[7]= Sqrt[3]

Out[8]= Sqrt[3]

Out[9]= Sqrt[3]

Out[10]= Sqrt[3]

Out[11]= Sqrt[3]

Out[12]= Sqrt[3]

Out[13]= Sqrt[3]

最近Mathematicaに関する事ばかりだな・・・

Mathematica Fourier変換について

バグを疑い始めたので,プログラムを純化させてみる.
In[1]:= IdentityTr[x_] := InverseFourier[Fourier[x]]

In[2]:= a = {0, 0, 0.33333333333333333, 0}
b = {0, 0, 0.333333333333333333, 0}
Precision[a[[3]]]
Precision[b[[3]]]

Out[2]= {0, 0, 0.333333, 0}

Out[3]= {0, 0, 0.33333333333333333, 0}

Out[4]= MachinePrecision

Out[5]= 17.5229

In[8]:= NestList[IdentityTr, a, 20]
NestList[IdentityTr, b, 20]

Out[8]= {{0, 0, 0.333333, 0}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 
  0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 
  0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 
  0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 
  0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 
  0.333333, 0.}, {0., 0., 0.333333, 0.}, {0., 0., 0.333333, 0.}}

Out[9]= {{0, 0, 0.33333333333333333, 0}, {0.*10^-18, 
  0.*10^-18 + 0.*10^-18 I, 0.33333333333333333, 
  0.*10^-18 + 0.*10^-18 I}, {0.*10^-17, 0.*10^-18 + 0.*10^-18 I, 
  0.3333333333333333, 0.*10^-18 + 0.*10^-18 I}, {0.*10^-17, 
  0.*10^-17 + 0.*10^-17 I, 0.333333333333333, 
  0.*10^-17 + 0.*10^-17 I}, {0.*10^-16, 0.*10^-16 + 0.*10^-16 I, 
  0.333333333333333, 0.*10^-16 + 0.*10^-16 I}, {0.*10^-15, 
  0.*10^-15 + 0.*10^-15 I, 0.33333333333333, 
  0.*10^-15 + 0.*10^-15 I}, {0.*10^-14, 0.*10^-15 + 0.*10^-15 I, 
  0.3333333333333, 0.*10^-15 + 0.*10^-15 I}, {0.*10^-14, 
  0.*10^-14 + 0.*10^-14 I, 0.333333333333, 
  0.*10^-14 + 0.*10^-14 I}, {0.*10^-13, 0.*10^-13 + 0.*10^-13 I, 
  0.33333333333, 0.*10^-13 + 0.*10^-13 I}, {0.*10^-12, 
  0.*10^-12 + 0.*10^-12 I, 0.33333333333, 
  0.*10^-12 + 0.*10^-12 I}, {0.*10^-11, 0.*10^-11 + 0.*10^-11 I, 
  0.3333333333, 0.*10^-11 + 0.*10^-11 I}, {0.*10^-10, 
  0.*10^-11 + 0.*10^-11 I, 0.333333333, 
  0.*10^-11 + 0.*10^-11 I}, {0.*10^-10, 0.*10^-10 + 0.*10^-10 I, 
  0.33333333, 0.*10^-10 + 0.*10^-10 I}, {0.*10^-9, 
  0.*10^-9 + 0.*10^-9 I, 0.33333333, 
  0.*10^-9 + 0.*10^-9 I}, {0.*10^-8, 0.*10^-8 + 0.*10^-8 I, 0.3333333,
   0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-8 + 0.*10^-8 I, 0.333333,
   0.*10^-8 + 0.*10^-8 I}, {0.*10^-7, 0.*10^-7 + 0.*10^-7 I, 0.33333, 
  0.*10^-7 + 0.*10^-7 I}, {0.*10^-6, 0.*10^-6 + 0.*10^-6 I, 0.3333, 
  0.*10^-6 + 0.*10^-6 I}, {0.*10^-5, 0.*10^-5 + 0.*10^-5 I, 0.3333, 
  0.*10^-5 + 0.*10^-5 I}, {0.*10^-4, 0.*10^-4 + 0.*10^-4 I, 0.333, 
  0.*10^-4 + 0.*10^-4 I}, {0.*10^-3, 0.*10^-4 + 0.*10^-4 I, 0.33, 
  0.*10^-4 + 0.*10^-4 I}}

2013年11月9日土曜日

Mathematicaでの任意制度計算

Mathematicaが手に入ったので,いろいろテストしてみました.
その結果, 任意精度計算を行うにあたって精度の保持はユーザーが監視しないと行けない事が分りました.
また,監視を怠ると正しい結果が得られない事も分りました.

簡単な例として離散Fourier変換(FFTか?)と逆変換を繰り返す演算をさせます.

Wolfram Mathematica Document Foureirの説明にある通り,
Inputのlistが厳密な数値である場合,
N関数を適用し精度はMachinePrecisionに落とされます.
そのためFoureir関数は厳密精度計算を行ってくれません.

しかしながらユーザーがInputのlistにN関数を適用し精度をユーザー側でコントロールしたい状況が発生しますが,安易な方法では正しい結果が得られません.

その結果が In[3]です.

In[1]:= func1[y_, iter_, dps_: 10] := Module[{i, x},
  x = y;
  For[i = 0, i < iter, i++;
   If[ Mod[i, 2] != 0, x = SetPrecision[Fourier[x], dps], 
    x = SetPrecision[InverseFourier[x], dps]];
   ];
  Return[x]]
func0[y_, iter_] := Module[{i, x},
  x = y;
  For[i = 0, i < iter, i++;
   If[ Mod[i, 2] != 0, x = Fourier[x], x = InverseFourier[x]];
   ];
  Return[x]]


In[3]:= a = {0, 0, 1/3, 0}
res = func0[a, 100]
Precision[res[[1]]]
$MachinePrecision

b = N[{0, 0, 1/3, 0}, 10]
res = func0[b, 100]
Precision[res[[1]]]


Out[3]= {0, 0, 1/3, 0}

Out[4]= {0., 0., 0.333333, 0.}

Out[5]= MachinePrecision

Out[6]= 15.9546

Out[7]= {0, 0, 0.3333333333, 0}

Out[8]= {0.*10^5, 0.*10^4 + 0.*10^4 I, 0.*10^5, 0.*10^4 + 0.*10^4 I}

Out[9]= 0.
listに厳密な数値を代入した場合,正しい結果がMachinPrecisionで帰ってきます.
一方listにN関数を作用させて作った初期条件の場合,全ての結果のPrecisionが0になっています.

これでは使い物になりません.1回のFourier関数の演算でどうやら10^2-10^10程度の桁落ちが発生し最終的にゼロが帰ってくる様です.

上記問題を回避する為にはSetPrecision関数を用いて,強制的に精度を維持する事です
(func1関数を参照).
ただしこの方法が正当な回避方法かどうか私には分りませんが,少なくとも正しい答えが望んだ精度で帰ってくる様です(In[10]以降).
In[10]:= a = {0, 0, 1/3, 0}
res = func1[a, 100, 10]
Precision[res[[3]]]
b = N[{0, 0, 1/3, 0}, 10]
res = func1[b, 100, 10]
Precision[res[[3]]]

Out[10]= {0, 0, 1/3, 0}

Out[11]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333, 
 0.*10^-21 + 0.*10^-21 I}

Out[12]= 10.

Out[13]= {0, 0, 0.3333333333, 0}

Out[14]= {0, 0.*10^-21 + 0.*10^-21 I, 0.3333333333, 
 0.*10^-21 + 0.*10^-21 I}

Out[15]= 10.

In[34]:= Exit[]

2013年10月27日日曜日

mpmath 速度比べ

import mpmath
import numpy as np
import time

sample=500
iter=500

x = np.array([1.0/3.0]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + 1.0/3.0
    else:
        x = x - 1.0/3.0
goal=time.time()
print("---- numpy array -----")
print("%.20f" % init)
print("%.20f" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 18
x = np.array([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath with numpy array (dps=%d) -----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 100
x = np.array([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath with numpy array (dps=%d) -----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 18
x = mpmath.matrix([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath matrix (dps=%d)-----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 100
x = mpmath.matrix([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath matrix (dps=%d)-----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)


x = [mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = [xx + mpmath.mpf("1.0")/mpmath.mpf("3.0") for xx in x]
    else:
        x = [xx - mpmath.mpf("1.0")/mpmath.mpf("3.0") for xx in x]
goal=time.time()
print("---- list -----")
print("%.20f" % init)
print("%.20f" % x[0])
print("time:", goal-start)




実行結果
$ python speed_test_array.py 
---- numpy array -----
0.33333333333333331483
0.33333333333332632042
('time:', 0.0022330284118652344)
---- mpmath with numpy array (dps=18) -----
0.333333333333333333
0.333333333333333326
('time:', 1.7795300483703613)
---- mpmath with numpy array (dps=100) -----
0.333333333333333333
0.333333333333333333
('time:', 2.18923282623291)
---- mpmath matrix (dps=18)-----
0.333333333333333333
0.333333333333333326
('time:', 5.815176010131836)
---- mpmath matrix (dps=100)-----
0.333333333333333333
0.333333333333333333
('time:', 6.364114046096802)
---- list -----
0.33333333333333331483
0.33333333333333331483
('time:', 12.696530103683472)

2013年10月15日火曜日

gccの4倍精度ライブラリー

gcc 4.6 から4倍精度計算ができる様になりました.昔は4倍精度で計算できるマシンがあったみたいだけどね.

今回は、c, (c++はやってない), fortranの言語で実装するとしてどの言語が最も実装しやすいか比較する.
gcc のversionは4.8 を用いる.結論としてはfortranが一番楽だと思う.
cおよびc++は
http://homepage3.nifty.com/maikaze/float128.html
で解説されているようにc(, c++)では幾分実装が不自然な所もあるので,一応比較してみる.
計算する内容は1/3, 4*arctan(1)=pi,sqrt(2)である.

まずfortranから
program fort_quadmath
 implicit none
 real (kind(0D0)) dx, dy, dpi,dsqrt2
 real (kind(0Q0)) qx, qy, qpi,qsqrt2
 dx = 1.0D0;
 dy = 3.0D0;
 dsqrt2 = sqrt(2.0D0)
 dpi = 4.0D0*atan(1.0D0)
 qx = 1.0Q0;
 qy = 3.0Q0;
 qsqrt2 = sqrt(2.0Q0)
 qpi = 4.0Q0*atan(1.0Q0)
 print *, "---- double precision ----"
 print *, "1/3     =", dx/dy
 print *, "pi      =", dpi
 print *, "sqrt(2) =", dsqrt2
 print *, "----  quad  precision ----" 
 print *, "1/3     =", qx/qy
 print *, "pi      =", qpi
 print *, "sqrt(2) =", qsqrt2
end program fort_quadmath 
コンパイルは
 gfortran -o fort_quadmath fort_quadmath.f03 
であり,実行結果は
$ ./fort_quadmath 
 ---- double precision ----
 1/3     =  0.33333333333333331     
 pi      =   3.1415926535897931     
 sqrt(2) =   1.4142135623730951     
 ----  quad  precision ----
 1/3     =  0.333333333333333333333333333333333317      
 pi      =   3.14159265358979323846264338327950280      
 sqrt(2) =   1.41421356237309504880168872420969798   
 pi           3.1415926535897932384626433832795028841971693993751
 sqrt(2)      1.4142135623730950488016887242096980785696718753769
厳密解と比較して32桁くらいまでは遜色無く一致している.
実装に関して変数の型の宣言と、リテラル(0Q0)さえ気をつけていれば,関数の型を意識する必要がない(というか関数が代表名で与えられているので,変数の型によって関数が勝手に陰的に変わるんだと理解している.)

一方c言語に関して
#include <stdio.h>
#include <quadmath.h>
#include <math.h>

int main(void)
{
 char c1[80],c2[80],c3[80],c4[80],c5[80];
 double dx,dy,dpi;
 long double lx,ly,lpi;
 __float128 x,y, pi1,pi2,sqrt1,sqrt2;
 dx = 1.0;
 dy = 3.0;
 dpi = 4.0*atan(1.0);
 lx = 1.0L;
 ly = 3.0L;
 lpi = 4.0L*atan(1.0L);
 x = 1.0Q;
 y = 3.0Q;
 pi1 = 4.0Q*atan(1.0Q);
 pi2 = 4.0Q*atanq(1.0Q); 
 sqrt1 = sqrt(2.0Q);
 sqrt2 = sqrtq(2.0Q);  
 
 printf("---- double precision ----\n");
 printf("1/3 = %.30lf\n", dx/dy);
 printf("pi  = %lf\n", dpi); 
 printf("---- long double prescision ----\n");
 printf("1/3 = %.30Lf\n", lx/ly);
 printf("pi  = %.30Lf\n", lpi);
 printf("---- quad precision ----\n");
 quadmath_snprintf(c1,80,"%.40Qf", x/y);
 quadmath_snprintf(c2,80,"%.40Qf", pi1); 
 quadmath_snprintf(c3,80,"%.40Qf", pi2);  
 quadmath_snprintf(c4,80,"%.40Qf", sqrt1);  
 quadmath_snprintf(c5,80,"%.40Qf", sqrt2);    
 printf("1/3 = %s\n", c1);
 printf("pi1  = %s\n", c2);
 printf("pi2  = %s\n",c3);
 printf("sqrt1  = %s\n",c4);
 printf("sqrt2  = %s\n",c5);  
 return 0;
}
のコンパイルは
gcc -o c_quadmath c_quadmath.c -lquadmath
で行う.実行結果は
$ ./c_quadmath 
---- double precision ----
1/3 = 0.333333333333333314829616256247
pi  = 3.141593
---- long double prescision ----
1/3 = 0.333333333333333333342368351437
pi  = 3.141592653589793115997963468544
---- quad precision ----
1/3 = 0.3333333333333333333333333333333333172839
pi1  = 3.1415926535897931159979634685441851615906
pi2  = 3.1415926535897932384626433832795027974791
sqrt1  = 1.4142135623730951454746218587388284504414
sqrt2  = 1.4142135623730950488016887242096981769402
 
pi     3.1415926535897932384626433832795028841971693993751
sqrt(2)  1.4142135623730950488016887242096980785696718753769
の厳密解と比較して pi1とsqrt1がdouble程度の精度しかない. これは代入する変数の型と関数の返り値の型が一致していないため,関数の方の型で丸められてしまったせいだと考えられる.(c++でも?)
既存のプログラムの精度を倍精度から4倍精度に書き直すにはc言語では関数も書き換えなければならない. (gccがversion upしたらなんとかなるんか?) 精度の仕様変更に関しては圧倒的にfortranの方がいじりやすいと個人的に思う.
 (fortran はあんまり使ってないけどね・・・)

2013年10月9日水曜日

gnuplot animation in C

 とある需要で,c言語でreal time animationを作る事になった.
gnuplotで描画させるのが一番簡単だろうなーと思ったので,基本的な操作方法をまとめる.
gnuplotを立ち上げて
gnuplot> plot "-"
と打つとデータの入力待ちになるので,
gnuplot> plot[-1:1][-1:1] "-" w p ps 7 pt 7
input data ('e' ends) > 0.0 0.0 
input data ('e' ends) > e
とかしてやると,データを手動で入力することができる.gnuplotを日常的に使っている人ならば直感的に,
gnuplot> plot[-1:1][-1:1] "-" w p ps 7 pt 7,"-" w p ps 7 pt 7
input data ('e' ends) > 0.5 0.0
input data ('e' ends) > e
input data ('e' ends) > -0.5 0.0
input data ('e' ends) > e
と二つ以上のデータ点を入力することができることに気がつくだろう.
結果は上図の様になる.つまり,cからgnuplotに上記のような命令でデータを入力してやれば,real time にデータを更新しながら描画することが可能になる.
(note:multiplotを利用するとデータ更新毎に再描画してしまうので,ちらつきの原因となる)
さて,ここからが本題.C言語からgnuplotに命令するため,ここではpopen関数をつかう.(多分やり方は色々有るはず)
FILE *gp;
gp = popen("gnuplot", "w");
でgnuplotにパイプを通すことができる.もちろんgnuplotのPATHが通ってなければならない.絶対PATHでも構わない. もし描画終了後にgnplotのウィンドウ画面が閉じるのが嫌であればgnuplotの起動オプション
gp = popen("gnuplot -persist", "w"); 
を付けてやれば良い.

以下簡単なサンプル
#include <stdio.h>
#include <stdlib.h>
#include <tmath.h>
#include <unistd.h>

void gnuplot_config(FILE *gp);
void plot2d(FILE *gp,const int num,const double *x,const double *y);
void move(const int num, double const t, double *x, double *y);
int main(void){
    FILE *gp;

    //gp = popen("gnuplot -persist", "w");    /*プログラムが終了後GnuplotのWindowを閉じない*/
    gp = popen("gnuplot", "w");               /*プログラムが終了後GnuplotのWindowを閉じる*/
    gnuplot_config(gp);  /* gnuplot configuration */


    int i;          // loop variable
    int step=100;
    int sample_num=10;
    double t;
    double dt=0.1;
    double x[sample_num], y[sample_num];
    
    printf("途中で終了する場合は,Ctrl-C をターミナルに入力\n");
    
    for (i=0;i<step;i++){
        fprintf(gp,"set title 'animated-%d/%d' font 'Times, 20'\n",i,step);    
        t = i*dt;
        move(sample_num, t, x, y);
        plot2d(gp,sample_num, x, y); 
        usleep(30000); //指定された(マイクロ?)秒プログラムを停止させます (Unix system only)
//        while(getchar() != '\n');       /*Enter keyが入力されるまで,計算停止*/        
    }
    printf("計算終了\n");


    return 0;
}

void gnuplot_config(FILE *gp)
{
/*
    fprintf(gp,"set terminal gif animate optimize size \n");
    fprintf(gp,"set output 'anime.gif'\n");
*/
    /*アニメーションをgif動画として保存する場合は,上のコメントを解除*/    
    fprintf(gp,"set xrange [-1.5:1.5]\n");
    fprintf(gp,"set yrange [-1.5:1.5]\n");
    fprintf(gp,"unset tics\n");
//    fprintf(gp,"set xtics 0.5\n");
//    fprintf(gp,"set ytics 0.5\n");
//    fprintf(gp,"set xlabel 'x' font 'Times, 20'\n");
//    fprintf(gp,"set ylabel 'y' font 'Times, 20'\n");
//    fprintf(gp,"set grid\n");
    fprintf(gp,"set border 0\n");
    fprintf(gp,"set size square\n");
}

void plot2d(FILE *gp, const int num, const double *x, const double *y)
{
    int i;
    fprintf(gp,"plot");
    for(i=0;i<num;i++){
        if (i!=num-1)
            fprintf(gp,"'-' with p pt 7 ps 10 lc %d tit '',", i);//plotデータの読み込み開始
        else
            fprintf(gp,"'-' with p pt 7 ps 10 lc %d tit ''\n", i);//plotデータの読み込み開始
    }
    for(i=0;i<num;i++){
        fprintf(gp,"%f %f\n",x[i],y[i]);
        fprintf(gp,"e\n"); //データ書き込み終了を知らせる.
    }
    fflush(gp); //バッファーに溜まったデーターを吐き出す.これをしないと,real time animation にならない.
}


void move(const int num, const double t, double *x, double *y)
{
    int i;
    double init_phase[num];
    for (i =0;i<num;i++){
        init_phase[i] =  i*2.0*M_PI/num;
    }

    for (i=0;i<num;i++){
        x[i] = cos(t+init_phase[i]);
        y[i] = sin(t+init_phase[i]);
    }
}
計算結果はこんな感じに
HTMLって鍵括弧(<>)を&lt &gtって入力しないと行けないのか・・・プログラムコピペする時めんどくさいな・・・

今回はspleep関数使ってるけど,FPSとかってどうやって制御すれば良いんだ?

Hello world

男は黙ってHello world.

print("Hello world")