Evaluating expressions in a local scope

[Jonas Rauch]

Hi,
I've tried to get the following to work and did not succed:
function f()
    local x=1;
    E=:(x+1)
    eval(E)
    x=2
    eval(E)
end
where the first eval should return 2 and the second one should return 3.   What I get from the above is "x not defined". I read in the manual, that eval works on the global scope. How can I include local variables here? I tried
function f()
    local x=1;
    E=:($x+1)
    eval(E)
    x=2
    eval(E)
end
which seems to include the value of $x in the expression. This will of course return 2 both times. When I change x to be an array
function f()
    local x=[1, 2];
    E=:($x+1)
    eval(E)
    x=2
    eval(E)
end
calling f() gives me a segmentation fault. Returning the Expr seems to work, but evaluation does not:
julia> function f()
           local x=[1, 2];
           :($x+1)
       end

julia> E=f()
+([1, 2],1)

julia> eval(E)
Segmentation fault
I suspect the last one is a bug? In addition, is there any way to do what I described at the beginning, possibly with a let ... block? I tried different things but did not succeed.

[Stefan Karpinski]

Eval always evaluates in the top-level context, not the context where it's called. Splicing x into an expression using $ will splice the value of x into the expression. Expressions are not like closures: they don't capture local variables. The last bit is definitely a bug. You should never get a segfault. It seems to me like what you want is just a closure:

function f()

    x = 1

    g() = x+1

    println(g())

    x = 2

    println(g())

end

julia> f()

2

3

[Jeff Bezanson]

- Bug fixed.
- In julia eval is just for making top-level objects. A good technique
is to build a function expression with arguments, eval that, then call
it.
- There are supposed to be rules about what can go in an AST, but
they're not enforced.

[Jonas Rauch]

Thanks for the clarification.


I was trying to use expressions instead of closures, because I thought about how to approach symbolic differentiation, as suggested by Douglas Bates in the ODE thread.
It might be possible to define a derivative operator
D(E::Expr, x::Symbol) -> ::Expr,
that computes the derivative of an expression w.r.t. a symbol. This should be fairly straightforward for most standard numerical functions and operators. In order to use the resulting expression to actually evaluate the derivative, it would be necessary to evaluate it in a local scope (I think).

In addition, it would be nice to be able to get the body of a function as an expression and vice versa to construct a function from an expression. So for a function
f(args...)=begin
...
end
we could do something like
df=(args...) -> eval(D(body(f), args[1]))
to define a function that computes the derivative of f w.r.t. its first argument.
This concept could also be used to compute partial derivatives of expressions that are calls to functions that don't have a known derivative. Putting everything together, it should be possible to symbolically compute a derivative of almost anything that has one.

Douglas: What do you think? As far as I know, symbolic differentiation in R works similarly, right?

[Jeff Bezanson]

The overall thinking is right. But instead of evaluating an expression
in another function's scope, it's better to wrap it with arguments:

julia> ex = :(2*x)
*(2,x)

julia> fex = :(x->($ex))
->(x,*(2,x))

julia> f = eval(fex)
#<function>

julia> f(3)
6

[Jonas Rauch]

Thanks, that helped a lot. I just could not figure out that syntax :-)

I just finished a first attempt at symbolic differentiation (see attachment). This is a proof of concept kind of thing, I hope someone can expand it to something fully functional, because I probably won't have the time. But I would be interested to hear what you guys think. Does this approach make sense? Right now I think it only works on simple functions defined as
f(args...)=...
The basic idea is as follows:

To get the body of a function as an Expr, I defined operators (only the basic ones for now) for Expressions and Symbols, so that for example
+(x::Expr, y::Expr) = :($x+$y)
now I can call a simple function
h(x,y)=x+2*y+x^2*y
with Symbol arguments:
julia> E=h(:x, :y)
+(+(x,*(2,y)),*(*(x,x),y))
Defining a derivative is pretty straightforward now:
D(E::Symbol, x::Symbol)=(E==x) ? 1 : 0
D(E::Expr, x::Symbol) implements the standard rules for differentiation
D(E, x::Symbol) = 0  #constants have derivative zero

Then, with the syntax you gave me in your last post,
dhdx=@eval (x,y)->$(D(E, :x))
will yield a function dhdx(x, y) that computes the partial derivative of h(x,y) w.r.t. x

More details in the attached source file. I think that the derivative function constructed this way will be JIT compiled at first evaluation, right? This would mean, that it would be very quick. Plus, it does not matter so much that some expressions are suboptimal, because the compiler should get rid of things like a*1, a*0 or b+0.

What is missing now is of course support for the simple mathematical functions (sin, exp, ...) and more advanced expressions. For example, the derivative of a call to
max(args...) would have to become a quoted form of something like
begin
  V, i=findmax(args)
  D(args, x)[i]
end

I have no idea at the moment how to handle assignments. Basically, for something like
a=...
we would need another line
Da=D(..., x)
so that later in an expression that involves a we can use Da for the chain rule and product rule. In particular, D(:a, :x)=Da  and not zero (as for other symbols). Maybe a new type, that contains an Expr and the derivative?

Sorry for the long post. I hope someone can make this work nicely (or propose something totally different!)  :-)

[Jonas Rauch]

Short update: the way it is now, assignments in a function body actually will just work on the expressions. What comes out of f(:x, :y) at the end will not have any assignments any more.

However, this amounts to rewriting the function in a purely functional way. For more complex algorithms including lots of temp variables and conditional assignments, transforming everything to an Expr would make things extremely complicated. So in this case a "code generation" kind of automatic differentiation would make sense, where the algorithm is changed so that for every variable it stores the derivative as well - similar to what I wrote about a and Da above. This would really require getting the body of a function to be able to manipulate it.

I think an optimal result would be a smart mix of symbolic differentiation and code generation, where loops and conditionals are treated via code generation and everything else is done symbolically. E.g. for
for i=1:n
x = sin(x*x)
end
the derivative of the RHS expression in the assignment can be computed symbolically to be cos(x*x)*2*x*x' and the loop can be written as
for i=1:n
(x, x') = (sin(x*x), cos(x*x)*2*x*x')
end
where of course x' has to be initialized in a sensible way. This should often be faster than my forward automatic differentiation, because in most cases x and x' will be (Arrays of) a bitstype, while the AD code I'm experimenting relies on composite types for every variable.

[david tweed]

Hi, I'm still getting to grips with Julia and haven't even looked at
which things are done by "julia itself" and which are pushed down into
the LLVM.

On Mar 19, 9:42 am, Jonas Rauch <jonas.ra...@googlemail.com> wrote:
> right? This would mean, that it would be very quick. Plus, it does not
> matter so much that some expressions are suboptimal, because the compiler
> should get rid of things like a*1, a*0 or b+0.

For what it's worth, as someone working with automatic differentiation
on a different project that goes via LLVM, expressions of various
_floating point_ types like a*0.0 don't get optimised away by LLVM
without using unsafe-math. (This is because things like this aren't
valid in cases where a or b are "special values" such as NaN.) I
looked at the possibility of writing LLVM code to apply those rules,
but personally decided to do the transfomation before lowering LLVM to
LLVM because it was clearer which 0's, 1's, etc, had arisen from
autodifferentiation as opposed to other constants in the code. (This
is because I'm trying to use the "mostly true" rule that applying
these simplifications will only produce a different derivative in
cases where the original function will contain a special value (NaN,
Inf, +0.0 vs -0.0, etc) so that you can be "mostly certain" that
everything's ok if you both check the original function evaluation and
the computed derivative, and it does lead to a very significant
reduction in "pointless" FP instructions.

Of course, it may be that Julia is applying those transformations
higher up, in which case this issue doesn't even occur.

Anyway, hope this helps,
David Tweed

> > On Sun, Mar 18, 2012 at 2:43 PM, Jonas Rauch <jonas.ra...@googlemail.com>

> > jonas.ra...@googlemail.com>

> > jonas.ra...@googlemail.com>

>  sd.jl
> 2KViewDownload

[Stefan Karpinski]

There are no fancy transformations applied at higher levels. We rely entirely on LLVM's floating point semantics, which are strictly IEEE and safe.