Extract intermediate data from inside a function

General Usage
1

Hi,
here comes a question which in any other language would be wishful, yet in Julia, I dare wish.

So I program a function
R = f(X)
think finite elements 101, given nodal displacements, the function returns nodal forces.

Of course, there are many intermediate results, say strain, and stress, not needed outside of my function f, in the logic of the FEA solver, but valuable for debugging and not least, for reporting.

My strategy is to not store intermediate results systematically during the analysis. This gobbles memory and CPU, and fills the code with boilerplate. What I want to do is store all the X (for all elements and steps…), and if I want, call f after the analysis is complete. And now the crux: while adding as little boilerplate as possible, I want to be able to “call f(X)” and add “give me the stresses”.

How would I do that?

  • store intermediate results in a struct, return the struct, and the solver just throws it away
  • metaprogramming? Transform readable code for function f into code with boilerplate for the reporting of intermediate results?
  • …???

With what strategy would you attack the problem? I am not asking for code, but pointers to the relevant language features that could be exploited here.

: )

Philippe

1 Like
2

Similar to the eigs function, you could return a tuple of results, with the most commonly needed results first, in the expectation that most callers will ignore the subsequent results. For example, x, = f(...) discards all results except for the first.

Or you could pass some optional keyword argument, with default value nothing—if the caller passes an empty array (or some other data structure with mutable contents) instead, it is mutated by your function to hold the requested data.

1 Like
3

Use a thunk?

function f(x)
    y = 2*x
    calculate_extra = () -> 3*y
    y, calculate_extra
end

# just the result
y, = f(1)

# extras calculated on demand
y, g = f(1)
g()
5 Likes
4

looks like the calculation happens regardless, so it would be more accurate to say assigned instead of calculated on demand, since the calculations are not affected

julia> function f(x)
           y = 2*x
           calculate_extra = (println("calculated some extra"); () -> 3*y)
           y, calculate_extra
       end
f (generic function with 1 method)

julia> y, = f(1)
calculated some extra
(2, getfield(Main, Symbol("##9#10")){Int64}(2))
5

just write it as

           calculate_extra = () -> (println("calculated some extra"); 3*y)

The calculations are indeed affected.

3 Likes
6

Depends on what you are trying to do, my point is that calculate_some extra is created either way, since it is an anonymous function you can defer some calculation for later if you want to, but this function is created regardless also if you don’t use it. I consider the creation of it as a calculation.

7

Do extra calculations on demand?

8

I’ve had a similar problem, and solved it by having the function return a tuple consisting of the solution and a struct with all intermediates. I then had a Bool flag passed to the function, and only calculated/populated the intermediates if the flag was true. (Otherwise, I left those fields empty.)

I’m not quite sure I follow the suggestion of calling f(X) a second time to get the intermediates, how would that avoid boilerplate code? Are you imagining a separate code path to calculate those?

Boilerplate code that makes the code less readable is of course an issue. A bigger issue, IMO, is repeated code. In a situation like this, I’d first try to reduce boilerplate code by putting the code that calculates the intermediates in a separate function, then put something like this somewhere in the main algorithm:

calculate_intermediates && calculate_intermediates!(X)

Hard to give more specific advice without seeing your code, perhaps your algorithm looks very different in the case that the intermediates are calculated?

9

Another possibility is to have a callback function that you call on the intermediate data, which can then decide whether to save it or do something else. For example, consider the following function that iteratively computes a recurrence relation (e.g. a discretized ODE or PDE):

function iterate_recurrence!(usercallback::Function, x, N)
    usercallback(x, 0) # initial data
    for n = 1:N
        some_recurrence!(x, n) # mutate x somehow, in-place
        usercallback(x, n)
    end
    return x # return final result only
end

# define a default callback that does nothing
iterate_recurrence!(x, N) = iterate_recurrence!((x,n)->nothing, x, N)

The default usercallback does nothing, and only the final result of the recurrence is returned. You certainly wouldn’t want to store all of the intermediate x values by default just in case the user might want them, because this could incur a huge storage cost (x could be a huge vector and N could be huge).

The advantage of a functional style (a usercallback function) is that then the caller has complete freedom to decide what to do. For example, you can use a callback to save every 10th x:

function doit(x, N)
    xsaved = [copy(x)] # array to store saved x's, initialized with first x
    iterate_recurrence!(x, N) do xâ‚™, n
        if n % 10 == 0 && n > 0
            push!(xsaved, copy(xâ‚™))
        end
    end
    return x, xsaved
end

Note that putting the callback as the first argument allowed me to use Julia’s nice do syntax. Also note that I pushed copy(xₙ), not xₙ, to the array xsaved because iterate_recurrence! mutates x in-place—I want to save N÷10 distinct vectors, not N÷10 references to the same vector.

You could imagine lots of other variations. For example, the callback could save only one component x[1] of the data. Or it could update a plot of the data. Or it could write data to the disk. It would not be possible to provide all possible behaviors efficiently just by passing flags to iterate_recurrence!, but a callback interface offers this flexibility with ease.

A much more sophisticated variation on this idea is implemented by the DifferentialEquations.jl package, which provides a wide variety of event-driven callback interfaces.

10 Likes
10

One technique which I employed in SyntaxTree.jl is to have intermediate helper functions which get called from a meta-function for organizing the shared allocations.

https://github.com/chakravala/SyntaxTree.jl/blob/master/src/exprval.jl

In this case, the expravg function computes some intermediate values needed in several later calculations:

function exprval(expr)
    val = expravg(expr)
    cal = callcount(expr)
    mal = sqrt(exprdev(expr,val[2],cal))
    cal*sqrt(abs(val[2])*mal)*val[4], cal, mal, val[2], val[4]
end

Then all the different values from the helper functions are combined in the meta-function. However, I only export the final meta-function and not the helper functions.

What’s special about this is that the helper function called expravg is both capable of recursively allocating itself to tally up calculations from an abstract syntactic tree and it can also be used for obtaining the intermediate calculations in the meta-function, thus splitting off that functionality neatly into a separate recursive function which gets called from and terminates at the meta-function needing intermediate result.

11

I might have thought that the multiplication would be optimized away. I don’t see how to change the value of y after the call returns.

So f(1) would return

2, () -> (println("calculated some extra");  6)

But, that’s not the case

julia> function f(x)
                  y = 2*x
                  calculate_extra = (println("calculated some extra"); () -> 3*y)
                  y, calculate_extra
              end
f (generic function with 1 method)

julia> z, g = f(1)
calculated some extra
(2, getfield(Main, Symbol("##5#6")){Int64}(2))

julia> @code_lowered g()
CodeInfo(
1 ─ %1 = (Core.getfield)(#self#, :y)
│   %2 = 3 * %1
└──      return %2
)
12

Hi Bennedich,

To explain my multiple pass idea: in the analysis, “f” is called for each element at each iteration of each step.

After the analysis, I want to allow the user to extract arbitrary “intermediaries”, interactively. But I want to avoid to store all intermediaries during the analysis, hence the extra calls to “f” “directly” by the user. So extra calls is not a strategy to solve the given problem, I make it part of my problem definition.

: )

Philippe

13

Neat.

Very neat. With, in addition, a little macro to hide the calls to the callback, a line in “f” could look like

function f(usercallback::Function,X, morestuff)
@remember intermediate_result = some_hack(X)
I will think about this very carefully, thank you!

: )

Philippe

14

…and a reply to all: I was right to ask, I got valuable input! In several cases, your answer points at features of Julia I am not well familiar with, so it will take me some time to understand in depth how you think.

But I am sure, that by the time I have done that, I will see one or more good solution to my question.

Thank you indeed!

1 Like