During the last few weeks I've been working on a library that would let the user do some differential-geometric calculations in Rust. By differential geometry I mean mostly the tensor calculus in curved spaces or space-times. I've already created something like that in C++, but I wanted to try and use some of the Rust features to create an improved version.
What could Rust do better?
The most convenient representation of tensors for doing calculations is in the form of arrays of numbers. The problem is that representing a tensor numerically requires choosing a coordinate system. Various operations, like for example addition of two tensors, only make sens when the tensors involved are expressed in the same coordinate system. The only possibility of enforcing this rule in C++ was to encode the coordinate system as a property of the tensor object and checking for compatibility in the operator code. This way any errors will be detected at runtime.
Ok, so the errors were detectable, so what could be done better? Well, for examples the tensors expressed in different coordinate systems could not only have a different value of some property, but be objects of different types. This way the error can be detected at compile time, before the program is even translated into an executable form. It wouldn't be very practical in C++, but the Rust type system allows to do it quite interestingly.
EDIT: It has been brought to my attention that C++'s templates also allow for this kind of thing. Nevertheless, doing it in Rust was a fun experiment :)
The array problem once again
So we want to check the compatibility of tensors with respect to the coordinate systems statically (that is, at compile time), but we could expand it a bit. For example, addition only makes sense for tensors of the same rank and variance (that is, the same "composition" of indices, that can be either covariant or contravariant). Adding, let's say, a vector and a covector (which are both of rank 1, but one is contravariant, and the other covariant) makes no sense. If the type system could also be used to detect this kind of errors, it would be awesome, but is it possible?
As it turns out, yes it is, but it's not that easy.
Dealing with tensors
Tensors are characterized basically by two properties: the dimension of the underlying space and their rank / variance. For example, as mentioned above, vectors and covectors have rank one. Matrix tensors (like linear transformations or bilinear forms), on the other hand, have rank 2. The difference lies in the variance - linear transformations have one contravariant and one covariant index, and bilinear forms - two covariant ones. The dimension along with the rank define the number of coordinates of the tensor - it is 
, where 
is the space dimension, and 
- the rank of the tensor. Since both these values will be known statically, it would also be nice if we could statically determine the size of the coordinates array. And this is where problems start.
As I already mentioned in the previous post, the Rust arrays can't be parametrized by their length. At the first glance it would be thus necessary to program each type of tensors separately with an array of an appropriate size. Fortunately, there is a solution which I also described in that post - it is the GenericArray struct. The types used by this struct to encode the array length have all of arithmetics defined in the typenum crate, so calculating values like shouldn't be a problem. Indeed, it's possible, but not really smooth.
It turns out that each operation on the number types from typenum causes all the guarantees assumed for those types to be lost. For example, if I defined the type representing the dimension and the type representing the rank to be applicable as the length of an array, this guarantee disappears upon creation of the type being the result of . This means that I need to define it separately, which bloats the trait requirements and leads to monster snippets like this one. The author of typenum, paholg has been of invaluable help in this matter. Thanks to him, I managed to write requirements that allowed me to compile a struct representing a tensor.
After these initial struggles everything went smoothly. Further operations on tensor types were easy to code, as I already knew what the compiler expects. This way I quickly implemented tensor addition and subtraction, multiplication by scalars and other tensors, and contraction (which means summing the tensor elements on a "diagonal" - in the case of a matrix it's just the trace).
The compiler gives up
Unfortunately, after I had everything coded and a few tests written, it turned out that the compiler couldn't take it anymore. The tensor multiplication was too much to bear and it got in an infinite loop.
My efforts aiming to identify and repair the problem were in vain. It is pretty obvious that the culprit is the tensor multiplication, as commenting out its code makes the error disappear, but it's impossible to determine how the error is being caused. In an act of despair I compiled the debug version of the compiler and generated some logs, but I only got about 500 MB of text which didn't help me much. The reddit community advised me to start an issue on GitHub, which I did. One of the developers manage to extract the essence of the offending code, but the cause of the problem is still unknown.
Nevertheless, using an alternate syntax (<Tensor<T,U> as Mul<Tensor<T, V>>::mul(tensor1, tensor2) instead of tensor1 * tensor2) allows to get the code to work and the library is functional. I have released version 0.1 yesterday.
Finally: an additional macro for GenericArray
I've also managed to improve one particular aspect of GenericArray. Until now, creation of such an array required the usage of the from_slice method, which a) bloats the code, b) only validates the slice length dynamically. Today I created an arr! macro, which not only makes the code shorter, but also checks everything statically:
|
1 2 3 4 5 6 7 8 9 |
// the old version let array1 = GenericArray::<u32, U3>::from_slice(&[1, 2, 3]); // the new version let array2 = arr![u32; 1, 2, 3]; // compiles, but panics at runtime let array3 = GenericArray<u32, U2>::from_slice(&[1, 2, 3]); // doesn't compile let array4: GenericArray<u32, U2> = arr![u32; 1, 2, 3]; |
The macro was uploaded to crates.io in version 0.1.1 of generic-array.
Summary and future plans
So, this was my recent coding activity. In the near future, I'm planning to add some support for metric tensors and Christoffel symbols to differential-geometry, which should be enough to create a clone of the aforementioned gr-engine. When I'm done, I'll write what became of it :)