2013-03-17

物理における確率とトポロジーと産業の関わりについて

作者: Claus Beisbart,Stephan Hartmann
出版社/メーカー: Oxford Univ Pr on Demand
発売日: 2011/12/15
メディア: ハードカバー
この商品を含むブログ (1件) を見る

が届いたので読んでいきたいですね．

あとのーてぃさんに教えてもらった

T.Kaczynski K.Mischaikow M.Mrozek,「Computational Homology」,Springer,2004
M.Mrozek B.Batko,「Coreduction Homology Algorithm」,Discrete Comput Geom(2009) 41;96-118
玉木大,「広がりゆくトポロジーの世界-言語としてのホモトピー論」, 現代数学社,2012

も読んでみたいかな？
ただ修士論文が終わるまでは修士論文に集中したいので読み始められるのは夏休みぐらいになるかもしれない．

のーてぃさんのトポロジーと産業の関わりについての面白い話はこちらから↓(pdf注意)
http://twitdoc.com/upload/conoughty/how-to-use-algebraic-topology.pdf

特にバウンダリー作用素を行列として扱うことには目から鱗が落ちました．これならコンピューターで再現できる！

2013-03-12

ビューティフルコード

Fortran

いつか読むはきっと読まない，ということで学校で『ビューティフル・コード』を借りてきました．

ビューティフルコード (THEORY/IN/PRACTICE)

作者: Brian Kernighan,Jon Bentley,まつもとゆきひろ,Andy Oram,Greg Wilson,久野禎子,久野靖
出版社/メーカー: オライリージャパン
発売日: 2008/04/23
メディア: 大型本
購入: 30人クリック: 617回
この商品を含むブログ (190件) を見る

20個ぐらいのキーワード（Python,Ruby,Java,信頼性，etc...）に分けて各専門家のエッセイを載せているというもの．その中で「第15章　美しいデザインの長期にわたる恩恵」が素晴らしすぎたので少しまとめてみる．

・CERNの数学ライブラリは読みやすいコード
→線形代数の部分は現在LAPACKライブラリにまとめられている．
　→この本では例としてLAPACKライブラリのSGBSVルーチンをあげている．
・ルーチンの使い方をコードの一番初めに書いておく．
　→マニュアルの説明とコード中のドキュメントを同じ内容にしておくことは良いこと
・入力に暗黙的な制限を作ると後から大変になる．
・コメントに引数の明示的な説明を書いておく．
　→また各引数が整合性を保っているか簡単なエラー処理をルーチンの中に書いておくのは良いこと．
・簡単なテストとサンプルをコード中に書いておく：重要

一番感動したのが最後の部分で，スモークテストでもいいからきちんと，Aという入力に対してA'という出力が出る，ということをサブルーチン単位で書くことは本当に大切なことだと思う．
このことを理解しない人がむちゃくちゃ多くて最近びっくりすることが多いので．

この章の内容だけでも科学的なコードを書く人が背景として持っておくと良いと思う．

ただ最後にも述べられているように，

美しい解法がいつまでも使えるわけではない

ということはコンピュータに携わる人は常に考えておかなければならない．

次は『リーダブルコード』を読みたい．

2013-03-11

離散フーリエ変換のコード

Fortran Math Physics

フーリエ変換については下記のURLが分かりやすい。
フーリエ級数展開をベクトルで直観的に理解するblog.physips.com

デジタル信号→離散フーリエ変換(DFT)→フーリエ級数(FS)→元の関数の出力(FA)という感じのコードです．
更新していきますが，とりあえずメモ代わり．
同一ディレクトリにDFT.f90,FS.f90,FA.f90って作って，

gfortran DFT.f90 FS.f90 FA.f90
./a.out

とかでいちおう動きます．

DFT.f90

Subroutine DFT(a,b_re,b_im,N,T)                                                                     
! This program calculates the descrete fourier transform.                                           
! This program transform tha array according to                                                     
! Xk = T / N * \sum^N-1_n=0 xn W^nk_N (k = 0,1,2,\cdots,N-1)                                        
! T and N is defined below.                                                                         
                                                                                                    
! Input                                                                                             
! array a   : the sampling data                                                                     
! double  T : sampling time                                                                         
! integer N : sampling frequency                                                                    
! Output                                                                                            
! array b_re,b_im : the transformed data                                                            
                                                                                                    
! Sample data                                                                                       
! Input                                                                                             
! N = 8, T = 1.d0                                                                                   
! a(i) = i / 8.d0                                                                                   
! Output                                                                                            
!     X(i)     b_re     b_im                                                                        
!     X(0)    0.438    0.000                                                                        
!     X(1)   -0.062    0.151                                                                        
!     X(2)   -0.062    0.063                                                                        
!     X(3)   -0.062    0.026                                                                        
!     X(4)   -0.062    0.000                                                                        
!     X(5)   -0.063   -0.026                                                                        
!     X(6)   -0.063   -0.062                                                                        
!     X(7)   -0.063   -0.151           
                                                                                                    
  implicit none                                                                                     
  real(8), dimension(0:N-1) :: a,b_re,b_im                                                          
  real(8) :: pi,T                                                                                   
  integer :: i,j,k,N                                                                                
                                                                                                    
  pi = 3.14159265358979                                                                             
                                                                                                    
  do i = 0,N-1                                                                                      
   b_re(i) = 0.d0                                                                                   
   b_im(i) = 0.d0                                                                                   
  end do                                                                                            
                                                                                                    
  do i = 0,N-1    
    do j = 0,N-1                                                                                    
      b_re(i) = b_re(i) + (dble(T) / dble(N)) * cos(2.d0*pi*j*i/8.d0) * a(j)                        
      b_im(i) = b_im(i) + (dble(T) / dble(N)) * (-1.d0) * sin(2.d0*pi*j*i/8.d0) * a(j)              
    end do                                                                                          
  end do                                                                                            
                                                                                                    
  write(*,'(3A9)') "X(i)","b_re","b_im"                                                             
  do i = 0,N-1                                                                                      
    write(*,'(A7,I1,A1,2F9.3)') "X(",i,")", b_re(i),b_im(i)                                         
  end do                                                                                            
end subroutine DFT

FS.f90

subroutine FS(b_re,b_im,x,N,T)                                                                      
! This program reverse the fourier transform data                                                   
  implicit none                                                                                     
  real(8), dimension(0:N-1) :: b_re,b_im                                                            
  integer :: N, it, i, j, k                                                                         
  real(8) :: pi, T, x, dt, rt                                                                       
                                                                                                    
  pi = 3.141592653589793238                                                                         
  dt = 0.01                                                                                         
  x  = 0.d0                                                                                         
                                                                                                    
  open(10,file='fourier.txt',status='unknown')                                                      
  write(*,'(2A10)') "t","x"                                                                         
                                                                                                    
  N = N / 2                                                                                         
  do it = 1,100                                                                                     
   rt = dt * it                                                                                     
   x = b_re(0)                                                                                      
   do i = 1, N                                                                                      
     x = x + (b_re(i)+b_re(8-N)) * cos(2 * pi * i * rt / T) + (-b_im(i)+b_im(8-i)) * sin(2 * pi * i * rt / T)
   end do                                                                                           
   write(*,'(2F10.2)') rt, x                                                                        
   write(10,'(2F10.2)') rt, x                                                                       
  end do                                                                                            
end subroutine FS

FA.f90

program FA                                                                                          
  implicit none                                                                                     
  real(8), allocatable, dimension(:) :: a, b_re, b_im                                               
  integer :: N,i                                                                                    
  real(8) :: T,x                                                                                    
                                                                                                    
  N = 8                                                                                             
  T = 1.d0                                                                                          
  allocate(a(0:N-1),b_re(0:N-1),b_im(0:N-1))                                                        
                                                                                                    
  do i = 0,7                                                                                        
   a(i) = dble(i) / dble(N)                                                                         
   b_re(i) = 0.d0                                                                                   
   b_im(i) = 0.d0                                                                                   
  end do                                                                                            
                                                                                                    
  call DFT(a,b_re,b_im,N,T)                                                                         
  call FS(b_re,b_im,x,N,T)                                                                          
end program FA

2013-03-07

上野旅行記

ちょっと前，東京に出張に行ったのでその記念に．
上野公園は実は初めてだったのでちょっとうきうき．

アンニュイな熊の像．恋人たちの聖地ってどこにでもありますよね．

一番びっくりしたのはスカイツリーのくっきり具合．こんなにはっきり見えるんだ．

で，本命はここ．国立西洋美術館．東京はしゃれてるなぁ．

この間に珈琲飲みの聖地，「○山珈琲店」に行ってきたんだけど写真をあげるのはあれかな，と思い割愛．
とりあえず，人生観が変わる一杯でした．美味しいとかそういうレベルじゃない．

理系学生らしく国立科学博物館に癒されて終了．一人で回る上野も中々楽しいものですね．

2013-03-01

Fortranのnint

Fortran

Fortranのnint(実数を四捨五入する関数)がどこで四捨五入してるのかなー，と思い調べました．
具体的に言うと，
1. 四捨五入する小数点以下の桁数
2. 整数部分は関係ないのかどうか？

入力

Program round
  implicit none
  integer :: i
  real(8),dimension(8) :: a,b

  a(1) = 0.33
  a(2) = 0.44
  a(3) = 0.49
  a(4) = 0.50
  a(5) = 0.51
  a(6) = 0.55
  a(7) = 0.60
  a(8) = 0.90

  write(*,'(3A10)') "i","a(i)","b(i)"
  do i = 1,8
    b(i) = nint(a(i))
    write(*,'(I10,2F10.2)') i,a(i),b(i)
  end do
end Program round

出力

i a(i) b(i)
1 0.33 0.00
2 0.44 0.00
3 0.49 0.00
4 0.50 1.00
5 0.51 1.00
6 0.55 1.00
7 0.60 1.00
8 0.90 1.00

入力

Program round
  implicit none
  integer :: i
  real(8),dimension(8) :: a,b

  a(1) = 1.33
  a(2) = 1.44
  a(3) = 1.49
  a(4) = 1.50
  a(5) = 1.51
  a(6) = 1.55
  a(7) = 1.60
  a(8) = 1.90

  write(*,'(3A10)') "i","a(i)","b(i)"
  do i = 1,8
    b(i) = nint(a(i))
    write(*,'(I10,2F10.2)') i,a(i),b(i)
  end do
end Program round

出力

i a(i) b(i)
1 1.33 1.00
2 1.44 1.00
3 1.49 1.00
4 1.50 2.00
5 1.51 2.00
6 1.55 2.00
7 1.60 2.00
8 1.90 2.00

というわけで，
1. 小数点以下1桁目を四捨五入している．
2. 整数部分は関係ない．
が答えですかねー．

2013-02-27

モンテカルロ法精度

Fortran Economy(経済、統計） Math Physics

モンテカルロ法において粒子に対する精度の話になったときに「えー，粒子数を $n$ としたとき精度は $\sqrt{n}$ でしょー」という話を聞いたときから，モンテカルロ法の精度が気になってはいたんですよね．いったい精度が $\sqrt{n}$ とはどういう状態を指すのか．
というわけで今回調べてみました．（モンテカルロ法自体については，itさんの素晴らしい解説があるのでそちらに丸投げで　ﾟωﾟ) つ○）

さて，精度が $\sqrt{n}$ というのは，粒子数が $n$ のモンテカルロ法の解を $a_{n}$ ，真の解を $a$ したとき，

$\left| a_{n} - a \right| \sim \sqrt{n}$

が成り立つことかなー，とか思ってたんですけど，これだと左辺はnに対して単調減少なのに右辺は単調増加でおかしなことになるわなー，ということで，こうか，と．

仮説：モンテカルロ法の精度が $\sqrt{n}$ であるとは，次の式が成り立つことである．
$\left| a_{n} - a \right| \sim \frac{1}{\sqrt{n}}$

つまり，Log10(粒子数)に対して誤差( $\left| a_{n} - a \right|$ )をプロットしたときに，

この曲線に乗るようにプロット点がくるのかなぁと．

さて，実際に粒子数を増やして $\pi$ を順次求めるプログラムを書いてみる．

Program MonteCarlo
  implicit none
  real(8) :: pi, seed, seed2, test, ans,error,fe,re,nr
  real(8),dimension(100000000) :: x,y
  integer :: i,j,k,ntest,nmax,sum,fntest

  pi    = 3.14159265358979d0

  call random_number(x)
  call random_number(y)

  write(*,'(A16,A9,A16,A16,A16)'),"squre of n ratio","ntest","answer","error","error ratio"
  do i = 1,8
     ntest = 10 ** i
     sum   = 0
     error = 0.d0
     do j = 1,ntest
       test = x(j) ** 2 + y(j) ** 2
       if (test < 1.d0) sum = sum + 1
     end do
     ans = 4.d0 * dble(sum) / dble(ntest)
     error = ans - pi
     error = abs(error)
     if(i == 1) then
       fe = error
       re = fe / fe
       fntest = ntest
     else
       re = error / fe
     end if
     nr = dble(ntest) / (fntest)
     nr = sqrt(nr)
     write(*,'(f16.4,I9,f16.8,f16.8,f16.8)'),nr,ntest,ans,error,re
  end do

end Program MonteCarlo

でこれに対する出力は次の通り．

うーん，良く分からんなー，ということでLog10(粒子数) - 誤差曲線を取ってみると，

となって，まぁ乗ってるのかなぁ，と思わなくもないけど粒子数が多くなると桁が小さくなって良く分からない．
しょうがないから誤差でも両対数（ただし底は全部10である）を取ってみる．

グラフの意味としては，正の方向に行けば行くほど誤差が大きいってことですー．
つまり，このグラフは粒子数が大きくなればなるほど，誤差が小さくなっていることを表しています．

で，仮説が正しければこれは次のグラフに乗るはず．

で，重ね合わせてみると，

お，けっこう綺麗に乗ってるじゃん，ということで，

仮説：モンテカルロ法において，粒子数を $n$ ，粒子数 $n$ での解を $a_{n}$ ，真の解を $a$ としたとき，精度が $\sqrt{n}$ であるとは，次の式が成り立つことである．
$\left| a_{n} - a \right| \sim \frac{1}{\sqrt{n}}$

がある程度正しいのかなーという結論に至りました．
多分ねー，もう少し拡張できて，

任意の粒子数 $n_{i},n_{j}$ に対して，
$\frac{\left| a_{n{i}} - a \right|}{\left| a_{n_{j}} - a \right|} \sim \frac{\sqrt{n{j}}}{\sqrt{n_{i}}}$

ぐらいが成り立つのではないか，と思ったりもする．

2013-02-25

Macの移行ユーティリティ

リアルでの引越しついでにMacの移行ユーティリティについて何点か

・送り先のOSが送り元のOSよりグレードが低いと移行のとき認識してくれない
・一番エラーが多いのはイラレとかの認証が必要な外部ソフト
・Mac製品に関してはあんまりエラーないような

とりあえず認証系に関しては面倒くさがらず送り元のMacの認証を解除して送り先で認証する．
こんな手間を考えてたら，認証が必要なソフトを入れていない場合は移行ユーティリティを使う．イラレとか入れてるんだったら新しく組みなおす．ぐらいが妥当かも知れませぬ