2016-04-02 8 views
2

丙よりも長い時間をかけてipyparallel使用して並列化、モンテカルロシミュレーション:連載

は私が散乱媒体で光子輸送のモンテカルロシミュレーションを作っています。私は並列化しようとしていますが、実行時の性能向上をシリアルシミュレーションと比較して観察するのが難しいです。

モンテカルロコードは以下にあります。クラスPhotonには、単一光子の輸送と散乱を計算するためのさまざまなメソッドが含まれています。クラスRunPhotonPackageは、散乱媒体の所定の厚さLに対して光子のシリーズNを実行します。それらは一瞬私の唯一の入力パラメータです:

import matplotlib.pyplot as plt 
import numpy as np 
from numpy.random import random as rand 


NPHOTONS = 100000 # Nb photons 
PI = np.pi 
EPS = 1.e-6 
L = 100. # scattering layer thickness 

class Photon(): 

    mut = 0.02 
    k = [0,0,1] 

    def __init__(self,ko,pos): 
     Photon.k = ko 
     self.x = pos[0] 
     self.y = pos[1] 
     self.z = pos[2] 

    def move(self): 
     ksi = rand(1) 
     s = -np.log(1-ksi)/Photon.mut 

     self.x = self.x + s*Photon.k[0] 
     self.y = self.y + s*Photon.k[1] 
     self.z = self.z + s*Photon.k[2]  
     zPos = self.z 
     return zPos 

    def exittop(self): 

     newZpos = 0 


    def exitbase(self): 
     newZpos = 0 


    def HG(self,g): 
     rand_teta = rand(1) 
     costeta = 0.5*(1+g**2-((1-g**2)/(1-g + 2.*g*rand_teta))**2)/g 

     return costeta 

    def scatter(self): 
     # calculate new angle of scattering 
     phi = 2*PI*rand(1)     
     costeta = self.HG(0.85) 
     sinteta = (1-costeta**2)**0.5 


     sinphi = np.sin(phi) 
     cosphi = np.cos(phi) 

     temp = (1-Photon.k[2]**2)**0.5 

     if np.abs(temp) > EPS:   

      mux = sinteta*(Photon.k[0]*Photon.k[2]*cosphi-Photon.k[1]*sinphi)/temp + Photon.k[0]*costeta 
      muy = sinteta*(Photon.k[1]*Photon.k[2]*cosphi+Photon.k[0]*sinphi)/temp + Photon.k[1]*costeta 
      muz = -sinteta*cosphi*temp + Photon.k[2]*costeta 

     else: 
      mux = sinteta*cosphi 
      muy = sinteta*sinphi 
      if Photon.k[2]>=0: 
       muz = costeta 
      else: 
       muz = -costeta 


     # update the new direction of the photon 
     Photon.k[0] = mux 
     Photon.k[1] = muy 
     Photon.k[2] = muz   


class RunPhotonPackage(): 

    def __init__(self,L,NPHOTONS): 
     self.L = L 
     self.NPHOTONS = NPHOTONS 

    def RunPhoton(self): 
     Dist_Pos = np.zeros((3,self.NPHOTONS)) 
     # loop over number of photon 
     for i in range(self.NPHOTONS): 

      # inititate initial photon direction 
      k_init = [0,0,1] 
      k_init_norm = k_init/np.linalg.norm(k_init) # initial photon direction. 
      # initiate new photon with initial direction 
      pos_init = [0,0,0] 
      newPhoton = Photon(k_init_norm,pos_init) 
      newZpos = 0. 

      # while the photon is still in the layer, move it and scatter it 
      while ((newZpos >= 0.) and (newZpos <= self.L)): 

       newZpos = newPhoton.move() 
       newscatter = newPhoton.scatter() 

      Dist_Pos[0,i] = newPhoton.x 
      Dist_Pos[1,i] = newPhoton.y 
      Dist_Pos[2,i] = newPhoton.z 


     return Dist_Pos 

私は、様々な層の厚さの長さに対する位置ヒストグラムと光子の与えられた数を記録するために、次のシリアルコードを実行します。

で、その後
import time 
tic = time.time() 
dictresult = {} 
for L in np.arange(10,100,10): 
    print('L={0} m'.format(L)) 
    Dist_Pos = RunPhotonPackage(L,10000).RunPhoton() 
    dictresult['{0}'.format(L)]=Dist_Pos 
toc = time.time() 
print('sec Elapsed: {0}s'.format(toc-tic)) 

この実行:

sec Elapsed: 26.425330162s 

私はipyparallel使用してコードを並列化しよう:

import ipyparallel 
clients = ipyparallel.Client() 
clients.ids 
dview = clients[:] 

dview.execute('import numpy as np') 
dview.execute('from numpy.random import random as rand') 
dview['PI'] = np.pi 
dview['EPS']= 1.e-6 

dview.push({"Photon": Photon, "RunPhotonPackage": RunPhotonPackage}) 

def RunPhotonPara(L): 
    LayerL = RunPhotonPackage(L,10000) 
    dPos = LayerL.RunPhoton() 
    return dPos 

tic = time.time() 
dictresultpara = [] 
for L in np.arange(10,100,10): 
    print('L={0}'.format(L)) 
    value = dview.apply_async(RunPhotonPara,L) 
    dictresultpara.append(value) 
    clients.wait(dictresultpara) 
toc = time.time() 
print('sec Elapsed: {0}s'.format(toc-tic)) 

それが実行されている:だから以上

sec Elapsed: 55.4289810658s 

ダブルタイム!私はこれをubuntu 32ビットで4つのコアで実行していて、ipcluster start -n 4を使って1つのコントローラと4つのエンジンをlocalhostで起動しています。パラレル化されたコードは、直列のコードを実行するのにかかる時間が1/4になると予想していましたが、明らかにそうではありません。

なぜそれを修正するのですか?

アドバイスありがとうございます。

グレッグ

+0

を。 1マシンで13秒、4コアで0.63秒、意味がありません。(Win10、Anacondaインストール) – roadrunner66

+0

それは私が考えていたものです。私はWindows 7、anacondaのインストール、8つのコアでそれを実行し、シリアルシミュレーションでは18秒、パラレル8コアでは28を得ました。 RunPhotonPackageをclass RunPhotonPackage(object)として定義するとエラーが発生しました。オブジェクトを削除するとコードが実行されます。私は編集し、私の質問を修正しました。 – gregory

+0

Hm、まだ連続12秒、4つのコアに平行0.3-0.5秒、すなわち速すぎる。しかし、私はあなたのコードを実際に読んだり、出力を見たりしませんでした。 – roadrunner66

答えて

0

例を単純化するためにいくつかの変更を加えました。私のMacではシリアルバージョンが〜18秒で実行され、4エンジンのパラレルバージョンは約半分で実行されます。これは、タスクの不均一な持続時間を考えると合理的です。

以前に設定された方法では、エンジンにエラーが発生していたため、高速に復帰しました。辞書を介してクラスを渡すことは適切ではないようです。代わりに、コードは各エンジンのクラスを定義するモジュールをインポートするようになりました。 私は単純にこの例ののsys.pathをハックしましたが、おそらく実稼働環境でこれを適切に処理するでしょう。

私はあなたがループ内に "待ち"したくないと思う。また、async_apply()メソッドよりもasync_map()メソッドが便利だったようです。

これを実行するには、ディレクトリを作成し、そのディレクトリの "photon.py"というファイルに次のコードをコピーし、空の "init .py"も作成します。 sys.pathに挿入するコード内の行を変更して、新しいディレクトリを参照します。変更があるディレクトリと「Pythonのphoton.py」を実行します。これは `あきれるほどparallel`なので、それはかなり簡単であるべき

# photon.py 

import ipyparallel 
import numpy as np 
from numpy.random import random as rand 
import time 

NPHOTONS = 100000 # Nb photons 
PI = np.pi 
EPS = 1.e-6 
L = 100. # scattering layer thickness 

class Photon(): 

    mut = 0.02 
    k = [0,0,1] 

    def __init__(self,ko,pos): 
     Photon.k = ko 
     self.x = pos[0] 
     self.y = pos[1] 
     self.z = pos[2] 

    def move(self): 
     ksi = rand(1) 
     s = -np.log(1-ksi)/Photon.mut 

     self.x = self.x + s*Photon.k[0] 
     self.y = self.y + s*Photon.k[1] 
     self.z = self.z + s*Photon.k[2]  
     zPos = self.z 
     return zPos 

    def exittop(self): 

     newZpos = 0 


    def exitbase(self): 
     newZpos = 0 


    def HG(self,g): 
     rand_teta = rand(1) 
     costeta = 0.5*(1+g**2-((1-g**2)/(1-g + 2.*g*rand_teta))**2)/g 

     return costeta 

    def scatter(self): 
     # calculate new angle of scattering 
     phi = 2*PI*rand(1)     
     costeta = self.HG(0.85) 
     sinteta = (1-costeta**2)**0.5 


     sinphi = np.sin(phi) 
     cosphi = np.cos(phi) 

     temp = (1-Photon.k[2]**2)**0.5 

     if np.abs(temp) > EPS:   

      mux = sinteta*(Photon.k[0]*Photon.k[2]*cosphi-Photon.k[1]*sinphi)/temp + Photon.k[0]*costeta 
      muy = sinteta*(Photon.k[1]*Photon.k[2]*cosphi+Photon.k[0]*sinphi)/temp + Photon.k[1]*costeta 
      muz = -sinteta*cosphi*temp + Photon.k[2]*costeta 

     else: 
      mux = sinteta*cosphi 
      muy = sinteta*sinphi 
      if Photon.k[2]>=0: 
       muz = costeta 
      else: 
       muz = -costeta 


     # update the new direction of the photon 
     Photon.k[0] = mux 
     Photon.k[1] = muy 
     Photon.k[2] = muz   


class RunPhotonPackage(): 

    def __init__(self,L,NPHOTONS): 
     self.L = L 
     self.NPHOTONS = NPHOTONS 

    def RunPhoton(self): 
     Dist_Pos = np.zeros((3,self.NPHOTONS)) 
     # loop over number of photon 
     for i in range(self.NPHOTONS): 

      # inititate initial photon direction 
      k_init = [0,0,1] 
      k_init_norm = k_init/np.linalg.norm(k_init) # initial photon direction. 
      # initiate new photon with initial direction 
      pos_init = [0,0,0] 
      newPhoton = Photon(k_init_norm,pos_init) 
      newZpos = 0. 

      # while the photon is still in the layer, move it and scatter it 
      while ((newZpos >= 0.) and (newZpos <= self.L)): 

       newZpos = newPhoton.move() 
       newscatter = newPhoton.scatter() 

      Dist_Pos[0,i] = newPhoton.x 
      Dist_Pos[1,i] = newPhoton.y 
      Dist_Pos[2,i] = newPhoton.z 

     return Dist_Pos 

def RunPhoton(L): 
    print('L={0}'.format(L)) 
    return RunPhotonPackage(L, 10000).RunPhoton() 

def serialTest(values): 
    print "Running serially..." 
    tic = time.time() 
    results = map(RunPhoton, values) 
    print results 
    toc = time.time() 
    print('sec Elapsed: {0}s'.format(toc-tic)) 

def parallelTest(values): 
    print "Running in parallel..." 
    client = ipyparallel.Client() 
    view = client[:] 

    view.execute('import sys') 

    # CHANGE THIS PATH TO REFER TO WHEREVER YOU PUT THIS CODE 
    view.execute('sys.path.insert(0, "/Users/rjp/ipp")') 
    view.execute('from photon import *') 

    tic = time.time() 
    asyncResults = view.map_async(RunPhoton, values) 
    print asyncResults.get() 
    toc = time.time() 
    print('sec Elapsed: {0}s'.format(toc-tic))  


if __name__ == "__main__": 
    values = np.arange(10, 100, 10) 

    serialTest(values) 
    parallelTest(values) 
+0

hei、それは完璧に動作します!シリアルバージョンよりも2倍高速のパラレルバージョンもあります。ありがとう。 – gregory

+0

Hei Rich、もう少し短い質問ですか?あなたが与えたソリューションは、 "python photon.py"を端末で完璧に動かしていましたが、parallelTest(値)をスパイダで実行しようとすると失敗します(serialTest(values)はうまく動作します)。スパイダーでは、asyncResults.get()コマンドは何も返さず、プログラムはハングアップしているようです。どうして? – gregory

関連する問題