Constructing a function of functions on the fly - approaches and best practices?

New to Julia
1

Think of something like R’s fivenum function (that returns the min, max, median, and upper/lower quartiles of a data vector). Now, instead of calculating all five numbers, say I wanted to be able to tell the function to just give me a particular subset of those statistics.

Is there a preferred or “Julian” way of handling stuff like this? The obvious (to me) thing to do would be to pass an argument to the function saying which statistics to calculate, and using a bunch of if…ifelse statements, checking each statistic against the arguments to determine whether to calculate it.

However, I’m trying to make this performant, and wasn’t sure whether a whole mess of if statements was a good idea. If nothing else, it makes for a very long and messy function.

An alternative that came to mind was something like the following (I want the results in a vector):

a(x) = [minimum(x), maximum(x), median(x)]

That works fine for me, and it keeps the code nice and clean, but it means hand-coding/constructing the function, and makes it harder for other users, who might be less comfortable writing code. It’s not the end of the world, but ideally, it’d be nice to do something like:

stats = ["maximum", "minimum"]
results = fivenum(x, stats)

especially if the resulting function is still performant.

To do something like that, I am stuck using a whole bunch of if statements? Is there a better way to (essentially) generate and call an array of expressions on a given input? In my use case, the “function of functions” is going to be defined once at runtime, and called many, many times as part of a larger MCMC algorithm.

I’m just learning Julia, so apologies if my questions are malformed, or I’m missing something obvious (I know R and a bit of C++). I’m guessing this is the sort of thing macros might be used for, but I frankly have trouble wrapping my head around macros in Julia…

Thanks!

2 
> stats = [maximum, minimum];
> fivenum(x, stats) = [f(x) for f in stats]
> results = fivenum([1,2,3,10], stats)
[10, 1]

?

1 Like
3

You can have the user pass you a function and set a default. Functions can be passed and used like variables.

You can also setup a type whose type-information denotes which analyses to run, and use the dispatch on booleans in its parameters. That would be a good use of “hyperpure” e.g. Base.@pure.

4

cstjean:

That syntax seems pretty nice and clean.

Some of my functions take multiple arguments (similar to the quantile function, for example); is the best approach to “pad” existing ones so every function takes the same arguments? E.g.,

maximum2args(x, y) = maximum(x)
minimum2args(x, y) = minimum(x)
stats = [maximum2args, minimum2args, quantile];
fivenum(x, y, stats) = [f(x, y) for f in stats]
results = fivenum([1,2,3,10], 0.1, stats)

Bit messier, but seems to work well enough, and my existing functions would only require a minimal change.

Now, some of my functions return a vector, rather than a scalar. In that case, the generator above gives me an “Any” vector, but the code below seems to work well enough:

(Note that I changed the “maximum” function so it returns a vector, allowing me to test that particular use case.)

maximum2args(x, y) = fill(maximum(x), 2)
minimum2args(x, y) = minimum(x)
function funclist(statlist, x, y)
  tempvec = fill(0.0, 0)
  for e in statlist
    append!(tempvec, e(x, y))
  end
  return tempvec
end

stats = [maximum2args, minimum2args, quantile]
results = funclist(stats, [1.0, 2.0, 3.0, 10.0], 0.1)

That’s not too complex, and I think it handles the use cases I’m concerned about right now.

ChrisRackauckas:

I think I get the gist what you’re saying, but my understanding of types is still really weak – right now, I don’t use them for much beyond defining types of variables that go into functions, and constructing some custom data structures (e.g., via composite types). Do you happen to have a link or example of the general sort of thing you’re talking about? Also, what is @pure? I can’t seem to find clear information about what it is/does.

5

@opera_malenky are you really sure you need this? I think that your very first idea:

a(x) = [minimum(x), maximum(x), median(x)]

is just fine. I think about it this way:

Let’s say a new user wants the maximum and minimum of their data. In the simple implementation, they have to know the name of the functions “maximum” and “minimum”, and they compute their statistics with:

mydata = [maximum(x), minimum(x)]

in the fivenum world, the user doesn’t have to know the function names, but they still have to know the names of the statistics (which are, actually, just the names of functions) plus they have to know how fivenum works (which isn’t something they will have encountered anywhere else in Julia) in order to do:

mydata = fiivenum(x, ["maximum", "minimum"])

which doesn’t actually save them from remembering anything at all. Plus, it requires you to create this complicated implementation and maintain it forever.

In essence, what you’re doing is creating a very very very simple programming language inside of Julia, in which users can pass a vector of function names and get a vector of results. But Julia’s already a pretty good language, so you might be better off just letting your users use Julia!

Also, your last example looks like the beginnings of an implementation of broadcasting. Fortunately, Julia has broadcasting built-in, and it’s amazing. For example, to do f(x, y) for each x in [1, 2, 3] and y = 5, you just define the scalar function:

f(x, y) = x + y

and call it with an extra dot:

f.([1, 2, 3], 5)
3 Likes
6

@rdeits

I appreciate the advice/comments.

I ran into some issues with the “a(x) = …” style when I wanted to add multiple function parameters and functions that return a vector, but (as you might expect) that was just PEBCAK. The following seems to work pretty well:

stats2(x, y) = vcat(fill(maximum(x), 2), minimum(x), quantile(x, y)
b = randn(10)
results = stats2(b, 0.5)

That’s more than good enough for my personal use, tbh.

Still, part of my reason for posting was to see what might be required to make my functions a bit more user-friendly down the line, and to get a little more comfortable using/manipulating/generating functions.

The reason I’m interested in that (aside form the learning part) is that in my actual algorithm, each statistic/function actually has a counterpart function. For every statistic T, there are actually two associated functions, T_f(x) and T_g(x). So, a model with statistics T and S actually requires two arrays of functions, [T_f(x), S_f(x)] and [T_g(x), S_g(x)].

Now, using the notation from the start of this post, it’s pretty easy to write something like:

f(x, y) = vcat(T_f(x), S_f(x))
g(x, y) = vcat(T_g(x), S_g(x))

Again, that would work perfectly well for me, and it’s much better than the way I have things hard-coded right now.

For future reference, though, I am curious to see if anyone has ideas about how one could take an array of functions like f(x) = [A(x), B(x), C(x, y), C(x, z)], and use that to generate some array of associated functions, e.g., g(x) = [A1(x, i), B1(x, i), C1(x, y, i), C1(x, z, i)].

(You can see, perhaps, why I mentioned the idea of using strings in my original post; as fiddly and annoying as they are, I can sort of wrap my head around using them for something like this…)

RE broadcasting: Terrific feature, definitely loving that so far. Array generators in Julia have also been fantastically useful in my work. TBH, the hardest part so far (that I still struggle with) has been figuring out when things are copied vs. passed by reference. I find myself getting paranoid sometimes and throwing deepcopy()'s around everywhere. :frowning:

7

Hm, I think I don’t understand what you mean by the associated functions. Do you have an example handy?

And you probably don’t need all those deepcopy()s :wink: The rules should hopefully be pretty simple: all function arguments are passed by reference, period. Immutable types might be copied when passed to a function, but since they’re immutable, there’s no difference between a copy and a reference. So the statement “all function arguments are passed by reference” is indistinguishable from the truth. The most common thing that quietly makes a copy is a slice a[1:3].

8

In this case, the “associated functions” I described are designed to calculate how the vector of selected statistics changes when one specific element of the data is mutated in some way (for example, increased or decreased). For example, if one of my statistics was mean(x), then I’d have a corresponding function Δmean(x, y, i) in my algorithm that tells me how much the mean changes if I were to change the value of x[i] to y.

@ChrisRackauckas mentioned trying to use types. No idea if I’m doing it right, but I took a stab at it, and the code below seems to do what I want. :slight_smile: In any case, it probably gives a clearer idea of what I was trying to do:

# Define a compound function that includes the statistic of 
# interest and it's associated "mutation"
type data_statistic
  statistic::Function
  Δstatistic::Function
end

# Replace x[i] with value y, and see how the mean changes
function Δmean(x::Array{Float64, 1}, y::Float64, i::Int64)
  return (x[i] - y)/size(x, 1)
end

# Replace x[i] with value y, and see how the variance changes
function Δvar(x::Array{Float64, 1}, y::Float64, i::Int64)
  z = copy(x)
  z[i] = y
  return (var(x) - var(z))
end

# Create the compound functions for the selected statistics
mean_statistic = data_statistic(mean, Δmean)
var_statistic = data_statistic(var, Δvar)

# Define an array of compound functions
statlist = [mean_statistic, var_statistic]

# Define functions to get numbers, given the list of compound functions
statistic_statlist(dat) = map(x -> x.statistic(dat), statlist)
delta_statlist(dat, y, i) = map(x -> x.Δstatistic(dat, y, i), statlist)

# Test
x = rand(10)
y = 10.0
i = 1
means = statistic_statlist(x)
deltas = delta_statlist(x, y, i)

I’m using the mean and variance here as minimal examples – my actual functions are non-standard and much more computationally intensive.

9

Cool! Thanks for the example!

Are you familiar with automatic differentiation like https://github.com/JuliaDiff/ForwardDiff.jl ? Rather than writing out your Δvar function by hand, you might be able to use ForwardDiff.derivative to compute the derivative of the variance. Forward-mode automatic differentiation like this will generally give the exact answer and is crazy fast and easy to use.

3 Likes