Possibly misunderstanding exit(0) statement.

[adrem7]

I have written my code to run a simulation of evenly distributing n points on a sphere. I have put a failsafe in place to end the code if there are any issues. My code is as follows:

#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,i,j,k,m,n,s;
  double p,Energy,energy,Distance,Length,DotProdForce,Forcemagnitude,ForceMagnitude[101],
         Force[101][4],E[10001],x[101][4],y[101][4];

  clock_t t1,t2;
  t1=clock();

  /*  Set the number of points */
  n=25;

  /* 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 these points */
  ofstream Bestrandompoints ("Bestrandompoints");
  Bestrandompoints << "Energy=" << Energy << "\n";
  for(i=1;i<=n;i++){
    Bestrandompoints << x[i][1] << " " <<   x[i][2] << " " << x[i][3] << "\n";
  }
  Bestrandompoints.close(); 

  /* Start doing gradient flow approach */
  a=1;
  s=0;
  f=0;
  p=0.05;
  Forcemagnitude=1.0;

  while(Forcemagnitude>0.00001){

    /* Reset initial force and energy change */
    for(i=1;i<=n;i++){ 
      Force[i][1]=0.0; 
      Force[i][2]=0.0; 
      Force[i][3]=0.0; 
    } 
    
    /* Calculate force on each particle */
    for(i=1;i<=n;i++){ 
      for(j=1;j<i;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));
        Force[i][1]=Force[i][1]+((x[i][1]-x[j][1])/(pow((Distance),3))); 
        Force[i][2]=Force[i][2]+((x[i][2]-x[j][2])/(pow((Distance),3))); 
        Force[i][3]=Force[i][3]+((x[i][3]-x[j][3])/(pow((Distance),3))); 
      } 

      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));
        Force[i][1]=Force[i][1]+((x[i][1]-x[j][1])/(pow((Distance),3))); 
        Force[i][2]=Force[i][2]+((x[i][2]-x[j][2])/(pow((Distance),3))); 
        Force[i][3]=Force[i][3]+((x[i][3]-x[j][3])/(pow((Distance),3)));
      }
    }

    /* Add the force to my points */
    for(i=1;i<=n;i++){

      DotProdForce=Force[i][1]*x[i][1]+Force[i][2]*x[i][2]+Force[i][3]*x[i][3];
      
      y[i][1]=x[i][1];
      y[i][2]=x[i][2];
      y[i][3]=x[i][3];
    
      Force[i][1]=Force[i][1]-DotProdForce*y[i][1];
      Force[i][2]=Force[i][2]-DotProdForce*y[i][2];
      Force[i][3]=Force[i][3]-DotProdForce*y[i][3];

      x[i][1] = y[i][1]+p*(Force[i][1]);
      x[i][2] = y[i][2]+p*(Force[i][2]);
      x[i][3] = y[i][3]+p*(Force[i][3]);

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

      for (j=1;j<=3;j++){
        x[i][j]=x[i][j]/Length;
      }
    }
      
    /* Calculate the new energy */
    E[a]=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));
        E[a]=E[a]+1.0/Distance;
      }
    }

    energy=E[a];

    cout << fixed << setprecision (10) << "Energy=" << Energy << "\n";
    cout << fixed << setprecision (10) << "energy=" << energy << "\n";

    /* Choose best energy and therefore best points */
    if (energy<Energy)
      Energy=energy,s=s+1;
    else
      energy=Energy,f=f+1;
  
    if(f=1){
      cout << fixed << setprecision (10) << "E=" << energy << "\n";
      cout << "Successes=" << s << " Failures=" << f <<"\n";
      cout << "Terminated \n";
      exit(0);
    }  
  
    cout << "Successes=" << s << " Failures=" << f <<"\n";
  
    for(i=1;i<=n;i++){
      ForceMagnitude[i]=pow((pow(Force[i][1],2)+pow(Force[i][2],2)
                             +pow(Force[i][3],2)),0.5);
    }

    for(i=1;i<=n-1;i++){
      if(ForceMagnitude[i]<ForceMagnitude[i+1])
        ForceMagnitude[i]=ForceMagnitude[i+1];
      else
        ForceMagnitude[i+1]=ForceMagnitude[i];
    }

    Forcemagnitude=ForceMagnitude[n];

    cout << "Forcemagnitude=" << fixed << setprecision(10) << Forcemagnitude << "\n";

    a=a+1;
    
  }

  /* Output new energy */
  cout << "Successes=" << s << " Failures=" << f <<"\n";
  cout << fixed << setprecision (10) << "Energy=" << Energy << "\n";
  cout << "Iterations=" << a-1 << "\n";

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

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

}

My failsafe is at line 213 and I want it to basically exit the code if energy I am calculating stops decreasing. It does this without fail, but my issue is that my outputs are confusing me.

A standard output gives me something like:

Energy=269.7821451992
energy=266.7549548756
E=266.7549548756
Successes=1 Failures=1
Terminated

As you can see, the Energy is higher than the energy and so this is a success, hence 1 success. However where has the failure come from?

The code states (I think) that on lines 204/205 I should get an output of my Energy and energy and then increment f if there is a failure. Since I only have one set of Energy and energy outputs, surely my successes and failures should add up to one. The clearly add to two and so therefore the while loop has run twice. Why is it only giving my one set of outputs?

Thanks for any assistance you can offer een if it is just links or extra reading. The only reason I can think of is either I have misunderstood exit(0) and so the code is exiting before I think it is, or I have misunderstood the order in which my code is being read and therefore expect more outputs than there should be.

Note: My code does get it right sometimes, but only if there is a fail on the first go. Then I get this output:

Energy=275.8325665982
energy=313.1149926643
E=275.8325665982
Successes=0 Failures=1
Terminated

(Obviously the Energy and energy values differ every time you run the code, but that is an example).

Thanks,

A.

[cire]

My failsafe is at line 213

You might want to use the comparison operator (==) as opposed to the assignment operator in that if condition.

[adrem7]

Excellent, such a stupid mistake. Thanks for that!