diffinvを使って
Rで数値積分をする方法を聞かれたので、ついでにここにも書いておくことにする。数値積分をするためには、基本的にはy*dxを足していけば良い。dxを求めるには、diff(x)とすれば良いが、数が一つ減ってしまうので、そのままではyをかけられない。そこで、最初の要素を除いてy[-1]diff(x)として、それを足し上げていくと数値積分となる。つまり、cumsum(y[-1]diff(x))とすれば良い。 この方法では少し誤差が大きいので、せめて台形積分にしてみよう。yのところは、隣り合ったyの値の平均にするために、(head(y,-1)+tail(y,-1))/2として、cumsum((head(y,-1)+tail(y,-1))/2diff(x))とすれば良い。この最初に0を加えるために、diffinv((head(y,-1)+tail(y,-1))/2diff(x))とすると、xと数が同じになって良いだろう。

Rで多項式やPade近似
Rで係数cnから多項式の値を求めるためには、例えば以下のようにすれば良い。

cn<- 1/gamma(1:11)
x<- -100:100/10
f<-0
for(i in 0:10)f<-f+cn[i+1]*x^i
plot(x,f,type="l")

しかし、値を求めるためにゼロで初期化したりforループを書いたりが面倒である。そこで、xのべき乗の行列を作って、行列の積を使うとスマートに書ける。

f<-outer(x,0:10,"^")%*%cn

これを使って、Pade近似をいじってみた。Pade近似というのは、関数を２つの多項式の比として近似する方法であり、n次の多項式をm次の多項式で割る場合に[n,m]次のPade近似という。例えば、Taylor展開は[n,0]次のPade近似ということもできる。sin(x)の10次のTaylor展開の係数cnは以下で表される。

cn<-c(0,1,0,-1)/gamma(1:11)

これを[4,6]次のPade近似で表してみよう。pracmaというlibraryを使うと係数が簡単に求められる。

library(pracma)
pq<-pade(rev(cn),,4,6)
pn<-rev(pq$r1)
qn<-rev(pq$r2)
f<-(outer(x,0:4,"^")%*%pn)/(outer(x,0:6,"^")%*%qn)

pracmaでは多項式の係数を逆順にしないといけないことに注意が必要である。libraryを使うのは反則かも知れないので、library無しで、係数を求めるプログラムも書いてみた。

mat<-cbind(diag(1,10,4),-sapply(0:5,function(i)c(rep(0,i),cn)[1:10] ))
pqn<-solve(mat)%*%cn[-1]
pn<-c(cn[1],pqn[1:4,1])
qn<-c(1,pqn[5:10,1])

関数によっては、Taylor展開よりもPade近似の方が、良く合う場合があり、 Rで簡単に使えるようになったので、今後は使ってみよう。

ウムラウトも大丈夫
ドイツ語の書類を読む必要があって、どうしようかと思ったが、スキャナーで取り込んで、OCRにかけて、google翻訳することにした。OCRソフトは、tesseractというのを使うと良いようなので、どのドイツ語対応のtesseract-ocr-deuをaptでインストールした。そして、スキャナで取り込んだpdfをpngに変換する。そして、tesseract -l deu+eng image.png tax-0 txtとしたら、それなりに認識してくれた。嬉しいのは、ウムラウトなどもきちんと認識してくれた点だ。ずいぶん前に別の方法でやった時には、ウムラウトやエスツェットが駄目で、手直しの必要があって苦労した覚えがある。画像が乱れているところは、うまく認識できない場合があったが、gimpで調整したら、読み取れるようになった。やはり、アルファベットは種類が少ないので、OCR向けかな。

rubyのclassとmodule
自作の装置で測定する場合には、rubyを用いて測定している。PCと機器で通信する場合には、インターフェースとしてGPIBやRS232CやUSBTMCを主に使っているが、これらを切り替える必要がある場合がある。当初は、機器に対してインターフェースを固定して、それ専用のプログラムを書いていたが、インターフェースを変えるために、別のプログラムを作らないといけないので面倒である。 現在ホームページに公開しているversionでは、インターフェースをclassとして定義して、機器はそれらのサブclassとして定義することによって、インターフェースを指定している。インターフェースを切り替えるためには、機器を記述しているプログラムのスーパーclassの部分を書き換えるだけで良い。スーパーclassは後から変えることは出来ないようなので、プログラムの書き換えが必須であった。 新たに装置を立ち上げるにあたって、いろいろと考えていたら、moduleを使えば、インターフェースと装置のプログラムを書き換えずに、測定プログラムからインターフェースを切り替えるという方法を思いついた。つまり、インターフェースはmoduleにして、測定プログラムの中から、機器のclassを再定義してincludeさせるのである。測定プログラムの先頭が少し長くなるが、そこでそれぞれの機器のインターフェースが明示されるのは、良い点とも考えられる。もう少しプログラムが固まってきたら、公開するかも知れません。

RとrglとwebGL

rglを使うと、Rで三次元の図を作って、マウスで回したりすることができる。webGLを利用して、htmlにすることもできる。

library('rgl')
open3d()
shade3d(cube3d(),col='red')
writeWebGL(filename='index.html')

とすると、index.htmlが出来て、これをブラウザーで表示させると、三次元の図が現れる。そのままの設定だとwebGLというフォルダが出来てしまうけど、まあいいか。このhtmlには、長いjavascriptが埋め込まれている。これを取り除いて、図形のデータだけにしたい場合は、

writeWebGL(filename='index.html',commonParts=FALSE)

とすればよい。しかし、そのままでは図形を表示できない。/usr/lib/R/site-library/rgl/htmlwidgets/lib/CanvasMatrix/CanvasMatrix.src.jsと/usr/lib/R/site-library/rgl/htmlwidgets/lib/rglClass/rglClass.src.jsを同じフォルダに置いて、"src."を除いたファイル名に変更して、htmlの最初の方に

<script src="CanvasMatrix.js" type="text/javascript"></script>
<script src="rglClass.js" type="text/javascript"></script>

を加える必要がある。また、同じhtmlに複数の図を置く場合には、div idを変更する必要がある。

<div id="div" class="rglWebGL"></div>
<script type="text/javascript">
	var div = document.getElementById("div"),

にある二つの"div"を重複しない名前に変更すれば良い。この方法がなかなか分からずに苦労した。

pdfとjpegの変換
ImageMagickを使って、pdfをjpegに変換しようとしたら、security policyのせいで出来ないと言われてしまった。昔はできたはずなんだけど。調べてみたら、回避方法が分かった。/etc/ImageMagick-6/policy.xmlでPDFの制限をしている行をコメントアウトすれば良いようである。実際にやってみたら、うまく動くようになった。

なかなか気づかなかった
lubuntu20.04でpdfを見たときには、使い勝手がこれまでと違うと思っていたら、実はevinceでは無いことにようやく気がついた。新たにevinceをインストールしたので、これからは同じ使い勝手で使うことができるでしょう。

思ったより簡単
最近lubuntu 20.04をインストールしたが、そのPCは時々しか使わなかったので、日本語入力が出来ないことになかなか気が付かなかった。fcitx-mozcを入れて、再ログインしたら、入力できるようになった。昔は言語設定とかしていた気がするけど、それは必要無くなったようだ。

一応復活
ディスプレイの交換作業をしたら、測定用のPCが立ち上がらなくなった。半年位前にXの部分でかなり苦労しながら立ち上げたものなのだが、おそらくその部分で立ち上がらなくなった。AMDのA4-7300というAPUを搭載したマシンなのだが、その時期のAPUはLinuxとの相性が良くないらしい。立ち上げる際は、kernelを4系統から5系統に置き換えたら、なんとかうまく行った。しかし、停電などで電源を落とした後には、いつも立ち上がるかヒヤヒヤものであった。予想通りと言うべきか、今回は立ち上がらなくなった。 grubのオプションでnomodesetを指定すると、X以外は立ち上がるので、その状態で折角なのでapt-get updateとupgradeをしたりして、kernelもより新しいものに入れ替えた。いろいろと試してみたが、なかなか立ち上がらなかった。 検索で同じマザーボードを使っている人のブログを見つけたが、そこではBIOSからover clockingの指定をdisableにしたら、うまく行ったという記述を見つけた。しかし、これも駄目だった。今から思い返すと、使っているPCのBIOSのVersionが古く、ブログのは最新のVersionを使って居たので、BIOSの更新をしたら、状況が変わっていたかも知れない。 さらにいろいろと試していたら、GRUBでquiet splash radeon.dpm=0と指定をして、DPM (Dynamic power management)を切ると、Xが立ち上がることを発見した。また立ち上がらなくなるのが怖いので、電源は落としていないが、これで安定してくれることを期待したい。 次に立ち上がらなくなったら、BIOSの更新をしてみよう。

Rのくせ

Rをいじっていたら、そのくせに苦労したので、その対処法をメモしておく。

まず苦労したのが、forループを使ってn回ループしようとして、n==0のときにうまく行かない点である。

for(i in 1:n){cat(i)}

とすると、n==0で実行すると1と0で二回実行されてしまうのである。いろいろと考えた結果、

for(i in (1:n)[0:n]){cat(i)}

とすると、うまく行くことが分かった。または、

for(i in seq_len(n)){cat(i)}

の方が素直かな。

次に苦労したのが、1行または1列の行列が自動的にベクトルに変換される点である。diag(3)[-1,-1]とすると2行2列の行列になるが、diag(2)[-1,-1]とすると、1行1列の行列になると思ったら、ベクトルになってしまうのである。これを行列にしたかったので、matrix(diag(2)[-1,-1],1)とした。美しくは無いが、仕方ないかな。

2021/4/6追記

for(i in (0:n)[-1]){cat(i)}

でも良いかな。