4

そこの。この特定の文字を使用して、MATLABで10^6個のランダムなポイントを生成します。ポイントはradious25の球の内側にある必要があり、3Dであるため、x、y、zまたはr、theta、phiがあります。各ポイント間には最小距離があります。まず、ポイントを生成してから距離を確認し、これらの条件がないポイントを除外することにしました。ただし、多くの点が省略される場合があります。もう1つの方法は、RSA(Random Sequential Addition)を使用することです。これは、ポイント間の最小距離でポイントを1つずつ生成することを意味します。たとえば、最初のポイントを生成し、次にポイント1からの最小距離からランダムに2番目のポイントを生成し、10^6ポイントに達するまで続けます。しかし、新しいポイントを探すのに時間がかかるので、時間がかかり、10^6ポイントに到達できません。

現在、私はこのプログラムを使用しています。

Nmax=10000; 
R=25; 
P=rand(1,3); 
k=1; 
while k<Nmax 
theta=2*pi*rand(1); 
phi=pi*rand(1); 
r = R*sqrt(rand(1)); 
% convert to cartesian 
x=r.*sin(theta).*cos(phi); 
y=r.*sin(theta).*sin(phi); 
z=r.*cos(theta); 
P1=[x y z]; 
r=sqrt((x-0)^2+(y-0)^2+(z-0)^2); 
D = pdist2(P1,P,'euclidean'); 
% euclidean distance 

    if D>0.146*r^(2/3) 
        P=[P;P1]; 
        k=k+1;
    end 
    i=i+1; 
end 
x=P(:,1);y=P(:,2);z=P(:,3); plot3(x,y,z,'.');

これらの条件で効率的にポイントを生成するにはどうすればよいですか?ありがとうございました。

4

3 に答える 3

4

私はあなたのアルゴリズムを詳しく調べて、それが機能する方法はないと結論付けました-少なくとも、その領域で本当に100万ポイントを獲得したい場合はそうではありません。理由を説明する簡単な図があります。これは、1つの追加の「適切な」ポイントを取得するために(RSAの手法を使用して)テストする必要があるポイントの数のプロットです。ご覧のとおり、これはわずか数千ポイントで漸近的になります(これを生成するために、200kポイントに対してわずかに高速なアルゴリズムを実行しました)。

テストされたポイントのグラフ

球を完全に配置したときに球に収まる理論上の点の数を計算しようとしたことがあるかどうかはわかりませんが、その数は1E6よりもかなり小さいと思い始めています。

これを調査するために使用した完全なコードと、それが生成した出力は、ここにあります。以前の回答で説明したテクニックまでは到達できませんでした...あなたが説明したセットアップでは、他にも多くのことが起こっていました。

編集:「完璧な」配置でも、100万ポイントに到達することは不可能かもしれないと私は考え始めました。私は次のように自分用の簡単なモデルを作成しました。

「外殻」(r = 25)から始めて、等距離の点を合わせようとしていると想像してください。「シェル」の面積を(半径r_sub_critの)1つの「除外ディスク」の面積で割ると、その距離にあるポイント数の(高い)推定値が得られます。

numpoints = 4*pi*r^2 / (pi*(0.146 * r^(2/3))^2) ~ 188 * r^(2/3)

次の「シェル」は、半径が0.146 * r ^(2/3)小さいはずですが、ポイントが非常に注意深く配置されていると考えると、少し近づくことができる場合があります。繰り返しになりますが、寛大に、シェルが基準よりも1 / sqrt(3)近くなる可能性があるとしましょう。次に、簡単なPythonスクリプトを使用して、外殻から始めて作業を進めることができます。

import scipy as sc
r = 25
npts = 0
def rc(r):
  return 0.146*sc.power(r, 2./3.)
while (r > rc(r)):  
  morePts = sc.floor(4/(0.146*0.146)*sc.power(r, 2./3.))
  npts = npts + morePts
  print morePts, ' more points at r = ', r
  r = r - rc(r)/sc.sqrt(3)

print 'total number of points fitted in sphere: ', npts

これの出力は次のとおりです。

1604.0  more points at r =  25
1573.0  more points at r =  24.2793037966
1542.0  more points at r =  23.5725257555
1512.0  more points at r =  22.8795314897
1482.0  more points at r =  22.2001865995
1452.0  more points at r =  21.5343566722
1422.0  more points at r =  20.8819072818
1393.0  more points at r =  20.2427039885
1364.0  more points at r =  19.6166123391
1336.0  more points at r =  19.0034978659
1308.0  more points at r =  18.4032260869
1280.0  more points at r =  17.8156625053
1252.0  more points at r =  17.2406726094
1224.0  more points at r =  16.6781218719
1197.0  more points at r =  16.1278757499
1171.0  more points at r =  15.5897996844
1144.0  more points at r =  15.0637590998
1118.0  more points at r =  14.549619404
1092.0  more points at r =  14.0472459873
1066.0  more points at r =  13.5565042228
1041.0  more points at r =  13.0772594652
1016.0  more points at r =  12.6093770509
991.0  more points at r =  12.1527222975
967.0  more points at r =  11.707160503
943.0  more points at r =  11.2725569457
919.0  more points at r =  10.8487768835
896.0  more points at r =  10.4356855535
872.0  more points at r =  10.0331481711
850.0  more points at r =  9.64102993012
827.0  more points at r =  9.25919600154
805.0  more points at r =  8.88751153329
783.0  more points at r =  8.52584164948
761.0  more points at r =  8.17405144976
740.0  more points at r =  7.83200600865
718.0  more points at r =  7.49957037478
698.0  more points at r =  7.17660957023
677.0  more points at r =  6.86298858965
657.0  more points at r =  6.55857239952
637.0  more points at r =  6.26322593726
618.0  more points at r =  5.97681411037
598.0  more points at r =  5.69920179546
579.0  more points at r =  5.43025383729
561.0  more points at r =  5.16983504778
542.0  more points at r =  4.91781020487
524.0  more points at r =  4.67404405146
506.0  more points at r =  4.43840129415
489.0  more points at r =  4.21074660206
472.0  more points at r =  3.9909446055
455.0  more points at r =  3.77885989456
438.0  more points at r =  3.57435701766
422.0  more points at r =  3.37730048004
406.0  more points at r =  3.1875547421
390.0  more points at r =  3.00498421767
375.0  more points at r =  2.82945327223
360.0  more points at r =  2.66082622092
345.0  more points at r =  2.49896732654
331.0  more points at r =  2.34374079733
316.0  more points at r =  2.19501078464
303.0  more points at r =  2.05264138052
289.0  more points at r =  1.91649661498
276.0  more points at r =  1.78644045325
263.0  more points at r =  1.66233679273
250.0  more points at r =  1.54404945973
238.0  more points at r =  1.43144220603
226.0  more points at r =  1.32437870508
214.0  more points at r =  1.22272254805
203.0  more points at r =  1.1263372394
192.0  more points at r =  1.03508619218
181.0  more points at r =  0.94883272297
170.0  more points at r =  0.867440046252
160.0  more points at r =  0.790771268402
150.0  more points at r =  0.718689381062
140.0  more points at r =  0.65105725389
131.0  more points at r =  0.587737626612
122.0  more points at r =  0.528593100237
113.0  more points at r =  0.473486127367
105.0  more points at r =  0.422279001431
97.0  more points at r =  0.374833844693
89.0  more points at r =  0.331012594847
82.0  more points at r =  0.290676989951
75.0  more points at r =  0.253688551418
68.0  more points at r =  0.219908564725
61.0  more points at r =  0.189198057381
55.0  more points at r =  0.161417773651
49.0  more points at r =  0.136428145311
44.0  more points at r =  0.114089257597
38.0  more points at r =  0.0942608092113
33.0  more points at r =  0.0768020649149
29.0  more points at r =  0.0615717987589
24.0  more points at r =  0.0484282253244
20.0  more points at r =  0.0372289153633
17.0  more points at r =  0.0278306908104
13.0  more points at r =  0.0200894920319
10.0  more points at r =  0.013860207063
8.0  more points at r =  0.00899644813842
5.0  more points at r =  0.00535025545232 

total number of points fitted in sphere:  55600.0

これは、あなたがどのように試みても、あなたが本当に百万に達することができないことを確認しているようです...

于 2013-02-09T14:50:57.577 に答える
3

プログラムを改善するためにできることはたくさんあります-アルゴリズムとコードの両方です。

コード側では、本当に遅くなっていることの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
....

ここからそれを取ることができるかどうかはわかりません。できない場合は、メモでそのように言ってください。今週末、これをより詳細にコーディングする時間があります。いくつかのベクトル化でそれをさらにスピードアップする方法がありますが、すぐに管理するのが難しくなります。

より多くの点をランダムに空間に配置し、上記を巨大なベクトル化されたカリングに使用するのが良い方法かもしれないと思います。ただし、最初に少し手順を実行することをお勧めします...上記をうまく機能させることができれば、さらに最適化できます(アレイサイズなど)。

于 2013-02-09T00:17:31.240 に答える
2

は参考文献を見つけました-「三次元セルオートマトンを使用したシミュレートされた脳腫瘍成長ダイナミクス」、Ansal et al(2000)。

私はそれが不可解であることに同意します-あなたが1つの重要なことに気付くまで。彼らは結果をで報告してmmいますが、あなたのコードはで書かれていますcm。それは取るに足らないように思えるかもしれませんが、「臨界半径」の式には、次元でrc = 0.146r^(2/3)ある定数、が含まれています。0.146次元はmm^(1/3)、ではなく、cm^(1/3)です。

可能なラティスサイトの数を評価するためにPythonコードにその変更を加えると、10倍にジャンプします。現在、彼らは0.38の「ジャミング制限」を使用していると主張しました。これは実際にはこれ以上サイトを見つけることができない数です。 。その制限を含めると、200kポイントしか見つからないと予測します。それでも、150万ポイントには達していませんが、それほどクレイジーではありません。

これについて彼らと話し合うために著者に連絡することを検討するかもしれませんか?会話に私を含めたい場合は、次のアドレスに電子メールを送信できます:SO(2文字のみ)my handle nameドットunited states。上記のリンクを投稿したのと同じドメイン...

于 2013-02-11T13:10:42.353 に答える