Risa2000によるグラフ表示

はじめに

Risa/Asirの標準のグラフ表示では与えられた多項式関数の零点の表示で、自己交差を含むグラフも綺麗に表示可能な点が挙げられる。この表示機能は他の数式処理システムと比較してもユニークな機能である。現在、ソース公開版のRisa2000ではグラフ表示で従来のrshを用いる方法から、OpenXMを用いる方法に移行しており、初期化の方法が異っている。ここではRISA/ASIR のソース公開版RISA2000に関して、そのグラフ表示の方法を簡単に述べる。

Risa2000は一見すると旧版のRisaとは大きな違いは無い様に見える。実際、通常の命令に大きな変更は特に無い様であるが、分散処理の処理の方法がOpenXMを利用するものに変更された為、グラフ表示の方法が大幅に変更され、従来のrshを用いる方法では表示が出来なくなっている。他には、利用しているpariのライブラリが従来の1.38から2.0.17に上っているが、こちらは内部的なもので、グラフ表示の方法が変更になった方が遊ぶには影響が大きい様に思える。

Risa/Asirが持つグラフ表示機能で遊ぶ話

RISA/ASIR2000でグラフ表示はOpenXMを介して行う。このOpenXM (Open message eXchange protocol)は分散計算に於ける通信プロコトルであり、次の手順で分散計算が行われる。

ClientがServerに対して計算実行依頼のメッセージを送る。
serverで計算を実行する。
clientかserverに対して計算結果送付依頼のメッセージを送る。
serverが計算結果を返し、clientは結果を受け取る。

このOpenXMによる分散計算を行う為には、serverを先ず立ち上げて通信を成立させる必要がある。

分散計算を開始する場合、一般的にはox_launch命令を用いる。このox_launchでは引数としてhost名、ディレクトリ、実行する命令を指定する。ここでhost名はasirが動作しているマシン上で命令を実行する場合は0を指定し、その際、ディレクトリに実行する命令が置かれているディレクトリを指定する。尚、asirの環境変数で設定されたディレクトリに置かれた命令を実行する場合、このディレクトリの指定は不要である。例えば、Asir組込のグラフ表示機能を使ってグラフ表示を行いたい場合には、単に、ox_launch(0,"ox_plot");を事前に実行するだけで十分である。ここで、Asirの動作環境がUNIXである場合、ログのエコーバック用にxtermが起動される。尚、Xが使えない環境では、ox_launch_noxを用いる。これらの引数はox_lauchと同じである。尚、ox_launchを実行すると、新規にxtermが起動されて、そこにメッセージが出力される。

Asir-Contrib xm version 20010310. Copyright 2000-2001, OpenX
M Committers.
help("keyword"); ox_help(0); ox_help("keyword"); ox_grep("ke
yword");
     for help messages. 
I'm an ox_plot, Version 20020301.
#0 Got OX_COMMAND SM_mathcap

#1 Got OX_COMMAND SM_popCMO
pop at 0

#2 Got OX_DATA -> data pushed
#3 Got OX_COMMAND SM_setMathcap
pop at 0

このox_launchによる初期化が行えた後は、以前のtcpinitによる初期化を行った後と同様に、ifplot等の表示命令を用いてグラフ表示が可能となる。尚、旧版のrcpを介したグラフ表示と比べ、.rhostの設定も不要であるが、迂濶にグラフ表示ウインドウを破棄した場合、ox_launchによる再初期化が必要な点は同じである。尚、内部的にはdraw_objcet等を用いてグラフ表示を行なっている。

尚、OpenXMを介してgnuplotを用いる事も可能である。この場合、最初に gnuplotを用いる為のパッケージを読込み、それから初期化を行う為に gnuplot_start();を実行する。すると、ox_launchを立ち上げた時の様にエコーバック専用のウインドウが開き、以降、gnuplotを用いたグラフ表示が可能となる。尚、ox_launch(id,"ox_plot")とgnuplot_start()は両立し得るもので、どちらかが使えなくなる様な排他的なものでは無い。

組込関数ifplot,conplot,plot,plotoverについて

最初にASIR単体が有する数式のグラフ表示について述べる。数式のグラフ表示を行う関数にはifplot,conplot,plot,plotoverがある。これらの関数では2次元グラフのみを表示するものである。ここでは簡単な説明と、その例を示す。尚、より原始的な命令にopen_canvas,draw_obj等があるが、これらをifplot等は用いている。この原始的な関数に関しては説明しない。

ifplot

ifplot( 2変数関数f(x,y) (,[w1,w2],[x,x0,x1],[y,y0,y1],id,"タイトル") );

ifplotは2変数多項式関数に対して、その解集合の表示を行う。()で括った部分はオプションである。先ず、[w1,w2]で表示ウインドウの大きさを指定する。何も指定しない場合には、300x300 の大きさのウインドウが出力される。次の[x,x0,x1],[y,y0,y1]で 2変数関数f(x,y)の表示領域を定め、idがサーバのid(ox_launchで指定した番号。localhostに表示するならば0で良い)を指定し、 "タイトル"にはウインドウのタイトルバーに表示するタイトルを指定する。

例

 [349] ifplot(x^3+x^2-y^2);

[350] ifplot(x^5-x*y^2+y^5,[x,-1,1],[y,-1,1],0,"Neko");

conplot

conplot(2変数関数f(x,y) (,[w1,w2],[x,x0,x1],[y,y0,y1],[level,l0,l1],id,"タイトル"));

conplotはifplotと同様に、与えられた2変数多項式関数に対して、その等高線表示を行う。つまり多項式関数f(x,y)が与えられると、 f(x,y)+level(デフォルト：-2<=level<=2、刻幅0.25)の ifplotを行う。levelはオプションで与える事が可能である。

例

[350] conplot(x^3+x^2-y^2);

[351] conplot(x^5-x*y^2+y^5,[x,-1,1],[y,-1,1],0,"Neko");

ここでウインドウの右側のスクロールバーをポインティングデバイスで上下に操作したり、バーを押したりすると、そのlevelに対応するグラフがハイライト表示されるので判別し易い。更に、右上のlevelにもその対応するlevelの値が表示される。

plot

plot(1変数関数f(x) (,[w1.w2],[x,x0,x1],id,"タイトル"));

plotは1変数関数のグラフの表示を行う。[w1,w2]は表示ウインドウの大きさを指定し、[x,x0,x1]で変域、idでサーバーのid、 "タイトル"で表示ウインドウのタイトルを指定する。

例

[355]  plot(sin(2*x)*exp(-x),[x,0,10]);

plotover

plotover(関数,id,"タイトル");

plotoverはplot,ifplot,conplotで予め表示したグラフに対し、新たに関数のグラフを書き加える場合に用いる。尚、与える関数は、plotとifplotやconplotで与えるものに合せる必要がある。つまり、plotで表示したグラフに対して追加する場合には、 plotoverでは1変数関数、ifplotやconplotの場合には2変数関数を表示する。

例

[360] ifplot(x^3-y^2+1);
0
[361] plotover(x^3-y^2,0,0);
0
[362]

例

下らない例を一つ。x^2*y^2+y^2*z^2+z^2*y^2-x*y*zの零点集合の平面x+y+z=kによる断面を表示する。この零点集合は次の様な曲面 (Steiner-Roman曲面)を構成している。

Steiner-Roman曲面(surfで作成(使い回し))

ifplotで、ベクトル[1,1,1]に直交し、[0,0,k]を含む平面x+y+z=k による上の曲面の断面を観察する。その為、z=k-x-yとして上記の式のzに代入したものを用いる。

関数定義

def f(K) {
Z=K-x-y;
P=x^2*y^2+y^2*Z^2+Z^2*x^2-x*y*Z;
return P;
}

尚、上記のプログラムではx,yを不定元とするが、Z,Pはこの様な不定元を含む式が代入される変数、即ち、プログラム変数である。risa/asir では、プログラム変数は大文字で開始する。

実行例

[368] for (I=-5;I<=10;I++)
ifplot(f(0.1*I),[x,-0.6,0.6],[y,-0.6,0.6]);
[369]

この処理を行うと、15枚のifplotによるグラフ表示を含むウインドウが続々と現われる。このグラフをepsファイルとして一々保存して、 ImageMagicのconvertでmpegに一気に変換するとアニメーションが楽しめる。但し、そのままでは白に黒なので面白く無い。そこで、epsファイルを弄ぶ事にする。その前にasirが生成するepsファイルの中身の抜粋を以下に示す。

出力ファイルの内容(抜粋)

%%BoundingBox: 100 100 400 400
%%Creator: This is generated by ifplot
%%Title: ifplot
%%EndComments:
0.1 setlinewidth
2 setlinecap
2 setlinejoin
/ifplot_putpixel {
    /yyy 2 1 roll def /xxx 2 1 roll def
    gsave newpath xxx yyy .5 0 360 arc
    fill grestore
} def
 137 198  ifplot_putpixel
 137 199  ifplot_putpixel
…略…
 323 344  ifplot_putpixel
 gsave [3] 0 setdash newpath
 249 100 moveto 249 400 lineto stroke
 stroke grestore
 gsave [3] 0 setdash newpath
 100 249 moveto 400 249 lineto stroke
 stroke grestore
showpage

ここでグラフ本体が} def以下のlifplot_putpixelが末尾にある行の個所で、 gsave以下にあるlineto strokeの個所がX,Y座標軸である。ここでPostScript ではsetrgbcolorで色の設定が行える。例えば、1 0 0 setrgbcolorで、それ以降の色は赤になる。その為、先頭にsetrgbcolorで赤を曲線に設定し、座標軸を黒とする。この為に以下の安易なスクリプトを用いる。

ifplotで生成したPSファイルのグラフの色を赤に変更するスクリプト

#! /bin/sh
cat $1| awk '{print $0}
             /} def/{print "1 0 0 setrgbcolor";}
             /gsave \[3\]/{print "0 0 0 setrgbcolor";}'>tmp;
mv tmp $1;

このスクリプトではawkを用いている。/パターン/{処理}でストリームにパターンが含まれていると処理を実行する。この場合、"} def"を含む行が入力ストリームに含まれていれば、その直後に"1 0 0 setrgbcolor" を含む行を追加し、gsave [3]が含まれていれば、"0 0 0 setrgbcolor"を含む行を追加する。こうする事で曲線のみを赤(1 0 0)とし、座標軸は黒(0 0 0)にする。このスクリプトによる変換例とMPEGファイルへの変換の応用例を以下に示す。尚、このスクリプトは兎に角動けば良いとの観点ででっちあげたもののであって、汎用性等は皆無である。

変換例

変換前

変換後

応用例

ponpoko@X20:~/asir > ls * | awk '/eps/{print "sh mkmk",$1}' >aa
ponpoko@X20:~/asir > sh aa; convert -density 100x100 *.eps asir.mpg

作成したMPEGファイル：asir.mpg

この様にして一気にmpegファイルに変換する。短か過ぎて白地が多い為に、ステップ毎に表示しないを意味不明かもしれない(少なくとも、X20では良く見えない)。こんな場合はgifアニメーションに変換した方が良いかもしれない。 Surf - Examples from SuSE Linux Manuals.に示したMPEGファイル作成の例と比べて派手さは無いが、これはこれで面白い。

尚、surfで生成したものを以下に示す。

Surfで生成した断面

この絵を生成するsurfのスクリプトも以下に示す。 surfに関しては Surf - Examples from SuSE Linux Manuals. を参照されたい。又、ImageMagicを用いてMPEGファイルを生成する話に関しては、上記のペイジと便利な話を参照されたい。

Surfのスクリプト

rot_x=-0.2;
rot_y=0.6;
rot_z=0.6;
origin_z=0;
scale_x=0.06;
scale_y=0.06;
scale_z=0.06;
illumination=ambient_light+
diffuse_light+reflected_light
+transmitted_light;
transparence=50;
clear_screen;
surface=x^2*y^2+x^2*z^2+y^2*z^2-x*y*z;
draw_surface;
plane=x+y+z-1;
curve_red=255;
curve_green=0;
curve_blue=0;
cut_with_plane;
plane=x+y+z-0.3;
curve_red=200;
curve_green=100;
curve_blue=0;
cut_with_plane;
plane=x+y+z;
curve_red=100;
curve_green=200;
curve_blue=100;
cut_with_plane;
plane=x+y+z+0.2;
curve_red=0;
curve_green=100;
curve_blue=255;
cut_with_plane;
plane=x+y+z+0.5;
curve_red=0;
curve_green=255;
curve_blue=0;
cut_with_plane;

GNUPLOTで遊ぶ話

RISA/ASIRはOpenXMを通して、gnuplotによる表示も可能である。この場合、最初にload("xm");でパッケージの読込を行い、それから、 gnuplot_start();を実行すれば、後はgnuplot命令で様々なグラフの表示が行える。以下にasirを立ち上げてからの様子を示す。

実行例

ponpoko@X24:~> asir
This is Risa/Asir, Version 20020301 (Kobe Distribution).
Copyright (C) 1994-2000, all rights reserved, FUJITSU LABORATORIES LIMITED.
Copyright 2000,2001, Risa/Asir committers, http://www.openxm.org/.
GC 5.3, copyright 1999, H-J. Boehm, A. J. Demers, Xerox, SGI, HP.
PARI 2.2.1(alpha), copyright (C) 2000,
     C. Batut, K. Belabas, D. Bernardi, H. Cohen and M. Olivier.
Asir-Contrib xm version 20010310. Copyright 2000-2001, OpenXM Committers.
help("keyword"); ox_help(0); ox_help("keyword"); ox_grep("keyword");
     for help messages. 
[716] load("xm");
1
Asir-Contrib xm version 20010310. Copyright 2000-2001, OpenXM Committers.
help("keyword"); ox_help(0); ox_help("keyword"); ox_grep("keyword");
     for help messages. 
[900] gnuplot_start();
0
[901]  gnuplot("plot sin(x**2);");
0
[902] gnuplot("splot (x**2+y**2)*sin(x**2+y**2)");
0
[903] gnuplot("set isosamples 50");
0
[904] gnuplot("splot (x**2+y**2)*sin(x**2+y**2)");
0
[905] gnuplot("splot x*y*sin(sqrt(x**2+y**2))");
0

実際にグラフ表示を行っているのが、gnuplot命令の個所である。 plotで曲線、splotで曲面を表示している。この上の例で一番下の x*y*sin(sqrt(x**2+y**2))のグラフを以下に示しておく。

このgnuplotの使い方は、OpenXM/lib/asir-contrib/gnulotのコメントの個所を参考にして、実際に動かしたものである。詳細はこのgnulotの個所を参照されたい。

povrayで遊ぶ話

上記の例では曲線や陰線処理を行った曲面表示で、surfの様な派手派手な曲面を描くものではない。Risa/Asirのパッケージの中には、povrayを用いて派手な曲面を描く為のデータを生成するものがある。但し、コメントを見る限りでは試験的なものの様である。このpovrayで遊ぶ為のパッケージはpovray.rrで、 z=f(x,y)の形式の曲面データの生成が行える。パッケージの読込は load("povray.rr");で行える。

先ず、povray.rrの中にある例題を以下に示す。

    povray_prolog("t.pov",0,0);
    povray_spheres(0,[ [0,3,0,3],[1,1,0,2] ]);
    povray_epilog(0);

最初のpovray_plologでpovray用のデータファイル名(ここではt.pov)を指定し、その指定したファイルを開く。次に、povray_spheresでリストで指定された二つの球面データを生成し、指定したファイルに書込み、最後のpovray_epilogでpovray_prologで開いたファイルを閉じている。これで生成したデータは指定ファイル(上の例ではt.pov)に書込まれるが、別にpovrayを起動してレンダリング結果が表示される訳では無い。この場合、単にpovray +It.povとすれば、t.pngファイルが生成される。以下にt.povをレンダリングした結果を示しておく。

上の例では自分でpovrayを起動する必要がある。このデータ生成に加え、 povrayを起動してレンダリングする事も可能である。この動作例は povray.rrのtest_mesh2等にある。こtest_mesh2では、レンダリングを povray_render("t.pov");で実行している。このtest_mesh2でレンダリングを行った結果も示しておく。

実際はもう少し簡単に曲面の表示も行える。例えば、povray_plot3d_mesh がある。この命令ではデータファアイルはt.povに固定され、XとYの領域も [-5,5]に固定される。更に、XY平面も同時にレンダリングされる為、多少 Z軸正方向に関数を上ておかないと隠れてしまう。以下に、povray_plot3d_mesh の入力と表示例を示す。

[744] povray_plot3d_mesh(x*y*sin((x^2+y^2)^(1/2))+25);
Z_bound=[13.6949,36.3051]
[[-5,-1],[5,5]]
 (t.pov) (w) file /asir_file_tmp set 

(出力データが長い為に省略)

@@@           @@@@@@@     ...........       ...@            @@@@@@@    @@@@@@
@@@@@@    @@@@@@@@@@@@        ....            ........    @@@@@@@@@@@@       
@@@@@         @@@@             @@@@.....    .............         @@@@       
@              @@@@@@@       @@@@@@@@@@@@.       .......              @@@@@@@
           @@@@@@@@@@@@@MM   MMMMMMMMMMM               @@@            @@@@@@@

Done Tracing
t.pov Statistics, Resolution 200 x 200
----------------------------------------------------------------------------
Pixels:           40000   Samples:           40000   Smpls/Pxl: 1.00
Rays:             67223   Saved:                 0   Max Level: 2/5
----------------------------------------------------------------------------
Ray->Shape Intersection          Tests       Succeeded  Percentage
----------------------------------------------------------------------------
Mesh                             97595            7394      7.58
Plane                            97595           31000     31.76
Sphere                           97595               0      0.00
Bounding Box                    571740          143115     25.03
----------------------------------------------------------------------------
Calls to Noise:                  0   Calls to DNoise:             10
----------------------------------------------------------------------------
Shadow Ray Tests:            91116   Succeeded:                 3569
Reflected Rays:              27223
----------------------------------------------------------------------------
Smallest Alloc:                 12 bytes   Largest:            26632
Peak memory used:           124484 bytes
----------------------------------------------------------------------------
Time For Trace:    0 hours  0 minutes   5.0 seconds (5 seconds)
    Total Time:    0 hours  0 minutes   5.0 seconds (5 seconds)

;
0
[746]

t.png

上の図で示す様に、安易にpovrayを実行すればレンダリングされる画像はややサムネイルに近い大きさのもの(200x200)である。しかし、生成したデータを用いて、povrayに画像の高さと幅を指定してやれば、幾らでも好きな大きさの画像が生成可能である。例えば、1000x700のサイズの画像が欲しければ、povray +It.png +W1000 +H700とすれば良い。勿論、大きなサイズを指定すればその分レンダリングに時間がかかるので、注意が必要になる。因に上の画像を押せば、1000x700の画像が表示される。

ここでpovray_renderはどの様な事をしているかと言えば、実際はshell命令でpovrayを立ち上げているだけである。以下にpovray_renderの定義内容を示しておく。

def povray_render(F)  {
   /* If FreeBSD, */
   shell("   povray +i"+F+" +L/usr/local/lib/povray3/include +D0 +P +H200 +W200&")$
}

ここで示す様に、povray_renderでは引数Fで指定したファイルに対し、 200x200でレンダリングを行うだけである。 povray.rrを参考にして自前のプログラムを作るのも楽しいだろう。

数式処理の目次に戻る

ponpoko

Last modified: Fri Jul 11 07:56:14 JST 2003