Previous ToC Up Next

15. A Complete Vector Class

15.1. Vector Multiplication

Carol: Which means it is time to move on to multiplication. But here we have another problem: there is multiplication and then there is multiplication.

Dan: You mean?

Carol: We can multiply a vector with a scalar number; for a vector 

  v = [1, 2, 3]
we would like to see multiplication by two giving us: 

  2 * v = [2, 4, 6]
But in addition (no pun intended) we want to form an inner product of two vectors. In particular, we would like to get the square of the length of a vector by forming the inner product with itself: 

  v * v = 2*2 + 4*4 + 6*6 = 56
Or more generally, for 

  w = [10, 10, 10]
we want to be able to take the inner product of v and w to give us: 

  v * w = 2*10 + 4*10 + 6*10 = 120
We could of course define different method names for these two operations, like multiply_scalar and inner_product, but something tells me that we will be happier using * for both.

Dan: Certainly I will be happier that way!

Carol: Well, how about this, in vector_try.rb: 

   def *(a)
     if a.class == Vector              # inner product
       product = 0
       self.each_index{|k| product += self[k]*a[k]}
     else
       product = Vector.new            # scalar product
       self.each_index{|k| product[k] = self[k]*a}
     end
     product
   end

Time for a workout: 

 |gravity> irb
 require "vector_try.rb"
 true
 v1 = Vector[1, 2, 3]
 [1, 2, 3]
 v2 = Vector[5, 6, 7]
 [5, 6, 7]
 v1 * 3
 [3, 6, 9]
 v1 * v1
 14
 v1 * v2
 38
 v1 * v2 * 3
 114
 v1 * 3 * v2
 114
 quit

15.2. An Unnatural Asymmetry

Dan: So far so good, but why do you put the number 3 only at the end and in the middle, why not in front? That would seem more natural! Let me try: 

 |gravity> irb
 require "vector_try.rb"
 true
 v1 = Vector[1, 2, 3]
 [1, 2, 3]
 v1 * 3
 [3, 6, 9]
 3 * v1
 TypeError: Vector can't be coerced into Fixnum
        from (irb):4:in `*'
        from (irb):4
        from :0
 quit
Hmmm, I guess not.

Carol: Yes, it would be more natural, but it doesn't work. Do you see why? Remember, when we write v2 * 3 we invoke the * method of the vector object v2. I understand that you would like to write 3 * v2, but in that case you have a problem: you would be trying to invoke the * method of the number object 3 . . . .

Erica: That's all and well, as a formal explanation, but if you can only write [1, 2, 3] * 3 and never 3 * [1, 2, 3], I'm not very happy with it. The whole point was to make our life easier, and to make the software notation more natural, more close to what you would write with pen and paper . . .

Dan: I agree. Half a solution is often worse than no solution at all. Perhaps we should just return to our earlier component notation.

Carol: I can't argue with that. But I'm not yet willing to give up. I wonder whether there is not some sort of way out. Let me see whether I can find something.

15.3. Vector Division

Erica: For now at least, let's finish the job we started . . .

Carol: . . . which means that we have to add a division method as well. In the case of division, we only have to deal with scalar division.

Erica: Ah, yes, of course, you have an inner product, but not an inner quotient.

Carol: Well, not completely `of course', there is something called `geometric algebra'. If you're curious, you can search for it on the internet.

Dan: I'm not curious, let's move on.

Carol: Okay, so I will punish attempts to divide two vectors by letting an error message appear. In Ruby you can use the command raise to halt execution and display an error message: 

   def /(a)
     if a.class == Vector
       raise
     else
       quotient = Vector.new           # scalar quotient
       self.each_index{|k| quotient[k] = self[k]/a}
     end
     quotient
   end

A quick check: 

 |gravity> irb
 require "vector_try.rb"
 true
 v1 = Vector[1, 2, 3]
 [1, 2, 3]
 v2 = Vector[5, 6, 7]
 [5, 6, 7]
 v1 / 3.0
 [0.333333333333333, 0.666666666666667, 1.0]
 v2 / 0.5
 [10.0, 12.0, 14.0]
 v1 / v2
 RuntimeError: 
        from ./vector_try.rb:41:in `/'
        from (irb):6
        from :0
 quit
You see, we got a run time error, because the raise statement had the effect of raising an error, as it is called in Ruby.

15.4. The Score: Six to One

Erica: Let me see whether I can sum up what we've learned. We've successfully implemented seven legal operations, unary/binary addition/multiplication, scalar/vector multiplication and scalar division.

Of these seven, six are okay. It is only with scalar multiplication that we encounter a problem, namely a lack of symmetry.

Dan: Can you define "okay" for the other six?

Erica: I will be explicit. I will use a notation where v1 and v2 are vectors, and s is a scalar number, which could be either an integer or a floating point number. I will call something "well defined" if it gives a specific answer, and not an error message.

Here is the whole list of the six operations that I am now considering "okay":

1) unary plus: +v1 is okay, because it is well defined.

2) unary minus: -v1 is okay, because it is well defined.

3) binary plus: v1 + v2 is okay, because both v1 + v2 and v2 + v1 are well defined, and both give the exact same answer, as they should.

4) binary minus: v1 - v2 is okay, because both v1 - v2 and v2 - v1 are well defined, and both give the exact opposite answer, as they should.

5) vector times: v1 * v2 is okay, because both v1 * v2 and v2 * v1 are well defined, and both give the exact same answer, as they should.

6) scalar slash: v1 / s is okay, because it is well defined. The alternative, s / v1 is not defined, and will give an error message. However, that's perfectly okay too, because we have decided that we don't allow division by vectors.

So far, the score is six to zero, and we seem to be winning. The problem is that we are losing in the case of number 7:

7) scalar times: v1 * s is okay, because it is well defined. The alternative, s * v1 should give the same answer too, but unfortunately, in our implementation it is not defined, and in fact, as we have seen, it gives an error message.

Dan: Thanks, Erica, that's a very clear summary. So all we have to do, to save the day, is to repair s * v1, so that it gives a well defined result, and in fact the same result as v1 * s.

15.5. A Solution

Carol: Ah, of course! How stupid, we should have thought about it right away!

Dan: About what?

Carol: Class augmentation! We can just augment the number classes, telling them how to multiply with vectors! Just as we have already augmented the array class, telling it how to convert an array into a vector.

Dan: I don't see that. What exactly are you trying to repair?

Carol: Let's go back to square one. We wanted to make our scalar-vector multiplication symmetric. The problem was, when you have a vector v and you want to write

3 * v

you are effectively asking the number 3 to provide a method that knows how to handle multiplication of 3 with a vector. But the poor number 3 has never heard of vectors! No wonder we got an error message. Let me show again explicitly what the error message says: 

 |gravity> irb
 require "vector_try.rb"
 true
 v = Vector[1, 2, 3]
 [1, 2, 3]
 v * 3
 [3, 6, 9]
 3 * v
 TypeError: Vector can't be coerced into Fixnum
        from (irb):4:in `*'
        from (irb):4
        from :0
 quit
You see: it tells us that the number 3 expects to make a multiplication with something like another number. It tries to convert a vector into such an other number, but doesn't succeed.

Erica: So, what we really would like to do is to augment the rules for the class Fixnum, which is the class of integers, to explain how to multiply with an object of class Vector. We can then do the same for the class Float, the class of floating point numbers.

15.6. Augmenting the Fixnum Class

Carol: Indeed. Well, there is no stopping us now to add to the behavior of the Fixnum class! Here we go, in vector.rb: 

 class Fixnum
   alias :original_mult :*
   def *(a)
     if a.class == Vector
       a*self
     else
       original_mult(a)
     end
   end
 end

Dan: Whoa! That's quite a bit more complicated than I had expected. You must have been reading Ruby manuals lately!

Carol: I admit, I did. Here is what happens. If we start with the number 3, and take a vector v, and write 

  3 * v
then the * method of the fixed point number 3 first checks whether v is a Vector. In our case, it is, so then it returns the result 

  v * self
which in the case of 3 is simply 

  v * 3
and that is something we have already defined, in the Vector class.

If, on the other hand we write anything else, such as 

  3 * 8.5
then the * method of the fixed point number 3 finds that 8.5 is not a Vector, so it applies the original definition of multiplication, to the number 8.5, as it should.

Dan: So the alias notion means that whatever the original * did is now assigned to original_mult instead?

Carol: Exactly. Writing alias :x :y means that x becomes an alias for y. Or so the theory goes. Let's see what the practice tells us: 

 |gravity> irb
 require "vector.rb"
 true
 v = [1, 2, 3].to_v
 [1, 2, 3]
 v * 3
 [3, 6, 9]
 (v * 3).class
 Vector
 3 * v
 [3, 6, 9]
 (3 * v).class
 Vector
 quit

15.7. Augmenting the Float Class

Dan: Not bad! Great! 3 * v now produces the exact same thing as v * 3. Congratulations! Another hope fulfilled. And I guess you should do the same thing for floating point numbers, right?

Carol: Right. All very similar: 

 class Float
   alias :original_mult :*
   def *(a)
     if a.class == Vector
       a*self
     else
       original_mult(a)
     end
   end
 end

This time my hope for success will be quite justified: 

 |gravity> irb
 require "vector.rb"
 true
 v = [1, 2, 3].to_v
 [1, 2, 3]
 v * 3.14
 [3.14, 6.28, 9.42]
 (v * 3.14).class
 Vector
 3.14 * v
 [3.14, 6.28, 9.42]
 (3.14 * v).class
 Vector
 quit
And indeed, it all comes out the way it should. So I declare victory: I'm very pleased with our vector class.

15.8. Vector Class Listing

Erica: Let's get the whole listing on the screen. Here it is. In fact, it is shorter than I thought, for all the things it does: 

 class Vector < Array
   def +(a)
     sum = Vector.new
     self.each_index{|k| sum[k] = self[k]+a[k]}
     sum
   end
   def -(a)
     diff = Vector.new
     self.each_index{|k| diff[k] = self[k]-a[k]}
     diff
   end
   def +@
     self
   end
   def -@
     self.map{|x| -x}.to_v
   end
   def *(a)
     if a.class == Vector              # inner product
       product = 0
       self.each_index{|k| product += self[k]*a[k]}
     else
       product = Vector.new            # scalar product
       self.each_index{|k| product[k] = self[k]*a}
     end
     product
   end
   def /(a)
     if a.class == Vector
       raise
     else
       quotient = Vector.new           # scalar quotient
       self.each_index{|k| quotient[k] = self[k]/a}
     end
     quotient
   end
 end
 
 class Array
   def to_v
     Vector[*self]
   end
 end
 
 class Fixnum
   alias :original_mult :*
   def *(a)
     if a.class == Vector
       a*self
     else
       original_mult(a)
     end
   end
 end
 
 class Float
   alias :original_mult :*
   def *(a)
     if a.class == Vector
       a*self
     else
       original_mult(a)
     end
   end
 end

15.9. Forward Euler in Vector Form

Dan: The whole point of this long exercise was to make our codes more readable. Here is one of the array versions we made of our first code, the forward Euler version, in euler_array_each_def.rb: 

 include Math
 
 def print_pos_vel(r,v)
   r.each{|x| print(x, "  ")}
   v.each{|x| print(x, "  ")}
   print "\n"
 end
 
 r = [1, 0, 0]
 v = [0, 0.5, 0]
 dt = 0.01
 
 print_pos_vel(r,v)
 1000.times{
   r2 = 0
   r.each{|x| r2 += x*x}
   r3 = r2 * sqrt(r2)
   a = r.map{|x| -x/r3}
   r.each_index{|k| r[k] += v[k]*dt}
   v.each_index{|k| v[k] += a[k]*dt}
   print_pos_vel(r,v)
 }

What would that look like in vector form?

Carol: To start with, we'll have to require the file vector.rb. We also have to convert the arrays into vectors, using .to_v. Then, in the loop, we can use our vector versions of addition, unary minus, multiplication and division. Let me write it all in file euler_vector.rb: 

 require "vector.rb"
 include Math
 
 def print_pos_vel(r,v)
   r.each{|x| print(x, "  ")}
   v.each{|x| print(x, "  ")}
   print "\n"
 end
 
 r = [1, 0, 0].to_v
 v = [0, 0.5, 0].to_v
 dt = 0.01
 
 print_pos_vel(r,v)
 1000.times{
   r2 = r*r
   r3 = r2 * sqrt(r2)
   a = -r/r3
   r += v*dt
   v += a*dt
   print_pos_vel(r,v)
 }

Dan: Wait a minute, we haven't define += yet!

Carol: Ah, you see, that's a nice feature of operator overloading in Ruby: once you redefine + the same redefinition applies to derivative expressions such as +=.

Dan: That's nice, if it really works. Let's compare the run with the old results: 

 |gravity> ruby euler_array_each_def.rb | tail -1
 7.6937453936572  -6.27772005661599  0.0  0.812206830641815  -0.574200201239989  0.0  

 |gravity> ruby euler_vector.rb | tail -1
 7.6937453936572  -6.27772005661599  0.0  0.812206830641815  -0.574200201239989  0.0
Wonderful. And I must say, the code does look a lot cleaner now.

Erica: Definitely. This is a lot more pretty. What a difference! None of all those ugly [k] occurrences anymore. I think it was worth all the work we put in, defining a Vector class.
Previous ToC Up Next