Introspecting into a function?

New to Julia
11

I’m actually not using a symbolic package, I wanted total control over evaluation and simplification so I rolled my own. Its still in the early stages but I’ll share some snippets here to show you whats going on but I think I’ll put the whole thing up on a github repo soon.

The idea is that I’ve added methods to all the basic mathematical functions such that

julia> :x + :y 
:(x + y)
julia> :x - 5 
:(x + -5)   #I'm not sure I'm super happy with this choice but it makes simplification much easier.
julia> sin(:x)
:(sin(x))

etc.

And then I made a simplification algorithm which looks for common mathematical patterns in the expression tree and applies replacement rules to them and then keeps on looping until the expression stops changing ie.

julia> :(x^1) |> expansion_loop
x
julia> :(0 * x) |> expansion_loop
0
julia> :(1 * x) |> expansion_loop
x
julia> :(0 + x) |> expansion_loop
x

etc.

Finally, to take derivatives of symbolic functions I made a struct

Mathy = Union{Number, Symbol, Expr}

struct Dual_Number
    real::Mathy
    infinitesimal::Mathy
end

function dual_number(a::Mathy, b::Mathy)
    if b == 0
        a
    else
        Dual_Number(a, b)
    end
end

dual_number(a) = dual_number(0, a)

and then screwed around with how dual numbers are displayed, ie.

function Base.show(io::IO, dual::Dual_Number)
    if (dual.infinitesimal isa Number) && (dual.infinitesimal < 0)
        op = "-"
    else
        op = "+"
    end
    real_string = dual.real == 0 ? "" : "$(dual.real) $op "
    infinitesimal_string = dual.infinitesimal == 0 ? "" :
                           abs(dual.infinitesimal) == 1 ? "ϵ" :
                           (dual.infinitesimal isa Expr) && (dual.infinitesimal.args[1] == :+ || dual.infinitesimal.args[1] == :-) ? "($(dual.infinitesimal))ϵ" :
                           "$(dual.infinitesimal)ϵ"
    print(io, "$real_string$infinitesimal_string")
end

Next, I’ve made methods for some standard mathematical functions for working with dual numbers that are agnostic about numbers versus symbols and expressions (because I can compose all of them). The idea here is that if I define


conj(a::Dual_Number) = dual_number(a.real, -(a.infinitesimal))

+(a::Dual_Number, b::Dual_Number) = dual_number(a.real + b.real, a.infinitesimal + b.infinitesimal)
+(a::Dual_Number, b::Mathy) = dual_number(a.real + b, a.infinitesimal)
+(a::Mathy, b::Dual_Number) = dual_number(a + b.real, b.infinitesimal)

-(a::Dual_Number, b::Dual_Number) = dual_number(a.real - b.real, a.infinitesimal - b.infinitesimal)
-(a::Dual_Number, b::Mathy) = dual_number(a.real - b, a.infinitesimal)
-(a::Mathy, b::Dual_Number) = dual_number(a - b.real, b.infinitesimal)

*(a::Dual_Number, b::Dual_Number) = dual_number(a.real * b.real, a.infinitesimal * b.real + a.real * b.infinitesimal)
*(a::Dual_Number, b::Mathy) = dual_number(a.real * b, a.infinitesimal * b)
*(a::Mathy, b::Dual_Number) = dual_number(a * b.real, a * b.infinitesimal)

/(a::Dual_Number, b::Dual_Number) = b.real != 0 ? (a * conj(b))/(b.real)^2 : Inf
/(a::Dual_Number, b::Mathy) = dual_number(a.real / b, a.infinitesimal / b)
/(a::Mathy, b::Dual_Number) = b.real != 0 ? (a * conj(b))/(b.real)^2 : Inf

^(a::Dual_Number, b::Mathy) = dual_number(a.real^b, b * a.real^(b-1) * a.infinitesimal)
^(a::Dual_Number, b::Integer) = dual_number(a.real^b, b * a.real^(b-1) * a.infinitesimal)
^(a::Mathy, b::Dual_Number) = dual_number(b^a.real, log(b) * a^b.real * a.infinitesimal)

log(a::Dual_Number) = dual_number(log(a.real), 1/a.real * a.infinitesimal)

sin(a::Dual_Number) = dual_number(sin(a.real), cos(a.real)*a.infinitesimal)
cos(a::Dual_Number) = dual_number(cos(a.real), -sin(a.real)*a.infinitesimal)

The idea here is that if I define ϵ = dual_number(0, 1) then ϵ^2 == 0 so ϵ is like a very small (infinitesimal) number and for any function f(x), f(x + ϵ) = f(x) + ϵ f'(x)

I then made a derivative function which takes a function f and returns a function which when evaluated will evaluate f(t + ϵ) and extracts out the infinitesimal part which will be the automatic derivative of f and then simplifies that result.

ϵ = dual_number(0, 1)
function D(f::Function)
   t::Mathy ->  f(t + ϵ).infinitesimal |> expansion_loop
end

With this I can do stuff like

julia> f(x) = log(sin((x^2)));
julia> D(f)(:x)
:((1 / sin(x ^ 2)) * (cos(x ^ 2) * (2x)))

which is correct(!!!) but could be simplified a bit more.

Sorry that turned into a lot of code but I was excited enough that I wanted to share what I had done. Once this is a bit more polished I’ll maybe make some blog posts about it if people are interested.

3 Likes