■
最近読んで面白かった話
High Energy Asymptotics of Multi--Colour QCD and Exactly Solvable Lattice Models
http://arxiv.org/abs/hep-th/9311037
High energy QCD as a completely integrable model
http://arxiv.org/abs/hep-th/9404173
High-Energy QCD as a Topological Field Theory
http://arxiv.org/abs/hep-ph/9807451
Integrability of High-Energy QCD
http://ci.nii.ac.jp/naid/110001208781
元々、高エネルギー極限(Regge極限)に於けるQCDを記述する有効理論として、Reggeon field theoryというのが70年代からあって、それは本質的に二次元の場の理論と同等らしい(このへんは、物理的な議論なので、あんまり理解してない)のだけど、更にlarge N極限を取ると、波動関数が正則部分と反正則部分に分かれるという、ありがちなことが起きて、正則部分だけ解けばよくなる。それで、正則部分のハミルトニアン(ハミルトニアンも正則部分と反正則部分に分解する)を見てやると、実は、XXX模型のハミルトニアンになっていることが分かる
但し、普通のHeisenberg模型では、各格子点に載ってるHilbert空間は、su(2)の二次元表現であるけども、今の場合は、各格子点にSL(2,C)の既約ユニタリ表現(主系列表現)が載っている。格子点の数は、丁度reggeon(reggeized gluonとかいうquasiparticle)の数に対応している。格子点に載ってるHilbert空間は違うのだけど、同じようにR-matrixとかL行列は決定できて、普通にXXX模型を解く手順を踏めば、固有値・固有ベクトルの計算はBaxter方程式に帰着し、この系は解ける(というのは"嘘"で、Baxter方程式が一般に解けないので、固有値とか厳密には計算できないのだけど、ここまで出来ることを可積分と呼んでいる)
最終的に、一次元格子模型に帰着させてはいるものの、元の理論はちゃんと4次元時空で定義されていて、超対称性を仮定したりもしてないし、最初から二次元時空上の理論であったりもしない点で、他の"解けるtoy model"と比べると、一番現実的な模型で一線を画している気がする。
コヒーレント状態とウェーブレット変換(3) Cauchy waveletとaffine coherent state
少し勘違いしてたことなどを補足
coherent state物語
http://d.hatena.ne.jp/m-a-o/20130303#p2
コヒーレント状態の有限群での類似物
http://d.hatena.ne.jp/m-a-o/20130917#p1
が(1)と(2)のつもり
Short-time Fourier変換(STFT)の場合、Gauss窓を使うのは、Weyl-Heisenberg群のコヒーレント状態(boson coherent stateと呼ぶことにする)を与え、この場合は、Gabor変換と呼ばれて、"数学的には"、最も自然なSTFTであると言える。連続ウェーブレット変換(CWT)の場合、背後にある群は、ax+b群であって、この場合は、コヒーレント状態に相当するものがないので、Gabor変換のような特別によいCWTは存在しないのだと思っていた。が、実は、Cauchy waveletというものが存在していて、これは、ax+b群のコヒーレント状態(affine coherent stateと呼ぶ)と考えられているらしい。つまり、
CWT : STFT = ax+b群 : Weyl-Heisenberg群 = Cauchy wavelet : Gabor変換
という対応関係がある。
コヒーレント状態の理論は、コンパクト半単純連結Lie群に対しては、もはや最小不確定性のような特徴づけは出来ないけれど、非常に綺麗に出来ている(コヒーレント状態は、余随伴軌道/旗多様体の小平埋め込みの具体的な実現を与え、それは最高ウェイトベクトルに群を作用させて得られる軌道と一致する)のに対して、一般の非コンパクト群については、少し怪しいものとなる。
非コンパクト群の場合、コヒーレント状態の集合は射影化した表現空間に埋め込むようなことはできないし、Weyl-Heisenberg群の場合、単純に真空ベクトルに群を作用させると、余計な位相因子が付いてきてコヒーレント状態をはみだす。アフィン・コヒーレント状態の場合は、コヒーレント状態の集合というのは一意には決まらない。それでも、あるベクトルの集合がコヒーレント状態と見なされるための強い条件として、幾何学的量子化によって対応する表現が得られる余随伴軌道と同型で、余随伴軌道には自然なケーラー構造が入ってるということがある。Weyl-Heisenberg群の場合、対応するケーラー多様体は複素平面で、ax+b群の場合は、上半平面(z=b+iaによって、座標を入れることができる)となる
#一般的に正しいのかどうか知らないけど、多分正しいこととして、以下のようなことがある。Lie群の表現を微分すればLie環の作用も得られることに注意すると、(コンパクトとも半単純とも限らない)Lig群Gの一般化コヒーレント状態|z>は、GのLie環の元Jに対して、
となって、
なので、座標(y_0,y_1,y_2)は二葉双曲面の連結成分と対応し(y_0>0に注意)、単位円板と同型である。更に、この余随伴軌道の幾何学的量子化によって、元の既約表現を復元できる(Kostant-Kirillovの軌道法)ので、コヒーレント状態は"最も古典的な状態"とか言われるけど、実際にある意味で量子化の逆を与える(まぁ、物理としては、状態のHilbert空間が、適当なLie群の(二乗可積分な)既約表現空間になっているというのは結構特殊な状況である気もするけど)。関数環の"量子化"を得たければ、普遍展開環の既約表現のannihilator(原始イデアルと呼ぶ)による商環を考えればよい？記号とか計算の詳細は、Perelomovの原論文
Coherent states for arbitrary Lie group
http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.cmp/1103858078
や
Some basics of su(1,1)
http://arxiv.org/abs/quant-ph/0407092
を参照。別のよく知られた一般化コヒーレント状態として、su(2)コヒーレント状態(スピン・コヒーレント状態)があるけども、この場合は、余随伴軌道は球面と同型で、全ての有限次元既約表現は、余随伴軌道の幾何学的量子化で得られる
geometric quantization of the 2-sphere
http://ncatlab.org/nlab/show/geometric+quantization+of+the+2-sphere
Cauchy waveletの具体的な形は
Time-frequency localisation operators - a geometric phase space approach: II. The use of dilations
http://www.math.duke.edu/~ingrid/publications/Inverse_Prob_4_p661.pdf
のAppendixに書いてある(ここに書いてあるのは、より一般的な何かで、特にγ=1の時が、Cauchy waveletで、βごとに異なるコヒーレント状態系が得られる)。Mother waveletの形は、Poisson分布の確率密度関数を思い出すけども、Poisson分布の時と、変数と定数が入れ替わってる
歴史的には、Klauderがaffine coherent stateの考案者らしい。
Path Integrals for Affine Variables
http://link.springer.com/chapter/10.1007%2F978-1-4615-7035-6_7
Unitary Representations of the Affine Group
http://jmp.aip.org/resource/1/jmapaq/v9/i2/p206_s1?isAuthorized=no
これは、Morletによるウェーブレット解析の発表より早い。Cauchy waveletという命名は、Holschneiderという人の
Wavelets: An Analysis Tool
http://www.amazon.co.jp/dp/0198505213
という本によるらしい(1995年)。何でCauchyなのかよく分からないのだけど、この名前は一般的になっているっぽい。
#別のアプローチとしてboson coherent stateの最小不確定性状態という定義を直接+b群に持ち込もうという試みが
The Affine Uncertainty Principle in One and Two Dimensions
http://www.sciencedirect.com/science/article/pii/0898122195001085
でなされているけど、多分あまりよい答えではないと、個人的には思う
ax+b群は、SU(1,1)に埋め込むことができ、ax+b群の既約ユニタリ表現は、SU(1,1)の離散系列表現をax+b群に制限することで得られる。Cauchy waveletとPerelomovによるSU(1,1)コヒーレント状態は実は、適当な関数空間で見ると、同じ形をしている
Characterization of SU(1,1) coherent states in terms of affine group wavelets
http://arxiv.org/abs/math-ph/0209006
それで数学的なことは、ともかく、実際にCauchy waveletはMorlet waveletとかよりよかったりするのかというと、謎。あんまり使ってる人いないし。
おまけ:というか、STFTでも、Gabor変換と他の窓関数使った場合で、どういう差があるのだろう
窓関数
http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/basic/chap9/index.htm
とか見ると、ガウス窓は、メインローブの幅がでかいけど、サイドローブの振幅が小さくて、"主成分の周波数分解能は劣りますが、小さい電力のスペクトルを検出するのに優れた特性を示します"とか書いてある。
最初何を言ってるか分からなかったけれど、メインローブの幅というのは、要するに単純な正弦波をSTFTした時の分散の大きさのことっぽい。これは矩形窓でやるのがベストに決まっていて、十分に大きな矩形窓を取れば、無限区間でFourier変換するのとほぼ同じ結果(デルタ関数的な何か)に近づいていく。完全に周期的な波では、サイドローブとか出てこないのだけど、現実の音とか信号は、減衰していくのが普通なので、これがサイドローブとか生み出す要因になるっぽい(直接的には、矩形関数のFourier変換がsinc関数なので、それがサイドローブの元であるらしい)。サイドローブがあると、複数個所にピークもどきが出てきてしまうので、あまり嬉しくない
Gauss窓だと、そんなにメインローブの幅が大きくなるのだろうかと思って、適当に実験して見た
#!/usr/bin/env python
import numpy as np
import scipy as sp
from scipy import signal
from math import sqrt,pi,sin
"""
import wave
def wavread(filename):
wf = wave.open(filename , "r")
fs = wf.getframerate()
x = wf.readframes(wf.getnframes())
if wf.getsampwidth()==2:
x = np.fromstring(x,dtype=np.int16)/32768.0
else:
x = np.fromstring(x,dtype=np.int8)/128.0
wf.close()
return x,fs
"""
fs = 44100
sinewave = np.array([sin(2*pi*n*400/fs) for n in xrange(fs*2)])
window_func = {"Hamming" : (lambda sz:sp.hamming(sz)) ,
"Hanning" : (lambda sz:sp.hanning(sz)) ,
"Rectangular" : (lambda sz:sp.kaiser(sz,0)) ,
"Kaiser(alpha=3)" : (lambda sz:sp.kaiser(sz,3)) ,
"Gaussian(window size=2 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/2.0)/(sqrt(2*pi)*sz/2.0)) ,
"Gaussian(window size=4 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/4.0)/(sqrt(2*pi)*sz/4.0)) ,
"Gaussian(window size=5 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/5.0)/(sqrt(2*pi)*sz/5.0)) ,
"Gaussian(window size=7 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/7.0)/(sqrt(2*pi)*sz/7.0)) ,
"Gaussian(window size=10 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/10.0)/(sqrt(2*pi)*sz/10.0)) ,
"Gaussian(window size=20 sigma)" : (lambda sz:signal.gaussian(sz,std=sz/20.0)/(sqrt(2*pi)*sz/20.0))}
if __name__=="__main__":
for (name , wfun) in window_func.iteritems():
start = int(0.5*fs)
framesize = int(0.04*fs) #-- 40msec
w = wfun(framesize)
xs = sp.fft(w*sinewave[start:start+framesize])
total = sum([abs(v)**2 for v in xs])
avg = sum([min(n,framesize-n)*abs(v)**2 for (n,v) in enumerate(xs)])/total
delta = sum([(min(n,framesize-n)-avg)**2*abs(v)**2 for (n,v) in enumerate(xs)])/total
print "avg=%f[Hz] delta=%f name=%s" % (avg*fs/framesize,delta,name)
やってることは、400Hzの正弦波に色んな窓関数でSTFTしてるだけ(窓幅は40ms固定)。出力は
avg=399.999999[Hz] delta=0.333712 name=Hanning avg=400.000949[Hz] delta=0.193046 name=Gaussian(window size=4 sigma) avg=400.000000[Hz] delta=1.266515 name=Gaussian(window size=10 sigma) avg=400.002197[Hz] delta=0.137814 name=Kaiser(alpha=3) avg=399.999999[Hz] delta=0.266531 name=Hamming avg=400.000000[Hz] delta=5.066059 name=Gaussian(window size=20 sigma) avg=400.000005[Hz] delta=0.620575 name=Gaussian(window size=7 sigma) avg=400.000658[Hz] delta=0.024861 name=Gaussian(window size=2 sigma) avg=400.000318[Hz] delta=0.314517 name=Gaussian(window size=5 sigma) avg=400.000000[Hz] delta=0.000000 name=Rectangular
とかなって、deltaが分散の大きさを表す(平方根とって適当にスケールすれば周波数の単位に変換できる)けど、正弦波相手だと、圧倒的に矩形窓最強っぽい。次によいのが、Gaussian(window size=2 sigma)というやつだけど、例えばsignal.gaussian(10,std=5)はarray([ 0.66697681, 0.78270454, 0.8824969 , 0.95599748, 0.99501248,0.99501248, 0.95599748, 0.8824969 , 0.78270454, 0.66697681])とかで、窓の端の値がピーク値の0.6倍くらいもあるので、ガウス窓といっても、積分区間が明らかに不十分。
signal.gaussian(10,std=2)だとarray([ 0.07955951, 0.21626517, 0.45783336, 0.7548396 , 0.96923323, 0.96923323, 0.7548396 , 0.45783336, 0.21626517, 0.07955951])で、scipy.hamming(10)がarray([ 0.08 , 0.18761956, 0.46012184, 0.77 , 0.97225861,0.97225861, 0.77 , 0.46012184, 0.18761956, 0.08 ])なので、これくらいだと、Hamming窓に近い感じがある。signal.gaussian(10,std=10.0/7)はarray([ 0.00700417, 0.04972487, 0.21626517, 0.57622907, 0.94058806, 0.94058806, 0.57622907, 0.21626517, 0.04972487, 0.00700417])でよさげではあるけど、このくらいになると、分散が他の窓関数より大きい
なんか別にGabor変換が殊更よいというわけでもないっぽい。。まぁ、そもそも何故Weyl-Heisenberg群やax+b群を選ぶのかという問いにも明確な答えがあるわけでもないと思うし
Huモーメント不変量による類似画像検索
Googleが"画像で検索する"機能をリリースして、暫く経つ。他の人はどうか知らないけど、自分の場合、実際に使うときって、"同じ"画像を探すのにしか使わない。類似性というのは、色んな視点がありえるので、正しい答えというのもなく、そういう意味で難しい問題であるけど、同じかどうかは人間ならほぼ100%判定できる。Googleは、この機能を実現するのにDeep何とかいう手法を使って沢山のマシンで学習させたとか、どっかに書いてたけど、同じ画像を探すという目的に限定すると、もっと手軽な手法がありそうなもんではある。
同じと言っても、フォーマットが違ったり、サイズが違うと、拡大縮小やら圧縮・符号化アルゴリズムによって、データとしては違うものになってしまうので、普通のハッシュ関数に通しても仕方ない。OpenCVのcv::matchTemplateのようなものを使うという手もあるけど、サイズを合わせないといけないので、2つのデータの類似度が対称にならないし、計算量も多い。というわけで、この問題もそれほど自明ではないけど、もし、画像データから、拡大縮小(スケール変換)、並進・回転の操作によって不変な量をいくつか計算することができれば、類似度/距離の計算を簡単化・高速化できるし、データ量も大幅に減らすことができる(可能性がある)
丁度いいことに、画像のスケール不変なモーメントの多項式として、並進・回転不変な7つのHuモーメント不変量というものがよく知られていて、OpenCVにも実装されているので、簡単に使える。
Image moment
http://en.wikipedia.org/wiki/Image_moment
#上Wikipediaの記事によると、Huのoriginalのセットは、独立でない元を含むばかりか、3次までの独立な元を全て含んでないとある。OpenCVの実装は、現時点では、Huのoriginalのセットに基づいているっぽい
HuMoments
http://docs.opencv.org/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html#humoments
Huモーメントから、データ間の類似度を測る方法は、特に自然なやり方とかはないけど、cv::matchShapesとかで使われてるやり方を参考にすればいいんだと思う(ただ、ここのやり方は高次モーメントの重みがやや高くなる傾向にあって、それほどよくないかもしれない)
matchShapes
http://docs.opencv.org/modules/imgproc/doc/structural_analysis_and_shape_descriptors.html#matchshapes
この方法では、同じ画像(拡大縮小したり、平行移動したり回転したもの)に対して、同じ不変量の値を取ることが保証されているけど、逆に、全然違う画像がたまたま近い不変量の値を取ってしまう可能性は排除されない。つまり、大規模なデータに対して、この方法で検索すると、再現率は保証できるけど、適合率はどうなるか分からない。適合率の低下を回避するには、更に多くの独立な不変量を持ってくる(Jan Flusserの論文では、4次までの独立な不変量が計算されている)とか、あるいは、別の方法(例えば、cv::matchTemplateとか)でad hocにフィルタリングしていくとか、色々やり方は考えられる。
一応、そもそも「(無限に)沢山の不変量を持ってきても原理的にすら区別不可能な"異なる"画像というのはありえるのか?」という問いがある。究極的には、定数だってスケール変換・回転・並進の不変量であるわけで、不変量を全部持ってきても、識別不能な画像が大量にあるようでは、これは使えない方法ということになってしまう。それで、まず問題となるのは、画像のモーメント(無限個ある)は、どれくらい画像を区別できるのかという話であると思う。わたしは正しい答えを知らないけど、異なる画像で全てのモーメントが等しくなることは滅多にないだろうとは思う(厳密には、モーメント問題とかいうテーマに属する問題な気がする。一次元の有限区間に対しては解かれているはず)。なので、十分に多くの不変量を持ってくれば、全然違う画像が原理的にすら区別できないという最悪の事態は避けられるのだと思う。まぁ、細かいことをいうと、不変量が、全ての"軌道"を区別可能であるとは限らないとか、色々あるし、全然違う画像が、全ての不変量で近い値を持ってしまうことがないかは別問題ではある
この方法は、3次元形状やボクセルデータにも拡張可能。動画にも使えるかどうか分からないけど、原理的には使えるはず。まぁ、理屈は純数学的なので、簡単で分かりやすい
一応、簡単なテスト
import cv2
from math import log
def hudist(f1 , f2):
def h2m(x):
if x==0:
return 0.0
elif x>0:
return 1.0/log(x)
else:
return -1.0/log(-x)
def hum(f):
img = cv2.imread(f , cv2.CV_LOAD_IMAGE_GRAYSCALE)
assert(img!=None),("failed to open %s" % f)
moments = cv2.moments(img)
return cv2.HuMoments(moments)
hu = hum(f1)
hu2 = hum(f2)
d = sum([abs(h2m(x1)-h2m(x2)) for (x1,x2) in zip(hu.T[0],hu2.T[0])])
return d
if __name__=="__main__":
imgs = ['Lenna.png','Lena3.jpg','sendo1.jpg','Sendou2.jpg','monnalisa.jpg','apple.jpg','Red_Apple.jpg']
for f1 in imgs:
print f1,
for f2 in imgs:
v = 100*hudist(f1,f2)
print ("%5.2f" % v),
print ""
結果
| Lenna.png | Lena3.jpg | sendo1.jpg | Sendou2.jpg | monnalisa.jpg | apple.jpg | Red_Apple.jpg | |
| Lenna.png | 0.00 | 0.04 | 5.93 | 5.87 | 4.19 | 2.43 | 4.83 |
| Lena3.jpg | 0.04 | 0.00 | 5.93 | 5.86 | 4.22 | 2.41 | 4.83 |
| sendo1.jpg | 5.93 | 5.93 | 0.00 | 0.09 | 7.82 | 4.55 | 1.19 |
| Sendou2.jpg | 5.87 | 5.86 | 0.09 | 0.00 | 7.86 | 4.47 | 1.11 |
| monnalisa.jpg | 4.19 | 4.22 | 7.82 | 7.86 | 0.00 | 3.97 | 7.88 |
| apple.jpg | 2.43 | 2.41 | 4.55 | 4.47 | 3.97 | 0.00 | 4.94 |
| Red_Apple.jpg | 4.83 | 4.83 | 1.19 | 1.11 | 7.88 | 4.94 | 0.00 |
テストに使用した画像は以下の通り。まぁ、そこそこいけそうな感じではある。
おまけ:Huモーメント不変量が代数的に独立でないことは、グレブナー基底を使って簡単に確認できる。具体的には、以下のようなRisa/Asirコードを書く
load("gr");
def relation(Exprs , Vars , NewVars){
/* assert(length(Exprs)==length(NewVars)); */
S = [];
for(I = 0 ; I < length(Exprs) ; I++)S = cons(Exprs[I] - NewVars[I] , S);
Basis = nd_gr(S , append(Vars , NewVars) , 0 , [[0, length(Vars)] , [1,length(NewVars)]]);
Com = [];
for(I = 0 ; I < length(Basis) ; I++){
Check = 1;
for(J= 0 ; J < length(Vars) ; J++){
if( deg(Basis[I] , Vars[J]) > 0 ){
Check = 0;
break;
}
}
if(Check == 1)Com = cons(Basis[I] , Com);
}
return nd_gr(Com , NewVars , 0 , 0);
}
/* Huモーメント不変量 */
I1=m20+m02;
I2=(m20-m02)^2+4*m11^2;
I3=(m30-3*m12)^2+(3*m21-m03)^2;
I4=(m30+m12)^2+(m21+m03)^2;
I5=(m30-3*m12)*(m30+m12)*((m30+m12)^2-3*(m21+m03)^2)+(3*m21-m03)*(m21+m03)*(3*(m30+m12)^2-(m21+m03)^2);
I6=(m20-m02)*((m30+m12)^2-(m21+m03)^2)+4*m11*(m30+m12)*(m21+m03);
I7=(3*m21-m03)*(m30+m12)*((m30+m12)^2-3*(m21+m03)^2)-(m30-3*m12)*(m21+m03)*(3*(m30+m12)^2-(m21+m03)^2);
relation([I1,I2,I3,I4,I5,I6,I7],[m11,m20,m02,m12,m21,m30,m03],[p1,p2,p3,p4,p5,p6,p7]);結果は、[-p4^3*p3+p5^2+p7^2]とかなって、I5^2+I7^2=I4^3*I3であることが分かる。これは、Flusserの論文
On the independence of rotation moment invariants
http://library.utia.cas.cz/prace/20000033.pdf
に書いてある結果(Flusserは、よりad hocな手法で、これを発見している)。ついでに、スケール不変なモーメントは無限個あるけども、有限次数で打ち切れば(例えば、Huモーメント不変量は三次までしか見てない)、それはただの(合同変換群による)不変式の計算に帰着して、やっぱりグレブナー基底を使って、独立な生成元を機械的に計算することができる(はず)。
合同変換群による不変式を計算するというのは、本質的には
三角形の合同条件と不変式論
http://d.hatena.ne.jp/m-a-o/20130104#p1
でやったのと同じようなもの(これ書いた当時は、不変式を計算するアルゴリズムを知らなかったけど、Derksenのアルゴリズムを使えば計算できる)。というわけで、これはComputer Vision分野におけるグレブナー基底のちょっと面白い応用事例でもあると思う
分子動力学での融点・沸点の計算
アルゴンの分子動力学計算
http://d.hatena.ne.jp/m-a-o/20130917#p3
の疑問について。
アルゴン原子集団を、固体から加熱していった時と、気体から冷却していった時で、融点・沸点がずれるというのを書いたけど、過冷却とか、そのへんの話らしい。過冷却？とかいうのは思ったんだけど、このへんの話を全然理解してなかった
バブルジェットプリンタの開発
http://www.nagare.or.jp/publication/nagare/archive/2005/6.html
http://www.nagare.or.jp/download/noauth.html?dd=assets/files/download/noauth/nagare/24-6/24-6-s03.pdf
『液体の沸点とは，液体の飽和蒸気圧が大気の圧力に等しくなる温度であり，この温度を境として，低温では液体状態の方が熱力学的に安定であり，高温では気体（蒸気）状態の方が安定である．従って，液体の温度が沸点を越えると，ただちに液体から気体への相転移すなわち沸騰がおこるように思われるかもしれない．しかしながら，液体中に蒸気相を作るためには，液体と蒸気の界面を形成するための余分のエネルギーを必要とするため，液体から蒸気への相変化の途中でポテンシャル障壁Wが現れる』
『液体から蒸気への相転移は，液体分子の熱運動によって密度に揺らぎが生じ，ある確率で臨界気泡の大きさ程度の蒸気領域が形成されることによって起こると考えることができる．これは，一次相転移における核生成現象の一例である．古典的な核生成理論によると，相転移の起こりやすさ(単位時間・単位体積当たりの核生成確率)は，N(σ/m)^{1/2} exp(W/kT) の程度と見積もられている．ここに，N は分子密度，m は分子質量，k はボルツマン定数である．』
#σは液体中にできた蒸気泡の表面張力で、Wはポテンシャル障壁
『一方，日常的な現象の中では，例えば，水を鍋で加熱すると，100℃で，簡単に沸騰してしまう．これは，空気などが鍋底のくぼみなどに潜んでいて，初めから気相が存在し，沸騰の核となるためであると考えられている．清浄なビーカーを用いて水を沸騰させると，沸点では沸騰せず，かなり高温になってから突沸することがあるのは，沸騰の核となるようなものがないからである』
まぁ、MDは、通常不純物とかないし、ある意味理想的な状態であるので、むしろ、こうなるのが普通であるらしい。一応、過冷却の限界は液滴サイズに依存するらしいけど、増加の仕方は相当緩やかで、上記のように巨視的な系でも、過熱・過冷却はある程だし、なかなか普通にMDやっても、(熱力学的な)融点・沸点は再現しないと。実験的には、相転移が一次か二次かの判定には、ヒステリシスの有無を確認するのが重要らしいので、前回の計算では、それが見えていると言えなくもない？
#どうでもいいけど、本文中にある"数pl の液滴を一秒あたり数万回噴射"って凄い。1pL=10^(-12)L=10^(-15) m^3なので、直径10^(-5)m=10μmオーダーと考えると、大体酵母のサイズくらい
Structural transformation in supercooled water controls the crystallization rate of ice
http://www.nature.com/nature/journal/v479/n7374/abs/nature10586.html
これは過冷却の話。32768個(32^3)の水分子を、色々考慮に入れて、MD(?)計算したところ、-48度まで凍らずに冷却できたという話。水の過冷却限界は、どこまで下がるのかみたいなことを調べたかったらしい
核生成と界面
http://ist.ksc.kwansei.ac.jp/~nishitani/Lectures/Kyoto/Solidification/nucleation.pdf
Gibbsに始まって、1930年代にBeckerさんという人が古典的核生成理論というのを集大成として、まとめたらしい。周期境界条件を課したMD計算だと、(表面とかないので)均質核生成が起きるという理解でいいんだろうかと思ったけど、古典的核生成理論は界面を形成するエネルギーがうんたら言ってることからも、原子100個程度の系に適応していいのか怪しい。まぁ、そもそも100原子程度で凝固とか沸騰とか出るってのもアレではあるが、Wikipediaに
古典理論の問題点
http://ja.wikipedia.org/wiki/%E6%A0%B8%E7%94%9F%E6%88%90#.E5.8F.A4.E5.85.B8.E7.90.86.E8.AB.96.E3.81.AE.E5.95.8F.E9.A1.8C.E7.82.B9
という項があって、『CNT(古典的核生成理論)は分子の巨視的性質を微視的な動きに適用できることを前提としているが、これは10分子程度からなる小さなクラスタの密度・表面張力・飽和蒸気圧などを扱う際に破綻する。また、核周辺での粒子の相互作用も考慮されていない。』と書いてあった。
過飽和気体における核生成過程の大規模分子動力学計算
- はじめて室内実験との直接比較を実現-
http://www.hokudai.ac.jp/news/130905_pr_lowtem.pdf
北大、高精度な過飽和気体の凝縮核生成モデルの構築手法を開発
http://news.mynavi.jp/news/2013/09/06/145/index.html
Large Scale Molecular Dynamics Simulations of Homogeneous Nucleation
http://arxiv.org/abs/1308.0972
LAMMPS を用いた気相からの凝縮核生成過程の大規模分子動力学計算
http://www.lowtem.hokudai.ac.jp/tech/report_2012/05tanaka.pdf
『核形成過程は物質が相変化する際に起きる物理現象であり、気象学(雲粒形成)や工学(ナノ粒子形成)などといった幅広い分野で研究対象となっており、その古典理論は核形成過程の巨視的な記述を与え広く用いられているが、同理論が与える核生成率は室内実験で得られる実験値から数桁以上もずれることが指摘されており、定量的に信頼できる理論モデルは未だ存在していなかった』
『同理論の問題点として、核生成の際に形成される臨界核の表面エネルギーの評価に巨視的な表面張力を用いるなど、分子レベルの効果を無視している点が挙げられており、その効果を明らかにするために、これまで分子の微視的な運動を計算することにより核生成現象を解析する分子動力学法が用いられてきたが、通常の計算規模は千から数十万分子程度であり、計算できる温度および圧力の範囲は狭い領域に限られるという課題があった』
『気体からの凝縮過程について，超並列計算機を用いて80 億分子までの大規模計算を行いました。その結果，室内実験条件と同様の低過飽和状態においてアルゴンガスの核生成を定量的に再現することに成功しました。核生成過程を決めるナノサイズの臨界核の形成エネルギーを分子動力学計算によって実測した』
80億分子とか、個人ではとても無理だが、まぁ、やっぱり古典的核生成理論は、定量的には色々ダメっぽいと。これ、見ると、上の水の過冷却を見たとかいうNatureの論文しょぼくね？という気がするのだけど。
あと、MD計算にLAMMPSというのを使ってるけど、
[lammps-users] LAMMPS vs GROMACS
http://lammps.sandia.gov/threads/msg05365.html
によれば、数コアしかないような計算機環境では、Gromacsの方が速いが、並列数を増加していった時には、LAMMPSの方がスケールするはず的なことが書いてある。まぁ、普通の計算機環境しかない人は、Gromacs使っとけば、とりあえず最速ではある？
このあたりの話題は、全然勉強したことなかったけど、工業的には、結構重要な話があるっぽいし、なかなか面白い
迷路ソルバー
ふと、A*-アルゴリズムって実装したことなかったな、と思ったので実装してみただけ。A*-アルゴリズムが正しいことは自明でないと思うのだけど、意外とどこにも証明とかなくて、結局、原論文を読めばいいという結論に至った。
A Formal Basis for the Heuristic Determination of Minimum Cost Paths"
http://fai.cs.uni-saarland.de/teaching/winter12-13/heuristic-search-material/Astar.pdf
#!/usr/bin/env python
#--breadth first
def route_bfs(M,start,end):
def next((x,y)):
if x>0 and M[y][x-1]!='*':yield (x-1,y)
if x+1<len(M[y]) and M[y][x+1]!='*':yield (x+1,y)
if y>0 and M[y-1][x]!='*':yield (x,y-1)
if y+1<len(M) and M[y+1][x]!='*':yield (x,y+1)
Q = [(start , [])]
while Q:
cp,rr = Q[0]
del Q[0]
if cp==end:return rr+[end]
for p in next(cp):
if p in rr:continue
Q.append( (p,rr+[cp]) )
return []
#-- depth first
def route_dfs(M,start,end):
def next((x,y)):
if x>0 and M[y][x-1]!='*':yield (x-1,y)
if x+1<len(M[y]) and M[y][x+1]!='*':yield (x+1,y)
if y>0 and M[y-1][x]!='*':yield (x,y-1)
if y+1<len(M) and M[y+1][x]!='*':yield (x,y+1)
def aux(M,s,e,rr):
if s==e:return [s]
routes = [aux(M , p , e , rr+[s]) for p in next(s) if not p in rr]
routes = [x for x in routes if len(x)>0]
if len(routes)>0:return [s]+min(routes,key=lambda x:len(x))
else:return []
return aux(M,start,end,[])
#-- Dijkstra
def route_dijk(M,start,end):
def next((x,y)):
if x>0 and M[y][x-1]!='*':yield (x-1,y)
if x+1<len(M[y]) and M[y][x+1]!='*':yield (x+1,y)
if y>0 and M[y-1][x]!='*':yield (x,y-1)
if y+1<len(M) and M[y+1][x]!='*':yield (x,y+1)
Q = [start]
paths = {start:[start]}
while Q:
cp = Q.pop()
cr = paths[cp]
for np in next(cp):
nr = paths.get(np , None)
if nr is None or len(nr)>=len(cr)+1:
paths[np] = cr+[np]
Q.append(np)
return paths.get(end , [])
#-- A*
def route_astar(M,start,end):
def h((x,y)): #--heuristic
return abs(x-end[0])+abs(y-end[1])
def next((x,y)):
if x>0 and M[y][x-1]!='*':yield (x-1,y)
if x+1<len(M[y]) and M[y][x+1]!='*':yield (x+1,y)
if y>0 and M[y-1][x]!='*':yield (x,y-1)
if y+1<len(M) and M[y+1][x]!='*':yield (x,y+1)
openSet = {start:h(start)}
closeSet = {}
parent = {}
while True:
if len(openSet)==0:return []
node,hval = min(openSet.items() , key=lambda x:x[1])
if node==end:
break
else:
closeSet[node] = openSet[node]
del openSet[node]
for np in next(node):
rval = (hval-h(node)) + (h(np) + 1)
if rval < openSet.get(np,rval):
openSet[np] = rval
parent[np] = node
elif rval < closeSet.get(np,rval):
openSet[np] = rval
parent[np] = node
del closeSet[np]
elif not closeSet.has_key(np) and not openSet.has_key(np):
openSet[np] = rval
parent[np] = node
else:
pass
paths = []
current = end
while True:
paths.append(current)
if current==start:return list(reversed(paths))
current = parent[current]
import time
if __name__=="__main__":
M = """*S*****************************
* * * *
*** ******* ***** * *** *** * *
* * * * * * * *
***** *** * *** * *** *** *****
* * * * * * *
*** * *** * * *** * * ******* *
* * * * * * * * *
*** ***** *** * * * *** * *****
* * * * * * * *
* *********** ********* ***** *
* * * * * * * * *
* *** * * *** * *** ***** *** *
* * * * * * *
* ********* ********* * ***** *
* * * * * * *
* * ******* *** ******* *** ***
* * * * * * * * * * *
*** * * *** * *** * * *** * * *
* * * * * *
*****************************G*"""
M=M.split('\n')
algorithms = [("bfs",route_bfs),("dfs",route_dfs),("Dijkstra",route_dijk),("A*",route_astar)]
for name,route in algorithms:
st = time.clock()
route(M,(1,0),(29,20))
et = time.clock()
print ("%s: %1.4f(sec)" % (name,et-st))
print "------------------"
M = [" "*100]
for name,route in algorithms:
st = time.clock()
route(M,(1,0),(80,0))
et = time.clock()
print ("%s: %1.4f(sec)" % (name,et-st))
print "-------------------"
M = [" "*6 for _ in xrange(6)]
for name,route in algorithms:
st = time.clock()
route(M,(0,0),(4,4))
et = time.clock()
print ("%s: %1.4f(sec)" % (name,et-st))
