• The Expression Problem in Rust

    Eli Bendersky’s website is fantastic and you should read it. In one series of posts he discusses the expression problem: a software design problem where in some languages it is easy to add new data types and in others it is easy to add new operations on data types. I like thinking about the expression problem because similar concepts appear in the design of so-called “deep learning” frameworks like MXNet and Tensorflow where operations and data types are abstracted out and treated symbolically.

    To summarize his excellent post, in object-oriented languages it is easy to add new data types. However, to add a new operation you need to modify the object interface and revisit each data type definition to implement the new operation. (The key word here, I believe, is revisit in the sense that the modification violates the open-closed principle.) Alternately, in functional programming languages the situation is reversed: adding a new data type means having to modify the definition of each operation, also violating the open-closed principle.

    I believe that Rust solves this problem using traits. In this post I will demonstrate how this is so. My knowledge of Rust is really elementary at the time of writing this so hopefully I haven’t violated any major Rust principles. My code may not be of top-level Rust design but it works.

    A Motivating Example

    Let’s replicate Eli’s example in Rust: a simple expression evaluator. Although he brings up this example in the context of compilers and interpreters expression construction and evaluation is also the bread and butter of deep learning frameworks.

    Eli’s examples use a base class (in OOP) or enum (in FP) called Expr from which all expression types are derived. From my reading of Rust such a notion may not be necessary, or perhaps even appropriate, if the goal is to perform operations on data types. Therefore, in my definition of Constant I will omit any “inheritance” from a parent type. See my final thoughts at the end of this post for more about design decision.

    We begin with the definition of a Constant data type,

    #[derive(Debug)]
    struct Constant<T> {
        value: T
    }
    

    We include the derive(Debug) macro to auto-magically make the data type printable for debugging purposes. Below is some code using this struct,

    fn main() {
        let c1 = Constant {value: 1.2};
        let c2 = Constant {value: 3.4};
    
        println!("c1 = {:?}", c1);  // output: c1 = Constant { value: 1.2 }
        println!("c2 = {:?}", c2);  // output: c2 = Constant { value: 3.4 }
    }
    

    Let’s now make Constant expressions “evaluatable”. To do so, we’ll create a new trait Eval which implements an eval() method. For Constant expressions eval() should simply return a copy of stored value. Note that we could implement eval() as a struct method but the end goal, here, is to have shared operation functionality across multiple data types.

    trait Eval<T> {
        fn eval(&self) -> T;
    }
    
    impl<T> Eval<T> for Constant<T>
        where T: Clone
    {
        fn eval(&self) -> T {
            self.value.clone()
        }
    }
    

    Because of Rust’s ownership model we return a copy/clone of Constant.value so that the struct can retain ownership of its value. Otherwise, the contents of value are moved out of the Constant struct to the caller. (Aside: is there a better way to do this?) Let’s test that (a) the eval() implementation works and (b) that the ownership properties are indeed satisfied.

    Continuing the main() function from above,

    fn main() {
        // --snip--
    
        println!("c1.eval() = {}", c1.eval());  // output: c1.eval() = 1.2
        println!("c1.eval() = {}", c1.eval());  // output: c1.eval() = 1.2
        println!("c2.eval() = {}", c2.eval());  // output: c2.eval() = 3.4
        println!("c2.eval() = {}", c2.eval());  // output: c2.eval() = 3.4
    }
    

    Excellent. If the value stored in c1 were moved then the second invocation of println!() would cause a compilation error. It’s incredible that Rust won’t let you compile if you have memory leaks!

    Adding a New Type

    Let’s add the second data type in Eli’s example.

    #[derive(Debug)]
    struct BinaryPlus<'a, 'b, U, V>
    where U: 'a,
          V: 'b,
    {
        left: &'a U,
        right: &'b V,
    }
    

    Again, I am not yet confident in my Rust skills but my understanding is that the BinaryPlus type should not actually own the data stored in its left and right attributes. In particular, we want the lifetime of BinaryPlus’s left attribute to last at least as long if not longer than the BinaryPlus instance, itself. Hence, the lifetime specifications. Regardless, trying to compile this data type without lifetime parameters results in a compile error.

    We now implement the Eval trait for BinaryPlus,

    impl<'a, 'b, T, U, V> Eval<T> for BinaryPlus<'a, 'b, U, V>
    where T: Add<Output=T>,
          U: Eval<T> + 'a,
          V: Eval<T> + 'b,
    {
        fn eval(&self) -> T {
            self.left.eval() + self.right.eval()
        }
    }
    

    Note the requirements on the generic data types. Whatever type T is must be add-able with output type also T. The left and right arguments must be Eval-able; after all, we call eval() on them. What’s interesting is that these restrictions are not present in the actual definition of the BinaryPlus type, above. Only if you want to use the Eval functionality do you have to narrow down the traits of the arguments.

    Here is an example usage of this new type and operation.

    fn main() {
        // --snip--
    
        let bp1 = BinaryPlus { left: &c1, right: &c2 };
        println!("bp1 = {:?}", bp1);              // output: bp1 = BinaryPlus { left: Constant { value: 1.2 }, right: Constant { value: 3.4 } }
        println!("bp1.eval() = {}", bp1.eval());  // output: bp1.eval() = 4.6
    
        // ownership check: make sure bp1 doesn't now own c1 and c2
        println!("c1.eval() = {}", c1.eval());    // output: c1.eval() = 1.2
        println!("c2.eval() = {}", c2.eval());    // output: c2.eval() = 3.4
    
        let bp2 = BinaryPlus { left: &c1, right: &bp1 };
        println!("bp2 = {:?}", bp2);              // output: bp2 = BinaryPlus { left: Constant { value: 1.2 }, right: BinaryPlus { left: Constant { value: 1.2 }, right: Constant { value: 3.4 } } }
        println!("bp2.eval() = {}", bp2.eval());  // bp2.eval() = 5.8
    

    Note that in adding a new type we did not have to modify and previous definitions of the Eval trait. Instead, we just specified how the Eval trait behaves with this new type. Therefore the FP version of the expression problem is resolved; we haven’t violated the open-closed principle.

    Adding a New Operation

    Time to add the second example operation to_string() as a method of the trait Stringify and implement it for both the Constant and BinaryPlus data types.

    trait Stringify {
        fn to_string(&self) -> String;
    }
    
    use std::fmt::Display;
    
    impl<T> Stringify for Constant<T>
    where T: Display
    {
        fn to_string(&self) -> String {
            format!("Constant({})", self.value)
        }
    }
    
    impl<'a, 'b, U, V> Stringify for BinaryPlus<'a, 'b, U, V>
    where U: Stringify,
          V: Stringify,
    {
        fn to_string(&self) -> String {
            format!("BinaryPlus({} + {})", self.left.to_string(), self.right.to_string())
        }
    }
    

    Note the type restrictions for Constant’s implementation is just that T implements the Display trait. As for BinaryPlus, we require that each constituent attribute is Stringify-able.

    Let’s see this new functionality in action:

    fn main() {
        println!("c1.to_string() = {}", c1.to_string());    // output: c1.to_string() = Constant(1.2)
        println!("c2.to_string() = {}", c2.to_string());    // output: c2.to_string() = Constant(3.4)
        println!("bp1.to_string() = {}", bp1.to_string());  // output: bp1.to_string() = BinaryPlus(Constant(1.2) + Constant(3.4))
        println!("bp2.to_string() = {}", bp2.to_string());  // output: bp2.to_string() = BinaryPlus(Constant(1.2) + BinaryPlus(Constant(1.2) + Constant(3.4)))
    }
    

    Nice. Again we see that adding the Stringify trait does not require modifications to the previously-written code for the Constant and BinaryPlus data types. Instead, we modify the capabilities of these types with new trait implementation code. Therefore, Rust also resolves the OOP version of the expression problem.

    I suppose the caveat, here, is if you wanted to modify the definition of trait Stringify to include a new method. Here is one way to do it without having to return to the definition of Stringify,

    trait StringifyExt {
        fn to_alternate_string(&self) -> String;
    }
    
    impl<T> StringifyExt for T
    where T: Stringify
    {
        fn to_alternate_string(&self) -> String {
            String::from("Hey!")
        }
    }
    
    
    fn main() {
        // --snip--
    
        println!("c1.to_alternate_string() = {}", c1.to_alternate_string());  // output: Hey!
    }
    

    Once again, no existing code was modified. I added the “extension trait” and implemented it for any type which implements Stringify. This is not quite the same pattern as I read in this article but, then again, Option<T> is a type (enum), not a trait. Regardless, this technique is useful if you want to extend traits that don’t belong to you.

    Final Thoughts

    Learning Rust has been educational. Just like when I played around with Haskell I came away looking at code design a little bit differently in other languages. (My post-Haskell Python code used more functional than iterative styles.)

    Earlier I mentioned that using an Expr enum seemed unnecessary for modeling the expression problem in Rust. Strictly speaking, to accurately model the problem I should have created an enum like so,

    enum Expr {
        ConstantExpr(Constant),
        BinaryPlusExpr(BinaryPlus),
    }
    

    The user could then create ConstantExpr and BinaryPlusExpr instances, both of which are of type Expr. The implementation details would be “hidden” underneath in the actual underlying structs, Constant and BinaryPlus, respectively. My thought is the use of an Expr enum does not add anything to resolving the expression problem: the only goal is to have a model where we can add data types and operations on those data types independently and easily such that the open-closed principle isn’t violated.

    That being said, Rust enums cannot be extended and therefore the open-closed principle is violated. However, there is a clunky workaround.

    enum MoreExpr {
        Expr(Expr),
        SomeNewExpr,
    }
    
    match my_expr_variable {
        Expr(ConstantExpr) => ...,
        Expr(BinaryPlusExpr) => ...,
        SomeNewExpr => ...
    }
    

    The role of the Expr interface in object-oriented programming languages is to establish a shared interface across all sub-types. The role of the Expr enum in functional programming languages is so that you can perform matching on various sub-types. With Rust traits the data types and operations on those data types are much more separated since you can just check if an input type implements a trait. Therefore, for the purposes of this problem I’m comfortable in omitting the Expr enum and calling the problem solved. Please prove me wrong because it seems too good to be true.


  • Announcement - Looking for New Abelfunctions Owners

    As much as I would like to keep hacking on Abelfunctions I doubt that I will have time to do so once I have a “real job”. Too many student projects evaporate into the ether. Therefore, I am hoping someone will take the mantle of this project and guide it to a better state.

    Please contact me if you are interested.

    Though not required, since I assume you’re a bright mathematician, but experience in the following would help in developing for Abelfunctions: (in order of importance)

    • Python / Sage
    • Object-Oriented programming
    • Cython
    • git / GitHub

    Of course, an appropriate mathematical background is needed. This project performs calculations that use the following:

    • algebraic geometry / planar curves
    • numerical analysis
    • computational geometry
    • special functions (in particular, the Riemann theta function)

    Help

    If you are interested in taking on this project I would be happy to meet with you via Facetime, Google Hangouts, Skype, or whatever your preferred video chat + screensharing method. I tried to document my code well but as some Ph.D. research projects the rush of arriving at certain checkpoints causes documentation and clarity take a hit while searching for correctness.

    Additionally, it’s not like I will cease to exist. I will be around to offer the occasional guidance and opinion.


  • Transitioning to Native Sage Datatypes

    Nils Bruin at Simon Fraser University brought up the idea / request to compute period matrices of plane curves to high precision. Sage has built-in capabilities to work with arbitrary-precision complex numbers via the CC ring. It also implements complex doubles in the the CDF ring.

    In the original version of Abelfunctions I decided to use Numpy + Cython for numerical performance. Now I’m wondering if it’s worth investigating using native Sage datatypes instead. Below, I perform some timings between different fast_callable declarations, where the domain is varied, and different input datatypes: numpy.complex, CDF, and CC.

    Setup

    from sage.all import fast_callable, QQ, CDF, CC, timeit
    import numpy
    import time
    
    R = QQ['x,y']
    x,y = R.gens()
    
    # functions:
    f = y**3 + 2*x**3*y - x**7
    f_complex = fast_callable(f, vars=[x,y], domain=complex)
    f_CDF = fast_callable(f, vars=[x,y], domain=CDF)
    f_CC = fast_callable(f, vars=[x,y], domain=CC)
    
    # inputs
    a = numpy.complex(1.000000000007 + 1.0000000000009j)
    b = numpy.complex(1.00000000005 - 2.0000000000008j)
    a_CDF = CDF(a)
    b_CDF = CDF(b)
    a_CC = CC(a)
    b_CC = CC(b)
    
    N = 100000

    Timings

    timeit('for _ in xrange(N): f(a,b)')
    timeit('for _ in xrange(N): f(a_CDF,b_CDF)')
    timeit('for _ in xrange(N): f(a_CC,b_CC)')
    5 loops, best of 3: 2.66 s per loop
    5 loops, best of 3: 2.49 s per loop
    5 loops, best of 3: 3.47 s per loop
    timeit('for _ in xrange(N): f_complex(a,b)')
    timeit('for _ in xrange(N): f_complex(a_CDF,b_CDF)')
    timeit('for _ in xrange(N): f_complex(a_CC,b_CC)')
    5 loops, best of 3: 234 ms per loop
    5 loops, best of 3: 271 ms per loop
    5 loops, best of 3: 275 ms per loop
    timeit('for _ in xrange(N): f_CDF(a,b)')
    timeit('for _ in xrange(N): f_CDF(a_CDF,b_CDF)')
    timeit('for _ in xrange(N): f_CDF(a_CC,b_CC)')
    5 loops, best of 3: 1.18 s per loop
    25 loops, best of 3: 20 ms per loop
    5 loops, best of 3: 88.3 ms per loop
    timeit('for _ in xrange(N): f_CC(a,b)')
    timeit('for _ in xrange(N): f_CC(a_CDF,b_CDF)')
    timeit('for _ in xrange(N): f_CC(a_CC,b_CC)')
    5 loops, best of 3: 2.4 s per loop
    5 loops, best of 3: 1.73 s per loop
    5 loops, best of 3: 1.18 s per loop

    The fastest result using only native Sage datatypes is f_CDF(a_CDF, b_CDF). This is only about x2 slower than directly during f_CDF.call_f(...) with native C complex datatypes. (Here, f_CDF is of type Wrapper_cdf which provides low-level access to the underlying JIT-compiled C function.) That is, by using native C get a x2 performance increase from the fastest native Sage version. Given that the pure-CDF version is already x100 faster than a non-CDF implementation I don’t think it’s worth the last bit of performance if it means easier to read code. (Not to mention, the ease of transitioning to using CC.)

    It will be challenging refactoring the numerics code for Sage. I will most likely not do anything about it until I’m done with my Ph.D. but it’s worth investigating in the future.