haskell + gnuplot + julia

Pretty Nine Zulu  Queens Ruled China.

美しい９人のズールー族の女王が中国を治めた

Pは素数、Nは自然数、Zは整数、Qは有理数、Rは実数、そして最後は複素数。数の包含関係を 覚えるのによいとの事。

素数は数の元素ってのが、こういう風に言われると激納得。自然数からスタートしてるような 本は疑った方がいいぞ。まてまて、元素が無限に有って良いのか？ はなはだ疑問。 科学界では、人口元素を作るのが流行っているけど、無限種類の元素って、理論的に作成 可能なの？

突然こんな事を書いたのは、 『算数宇宙の冒険　アリスメトリック！』（実業之日本社） なんてのを読んだからだ。この本を読んで、Rで数学をやってみる まで、足を延ばしておられるな。

上記の数の包含関係をきっちり実現してる言語って何があるだろう？ Schemeがまず頭に浮かぶな。 続いてHaskellかな。数学に強いと言うJuliaはどうかな？

さ迷っていたら、 Euler Math Toolboxなんてのに遭遇。このページの 所の Demo - Plots in Euler に面白いグラフが有るぞ。（はあと）

Haskellとグラフ

色々やってもしょうがないんで、取り合えずHaskellです。Haskellでグラフしたい。さて どうする。

そりゃ、Gnuplotと組み合わせるに決まってる。 HaskellからGnuplotを呼び出してチャートを描く とか、 gnuplotに目覚めた方とか。

でも、Gnuplot以外の選択は？ Haskell Chart パッケージでグラフ描画 なんて方もおられます。多様性は善。

取りあえずThe gnuplot packageに行って デモ付きでcabalしてみました。そして、試運転。

Loading package array-0.5.0.0 ... linking ... done.
Loading package deepseq-1.3.0.2 ... linking ... done.
Loading package containers-0.5.5.1 ... linking ... done.
Loading package stm-2.4.2 ... linking ... done.
Loading package filepath-1.3.0.2 ... linking ... done.
Loading package transformers-0.3.0.0 ... linking ... done.
Loading package mtl-2.1.3.1 ... linking ... done.
Loading package data-accessor-0.2.2.7 ... linking ... done.
Loading package data-accessor-transformers-0.2.1.7 ... linking ... done.
Loading package time-1.5.0.1 ... linking ... done.
Loading package bytestring-0.10.4.0 ... linking ... done.
Loading package unix-2.7.1.0 ... linking ... done.
Loading package directory-1.2.5.0 ... linking ... done.
Loading package process-1.2.3.0 ... linking ... done.
Loading package pretty-1.1.1.1 ... linking ... done.
Loading package template-haskell ... linking ... done.
Loading package transformers-compat-0.5.1.4 ... linking ... done.
Loading package exceptions-0.8.2.1 ... linking ... done.
Loading package temporary-1.2.0.3 ... linking ... done.
Loading package utility-ht-0.0.11 ... linking ... done.
Loading package gnuplot-0.5.4 ... linking ... done.

色々な物が呼び出されるのは、デモの為？ .cabal/bin/gnuplot-demoで美しい例が 観賞出来る。

今度は、例のやつをrunghcすると、

[sakae@fedora plt]$ runghc bb.hs

bb.hs:4:6:
    No instance for (Fractional a0) arising from a use of ‘linearScale’
    The type variable ‘a0’ is ambiguous
    Relevant bindings include xs :: [a0] (bound at bb.hs:4:1)
    Note: there are several potential instances:
      instance Fractional Double -- Defined in ‘GHC.Float’
      instance Fractional Float -- Defined in ‘GHC.Float’
      instance Integral a => Fractional (GHC.Real.Ratio a)
        -- Defined in ‘GHC.Real’
      ...plus one other
    In the expression: linearScale 1000 (0, 2 * pi)
    In an equation for ‘xs’: xs = linearScale 1000 (0, 2 * pi)

bb.hs:4:24:
    No instance for (Num a0) arising from the literal ‘0’
    The type variable ‘a0’ is ambiguous
    Relevant bindings include xs :: [a0] (bound at bb.hs:4:1)
    Note: there are several potential instances:
      instance Num Double -- Defined in ‘GHC.Float’
      instance Num Float -- Defined in ‘GHC.Float’
      instance Integral a => Num (GHC.Real.Ratio a)
        -- Defined in ‘GHC.Real’
      ...plus 12 others
    In the expression: 0
    In the second argument of ‘linearScale’, namely ‘(0, 2 * pi)’
    In the expression: linearScale 1000 (0, 2 * pi)
     :

こんな具合に盛大なエラーですよ。ghciからはちゃんと動くのにね。何と言う事で しょう。ソース嫁って事で、ご対面。Utility.hsに目当ての物が有りました。

linearScale :: Fractional a => Integer -> (a,a) -> [a]
linearScale n (x0,x1) =
   map (\m -> x0 + (x1-x0) * fromIntegral m / fromIntegral n) [0..n]

始点と終点がタプルで与えられて、それをｎ等分した等差数列を返す関数なんだな。 等比数列とかログスケールは有るのかな、なんて思っちゃうぞ。

この定義を見て、はたと気が付いた。タプル内の初項が、何とでも解釈出来ちゃうから 厳格な数学者が怒り心頭で文句を言ってきたのね。

import Graphics.Gnuplot.Simple

xs = linearScale 1000 (0::Float, 2 * pi)
ys = fmap sin xs
points = zip xs ys
main = do
        plotPath [(Title "hello")] points

これで、ちゃんと動いた。数学者は、こうしたら良いよって、ちゃんと英語でサゼッション してくれていたのに、無視してたオイラーは島国住人。ちょいと恥ずかしいぞ。

import Graphics.Gnuplot.Simple

xs = linearScale 1000 ((-2::Float)*pi, 2*pi)
ys = map (\x -> sin(x)/x) xs
points = zip xs ys
main = do
        plotPath [(Title "hello")] points

要領が分かれば、ちょいと変形も可能。学ぶは真似るから来ているらしいですから。

import Graphics.Gnuplot.Simple

pnts :: [Float]
pnts = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0]
doPlot = plotPathStyle [ ( Title "plotting dots" )]
    (PlotStyle Points (CustomStyle []))  (zip [0..] pnts)

main = do
    doPlot

ドットスタイルの例。

import Graphics.Gnuplot.Simple 
Prelude Graphics.Gnuplot.Simple> plot
plotDots        plotFuncs       plotLists       plotParamFunc   plotPathStyle   plotType
plotFunc        plotList        plotListsStyle  plotParamFuncs  plotPaths
plotFunc3d      plotListStyle   plotMesh3d      plotPath        plotPathsStyle

plotまで叩いて、続いてTABキーを叩くと、ghciが補間してくれて、どんなのが有るか 教えてくれます。詳細を知りたい場合は

Prelude Graphics.Gnuplot.Simple> :t plotPath
plotPath
  :: Graphics.Gnuplot.Value.Tuple.C a =>
     [Attribute] -> [(a, a)] -> IO ()

Prelude Graphics.Gnuplot.Simple> :t plotPathStyle
plotPathStyle
  :: Graphics.Gnuplot.Value.Tuple.C a =>
     [Attribute] -> PlotStyle -> [(a, a)] -> IO ()

こんな具合に調べられます。

Prelude Graphics.Gnuplot.Simple> :t plotFunc
plotFunc
  :: (Graphics.Gnuplot.Value.Atom.C a,
      Graphics.Gnuplot.Value.Tuple.C a) =>
     [Attribute] -> [a] -> (a -> a) -> IO ()
Prelude Graphics.Gnuplot.Simple> plotFunc [] [0,0.01..20] sin

これで、簡単にサイン波を描画出来ましたよ。

Attributeとかの調べ方

上でtypeを調べた時、AttributeとかPlotStyleとかが出てきた。マニュアルも無いし、 不親切だと思うぞ。大体、どのパッケージにもマニュアルが付いているのを見た事が無い。 それで文句が来ないのは、何処かに抜け道が有るに違いない。

そりゃ、ソースが有るんだから読め。以上、説明終了。数学者らしい。

諦めきれないオイラーは、ghciの使えそうなコマンドをいろいろ試してみたよ。 そして、分かったのは、

Prelude Graphics.Gnuplot.Simple> :browse!
  :
data Attribute
  = Custom String [String]
  | EPS FilePath
  | PNG FilePath
  | Terminal gnuplot-0.5.4:Graphics.Gnuplot.Private.Terminal.T
  | Grid (Maybe [String])
  | Key (Maybe [String])
  | Border (Maybe [String])
  | XTicks (Maybe [String])
  | YTicks (Maybe [String])
  | Size Size
  | Aspect Aspect
  | BoxAspect Aspect
  | LineStyle Int [LineAttr]
  | Title String
  | XLabel String
  | YLabel String
  | ZLabel String
  | XRange (Double, Double)
  | YRange (Double, Double)
  | ZRange (Double, Double)
  | Palette [(Double, (Double, Double, Double))]
  | ColorBox (Maybe [String])
  | XTime
  | XFormat String
data PlotType
  = Lines
  | Points
  | LinesPoints
  | Impulses
  | Dots
  | Steps
  | FSteps
  | HiSteps
  | ErrorBars
  | XErrorBars
  | YErrorBars
  | XYErrorBars
  | ErrorLines
  | XErrorLines
  | YErrorLines
  | XYErrorLines
  | Boxes
  | FilledCurves
  | BoxErrorBars
  | BoxXYErrorBars
  | FinanceBars
  | CandleSticks
  | Vectors
  | PM3d

ダーと出てくる表示から眼grepして抜き出すっていう特殊技能を身に付けるのが一番。 そんなの無理だって人は、擬似grepもどきの

Prelude Graphics.Gnuplot.Simple> :info Attribute
data Attribute
  = Custom String [String]
  | EPS FilePath
  | PNG FilePath
  | Terminal gnuplot-0.5.4:Graphics.Gnuplot.Private.Terminal.T
  | Grid (Maybe [String])
  | Key (Maybe [String])
  | Border (Maybe [String])
  | XTicks (Maybe [String])
  | YTicks (Maybe [String])
  | Size Size
  | Aspect Aspect
  | BoxAspect Aspect
  | LineStyle Int [LineAttr]
  | Title String
  | XLabel String
  | YLabel String
  | ZLabel String
  | XRange (Double, Double)
  | YRange (Double, Double)
  | ZRange (Double, Double)
  | Palette [(Double, (Double, Double, Double))]
  | ColorBox (Maybe [String])
  | XTime
  | XFormat String
        -- Defined in ‘Graphics.Gnuplot.Simple’

定義場所の教えてくれていて、便利だよ。後はgnuplotをどれぐらい知ってるかだな。

参考までに該当箇所のソースを当ると

data Attribute =
     Custom String [String]  -- ^ anything that is allowed after gnuplot's @set@ command
   | EPS    FilePath
   | PNG    FilePath
   | Terminal Terminal.T     -- ^ you cannot use this, call 'terminal' instead
   | Grid   (Maybe [String])
:

さすがに、コメントまではghciからは見られない仕様（コンパイル時に削除）になってました。 ソースを読むと面白い。

data Aspect =
     Ratio Double
   | NoRatio

{- |
Be careful with 'LineTitle'
which can only be used as part of gnuplot's @plot@ command
but not as part of @set@.
  :

冒頭には、ghciでもhugsでも動きますって書いてあったぞ。本当にhugsで動くかやってミレー。 やってみたけど、hugsに無いモジュール多数。そのモジュール内で使ってるものをコピペ移植 を始めたけど、きりがないので断念。根性あれば、動くだろうよ。

プラス アルファ

科学技術計算たら、行列が出てくる事になっている。世のスパコンも行列を高速で計算出来る ような作りになってる。

Haskellで行列、それにグラフを合わせれば、鬼に金棒。juliaといい勝負出来るぞ。 haskellからgnuplotでグラフ出力してみた なんてのにも遭遇。かっこいい。オライーもやってみたい。

Numeric.LinearAlgebra なんてのをインポートしてるけど、どんなパッケージを入れたら よいの？ そういう時は、cabalに聞くなり、Hayoo! に聞くなりします。

hmatrixを入れればよいと分かったんで開始。有名な土台が足りないと言われたぞ。 fedoraのdnfで関係者を入れた。

インストール:
  blas.i686 3.5.0-12.fc23              blas-devel.i686 3.5.0-12.fc23
  gcc-gfortran.i686 5.3.1-2.fc23       lapack.i686 3.5.0-12.fc23
  lapack-devel.i686 3.5.0-12.fc23      libgfortran.i686 5.3.1-2.fc23
  libquadmath.i686 5.3.1-2.fc23        libquadmath-devel.i686 5.3.1-2.fc23

これって、スパコンの京にも入っているのだろうか。fortranは偉大なり。

juliaでグラフ

暫くぶりにjuliaでもと思って、サンプルを探してきた。 Julia By Example いろいろな言語のいいとこ取りしてるそうなんで、すらすらと読める。途中からRっぽく なってきて、グラフが色々と紹介されてた。

色々なグラフはパッケージから入れるんだけど、要約がまとまっていた。 Graphics in Julia ありがたい事です。そして、ぐぐる様のチャートも使えるようです。 Juliaできれいなグラフを描こう

取りあえず、Windows7に安定版のjulia 0.4.3を入れました。IDE付きは大きくて重いので 入れてはみたものの、素のjuliaにしちゃいました。(後で気が付いたけど、素のやつでも、CPU リソース喰いまくり。結局、ごみ箱へ棄てました。)

本家のページで紹介されてた、Gadflyは結果をブラウザーに表示するって事で期待したんですが、 結果が表示されるまでの待ち時間が長すぎていらいら。

Winstonは、表示をTkに振っているようで、そこそこの使い心地。但し、色々なpkgを 併用するので、どうかなと。その点、gnuplotとの橋渡しをするだけのGastonは、余計な ものが入らずにすっきりです。

julia> Pkg.dependents("Winston")
5-element Array{AbstractString,1}:
 "ImageView"
 "EEG"
 "OpenStreetMap"
 "ParallelAccelerator"
 "RobustStats"

julia> Pkg.dependents("Gaston")
0-element Array{AbstractString,1}

Gaston

ASCIIPlots

ちょいとあらましだけを表示するならASCIIPlotsが良いかと。対抗馬でTextPlotも あるようですが、遠慮です。

fedora23にjuliaを入れて、Gastonを使おうとすると、

WARNING: Base.String is deprecated, use AbstractString instead.
  likely near /home/sakae/.julia/v0.4/Gaston/src/gaston_types.jl:117
ERROR: LoadError: LoadError: UndefVarError: Range1 not defined
while loading /home/sakae/.julia/v0.4/Gaston/src/gaston_types.jl, in expression starting on line 161
while loading /home/sakae/.julia/v0.4/Gaston/src/Gaston.jl, in expression starting on line 26

なにやらRange1ってのが未定義のようでエラーです。該当ファイルを開くとRangeとRange1が 使われているけど、Rangeも定義無し。この不平等さは何よ。それにRange1が使われている 関数も気になる。

こういう時は、replから聞いてみます。replの冒頭で？って叩くと、help?>に なるので、そこで質問。replには、もう一つshellを呼び出す ； が用意されてて、なかなか お洒落。

help?> Range
search: Range range RangeIndex nzrange linrange UnitRange StepRange histrange

  No documentation found.

  Summary:

  abstract Range{T} <: AbstractArray{T,1}

  Subtypes:

  FloatRange{T<:AbstractFloat}
  LinSpace{T<:AbstractFloat}
  OrdinalRange{T,S}

help?> isa
search: isa isascii isalpha isalnum isapprox isabspath isassigned disable_sigint

  isa(x, type) -> Bool

  Determine whether x is of the given type.

Quick hackで、Range1をRangeに変更。（gaston_midlvl.jlとgaston_types.jl）

これで、どうやら、パッケージをロード出来た。いざ実行しようとすると、

julia> histogram(y,"bins",25,"norm",1,"color","blue","title",
           "Rayleigh density (25 bins)")
WARNING: [a,b] concatenation is deprecated; use [a;b] instead
 in depwarn at deprecated.jl:73
 in vect at abstractarray.jl:38
 in figure at /home/sakae/.julia/v0.4/Gaston/src/gaston_hilvl.jl:99
 in histogram at /home/sakae/.julia/v0.4/Gaston/src/gaston_hilvl.jl:250
while loading no file, in expression starting on line 0

gnuplot> set term wxt 1
                  ^
         line 0: unknown or ambiguous terminal type; type just 'set terminal' for a list

今度は、ターミナルタイプが違うって、fedoraのgnuplotはデフォで qt ってなってたので、 gaston_types.jlとgaston_aux.jlで該当箇所を書き換え。これで、ワーニングは出るものの、 グラフを表示してくれた。やれやれ。そしてhackの御褒美は gaston_demo() これ、圧巻 だぞ。

using Gaston

t = 0:0.0001:.15
carrier = cos(2pi*t*200)
modulator = 0.7+0.5*cos(2pi*t*15)
am = carrier.*modulator
plot(t,am,"color","blue","legend","AM DSB-SC","linewidth",1.5,
    t,modulator,"color","black","legend","Envelope",
    t,-modulator,"color","black","title","AM DSB-SC example",
    "xlabel","Time (s)","ylabel","Amplitude",
    "box","horizontal top left")
set_filename("plotex.pdf")
set_print_linewidth(2)
set_print_fontsize(12)
set_print_size("6.5in,3.9in")
printfigure()

上図を実行した時の変調度を求めよって、酷試にでるぞ。

最近、友人が、50Mhz帯AM波を出すんだとか言って、 FT-991シリーズを衝動買い したって言ってた。きっかけは、昔の無線仲間が退職して暇になったので、遊びにきたからとか。昔の電波少年が、 電波老人になった訳ね。頑張ってラグチューしてください。

ボタンが沢山ある上、長押しとかダブルタップとか、覚えるの大変との事。 そりゃ、中身はラズパイで、周辺装置に、送信機と受信機が付属したものですから当然の 帰結。