Feeling Guilty

Today on #perl6, masak commented on the many failing tests in the various projects in proto. Unfortunately, the Vector project is one of those that fails a test. The core Vector stuff passes with flying colors (at least on my Rakudo build here) but my most recent test file on this side project fails. It's not even in failing in my code -- it's just a test I wrote for something very basic in Perl 6 which doesn't seem to work in this test. I've no idea why, it seems to work fine stand alone.


use v6;
use Test;

plan *;

my @array = (1, 2, 2, 3, 4, 5, 5, 5, 5, 6, 7, 8);

ok([<=] @array, "array is sorted properly");

done_testing;


If anyone has a notion why this fails, I'm all ears. I'd love to clear this up.

Polynomial and <<+>>

Defining a one-variable polynomial class is actually a pretty straightforward thing. We use a standard Perl 6 array to hold the polynomial's coefficients, with the Nth element of the array representing the coefficient of x^N. You can check out the complete source if you'd like; for the most part, Perl 6's nifty hyperoperators and list operations make the implementation simple.

One example where that isn't true is simple polynomial addition. In principle, this is the simplest thing possible; you just do a pairwise add on all the coefficients of the two polynomials. The only catch is when the polynomials have different degrees, so the arrays have different sizes. In that case we'd like to pad the array with zeros. Unfortunately, <<+>> pads with the last element of the shorter array, with is completely mathematically unsuitable in this case. So here was my first stab at implementing Polynomial addition:


That sort of thing works, but seems more like it belongs in the ugly world of C++ than the elegant world of Perl 6. The subject came up on the #perl6 channel yesterday, and pmichaud++ came up with a nice workaround: $a.coefficients, 0 <<+>> $b.coefficients, 0. As it is, this adds an additional zero coefficient each time it is called, which in the long run will be trouble. But it can be counteracted by adjusting the Polynomial constructor to remove "trailing" zero coefficients. (This is how Polynomial.pm is currently implemented.) (Errr... perhaps I spoke too soon; I'm having difficulty successfully deleting the trailing zeros. Sigh.)

On the IRC chat, colomon suggested using an adverb to tell <<+>> how to extend. His suggestion of :extend-with seems like a weak name to me, but I do think it seems like the proper Perl 6ish way to handle this issue.

The Best Laid Schemes...

"The best laid schemes o' Mice an' Men, Gang aft agley" -- Robert Burns

I've hit on one of those programming conundrums where a quick diversion becomes a huge project. Let me explain what I'm going for; it will probably take three or four blog posts to cover all the details.

My basic idea for a cool Perl 6 twist on knot vectors was evaluating the basis function to create a polynomial rather than a set value. The basis function defines a piecewise polynomial, with a continuous polynomial between each pair of distinct neighboring knots. So if the knot vector is (0, 0, 0, 1, 2, 3, 3, 3), that defines three polynomials: from 0 to 1, from 1 to 2, and from 2 to 3.

For any reasonably complicated knot vector, those polynomials are pretty tricky to calculate by hand. The NURBS Book approach is a change of basis, which works well but is not particularly mathematically enlightening. But the knot vector basis calculation function N we've defined can already do this: if it is passed a reasonably capable polynomial object for the $u parameter. This is the payoff for not forcing $u to be a Num; we can pass a polynomial instead and get a polynomial out.

Well, almost. We have two KnotVector.N functions, one for the $p == 0 case and one for $p > 0. The latter will work perfectly when passed a polynomial for $u. (At least in theory; there may be issues with Rakudo bugs, of course.) The former, however, actually needs to know what the numeric value of $u is. Originally I planned to somehow overload the polynomial class to have a value for this purpose. But as far as I know, to make that work, I'd have to overload prefix:<+> -- and last time I checked, that didn't work in Rakudo.

You could fix that by changing the second N to take an impulse array, instead of calculating it using the value of $u. But once you're thinking about doing that, it's time to face the facts: even the hypothetical "cool" implementation of N we have is grotesquely inefficient for reasonably complicated knot vectors.



To calculate the basis efficiently, you should only calculate the values that have a chance of being non-zero. And you should get the impulse value by doing something like a binary search.

So before I can get where I want to go with knot vectors, I need to have a working polynomial class and a working binary search. Plus I think it might be worthwhile to have a short post discussing the .Str and .perl functions for KnotVector and Nubs...

NURBS Knot Vectors In Perl 6 (part 3)

Well, instead of being a slick demonstration of how powerful Perl 6 is, it has turned into a slog to find something the current Rakudo is happy with. The beautiful hyper-equation of my last post has turned into a dull and confusing series of little hyperoperations (and one which got expanded because the hyperoperator just wouldn't work).


I think all the bugs I've encountered here have already been reported, but I'm not sure. If anyone has suggestions on how to ably search RT, I'm all ears. My attempts so far have been painfully clumsy, to the point where my inclination is either to not report if I think it's likely to have been reported so far, or just report it without searching.

At any rate, let me finish by quickly sketching out a NUBS (Non-Uniform BSpline) curve class. This isn't quite as powerful a tool as a NURBS (Non-Uniform Rational BSpline), but it's still pretty powerful, and it's dead easy to implement the basic evaluate function with the tools we have so far.


Combined with the Vector class, this allows you to easily define N-dimensional NUBS curves. But we don't force the control points to be Vectors -- you can get a NUBS representation of polynomials over reals or complexes easily using this too, and if you have your own Vector class you'd prefer, that should work too, as long as you have scalar multiplication and vector addition implemented using * and +.

I do have one interesting twist left for another post on the KnotVector class, something I've always wanted to do that wasn't very practical in C++. And then I will work on expanding Nubs, and implmenting a proper Nurbs class, and looking at more complicated structures for each (surfaces, etc).

NURBS Knot Vectors In Perl 6 (part 2)

This is the sort of frustrating post you occasionally get when working with a compiler in the process of development like Rakudo. Because the cool thing I want to do is legal Perl 6 (I'm pretty sure), but it doesn't come close to working in Rakudo.

When we last left knot vectors, we had defined a recursive function for calculating the knot vector basis. The problem with that implementation is each level calls the prior level twice. In practice, this means most of values are calculated repeatedly. For higher order NURBS objects, this is massively inefficient. So ideally we'd like to calculate each value of N exactly once. This can be done by calculating all the degree 0 values of N, then use those results to calculate all the degree 1 values, then use those to calculate the degree 2 values, and so on. My big idea here was that with Perl 6's hyper operators, each level can be calculated without even using a loop! So given an array of knots named @.knots, the prior degree's array of N values called @N_prev, and $end = @.knots.end, the calculation looks like this:


I think that is correct. But it makes Rakudo very unhappy: "Non-dwimmy hyperoperator cannot be used on arrays of different sizes or dimensions." Mind you, that's happening in $u «-« @.knots[0..($end - $p - 1)], which is very definitely a dwimmy hyperoperator. And that's just the tip of the iceberg -- in fact, I have yet to simplify this calculation enough to make Rakudo happy with it. It seems Rakudo is really not ready to handle deeply nested hyperoperators yet. You can work around it in this case by using the Texas version of the operator. It also seems there is some odd glitch the second time you call »O/«: you get the error message "ResizablePMCArray: Can't pop from an empty array!" even simplifying the use until it is very obviously right.

I believe both bugs have been reported already, so there's nothing to do but dumb down the code so that it will work with the current Rakudo. Sigh.

NURBS Knot Vectors In Perl 6 (part 1)

I'm not going to provide a lot of background for the wheres and whys of Non-Uniform Rational B-Spline (NURBS) knot vectors yet. Basically, a NURBS object is one or more knot vectors which define the parametrization and polynomials of the object, along with a bunch of Vectors to define the frame of the shape. I'm planning on spending a couple of posts working on getting a good knot vector class built up.

Confusingly enough, knot vectors have nothing directly to do with the sort of Vector we have been working with so far. From, say, a Perl 6 perspective, it would probably be more correct to call them knot arrays or knot lists. But knot vector is the traditional term, so that is what we will use.

The knot vector is a series of numbers in non-descending order, some of them repeated, that define the parameter range of the NURBS object we are looking at, and where it does interesting things. That is to say, adding a knot to the middle of the knot vector essentially splits one polynomial span the NURBS object is defined over into two joined polynomials. It's a fairly simple basis for a lot of powerful geometry.

Instead of starting by defining a knot vector class which can be the basis for NURBS objects, we're going to start by implementing the definition of the B-spline basis functions directly. (Note that the Wikipedia version there specifies it using different variable names -- our notation comes from The NURBS Book.) This would be horribly inefficient for real-world calculations, but we're just going to use it to establish good tests for our real class when that is ready.

Here's a straight forward implementation of the basis function formulas.



There are two tricky bits here. The first is that the standard definition can result in dividing zero by zero, which is defined to result in zero in this case. To keep the function implementation clean and simple, we define a "O/" operator which implements this exception to the standard rules of division.

The second wrinkle is that the standard definition has one curious blind spot. The equation does not work out when you reach the very last point! That is to say, the B-spline basis defined the standard way on a knot vector whose first value is A and last value is B covers the range [A, B). (Actually, I'm using first and last rather loosely there -- that's only true with the most common kind of knot vectors.) We define a KnotBasisDirection argument to N allows you to switch to the symmetrical definition, which is defined over (A, B]. The two definitions are identical over the interior range (A, B). The KnotBasisDirection enum has two values, Left and Right, corresponding to these choices. That's the first time I've defined an enum in Perl 6, but it seems straightforward enough.

Thoughts about numeric roles in Perl 6

There was a big discussion about potential numeric classes/roles in Perl 6 on #perl6 a few days ago. Most of it went right over my head, I fear. But I've kept on thinking about it in the days since, and it seems to me there is a pretty simple and mathematically sound layout for the roles. (Instead of worrying about coming up with beautiful names, I'm just going to stick "Role" onto the name of each mathematical concept.) This is a quick high level sketch which I have made no attempt to implement.

The biggest role would be ComplexRole, capturing the idea of numbers with real and imaginary parts (either or both of which might be zero), and implementing a slew of standard mathematical functions that work on complex numbers.

Then there would be a subset RealRole of ComplexRole where { $^z.im == 0.0 }. The most notable added method would be Num, but there would probably also be a lot of RealRole-specific mathematical functions, as they can be a lot simpler than the complex equivalents.

I don't know how to express is as a where, but the RationalRole subset of RealRole would mostly just add numerator and denominator methods. And finally we would have subset IntegerRole of RationalRole where { $^q.denominator == 1 }.

The idea would be that the classes would be built up in the opposite order than this (or not related at all), using the logical computer programming approach, with these roles applied on top to give a mathematical structure. So you'd define a rational class as two integers, then define numerator and denominator methods for the RationalRole, a Num method for the RealRole, and an imaginary method for the ComplexRole.

I'm not sure how practical this is to actually implement in Perl 6, but it seems like it would be a nice way to capture these notions. I'm presuming most of the mathematical functions would be implemented on ComplexRole and RealRole, so this would provide a simple way of plugging any new math types you defined into a full library of math functions.

Update: I've started work on implementing B-spline knot vectors and basis functions in Perl 6, but it will be a few more days before I have something worth posting.

Potential Cool Use (Extension?) For Range

While working on the tests for my next piece of code to talk about on this blog, I ran into one of those perennial nuisances of numeric programming. It is very common to need a list of numbers from A to B (inclusive) with either N numbers or split into S segments. So perhaps I want ten numbers from 1.0 to 2.0; then the numbers are 1.0, 1.11111, 1.22222, 1.33333, 1.44444, 1.55555, 1.66666, 1.77777, 1.88888, 2.0. There are two issues here. The first is the potential for fencepost error: if I think ten numbers from 1.0 to 2.0, my first thought is always that they go up by 0.1, when they actually go up by 1/9th. The second issue is that if you do the obvious loop ($x = 1.0; $x <= 2.0; $x += (1/9)), it's a crapshoot whether you get the "2.0" or not, and even if you do, it might actually be 1.99998 or something like that. Usually in these cases getting the exact ending value is highly desirable (as it is an edge case).

(Before getting into potential solutions for this, let me note that this particular example is already solvable in the latest Rakudo builds using Rats. loop ($x = 1; $x <= 2; $x += (1/9)) should give the exact list 1, 10/9, 11/9, 12/9, ... 17/9, 2. But that doesn't help with the fencepost, and the general problem is very ugly to solve this way!)

My first thought was that Range with :by might do the trick. Unfortunately, I don't see Perl 6 specs for how :by actually works in this sort of case, and as far as I know, it's not actually implemented yet to test. At any rate, I'd like to suggest that :by does actually handle this case, if possible (ie if an element of the range is within epsilon of the to value, then the to value is returned instead). And I'd like to suggest two more modifiers, something like :size(N) to specify that the Range should return N values, with the :by value calculated as appropriate, and :sections(N) to specify that the Range should return N+1 values.

In the meantime, here's a quick implementation of what I think the Range :size function should do. I believe that it is properly set up to be lazy, while laziness is available. Note that if you use integer parameters, it will return Rats!