プログラムを改善するためにできることはたくさんあります-アルゴリズムとコードの両方です。
コード側では、本当に遅くなっていることの1つは、for
ループ(遅い)を使用するだけでなく、行でループを使用しているという事実です。
P = [P;P1];
配列に要素を追加します。そのたびに、Matlabはデータを配置する新しい場所を見つけ、プロセス内のすべてのポイントをコピーする必要があります。これはすぐに非常に遅くなります。でアレイを事前に割り当てる
P = zeros(1000000, 3);
これまでに見つけたポイントの数Nを追跡し、距離の計算を次のように変更します。
D = pdist2(P1, P(1:N, :), 'euclidean');
少なくともそれに対処します...
もう1つの問題は、以前に見つかったすべてのポイントに対して新しいポイントをチェックすることです。つまり、100ポイントの場合、約100x100をチェックし、1000の場合は1000x1000になります。この場合、このアルゴリズムは少なくともO(N ^ 3)であることがわかります...密度が上がるにつれて、より多くの「ミス」が発生するという事実は考慮されていません。N = 10E6のAO(N ^ 3)アルゴリズムには、少なくとも10E18サイクルかかります。比較ごとに1クロックサイクルの4GHzマシンがある場合、2.5E9秒=約80年が必要になります。並列処理を試すことはできますが、それは単なるブルートフォースです-誰がそれを望んでいますか?
問題をより小さな部分に分割することを検討することをお勧めします(文字通り):たとえば、球を最大距離とほぼ同じサイズの小さなボックスに分割し、各ボックスについて、どのポイントが入っているかを追跡する場合次に、「this」ボックスとそのすぐ隣のポイント(全部で27個のボックス)をチェックするだけで済みます。ボックスの幅が2.5mmの場合、100x100x100=1Mのボックスになります。これは大変なことのように思えますが、(アルゴリズムの終わりまでに)ボックスごとに平均1ポイントしかないため、計算時間が大幅に短縮されます...もちろん、使用している距離基準を使用すると、中心近くにより多くのポイントがありますが、それは詳細です。
必要なデータ構造は100x100x100のセル配列であり、各セルには「そのセル内」でこれまでに見つかった良い点のインデックスが含まれています。セル配列の問題は、ベクトル化に適していないことです。代わりにメモリがある場合は、セルごとに10ポイント以下であると仮定して、10x100x100x100の4D配列として割り当てることができます(そうする場合は、個別に処理する必要があります。ここで作業してください...) 。-1
まだ見つからないポイントのインデックスを使用する
次に、チェックは次のようになります。
% initializing:
bigList = zeros(10,102,102,102)-1; % avoid hitting the edge...
NPlist = zeros(102, 102, 102); % track # valid points in each box
bottomcorner = [-25.5, -25.5, -25.5]; % boxes span from -25.5 to +25.5
cellSize = 0.5;
。
% in your loop:
P1= [x, y, z];
cellCoords = ceil(P1/cellSize);
goodFlag = true;
pointsSoFar = bigList(:, cellCoords(1)+(-1:1), cellCoords(2)+(-1:1), cellCoords(3)+(-1:1));
pointsToCheck = find(pointsSoFar>0); % this is where the big gains come...
r=sum(P1.^2);
D = pdist2(P1,P(pointsToCheck, :),'euclidean'); % euclidean distance
if D>0.146*r^(2/3)
P(k,:) = P1;
% check the syntax of this line...
cci = ind2sub([102 102 102], cellCoords(1), cellCoords(2), cellCoords(3));
NP(cci)=NP(cci)+1; % increasing number of points in this box
% you want to handle the case where this > 10!!!
bigList(NP(cci), cci) = k;
k=k+1;
end
....
ここからそれを取ることができるかどうかはわかりません。できない場合は、メモでそのように言ってください。今週末、これをより詳細にコーディングする時間があります。いくつかのベクトル化でそれをさらにスピードアップする方法がありますが、すぐに管理するのが難しくなります。
より多くの点をランダムに空間に配置し、上記を巨大なベクトル化されたカリングに使用するのが良い方法かもしれないと思います。ただし、最初に少し手順を実行することをお勧めします...上記をうまく機能させることができれば、さらに最適化できます(アレイサイズなど)。