Pharoからgnuplotでsin(x)

UnifiedFFIを使ってC言語の関数を呼べたので、 Pharoからgnuplotを呼び出してみた。

呼び出す関数は、以前書いたgnuplotへのpipeを開いて コマンドを書き込む以下の関数。

/* gpdo.c */
#include <stdio.h>

FILE *gpstart(void);
void gpstop(FILE *p);
void gpdo(FILE *p, char *s);

FILE *gpstart(void) {
    /* start gnuplot process */
    return popen("gnuplot.exe --persist", "w");
}

void gpstop(FILE *p) {
    /* stop gnuplot process */
    gpdo(p, "quit");
    pclose(p);
}

void gpdo(FILE *p, char *s) {
    /* send command string to gnuplot process */
    fprintf(p, s);
    fprintf(p, "\n");
    fflush(p);
}

Windowsのgcc(MinGW)で共有ライブラリにします。

gcc -shared -o gpdo.so gpdo.c

出力されたgpdo.soをPharo.exeがあるフォルダに配置しました。

前と同様に、Pharoを起動してFFILibraryを継承したライブラリファイルを示す クラスGnuplotDoLibを作成しました。そのクラスに、win32ModuleNameメソッド を追加し、共有ライブラリのファイル名を返すようにします。

win32ModuleName
    ^ 'gpdo.so'

次は、関数呼び出し用のクラスを作ります。Objectクラスを継承してGnuplotDoクラスを 作りました。メソッドは関数を呼び出すだけなので、クラス側のメソッドにします。 これで、呼び出すときにわざわざインスタンスを生成する必要がなくなります。 三つの関数についてそれぞれメソッドを追加しました。 FILE構造体へのポインタ(FILE *)については、 今回はPharo側で解釈する必要がない（メンバーにアクセスしたりしない）ため、(void *)とします。

gpOpen
    ^ self ffiCall: #(void * gpstart ()) module: GnuplotDoLib.

gpClose: aHandle
    self ffiCall: #( void gpstop (void * aHandle)) module: GnuplotDoLib

gpHndl: aHandle gpCmd: aString
    self ffiCall: #( void gpdo (void * aHandle, String aString) ) module: GnuplotDoLib

これで準備完了です。Playgroundを開いて動かしてみます。

hndl := GnuplotDo gpOpen.
GnuplotDo gpHndl: hndl gpCmd: 'set xlabel "x"'.
GnuplotDo gpHndl: hndl gpCmd: 'set ylabel "y"'.
GnuplotDo gpHndl: hndl gpCmd: 'plot sin(x)'.
GnuplotDo gpClose: hndl.

とりあえず、これでPharoからgnuplotを起動してラベル設定してsin(x)を プロットする手続きが行えました。 エラーが起こるようなコマンドをgnuplotに渡すとgnuplotがクラッシュ してしまいますが、よしとします。

pythonでSTL形式のデータを作ってみた

gnuplotでは、以下のようなスクリプトで球面を描画したりします。

set parametric # 媒介変数モード
set urange[0:2*pi] # uの範囲を設定
set vrange[0:2*pi] # vの範囲を設定
set ticslevel 0
set hidden3d
set isosample 50
set view equal xyz
splot cos(u/2.0)*cos(v),sin(u/2.0)*cos(v),sin(v) w l

このようなものを3Dモデルとして出力したいということで、 pythonで(u, v)の関数をSTL形式にして出力するスクリプトを書いてみました。 ほかの言語で書き換えて遊ぶのにも、ちょうどよい内容と長さかと思います。

import math

def Sphere(u, v):
    x = math.cos(u/2.0)*math.cos(v)
    y = math.sin(u/2.0)*math.cos(v)
    z = math.sin(v)
    return x, y, z

def VectorSub(v1, v2):
    return [e1-e2 for e1, e2 in zip(v1, v2)] 

def VectorDiv(v, a):
    return [e/a for e in v] 

def VectorCross(v1, v2):
    x1, y1, z1 = v1
    x2, y2, z2 = v2
    return [y1*z2-z1*y2, z1*x2-x1*z2, x1*y2-y1*x2]

def VectorAbs(v):
    return math.sqrt(sum([e**2 for e in v]))

def NormalVector(p1, p2, p3):
    v1 = VectorSub(p2, p1)
    v2 = VectorSub(p3, p2)
    vc = VectorCross(v1, v2)
    a = VectorAbs(vc)
    v = VectorDiv(vc, a) if a!=0.0 else False
    return v

def STLfacet(apolygon, inverse=False):
    if inverse:
        points = apolygon[::-1]
    else:
        points = apolygon[:]

    nv = NormalVector(*points)

    if not nv:
        return ""

    l = "facet normal %g %g %g\n" % tuple(nv)
    l = l + "outer loop\n"
    for p in points:
        l = l + "vertex %g %g %g\n" % tuple(p)
    l = l + "endloop\n"
    l = l + "endfacet\n"
    return l

def PolygonsToSTL(name, polygons, inverse=False):
    l = "solid %s\n" % name
    for p in polygons:
        l = l + STLfacet(p, inverse)
    l = l+"endsolid\n"
    return l

def Polygons_of_Function_UV(func):
    nu = 90 
    nv = 90 
    du = 2.0 * math.pi / nu
    dv = 2.0 * math.pi / nv

    u = [du * i for i in range(nu+1)]
    v = [dv * i for i in range(nv+1)]

    polygons = []
    for i in range(nu):
        for j in range(nv):
            p1 = func(u[i], v[j])
            p2 = func(u[i+1], v[j])
            p3 = func(u[i+1], v[j+1])
            p4 = func(u[i], v[j+1])
            polygons.append((p1, p2, p4))
            polygons.append((p4, p2, p3))
    return polygons

def OutputSTL(name, polygons):
    with open('%s.stl' % name, "w") as fp:
        fp.write(PolygonsToSTL(name, polygons, False))

if __name__ == '__main__':
    OutputSTL("Sphere", Polygons_of_Function_UV(Sphere))

出力されたSTLデータを見るには、GLC Playerが便利ですね。また、Windows10には、3D Builderがついているので、そのままでSTLデータを見ることができます。GLC Playerでスナップショットを取ったものが下の図です。

関数を次のものにしてみます。

def Dounut(u, v):
    x = math.cos(u)*(1+0.5*math.cos(v))
    y = math.sin(u)*(1+0.5*math.cos(v))
    z = 0.5*math.sin(v)
    return x, y, z

ドーナツの形になりました。

def Helical(u, v):
    n = 10.0
    a = 0.6
    b = 0.3
    rm = (lambda v: a*b/math.sqrt(b*b*math.cos(v)**2 + a*a*math.sin(v)**2))
    r = (lambda u, v: 1 + rm(v)*math.cos(v+n*u/2.0))
    x = r(u, v) * math.cos(u)
    y = r(u, v) * math.sin(u)
    z = rm(v) * math.sin(v + n*u/2.0)
    return x, y, z

楕円をぐるぐる回転させていくと、ねじれた形状になりました。

Windowsのgnuplotで、plot sin(x)をしてみた

Javaは知らないけど、便利なライブラリがあるようなので、 Clojureを触ってみようということでやってみました。

数値を得たら、プロットしてみるということはしばしばあります。 pythonだとmatplotlibを使ったりしますが、 gnuplotを呼び出すことができれば、 グラフ描画をgnuplotにまかせることもできそうです。

gnuplotは、gp463-win32-setup.exeを使ってインストールしたも のを利用しています。(ver.4.6.5でも大丈夫です) コマンドラインでgnuplotが実行できるように、インストールされたgnuplot.exeがあるディレクトリへあらかじめPathを通しておきます。

Clojureは、1.6.0をダウンロードしました。 Getting Startedにあるように、

java -cp clojure-1.6.0.jar clojure.main

で、REPLを起動して利用しています。

以下が、そのスクリプトですが、 プロセスを開始するgpstart、gnuplotのコマンドを送るgpdo、 終了するgpstopという関数を作っています。 sin(x)をプロットするのは、demo関数で行います。

(defrecord GnuplotProc [proc out in])

(defn gpstart
  "Start gnuplot process"
  []
  (let [proc (.start (doto (ProcessBuilder. '("gnuplot" "--persist"))
                       (.redirectErrorStream true)))
        out (clojure.java.io/writer (.getOutputStream proc))
        in (clojure.java.io/reader (.getInputStream proc))]
    (GnuplotProc. proc out in)))

(defn gpstop
  "Stop gnuplot process"
  [proc]
  (.destroy (get proc :proc)))

(defn gpdo
  "send command string to gnuplot process"
  [proc s]
  (let [w (get proc :out)]
       (.write w (str s))
       (.newLine w)
       (.flush w)))

(defn demo
  "Demo function"
  []
  (def p (gpstart))
  (gpdo p "set xlabel \"x\"")
  (gpdo p "set ylabel \"y\"")
  (gpdo p "plot sin(x)")
  (read)
  (gpstop p))

(demo)

clispと、pythonでも書いてみました。

(defun gpstart 
  ;Start gnuplot process
  ()
  (make-pipe-io-stream "gnuplot --persist" :buffered t))

(defun gpstop
  ;Stop gnuplot process
  (proc)
  (let ()
    (gpdo proc "quit")
    (close proc)))

(defun gpdo
  ;send command string to gnuplot process
  (proc s)
  (let ()
    (format proc "~A~%" s)
    (force-output proc)))

(defun demo
  ;demo function
  ()
  (let ((p (gpstart)))
    (gpdo p "set xlabel \"x\"")
    (gpdo p "set ylabel \"y\"")
    (gpdo p "plot sin(x)")
    (read)
    (gpstop p)))

(demo)

Pythonで書いたスクリプトは以下の通り。

import subprocess

def gpstart():
    """ Start gnuplot process """
    proc = subprocess.Popen(["gnuplot", "--persist"], 
                            stdin = subprocess.PIPE,
                            stdout = subprocess.PIPE,
                            stderr = subprocess.PIPE)
    return proc

def gpstop(proc):
    """ Stop gnuplot process """
    gpdo(proc, "quit")
    proc.kill()

def gpdo(proc, s):
    """ send command string to gnuplot process """
    proc.stdin.write(s)
    proc.stdin.write("\n")
    proc.stdin.flush()

def demo():
    """ Demo function """
    p = gpstart()
    gpdo(p, 'set xlabel "x"')
    gpdo(p, 'set ylabel "y"')
    gpdo(p, 'plot sin(x)')
    raw_input()
    gpstop(p)

if __name__ == '__main__':
    demo()

この程度なら、行数はほぼ一緒ですね。
clispやClojureを使ってもgnuplotを操作できそうなので、ひと安心です。

Rubyでもやってみました。

def gpstart
    # Start gnuplot process
    IO.popen("gnuplot --persist", "r+")
end

def gpstop(io)
    # Stop gnuplot process
    io.close
end

def gpdo(io, s)
    # send command string to gnuplot process
    io.puts(s)
end

def demo
    io = gpstart
    gpdo(io, "set xlabel \"x\"")
    gpdo(io, "set ylabel \"y\"")
    gpdo(io, "plot sin(x)")
    gets
    gpstop(io)
end

if __FILE__ == $0
    demo
end

シンプルですね。
Windows環境でirbを利用する人は少ないのかな。 ActiveScriptRuby 2.1.2-p95を利用してみたのですが、 irbのインタラクティブな環境で、Backspaceがちゃんと動作しませんでした。 検索してみたら、1.9.2での同様な症状と対処法がこちらにありました。 irbに--noreadlineオプションをつけて起動すると、正常にBackspaceできるようになりました。


2016/07/24追記
同じものをC言語で。Windows Vista 32bit MinGWのgcc-4.9.3で確認。

#include <stdio.h>

FILE *gpstart(void);
void gpstop(FILE *p);
void gpdo(FILE *p, char *s);

FILE *gpstart(void) {
    /* start gnuplot process */
    return popen("gnuplot.exe --persist", "w");
}

void gpstop(FILE *p) {
    /* stop gnuplot process */
    gpdo(p, "quit");
    pclose(p);
}

void gpdo(FILE *p, char *s) {
    /* send command string to gnuplot process */
    fprintf(p, s);
    fprintf(p, "\n");
}

void demo(void) {
    /* demo function */
    FILE *p;
    p = gpstart();
    gpdo(p, "set xlabel \"x\"");
    gpdo(p, "set ylabel \"y\"");
    gpdo(p, "plot sin(x)");
    gpstop(p);
}

void main(void) {
    demo();
}

カラムチェッカー

pythonで書いたスクリプトの追加です。
グラフを描画するソフトはいろいろありますが、その中でもgnuplotはおすすめです。gnuplotは多くのプラットホーム上で動作し、誰でも利用できるからです。また、ちょっとデータをグラフにして確認したいときに、コマンド入力による軽快な操作でプロットすることができます（もちろんコマンドを知っていないといけませんが）。ただ、データをプロットするときに、どのカラム（列）のデータをプロットするのかを数字で指定しなくてはなりません。カラム番号があらかじめ分かっていればよいのですが、多くのカラムがある場合、プロットしたいデータが何番目のカラムにあるのかを知るのは大変です。
そのようなときに便利なスクリプトが、今回追加したcColumn.pyです。
ここからダウンロードできます。

Windows上で利用するために作成しましたが、ubuntu10.04のPC上でも動作しました。
python2.6以降とwxPython2.8以降があれば動くのではないでしょうか。