Difference between revisions of "Tutorials:Code-01 Python/ja"

From ALPS
Jump to: navigation, search
(Created page with "{{Languages|Tutorials:Code-01_Python}} '''Pythonを利用したALPS入門''' <!--In this tutorial we will show how a simulation can be written in a few lines of code using pyth…")
 
Line 3: Line 3:
 
'''Pythonを利用したALPS入門'''
 
'''Pythonを利用したALPS入門'''
  
<!--In this tutorial we will show how a simulation can be written in a few lines of code using python-ALPS. We will look at the "hello world"-example in the world of physics simulations and perform a Monte Carlo simulation of the classical 2D Ising model with local updates.-->
 
 
このチュートリアルでは、どのように記述しシミュレーションをおこなうかPython-ALPSの簡単な例を紹介します。"hello world"の例として局所更新古典2次元イジングモデルのモンテカルロシミュレーション法をとりあげます。
 
このチュートリアルでは、どのように記述しシミュレーションをおこなうかPython-ALPSの簡単な例を紹介します。"hello world"の例として局所更新古典2次元イジングモデルのモンテカルロシミュレーション法をとりあげます。
<!--
+
 
A skeleton code outlining the typical structure of a Monte Carlo simulation is provided in the file [http://alps.comp-phys.org/static/tutorials2.0.0/code-01-python/ising-skeleton.py ising-skeleton.py] and will be discussed step by step below:-->
 
 
モンテカルロシミュレーションの概要が記述してあるコードは[http://alps.comp-phys.org/static/tutorials2.0.0/code-01-python/ising-skeleton.py ising-skeleton.py]を参照してください。
 
モンテカルロシミュレーションの概要が記述してあるコードは[http://alps.comp-phys.org/static/tutorials2.0.0/code-01-python/ising-skeleton.py ising-skeleton.py]を参照してください。
  
<!--First we import the required python modules-->
 
 
最初に、必要なPythonモジュールをインポートします。
 
最初に、必要なPythonモジュールをインポートします。
  
Line 17: Line 14:
 
  import pyalps.pytools as alpstools
 
  import pyalps.pytools as alpstools
  
<!--We will start by implementing a Simulation class, which will contain the methods for initializing a simulation, running it and storing the measurements into a HDF5 file. -->
 
 
まず、シミュレーションクラスを実装します。次に、シミュレーションの初期化、実行をおこない、計算結果をHDF5ファイルに格納します。
 
まず、シミュレーションクラスを実装します。次に、シミュレーションの初期化、実行をおこない、計算結果をHDF5ファイルに格納します。
  
Line 41: Line 37:
 
         self.abs_magnetization = alpsalea.RealObservable('|m|')
 
         self.abs_magnetization = alpsalea.RealObservable('|m|')
  
<!--The <tt>__init__</tt> method defines how an object of this Simulation class will be instantiated. As arguments the lattice size <math>L</math> and the inverse temperature <math>\beta</math> will be passed. Based on these parameters we deduce the possible Boltzmann weights and initialize a square lattice of Ising spins in a random configuration. Furthermore we also initialize the observables we are going to measure. Here we make use of the python-ALPS framework in order to let the ALPS alea library handle the evaluation of observables for us. Within the class we have also initialized and seeded a random number generator. As an engine we are using the Mersenne Twister MT19937, whose long period and statistical properties make it a good choice for Monte Carlo simulations.-->
+
The <tt>__init__</tt>メソッドはSimulationクラスのオブジェクトがどのようにインスタンス化するのか定義します。引数として格子サイズの<math>L</math>と逆温度<math>beta</math>が渡されます。これらのパラメータに基づいてボルツマンの重みを推定し、ランダム配置でのイジングスピンの正方格子を初期化します。さらに求める観測値の初期化をおこないます。ここではALPS aleaライブラリが観測値の解析をおこなうようにpython-ALPSフレームワークを使用しています。クラス内では乱数種発生の初期化もおこなわれます。乱数の発生にはメルセンヌ・ツイスタMT19937を使用しています。
The <tt>__init__</tt>メソッドはSimulationクラスのオブジェクトがどのようにインスタンス化するのか定義します。引数として格子サイズの<math>L</math>と逆温度<math>beta</math>が渡されます。これらのパラメータに基づいてボルツマンの重みを推定し、ランダム配置でのイジングスピンの正方格子を初期化します。さらに求める観測値の初期化をおこないます。ここではALPS aleaライブラリが観測値の解析をおこなうようにpython-ALPSフレームワークを使用しています。クラス内では乱数種発生の初期化もおこなわれます。乱数の発生にはメルセンヌ・ツイスタMT19937を使用している。
 
  
 
  def save(self, filename):
 
  def save(self, filename):
Line 50: Line 45:
 
         self.magnetization.save(filename)
 
         self.magnetization.save(filename)
  
<!--The <tt>save</tt> method stores the simulation parameters and results in a HDF5-file.-->
 
 
<tt>save</tt>メソッドはシミュレーションパラメータと結果をHDF5ファイルにストアします。
 
<tt>save</tt>メソッドはシミュレーションパラメータと結果をHDF5ファイルにストアします。
  
Line 70: Line 64:
 
         print 'm:\t', self.magnetization.mean, '+-', self.magnetization.error, ',\t tau =', self.magnetization.tau
 
         print 'm:\t', self.magnetization.mean, '+-', self.magnetization.error, ',\t tau =', self.magnetization.tau
  
<!--The <tt>run</tt> method manages the Monte Carlo updates defined in the <tt>step</tt> routine and the measurement of the observables in the <tt>measure</tt> function. While the system is thermalizing we do not perform any measurements. Once all steps are done we print mean, error and autocorrelation time of the observables.-->
+
<tt>run</tt>メソッドはモンテカルロ計算の更新を<tt>step</tt>で定義し、<tt>measure</tt>関数で観測値について設定します。システムが熱化の計算をしている間は他の計測は実行しないでください。すべての計算が完了したら、平均値、誤差、相関時間の出力をおこないます。
<tt>run</tt>メソッドはモンテカルロ計算の更新を<tt>step<tt>で定義し、<tt>measure</tt>関数で観測値について設定します。システムが熱化の計算をしている間は他の計測は実行しないでください。すべての計算が完了したら、平均値、誤差、相関時間の出力をおこないます。
 
  
 
  def step(self):
 
  def step(self):
Line 82: Line 75:
 
       ...
 
       ...
  
<!--The Monte Carlo sweeps are done in the <tt>step</tt> method. In the Metropolis algorithm  a spin is a randomly picked and flipped with probability <math>p_{accept} = min(1,e^{-\beta \Delta E})</math>, <math>\Delta E</math> being the energy difference of the initial and proposed configuration. This procedure is repeated <math>L^2</math>times. The implementation of the Metropolis algorithm is left to you as an exercise. You can make use of the <tt>randint</tt> function defined below:-->
 
 
モンテカルロの掃引は<tt>step</tt>メソッドでおこなわれます。メトロポリスアルゴリズムでは、スピンはランダムに選択され、確率<math>p_{accept} = min(1,e^{-\beta \Delta E})</math>、で放出します。<math>\Delta E</math>は初期構造と提案された構造でのエネルギー差です。この手順は<math>L^2</math>回繰り返されます。メトロポリスアルゴリムの実装は良い練習になるでしょう。次のように記述し<tt>randint</tt>関数の利用が可能です。
 
モンテカルロの掃引は<tt>step</tt>メソッドでおこなわれます。メトロポリスアルゴリズムでは、スピンはランダムに選択され、確率<math>p_{accept} = min(1,e^{-\beta \Delta E})</math>、で放出します。<math>\Delta E</math>は初期構造と提案された構造でのエネルギー差です。この手順は<math>L^2</math>回繰り返されます。メトロポリスアルゴリムの実装は良い練習になるでしょう。次のように記述し<tt>randint</tt>関数の利用が可能です。
  
Line 88: Line 80:
 
         return int(max*self.rng())
 
         return int(max*self.rng())
  
<!--The measurements of our chosen observables are going to be implemented in the method <tt>measure</tt>:-->
 
 
選択した測定値は<tt>measure</tt>メソッドで実施されます。
 
選択した測定値は<tt>measure</tt>メソッドで実施されます。
  
Line 103: Line 94:
 
     self.abs_magnetization << abs(M)/(self.L*self.L)
 
     self.abs_magnetization << abs(M)/(self.L*self.L)
  
<!--The values of the energy and magnetization are determined for the given spin configuration and added to the ALPS observable. The implementation is again left to you as an exercise.-->
 
 
エネルギーと磁化の計算値はスピン配置によって決定され、ALPSobservableに加えられます。試してみてください。
 
エネルギーと磁化の計算値はスピン配置によって決定され、ALPSobservableに加えられます。試してみてください。
  
<!--Once you have completed the implementation of the observable measurements and Metropolis update you can run the simulation using the <tt>vispython</tt> or <tt>alpspython</tt> python interpreter. In this example we will do a scan over different values of <math>\beta = 1/k_B T</math>. The "main" program is given below:-->
+
次に、<tt>vispython</tt>や<tt>alpspython</tt>を使ったシミュレーションの実行方法をやってみてください。
次に、<tt>vispython</tt>や <tt>alpspython</tt>を使ったシミュレーションの実行方法をやってみてください。
 
 
次の例では、<math>\beta = 1/k_B T</math>の異なる値でのスキャンをおこないます。
 
次の例では、<math>\beta = 1/k_B T</math>の異なる値でのスキャンをおこないます。
  
Line 122: Line 111:
 
     sim.save('ising.'+str(beta)+'.h5')
 
     sim.save('ising.'+str(beta)+'.h5')
  
<!--A nice thing about python-ALPS is that if you evaluate composite observables, e.g. of the form <math>U = <A>/<B></math>, a jackknife analysis will be performed and you will automatically obtain correctly evaluated values for the mean and error.
 
As an example we are going to extend our Ising simulation and add measurements of <math>m^2</math> and <math>m^4</math>, from which we determine the Binder cumulant <math>U_4=< m^4> /< m^2>^2</math>. Since the Binder cumulant is usually used to determine the critical temperature using finite size scaling, we are also going to simulate <math>L=6,8</math> in addition.-->
 
python-ALPSを使用する利点は、もし、<math>U = <A>/<B></math>のような観測値は評価する場合、ジャックナイフ法が実行されると、自動的に正しい平均値や誤差が含まれることにあります。
 
例では、イジングモデルを拡張し、Binder cumulant<math>U_4=< m^4> /< m^2>^2</math>から定義した<math>m^2</math>と<math>m^4</math>を観測値として加えます。通常、Binder cumulantは有限サイズスケーリングを用いて臨界温度を求めるのに使用します。また、<math>L=6,8</math>の場合でもシミュレーションしてみます。
 
  
<!--Assuming that you have succesfully implemented these two additional observables and run the simulation we will first load the data and store for each system size <math>L</math> the two observables <math>m^2</math> and <math>m^4</math> as a function of inverse temperature  <math>\beta </math>:-->
+
python-ALPSを使用する利点は、もし、<math>U = <A>/<B></math>のような観測値を評価する場合、ジャックナイフ法が実行されると、自動的に正しい平均値や誤差が含まれることにあります。
 +
例では、イジングモデルを拡張し、バインダーキュムラント<math>U_4=< m^4> /< m^2>^2</math>から定義した<math>m^2</math><math>m^4</math>を観測値として加えます。通常、バインダーキュムラントは有限サイズスケーリングを用いて臨界温度を求めるのに使用します。また、<math>L=6,8</math>の場合でもシミュレーションしてみます。
 +
 
 
これらの2つの系での計算が正常に実行されていると仮定して、まず最初にデータをロードし、各システムサイズ<math>L</math>に格納します。データは逆温度<math>\beta </math>の関数として <math>m^2</math>と<math>m^4</math>の観測値です。
 
これらの2つの系での計算が正常に実行されていると仮定して、まず最初にデータをロードし、各システムサイズ<math>L</math>に格納します。データは逆温度<math>\beta </math>の関数として <math>m^2</math>と<math>m^4</math>の観測値です。
  
Line 134: Line 121:
 
  m4=pyalps.collectXY(data,x='BETA',y='m^4',foreach=['L'])
 
  m4=pyalps.collectXY(data,x='BETA',y='m^4',foreach=['L'])
  
<!--Now we can calculate the Binder cumulant, which we will store in a list of datasets:-->
+
バインダーキュムラントの計算が実行され、datasetのリストに格納されます。
Binder cumulantの計算が実行され、datasetのリストに格納されます。
 
  
 
  u=[]
 
  u=[]
Line 146: Line 132:
 
     u.append(d)
 
     u.append(d)
  
<!--You can plot the Binder cumulant using the commands:-->
 
 
次のコマンドで、Binder Cumulantをプロットできます。
 
次のコマンドで、Binder Cumulantをプロットできます。
  
Line 158: Line 143:
 
  plt.show()
 
  plt.show()
  
<!--What do you observe? If you want to learn more about phase transitions and finite size scaling methods, check out the tutorial
 
[[ALPS_2_Tutorials:MC-07_Phase_Transition |  MC-07 Phase transition in the Ising model ]] .-->
 
  
 
何が観測されましたか?もしさらに相転移や有限サイズスケーリングについて学びたいのであれば、[[ALPS_2_Tutorials:MC-07_Phase_Transition/ja |  MC-07 Phase transition in the Ising model ]]を参照してください。
 
何が観測されましたか?もしさらに相転移や有限サイズスケーリングについて学びたいのであれば、[[ALPS_2_Tutorials:MC-07_Phase_Transition/ja |  MC-07 Phase transition in the Ising model ]]を参照してください。

Revision as of 06:47, 16 March 2012

Pythonを利用したALPS入門

このチュートリアルでは、どのように記述しシミュレーションをおこなうかPython-ALPSの簡単な例を紹介します。"hello world"の例として局所更新古典2次元イジングモデルのモンテカルロシミュレーション法をとりあげます。

モンテカルロシミュレーションの概要が記述してあるコードはising-skeleton.pyを参照してください。

最初に、必要なPythonモジュールをインポートします。

import math
import pyalps
import pyalps.alea as alpsalea
import pyalps.pytools as alpstools

まず、シミュレーションクラスを実装します。次に、シミュレーションの初期化、実行をおこない、計算結果をHDF5ファイルに格納します。

class Simulation:
# Seed random number generator: self.rng() will give a random float from the interval [0,1)
   rng = alpstools.rng(42)
   
   def __init__(self,beta,L):
       self.L = L
       self.beta = beta
       
       # Init exponential map
       self.exp_table = dict()
       for E in range(-4,5,2): 
         self.exp_table[E] = math.exp(2*beta*E)
       
       # Init random spin configuration
       self.spins = [ [2*self.randint(2)-1 for j in range(L)] for i in range(L) ]
       
       # Init observables
       self.energy = alpsalea.RealObservable('E')
       self.magnetization = alpsalea.RealObservable('m')
       self.abs_magnetization = alpsalea.RealObservable('|m|')

The __init__メソッドはSimulationクラスのオブジェクトがどのようにインスタンス化するのか定義します。引数として格子サイズのLと逆温度betaが渡されます。これらのパラメータに基づいてボルツマンの重みを推定し、ランダム配置でのイジングスピンの正方格子を初期化します。さらに求める観測値の初期化をおこないます。ここではALPS aleaライブラリが観測値の解析をおこなうようにpython-ALPSフレームワークを使用しています。クラス内では乱数種発生の初期化もおこなわれます。乱数の発生にはメルセンヌ・ツイスタMT19937を使用しています。

def save(self, filename):
       pyalps.save_parameters(filename, {'L':self.L, 'BETA':self.beta, 'SWEEPS':self.n, 'THERMALIZATION':self.ntherm})
       self.abs_magnetization.save(filename)
       self.energy.save(filename)
       self.magnetization.save(filename)

saveメソッドはシミュレーションパラメータと結果をHDF5ファイルにストアします。

def run(self,ntherm,n):
       # Thermalize for ntherm steps
       self.n = n
       self.ntherm = ntherm
       while ntherm > 0:
           self.step()
           ntherm = ntherm-1
       # Run n steps
       while n > 0:
           self.step()
           self.measure()
           n = n-1
       # Print observables
       print '|m|:\t', self.abs_magnetization.mean, '+-', self.abs_magnetization.error, ',\t tau =', self.abs_magnetization.tau
       print 'E:\t', self.energy.mean, '+-', self.energy.error, ',\t tau =', self.energy.tau
       print 'm:\t', self.magnetization.mean, '+-', self.magnetization.error, ',\t tau =', self.magnetization.tau

runメソッドはモンテカルロ計算の更新をstepで定義し、measure関数で観測値について設定します。システムが熱化の計算をしている間は他の計測は実行しないでください。すべての計算が完了したら、平均値、誤差、相関時間の出力をおこないます。

def step(self):
    for s in range(self.L*self.L):
      # Pick random site k=(i,j)
       ...
      # Measure local energy e = -s_k * sum_{l nn k} s_l
       ...		
      # Flip s_k with probability exp(2 beta e)
      ...

モンテカルロの掃引はstepメソッドでおこなわれます。メトロポリスアルゴリズムでは、スピンはランダムに選択され、確率p_{accept} = min(1,e^{-\beta \Delta E})、で放出します。\Delta Eは初期構造と提案された構造でのエネルギー差です。この手順はL^2回繰り返されます。メトロポリスアルゴリムの実装は良い練習になるでしょう。次のように記述しrandint関数の利用が可能です。

def randint(self,max):
       return int(max*self.rng())

選択した測定値はmeasureメソッドで実施されます。

def measure(self):
   E = 0.	# energy
   M = 0.	# magnetization
   for i in range(self.L):
       for j in range(self.L):
           E -= ...
           M += ...
   # Add sample to observables
   self.energy << E/(self.L*self.L)
   self.magnetization << M/(self.L*self.L)
   self.abs_magnetization << abs(M)/(self.L*self.L)

エネルギーと磁化の計算値はスピン配置によって決定され、ALPSobservableに加えられます。試してみてください。

次に、vispythonalpspythonを使ったシミュレーションの実行方法をやってみてください。 次の例では、\beta = 1/k_B Tの異なる値でのスキャンをおこないます。


L = 4	# Linear lattice size
N = 5000	# of simulation steps
print '# L:', L, 'N:', N
# Scan beta range [0,1] in steps of 0.1
for beta in [0.,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.]:
   print '-----------'
   print 'beta =', beta
   sim = Simulation(beta,L)
   sim.run(N/2,N)
   sim.save('ising.'+str(beta)+'.h5')


python-ALPSを使用する利点は、もし、U = <A>/<B>のような観測値を評価する場合、ジャックナイフ法が実行されると、自動的に正しい平均値や誤差が含まれることにあります。 例では、イジングモデルを拡張し、バインダーキュムラントU_4=< m^4> /< m^2>^2から定義したm^2m^4を観測値として加えます。通常、バインダーキュムラントは有限サイズスケーリングを用いて臨界温度を求めるのに使用します。また、L=6,8の場合でもシミュレーションしてみます。

これらの2つの系での計算が正常に実行されていると仮定して、まず最初にデータをロードし、各システムサイズLに格納します。データは逆温度\beta の関数として m^2m^4の観測値です。

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='ising.L*'),['m^2', 'm^4'])
m2=pyalps.collectXY(data,x='BETA',y='m^2',foreach=['L'])
m4=pyalps.collectXY(data,x='BETA',y='m^4',foreach=['L'])

バインダーキュムラントの計算が実行され、datasetのリストに格納されます。

u=[]
for i in range(len(m2)):
    d = pyalps.DataSet()
    d.propsylabel='U4'
    d.props = m2[i].props
    d.x= m2[i].x
    d.y = m4[i].y/m2[i].y/m2[i].y
    u.append(d)

次のコマンドで、Binder Cumulantをプロットできます。

import pyalps.plot as plt 
plt.figure()
plt.plot(u)
plt.xlabel('Inverse Temperature $\beta$')
plt.ylabel('Binder Cumulant U4 $g$')
plt.title('2D Ising model')
plt.show()


何が観測されましたか?もしさらに相転移や有限サイズスケーリングについて学びたいのであれば、 MC-07 Phase transition in the Ising model を参照してください。