octave

1Kgの重さの原器から新しい方式へと移行するそうな。

(1)  E = hf
(2)  E = mc^2
(3)  m = hf/c^2

1式は、電磁波のエネルギーは、周波数(f)とプランク定数(h)の積だよと。2式はアインシュタインの有名な式。エネルギーは、質量(m)と光速(c)の2乗との積だよと。この２つの 式から、3式が導き出せる。

周波数は１秒間に振動する回数。１秒は、日本の誇る光子時計で十分に精度よく決定出来ている。光の速さも十分な精度が得られている。後は物理定数のプランク定数を精度よく決定出来れば、質量もおのずと未来永劫に精度を担保出来る。

質量の単位「キログラム」の新たな基準となるプランク定数の決定に貢献

ワット天秤

国際単位系（SI）kg再定義の舞台裏

アボガドロ定数を正確に求める方法もあるんだな。

setup octave

octaveは、debian機に入っている。Windowsなマシンの仮想PCにも入れておこ うと思ったんだ。CentOSに入れたらセグフォ。これはまずいな次はウブかな。

plotが動かない、ってか画面をコピペする様な変な動作。ならばFreeBSDだな。 随分と余計な物もインストールさせられた。(Java,lua,python,gcc,clang,gnuplot)

$ octave
GNU Octave, version 4.4.1
 :
octave:1> hold on
octave:2> libGL error: No matching fbConfigs or visuals found
libGL error: failed to load driver: swrast

何これ？ ググるも解決法には届かず。でも、ヒントがあった。mesaのライブ ラリィーに足りない物があるから、WSLなウブな場合はどうしろ、こうしろと。

CUIな端末からXを飛ばしているから、Windows側のVcXsrvが満足しないんだろう。 ならば、twmを起動して、そこで試してみよう。ビンゴでしたね。ちゃんとグ ラフが描かれた。これにて、一件落着でもいいんだけど、いまいち動きが遅い。

それで思いついたんだ。昔盛んに使っていたMobaXtermは、どうかと。ビンゴ でしたね。これで、CUIな端末から操作してて、グラフ画面だけを別に表示で きる。操作履歴も容易にコピペ出来るってもんです。

$ echo $DISPLAY
localhost:10.0

MobaXtermの場合、起動すると勝手に環境変数がセットされるんで何も設定す る必要無し。

余談だけど、NetBSDの注意書きに面白い事が書いてあった。クラッシュを体験 なされたなら、、、、ってね。

$NetBSD: MESSAGE,v 1.2 2017/01/24 12:31:18 maya Exp $

If you are experiencing crashes while plotting, you can try adding the
following line to ~/.octaverc:
        graphics_toolkit("gnuplot");

まくまくOctaveノート

octave:6> x = [1; 0]
x =

   1
   0
octave:7> c = [3 2; 1 4]
c =

   3   2
   1   4
octave:8> c * x
ans =

   3
   1

octave:9> [x,y] = meshgrid(-2:2)
x =

  -2  -1   0   1   2
  -2  -1   0   1   2
  -2  -1   0   1   2
  -2  -1   0   1   2
  -2  -1   0   1   2

y =

  -2  -2  -2  -2  -2
  -1  -1  -1  -1  -1
   0   0   0   0   0
   1   1   1   1   1
   2   2   2   2   2

ここからは、少しコードを書いてみる。線形変換行列 c を受け取って変換す る奴。

$ cat xx.m
function xx(c)
  r = -3:3
  for x = r
    for y = r
      o = [x;y]
      z = c * o
      x2 = z(1)
      y2 = z(2)
      quiver(0,0,x2,y2)
    end
  end
end

演算過程が見える様に、式の後ろにセミコロンを付ていない。

octave:16> c = [2 1.5; -0.8 1.2]
c =

   2.00000   1.50000
  -0.80000   1.20000

octave:17> hold on
octave:18> xx(c)
r =

  -3  -2  -1   0   1   2   3

  :

o =

   3
   3

z =

   10.5000
    1.2000

x2 =  10.500
y2 =  1.2000

octave-mode

気を良くしたんで、emacsから編集しようとしたんだ。そしたら、あろう事か ObjCを編集するモードになったぞ。アプル御用達なサフィックスですかい？ m ってのは。強制的にoctave-modeになる様に下記の設定をした。

(setq auto-mode-alist
       (append '(("\\.m$" . octave-mode))
           auto-mode-alist))

M-x run-octaveで画面が割れてreplが動き出すのはお約束。編集して実行が楽 で良いわい。

function xx(c)
  hold on
  r = -1:1;
  for x = r
    for y = r
      o = [x;y];
      z = c * o;
      x2 = z(1);
      y2 = z(2);
      quiver(x,y,x2,y2, 0,  "b")
    end
  end
  disp("Hit any key to hide graph ...")
  pause()
  close
end

始点を線形変換した先を矢印で示すようにチョイと変更。更にグラフが自動で 作られる様にした。

octave> xx([2 -0.8; 1.2 0.5])
Hit any key to hide graph ...
q
octave>

こんな風にして、ワープを実感してくれ。それにしても計算に時間がかかり過 ぎのように思える。なんでかな？ 実行時間を測るには、どうする？ システムに聞くのが一番だな。

octave> lookfor time
 :
tic                 Initialize a wall-clock timer.
toc                 Measure elapsed time on a wall-clock timer.
speed               Determine the execution time of an expression (F) for vario
                    us input values (N).

使えそうな物を見繕ってみた。そして、speedの例。

octave> speed ("sum (x)", "", [10000, 100000], ...
                 "v = 0; for i = 1:length (x), v += x(i); endfor")
> testing sum (x)
init: x = randn (n, 1)
n1 = 10000  n2 = 11788  n3 = 13895  n4 = 16379  n5 = 19307  n6 = 22758  n7 = 26\
827  n8 = 31623  n9 = 37276  n10 = 43940  n11 = 51795  n12 = 61054  n13 = 71969\
  n14 = 84834  n15 = 100000

Mean runtime ratio = 235 for 'v = 0; for i = 1:length (x), v += x(i); endfor' v\
s 'sum (x)'

For sum (x):
  asymptotic power: O(n^1.2)
  approximate time per operation: 9 ns

ふむ、乱数な配列のサイズを10000から100000の間で発生させて、その合計を 求める計算をさせているんだな。そして、配列のサイズによって、計算時間 がどう変化してるか求めている。リニアになるかと思たら、1.2乗になりま したとな。良いアルゴリズムを求めてくださいって趣旨がよく分る。

今回はもっと大雑把な奴が欲しい。tic,tocを試してみる。

octave> tic
octave> toc
Elapsed time is 3.52813 seconds.
octave> toc
Elapsed time is 14.688 seconds.

ticがストップウォッチを初期化して、tocで経過時間ね。

octave> xx([2 -0.8; 1.2 0.5])
Elapsed time is 1.62495 seconds.
Hit any key to hide graph ...

待ってる時間って長く感じるものなのね。これって作ったデータをその都度 gnuplotに渡しているのかな。

octave> demo quiver
  ：
quiver example 3:
 clf;
 [x,y] = meshgrid (1:2:20);
 h = quiver (x,y, sin (2*pi*x/10), sin (2*pi*y/10));
 set (h, "marker", "o");
 title ("quiver() plot w/origin markers and arrowheads");

Press <enter> to continue:

この例をファイルにコピペして、実行時間を計測。

octave> aa
Elapsed time is 0.127378 seconds.

オイラーが書いたコードより、よっぽど複雑なプロットなのに、断トツで実行 時間は短かい。多分gnuplotには一回しかデータを渡していないのだろうね。

/usr/local/share/octave/4.4.1/m/plot/draw/quiver.m

function h = quiver (varargin)
  [hax, varargin, nargin] = __plt_get_axis_arg__ ("quiver", varargin{:});
    :
  unwind_protect
    hax = newplot (hax);
    htmp = __quiver__ (hax, false, varargin{:});

    ## FIXME: This should be moved into __quiver__ when problem with
    ##        re-initialization of title object is fixed.
    if (! ishold ())
      set (hax, "box", "on");
    endif

主要部分はこれっぽい。__quiver__ が、プリミティブで、演算とかgnuplotと のやり取りをしてるんだろうね。wnwind_protect 失敗するか も知れない演算を表わすlisp風の防御関数だから。

octave-forge

好き者が作った、Community packagesが集積した場所がある。

octave-forge

色々あるなあ。pkg -forge とかやると、インストールできるようだ。 でも、FreeBSDには、色々と用意されてるぞ。

$ sudo pkg install octave-forge-ltfat
  ：
The following 3 package(s) will be affected (of 0 checked):

New packages to be INSTALLED:
        octave-forge-ltfat: 2.3.1_2,1
        cblas: 1.0_9
        portaudio: 19.6.0,1

Number of packages to be installed: 3

The process will require 8 MiB more space.
7 MiB to be downloaded.

後は、これをどう使うかだなあ。

gnuplotの、あれ

あれとかそれとか言うようになると、ボケの初まりらしい。 ちゃんと、vectors って言えよ。gnuplotのヘルプによると

gnuplot> help vec

 The 2D `vectors` style draws a vector from (x,y) to (x+xdelta,y+ydelta).
 The 3D `vectors` style is similar, but requires six columns of basic data.
 In both cases, an additional input column (5th in 2D, 7th in 3D) may be used
 to provide variable (per-datapoint) color information.
 (see `linecolor` and `rgbcolor variable`).
 A small arrowhead is drawn at the end of each vector.

      4 columns:  x  y  xdelta  ydelta
      6 columns:  x  y  z  xdelta  ydelta  zdelta

 The keywords "with vectors" may be followed by an inline arrow style
 specifications, a reference to a predefined arrow style, or a request to read
 the index of the desired arrow style for each vector from a separate column.
 Note: If you choose "arrowstyle variable" it will fill in all arrow properties
 at the time the corresponding vector is drawn; you cannot mix this keyword with
 other line or arrow style qualifiers in the plot command.

      plot ... with vectors filled heads
      plot ... with vectors arrowstyle 3
      plot ... using 1:2:3:4:5 with vectors arrowstyle variable

 Example:
       plot 'file.dat' using 1:2:3:4 with vectors head filled lt 2
       splot 'file.dat' using 1:2:3:(1):(1):(1) with vectors filled head lw 2

/usr/local/share/examples/gnuplot/vector.dem あたりを参考にすれば良い のかな。何やらダイポールアンテナの電界分布を表示する例だ。こういうのは、 CQ誌のあの人の得意分野。マクスウェルの電磁方程式を計算してるんだな。

#set table $field2xy
set table "tab.csv" separator comma
splot vtot(x,y) w l
unset table
pause 0

tableなんていうコマンドが用意されてて、(s)plotコマンドで描画されるデー タをin-memoryで保持できるみたい。それを横取りして、内容を覗いてみた。3 次元の世界で電波は拡がっていくんだなあ。(あったり前よ。まて、時間の次 元が抜けていないか。幾ら電波伝搬が速いと言っても、有限だぞ。過渡応答は 難しいので、定常状態の解析なんだな)

$ head tab.csv

# Surface 0 of 1 surfaces

# Curve title: "vtot(x,y)"

# IsoCurve 0, 100 points
# x y z type
-10 -10 -0.00706214 i
-9.79798 -10 -0.00713123 i
-9.59596 -10 -0.00719778 i

gnuplot-mode

gnuplotもある意味、スクリプトを書いて動かす、インタープリタだ。ならば、 emacsから使えるようになっているはず。

Emacs の gnuplot-mode で gnuplot を快適に

やっぱりあったな。スクリプトを書いてから、C-c C-bでバッファーの内容を gnuplotへ送りつける事が出来る。(あるいは、C-c C-fでファイルを喰わせる 事もできる)

test.gp

set xr [0:8]
set yr [0:8]
set grid

plot '-' using 1:2:3:4 with vectors head filled lt 2, \
     '-' using 1:2 w points pt 'o'

2 3 1 -2
3 2 2 2
4 4 0 1
end
2 3 1 -2
3 2 2 2
4 4 0 1
end

pause -1 "Hit return to continue"

ベクタープロット・データ付きなスクリプト。shellからgnuplot test.gp し ても良いし、replな場面から、load "test.gp" しても良い。

ベクタープロットで始点にmakrする方法が見付からなかったので、グラフの重 ね書きをしてる。その為、必要なデータを二回繰りかえしてしまった。 本来はデータは、ファイルに用意するんで、問題は無かろう。

maxima

octaveを入れた時に、随分と色々な物が付属してきた。ここまで数学寄りな奴 が入っているなら、maximaを入れてもCommon Lispぐらいしか必要ないだろう。

試しに入れたら、sbclとtkだけが足りないから入れるよと言われて、完了した。 これで、現在のFreeBSDは、数学学習の為のOSになったな。まるでfedoraのあ のセットみたい。

wxmaximaなんてのも有ったなと思って入れてみたけど、エラーになってしまっ て使えなかった。よって、emacsからだな。

$ locate maxima.el
/usr/local/share/emacs/26.1/lisp/org/ob-maxima.el.gz
/usr/local/share/emacs/26.1/lisp/org/ob-maxima.elc
/usr/local/share/maxima/5.41.0/emacs/emaxima.el
/usr/local/share/maxima/5.41.0/emacs/imaxima.el
/usr/local/share/maxima/5.41.0/emacs/maxima.el

こんなのが付属してたんで、setupした。

;; maxima
(autoload 'maxima-mode "maxima" "Maxima mode" t)
(autoload 'maxima "maxima" "Maxima interaction" t)
(setq auto-mode-alist (cons '("\\.mac" . maxima-mode) auto-mode-alist))

maxima業界のサフックスは、.macなのね。知らんかったわい。

lisp繋りで、面白い記事を発見。

Lispの車窓から見た人工知能

機械学習は、Pythonが覇権を握ってしまったけど、昔は頑張っていたんだよ。 たまには、素のLispと戯れてみるか。