Category Archives: Julia

Glue AD for Full Language Differentiable Programming

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/glue-ad-for-full-language-differentiable-programming/

No design choice will be the best choice for all possible users. That’s a statement that is provocative but at the same time I think everyone would easily agree with it. But that should make us all question whether it’s a good idea to ever try and make all users happy with one piece of code. Under the differentiable programming mindset we are trying to make all code in the entire programming language be differentiable, but why would we think that a single system with a single set of rules and assumptions would be the best for everyone?

Optimized Algorithms Across Scientific Computing and Machine Learning

Differentiable programming is a subset of modeling where you model with a program where each of the steps are differentiable, for the purpose of being able to find the correct program with parameter fitting using said derivatives. Just like any modeling domain, different problems have different code styles which must be optimized in different ways. Traditional scientific computing code makes use of mutable buffers writing out nonlinear scalar operations and avoid memory allocations in order to keep top performance. On the other hand, many machine learning libraries allocate a ton of temporary arrays due to using out of place matrix multiplications, [which is fine because dense linear algebra costs grow much faster than the costs of the allocations](https://www.stochasticlifestyle.com/when-do-micro-optimizations-matter-in-scientific-computing/). Some need sparsity everywhere, others just need to fuse and build the fastest dense kernels possible. Some algorithms do great on GPUs, while some do not. This intersection between scientific computing and machine learning, i.e. scientific machine learning and other applications of differentiable programing, is too large of a domain for one approach to make anyone happy. And if an AD system is unable to reach top notch performance for a specific subdomain, it’s simply better for the hardcore package author to not use the AD system and instead write their own adjoints.

Even worse is the fact that mathematically there are many cases where you should write your own adjoints, since differentiating through the code is very suboptimal. Any iterative algorithm is of this sort, where the derivative of a nonlinear solve f(x)=0 may use Newton’s method to get f(x*)=0, but the adjoint is only defined at x* with f'(x*), so there’s no need to ever differentiate through Newton’s method. So we should all be writing adjoints! Does this mean that the story of differentiable programming is destroyed? Is it just always better to not do differentiable programming, so any hardcore library writer will ignore it?

Is There A Common Ground?

Instead of falling into that despair, let’s just follow down that road with a positive light. Let’s assume that the best way to do differentiable programming is to write adjoints on library code. Then what’s the purpose of a differentiable programming system? It’s to help your adjoints get written and be useful. It’s just matrix multiplications in machine learning. If the majority of the code is in some optimized kernel, then you don’t need to worry about the performance of many other aspects rest: you just want it to work. With differentiable programming, if 99% of the computation is in the DiffEq/NLsolve/FFTW/etc. adjoints, what we need from a differentiable programming system is something that will get the rest of the adjoint done and be very easy to make correct. The way to facilitate this kind of workflow would be for the differentiable programming system to:

  1. Have very high coverage of the language. Sacrifice some speed if it needs to, that’s okay, because if 99% of the compute time is in my adjoint, then I don’t want that 1% to be hard to develop. It should just work, however it works.
  2. Be easy to debug and profile. Stacktraces should point to real code. Profiles should point to real lines of code. Errors should be legible.
  3. Have a language-wide system for defining adjoints. We can’t have walled gardens if the way to “get good” is to have adjoints for everything: we need everyone to plug in and distribute the work. Not to just the developers of one framework, and not just to the users of one framework, but to every scientific developer in the entire programming language.
  4. Make it easy to swap out AD systems. More constrained systems may be more optimized, and if I don’t want to define an adjoint, at least I can fallback to something that (a) works on my code and (b) matches its assumptions.

Thus what I think we’re looking for is not one differentiable programming system that is the best in all aspects, but instead we’re looking for a differentiable programming system that can glue together everything that’s out there. “Differentiate all of the things, but also tell me how to do things better”. We’re looking for a glue AD.

If that’s where we need to go, how do we get there?

Zygote is surprisingly close to being this perfect glue AD. It’s stacktraces and profiling are fairly good because they point to the pieces generating the backpasses. It just needs some focus on this goal if it wants to obtain it. For (1), it would need to get higher coverage of the language, focusing on its expanse moreso than doing everything as fast as possible. Of course, it should do as well as it can, but for example, if it needs to sacrifice a bit of speed to get full performance in mutability today, that might be a good trade-off if the goal is to be a glue AD. Perfect? No, but if that’s that would give you the coverage to then tell the user that if they need more on a particular piece of code, seek out more. To seek out more performance, users could just have Zygote call ReverseDiff.jl on a function and have that compile the tape (or other specialized AD systems which will be announced more broadly soon), or may want to write a partial adjoint.

So (4) is really the kicker. If I was to hit a slow mutating code today inside of a differential equation, it would probably be something perfect for ModelingToolkit.jl to handle, so the best thing to do is to build hyper-optimized adjoints of that differential equation using ModelingToolkit.jl. At that level, I can symbolically handle it and generate code that a compiler because I can make a lot of extra assumptions, like cos^2(x) + sin^2(x) = 1 in my mathematical context. I can move code around, auto-parallelize it, etc. easily because of the simple static graph I’m working on. Wouldn’t it be a treat to just `@ModelingToolkitAdjoint f` and bingo now it’s using ModelingToolkit on a portion of code? `@ForwardDiffAdjoint f` to tell it that it “you should forward mode mere”. Yota.jl is a great reverse mode project, so `@YotaAdjoint f` and boom that could be more optimized than Zygote on some cases. `@ReverseDiff f` and let it compile the tape and it’ll get fairly optimal on the places where ReverseDiff.jl is applicable.

Julia is the perfect language to develop such a system for because its AST is so nice and constrained for mathematical contexts that all of these AD libraries do not work on a special DSL language like TensorFlow graphs or torch.numpy, but instead work directly on the language itself and its original existing libraries. With ChainRules.jl allowing for adjoint overloads that apply to all AD packages, focusing on these “Glue AD” properties could really open up the playing field, allowing Zygote to be at the center of an expansive differentiable programming world that works everywhere, maybe makes some compromises to do so, but then gives a system for other developers to make assumptions and define easily define adjoints and plug alternative AD systems into the whole game. This is a true mixed mode which incorporates not just forward and reverse, but also different implementations with different performance profiles (and this can be implemented just through ChainRules overloads!). Zygote would then facilitate this playing field with just a solid debugging and profiling experience, along with a very high chance of working on your code on your first try. That plus buy-in by package authors would be a true solution to differentiable programming.

The post Glue AD for Full Language Differentiable Programming appeared first on Stochastic Lifestyle.

Solving Sicherman dice puzzle using DataFrames.jl

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2020/06/14/Sicherman-dice.html

Sicherman dice puzzle

In many games players roll two normal dice to get a result that is then later
used do decide the course of play.

By normal die we understand a 6-sided die with faces numbered from 1 do 6.

Now the puzzle is to check if there exists other pairs of two 6-sided dice
with faces numbered with positive integers and have the same probability
distribution for the sum as normal dice.

A standard approach to answer this question is to use generating functions to
show that it is actually possible, see e.g. Wikipedia article about
Sicherman dice.

However, in this post we will want to use DataFrames.jl to enumerate the
possible solutions and find the feasible ones. The exercise mainly showcases
the new API for the filter function.

The code was tested under Julia 1.4.2 and DataFrames.jl 0.21, so please make
sure you have a their proper versions (the examples should work under any
post 1.0 release of Julia, but require DataFrames.jl to be at least 0.21).

Getting the reference probability distribution

First let us get the distribution of outcomes on a pair of normal dice.
We start with defining the getdist function that takes the numbers on
faces of the dies and returns their distribution:

function getdist(d1, d2)
    min1, max1 = extrema(d1)
    min2, max2 = extrema(d2)
    @assert min1 > 0 && min2 > 0
    d = zeros(Rational{Int}, max1 + max2)
    for p1 in d1, p2 in d2
        s = p1 + p2
            d[s] += 1
    end
    d .//= length(d1) * length(d2)
    return d
end

In the function we assume that sides of both dies are numbered with positive
integers. As we are working with a finite probability space all probabilities
are rationals so we store them as Rational{Int} type to avoid rounding. We
could have used Float64 (or even just Int without doing normalization), but
as we will soon learn our code will be still fast enough so no such
approximation is required.

To test the code let us get a distribution for two normal dice and store it in
the NORMAL_DIST constant (we will use it later):

julia> const NORMAL_DIST = getdist(1:6, 1:6)
12-element Array{Rational{Int64},1}:
 0//1
 1//36
 1//18
 1//12
 1//9
 5//36
 1//6
 5//36
 1//9
 1//12
 1//18
 1//36

julia> sum(NORMAL_DIST)
1//1

In the second line we have checked that we actually have a probability
distribution, as all its entries add up to 1.

Generating all possible dice

In the next step we generate all possible 6-sided die that possibly could be
used to generate the NORMAL_DIST distribution. In order to use some new features
of DataFrames.jl let us make one observation. Note that the sum equal to 2
is obtained with 1/36 probability. This means that each dice must have 1 on
its face exactly once.

Now what is the maximal possible value on the face? We see that the sum of two
maximal values is 12 and is obtained with 1/36 probability. This means that
the maximal value must be unique on both dice. But this implies that it must be
at most 8:

  • if it were 11, then the other die would have to have only 1 on all sides
    which is not possible;
  • if it were 10, then the other die would have to have one 1 and five 2s
    on its sides, which again is not allowed;
  • if it were 9, then the other die would have to be (1,2,2,2,2,3), but this
    would mean that the probability of rolling 11 would be at least 1/9 and it
    must be equal to 1/18.

So let us create a data frame, call it df1, with six columns, where each
column represents a single side of the die, and rows represent possible numbers
on its sides:

julia> using DataFrames

julia> df1 = DataFrame(Iterators.product((2:8 for i in 1:5)...));

julia> insertcols!(df1, 1, "0" => 1);

julia> show(df1, eltypes=false)
16807×6 DataFrame
│ Row   │ 0 │ 1 │ 2 │ 3 │ 4 │ 5 │
├───────┼───┼───┼───┼───┼───┼───┤
│ 1     │ 1 │ 2 │ 2 │ 2 │ 2 │ 2 │
│ 2     │ 1 │ 3 │ 2 │ 2 │ 2 │ 2 │
│ 3     │ 1 │ 4 │ 2 │ 2 │ 2 │ 2 │
│ 4     │ 1 │ 5 │ 2 │ 2 │ 2 │ 2 │
│ 5     │ 1 │ 6 │ 2 │ 2 │ 2 │ 2 │
│ 6     │ 1 │ 7 │ 2 │ 2 │ 2 │ 2 │
│ 7     │ 1 │ 8 │ 2 │ 2 │ 2 │ 2 │
│ 8     │ 1 │ 2 │ 3 │ 2 │ 2 │ 2 │
⋮
│ 16799 │ 1 │ 7 │ 7 │ 8 │ 8 │ 8 │
│ 16800 │ 1 │ 8 │ 7 │ 8 │ 8 │ 8 │
│ 16801 │ 1 │ 2 │ 8 │ 8 │ 8 │ 8 │
│ 16802 │ 1 │ 3 │ 8 │ 8 │ 8 │ 8 │
│ 16803 │ 1 │ 4 │ 8 │ 8 │ 8 │ 8 │
│ 16804 │ 1 │ 5 │ 8 │ 8 │ 8 │ 8 │
│ 16805 │ 1 │ 6 │ 8 │ 8 │ 8 │ 8 │
│ 16806 │ 1 │ 7 │ 8 │ 8 │ 8 │ 8 │
│ 16807 │ 1 │ 8 │ 8 │ 8 │ 8 │ 8 │

After loading the DataFrames.jl package, we first create a df1 data frame
with sides not equal to 1 (remember that there is exactly one such side).
For this we use the Iterators.product function, which gives us an iterator
that is a product of passed iterators. As the side with value 1 on it is
excluded we used five such iterators. Also note that by default DataFrame
constructor treats the values produced by the iterator as rows of the produced
data frame.

Next we use insertcols! function to insert a column containing only 1 in the
first position and name it "0" (the reason will be soon seen – the
DataFrame constructor by default has named the other columns as "1", "2",
etc.). Note that the column_name => value syntax of insertcols! performs
automatic broadcasting of single values if needed (similar behavior is
implemented in DataFrame constructor, select, transform, and combine).

Note that we could have just written:

julia> DataFrame(Iterators.product(1, (2:8 for i in 1:5)...))
16807×6 DataFrame
│ Row   │ 1     │ 2     │ 3     │ 4     │ 5     │ 6     │
│       │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │
├───────┼───────┼───────┼───────┼───────┼───────┼───────┤
│ 1     │ 1     │ 2     │ 2     │ 2     │ 2     │ 2     │
│ 2     │ 1     │ 3     │ 2     │ 2     │ 2     │ 2     │
│ 3     │ 1     │ 4     │ 2     │ 2     │ 2     │ 2     │
│ 4     │ 1     │ 5     │ 2     │ 2     │ 2     │ 2     │
│ 5     │ 1     │ 6     │ 2     │ 2     │ 2     │ 2     │
│ 6     │ 1     │ 7     │ 2     │ 2     │ 2     │ 2     │
│ 7     │ 1     │ 8     │ 2     │ 2     │ 2     │ 2     │
│ 8     │ 1     │ 2     │ 3     │ 2     │ 2     │ 2     │
⋮
│ 16799 │ 1     │ 7     │ 7     │ 8     │ 8     │ 8     │
│ 16800 │ 1     │ 8     │ 7     │ 8     │ 8     │ 8     │
│ 16801 │ 1     │ 2     │ 8     │ 8     │ 8     │ 8     │
│ 16802 │ 1     │ 3     │ 8     │ 8     │ 8     │ 8     │
│ 16803 │ 1     │ 4     │ 8     │ 8     │ 8     │ 8     │
│ 16804 │ 1     │ 5     │ 8     │ 8     │ 8     │ 8     │
│ 16805 │ 1     │ 6     │ 8     │ 8     │ 8     │ 8     │
│ 16806 │ 1     │ 7     │ 8     │ 8     │ 8     │ 8     │
│ 16807 │ 1     │ 8     │ 8     │ 8     │ 8     │ 8     │

to get a similar result (but with different column names), but I wanted to show
the use of the insertcols! function.

Finally we show our data frame, but pass eltypes=false keyword argument to
avoid printing the column type information, as we do not need it.

We immediately notice that some rows in our df1 data frame are duplicates. For
example row 2 and row 8 represent the same die (permutation of numbers on sides
does not affect the distribution of outcomes). We get rid of the duplicates
in-place by requiring that the numbers on sides are sorted in the filter!
function:

julia> filter!(AsTable(:) => issorted, df1)
462×6 DataFrame
│ Row │ 0     │ 1     │ 2     │ 3     │ 4     │ 5     │
│     │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 2     │ 2     │ 2     │ 2     │
│ 2   │ 1     │ 2     │ 2     │ 2     │ 2     │ 3     │
│ 3   │ 1     │ 2     │ 2     │ 2     │ 3     │ 3     │
│ 4   │ 1     │ 2     │ 2     │ 3     │ 3     │ 3     │
│ 5   │ 1     │ 2     │ 3     │ 3     │ 3     │ 3     │
│ 6   │ 1     │ 3     │ 3     │ 3     │ 3     │ 3     │
│ 7   │ 1     │ 2     │ 2     │ 2     │ 2     │ 4     │
│ 8   │ 1     │ 2     │ 2     │ 2     │ 3     │ 4     │
⋮
│ 454 │ 1     │ 6     │ 7     │ 8     │ 8     │ 8     │
│ 455 │ 1     │ 7     │ 7     │ 8     │ 8     │ 8     │
│ 456 │ 1     │ 2     │ 8     │ 8     │ 8     │ 8     │
│ 457 │ 1     │ 3     │ 8     │ 8     │ 8     │ 8     │
│ 458 │ 1     │ 4     │ 8     │ 8     │ 8     │ 8     │
│ 459 │ 1     │ 5     │ 8     │ 8     │ 8     │ 8     │
│ 460 │ 1     │ 6     │ 8     │ 8     │ 8     │ 8     │
│ 461 │ 1     │ 7     │ 8     │ 8     │ 8     │ 8     │
│ 462 │ 1     │ 8     │ 8     │ 8     │ 8     │ 8     │

Notice that we have significantly reduced the number of possibilities this way
– from 16807 to 462.

In the filter! call note that we used AsTble(:) => issorted predicate
specifier. It means that each row of a DataFrame is converted to a
NamedTuple before being passed to the issorted function.

Going from one die to two dice

In df1 we have all possible configuration of one die. Now let us generate
all possibilities for two dice:

julia> df2 = crossjoin(df1, df1, makeunique=true);

julia> rename!(df2, [Symbol(die, side) for die in ["l", "r"] for side in 1:6]);

julia> show(df2, eltypes=false)
213444×12 DataFrame
│ Row    │ l1 │ l2 │ l3 │ l4 │ l5 │ l6 │ r1 │ r2 │ r3 │ r4 │ r5 │ r6 │
├────────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┤
│ 1      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │
│ 2      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 2  │ 2  │ 2  │ 2  │ 3  │
│ 3      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 2  │ 2  │ 2  │ 3  │ 3  │
│ 4      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 2  │ 2  │ 3  │ 3  │ 3  │
│ 5      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 2  │ 3  │ 3  │ 3  │ 3  │
│ 6      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 3  │ 3  │ 3  │ 3  │ 3  │
│ 7      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 2  │ 2  │ 2  │ 2  │ 4  │
│ 8      │ 1  │ 2  │ 2  │ 2  │ 2  │ 2  │ 1  │ 2  │ 2  │ 2  │ 3  │ 4  │
⋮
│ 213436 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 6  │ 7  │ 8  │ 8  │ 8  │
│ 213437 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 7  │ 7  │ 8  │ 8  │ 8  │
│ 213438 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 2  │ 8  │ 8  │ 8  │ 8  │
│ 213439 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 3  │ 8  │ 8  │ 8  │ 8  │
│ 213440 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 4  │ 8  │ 8  │ 8  │ 8  │
│ 213441 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 5  │ 8  │ 8  │ 8  │ 8  │
│ 213442 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 6  │ 8  │ 8  │ 8  │ 8  │
│ 213443 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 7  │ 8  │ 8  │ 8  │ 8  │
│ 213444 │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │ 1  │ 8  │ 8  │ 8  │ 8  │ 8  │

Using crossjoin we generate all possible combinations of both dice. We use
makeunique=true as we pass df1 as the left and right data frame in cross
join. Therefore we next rename! the data frame that we got to properly name
its columns so that they clearly show if we are considering left or right data
frame and which side of it (note that now we number sides from 1 to 6).

Final step – finding Sicherman dice

So now we want to find if our df2 data frame contains any pair of dice that
produces the same probability distribution as NORMAL_DIST we have computed
above. For this we define a helper function:

function test_dice(x...)
    d1 = ntuple(i -> x[i], 6)
    d2 = ntuple(i -> x[i + 6], 6)
    return d1 <= d2 && getdist(d1, d2) == NORMAL_DIST
end

The function assumes it is passed twelve positional arguments (soon these will
be values sored in a single row of our data frame). It constructs d1 and d2
tuples from them, to represent the dice. First we check if d1 <= d2 to avoid
permuted duplicates in the results, and if this test passes we check if the
probability distribution produced by our dice is the same as NORMAL_DIST.

Let us run the filter function to find the solution of the Sicherman dice
puzzle:

julia> show(filter(All() => test_dice, df2), eltypes=false)
2×12 DataFrame
│ Row │ l1 │ l2 │ l3 │ l4 │ l5 │ l6 │ r1 │ r2 │ r3 │ r4 │ r5 │ r6 │
├─────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┤
│ 1   │ 1  │ 2  │ 2  │ 3  │ 3  │ 4  │ 1  │ 3  │ 4  │ 5  │ 6  │ 8  │
│ 2   │ 1  │ 2  │ 3  │ 4  │ 5  │ 6  │ 1  │ 2  │ 3  │ 4  │ 5  │ 6  │

Indeed we get that there exists only one pair of dice different from two normal
dice that meets our requirements, namely (1,2,2,3,3,4) and (1,3,4,5,6,8).
A surprising finding indeed!

Note that this time in filter we have used All() => test_dice predicate.
It means that test_dice for each row of a data frame is passed all its column
as positional arguments.

Concluding remarks

I hope you found the examples interesting and giving you some insight how
filtering in DataFrames.jl can be used efficiently.

Note that the proposed codes are not only relatively terse but quite fast:

julia> @time begin
           df1 = DataFrame(Iterators.product((2:8 for i in 1:5)...))
           insertcols!(df1, 1, "0" => 1)
           filter!(AsTable(:) => issorted, df1)
           df2 = crossjoin(df1, df1, makeunique=true)
           rename!(df2, [Symbol(die, side) for die in ["l", "r"] for side in 1:6])
           @time filter(All() => test_dice, df2)
       end
  0.306816 seconds (307.21 k allocations: 62.547 MiB, 2.62% gc time)
2×12 DataFrame
│ Row │ l1    │ l2    │ l3    │ l4    │ l5    │ l6    │ r1    │ r2    │ r3    │ r4    │ r5    │ r6    │
│     │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 2     │ 3     │ 3     │ 4     │ 1     │ 3     │ 4     │ 5     │ 6     │ 8     │
│ 2   │ 1     │ 2     │ 3     │ 4     │ 5     │ 6     │ 1     │ 2     │ 3     │ 4     │ 5     │ 6     │

As you can see we are able to solve the puzzle in sub-second time.

I have not attempted to reproduce these examples using data frames in R or
Python, as I do not feel competent enough to write exemplary codes for these
environments. However, I would be interested to see how to reproduce the steps
I have shown here and how fast they would run.

If someone would be interested to make such an implementation and its benchmark
please contact me with your proposal and I will update this post below,
giving a solution and a credit to the submitter.