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 之间得到它。我知道这相对较好，但我希望它更准确（即使确实需要一段时间）。如果可能的话，我还希望提高成功率（我的整体成功率绝对令人震惊）。

如果有人有任何想法，我将非常感谢您的帮助！

谢谢。

（注意：如果你想运行代码，你必须把双 * 去掉。有四个部分用双 * 包围它们）。

4

1 回答 1

3

首先，你看起来像一个聪明的科学家/数学家，正在尝试做一些编程。我是一名物理学家，根据我的经验，这些人是一些最糟糕的程序员；如果可能的话，请从经验丰富的编码人员那里获得一些帮助。

其次，看这段代码（重复，见First）：

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

您将所有三个坐标都修改了相同的量，这意味着您始终沿 (1,1,1) 射线移动一个点。如果一次修改一个坐标，结果会得到改善。

第三，在最后一个循环中（这是花费大部分时间的循环），您的逻辑有点混乱——您修改x，然后使用y计算能量。结果仍然相当不错，因为在循环结束时您还对 x 和 y 进行了转置，但纠正这一点可以提高结果的准确性。

第四，这是一个大问题，当你扰动一个点然后重新计算能量时，你重新计算所有点的贡献；只有一个点发生了变化，这意味着大部分点对没有变化，不需要重新计算。相反，在您选择一个点后，您可以使用以下方式计算该点的贡献：

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 回答