Mathematicaにおけるコンパイラの効用
Mathematicaにおけるコンパイラの効用
松田裕幸
計算環境 Mathematica 14.0、 MacOS Sonoma 14.3.1, CPU Apple M1, メモリ 16GB
Mathematica関数 Compile, AbsoluteTiming, Clear
計算環境 Mathematica 14.0、 MacOS Sonoma 14.3.1, CPU Apple M1, メモリ 16GB
Mathematica関数 Compile, AbsoluteTiming, Clear
はじめに
はじめに
Mathematica にコンパイル機能があることは昔からよく知られていますが、いま一つ使い方がわからず、あまり利用されていません。本稿では具体的な例 https://reference.wolfram.com/language/ref/Compile.html.ja 内[おもしろい例題]を参考に、コンパイルすると高速化(この例の場合は約30倍に高速化)されることを示します。
基本事項の確認
基本事項の確認
コンパイラを知らない人はいないと思いますが、コンパイルの基本がメモリ参照をCPU内メモリ(普通はレジスタ)に切り替えると理解している方は少ないかもしれません。メモリアクセス速度を低い順に並べると ファイル << メインメモリ < CPUレジスタとなります。Mathematica は通常、インタプリタとして機能し、変数が登場するたびにメインメモリへのアクセスが発生します。これに対し、コンパイラは変数をCPUレジスタに割り付けるために、参照速度が上がります。なおコンパイルにおいて当該データが何バイト使用するか事前に知る必要があり、Mathematica のコンパイル関数 Compile では引数に必ず型を加えます。たとえば、} 「引数 の型は実数」はその例です。したがってコンパイルの対象変数は実行時、整数や実数等の数値でなければなりません。ただし、それらの変数が Sin 等の組み込み関数の引数あるいはリスト構造の一部であっても、コンパイルの効果はあります。
{x,_Real
x
コンパイル関数 Compile の使い方
コンパイル関数 Compile の使い方
式 を本体に持つ関数 cf を定義します。
Sin[x]+x^2-1/(1+x)
In[]:=
Clear[cf]
In[]:=
cf[x_]:=Sin[x]+x^2-1/(1+x);cf[N[Pi]]
Out[]=
9.62815
この関数の本体部分をコンパイルするには引数部分 (二重括弧は複数の引数を想定)を含め、本体を Compile で囲みます。結果の CompileFunciton が代入された変数 cf は通常の関数名として使用できます: 。なお、引数が数値を仮定しているため、シンボル を によって強制的に数値化する必要はありません。
{{x,_Real}}
\bcf[Pi]
Pi
N
In[]:=
cf=Compile[{{x,_Real}},Sin[x]+x^2-1/(1+x)]cf[Pi]
Out[]=
CompiledFunction
Out[]=
9.62815
ここでコンパイルの効果をみるために、かなり不自然な関数(同じ計算を100万回
繰り返す)を定義します。現在、Mathematica の Do 関数はかなり最適化されていますが、それでもコンパイルするのとしないとでは約10倍の速度差が出ることがわかります。
繰り返す)を定義します。現在、Mathematica の Do 関数はかなり最適化されていますが、それでもコンパイルするのとしないとでは約10倍の速度差が出ることがわかります。
In[]:=
Clear[cf]
In[]:=
cf[x_Real]:=Do[Sin[x]+x^2-1/(1+x),{1000000}]AbsoluteTiming[cf[N[Pi]];]
Out[]=
{1.1818,Null}
In[]:=
cf=Compile[{{x,_Real}},Do[Sin[x]+x^2-1/(1+x),{1000000}]];AbsoluteTiming[cf[N[Pi]];]
Out[]=
{0.10731,Null}
実例
実例
以下の各小節に「原コード」と「コンパイル化コード」を挙げますが、その前に関数引数部分に関し、コンパイル化前と後の比較を行います。なお With は一時的な値を設定するのに使用し、Block のように内部の変数を局所化しません。つまり With 内で関数は大域的に定義されます。
比較1.
fade[t_]:=t*t*t*(t*(t*6.0-15.0)+10.0);fade=Compile[{{t,_Real}},t*t*t*(t*(t*6.0-15.0)+10.0)];
比較2.
lerp[x_,y_,t_]:=(1.0-t)*x+t*y;lerp=Compile[{{x,_Real},{y,_Real},{t,_Real}},(1.0-t)*x+t*y];
比較3.
With[{grad={{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}}},dot[gradIdx_,x_,y_,z_]:=grad[[gradIdx+1]][[1]]*x+grad[[gradIdx+1]][[2]]*y+grad[[gradIdx+1]][[3]]*z];dot=With[{grad={{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}}},Compile[{{gradIdx,_Integer},{x,_Real},{y,_Real},{z,_Real}},grad[[gradIdx+1]][[1]]*x+grad[[gradIdx+1]][[2]]*y+grad[[gradIdx+1]][[3]]*z]];
比較4.
With[{permutations=Join[(permutations=RandomSample[Range[0,255]]),permutations]},signedNoise[x0_,y0_,z0_]:=(*以降省略*)]signedNoise=With[{permutations=Join[(permutations=RandomSample[Range[0,255]]),permutations]},Compile[{{x0,_Real},{y0,_Real},{z0,_Real}},(*以降省略*)]]
比較5.
With[{octaves=8},classicPerlin[xIndex_,yIndex_,amplitude_,frequency_,gain_,lacunarity_,scale_,increment_,width_,height_]:=(*以降省略*)]classicPerlin=With[{octaves=8},Compile[{{xIndex,_Integer},{yIndex,_Integer},{amplitude,_Real},{frequency,_Real},{gain,_Real},{lacunarity,_Real},{scale,_Real},{increment,_Real},{width,_Integer},{height,_Integer}},(*以降省略*)]]
比較4.比較5.ではコンパイラオプション
CompilationOptions{"InlineExternalDefinitions"True},"CompilationTarget""WVM"
が指定されています。
CompilationTarget はデフォルトでは Cコードを想定します。ただし signedNoise のように計算よりリスト処理が中心の場合、Mathematicaの仮想マシン用コード WVM を指定し、高速化します。
InlineExternalDefinitions を指定した場合、signedNoise 内で使用される外部関数(たとえば dot 等)がインライン展開され、コンパイルされます。これによりコンパイルコードは長くなりますが、関数呼び出しが減る分、処理時間は高速化されます。
原コード
原コード
In[]:=
Clear[dot,fade,lerp,signedNoise,classicPerlin]
In[]:=
With[{grad={{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}}},dot[gradIdx_,x_,y_,z_]:=grad[[gradIdx+1]][[1]]*x+grad[[gradIdx+1]][[2]]*y+grad[[gradIdx+1]][[3]]*z];fade[t_]:=t*t*t*(t*(t*6.0-15.0)+10.0);lerp[x_,y_,t_]:=(1.0-t)*x+t*y;With[{permutations=Join[(permutations=RandomSample[Range[0,255]]),permutations]},signedNoise[x0_,y0_,z0_]:=Module[{x,y,z,ix,iy,iz,g000,g001,g010,g011,g100,g101,g110,g111,n000,n100,n010,n110,n001,n101,n011,n111,u,v,w,nx00,nx01,nx10,nx11,nxy0,nxy1,nxyz},ix=IntegerPart[x0];iy=IntegerPart[y0];iz=IntegerPart[z0];x=x0-ix;y=y0-iy;z=z0-iz;ix=Mod[ix,255]+1;iy=Mod[iy,255]+1;iz=Mod[iz,255]+1;g000=Mod[permutations[[ix+permutations[[iy+permutations[[iz]]]]]],12];g001=Mod[permutations[[ix+permutations[[iy+permutations[[iz+1]]]]]],12];g010=Mod[permutations[[ix+permutations[[iy+1+permutations[[iz]]]]]],12];g011=Mod[permutations[[ix+permutations[[iy+1+permutations[[iz+1]]]]]],12];g100=Mod[permutations[[ix+1+permutations[[iy+permutations[[iz]]]]]],12];g101=Mod[permutations[[ix+1+permutations[[iy+permutations[[iz+1]]]]]],12];g110=Mod[permutations[[ix+1+permutations[[iy+1+permutations[[iz]]]]]],12];g111=Mod[permutations[[ix+1+permutations[[iy+1+permutations[[iz+1]]]]]],12];n000=dot[g000,x,y,z];n100=dot[g100,x-1,y,z];n010=dot[g010,x,y-1,z];n110=dot[g110,x-1,y-1,z];n001=dot[g001,x,y,z-1];n101=dot[g101,x-1,y,z-1];n011=dot[g011,x,y-1,z-1];n111=dot[g111,x-1,y-1,z-1];u=fade[x];v=fade[y];w=fade[z];nx00=lerp[n000,n100,u];nx01=lerp[n001,n101,u];nx10=lerp[n010,n110,u];nx11=lerp[n011,n111,u];nxy0=lerp[nx00,nx10,v];nxy1=lerp[nx01,nx11,v];nxyz=lerp[nxy0,nxy1,w];nxyz]];With[{octaves=8},classicPerlin[xIndex_,yIndex_,amplitude_,frequency_,gain_,lacunarity_,scale_,increment_,width_,height_]:=Module[{noiseVal=0.0,x,y,z,freq=frequency,amp=amplitude},x=xIndex*frequency/scale;y=yIndex*frequency/scale;z=1.0*frequency/scale;Do[noiseVal+=signedNoise[x*freq,y*freq,z*freq]*amp;freq*=lacunarity;amp*=gain,{octaves}];Min[Max[noiseVal,0.0],1.0]]];createImage[img_,r_,g_,b_]:=Image[{r*img,g*img,b*img},InterleavingFalse]perlin[width_Integer,height_Integer]:=Table[classicPerlin[ii,jj,amplitude,frequency,gain,lacunarity,scale,increment,width,height],{ii,0,width},{jj,0,height}];
コンパイル化コード
コンパイル化コード