C++

Forum

 
Gradient Flow/Steepest Descent method for mutually repelling points on a sphere

adrem7 (41)
I have written some code to give me an output of the energy associated with n points constrained to a sphere all repelling each other (much like electrons around an atom).

The issue is for certain points, the Forcemagnitude does not tend to zero after a a certain amount of run time. Does anyone know why? This issue means the code just runs indefinitely and will never cease. I can't think why it would do this. I can get a "close" answer if I alter the 0.00001 bit in the while loop, but this shouldn't have to be done so confuses me. Also I need a piece of code where you cn plug any n from 2 to 100 in and so altering the Forcemagnitude bit for a few points (n=7, 13 etc.) isn't really a solution.

The code is as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
 

#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 Forcemagnitude,ForceMagnitude[101],DotProdForce,Force[101][4],p,q[101],b[101],c[101],d[101],l[101],energy,Energy,Distance,E[10001],F[101][101][4],y[101][4],x[101][101][4],z[101][4],length[101],distance[101];

clock_t t1,t2;
    t1=clock();

  /*  set the number of points */
  n=72;

  /* check that there are no more than 100 points */
  if(n>100){
    cout << n << " is too many points for me :-( \n";
  exit(0);
  }
  
  /* Loop for random points m times*/
  m=10;
  
  if (m>100){
  cout << "The m="<< m << " loop is inefficient...lessen m \n";
  exit(0);
  }
  
  /* reset the random number generator */
  srand((unsigned)time(0));  
  
   for (a=1;a<=m;a++){

  /* assign random points */
  for (i=1;i<=n;i++){
      x[a][i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[a][i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[a][i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      length[a]=sqrt(pow(x[a][i][1],2)+pow(x[a][i][2],2)+pow(x[a][i][3],2));
      for (k=1;k<=3;k++){x[a][i][k]=x[a][i][k]/length[a];}
    }

  /* calculate the energy */
  E[a]=0.0;
  
  for(i=1;i<=n;i++){
  for(j=i+1;j<=n;j++){
      distance[a]=sqrt(pow(x[a][i][1]-x[a][j][1],2)+pow(x[a][i][2]-x[a][j][2],2)
		    +pow(x[a][i][3]-x[a][j][3],2));
      E[a]=E[a]+1.0/distance[a];
    }}

}

/* Find the lowest of these energies */

  for(a=1;a<m;a++){
  if (E[a]<E[a+1])
  for(i=1;i<=n;i++){
  E[a+1]=E[a],x[a+1][i][1]=x[a][i][1],x[a+1][i][2]=x[a][i][2],x[a+1][i][3]=x[a][i][3];
  }
  else
  for(i=1;i<=n;i++){
  E[a]=E[a+1],x[a][i][1]=x[a+1][i][1],x[a][i][2]=x[a+1][i][2],x[a][i][3]=x[a+1][i][3];
  }}
  
  Energy=E[a];

/* Save the points to separate variable */

  for(i=1;i<=n;i++){
  y[i][1]=x[a][i][1];
  y[i][2]=x[a][i][2];
  y[i][3]=x[a][i][3];
  }
  
  /* Output these points */
  
  ofstream Bestrandompoints ("Bestrandompoints");

      Bestrandompoints << "E=" << E[a] << "\n";
  for(i=1;i<=n;i++){
      Bestrandompoints << y[i][1] << " " <<   y[i][2] << " " << y[i][3] << "\n";
      }

   Bestrandompoints.close(); 

/* Start doing gradient flow approach */
r=1;
s=0;
f=0;
p=0.01;
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 (j=1;j<=n;j++){ 
for (i=1;i<j;i++){

Distance=sqrt(pow(y[j][1]-y[i][1],2)+pow(y[j][2]-y[i][2],2)+pow(y[j][3]-y[i][3],2));

Force[j][1]=Force[j][1]+((y[j][1]-y[i][1])/(pow((Distance),3))); 
Force[j][2]=Force[j][2]+((y[j][2]-y[i][2])/(pow((Distance),3))); 
Force[j][3]=Force[j][3]+((y[j][3]-y[i][3])/(pow((Distance),3))); 
} 

for (i=j+1;i<=n;i++){ 

Distance=sqrt(pow(y[j][1]-y[i][1],2)+pow(y[j][2]-y[i][2],2)+pow(y[j][3]-y[i][3],2));

Force[j][1]=Force[j][1]+((y[j][1]-y[i][1])/(pow((Distance),3))); 
Force[j][2]=Force[j][2]+((y[j][2]-y[i][2])/(pow((Distance),3))); 
Force[j][3]=Force[j][3]+((y[j][3]-y[i][3])/(pow((Distance),3)));
}}

/* Add the force to my points */

for(i=1;i<=n;i++){

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

b[i]=y[i][1];
c[i]=y[i][2];
d[i]=y[i][3];

Force[i][1]=Force[i][1]-DotProdForce*b[i];
Force[i][2]=Force[i][2]-DotProdForce*c[i];
Force[i][3]=Force[i][3]-DotProdForce*d[i];

y[i][1] = b[i]+p*(Force[i][1]);
y[i][2] = c[i]+p*(Force[i][2]);
y[i][3] = d[i]+p*(Force[i][3]);

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

for (j=1;j<=3;j++){
y[i][j]=y[i][j]/length[i];
}}

/* Calculate the new energy */
  E[r]=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));
E[r]=E[r]+1.0/Distance;
    }}

energy=E[r];

/* Choose best energy and therefore best points */
if (energy<Energy)
  Energy=energy,s=s+1;
  else
  energy=Energy,f=f+1;
  
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];

r=r+1;
}

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

ofstream mypoints ("mypoints");

for(i=1;i<=n;i++){
mypoints << y[i][1] << " " <<   y[i][2] << " " << y[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;

}


If you can see any issues please let me know, this is my first post on the forum and I have only been coding from about half way into term (about 6-8 weeks) so I wouldn't be surprised if I have done something wrong.

Thanks for any help in advance,

A.
Last edited on
TheIdeasMan (6781)
One potential problem is the direct comparison of doubles. You need an absolute value compared with a scaled epsilon value to do this correctly. There is a numeric_limits<double>epsilon Google C++ floating point comparison I am not sure whether this will fix your problem.

By far the easiest way to find out why things don't work, is to use a debugger. Should be easy if you are using an IDE. You can make a watch list of variables and track their values as you step through the code 1 line at a time, thus deducing where the problem is. There are command line ones like gdb, but these are a bit trickier to learn. There is an article on how to use it in the articles section on this site.

Some other comments:

Make more use of functions, rather than having 200 LOC in main. Declare the functions before main & put the definitions after main. Put some comments to describe what they do, and / or give them meaningful names.

You could then make more use of variables local to the function, with out having so many in main. The other thing about this, is wondering whether they are all initialised.

If you are going to use single char variable names, explain what they mean with comments at declaration. Edit - I can see that a lot of them a used in for loops as counters - in that case declare them just before the loop.

Variable names like Energy & energy can cause confusion. And Distance & distance - see below.

Avoid using namespace std; - your were lucky with your variable Distance, because there is a std::distance which is a very different thing. I have just seen you have a distance variable as well even though it is an array, it might cause problems. I try to avoid naming my variables the same as something in the std namespace. You could create your own namespace, but I think you could get away with slightly more imaginative variable names.

So delete line 10, and put std:: before each std thing, as in std::cout. If you look in the reference section, you can see all the STL things which require this. You can do using std::string; or using std::vector; say for things which you might use a lot.

Line 16: Avoid putting so many declarations on 1 line.

With your for loops, you have a lot of these:

for(i=1;i<=n;i++){

IMO, you should try to get used to starting at zero. the idiom for doing something 10 times is:

for (int a = 0 ; a <10; a++){}

Sorry to point out something so basic, but the use of <= of ten causes the overrun of arrays.

Also be wary of divisions - you should check with code that the divisor is not zero (make sure to use the epsilon)

Hope this helps. Cheers


Last edited on
Lowest0ne (1536)
Your loop (line 104-132) never updates ForceMagnitude.

As far as entering a number, the simplest thing is:
1
2
3
 
int i;
cout << "Enter a number: ";
cin >> i;


You will need/want to ask for a number for as long as the number is not within your range. You can also break cin this way, but this does seem like something you are doing for fun, so maybe who cares.
Last edited on
adrem7 (41)
Your loop (line 104-132) never updates ForceMagnitude.


Why should it need to? "ForceMagnitude" is defined later on as an equation of "Force" and Force updates in that loop so it should update just from that no?

Like I said I am new to C++ (relatively) and haven't actually done any coding before this project so things like arrays starting at 1 instead of zero are just me being too mathsy brained rather than thinking as a programmer.

I don't think these little things are the issue though (although I appreciate the code is messy and shall clean it up in the way you have suggested. - Thanks for all the tips!).

One potential problem is the direct comparison of doubles. You need an absolute value compared with a scaled epsilon value to do this correctly. There is a numeric_limits<double>epsilon Google C++ floating point comparison I am not sure whether this will fix your problem.


I shall look into the direct comparison of doubles thing and see if that makes it any better. Are you meaning the lines 181 to 186 for this?

Just a couple more questions:
Why can I not directly compare doubles?
Does either a Mac (using a MacBook Pro) or C++ only count to a certain number of decimal places?

If so comparing doubles directly like you say might mean that if the Forcemagnitude doesn't converge quickly enough it might blow up in my face like it has been doing.

Thanks for the help!

A.

Edit: Found this online for anyone who stumbles across this in future.
http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
Scroll down to the system aspects bit and from "Incidentally..." it gives a nice summary.
Also if the code works for most points why does it cease to work at n=30, surely the double comparison thing would be an issue for all n not just some of them. Since the code doesn't blow up in my face every time, does this not mean the comparisons I have put in place are acceptable?
Last edited on
TheIdeasMan (6781)
Are you meaning the lines 181 to 186 for this?


Yes, lines 172 & 182. You use < which might be OK, but see below where this might not be what you want.

That article was quite technical, I will explain in laymans terms.

Because FP numbers are represented as binary fractions, they cannot represent all real numbers exactly. So you finish up with the final value being the target value plus or minus a small amount. For a double of 1.0 , this small amount is about 1E-16 and is known as epsilon. It can be thought of as the 'distance' to the next representable number. It has to be scaled (because it varies depending on the size of the number)- if your number is 1000.0 then the scaled epsilon would be 1000.0 * std::numeric_limits<double>epsilon . You can use an absolute function to see if the difference between your number and the target number is less than the scaled epsilon. If so then they are 'equal'.

So all this becomes a problem when numbers are nearly equal, and not for all the values. It is a problem when you have a loop end condition or if statement depending on it. Consider this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
float a = 0.1; // a is 0.09999997
float b = 10 * a; // b is 0.9999997

if (b == 1.0)  //false but should really be true
if (a < 1.0)  //true for this number but not what you want
                    //may be false for other numbers because it is slightly larger

double MyDecPlaces = 3;
double Epsilon = std::numeric_limits<double>epsilon() ;
double MyNumber = 1000.0;

double MyEpsilon = MyNumber * Epsilon * MyDecPlaces;

b = 100.0 * 10.0;
if (std::abs(b - MyNumber) < MyEpsilon )  //true within MyDecPlaces


I used floats in the example to show what happens to the actual value, and I am saying changing them to double doesn't help.

Check out this:

http://en.cppreference.com/w/cpp/types/numeric_limits/epsilon


My example is similar, but without the template.

How did you get on with the other things I suggested - especially the debugger?
ne555 (10692)
I would blame `catastrophic cancellation'
1
2
3
 
Force[i][1]=Force[i][1]-DotProdForce*b[i];
Force[i][2]=Force[i][2]-DotProdForce*c[i];
Force[i][3]=Force[i][3]-DotProdForce*d[i];
In order to exit the loop, the `Force' maximum magnitude should be really low.
That means that you are subtracting really close numbers.

> Your loop (line 104-132) never updates ForceMagnitude.
It does update it. The loop actually ends in line 191

> Why should it need to? "ForceMagnitude" is defined later on as an equation of "Force" and Force updates in that loop so it should update just from that no?
There are no equations in c++, you've got assignments.
The code is executed from top to bottom

So please realize the importance of indentation.


> I try to avoid naming my variables the same as something in the std namespace.
Then you don't understand namespaces.
TheIdeasMan (6781)
Then you don't understand namespaces.


Well I do, it's just that I try to be a little imaginative with my variable names.
adrem7 (41)
I would blame `catastrophic cancellation'
Force[i][1]=Force[i][1]-DotProdForce*b[i];
Force[i][2]=Force[i][2]-DotProdForce*c[i];
Force[i][3]=Force[i][3]-DotProdForce*d[i];
In order to exit the loop, the `Force' maximum magnitude should be really low.
That means that you are subtracting really close numbers.


Could you please expand on this in more layman terms. Sorry to be slow but I am not as good as I should be at coding and don't really understand your point.


> Your loop (line 104-132) never updates ForceMagnitude.
It does update it. The loop actually ends in line 191

> Why should it need to? "ForceMagnitude" is defined later on as an equation of "Force" and Force updates in that loop so it should update just from that no?
There are no equations in c++, you've got assignments.
The code is executed from top to bottom


Sorry that's my maths brain speaking. I wanted to say that on line 178 the ForceMagnitude is a function of Force (mathematically speaking) but I didn't want to use the word function as I haven't used functions (coding wise) in my piece of code. Therefore as the Force updates for each time the loop runs, so must the ForceMagnitude. Like you say the loop ends at line 191. Therefore it does update and this is not the problem.


So please realize the importance of indentation.


Having been fiddling with this for a while now I know where everything is, but yes you are right, I should use indentation. Especially when asking for help like here. Sorry for my sloppy presentation!

I personally am stumped with the issue as you can probably tell. Why would comparing doubles only be an issue for some values of n? (Line 22 I define n). If comparing doubles was an issue, surely it would be a constant issue, but the code works for some n and not others. This is what really confuses me. (Note: I need the code to work for all numbers from 2 to 100 inclusive).

You might be right about the 'Catastrophic Cancellation' but I shall look into it as right now I don't really get why it would be an issue.

Thanks, A.

Edit:
If I wanted to compare my doubles say lines 181 to 188, how would I do this? Specifically for my problem I mean. I could take x=ForceMagnitude[i], y=ForceMagnitude[i+1]. and compare by seeing if abs(x-y)<epsilon for some epsilon, but my problem is now, surely my epsilon will have to change for different values of n at the start of the code... Am I right in thinking this (Note: going to give it a go and see where I get to).
Last edited on
ne555 (10692)
here's an example http://en.wikipedia.org/wiki/Loss_of_significance

My guess is that when you try to calculate the force, the computation is meaningless because of the error.
So you miss the target a little, and you continue from that bad position.
Keep cycling near the answer, but can't approach with such precision.

> Sorry for my sloppy presentation!
Don't be, just fix it.
adrem7 (41)
Is this better (just presentation wise I mean. I am still thinking about the error I get).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
 

#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 Forcemagnitude,ForceMagnitude[101],DotProdForce,Force[101][4],p,q[101],b[101],c[101],d[101],l[101],energy,Energy,Distance,E[10001],F[101][101][4],y[101][4],x[101][101][4],z[101][4],length[101],distance[101];

  clock_t t1,t2;
  t1=clock();

  /*  set the number of points */
  n=72;

  /* check that there are no more than 100 points */
  if(n>100){
    cout << n << " is too many points for me :-( \n";
    exit(0);
  }
  
  /* Loop for random points m times*/
  m=10;
  
  if (m>100){
    cout << "The m="<< m << " loop is inefficient...lessen m \n";
    exit(0);
  }
  
  /* reset the random number generator */
  srand((unsigned)time(0));  
  
  for (a=1;a<=m;a++){

    /* assign random points */
    for (i=1;i<=n;i++){
      x[a][i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[a][i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[a][i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      length[a]=sqrt(pow(x[a][i][1],2)+pow(x[a][i][2],2)+pow(x[a][i][3],2));
      
      for (k=1;k<=3;k++){
        x[a][i][k]=x[a][i][k]/length[a];
      }
    }

    /* calculate the energy */
    E[a]=0.0;
  
    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        distance[a]=sqrt(pow(x[a][i][1]-x[a][j][1],2)+pow(x[a][i][2]-x[a][j][2],2)+pow(x[a][i][3]-x[a][j][3],2));
        E[a]=E[a]+1.0/distance[a];
      }
    }
  }

  /* Find the lowest of these energies */

  for(a=1;a<m;a++){
    if (E[a]<E[a+1])
      for(i=1;i<=n;i++){
        E[a+1]=E[a],x[a+1][i][1]=x[a][i][1],x[a+1][i][2]=x[a][i][2],x[a+1][i][3]=x[a][i][3];
      }
    else
      for(i=1;i<=n;i++){
        E[a]=E[a+1],x[a][i][1]=x[a+1][i][1],x[a][i][2]=x[a+1][i][2],x[a][i][3]=x[a+1][i][3];
      }
  }
  
  Energy=E[a];

  /* Save the points to separate variable */

  for(i=1;i<=n;i++){
    y[i][1]=x[a][i][1];
    y[i][2]=x[a][i][2];
    y[i][3]=x[a][i][3];
  }
  
  /* Output these points */
  
  ofstream Bestrandompoints ("Bestrandompoints");
  Bestrandompoints << "E=" << E[a] << "\n";
  for(i=1;i<=n;i++){
    Bestrandompoints << y[i][1] << " " <<   y[i][2] << " " << y[i][3] << "\n";
  }
  Bestrandompoints.close(); 

  /* Start doing gradient flow approach */
  r=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(j=1;j<=n;j++){ 
      for(i=1;i<j;i++){
        
        Distance=sqrt(pow(y[j][1]-y[i][1],2)+pow(y[j][2]-y[i][2],2)+pow(y[j][3]-y[i][3],2));
        Force[j][1]=Force[j][1]+((y[j][1]-y[i][1])/(pow((Distance),3))); 
        Force[j][2]=Force[j][2]+((y[j][2]-y[i][2])/(pow((Distance),3))); 
        Force[j][3]=Force[j][3]+((y[j][3]-y[i][3])/(pow((Distance),3))); 
      } 

      for (i=j+1;i<=n;i++){ 

        Distance=sqrt(pow(y[j][1]-y[i][1],2)+pow(y[j][2]-y[i][2],2)+pow(y[j][3]-y[i][3],2));
        Force[j][1]=Force[j][1]+((y[j][1]-y[i][1])/(pow((Distance),3))); 
        Force[j][2]=Force[j][2]+((y[j][2]-y[i][2])/(pow((Distance),3))); 
        Force[j][3]=Force[j][3]+((y[j][3]-y[i][3])/(pow((Distance),3)));
      }
    }

    /* Add the force to my points */

    for(i=1;i<=n;i++){

      DotProdForce=Force[i][1]*y[i][1]+Force[i][2]*y[i][2]+Force[i][3]*y[i][3];
      
      b[i]=y[i][1];
      c[i]=y[i][2];
      d[i]=y[i][3];

      Force[i][1]=Force[i][1]-DotProdForce*b[i];
      Force[i][2]=Force[i][2]-DotProdForce*c[i];
      Force[i][3]=Force[i][3]-DotProdForce*d[i];

      y[i][1] = b[i]+p*(Force[i][1]);
      y[i][2] = c[i]+p*(Force[i][2]);
      y[i][3] = d[i]+p*(Force[i][3]);

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

      for (j=1;j<=3;j++){
        y[i][j]=y[i][j]/length[i];
      }
    }

    /* Calculate the new energy */
    
    E[r]=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));
        E[r]=E[r]+1.0/Distance;
      }
    }

    energy=E[r];

    /* Choose best energy and therefore best points */
    if (energy<Energy)
      Energy=energy,s=s+1;
    else
      energy=Energy,f=f+1;
  
    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=" << Forcemagnitude << "\n";
    cout << "r=" << r << "\n";

    r=r+1;
  }

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

  ofstream mypoints ("mypoints");
  for(i=1;i<=n;i++){
    mypoints << y[i][1] << " " <<   y[i][2] << " " << y[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;

}


Also lets say that you have identified the error there. How do I go about fixing this? (I am really new to C++ so please bare with me. I appreciate I am very slow!)

Edit: Is there a way to get the code on the next line even if I want it to read as if it were on the same one? Like on line 17. Just to make the code less wide. Can you just hit return and put it on the next line? Will it read in the same way?

Thanks, A.
Last edited on
ne555 (10692)
I don't know how to fix it.
You may try to increase the bits in the representation (long double),
change the computation, in order to not do that subtraction (dunno how)
assume that you can't have so much precision.

> Like on line 17. Just to make the code less wide.
Yes, whitespace is free.
1
2
3
4
5
6
7
8
9
 
	double Forcemagnitude, ForceMagnitude[101], DotProdForce, Force[101][4],
		p,
		q[101], b[101], c[101], d[101], l[101],
		energy, Energy,
		Distance,
		E[10001],
		F[101][101][4],
		y[101][4], x[101][101][4], z[101][4],
		length[101], distance[101];


Still it's quite obfuscated (use meaningful names, limit scope, use functions)
1
2
3
4
5
6
7
8
9
10
11
 
	/* Find the lowest of these energies */
	for(a=1;a<m;a++){
		if (E[a]<E[a+1])
			for(i=1;i<=n;i++){
				E[a+1]=E[a],x[a+1][i][1]=x[a][i][1],x[a+1][i][2]=x[a][i][2],x[a+1][i][3]=x[a][i][3];
			}
		else
			for(i=1;i<=n;i++){
				E[a]=E[a+1],x[a][i][1]=x[a+1][i][1],x[a][i][2]=x[a+1][i][2],x[a][i][3]=x[a+1][i][3];
			}
	}
¿what's that supposed to do? Because find the lowest require just one loop.
adrem7 (41)
¿what's that supposed to do? Because find the lowest require just one loop.


It is supposed to find the lowest energy. Here energy is denoted E[a] and I want to find the lowest of every E[a]. I do not know how to do this in one loop and so did it like this. This is quite inefficient I know as essentially you look at E[1] and E[2], see which one is lowest, E[2] and E[3], see which one is lowest etc. etc. unitl you have the lowest one overall. If this can be done in one loop please let me know how.

Thanks, A.

Edit: Excuse last edit, I was just putting loops inside loops and getting crap load of stuff.
Last edited on
ne555 (10692)
adrem7 (41)
Thanks for the link. To be honest lines 31-97 could be replaced with something like:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
 
  /* 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;
  }


It is longer, but less messy I think. Obviously the minimum link you just gave me will be better still so I shall try to get my head round that!
(Honestly I didn't know that existed so I just came up with something to do the job from my limited knowledge).

A.
adrem7 (41)
Going back to your original point of why my code might stop working. You said it could be Catastrophic Cancellation.

Say you are right (for arguments sake). Why would this make my code run indefinitely.

1
2
3
 
Force[i][1]=Force[i][1]-DotProdForce*b[i];
Force[i][2]=Force[i][2]-DotProdForce*c[i];
Force[i][3]=Force[i][3]-DotProdForce*d[i];


Surely if they cancelled I would just get a Force of 0.0000 at some point. What is wrong with that?

There are two other places I use Force after this before it is set back to 0.0 at line 110. These are:

Line 149:
1
2
3
 
y[i][1] = b[i]+p*(Force[i][1]);
      y[i][2] = c[i]+p*(Force[i][2]);
      y[i][3] = d[i]+p*(Force[i][3]);


Line 182:
ForceMagnitude[i]=pow((pow(Force[i][1],2)+pow(Force[i][2],2)+pow(Force[i][3],2)),0.5);

If either of these were taking the Force value to be 0.0000 then I don't see an issue. The first would just give:
y[i][1] = b[i]+0.00;

or c[i] or d[i] respectively, and the second would just give:
ForceMagnitude[i]=pow((pow(0.00,2)+pow(Force[i][2],2)+pow(Force[i][3],2)),0.5);

Which is just:
ForceMagnitude[i]=pow((0.00+pow(Force[i][2],2)+pow(Force[i][3],2)),0.5);

Which is just:
ForceMagnitude[i]=pow((pow(Force[i][2],2)+pow(Force[i][3],2)),0.5);

(or interchange this 0 for the second or third term respectively).
Neither of which I would consider an issue. I am still trying to see if this is the issue, but it takes a while as I have to check so many numbers.

Thanks for all the help so far. I really appreciate everything you have done!

A.
adrem7 (41)
I have fixed the issue by looking at it in a more mathematical manner. The changes have been made and the code is as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
 
#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=73;

  /* 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;
  Forcemagnitude=1.0;

  if(n>80)
    p=0.005;
  else if(n>60)
    p=0.01;
  else if(n>40)
    p=0.02;
  else
    p=0.05;

  cout << "p=" << p << "\n";

  while(Forcemagnitude>0.00005){

    /* 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 << 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,p=(9.9*p)/10;
 
    cout << "p=" << p << "\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 << "FM=" << Forcemagnitude << "\n";
    
    a=a+1;
  }

  /* Output new energy */
  cout << "Successes=" << s << " Failures=" << f <<"\n";
  cout <<fixed;
  cout << setprecision (20) << "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;

}


There are issues that need to be solved for n>85ish, but these should be able to be done by fiddling with values of p on lines 134 and 244, as well as changing required accuracy on line 145.

Thanks for all the posts and help. I appreciate that the actual coding aspect has a way to go really (functions, not comparing floating point numbers (use of an epsilon)) and I shall implement these in due time, asking for help on new thread with more specific questions if I need it.

Thanks, A.
Last edited on
Topic archived. No new replies allowed.