macos — метод Монте-Карло (возможно, имитация отжига?) для N взаимно отталкивающихся точек при переполнении стека единичной сферы

Мне нужно создать алгоритм в 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. Я знаю, что это относительно хорошо, но я хочу, чтобы это было более точным (даже если это займет некоторое время). Я также хотел бы увеличить вероятность успеха, если это возможно (моя общая вероятность успеха абсолютно ужасна).

Если у кого-то есть идеи, буду очень признателен за вашу помощь!

Спасибо.

(Примечание: если вы хотите запустить код, вы должны убрать двойные *. Есть четыре секции с двойными *, окружающими их).