Aitken加速[programming][lisp][python][rust]
概要
数列のAitken加速というのを説明している http://cympfh.cc/aiura/_/aitken.html を見て、真似た。昔授業で習ったはずだが、すでに記憶の彼方だったので。あと最近Rustについて調べてみてるのでその練習も兼ねて幾つかの言語で書き比べてみようとした(python, rust, common lisp)。
数列の加速
http://cympfh.cc/aiura/_/aitken.html の通り、
sに収束する数列が
(|λ|<1)というふうに収束していくとき、隣り合う3項を使ってβとλを消去してやるとsに非常に近い値が得られ、結果として元の数列より早く極限値に収束する数列となるというもの。
解いた結果としては
となる。こうやって数列から数列
を得る手続きをAitken加速という。さらにこの手続きを
にもう一度行うことで、2段のAitken加速となる。k回繰り返したものをk段のAitken加速列
という。特に
。
が有限個しかないと、Aitken加速で数列の長さが2ずつ小さくなるので、N項の数列ではN/2段程度までしか行えない。
参考リンク2つ目の内容を読んだメモ:
数列がsにp次収束するとは、
のときいう。またsに収束する数列が極限値
が存在するときは|λ|<1になるのだが、λ=1の場合を対数収束、λ=0を超線形収束、 を線形収束という。
超線形収束の場合は加速法は普通はなくても実用上問題ない。対数収束する数列の加速は一般に難しい。
Aitken加速は任意の線形収束する数列を加速する。数値積分とかに使われる。
実装
python(2.7), rust, common lisp (roswell script)で多段Aitken加速を実装してみた。
参考リンク1つ目と同じく、 という数列の最初の20個をとってきて、順次加速している。
どれも実行すると
xs^0: 1.000000, 1.250000, 1.361111, 1.423611, 1.463611, 1.491389, 1.511797, 1.527422, 1.539768, 1.549768, 1.558032, 1.564977, 1.570894, 1.575996, 1.580440, 1.584347, 1.587807, 1.590893, 1.593663, 1.596163 xs^1: 1.450000, 1.503968, 1.534722, 1.554520, 1.568312, 1.578464, 1.586246, 1.592399, 1.597387, 1.601510, 1.604977, 1.607931, 1.610479, 1.612698, 1.614650, 1.616378, 1.617920, 1.619304 xs^2: 1.575465, 1.590296, 1.599981, 1.606776, 1.611798, 1.615658, 1.618716, 1.621197, 1.623250, 1.624977, 1.626449, 1.627720, 1.628827, 1.629801, 1.630664, 1.631434 xs^3: 1.618209, 1.622752, 1.626025, 1.628473, 1.630368, 1.631876, 1.633103, 1.634120, 1.634977, 1.635709, 1.636341, 1.636892, 1.637377, 1.637807 xs^4: 1.634458, 1.635750, 1.636852, 1.637743, 1.638467, 1.639062, 1.639559, 1.639980, 1.640341, 1.640653, 1.640925, 1.641165 xs^5: 1.643242, 1.641498, 1.641588, 1.641833, 1.642077, 1.642299, 1.642495, 1.642660, 1.642825, 1.642911 xs^6: 1.641584, 1.641444, 1.444996, 1.644407, 1.644065, 1.643488, 1.618429, 1.643004 xs^7: 1.641584, 1.543955, 1.644065, 1.644906, 1.644078, 1.630836 xs^8: 1.593382, 1.644913, 1.644489, 1.644961 xs^9: 1.644492, 1.644712
という出力になる。この数列はpi^2/6 = 1.644934...に収束する。
実行速度は Rust > Common lisp > Python という感じ。しかしRustはコンパイラの警告が厳しくて*1、あんまり気楽に書けないなあ。pythonは全体的に簡潔に書けてやっぱりうまくできた言語だと感じる。でも少しでも変則的な accumulation になるとlisp(というかloopマクロか)の方がずっと楽。ここのところPythonをずっと書いていたけどそこはずっと不便に感じていた。lispはlispで配列参照と数式を書くのが僕にはちょっと気が重いんだけど*2。
Rust
あとPython版については ./python aitken.py plot として起動すると加速した数列を重ね書きplotするようにしてみた。
加速1段ごとに最初から収束値(点線)に近いところに行っているのが見える。6段くらいで不安定になっている。要するにが小さくなっているのだが、多段で加速しているうちに誤差が積もるかというとまだdoubleだったら余裕な気がするので、あんまり線形収束的な振る舞いをしていない部分ではそうなるということに思われるが。
参考
*1:その割にはindexが範囲外とかのエラーでは落ちるし、いやセキュリティを気にする用途なら落ちて正解なんだろうが
*2:http://qiita.com/y2q_actionman/items/0a7737710ba647697832 あたりを参考になんか作るといいのかもしれない
iphoneからSSHのsocksプロキシを使う(macbookがあるとき)
iphoneとmacbookが同じwifiネットワーク上にあるとき、macbookからのSSH接続をiphoneからsocksプロキシとして使うメモ。
需要的には海外出張でキャリア回線が使えない上にどうにか日本からしか見られないコンテンツを見たい、自宅サーバーにSSHでトンネルすれば…みたいなものすごいニッチなのしかなさそうだけど、単純におもしろい気もするので書いておく。
やってることは
なので、macbookである必然性はそこまでない。windowsでもcygwinとかで頑張れば行けそうな気がする。
環境:
macbook: OSX 10.9.5 (Mavericks)
iphone: iOS 9.0.2
手順
1. macbookのwebserver機能を有効にする。
Mavericksならターミナルから
$ sudo apachectl start
でよい。それ以外なら参考URLの一番上でいけるはず。
とりあえずmacbookのブラウザから http://localhost にアクセスして it Works! というページが出るのを確認する。
iPhoneから繋がることを確認する。ターミナルで
$ hostname
として nos-MacBook と出た場合iphoneのsafariで http://nos-macbook.local と打つと it Works! が出るはず。出なければipアドレスを使うなりすればいいはず。
2. iPhone用のプロキシ設定ファイルを置く
/Library/WebServer/Documents/proxy.pac として
function FindProxyForURL(url, host) { return "SOCKS nos-macbook.local:XXXX"; }
を置く*1。アドレスは手順1.で確認したやつ。XXXXは9999とか適当な数字にする。あんまり小さかったり大きかったりすると多分失敗する。
3. iphoneのプロキシを設定
設定 → WiFi → 接続のiマークを押すと詳細設定に行く。一番下のHTTPプロキシのところで自動を選ぶ。アドレスの欄に http://nos-macbook.local/proxy.pac を入れる。
手順4.でトンネルを用意する前のこの時点でも(プロキシなしで)ネットに繋がるが、これが"正しい"挙動なのかはよくわからない。
参考
*1:実体はjavascriptらしい
common lispのプログラムにemacsでアタッチ
cursesっていうターミナル表示の制御ライブラリ*1を使ってテトリスを作ってみてるのだけど、途中でエラーが起こったときにデバッグ表示が乱れる。特にマルチスレッドだと目も当てられない。そこでemacsでアタッチしてデバッグすると便利。
ちなみに↓がテトリスの動作風景。ソースはもう少し整理してから公開したいと思う。
開発してた時の手順
処理系はsbcl*2。quicklispを使っている。slimeもquicklispでインストールしてある。
環境は Mac (Mavaricks)、sbcl 1.2.11。
~/quicklisp/local-projects 以下にあるプロジェクトは設定なしに quicklisp で読めるので、
quickproject でこの下に作っている。mainという関数を作って export しておき、
$ sbcl --eval '(ql:quickload :tetris :silent t)' --eval '(tetris:main)' --quit
という感じで走らせる。もしくは
(defun make-tetris-command (&optional (name "tetris")) (sb-ext:save-lisp-and-die name :compression t :toplevel #'main :executable t))
というのを作っておくと
$ sbcl --eval '(ql:quickload :tetris :silent t)' --eval '(tetris:make-tetris-command)'
でカレントディレクトリに tetris という実行ファイルを作れるので
$ ./tetris
で走らせられる。
デバッグ
ターミナルから
$ sbcl --eval '(ql:quickload :tetris :silent t)'
で起動して
* (ql:quickload :swank)
To load "swank":
Load 1 ASDF system:
swank
; Loading "swank"
.
(:SWANK)
* (swank:create-server :port 5555 :style :fd-handler :dont-close t)
;; Swank started at port: 5555.
5555とする*3。ポートがかぶったりしてたら適当に変えてください。
emacs側で M-x slime-connect localhost 5555 とすると繋がる。
あとは emacs 側のreplで
CL-USER> (tetris:main)
と打てば、コンソール側でプログラムが走り出す。何か問題が起きたら emacs 側でデバッガが立ち上がる。M-n や M-p でスタックトレースを行き来したり、ソースファイルを開いてから修正して C-c C−c でコンパイルして実行再開とかできて便利。
しかし、ゲームを作るというのは、特に曲がりなりにも動くようになってからゲームとして成立しだす時というのはなかなか面白いね。
raspberry piの音周り
ラズベリーパイにしゃべらせたりairplayサーバーにしたりして遊んでいたメモを放出。
モデルはRaspberrry pi type B 512MB。
イヤホンジャック/HDMIから出力
参考: http://elinux.org/R-Pi_Troubleshooting#Sound
$ amixer cset numid=3
n はヘッドフォン端子が1、HDMIでは2、0で自動判別
テスト
$ aplay /usr/share/sounds/alsa/Noise.wav
ボリューム調節
$ amixer set PCM 20%
USBスピーカー
USB経由での方が音がよいという噂なので
サンワサプライ USBスピーカー ブラック MM-SPU6BK
- 出版社/メーカー: サンワサプライ
- 発売日: 2012/05/01
- メディア: Personal Computers
- この商品を含むブログを見る
電力を食うのか、繋いだ瞬間にraspberry piが落ちる。電源付きのUSBハブ越しでも落ちる。wifiのドングルはハブ越しなら平気だったのに…。
ハブにつないだまま起動すれば問題はない。それだけで認識はする(alsamixerコマンドで選択出来るようになっていた気がする)が、
$ sudo vi /etc/modprobe.d/alsa-base.conf
で
options snd-usb-audio index=-2
の行をコメントアウトして再起動するとデフォルトのサウンドデバイスになる。
音質はヘッドフォン端子から出すよりはよい、がまあそれなり、という感じ。ゲインを0にするとよくなるという噂を見てやってみたら(要は上でvolumeを100% にする。本体にも物理的なボリュームがあるのでそっちで調節)少しよくなった気がする?
あと1ヶ月くらい使ってるうちに音質悪くなくなってきた気がする(耳が慣れただけだ)。さすがにiphoneよりは良い音がする。
ゆっくり
ゆっくりの名で有名な AquasTalk のラズベリーパイ版。商用だが、個人・非営利利用ならライセンス料はかからない。http://www.a-quest.com/products/aquestalkpi.html からいただいたパッケージをどこかに解凍するだけ。解凍したディレクトリに行って
$ ./AquesTalkPi ゆっくりしていってね | aplay
とするとしゃべる。英単語もある程度読める。
オープンソースの OpenJTalk というのも試したのだが、商用だけあってこっちの方が読み上げが自然だし動作も軽快。
このままだとなんかエラーが発生したとき aplay に変なデータが流れ込んでザーッと言う音になって耳障りなので、シェルスクリプトを作った。
atalk:
#!/bin/bash aquestalkpi=/home/pi/build/aquestalkpi/AquesTalkPi var=`$aquestalkpi "$@" | base64; echo ":${PIPESTATUS[0]}"` ret=(${var##*:}) data=${var%:*} if [ $ret -eq 0 ]; then echo $data | base64 --decode --ignore-garbage | aplay -q else echo $data | base64 --decode --ignore-garbage exit $ret fi
というのを作って実行権限を付けてパスの通ったところに置いて、
$ atalk -g 15 ほげ
とする。オプションはもともとのAquesTalkと同じで、例えば上の -g は音量(100が最大)。なんも指定しないで走らせるともとの AquesTalkPi のヘルプを表示するのでそれを参照。
$ atalk NAME AquesTalkPi - Raspberry Pi用テキスト音声合成コマンド Ver.1.00 SYNOPSIS AquesTalkPi [-h] [-s speed] [-g volume] [-b] [-v f1|f2] [-k] [-t] [-o out.wav] [-f file | string] DESCRIPTION 日本語のテキスト音声合成を行います。 標準出力にWAV形式の音声データを出力します。 OPTIONS string 発声する文字列を指定します(UTF-8)。漢字も読めます。 音声記号列での指定も可能です(-k オプション) 音声記号列の詳細は、AquesTalk音声記号列仕様を参照。 スペース等を含む場合は ""で囲んで指定してください。 -f file 発声する文字列をファイルで指定するときに指定します。 先頭行だけしか処理しません。 stringを指定せず、且つ fileに - を指定したときは、 標準入力からの入力になります(pipe使用可能)。 -o file WAVファイルとして出力するときにファイル名を指定します。 指定しないときは標準出力に出力されます。 -t WAV形式音声データの代わりに音声記号列を出力します。 言語処理の結果を返します。 -k 発声する文字列が音声記号列の場合に指定します。 -v f1 | f2 声種を指定します。 f1:女声１(default) f2:女声２ -b 棒読み(アクセントが平板）になります。 -g volume 音量を指定します。 (0-100) default:100 -s speed 発話速度を指定します。 (50-300) default:100 -h このメッセージを表示します。 EXAMPLE $ ./AquesTalkPi 漢字も読めます。 | aplay $ echo ゆっくりしていってね？ | ./AquesTalkPi -b -f - | aplay $ ./AquesTalkPi -s 150 -v f2 -k -o out.wav "ファイルニ、シュツ'リョクシマ_ス。" *実行時にはAquesTalkPiと同じディレクトリに付属の/aq_dic以下が必要です。 LICENCE, etc. 商用利用、再配布には別途ライセンス契約(有償)が必要です。 ライセンス、その他につきましては、下記サイトを参照ください。 http://www.a-quest.com/products/aquestalkpi.html ---- COPYRIGHT 2013 AQUEST Corp. ----
AirPlay
参考: http://www.lifehacker.jp/2013/03/130306raspberry_piairplay.html
sudo apt-get install git libao-dev libssl-dev libcrypt-openssl-rsa-perl libio-socket-inet6-perl libwww-perl avahi-utils libmodule-build-perl git clone https://github.com/njh/perl-net-sdp.git cd perl-net-sdp perl Build.PL sudo ./Build sudo ./Build test sudo ./Build install cd .. git clone https://github.com/hendrikw82/shairport.git cd shairport make # ./shairport.pl -a AirPi sudo make install sudo cp shairport.init.sample /etc/init.d/shairport sudo chmod 755 /etc/init.d/shairport sudo update-rc.d shairport defaults sudo vi /etc/init.d/shairport #DAEMON_ARGS="-w $PIDFILE" #DAEMON_ARGS="-w $PIDFILE -a AirPi" sudo service start shairport
で、同じネットワークのiphoneとかから AirPi という名前で見えるようになるはず。
youtube
$ sudo apt-get install mplayer youtube-dl $ sudo apt-get install rtmpdump swftools libxml2-utils $ mplayer $(youtube-dl -g www.youtube.com/watch?v=Y6ljFaKRTrI)
NHKラジオ
第1
$ mplayer -playlist http://mfile.akamai.com/129931/live/reflector:46032.asx
時報・天気予報をしゃべらせる
応用としてcrontabに入れて定時でしゃべるようにしてみる。
時報
上で作った atalk が ~/bin 以下に置いてあれば
0 */1 * * * /home/pi/bin/atalk -g 15 `date +\%-H`'時だよ'
としておくと毎時時報をしゃべる。crontabだと % をバックスラッシュでエスケープしなければいけないらしいのでそこは注意。
天気予報
コード参考(というかほぼそのまま):
- http://nwpct1.hatenablog.com/entry/2013/10/18/175913
- http://www.programming-magic.com/20080820002254/
pythonでtenki.jpから天気予報の情報を読んで、atalkでしゃべらせられるようなテキストにする。
tenki:
#!/usr/bin/env python #coding: utf-8 import urllib2 import re import socket import htmlentitydefs def unescape_html_entity(text): reference_regex = re.compile(u'&(#x?[0-9a-f]+|[a-z]+);', re.IGNORECASE) num16_regex = re.compile(u'#x\d+', re.IGNORECASE) num10_regex = re.compile(u'#\d+', re.IGNORECASE) result = '' i = 0 while True: match = reference_regex.search(text, i) if match is None: result += text[i:] break result += text[i:match.start()] i = match.end() name = match.group(1) if name in htmlentitydefs.name2codepoint.keys(): result += unichr(htmlentitydefs.name2codepoint[name]) elif num16_regex.match(name): result += unichr(int(u'0'+name[1:], 16)) elif num10_regex.match(name): result += unichr(int(name[1:])) return result def tenki(url): socket.setdefaulttimeout(10.0) # 正規表現のパターンを定義 - タグ消し remove_tag = re.compile(r'<.*?>') try: htmldata = urllib2.urlopen(url) except urllib2.HTTPError as err: print('HTTPError') print(err) except urllib2.URLError as err: print('URLError') print(err) if isinstance(err.reason, socket.timeout): print('timeout') else: # print('Get HTML') pass content = htmldata.read() htmldata.close() content_list = content.split('\n') flag = False result = '' for line in content_list: if flag==True and not "</div>" in line: #正規表現でタグ消しした後にstripメソッドでスペースを除去。後ろに,をつけると改行しない result += remove_tag.sub("", line).strip() #print line #filewrite(line) if flag==True and "</div>" in line: result += remove_tag.sub("", line).strip() #print line #filewrite(line) flag = False elif "weatherCountryDescriptionBody" in line: result += remove_tag.sub("", line).strip() #print line #filewrite(line) flag = True print unescape_html_entity(result.decode('utf-8')).encode('utf-8') # print result if __name__ == "__main__": url = "http://tenki.jp/forecast/pref-14.html" tenki(url)
上は埼玉県の天気を出力するが、一番下のurlを変えたら他の県もいけるはず。
実行権限をつけて ~/bin において、
$ tenki
で埼玉県の天気がテキストで出力されるのを確認する。
そしたらcrontabに
30 8 * * * /home/pi/bin/atalk -g 15 '天気予報です、'`/home/pi/bin/tenki`
と書いておくと毎朝8時半に天気予報をしゃべってくれる。
pLaTeXのフォント埋め込み設定(Mac)
kindleは日本語のpdfをうまく表示出来ない(日本語のところだけ空白になる)ようで、フォントを埋め込んでやる必要がある。
ヒラギノを埋め込む手順が http://osksn2.hep.sci.osaka-u.ac.jp/~taku/osx/embed_hiragino.html にあったのでメモ。
現在の埋め込み設定の確認
$ kanji-config-updmap status
ヒラギノフォントの設定
sudo mkdir -p /usr/local/texlive/texmf-local/fonts/opentype/public/hiragino/ cd /usr/local/texlive/texmf-local/fonts/opentype/public/hiragino/ sudo ln -s "/Library/Fonts/ヒラギノ明朝 Pro W3.otf" HiraMinPro-W3.otf sudo ln -s "/Library/Fonts/ヒラギノ明朝 Pro W6.otf" HiraMinPro-W6.otf sudo ln -s "/Library/Fonts/ヒラギノ角ゴ Pro W3.otf" HiraKakuPro-W3.otf sudo ln -s "/Library/Fonts/ヒラギノ角ゴ Pro W6.otf" HiraKakuPro-W6.otf sudo ln -s "/Library/Fonts/ヒラギノ角ゴ Std W8.otf" HiraKakuStd-W8.otf sudo ln -s "/Library/Fonts/ヒラギノ丸ゴ Pro W4.otf" HiraMaruPro-W4.otf sudo ln -s "/System/Library/Fonts/ヒラギノ明朝 ProN W3.otf" HiraMinProN-W3.otf sudo ln -s "/System/Library/Fonts/ヒラギノ明朝 ProN W6.otf" HiraMinProN-W6.otf sudo ln -s "/System/Library/Fonts/ヒラギノ角ゴ ProN W3.otf" HiraKakuProN-W3.otf sudo ln -s "/System/Library/Fonts/ヒラギノ角ゴ ProN W6.otf" HiraKakuProN-W6.otf sudo ln -s "/Library/Fonts/ヒラギノ角ゴ StdN W8.otf" HiraKakuStdN-W8.otf sudo ln -s "/Library/Fonts/ヒラギノ丸ゴ ProN W4.otf" HiraMaruProN-W4.otf sudo mktexlsr sudo updmap-sys --setoption kanjiEmbed hiragino kanji-config-updmap hiragino
ヒラギノを埋め込む
$ kanji-config-updmap hiragino
埋め込みの停止
$ kanji-config-updmap nofont
mpcを使ってみた
mpc*1はcommon lisp用のパーサコンビネータライブラリ。
Common Lisp で Lispインタプリタを作ってみた – さくらんぼの技術備忘録などを見て気になったので使ってみた。
使用例
以下はマニュアルにあるメールアドレス < user >@< host > をパースする例。
まずパッケージを定義。
(defpackage simple-address (:use :cl :mpc :mpc.characters)) (in-package :simple-address)
で、userやhostの部分に入っていい文字を定義する。
(defun =address-character () (=or (=satisfies #'alphanumericp) (=one-of '(#\- #\_ #\. #\+)))) ;; アルファベット、数字、-、_、.、+ が許される
これを使って求めるパーサーを実装する。
(defun =simple-address () (=let* ((user (=string-of (=address-character))) (_ (=character #\@)) (host (=string-of (=address-character)))) (=result (list user host)))) ;; =let* 中の _ は文字@を無視するという意味。結果を =result を使って返す。
あとは run 関数で走らせる事が出来る:
(run (=simple-address) "foo@example.com") ⇒ ("foo" "example.com") (run (=simple-address) "!!!@@@.com") ⇒ NIL
ごく簡単な電卓
練習として整数だけ(分数も結果的に扱えるのだが)対応のごく簡単な4則演算+括弧の電卓を作ってみた。
simple-math-parser.lisp
(defpackage simple-math (:use :cl :mpc :mpc.characters :mpc.numerals)) (in-package :simple-math) (defun =add-op () (=let* ((op (=one-of '(#\+ #\-)))) (=result (case op (#\+ '+) (#\- '-))))) (defun =mul-op () (=let* ((op (=one-of '(#\* #\/)))) (=result (case op (#\* '*) (#\/ '/))))) (defun =primary () (=or (=let* ((_ (=character #\()) (p (=add-expr)) (_ (=character #\)))) (=result p)) (=integer-number))) (defun =mul-expr () (=let* ((result (=primary)) (rights (=zero-or-more (=list (=skip-whitespace (=mul-op)) (=skip-whitespace (=primary)))))) (=result (loop for r in rights do (setf result (list (car r) result (cadr r))) finally (return result))))) (defun =add-expr () (=let* ((result (=mul-expr)) (rights (=zero-or-more (=list (=skip-whitespace (=add-op)) (=skip-whitespace (=mul-expr)))))) (=result (loop for r in rights do (setf result (list (car r) result (cadr r))) finally (return result))))) (defun =simple-math () (=prog1 (=add-expr) (=end-of-input))) ;; (run (=simple-math) "1+2-3*4/4") (defun main () (loop (princ "> ") (princ (eval (run (=simple-math) (read-line)))) (fresh-line))) (main)
こんな感じ↓
$ clisp
(略)
[1]> (ql:quickload :mpc)
To load "mpc":
Load 1 ASDF system:
mpc
; Loading "mpc"
(:MPC)
[2]> (load "simple-math-parser.lisp")
;; Loading file /User/nos/simple-math-parser.lisp ...
> 2*(3+4)
14
> 1+2/3
5/3
使ってみて
ちょっとはまったのは、=ifや=and、=prog2などがマクロでなく関数なので、引数が全て評価されてしまうという事。これが原因で再帰的なパーサーを書こうとすると意図しないスタックオーバーフローが発生する。例えば上の電卓で =primary を
(defun =primary () (=if (=character #\() (=prog2 (=character #\() (=add-expr) (=character #\))) (=integer-number)))
と書くとrunした際にスタックオーバーフローする。この=ifの条件部分が成功しようがしまいが"(=add-expr)"が評価されてしまうので、関数の間で呼び出し関係が循環しているとスタックがあふれるまで停止しないわけだ。 =let* がショートサーキットに振る舞うマクロなのでこっちを使えば回避出来るのだけど…。これマクロに置き換えて不都合あるかなあ。ちょっと暇なときに調べてみたい。
しかしその欠点を差し引いてもよくできたライブラリと思う。使っていてあとソースを見てみると思いのほか短くて驚く。
*1: http://mr.gy/maintenance/mpc/
*2:fork元のsmugはできないというのに…
