workbench
Emacs上で動くVimを実装する。
（宗教戦争の解決）
CやC++ではなく、C#やJavaを使う。
（ゴミ問題の解決）
PHPをdisらない。
（差別問題の解決）
ナンやライスにHaskellをつけて食べる。
（食糧問題の解決）
適切なクラス設計を行う。
（資源の再利用）
Objectのprototypeを拡張しない。
（環境汚染の回避）
マクロを作らない。
（雇用減少の防止）
松岡修造botを作らない。
（地球温暖化の抑制）
Googleの背景色を#000にするGreasemonkeyスクリプトを書く。
（電力問題の解決）
オライリーのPDF版を購入する。
（資源の節約）
プログラミングをやめる。
（書くの放棄）
このブログは証明できない。 こちらも 面白いぞ。久々に腹が痛くなったよ。
random
前回ちょと乱数の事について書いた。んでもって、ふと次のような疑問が浮かんできたんだ。
random() + random()のように乱数同士を足し算すると均等?
random()とrandom()×random()のどっちがランダムなの？
乱数同士を足し算したら、２倍ばらばらになるんじゃなかろうか。更に掛け算しちゃったら、それはもう ２乗で乱数度が向上するんではなかろうか？ 頭の中で、しばしシュミレーションしてみてください。 答えが分からなかったら、下記を走らせてみましょう。
import random
from pylab import *
from numpy import *
n = 200
ans = n*[0]
s = 1.0/n
def run():
for n in xrange(1000000):
a = random.random()
b = random.random()
#r = a
r = (a + b) / 2
#r = a * b
ans[int(r/s)] += 1
plot(ans)
show()
run()
０以上１未満の乱数を１００万回発生させて、その度数をグラフにするスクリプトです。適時ｒの所の コメントを付けたり外したりして、走らせてみてください。
結果は意外なものでした。思考力が落ちているんですかねぇ？ 更に面白いのが、上記では乱数を２回 しか足していないけど、これを３，４、と増やして行くと、だんだん あれ の形に似てくる。こんなの 想像の範囲外でしたよ。こういうの、お馬さんと勝負してる人は強いんですかね？ いやいや、半丁とか のサイコロの世界の人の方が強かったり。
timeit
で、調子に乗って、足す回数を増やして行くと、結果が表示されるまでの時間がどんどん増えて行く。腹時計では 正確に実行時間が計測出来ないので、Pythonの助けを借りて測ってみた。（事前にplotとshowはコメントアウト しておきます）
>>> from timeit import Timer
>>> Timer('run()', 'from %s import run' % __name__).timeit(1)
1.5477454441102978
あるいは、もっと素直に
>>> import time >>> def meas(): t1 = time.time() run() t2 = time.time() print t2 - t1 >>> meas() 1.59099984169
ちょいと一つの式の実行時間を計りたいような場合は
>>> import timeit
>>> timeit.Timer('[n for n in xrange(1000000)]').timeit(1)
0.11548197591810094
これを見ると、１００万回のループの空回しが、０．１秒強って事が分かる。run()は内部に100万回回るループ を内包してたからいいけど、単独の手続きだったらどうしよう。こんな場合でもtimeitを使うと、１００万回 ループさせて計測出来るようになっている。
>>> import timeit
>>> setup='''
def foo(x):
return x*x
'''
>>> stmt = 'foo(1)'
>>> t = timeit.Timer(stmt, setup)
>>> t.timeit()
0.33254254381492616
>>>
>>> t.repeat(3, 1000000)
[0.33528088067059088, 0.33709172534497895, 0.35431430530975234]
repeatを使うと、ｎ回測定した結果を表示してくれるので、プチ嬉しかったりします。
んでもって、timeitってエバッてるなと思ったんだ。コードを覗いてみる鹿。
template = """
def inner(_it, _timer):
%(setup)s
_t0 = _timer()
for _i in _it:
%(stmt)s
_t1 = _timer()
return _t1 - _t0
"""
:
src = template % {'stmt': stmt, 'setup': setup}
:
code = compile(src, dummy_src_name, "exec")
exec code in globals(), ns
eval使ってなかった。python流のマクロだな。テンプレートに実行したい命令を埋め込み、それを コンパイルしてから、実行。あれ？ 名前の衝突はどう処理してるのかな？
仮にstmtが、_t0 = 'abc' なんてのだったりすると、時間から文字列を引き算する事になっちゃって、 エラーになってしまうと思うんだけど。。。 一般ユーザーは、アンダーバーから始まる変数は使うな と言う、暗黙の了承事項を受け入れなければいけないのかな？ それとも、隔離した環境で動かす ような仕組みがあって、その中でstmtを動かすから、衝突は無いのかな。もう少し精査しないと 分からんな。
workbench
毎度の事で恐縮だけど、前回ワークベンチと称して
>>> w = wave(1, 1) >>> f = fft(w) >>> t = gett(f) >>> plot(t) >>> show()
こんな事をやった。波形を生成して、それをFFTして、各周波数における位相を求めて、それを表示させる。 いちいち変数を引き回していて、めんどうくさい。
これを、何とか数珠繋ぎのように出来ないか？ 実世界の実験台の上には、発振器があり自作のFFT演算器が あり、位相の計算機があり、パソコンなりの表示機がある。それらを配線して、一つの機能を実現する。
フィルターの伝達特性を見たかったら、発振器（ホワイトノイズ発生器の方がいいかも）とFFT演算器の間にフィルターを入れてやれば良い。それらを 繋ぐのは、コードと言うかケーブルと言うかワイヤーと言うか。要するに信号が流れる線だな。
この発想で直ぐに思いつくのは、unixのpipe機能だ。イメージはこんな感じかな。
wave(1,1) | fft() | gett() | plot() | show()
信号が流れて行くってのを強調したいのなら、＞ でも －＞ でも、Haskell風に ＞＞＝ でもいいけど、 やっぱり、日ごろから、慣れ親しんでいる、バーチカルバーがいいな。
で、この機能をpythonでどうやって実現する？ 各装置間を結ぶ線の代わりに変数を使おう。 evalが役に立つか、ちょっと調べてみる。
>>> eval('x = 3 + 4')
Traceback (most recent call last):
File "<pyshell#32>", line 1, in <module>
eval('x = 3 + 4')
File "<string>", line 1
x = 3 + 4
^
SyntaxError: invalid syntax
evalの中では、変数のbindは出来ないのね。しゃーない、これはどうだ？
>>> z = eval('3 * 4')
>>> z
12
evalの評価値を変数にbindするのは、OKだった。 これを踏まえて、簡単なパイプを書いてみよう。ああ、パイプとは図々しいので、土管ぐらいかな。 だったら、頭文字を取って、パイプの起動コマンドは、do ぐらいでいいか。haskellにも、それっぽい のが有ったようなので、止めておこう。ren（連）でもいいけど、無難に run()ぐらいかな。
from pylab import *
from numpy import *
from numpy.fft import *
from math import *
from cmath import phase
SMPL= 4096 # sampling points, must be 2**n
r2d = lambda r: 180 * r / pi # radian to degree
d2r = lambda d: d * pi / 180 # degree to radian
th = lambda z: r2d(phase(z)) # angle
getr = lambda a: array([x.real for x in a]) # real
geti = lambda a: array([x.imag for x in a]) # imag
geta = lambda a: array([abs(x) for x in a]) # power
gett = lambda a: array([th(x) for x in a]) # phase
def dt(a): # display time domain
plot(a)
show()
def df(a): # display freq domain
bar(arange(SMPL/2), a[:SMPL/2])
show()
def st(fs1=3, g=1.0): # single tone (for audio)
wave = []
for n in range(SMPL):
wave.append(g * cos(2 * pi * fs1 * n / SMPL) )
return array(wave)
def tt(fs1=3, fs2=5, g=0.5): # two tone (for audio)
wave = []
for n in range(SMPL):
wave.append(g * (cos(2 * pi * fs1 * n / SMPL) +
sin(2 * pi * fs2 * n / SMPL)) )
return array(wave)
def am(fs, fc=300, m=0.5): # amplitude modurator
for n in range(SMPL):
fs[n] = (1 + m * fs[n] ) * cos(2 * pi * fc * n / SMPL)
return fs
def run():
try:
exlude = set(['show()',])
cmds = wire.split('|')
for n, cmd in enumerate(cmds):
cmd = cmd.strip()
if (n == 0) or cmd in exlude:
do_cmd = (cmd if cmd[-1] == ')' else (cmd + '()' ))
else:
do_cmd = (cmd.replace('(', '(val,', 1) if cmd[-1] == ')'
else (cmd + '(val)' ))
print(do_cmd)
val = eval(do_cmd)
except:
print("please define 'wire' before run(), such as ...")
print(" wire = 'tt | am(m=0.7) | fft | geta | df'")
使い方は、こんな感じ。
>>> run() please define 'wire' before run(), such as ... wire = 'tt | am(m=0.7) | fft | geta | df' >>> >>> wire = 'tt | am(m=0.7) | fft | geta | df'
いきなりrun()したので、怒られました。まず最初に実体配線図を 定義します。 そして通電。
>>> run() tt() am(val,m=0.7) fft(val) geta(val) df(val)
マイクの代わりにオーデォ帯のツートーン発信器を用意し、それをAM送信機に接続しました。送信機は、 搬送周波数が、標準で300、変調度50%になってますが、変調度だけ変更しました。そして、送信機の 出力をFFT解析し、そこからパワースペクトラムを取り出して、周波数ドメイン表示してます。
尚、上記のようにデフォルト値から外れた状態で動かしたい場合のみ括弧で括って設定します。デフォルト値で 良い場合は、括弧は不要な、もの草仕様になってます。
run()時の信号の流れが表示されてます。定義した内容が一部書き換えられてるって取ると、runは マクロと思われます。自動的にmacro-expand1をやってますしね。また、run()をインタープリターと理解 すれば、実行トレースが表示されてると思っても良いでしょう。
>>> wire = 'st | am | getb | dt'; run() st() am(val) getb(val) please define 'wire' before run(), such as ... wire = 'tt | am(m=0.7) | fft | geta | df' >>> wire = 'st | am | geta | dt'; run() st() am(val) geta(val) dt(val)
今度は、シングルトーンで変調をかけて、その波形を検波してからオシロスコープで眺めようと算段です。最初は 名前を間違えてしまい最後まで実行されませんでした。例外の扱いをもう少し丁寧にしておくべきだな。まあ、ともかく波形が 表示されたら、マウスでぐりぐりやって、波形を拡大する事も出来ます。
アナログ・オシロスコープの高級なやつには、B掃引が付いていて、波形の一部を拡大表示できたけど、 今なら、こんなの朝飯前なのね。写真撮影機能も付いてるし。。。昔は証拠写真をポラロイドカメラで 取ったものです。
これをどんどんと発展させていくと、去年の３月頃試してみた 'puredata'に行き着く訳です。こちらは、 本格的で、マウスでぐるぐりと配線したり、リアルタイムにパラメータを変えるなんて事が容易に 出来、GUIな人には嬉しい事でしょう。
環境違い？
上記のお遊びは、IDLEの上でやってたんだけど、ipythonの上へ持ってくると、怒られた。
In [1]: from wb import * In [2]: wire = 'st | dt' In [3]: run() please define 'wire' before run(), such as ... wire = 'tt | am(m=0.7) | fft | geta | df' In [4]: wire Out[4]: 'st | dt'
ちゃんと配線してるのに、配線してないと言い張っている。読み込んだ時の状態が 凍結されてるっぽい。おかしいなあ。で、ちょと実験。
>>> def hoge(): ... print(aa) ... >>> hoge() Traceback (most recent call last): File "<stdin>", line 1, in <module> File "<stdin>", line 2, in hoge NameError: global name 'aa' is not defined >>> aa = 'fuga' >>> hoge() fuga >>> aa = 'foo' >>> hoge() foo
適当な関数を定義した後、走らせると、そりゃ、関数内で使ってる変数が無いので、上位を見に行く けど、無いのでエラー。
次に、関数の外で変数を定義しておけば、ちゃんと動く。どこが違うんだろう？
上の実験に習って、run()をコピペして再定義したら(REPL側での定義ですな)、その時点から動き始めた。どうやら、importを使って 読み込んだやつは、安全のために、凍結されちゃうのね。動きからクロージャだな。で、回避策は、
ipython -i wb.py
なのね。ipython上で使うと、Ctrl-p で、過去の配線を呼び出して、ちょっとパラメータを変更して実行 なんてのが、手軽に出来るので、とっても便利です。