c++ - 单位球体 C++ 上 N 个相互排斥点的蒙特卡罗(可能是模拟退火?)方法

标签 c++ macos montecarlo simulated-annealing

我需要在 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。

(注意:如果你想运行代码,你必须把双 * 去掉。有四个部分被双 * 包围)。

最佳答案

首先,您看起来像是一位聪明的科学家/数学家,正在尝试进行一些编程。我是一名物理学家,根据我的经验,这些人会成为最糟糕的程序员;如果可能的话,从有经验的编码人员那里获得一些帮助。

其次,看这段代码(有重复,见第一):

/* 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(n2) 到 O(n),这使得它更快。

当我进行这些修改时(并使 p 收敛速度快 10 倍,因为我不是很有耐心)我的能量在 49.1652530576 处出现。

关于c++ - 单位球体 C++ 上 N 个相互排斥点的蒙特卡罗(可能是模拟退火?)方法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/14084545/

相关文章:

python - Tensorflow 可以计算积分近似的梯度吗?

c - OpenMP Monte_Carlo 模拟实现目标与 PI 的接近度

c++ - 使用我的结构作为返回类型的问题

c++ - 为什么要打扰预处理器指令?

java - Minecraft Coder Pack终端内存错误伪造[Mac]

swift - 快速连续点击不同的按钮

c++ - 在声明前使用静态常量成员

android - "fatal error: Eigen/Dense: No such file or directory"//特征库

macos - 如何检测我的应用程序是否已沙盒化?

matlab - 使用蒙特卡罗方法将矩阵元素扰动到其标称值周围约 10%