1

C++ でアルゴリズムを作成して、モンテカルロ法を使用して球面上の相互に反発する点をシミュレートする必要があります。これまでのところ、私が持っているのはこれです:

#include <stdio.h> 
#include <string.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <time.h>
#include <stdlib.h>
using namespace std;

int main()
{

  int a,f,g,n,m,i,j,k,r,s;
  double p,q,Energy,energy,y[101][4],x[101][4],Length,Distance;

 clock_t t1,t2;
  t1=clock();

  /*  set the number of points */
  n=12;

  /* check that there are no more than 100 points */
  if(n>100){
    cout << n << " is too many points for me :-( \n";
    exit(0);
  }

  /* reset the random number generator */
  srand((unsigned)time(0));  

  for (i=1;i<=n;i++){
    x[i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
    x[i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
    x[i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;

    Length=sqrt(pow(x[i][1],2)+pow(x[i][2],2)+pow(x[i][3],2));

    for (k=1;k<=3;k++){
      x[i][k]=x[i][k]/Length;
    }
  }

  /* calculate the energy */
  Energy=0.0;

  for(i=1;i<=n;i++){
    for(j=i+1;j<=n;j++){
      Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                    +pow(x[i][3]-x[j][3],2));

      Energy=Energy+1.0/Distance;
    }
  }

  /* Save Original Points */
  for(i=1;i<=n;i++){
    y[i][1]=x[i][1];
    y[i][2]=x[i][2];
    y[i][3]=x[i][3];
  }

  /* Loop for random points m times*/
  m=10;

  if (m>100){
    cout << "The m="<< m << " loop is inefficient...lessen m \n";
    exit(0);
  }

  a=1;

  while(a<m){

    /* assign random points */
    for (i=1;i<=n;i++){
      x[i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;

      Length=sqrt(pow(x[i][1],2)+pow(x[i][2],2)+pow(x[i][3],2));

      for (k=1;k<=3;k++){
        x[i][k]=x[i][k]/Length;
      }
    }

    /* calculate the energy */
    energy=0.0;

    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                      +pow(x[i][3]-x[j][3],2));

        energy=energy+1.0/Distance;
      }
    }

    if(energy<Energy)
      for(i=1;i<=n;i++){
        for(j=1;j<=3;j++){
          Energy=energy;
          y[i][j]=x[i][j];
        }
      }
    else
      for(i=1;i<=n;i++){
        for(j=1;j<=3;j++){
          energy=Energy;
          x[i][j]=y[i][j];
        }
      }

    a=a+1;
  }

  /* Output the best random energy */
  cout << "Energy=" << Energy << "\n";

  m=10;
  a=1;

  while(a<m){
    /* Choose random point to move */
    g=(rand() % n)+1;

    /* Choose a p small to give q in a range -p <= q <= p */
    p=0.1;

    /* q is how much I am moving the random point by */
    q=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0*p;

    /* Move the point by q */
    for(j=1;j<=3;j++){
      x[g][j]=((x[g][j])+q);
    }

    /* Bring it back onto sphere */
    Length=sqrt(pow(x[g][1],2)+pow(x[g][2],2)+pow(x[g][3],2));

    for (k=1;k<=3;k++){
      x[g][k]=x[g][k]/Length;
    }

    /* Calculate the new energy */
    energy=0.0;

    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                         +pow(x[i][3]-x[j][3],2));

        energy=energy+1.0/Distance;
      }
    }

    /* Choose best energy and therefore best point */
    if (energy<Energy)
      Energy=energy,x[g][1]=y[g][1],x[g][2]=y[g][2],x[g][3]=y[g][3];
    else
      energy=Energy,y[g][1]=x[g][1],y[g][2]=x[g][2],y[g][3]=x[g][3];

    a=a+1;  

  }

   /* Output the best single shift energy */
  cout << "Energy=" << Energy << "\n";

  /* Set fail count to 0 */
  s=0;
  f=0;
  r=1;
  **p=0.1;**

  /* Maximum distance to move the random point */

  while (**p>0.00001**) {

    /* Number of loops to do */

    while (**r<3000**) {

      g=(rand() % n)+1;

      /* q is how much I am moving the random point by -p<=q<=p*/
      q=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0*p;

      /* Move the point by q */
      for(j=1;j<=3;j++){
        x[g][j]=((x[g][j])+q);
      }

      /* Bring it back onto sphere */
      Length=sqrt(pow(x[g][1],2)+pow(x[g][2],2)+pow(x[g][3],2));

      for (k=1;k<=3;k++){
        x[g][k]=x[g][k]/Length;
      }

      /* Calculate the new energy */
      energy=0.0;

      for(i=1;i<=n;i++){
        for(j=i+1;j<=n;j++){
          Distance=sqrt(pow(y[i][1]-y[j][1],2)+pow(y[i][2]-y[j][2],2)
                        +pow(y[i][3]-y[j][3],2));
          energy=energy+1.0/Distance;
        }
      }

      /* Choose best energy and therefore best point */
      if (energy<Energy)
        Energy=energy,x[g][1]=y[g][1],x[g][2]=y[g][2],x[g][3]=y[g][3],s=s+1;
      else
        energy=Energy,y[g][1]=x[g][1],y[g][2]=x[g][2],y[g][3]=x[g][3],f=f+1;

      r=r+1;

    }

    **/* Calculate percentage fails */

    if ((100.0*(f/r))>50.0)
      p=(p-0.00001);
    else
      p=p;**

    r=0;  

  }

  cout << "Overall Success Rate = " << ((s*1.0)/((s+f)*1.0))*100 << "%" << "\n";
  cout << "Energy=" << fixed << setprecision(10) << Energy << "\n";


  ofstream Bestpointssofar ("Bestpointssofar");
  for(i=1;i<=n;i++){
    Bestpointssofar << y[i][1] << " " <<   y[i][2] << " " << y[i][3] << "\n";
  }
  Bestpointssofar.close(); 

  t2=clock();
    float diff ((float)t2-(float)t1);
    float seconds = diff / CLOCKS_PER_SEC;
    cout << fixed << setprecision(2) << "Run time: " << seconds << "(s)" << "\n";
    return 0;

}

私は大丈夫だと思います（私は基本的にエネルギー関数を最小化しようとしていることに注意してください）が、より正確にしたり、より速く実行したりしたいと考えています。そのためには、p の値、while ループ条件、またはコードの最後で p を変更する方法を変更する必要があると思います。(これらはすべて * ... *に含まれています。これは、私が意味するところを明確にするためにそれらを大胆にしようとしていたためです。コードの約 3/4 まで)。私はこれらの条件を変えようと何時間も座っていましたが、何も機能していません. n=12 (球上の 12 ポイント) の場合、私のエネルギーは 49.16525306 になるはずですが、実際には 50.5 から 54.0 の間しか得られません。これが比較的良いことはわかっていますが、もっと正確にしたいです (時間がかかるとしても)。また、可能であれば成功率を上げたいと思っています (私の全体的な成功率は絶対に恐ろしいものです)。

誰かが何かアイデアを持っているなら、私はあなたの助けにとても感謝しています!

ありがとう、A.

(注: コードを実行する場合は、二重の * を削除する必要があります。二重の * で囲まれた 4 つのセクションがあります)。

4

1 に答える 1

3

まず、あなたはプログラミングをしようとしている知的な科学者/数学者のように見えます。私は物理学者であり、私の経験では、そのような人々は最悪のプログラマーです。可能であれば、経験豊富なコーダーの助けを借りてください。

次に、このコードを見てください (これは繰り返されます。最初を参照してください)。

/* Move the point by q */
for(j=1;j<=3;j++){
  x[g][j]=((x[g][j])+q);
}

3 つの座標すべてを同じ量だけ変更します。つまり、常に (1,1,1) 光線に沿ってポイントを移動します。一度に 1 つの座標を変更すると、結果が改善されます。

3 番目に、最後のループ (最も時間がかかるループ) では、ロジックが少し厄介です。 xを変更しますが、 yを使用してエネルギーを計算します。ループの最後で x と y が転置されているため、結果は依然としてかなり良好ですが、これを修正すると結果の精度が向上します。

第 4 に、これは大きな問題です。ポイントを摂動してからエネルギーを再計算すると、すべてのポイントの寄与が再計算されます。変更されたポイントは 1 つだけです。つまり、ほとんどのポイント ペアは変更されておらず、再計算する必要はありません。代わりに、ポイントを選択した後、次のようにしてそのポイントの寄与を計算できます。

double oldEnergy = 0.0;
  for(i=1;i<=n;i++)
    {
      if(i!=g)
        {
          Distance=myDistance(x[i], x[g]);
          oldEnergy += 1.0/Distance;
        }
    }

次に、摂動後に再度計算し、比較します。これにより、O(n 2 ) から O(n) への計算が行われるため、計算が大幅に高速化されます。

これらの変更を行うと (そして p を 10 倍速く収束させます。なぜなら、私はあまり忍耐強くないからです) 、私のエネルギーは 49.1652530576 になります。

于 2013-01-01T02:14:44.720 に答える