plot
図書館から、『数学ガールの秘密ノート 式とグラフ編』なんて本を借りてきた。 著者はもちろんあの人。懐かしくなって手にした次第。それもあるんだけど、前回から やってるmaximaやgnuplotでグラフを書くのに触発された次第。
x^3 - x x^4 - 5x^2 + 4
等の、高次関数で、山（谷）が幾つあるかは、オイラー先生が盛んに研究してたとか。 眼の付け所が違うな。
ああ、眼の付け所と言えば、アマゾンの親分さんは、双曲線のグラフを見て、ロングテールに 商機を見出したんでしょうか。反比例のグラフね。
y = a / x
で、aは何になるでしょう？ と言う問い、恥ずかしながら答えられなかった。中学生にも 劣るな。
もう一冊、数学と言うか物理の本を借りてきた。『今日から使える量子工学』。ここに出て来る、 シュレーディンガー方程式ってのに惹かれたから。これ、よく名前を聞く。難しい、チンプンカンプンと言う代物。怖いもの見たさ、ミーハーの極み。
大体、今日から使えるって、そんな日常的なもの？ まあ、量子コンピューターが当たり前に 使えるようになれば、今日から使えるってのも頷けるけど、、、いや頷けるないか。今日使ってる普通のコンピュータでも、量子力学が応用されてるぞ。
例えば、オイラーのパソコンにはSSDが搭載されてる。これって、フラッシュメモリーの応用だよね。絶縁体の 中に電子を閉じ込めるか閉じ込めないかで、電子の流れを邪魔させたり（回路がOFFの状態）、 通したり（回路がON）出来るようになってる。いわゆるROMなんだけど、電子を絶縁体の中に 注入するって、量子力学のトンネル効果の応用だ。
不思議な世界。大体絶縁体の中にどうやって電子を注入したり、電子が入っていない状態を 作り出すのか？ フェルミ粒子とかボーズ粒子とかいう、とっても細かい粒子が動き回って いる世界が相手。やはり、凡人には意味不な世界でしたよ。
いやいや、一つ分かった事がある。物理の世界で色々実験したり考察したりして、仮定を立てる。それを式（微分方程式）で表す。後はそれを数学の世界に持って行って解く。その結果を 物理の世界に引き戻して、検討修正してく。この繰り返し。だから、数学は物理にとっては 単なる道具（ちと、言い過ぎか）。
同じ事が、数学ガールにも書かれていた。簡単なつるかめ算を解くのに、式を立て、数学の世界で連立方程式を解く。結果を実世界に戻して、答えとする。その方法が、算数とは違うんだよとな。
ともかく、シュレーディンガー方程式では討ち死にしたので、関連なトピックスを挙げておく。
plot2d, plot3d
上で出て来た、反比例のグラフ、双曲線を描いてみる。ああ、maximaが勝手に書くから、 オイラーは指示するだけだけど。
(%i1) plot2d(1/x, [x,-3,3]);
ｘの範囲は、-3から3までと、適当に決めた。そうすると分母がZEROになる時が有って、ZEROに よる割り算って事で、落ちないかひやひや。落ちなかった。でもY軸の範囲が勝手に選ばれて、 ±800になってたぞ。オートスケールで、気にしなくても良いので、楽でええわ。
気になるZERO近辺のデータをログから調べておくか。
-0.0323275862068965469e-1 -309.33333333333337 -0.1616379310344827347e-2 -618.6666666666667 0.1616379310344827347e-2 618.6666666666667 0.0323275862068965469e-1 309.33333333333337
見事に避けて通ってますなあ。これって、maximaが配慮してくれてるんだよね。あんた、数式の 事良く分かってるじゃん。
(%i2) plot2d(1/x, [x,-3,3], [y,-5,5]);
Y軸の範囲も、上記のようにすると、決めておける。この場合は、クリッピングの処理をgnuplot側でやってるのかな？ 一応確認。
-0.20366379310344823e+0 -4.910052910052912 -0.2004310344827586e+0 -4.989247311827958 0.2068965517241379e+0 4.833333333333334 0.21012931034482757e+0 4.758974358974359
maximaが、枠内に収まる範囲だけを計算してるよ。上からの類推がオイラーの頭に働かないのは 、回転が鈍ってきた証拠だな。
(%i3) plot2d([1/x, 2/x], [x,-3,3], [y,-5,5]);
複数のグラフを同時に描く場合は、リスト中に指定するようだ。
(%i4) plot2d(1/x, [x,-3,3], [y,-5,5]); plot2d: expression evaluates to non-numeric value somewhere in plotting range. plot2d: some values were clipped. (%o4) [/tmp/maxout76515.gnuplot_pipes]
ああ、メッセージを良く見たら、幾つかのデータはクリッピングしたって、親切に教えて くれていたよ。まったくもう、モウロクしてきたな。
(%i5) plot3d([cos(x)*(3+y*cos(x/2)), sin(x)*(3+y*cos(x/2)), y*sin(x/2)],
[x,-%pi,%pi], [y,-1,1], ['grid,50,15]) ;
これ、メビウスの輪らしいです。色が付いてるけど、勝手につくのかな。
ちょっとグラフ描きを調べてみたけど、まあこれぐらい分かれば困る事は無さそうだな。 それより、divide by zero のエラーが出るかの方が心配なんで、実験しておく。
(%i11) 1 / 0 ; expt: undefined: 0 to a negative exponent. -- an error. To debug this try: debugmode(true); (%i12) 1 / 1e-300 ; (%o12) 9.999999999999999e+299 (%i13) 1 / 1e-1000 ; expt: undefined: 0 to a negative exponent. -- an error. To debug this try: debugmode(true);
式の結果が未定義なら理解出来るけど、ZEROがダメとはねぇ。浮動小数で、うんと小さい数を 設定してみると、何処かに極限が有りそうな応答でしたよ。まあ、普通にプランク定数ぐらいは 扱えますよって事にしておこう。(6.62607004e-34 Js)
ここで、ふと疑問が出てきた。素のgnuplotを使って、上記のような 1 / 0 に遭遇する ようなグラフを書かせたら、どうなるかっての。（ええ、重箱の隅は重々承知してますよ） それには、プロットされたデータをファイルに落とせればいいんだな。果たして、そんな おもてなしをgnuplotがしてくれるか、それが問題だ。
これで調べろとな。set tabke "hoge.txt" とかやってから、グラフを書く。後はunset tableしとけとな。
gnuplot> plot [0:0.1] 1/x
に対して
# Curve 0 of 1, 100 points
# Curve title: "1/x"
# x y type
-6.67761e+153 -6.67761e+153 u
0.0010101 990 i
0.0020202 495 i
:
こんなデータが保存された。XはZEROに近い何かって事で、-6.67761e+153が選ばれているけど、 指数の極性間違ってる（よね）。これじゃ、とてつもなく大きな負数じゃん。
Debian9に入っている、GNUPLOT Version 5.0 patchlevel 5 では
0 0 u
0.010101 99 i
0.020202 49.5 i
:
これが正解。xrangeの初期値もステップ値も満足する値を刻んでいる。上のおかしなやつは、 OpenBSDに入っていたやつなんだけど(Version 4.6 patchlevel 6)、ZERO点意外でも、ステップ値がおかしいね。あれ？ 日を置いて、同じ事をやったら今度は正常になったぞ。イニシャライズ関係がおかしいのかね。
draw for haskell
色々あるけど、面白そうな『関数画家』を選んでみた。昔の事なので、使われていたのは、runhugsだった。和田先生に敬意を表して、hilbertを試してみた。お絵描きはポストスクリプトなのね。 何たるかを知らないので、ちと学習。レーザープリンターとかに使われている言語みたい。
これで、概要は分かった。紙に描くにはポストスクリプトを使う。WebのキャンパスにはSVGを 使うって事でよいかな。
それはそうと、現代の主流、ghcでも動くかしら？ 試してみる。
$ ghci hilbert.hs GHCi, version 8.2.1: http://www.haskell.org/ghc/ :? for help [1 of 1] Compiling Main ( hilbert.hs, interpreted ) hilbert.hs:12:4: error: Parse error in pattern: n + 1 | 12 | x (n+1) (x0,y0) = do y n (x1,y1); f (x1,y1); x n (x0,y0); f (x0,y0) | ^^^ Failed, 0 modules loaded.
見事にエラーだ。でも、この原因知ってるぞ。パターンマッチで、n + 1 と言う未来には マッチするのを禁止にしようという運動が有り、ghcでは使えないようになってしまったのだ。 和田先生みたいに、甘いも辛いも分かる方が使う分には、記事にも説明が有るように便利 だと思うんだけどなー。抜け穴スイッチでも用意されてるのかな？ 調べるのも面倒なので、 現代風に書き換えてみた。
import Numeric
left,right :: (Int,Int) -> (Int,Int)
left (dx,dy) = (-dy, dx) -- 左折
right (dx,dy) = (dy, -dx) -- 右折
f :: (Int,Int) -> IO () -- Fに対応するrlinetoは相対位置まで線を引く
f (dx,dy) = putStrLn ((show dx) ++ " " ++ (show dy) ++ " rlineto")
z :: Num a => a -> a
z = flip (-) 1
x,y :: Int -> (Int,Int) -> IO () -- X,Yに対応
x 0 (_,_) = putStr "" -- 0次ならなにもしない. 空を出力
x n (x0,y0) = do y (z n) (x1,y1); f (x1,y1); x (z n) (x0,y0); f (x0,y0)
x (z n) (x0,y0); f (x3,y3); y (z n) (x3,y3)
where (x1,y1) = left (x0,y0); (x3,y3) = right (x0,y0)
y 0 (_,_) = putStr ""
y n (x0,y0) =do x (z n) (x3,y3); f (x3,y3); y (z n) (x0,y0); f (x0,y0)
y (z n) (x0,y0); f (x1,y1); x (z n) (x1,y1)
where (x1,y1) = left (x0,y0); (x3,y3) = right (x0,y0)
hilbert :: Int -> Int -> IO () -- ドライバ
hilbert size n = do putStrLn ((show o) ++ " " ++ (show o) ++ " moveto")
x n (x0,y0)
putStrLn "stroke"
where x0 = size `div` (2 ^ n) -- セグメント長 (= 移動距離)
y0 = 0
o = x0 `div` 2 -- 原点から出発位置までの距離
main = hilbert 256 5
-- Usage: runghc hilbert.hs | gv -
n+1のパターンを止めて、nにし、本体側で、n-1 するように、関数zを定義した。あちこちに、 この関数をふりまかなければならず、非常な面倒さが実感出来るよ。
$ runghc hilbert.hs 4 4 moveto 0 8 rlineto 8 0 rlineto 0 -8 rlineto : 8 0 rlineto 0 -8 rlineto stroke
調子に乗ってfish.hsもghcで動かそうとしたら、上記のエラーの他に、rotが無いとかPainterが無いとか、色々言われた。必要な者は、functionalprogrammer.hsと squarelimit.hsに定義されてるので、適時持ってくればよい。
このfish.hsが吐き出す描画情報は、12948行にもなった。印刷具合確認ソフトも、これだけの 線画を受け取って、右往左往しながら表示してたよ。結構太い線が使われているので、細かい 所は、線が潰れていた。線を細く指定出来るのかな？
postscript
思わぬ所からpostscriptが出てきた。でも、今ならWebでSVGでしょ。その昔に C to GOなんてタイトルでやってた。
和田先生のps出力をSVG用に書き換えてみようかと思ったけど、ごちゃごちゃした記法に嫌気が さしてきたぞ。SVGと言えども、HTMLの一種。説明的なタグを付けるのは親譲りで、まぬがれる 事は出来ない。
だったら、psでいいじゃん。この形式をブラウザーは喰ってくれるのかと思って、喰わせてみた。だって、表示器のgvなんて、そんじょそこらに転がってはいないからね。ああ、そんじょそこらってのは、Windowsの事です。debianにも*BSDにも普通に有りますけどね。
で、firefoxの言う事には、自前で表示出来ないから、専用のアプリを呼ぶわと言う事でした。 PDFビューワーね。アドビなやつは重いので使っていない。軽いやつが起動したんだけど、 エラーだよ。
で、上で調べておいたpsファイルの書き方を見たら、psデータですよって、シェバング相当を 書くのが作法らしい。（gvは、そんなの無くても、ちゃんと解釈して表示する）
ってな事で、psデータファイルの冒頭にちょっと細工をした。
%! %% ps-file mark (Must be need) 2 2 scale %% x,y scale (option) 0.5 setlinewidth %% line width (option)
なお、左下が原点(0,0)になってて、X,Y座標の第一面の扱い。特別な設定をしない限り、 72DPIの分解能で扱われる。言語体系はForth語だ。通な仕様だな。
和田先生のお勧め、SICPも見ておくか。
やや、maximaの初歩もあるぞ。昔やったはずなんだけど、すっかり忘却の彼方。ここで会ったが 100年目、しっかり再勉強しとけ。
maximaは、MIT製だ。SICPも、勿論MIT製。さすが、横の連携が出来ているな。いや、SICPを 難なくこなす学生をmaximaプロジェクトに取り込んで発展させようという野望が有るに 違いない。
こういう下種の感ぐりはやめれ。折角postscriptに出会ったんだから、postscriptを楽しめ。 ちょいと調べてみると、postscriptは、OpenBSDのportsでは、print系に分類されてた。
ghostscript,ghostview,gvと言った所がその主体になってる。オイラーは何も考えずに、和田 先生の指示（記事）に従って、gvを入れたんだけど、これって、ghostviewの進化形なのね。 じゃ、ghostscriptは何かと言うと、gsってのを提供してるのね。これ、forthに表示系を バンドルしたインタープリタだ。
haskellになぞらえると、ghciに相当するのがgsで、runghcに相当するのがgvって事になるかな。そんな講釈はいいから、gsを使ってみれ。
$ gs GPL Ghostscript 9.07 (2013-02-14) Copyright (C) 2012 Artifex Software, Inc. All rights reserved. This software comes with NO WARRANTY: see the file PUBLIC for details. GS>0 0 moveto GS>200 200 lineto GS>stroke
起動すると、ちょいと長めの真っ白なキャンパスが現れました。いわゆるレターサイズの紙 なんですかね。こんな所にもアメリカ独自の規格が出てくるんか。国際規格？のA4サイズの 紙を指定するには、どうするんだろう？ まあいいや、幼児宜しく、そこに紙があれば、いたずらしたくなります。
movetoで原点に移動。そこから、座標軸200,200に向かって線と言うか、軌跡を書きます。 この時点では、まだ透明になってて線は見えない。次にstrokeって呪文を唱えると、炙り出し 宜しく、その軌跡が見えるようになるとな。そうか、strokeって呪文を知らなければ、秘密の 文を書いておけるのだな。
と、postscript対応のハロワをやった後は、上にあげた実習資料で、forthを十分に楽しみましょう。gvを入れれば、gsはつられて入ってきますからね。簡単に環境を作れますよ。
GS>3 4 GS<2>mul GS<1>== 12 GS>== Error: /stackunderflow in --exch-- Operand stack: --nostringval-- Execution stack: %interp_exit .runexec2 --nostringval-- --nostringval-- --nostringval-- 2 %stopped_push --nostringval-- --nostringval-- %loop_continue --n ostringval-- --nostringval-- false 1 %stopped_push .runexec2 --nostringval-- --nostringval-- --nostringval-- 2 %stopped_push --nostringval -- --nostringval-- --nostringval-- Dictionary stack: --dict:1160/1684(ro)(G)-- --dict:0/20(G)-- --dict:77/200(L)-- Current allocation mode is local Last OS error: 35 Current file position is 3 GS<1>== -file-
自然主義入門なんて言う哲学書を読んでいるんだ。そこに言語はどうやって獲得するなんていう 考察が出てた。チョムスキーとか言う高名な先生によると、文法は生まれつき持っているそうな。周囲（父母等）の観察により、その文法を母語に合わせて自動調整してくとか。そう解釈せざるを得ないそうだ。
よって、オイラーもそれに倣って、色々な事をやるのさ。プロンプト中に表示されてる数字は、 スタックの深さを教えてくれているんだなと、幼児並みの観察をしましたよ。それから、GSを 怒らせたら、腹のうちを示してくれたね。スタックが３本使われているのか。
その哲学書には、生まれつき持ってるって主張する派閥と、それは違うだろう経験で学習する んだよ派が、血を血で洗う対立してるって書いてあった。
経験派の言う事には、最近のAIを見てごらんなさい。深層学習とやらで、色々と知恵を付けて いるじゃんってのがある。大量にあるデータから統計によって学習してくといいんだよ。 なるほどねぇ。こういう分野も面白いな。
haskell + gnuplot
Haskellライブラリ所感2016 なんてのが、今流行のライブラリィーなのね。流行り廃りがあって大変そう。インターフェースが落ち着かないものは使うなが、おじいちゃんの遺訓だな。
そんな事を思いながら、gnuplotをハスけるから使うのは苦労しないと思われ。 The gnuplot package へ行って眺めてみた。
その前に、昔やってたな。stackから難なく入った。枯れてると言う事かな。
deb9:tmp$ stack ghci Configuring GHCi with the following packages: GHCi, version 8.0.2: http://www.haskell.org/ghc/ :? for help Loaded GHCi configuration from /tmp/ghci1365/ghci-script Prelude> import Graphics.Gnuplot.Simple Prelude Graphics.Gnuplot.Simple> let xs = linearScale 1000 (0, 2*pi) Prelude Graphics.Gnuplot.Simple> let ys = fmap sin xs Prelude Graphics.Gnuplot.Simple> let points = zip xs ys Prelude Graphics.Gnuplot.Simple> plotPath [(Title "hello")] points Prelude Graphics.Gnuplot.Simple>
リポジトリーからソースも見られるはずなんだけど、今回はtar玉を落としてきて眺めている。
Simple.hsの中で、上で実行したのがどう定義されてるか確認。
plotPath :: (Tuple.C a) => [Attribute] -> [(a,a)] -> IO () plotPath attrs = plot2d attrs . list
モジュールのリンクを辿って、 Graphics.Gnuplot.Simpleと見比べしてる。全体を把握するなら、htmlの方が便利。ソースは、心の平安の為に。
id
図書館から借りてきた本を全部読み終わって、本の交換に行こうとしたら、女房が化粧を始めた。あーあ、暫く時間がかかるなと思って、手元に有った『プログラミングHaskell』をパラパラ してたんだ。
関数合成演算子の所に、申し訳ない程度に、関数合成は単位元を持つ。それは id ってのが 載ってた。idは引数を変更することなく返す。。。。また、複数の関数を合成する際の初期値と して用いる事が出来るって説明が出てた。さりげなく書いてあるんで、見落としていた。
前回やったcat.hsにもidが使われていて、分かった積りでいたけど、下記のコードで より明確になったよ。
*Main> let flg = [lineCount , tabConv , lineSepConv] *Main> :t flg flg :: [FilterProgram] *Main> let bar = foldr (.) id flg *Main> :t bar bar :: String -> String *Main> putStr $ bar "aa\n\tbb\n" 1 aa$ 2 ^Ibb$
合成する関数のリストを用意。そいつを使ってfoldrで合成。この時に種として使う。 合成子が(+)だと、その単位元は０。関数ならidって事。
flgを畳み込むと、関数達が合成される。そしてそれを実行。 idはfoldrの基底として使われている。
foldr :: (a -> b -> b) -> b -> [a] -> b foldr f z [] = z foldr f z (x:xs) = f x (foldr f z xs)
これがfoldrの定義なので、基底だけが選ばれるように、NULのリストを渡してみる。
*Main> let zz = foldr (.) id []
*Main> :t zz
zz :: b -> b
*Main> putStr $ zz "aa\n\tbb\n"
aa
bb
この時の型は、受け取ったものをそのまま返す、idの定義を示している。これで納得 しました。