Introduction to metaprogramming: "Code that creates code"

This is pretty much borrowed from David Sanders Julia 2016 talk.

Julia has strong metaprogramming capabilities. What does this mean?

meta: something on a higher level

metaprogramming = "higher-level programming"

i.e. writing code (a program) to manipulate not data, but code (that itself manipulates data)

Motivating example

Metaprogramming has many different uses, several of which we will explore. It is often a way to add, or change, the syntax of Julia, to write a "domain-specific language" (DSL), that converts a simple syntax (but that is not standard Julia syntax) into true Julia code.

As a motivating example, we might like to be able to write what is called Wilkinson-type polynomials. The Wilkinson polynomial is

p20(x):=(x−1)⋅(x−2)⋅⋯⋅(x−20)=∏i=120(x−i).p_{20}(x) := (x-1) \cdot (x-2) \cdot \cdots \cdot (x-20) = \prod_{i=1}^{20} (x-i).

[Polynomials like this are interesting, since the eigenvalues of the associated "companion matrix", which are used to find the roots of the polynomial, are very sensitive to perturbations of the coefficients of the polynomial.]

Suppose we wish to define this polynomial in Julia. The simple way would be to write it out explicitly:

p_5(x) = (x-1) * (x-2) * (x-3) * (x-4) * (x-5)

p10p_{10} is already a pain to type by hand, p20p_{20} more so, and p100p_{100} is basically impossible. But this is just a case of repetition, and computers are designed for that. One possible definition uses a for loop and with the help of anonymous function to define the function pnp_n:

function wilkinson(n, x)
    result = x - 1
    
    for i in 2:n
        result *= (x - i)
    end
    
    result
end

p(n) = x -> wilkinson(n, x)

p(4)(3.1)

It seems, though, that it should be possible to use the original definition of the function to write the equivalent of the following for n=100n=100.

p_5(x) = (x-1)*(x-2)*(x-3)*(x-4)*(x-5)

In other languages, this would often be accomplished by manipulating strings that represent the "surface syntax", i.e. the string of characters that you would actually type, and then evaluate this string. However, in Julia

we never use strings for this, since Julia has a way to refer to Julia code objects within Julia. So Julia codes are objects and you can also manipulate them like manipulate any other things in Julia.

Expressions

Let's take the simplest polynomial, that we write as (x-1) * (x-2). We can view this as a piece of Julia code, called an expression, and can tell Julia this as follows:

ex = quote   
        (x - 1) * (x - 2)
    end

# or direct way

ex = :( (x-1) * (x-2) ) 

typeof(ex) #What does Julia think this is?

ex

We see that Julia has a type Expr that represents an expression object, that we can think of as a "Julia code snippet". Now if we would like to be able to execute this code; we can do so with eval ("evaluate"):

eval(ex) # We see that `x` is not defined,

If we give x a value first, then it works:

x = 3.5
(x-1) * (x-2) # works
eval(ex)

For example, we can try to write a simple formula evaluator:

formula = "3x^2 + 4x - 3sin(x)"
eval(formula) #This does not do what we expect, since we have a string, not a Julia expression:
typeof(formula)
formula2 = Meta.parse(formula)  # To convert this string into an expression, we use `parse`:
eval(formula2)

Note that the really hard work, that of parsing the expression, i.e. converting it into an internal representation in terms of a tree, has already been done for us by Julia. It is thus, happily, usually not necessary for us to write our own parser; we just leverage Julia's own! This is one of many ways in which Julia minimizes the work that we need to do, compared to other languages.

Structure of Expressions

An expression object of type Expr has an internal structure that corresponds to the abstract syntax tree (AST) of the expression, i.e. a tree, whose nodes are operations, and whose children are subtrees. We can see this in two ways, using dump:

ex = :((x-1)*(x-2))
dump(ex)


## print the tree with the nice package
using AbstractTrees
print_tree(ex)



julia> dump(ex)
Expr
head: Symbol call
args: Array{Any}((3,))
    1: Symbol *
    2: Expr
    head: Symbol call
    args: Array{Any}((3,))
        1: Symbol -
        2: Symbol x
        3: Int64 1
    3: Expr
    head: Symbol call
    args: Array{Any}((3,))
        1: Symbol -
        2: Symbol x
        3: Int64 2

julia> using AbstractTrees

julia> print_tree(ex)
Expr(:call)
├─ :*
├─ Expr(:call)
│  ├─ :-
│  ├─ :x
│  └─ 1
└─ Expr(:call)
├─ :-
├─ :x
└─ 2

which shows everything [up to a pre-determined depth] in detail, or

Meta.show_sexpr(ex)  ## the Lisp way, but let's ignore this

(:call, :-, (:call, :+, (:call, :*, 3, (:call, :^, :x, 2)), (:call, :*, 4, :x)), (:call, :*, 3, (:call, :sin, :x)))

Let's stick with dump

ex = :((x-1) * (x-2))
dump(ex)

Expr 
  head: Symbol call
  args: Array(Any,(3,))
    1: Symbol *
    2: Expr 
      head: Symbol call
      args: Array(Any,(3,))
        1: Symbol -
        2: Symbol x
        3: Int64 1
      typ: Any
    3: Expr 
      head: Symbol call
      args: Array(Any,(3,))
        1: Symbol -
        2: Symbol x
        3: Int64 2
      typ: Any
  typ: Any

Now we can also access what is inside the type, using our .fieldnames syntax

dump(ex)

ex.head
ex.args

julia> ex.head
:call  ## a function call

julia> ex.args
3-element Array{Any,1}:
 :* 
 :(x - 1)
 :(x - 2)

We can also go one level further. The first element of the array ex.args is the function to be called:

ex.args[1]
ex.args[2]
typeof(ans) # that is, it is itself an `Expr`ession:

Thus the structure of Expressions is fundamentally recursive, and so lends itself naturally to recursive algorithms. We can dive, in turn, into the structure of each element:

ex2 = :(a = b + c)
ex2.head

ex.args[2]


ex.args[2].args


ex.args[2].args[1]

Also with the example we started

ex.head
ex.args
ex.args[2]
ex.args[2].args[1]
ex.args[2].args[2]
typeof(ex.args[2].args[2])

We see that it has a special Symbol type, which represents the atoms (smallest parts) of an expression.

Modifying expressions

Now, what happens if we modify this? Let's make a copy of the expression first so that we can compare the two. Wht copy? Because we just created a different copy (deepcopy is not necessary!)

ex = :((x - 1) * (x - 2))
ex2 = copy(ex)

ex2.args[2] = :z
ex2
ex
ex2.args[2].args[2]
ex2.args[2].args[2] = :z

We have changed something inside the object ex2, so let's look at it:

ex2
ex

The original expression has, indeed, changed! That is, we have taken a piece of Julia code, and used Julia itself to manipulate it into a different piece of Julia code. This is one of the simplest examples of metaprogramming. Now we can define x and z and evaluate the expressions:

x = 3.5; z = 4.5
eval(ex), eval(ex2)  ## And of course they will be different

Walking a syntax tree

What if we wanted to replace all of the xs in a given expression?

The problem is that they may be buried arbitrarily deeply in the nested syntax tree. So we have to traverse the tree, in order to visit each piece of it and check if it is an :x. So I need to have an algorithm which will go to each nodes of the tree and checks whether this is an x if it is, then it is going to change it to z. This is maybe possible, right? Important thing is I have to visit each node in some way and check whether I have x

Here is a function that takes an expression object and replaces all of the :xs by :zs.

function traverse!(ex::Expr, find, replace)
    for (i, arg) in enumerate(ex.args)
        
        if arg == find
            ex.args[i] = replace   # we cannot use arg here, since it is a copy, not a reference
        end
        
        if isa(arg, Expr)
            traverse!(arg, find, replace)   # recursive
        end
    end
    
    ex
end

It's calling the function recursively and replacing. Typical example of a recursion

function recur_factorial(n)
   if n == 1
       return n
   else
       return n*recur_factorial(n-1)
    end

end

recur_factorial(10)
10*9*8*7*6*5*4*3*2


traverse!( :( (x-1) * (x-2) ), :(x), :(z) )

This is, of course, not general enough - for example, we cannot replace something of the form :(x - a) with :(x - 2a) with the current version; this would require more sophisticated pattern matching capabilities - but it demonstrates the basic idea.

Code generation

Let's return to the original question: how to create a long polynomial expression. So far, we have not seen how to add code, only change code that is already there.

A natural idea is to build the code up step by step; this is known as code generation.

Let's start again from the first term:

ex = :(x-1)
typeof(ex)
ex_new = :(ex * (x-2)) ## not correct

We would like to take what we have and create code that is "what we have multiplied by :(x-2)". Let's try. This doesn't work, since ex is treated as a symbol, whereas we need the value contained in the variable called ex. This is obtained using the $ operator, a procedure called interpolation. (Compare this to string interpolation.)

name = "David"
s = "Hello $name"  ## string interpolation

ex = :( $ex * (x-2) )

# Now we can continue:
ex = :( $ex * (x-3) )

Finally we see how to construct our loop:

n = 10
ex = :(x-1)

for i in 2:n
    ex = :( $ex * (x - i) )  
end

julia> ex
:((((((((((x - 1) * (x - i)) * (x - i)) * (x - i)) * (x - i)) * (x - i)) * (x - i)) * (x - i)) * (x - i)) * (x - i)) 

## So I get x-i everytime.

This did not work, since once again we did not want "the code 'i'", but rather the value of the variable i. Now let's try this one,

n = 10

ex = :(x-1)

for i in 2:n
    ex = :( $ex * (x - $i) )
end

ex

Pretty neat (Not because I did it, the way it shows!). This is almost what we would write by hand, except for the large number of parentheses. Now let's do n=100n = 100

n = 100

ex = :(x-1)

for i in 2:n
    ex = :( $ex * (x - $i) )
end

ex

Now we have the right hand side of the function, but we need the left hand side, which is the name of the function, so I can have the name of the function as,

n = 10
name = Symbol("p_", n) ## note if you give string, it keeps string, if you give variable it evaluates

The code is then

code = :( $name(x) = $ex )

We can wrap this inside a function,

function make_wilkinson(n)
    
    ex = :(x-1)

    for i in 2:n
        ex = :( $ex * (x - $i) )
    end
    
    name = Symbol("p_", n) 
    code = :( $name(x) = $ex )
    
    eval(code)  ## this is a bit dirty, using eval, eval is always in gloabal... that is why it is 
    ## discouraged
end 

x = 2
function check_eval()
     eval(:(y = x))
end  

check_eval()

Finally we evaluate this:

make_wilkinson(10)  ##  it makes the function p_10
make_wilkinson(100)

p_10(4)
p_100(2)

This creates a function with the name p_10 that does what we would write by hand.

Let's compare the two options by evaluating the function on a grid of values:

function f1(range)
    total = 0.0
    for x in range
        total += p_100(x)
    end
    return total
end

function f2(range)
    total = 0.0
    for x in range
        total += wilkinson(100, x)
    end
    return total
end

range = -10:0.000001:10
t1 = @elapsed f1(range);
t2 = @elapsed f2(range);
t2 / t1  ## around 20% gain in performance

We see that the generated code with the unrolled loop is 10% faster than the naive loop with 100 terms.

Starting from empty

In some cases, it is useful to start from something empty and add statements to it. To do this, we can do, for example,

code = quote end   # empty
dump(code)

We can add statements using the standard Julia push!:

Base.push!(code.args, :(a = 1))
Base.push!(code.args, :(b = a + 1))
code

eval(code)

dump(:(a * b * c))

Expr 
  head: Symbol call
  args: Array(Any,(4,))
    1: Symbol *
    2: Symbol a
    3: Symbol b
    4: Symbol c
  typ: Any

Exercise: Build up the Wilkinson example using this technique.

Repetitive code

Code generation is used frequently in Julia code when repetitive code is required. For example, let's return to the idea of wrapping a type:

struct OurFloat
    x::Float64
end

We can generate objects of this type:

a = OurFloat(3)
b = OurFloat(4)

But arithmetic operations are not defined:

a + b

We can define them in the natural way:

import Base: +, -, *, /
+(a::OurFloat, b::OurFloat) = a.x + b.x
-(a::OurFloat, b::OurFloat) = a.x - b.x

But this will quickly get dull, and we could easily make a mistake. As usual, whenever we are repeating something more than twice, we should try to automate it.

We have some code of the form

*op*(a::OurFloat, b::OurFloat) = a.x *op* b.x

where *op* is supposed to represent the operator. Julia allows us to do this almost literally; we just need to substitute in the value of the variable op!

@show x

## prints, name = value

Following is just for checking

for op in (:+, :-, :*, :/)
    @show op  ## name 
    code = :( $(op)(a::OurFloat, b::OurFloat) = $(op)(a.x, b.x) )
    @show code
end

op = :+
code = :(a::OurFloat + b::OurFloat = begin  # In[127], line 3:
            a.x + b.x
        end)
op = :-
code = :(a::OurFloat - b::OurFloat = begin  # In[127], line 3:
            a.x - b.x
        end)
op = :*
code = :(a::OurFloat * b::OurFloat = begin  # In[127], line 3:
            a.x * b.x
        end)
op = :/
code = :(a::OurFloat / b::OurFloat = begin  # In[127], line 3:
            a.x / b.x
        end)

Finally we need to evaluate the code. The combination of eval and :(...) that we used above can be abbreviated to @eval:

for op in (:+, :-, :*, :/)
    @eval $(op)(a::OurFloat, b::OurFloat) = $(op)(a.x, b.x) 
end

a*b
@which a*b

@which sin(3.5)