Previous ToC Up Next

2. Distance Measures

2.1. Phase Space

Bob: Now that we have a shared time step code nbody_sh1.rb, I'd like to make sure that it really gives the same results as the previous code nbody_cst1.rb that used constant time steps. Having good energy accuracy is one thing, but it would be much better to check the actual particle positions.

Alice: How about computing the distance between the results of two runs, in 6N-dimensional phase space? A whole N-body system can be viewed as a single point in a space that spans all the degrees of freedom of the positions and velocities of all particles.

More generally, if we model a system in D dimensions, a single particle will have D components for its position as well as its velocity, or 2D in total. N particles will have 2DN such degrees of freedom. For planar N-body systems, we need a 4N-dimensional phase space, and for N-body systems in the real three-dimensional world we are dealing with a 6N-dimensional phase space.

Bob: In practice this means, I presume, that you use Pythagoras to compute the distance in such a large space. In other words you take the square root of the sum of the squares of all the differences, right?

Alice: Indeed.

Bob: So another way of looking at your suggestion is that you compute the average distance between corresponding particles in the usual 6-dimensional phase space. Apart from a normalization factor N, the result would be the same.

Alice: Yes, that's right. I like the more abstract picture, but your description boils down to the same thing.

Bob: I find it a lot easier to think about 2D dimensions than about 2ND dimensions. In either case, the practical procedure is that we have to subtract two N-body systems, point by point, and then we have to compute the norm or `size' of the resulting system, a type of absolute value.

Alice: Ah, so you do like abstract thinking, after all! How about defining a minus operator, that allows us to subtract two N-body systems, nb1 and nb2, just by writing nb1-nb2? That would give us a nicely compact notation. And the 2DN-dimensional distance between the two systems could then be given by (nb1-nb2).abs if we define the correct absolute value operator abs for the NBody class.

Bob: Good idea! However, there is no guarantee that two N-body systems will show their particles in the output in the exact same order. In other words, what we need is a way to identify which particle in the one system corresponds to which particle in the other system.

Alice: I guess we have to number them, by giving each body a unique identifier.

2.2. Numbering Bodies

Bob: Given that every object in Ruby already has an identifier called object_id, a natural name for such an ID would be body_id. Now that is easy to implement. How about writing a nifty tool like this, in file nbody_set_id.rb:

 #!/usr/local/bin/ruby -w
 
 require "nbody.rb"
 
 options_text= <<-END
 
   Description: Takes an N-body system, and gives each body a unique ID
   Long description:
     This program accepts an N-body system, and assigns a number to each
     body consecutively, as an instance variable @body_id which takes on
     integer values.
 
     (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org
 
     example:
                kali #{$0} -n 0
 
 
   Short name:          -n
   Long name:            --starting_number
   Value type:           int
   Default value:        1
   Description:          value of @body_id for 1st body
   Variable name:        n
   Long description:
     This option allows the user to start numbering the particles at a
     value different from the default value (0, say, or 10, or whatever).
 
 
   END
 
 c = parse_command_line(options_text)
 
 class Body
   attr_accessor :body_id
 end
 
 nb = ACS_IO.acs_read(NBody)
 i = c.n
 nb.body.each do |b|
   b.body_id = i
   i += 1
 end
 nb.acs_write($stdout, false, c.precision)

2.3. A Test Count

You see, I even allowed for taste: the default is to start numbering at one, but you can start at zero as well, or at any other place you like. Let's test it with our Plummer Model:

 |gravity> kali nbody_set_id.rb -h
 Takes an N-body system, and gives each body a unique ID
 -n  --starting_number  : value of @body_id for 1st body [default: 1]
 --verbosity            : Screen Output Verbosity Level  [default: 1]
 --acs_verbosity                : ACS Output Verbosity Level     [default: 1]
 --precision            : Floating point precision       [default: 16]
 --indentation          : Incremental indentation        [default: 2]
 -h  --help             : Help facility
 ---help                        : Program description (the header part of --help)
 |gravity> kali mkplummer.rb -n 2 | kali nbody_set_id.rb --precision 5
 ==> Plummer's Model Builder <==
 Number of particles            : N = 2
 pseudorandom number seed given : 0
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 ==> Takes an N-body system, and gives each body a unique ID <==
 value of @body_id for 1st body : n = 1
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 5
 Incremental indentation                : add_indent = 2
              actual seed used  : 252462866881473732163958989257426001914
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           1
         Float mass
              5.00000e-01
         Vector pos
              1.94755e-01   9.73956e-02  -1.22820e-01
         Vector vel
              1.78002e-01   5.49335e-01  -4.08101e-01
       Body body[1]
         Fixnum body_id
           2
         Float mass
              5.00000e-01
         Vector pos
             -1.94755e-01  -9.73956e-02   1.22820e-01
         Vector vel
             -1.78002e-01  -5.49335e-01   4.08101e-01
     Float time
          0.00000e+00
 SCA
 |gravity> kali mkplummer.rb -n 2 | kali nbody_set_id.rb --precision 5 -n 42
 ==> Plummer's Model Builder <==
 Number of particles            : N = 2
 pseudorandom number seed given : 0
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 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  : 210733041374937108452252448370138453664
 ==> Takes an N-body system, and gives each body a unique ID <==
 value of @body_id for 1st body : n = 42
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 5
 Incremental indentation                : add_indent = 2
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           42
         Float mass
              5.00000e-01
         Vector pos
              1.16550e-01  -4.85085e-02   2.15785e-01
         Vector vel
             -5.72883e-02   3.45794e-01   6.14121e-01
       Body body[1]
         Fixnum body_id
           43
         Float mass
              5.00000e-01
         Vector pos
             -1.16550e-01   4.85085e-02  -2.15785e-01
         Vector vel
              5.72883e-02  -3.45794e-01  -6.14121e-01
     Float time
          0.00000e+00
 SCA

2.4. Defensive Programming

Bob: Now, with identifiers in place, we can extend our nbody.rb file by adding a subtraction operator:

   def -(other)
     if other.class != self.class
       raise "other.class.name = #{other.class.name} != #{self.class.name}"
     end
     if (n = other.body.size) != body.size
       raise "other.body.size = #{other.body.size} != #{body.size}"
     end
     nb = NBody.new
     body.each_index do |i|
       if (id = body[i].body_id) == nil
         raise "body[#{i}].body_id == nil"
       end
       ob = other.body.find{|oi| oi.body_id == id}
       if ob == nil
         raise "body_id = #{body[i].body_id} not found in other N-body system"
       end
       nb.body[i] = Body.new
       nb.body[i].pos = body[i].pos - ob.pos
       nb.body[i].vel = body[i].vel - ob.vel
     end
     nb
   end

This includes proper testing of particle identities.

Alice: And I see that you are returning the newly constructed nbody system nb at the end of the method. Also, you're getting quite defensive in your programming!

Bob: This is a situation where defensive programming is called for. In general, I don't like to litter my code with many error checking lines all over the place, but in this particular case, it is quite likely that the difference operator will be invoked for incompatible N-body systems. So I'm checking whether the other object is really of NBody class, whether it has the same number of particles as the calling NBody object, and most importantly, I am then checking whether each particle in the calling system has a counter part in the other system, with the same body_id.

Alice: Applause, applause! As you know, I like to check things, as a matter of habit, in order to prevent future surprises that may be hard to debug. But I agree, you can overdo with anything, and there is a lot to say for keeping code small. In this case, though, the whole code is so small that it doesn't distract to have these extra safeguards built in.

2.5. Putting it together

Bob: With our new tools in hand, it now becomes really trivial to check for the actual size of the difference between two N-body systems, by defining the absolute value method, also within the file nbody.rb:

   def abs
     a = 0
     body.each{ |b| a += b.pos*b.pos + b.vel*b.vel }
     sqrt a
   end

Alice: Sometimes we may want to look only at the positions, restricting ourselves to configuration space, rather than to phase space.

Bob: That's easy. Here is a position-only version of abs:

   def abs_pos
     a = 0
     body.each{ |b| a += b.pos*b.pos }
     sqrt a
   end

Alice: Let me try to implement the module that gives the difference between two N-body systems, as the 6N-dimensional distance in phase space, or the 2DN-dimensional distance for the general case of D spatial dimensions. It's time to practice my Ruby skills again. Let's call the file nbody_diff.rb:

 #!/usr/local/bin/ruby -w
 
 require "nbody.rb"
 
 options_text= <<-END
 
   Description: 6N-dimensional phase space distance between two N-body systems
   Long description:
     This program accepts two N-body systems, and computes the distance between
     them in the 2DN-dimensional phase space of all positions and velocities,
     if they are given in D spatial dimensions.  The values of the masses of
     the particles are not used, and each particle is given equal weight.  In
     order to compare corresponding particles in the two systems, the body_id
     instance variable of each body is used as an identifier.
 
     (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org
 
     example:
     cat run1.out run2.out | kali #{$0}
 
 
   Short name:          -r
   Long name:            --positions_only
   Value type:           bool
   Description:          Using only positions, not velocities
   Variable name:        position_only_flag
   Long description:
     This option allows the user to compute the distance between two N-body
     systems using only particle information; velocity information is not used.
 
 
   END
 
 c = parse_command_line(options_text)
 
 class Body
   attr_accessor :body_id
 end
 
 nb1 = ACS_IO.acs_read(NBody)
 nb2 = ACS_IO.acs_read(NBody)
 nb = nb1 - nb2
 d = nb.body[0].pos.size
 n = nb.body.size
 print "#{2*d}N-dim. phase space dist. for two #{n}-body systems:  "
 if c.position_only_flag == false
   printf(" %.#{c.precision}e\n", nb.abs)
 else
   printf(" %.#{c.precision}e\n", nb.abs_pos)
 end

How about this?

2.6. Checking it out

Bob: Looks fine, as far as I can see. You put in options to allow for more precision, without burdening the reader with all those extra digits in the general case. And you also gave an option for positions-only comparisons.

Time to test it! We still have the initial conditions lying around for a circular binary, I believe. Ah yes, here it is, in circular_binary.in:

 ACS
   NBody 
     Array body
       Body body[0]
         Float mass
              4
         Vector pos
              1  0
         Vector vel
              0  1
       Body body[1]
         Float mass
              4
         Vector pos
              -1  0
         Vector vel
              0  -1
 SCA
To start with, let's see whether two identical systems are assigned a distance of zero correctly:

 |gravity> kali nbody_set_id.rb < circular_binary.in > tmp0.in
 ==> Takes an N-body system, and gives each body a unique ID <==
 value of @body_id for 1st body : n = 1
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 |gravity> cat tmp0.in tmp0.in | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 4N-dim. phase space dist. for two 2-body systems:   0.0000000000000000e+00
Alice: Good! Let's perturb our binary slightly. If we use the most popular Pythagorean triangle, we can immediately check whether our Cartesian distances in phase space come out correctly. Here is a file pert_circ_binary.in:

 ACS
   NBody 
     Array body
       Body body[0]
         Float mass
              4
         Vector pos
              1  0.03
         Vector vel
              0  1
       Body body[1]
         Float mass
              4
         Vector pos
              -1  0
         Vector vel
              0  -1.04
 SCA
And here is our first real result:

 |gravity> kali nbody_set_id.rb < pert_circ_binary.in > tmp1.in
 ==> Takes an N-body system, and gives each body a unique ID <==
 value of @body_id for 1st body : n = 1
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 |gravity> cat tmp0.in tmp1.in | kali nbody_diff.rb
 ==> 6N-dimensional phase space distance between two N-body systems <==
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 4N-dim. phase space dist. for two 2-body systems:   5.0000000000000031e-02
Just what it should be! Now for positions only, and a higher precision:

 |gravity> cat tmp0.in tmp1.in | kali nbody_diff.rb -r --precision 5
 ==> 6N-dimensional phase space distance between two N-body systems <==
 Using only positions, not velocities
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 5
 Incremental indentation                : add_indent = 2
 4N-dim. phase space dist. for two 2-body systems:   3.00000e-02
Perfect. I think we're ready to move on.
Previous ToC Up Next