Transforming a function into a sequence of mutable operations on objects

The following example comes up in the Taylor method for integrating differential equations, but the question keeps coming up over and over again.

Suppose I have something like the following, which is over-simplified, in particular with the wrong formula for * for simplicity, but contains the kernel of the question:

struct Taylor
    coeffs::Vector{Float64}
end 

add!(z::Taylor, x::Taylor, y::Taylor, n) = z.coeffs[n] = x.coeffs[n] + y.coeffs[n]
mul!(z::Taylor, x::Taylor, y::Taylor, n) = z.coeffs[n] = x.coeffs[n] * y.coeffs[n]

x = Taylor(zeros(5))
y = Taylor(zeros(5))

Now if a user specifies a polynomial function like
f(x, y) = x*x + y,
we need to transform it into code that produces finally a new Taylor object, but in order to do so, it must calculate the coefficients of the new object element by element.

Currently in TaylorIntegration.jl, there is a macro @taylorize that does this by direct and painful code generation.
My question is if there is a better method, e.g. using operator overloading or creating an Intermediate Representation, e.g. using ModelingToolkit.jl, and doing something clever with that?

z1 = copy(x)
z2 = copy(y)
for n in 1:length(coeffs.x)
    mul!(z1, x, x, n)
    add!(z2, z1, y, n)
end

We really do need to do this element by element because of how the Taylor method works. Sowe need some way of persistently recording which intermediate objects get created along the way and modifying those later.

I have a rather different approach outlined in https://github.com/dpsanders/LazyTaylorSeries.jl, but I think it would be very useful to see any suggestions that people can make about this kind of coding pattern.

cc @lbenet

I am certainly not a julia expert, but maybe using an aproach like forwardDiff can work, they have a struct that stores a number <:Real and the derivative, and they overload the base operations with this struct:
(the following code is a modification of the dual number implementation in forwardDiff.jl

struct TaylorNumber{V,N<:Real} <: Real
number::N
coeffs::Vector{V}

and then import the base operators:

import Base: *
function *(x::taylorNumber,y::taylorNumber)
return taylorNumber(x.value+y.value,taylorMul(x.coeff,y.coeff)
end 

doing this, using a function of the form f(x,y) where x,y are Real, can automatically generate such coefficients, please correct me if am wrong.

I totally recomend reading that source code, they know what they are doing, i am not

Thanks, I am aware of ForwardDiff.jl; actually I think I might need something more like ReverseDiff.jl, which builds up a tape of the operations that you can play back later. (Or maybe even Cassette.jl…)

Unfortunately the difficulty here is that I need to refer back to the same intermediate result later to modify it.

This sounds a bit like the LazyExpressions from Parametron.jl: Parametron.jl/lazyexpression.jl at master · tkoolen/Parametron.jl · GitHub . I wonder if there is some code (or at least some ideas) you could adapt from there.

2 Likes

That looks very interesting - I’ll check it out, thanks!

This also looks similar to the macro rewritting done in JuMP in parse_expr.jl.
When you write @constraint(model, x + u + y * z), it rewrite the expression to JuMP.destructive_add!(JuMP.destructive_add!(copy(x), u), y, z).
We plan to make this rewriting independent of JuMP as a GSoC project as this would allow any object used in a @constraint expression to take advantage of the rewritting (I have a particular application in minds which is polynomials for Sum-of-Squares).
But the package would also have a, say @expression, macro and if your object is a subtype of, say AbstractMutable, and implements the mutable API, it could be constructed with @expression and we can share the rewritting code and not duplicate the work in JuMP, Parametron, LazyTaylorSeries, …

1 Like