ご指摘のとおり、問題は 2 つあります。一方では、非常に多くの試行を行う必要があるため、記憶の問題があります。一方で、これらすべての試行を処理する必要があるため、パフォーマンスの問題があります。
各問題の解決策は、多くの場合、他の問題に悪影響を及ぼします。私見、最善のアプローチは妥協点を見つけることです。
ベクトル化に必要な膨大な配列を取り除き、別の戦略を使用してループを実行する場合にのみ、より多くの試行を行うことができます。おそらく最適なパフォーマンスを犠牲にして、より多くのトライアルを使用する可能性を優先します。
Matlab プロファイラーでコードをそのまま実行すると、すべての変数の初期メモリ割り当てに時間がかかることがすぐにわかります。plotまた、 コマンドとhold allコマンドがすべての行の中で最も時間のかかる行であることも示しています。さらに試行錯誤を重ねた結果、エラーが発生し始めるtrials前に実行できるの最大値が残念なほど低いことがわかりました。OUT OF MEMORY
Matlab での制限についていくつか知っていれば、ループを大幅に高速化できます。Matlab の古いバージョンでは、「ベクトル化された」コードを優先して、ループを完全に回避する必要があることが真実でした。最近のバージョン (R2008a 以降と思われます) では、Mathworks は、実行中に M コードをその場で機械語に変換する JIT アクセラレータ (ジャストインタイム コンパイラ) と呼ばれるテクノロジを導入しました。簡単に言えば、JIT アクセラレータを使用すると、コードが Matlab のインタープリターをバイパスし、基盤となるハードウェアとより直接的にやり取りできるようになり、時間を大幅に節約できます。
Matlab ではループを避けるべきであるというアドバイスをよく耳にしますが、一般的には当てはまりません。ベクトル化には依然として価値がありますが、ベクトル化されたコードのみを使用して実装されたかなり複雑な手順は、多くの場合、判読しにくく、理解しにくく、変更しにくく、維持するのが困難です。ループを使用する同じ手順の実装には、多くの場合、これらの欠点がなく、さらに、より高速で必要なメモリが少なくて済みます。
残念ながら、JIT アクセラレータには厄介な (そして私見ですが、不必要な) 制限がいくつかあります。
そのようなことの 1 つですplot。一般に、ループにデータの収集と操作以外のことをさせず、プロット コマンドなどをループの後まで遅らせることをお勧めします。
別のそのようなことは次のとおりですhold。hold関数はMatlab の組み込み関数ではありません。つまり、M 言語で実装されています。Matlab の JIT アクセラレータは、ループ内で使用される場合、非組み込み関数を高速化できません。つまり、ループ全体が機械語の速度ではなく、Matlab の解釈速度で実行されます! したがって、ループの後までこのコマンドも遅らせてください:)
ご参考までに、この最後のステップは大きな違いを生む可能性があります。関数本体をコピーして上位レベルのループに貼り付けると、パフォーマンスが 1200 倍向上したケースを私は知っています。実行時間の日数は数分に短縮されました!)。
ループには実際には別の小さな問題があります (これは非常に小さく、かなり不便です。すぐに同意します)。ループ変数の名前はi. 名前iはMatlabの虚数単位の名前であり、名前解決も各反復で不必要に時間を消費します。小さいですが、無視できません。
さて、これらすべてを考慮して、次の実装にたどり着きました。
function [hot, cold, h] = MuonTracks(tracks)
% NOTE: no variables larger than 1x1 are initialized
width = 1e-4;
height = 2e-4;
% constant used for Laplacian noise distribution
bL = 15 / sqrt(2);
% Loop through all tracks
X = [];
hot = 0;
ii = 0;
while ii <= tracks
ii = ii + 1;
% Note that I've inlined (== copy-pasted) the original laprnd()
% function call. This was necessary to work around limitations
% in loops in Matlab, and prevent the nececessity of those HUGE
% variables.
%
% Of course, you can still easily generalize all of this:
% the new data
u = rand-0.5;
Ystart = 15;
Xstart = 800*rand-400;
Xend = Xstart - bL*sign(u)*log(1-2*abs(u));
b = (Ystart*Xend)/(Xend-Xstart);
% the test
if ((b < height && b > 0)) ||...
(Xend < width/2 && Xend > -width/2)
hot = hot+1;
% growing an array is perfectly fine when the chances of it
% happening are so slim
X = [X [Xstart; Xend]]; %#ok
end
end
% This is trivial to do here, and prevents an 'else' in the loop
cold = tracks - hot;
% Now plot the chosen ones
h = figure;
hold all
Y = repmat([15;0], 1, size(X,2));
plot(X, Y, 'r');
end
この実装では、次のことができます。
>> tic, MuonTracks(1e8); toc
Elapsed time is 24.738725 seconds.
完全に無視できるメモリフットプリントで。
プロファイラーは、コードに沿った労力の適切で均一な分布も表示するようになりました。メモリの使用やパフォーマンスのために実際に目立つ行はありません。
これは可能な限り最速の実装ではない可能性があります (明らかな改善が見られた場合は、自由に編集してください)。しかし、あなたが待つ気があるなら、あなたはできるでしょうMuonTracks(1e23)（またはそれ以上:)
また、Matlab MEX ファイルにコンパイルできる C での実装も行いました。
/* DoMuonCounting.c */
#include <math.h>
#include <matrix.h>
#include <mex.h>
#include <time.h>
#include <stdlib.h>
void CountMuons(
unsigned long long tracks,
unsigned long long *hot, unsigned long long *cold, double *Xout);
/* simple little helper functions */
double sign(double x) { return (x>0)-(x<0); }
double rand_double() { return (double)rand()/(double)RAND_MAX; }
/* the gateway function */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int
dims[] = {1,1};
const mxArray
/* Output arguments */
*hot_out = plhs[0] = mxCreateNumericArray(2,dims, mxUINT64_CLASS,0),
*cold_out = plhs[1] = mxCreateNumericArray(2,dims, mxUINT64_CLASS,0),
*X_out = plhs[2] = mxCreateDoubleMatrix(2,10000, mxREAL);
const unsigned long long
tracks = (const unsigned long long)mxGetPr(prhs[0])[0];
unsigned long long
*hot = (unsigned long long*)mxGetPr(hot_out),
*cold = (unsigned long long*)mxGetPr(cold_out);
double
*Xout = mxGetPr(X_out);
/* call the actual function, and return */
CountMuons(tracks, hot,cold, Xout);
}
// The actual muon counting
void CountMuons(
unsigned long long tracks,
unsigned long long *hot, unsigned long long *cold, double *Xout)
{
const double
width = 1.0e-4,
height = 2.0e-4,
bL = 15.0/sqrt(2.0),
Ystart = 15.0;
double
Xstart,
Xend,
u,
b;
unsigned long long
i = 0ul;
*hot = 0ul;
*cold = tracks;
/* seed the RNG */
srand((unsigned)time(NULL));
/* aaaand start! */
while (i++ < tracks)
{
u = rand_double() - 0.5;
Xstart = 800.0*rand_double() - 400.0;
Xend = Xstart - bL*sign(u)*log(1.0-2.0*fabs(u));
b = (Ystart*Xend)/(Xend-Xstart);
if ((b < height && b > 0.0) || (Xend < width/2.0 && Xend > -width/2.0))
{
Xout[0 + *hot*2] = Xstart;
Xout[1 + *hot*2] = Xend;
++(*hot);
--(*cold);
}
}
}
でMatlabでコンパイル
mex DoMuonCounting.c
( :)を実行した後mex setup、次のような小さな M-wrapper と組み合わせて使用します。
function [hot,cold, h] = MuonTrack2(tracks)
% call the MEX function
[hot,cold, Xtmp] = DoMuonCounting(tracks);
% process outputs, and generate plots
hot = uint32(hot); % circumvents limitations in 32-bit matlab
X = Xtmp(:,1:hot);
clear Xtmp
h = NaN;
if ~isempty(X)
h = figure;
hold all
Y = repmat([15;0], 1, hot);
plot(X, Y, 'r');
end
end
それは私がすることを可能にします
>> tic, MuonTrack2(1e8); toc
Elapsed time is 14.496355 seconds.
MEX 版のメモリ フットプリントはわずかに大きくなりますが、心配する必要はないと思います。
私が見る唯一の欠陥は、ミューオンカウントの固定最大数です（の初期配列サイズとして10000としてハードコーディングされていXoutます;標準Cには動的に成長する配列がないため必要です）...心配している場合、この制限は壊れている場合は、単純に増やすか、 の一部に等しくなるように変更するtracksか、よりスマートな (しかしより苦痛な) 動的な配列拡張のトリックを実行します。
