Monte Carlo (Possibly Simulated Annealing?) Method For N Mutually Repelling Points on a Unit Sphere C++

Question

I need to create an algorithm in C++ to simulate mutually repelling points on a sphere using a Monte Carlo method. So far what I have is this:

#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;

}

Which I think is ok (note I am essentially trying to minimise the energy function), but I want to make it more accurate/make it run quicker. To do so I think I should change my value of p, the while loop conditions or how to alter p at the end of the code. (All of these are in *... * as I was trying to embolden them to make it clear to you where I mean. About 3/4 of the way through the code). I have been sitting for hours trying to alter these conditions but nothing is working. For n=12 (12 points on the sphere) my energy should come out at 49.16525306, but I can only get it between 50.5 and 54.0 really. I know this is relatively good, but I want it more accurate (even if it does take a while). I would alsolike the success rate to increase if possible (my overall success rate it absolutely appalling).

If anyone has any ideas, I would be very grateful for your help!

Thanks, A.

(Note: If you want to run the code you must take the double *'s out. There are four sections with double *'s surrounding them).

Answer 1

First, you seem like an intelligent scientist/mathematician who is trying to do some programming. I'm a physicist, and in my experience such people make some of the worst programmers; if at all possible, get some help from an experienced coder.

Second, look at this code (which is repeated, see First):

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

You are modifying all three coordinates by the same amount, which means you always move a point along the (1,1,1) ray. The results improve if you modify one coordinate at a time.

Third, in the final loop (which is the one that takes most of the time) your logic is a little screwy-- you modify x, but then calculate energy using y. The results are still pretty good, because you also have x and y transposed at the end of the loop, but correcting this improves the accuracy of the results.

Fourth, and this is a big one, when you perturb a point and then recalculate energy, you recalculate the contributions of all points; only one point has changed, which means that most of the point pairs have not changed and need not be recalculated. Instead, after you choose a point, you can calculate the contribution of that point with something like this:

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

Then calculate it again after the perturbation, and compare. This takes the calculation from O(n²) to O(n), which makes it a lot faster.

When I make these modifications (and make p converge 10 times faster, because I'm not very patient) my energy comes out at 49.1652530576.