Previous ToC Up Next

8. Quartile Checks

8.1. Quartiles

Bob: Hi Alice, what's next? I can think of a few further improvements to my Plummer code, but do you have something specific, you'd like to add?

Alice: I'm happy that we have done this energy checking, which seems to indicate that your code probably does produce Plummer's model correctly. However, I'd like to be even more sure, since it would really be terrible if we start our future simulations with the wrong initial conditions. That would make debugging our evolution codes quite confusing, to say the least.

Bob: There are many ways to check things. The question is, what would be the simplest way. Ah, how about just looking at mass shells, shells with a certain mass fraction, and check that they have the correct radius? After all, we computed the cumulative mass as a function of the radius already.

Alice: That is a good idea. For one thing, the half-mass radius should come out with the correct value. And for good measure, we may as well compute the quarter-mass radius and the three-quarter-mass radius. In other words: the radii of the three quartiles, within which 1/4, 1/2, and 3/4 of the total mass are enclosed.

Bob: That we should be able to do from scratch. We can start again with minimal versions of the Body and Nbody classes, just enough to read a snapshot in, and initialize everything properly. The only thing to do is add a method quartiles that prints out the three quartiles.

Fortunately, Ruby has a handy sort method, that comes with the Array class. As usual, it has an associated sort! method that affects the array itself, unlike sort that returns a new array. Let's see. We have to sort the particles in radial order, using , the distance from the center. It will be easier to use the square instead, . This can be done in three lines: 

   def order_squared_radii
     a = []
     @body.each{|b| a.push b.pos*b.pos}
     a.sort!
   end

The first line defines a to be an empty array. The second line fills the array with the squares of the radial positions of all particles. Then the third line sorts that array. Couldn't be simpler!

8.2. Coding

Alice: Can I write the quartiles method? I'd like to get some more experience with Ruby.

Bob: Here is the keyboard!

Alice: Since you have defined order_squared_radii as a member function, we may as well define quartiles as a member function too. We can then immediately invoke order_squared_radii to create our order list of particle positions. All we have to do then is to cut the list in four equal parts, and document the places where we cut the list. Hey, that is too simple to learn much about Ruby! This should do the job: 

   def quartiles
     a = order_squared_radii
     n = a.size
     r_1 = a[(n/4.0).round - 1]
     r_2 = a[(n/2.0).round - 1]
     r_3 = a[(n*3/4.0).round - 1]
     print "The values of the three quartiles for r(M) are:\n"
     print "  r(1/4) = "
     printf("%.4g\n", r_1)
     print "  r(1/2) = "
     printf("%.4g\n", r_2)
     print "  r(3/4) = "
     printf("%.4g\n", r_3)
   end

Bob: I bet it does. I saw you looking up round: indeed, that is a predefined method that rounds a floating point number to an integer. In case of a number exactly in between, such as 2.5, it rounds up, to 3.0. You could have used integer division directly, instead of floating point division, but I agree that what you wrote might be easier to understand.

Alice: We could quibble about exactly where to choose the boundaries for the quartiles. For example, for a four-particle system, things come out just right, but for a five particle system, the half mass radius encloses 3/5 of the mass. We could have take the average of the two-particle and three-particle distance, but that would have been overkill, in my opinion. Anyway, the fluctuations in the positions of the particles, in an N-body system will be of order , so we don't have to worry about enclosing one particle more or less, at least for .

8.3. Code

Bob: Let me print out the whole program, and write it out to a file called quartiles1.rb: 

 require "acs"
 
 class Body
 
   attr_accessor :mass, :pos, :vel
 
   def initialize(mass = 0, pos = Vector[0,0,0], vel = Vector[0,0,0])
     @mass, @pos, @vel = mass, pos, vel
   end
 
 end
 
 class NBody
 
   attr_accessor :time, :body
 
   def initialize
     @body = []
   end
 
   def order_squared_radii
     a = []
     @body.each{|b| a.push b.pos*b.pos}
     a.sort!
   end    
 
   def quartiles
     a = order_squared_radii
     n = a.size
     r_1 = a[(n/4.0).round - 1]
     r_2 = a[(n/2.0).round - 1]
     r_3 = a[(n*3/4.0).round - 1]
     print "The values of the three quartiles for r(M) are:\n"
     print "  r(1/4) = "
     printf("%.4g\n", r_1)
     print "  r(1/2) = "
     printf("%.4g\n", r_2)
     print "  r(3/4) = "
     printf("%.4g\n", r_3)
   end
 
 end
 
 include Math
 
 nb = ACS_IO.acs_read(NBody)
 nb.quartiles

8.4. Testing

Alice: Let's try it out: 

 |gravity> kali mkplummer2.rb -n 100 | kali quartiles1.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 100
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 157807557255272290078302474298427494808
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.642
   r(1/2) = 1.737
   r(3/4) = 4.359

 |gravity> kali mkplummer2.rb -n 100 | kali quartiles1.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 100
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 251979634788701045369825442993623037445
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.5245
   r(1/2) = 1.511
   r(3/4) = 3.524

 |gravity> kali mkplummer2.rb -n 100 | kali quartiles1.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 100
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 116741143007870386439078079644181238717
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.5138
   r(1/2) = 1.192
   r(3/4) = 3.349
Quite noisy. Not too surprising: for the inner quartile we are dealing with only 25 particles, so we should expect noise. Indeed, there are fluctuations of about twenty percent in the first quartile.

Bob: Let's try 10,000 particles instead: 

 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles1.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 195924807145294240096344655914317919224
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.6504
   r(1/2) = 1.665
   r(3/4) = 4.717

 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles1.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 42611728057165554375136236111523795370
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.6226
   r(1/2) = 1.655
   r(3/4) = 4.727

 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles1.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 323527003508540608555251954323353747057
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.6753
   r(1/2) = 1.724
   r(3/4) = 4.761
Not bad! Fluctuations of a few percent, at most, as it should be. Encouraging.

8.5. Checking the Math

Alice: Now let's check the numbers. We have used already a few times eq. (26):

Or in our units, with :

While we walked through your code, we reproduced the inverted form, expressing , in eq. (20) as:

Bob: I have one problem with this kind of testing. First we write a code, using the above equation. Then we test the code, using the same equation. This will catch implementation bugs alright, but it won't catch any mistake we might have made in the above equation. In this case the equation is simple, but as a general procedure, I'm not sure how much this really buys us.

Alice: Excellent point! And there is an easy check: Let us insert eq. (96) into the right-hand side of eq. (95), to check whether we get the original amount of mass back. Let us call the left-hand side of eq. (95) , so that we can distinguish it from the mass we start with, in eq. (96). Eq. (95) then reads:

This was an almost trivial exercise. But I'm glad we checked.

Bob: Yes, many bugs turn out to be `almost trivial' as well, and yet they can waste many hours of your time while you're chasing them. By the way, wouldn't it have been much faster to substitute eq. (96) into the right-hand side of eq. (94)? Calling the used in eq. (94) , we would have gotten:

Alice: You're right, that is faster. And if we were to write a text book, we would certainly hide the fact that we did it first the hard way around. But in practice, who cares? The important thing is to check that things are correct. And they are, we now know.

8.6. Checking the Code

Bob: Let's compute the expected quartile radii, from eq. (96):

Alice: Those numbers look completely different from what we got from our code! And the strange thing is, the first number comes out too small, while the next ones are too large, especially the last one. So it is not a matter of having overlooked a scale factor, or something like that.

Bob: What a pity! I was more or less convinced that my Plummer code was correct. Obviously there is something seriously wrong with the way particles are sprinkled into space.

Alice: What is really strange is the fact that we got the energy correct. First of all, the virial theorem did hold, something that is typically violated if you make a random mistake in the energies, since it is unlikely that you will make the same mistake in the kinetic as in the potential energy. But secondly, we actually computed the energy, and we found that your code produced it correctly.

Bob: Yes, that is puzzling. Unless . . .

Alice: . . . unless my quartile code is wrong. But that is also very unlikely. How much can there be wrong in just a few lines of code?

Bob: It wouldn't be the first time, to find a nasty bug in a small code. I guess that's why they are called bugs! Bugs are small, and can hide even in a code of a single line. Let's have a look.

Alice: Yes, the good thing about a small code is that it makes sense to go through each line, one by one. And there are less than ten lines here where actual calculations are done.

8.7. Checking the Code-Checking Code

Bob: Let's be systematic. The driver part calls quartiles, and that method in turn starts by calling order_squared_radii. You like long names, don't you? You could have called it place_an_order_to_order_squared_radii.

Alice: Yes, I do like long names, since that way I will remember later what was doing what. I hate that old style in which it probably would have been called srtrsq for `sorting r squared'. Forty years ago there was a reason, when single memory locations were expensive, but that excuse has long since gone.

Bob: Well, I can't for the life of me see anything wrong with order_squared_radii. Back to quartiles. In the second line you find the number of particles. Nothing wrong with that either. In the third line you cut the list at the place of the first quarter. Indeed, in that way you are bound to find the mass shell that contains one quarter of the mass, and precisely the inner quarter.

Alice: Iron logic, and I, too, cannot see anything wrong with that.

Bob: So we both agree that a[(n/4.0).round - 1] returns what you have ordered: the squared radius of the inner quarter of . . .

Alice: . . . the squared radius!!

Bob: Ah, yes, that's the problem! You had forgotten that you ordered the squares of the radial distances, for convenience. Instead, you just assigned the values to the quartile radii themselves.

Alice: That must have been the problem. And that immediately explains why the errors grew bigger for larger values of r. Okay, that is easy to fix. Let us rewrite these three lines, and put the result in a file called quartiles.rb: 

     r_1 = sqrt(a[(n/4.0).round - 1])                                         
     r_2 = sqrt(a[(n/2.0).round - 1])                                         
     r_3 = sqrt(a[(n*3/4.0).round - 1])

And let's test it right away: 

 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 201299017348521780587188254399647312211
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.8137
   r(1/2) = 1.292
   r(3/4) = 2.176
 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 42982039013636173889503342183695809112
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.8216
   r(1/2) = 1.315
   r(3/4) = 2.181
 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 198688716986344613646322407980102362651
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.8141
   r(1/2) = 1.307
   r(3/4) = 2.186
Bob: Wonderful! So my Plummer code was correct after all. But I'm glad we checked so thoroughly.

Alice: You may laugh, and perhaps I'm getting to much addicted to testing, but let me do one final little check. Having just made a really silly mistake determining the quartile radii, I feel I cannot take anything for granted. Let me see whether the quartile radii are really quartile radii. Starting with:

We have

Similarly with the half-mass radius , and:

Bob: Now that is what I would call overkill. But you're right, it doesn't hurt, and better safe than sorry.

Alice: Indeed!
Previous ToC Up Next