ME差 (2)

やませが吹いて、今年の梅雨は寒い。農家の方は、さぞヤキモキされている事でしょう。でも季節は順当に巡ってきて夏野菜が豊富に出回るようになりました。

近くのマルシェ（気取ってみました）へ行ったら、今まで盛んに並んでいたブロッコリーが品薄。その代わりにキュウリが７本で100均価格、丸茄子２個100均。長茄子4個のセットとか、モロヘイヤとかズッキーニとかね。あれ、ぬるぬる系のオクラはまだだなあ。

人気は大玉の梅と赤紫蘇かな。おばちゃん連中が群がっていた。紫蘇ジュースはおばちゃんパワーに圧倒されて、原材料を入手出来ず。じゃ、トマトはどうよ？リコピンと何とかいうのが美と若返りの元ってTVで盛んに宣伝してたけど。何故か高価なんだよな。ってか、いつも高止まりでした。

季節外れっぽいイチゴが１パック150円。食べ収めと思って２パック入手。果物大好き女房は、酸っぱいからってんで、自家製のジャムに変身させてた。

今日は夏野菜のトマト煮。トマトはイタリアから輸入のホールトマトだった。EU産の方が国産より安いな。明日は、夏野菜のカレー煮かな。モロヘイヤはミキサーしてスープにするのが一番美味しいと思うけど、どうよ？

もう少しすると、えだ豆のシーズンだなあ。最近は農業のIoT化で、枝豆の発育状況と画像判定してるとか。美味しい枝豆の影にAI有りだな。

debian dist-upgrade

某みずほ銀行に倣って、３連休中の早朝にシステム更新を実施する。休日のしかも早朝を狙ったのは、世間様にご迷惑をかけない為。

debian機(i386版）を、9.9から10.0に上げる作業。今までは新しいのが出ると、古いのを消して新しいのに移るって事をやってたけど、実機なのでそれは避けたい。よって、こういう作業は初体験です。

# sed -i 's/stretch/buster/g' /etc/apt/sources.list

まずは、こんなコマンドを叩いて、パッケージの取り寄せ先を変更します。同じ敷地で倉庫名を変えるという、ことなかれ主義。パッケージダウンロードに利用するミラーサイトについてなんてのを見て、あちこち彷徨うのもいいかも知れません。

後は、 apt(update,upgrade,dist-upgrade) 左式を展開実行するだけ、と軽く考えていたんだけど、人生は躓きの連鎖って法則が適用されたぞ。

upgradeの所で、お前の所に入ってるやつは10.0でパッケージになってないものが有るから、諦めろと言う。それにWiFi用に入ってるping系がコンフリクトだよとも言われたぞ。そんな話何処にも出てなかったぞ。LXDEバージョンで、WiFi使用ってどこにもありそうな構成だと思うんだけどな。

先行者のUさんからのレポートが参考になった。(Tnx Uさん）

# apt-get update
# apt-get upgrade --fix-missing
( # apt-get -s dist-upgrade は好みで)
# apt-get dist-upgrade
#
# reboot

 再起動前の "apt-get update～apt-get upgrade"時に最初エラーが出て、そのメッセー
  ジの中で --fix-missing オプションが推奨されたのでそれに従いました。update 時点
  でエラーが出なかったら --fix-missing は要らないかもしれません。

  また、"apt-get dist-upgrade"の最中にいくつか対話式に確認が求められましたが、わ
  たしの場合は全部デフォルト(の状態のまま)で返しました。

オイラーの所では、つれなくされたんだけどな。人をみるのか？

ぶつくさ言いながらも、５時に工事を始めて７時半に終了。途中、sshでモニターしてた回線は切れる（別IPになった為）は、GUIの画面が１分程消えるはでひやひやしたよ。でも無事に再起動出来た。

upgradeに伴って、 virtualbox-5.2が破損扱いになってた。i386版のやつはもうディスコンなんだな。６系は勿論amd64版のみの提供でしたよ。しょうがないな。

ついでに、Windows10のWSLに入ってるdebianもupgradeした。こちらは何のトラブルも無く完了したよ。

また、FreeBSDは３ケ月毎に、portsが一斉更新されるんで、pkg upgradeしといた。

fsolve

連立1次方程式とか、N次方程式の解法は、前回調べたけど、これらが合体した方程式は、どうやって解けばいいの？そんなの、fsolveを使えらしい。ただ、一発では解けず、ちょっと汗を流さないといけないらしい。以下は、fsolveのマニュアルから抜粋。非線型連立方程式(日本語解説)

   The following is a complete example.  To solve the set of equations

     -2x^2 + 3xy   + 4 sin(y) = 6
      3x^2 - 2xy^2 + 3 cos(x) = -4

you first need to write a function to compute the value of the given
function.  For example:

     function y = f (x)
       y = zeros (2, 1);
       y(1) = -2*x(1)^2 + 3*x(1)*x(2)   + 4*sin(x(2)) - 6;
       y(2) =  3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4;
     endfunction

   Then, call ‘fsolve’ with a specified initial condition to find the
roots of the system of equations.  For example, given the function ‘f’
defined above,

     [x, fval, info] = fsolve (@f, [1; 2])

results in the solution

     x =

       0.57983
       2.54621

     fval =

       -5.7184e-10
        5.5460e-10

     info = 1

上記を見つける為に、彷徨っていたら、シミュレーションこんな、関連のありそうな物も発見したよ。また、非線形方程式の解法なんてのもあった。

そして、超力作 Octaveによる数値計算入門

近似

よく使いそうなので、コピペしたよ。関連が次々と出てくるな。

・ 最小二乗近似
>>  x=[  1    2    3    4    5]
>>  y=[  4    5    7    9    12]
>> polyfit(x,y,1)
ans =
  2.00000 --- 傾き
  1.40000 --- 切片   y = 2x + 1.4

・ 二次関数近似
>> polyfit(x,y,2)
ans =
  0.28571
  0.28571
  3.40000            y = 0.28571x^2 + 0.28571x + 3.4

shadows

[ham@cent tmp]$ octave -q
warning: function ./speed.m shadows a core library function

これ、たまたま前回のME差プログラムしてた時に遭遇した警告。プログラムの実行が速くなるには、速い名前を付けるに限るって事で、speed.mって名前にしたんだ。

この名前がコアライブラリィー名と被ってしまって、大事な方を覆い隠しちゃったよっていう重大警告だ。

octave:1> help speed
'speed' is a script from the file /tmp/speed.m

 me blood

自作の奴を脇に避けてから起動してhelpしてみると

'speed' is a function from the file /app/share/octave/5.1.0/m/testfun/speed.m

 -- speed (F, INIT, MAX_N, F2, TOL)
 -- [ORDER, N, T_F, T_F2] = speed (...)

     Determine the execution time of an expression (F) for various input
     values (N).

オリジナルのやつが見えてきた。current-dirが最初の検索pathになって、ダブりをチェックしてくれてるのね。

octave:9> sin(pi)
ans =    1.2246e-16
octave:10> sin = 1212
sin =  1212
octave:11> sin(pi)
error: sin(3.14159): subscripts must be either integers 1 to (2^63)-1 or logicals (note: variable 'sin' shadows function)

sin関数。同名な名前に数値をバインディング。何の警告もなく実行された。もうだめだ。シャドーとか注意深くみるとミスった事は分かるけど。

octave:11> help cos
'cos' is a built-in function from the file libinterp/corefcn/mappers.cc

根幹の部類なんだけどな。

octave:12> cos(0.2)
ans =  0.98007
octave:13> a = @cos
a = @cos
octave:14> a(0.2)
ans =  0.98007

別名で保存しておけとな。これらから、octaveはLisp-1の挙動をしてるんですな。

dict

前回書いた、me.mなんだけど、期待したスピードは出ていなかった。forの遅いのは諦めるとして、気になるのはlookupを使ってる事。検索するテーブルはソートされてる事って条件から、容易にバイナリーサーチしてる事が想像出来る。

普通そんなのは、連想配列を使うしょ。octaveにはそういう構造は無いのか？前から気になってた、 matlab: セル配列で困っている人へを調べるも、あくまで配列の扱いだ。

ならばってんで、octave key valueあたりで検索すると、 containers.Mapがヒットする。ぴったりじゃんと思うのもつかの間、まだ実装されていない。（未実装品をマニュアルに載せないでください。＞octaveの偉い人）

２番目のヒットは、 dict (keys, values)と言う奴。generalって言うパッケージに収録されてるとな。将軍じゃなくても一般人が使いたそうなやつである。でも、これを入れなくても、何とかなりそうな雰囲気を感じる。ちょっと実験。

octave:47> qq.(123)=456
error: wrong type argument 'scalar'
error: dynamic structure field names must be strings
octave:47> qq.('123')=456
qq =

  scalar structure containing the fields:

    123 =  456
octave:49> qq.('123')
ans =  456
octave:50> qq.('4395')
error: structure has no member '4395'

キーには数値を使えません。文字列でお願いします。登録されていないキーに対してアクセスすると、エラーになると言う厳しい制約があるようだ。それから、内部表現が構造体って事で、多分スピードは期待出来ないと予測。でも、一応使えるか検証してみたいぞ。

tic();
bld = csvread('current.csv');
am  = bld( rem(bld(:,1),100)<12, [1 2] ); % [ [ymdh hi lo ps] ... ]
pm  = bld( rem(bld(:,1),100)>12, [1 2] );
am  = [fix(am(:,1)/100) am(:,2)];         % [ [ymd hi] ... ]
pm  = [fix(pm(:,1)/100) pm(:,2)];
global qq
for i=1:length(pm)
  k = num2str(pm(i,1)); v = pm(i,2);
  qq.(k) = v;
endfor

function v = pick(x)
  global qq
  try
    k = num2str(x);
    v = qq.(k);
  catch
    v = 0;
  end
endfunction

toc()
res = zeros(2,2);
for i=1:length(am)
   x = pick(am(i,1)) ;
   if (x)
     av = (am(i,2) + x) / 2.0;
     dl = am(i,2) - x;
     if (av >  135 && dl > 15)   res(1,2) +=1; endif
     if (av <= 135 && dl > 15)   res(1,1) +=1; endif
     if (av >  135 && dl<= 15)   res(2,2) +=1; endif
     if (av <= 135 && dl<= 15)   res(2,1) +=1; endif
   endif
endfor
res
sum(sum(res))
toc()

結果は勿論、前回書いたものと一致した。

で、気になる実行スピードだか、激遅。参考までに、下段に初回に書いたもののスピードを掲載。

Elapsed time is 2.41655 seconds.
Elapsed time is 5.08257 seconds.

Elapsed time is 0.010308 seconds.
Elapsed time is 0.129669 seconds.

何の得にもならなかったじゃんと言う否定的な見解は伏せておこう。try catch の実例を試せたと前向きに受け取っておく。

speed

間違って上書きしたコードのオリジナルは、スピード評価用の試験機だった。

下記は実行例。

octave:6> 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  ...  n14 = 84834  n15 = 100000

Mean runtime ratio = 1.84e+03 for 'v = 0; for i = 1:length (x), v += x(i); endfor' vs 'sum (x)'

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

randn()で任意長の乱数を発生させ、その合計を組み込み関数sumでやるか、forで回すかの対決させたのも。結果は３桁違う。グラフも出てくるので参考に。

demoでは、こんなのが出てきた。

-----------------------
function x = build_orig (n)
  for i=0:n-1, x([1:100]+i*100) = 1:100; endfor
endfunction
-----------------------
function x = build (n)
  idx = [1:100]';
  x = idx(:,ones(1,n));
  x = reshape (x, 1, n*100);
endfunction
-----------------------
  :
Mean runtime ratio = 27.8 for 'build_orig (n)' vs 'build (n)'

初心者が書くやつに比べて、プロが書くやつは、27倍も速かったよと。forをいかに撲滅して組み込みだけでやるか、クイズを出されているみたいですなあ。

example('speed')とかすると、デモのソースを見れるぞ。使い方が分かって面白い。

containers.Map

上で無い々と騒いだcontainers.MAPだが、octave 5.1.0にはちゃんと実装されてた。ひょっとしてOpenBSDで無いと誤解してたのかも知れない。と言う事で、上の差分。

--- sme.m       2019-07-15 14:20:56.946463147 +0900
+++ cme.m       2019-07-15 16:48:24.944187356 +0900
@@ -5,16 +5,12 @@
 am  = [fix(am(:,1)/100) am(:,2)];         % [ [ymd hi] ... ]
 pm  = [fix(pm(:,1)/100) pm(:,2)];
 global qq
-for i=1:length(pm)
-  k = num2str(pm(i,1)); v = pm(i,2);
-  qq.(k) = v;
-endfor
+qq = containers.Map(pm(:,1),pm(:,2));

 function v = pick(x)
   global qq
   try
-    k = num2str(x);
-    v = qq.(k);
+    v = qq(x);
   catch
     v = 0;
   end

登録は一発で出来たけど、アクセスする時は、やはりtry catchで防御が必要。キーが無かったら、適当な値を返すようにしておいてくれれば、、、気になるスピード

Elapsed time is 0.0478568 seconds.
Elapsed time is 4.33144 seconds.

やっぱり、遅いね。