Category Archives: Julia

Building a VCR Clone in 28 Lines of Julia

By: julia on Chris de Graaf

Re-posted from: https://cdg.dev/post/vcr/

The Ruby gem VCR is a tool that allows you to “[r]ecord your test suite’s HTTP interactions and replay them during future test runs for fast, deterministic, accurate tests”. It’s really useful for testing things like web API wrappers. It also contains 3,000 lines of code. Let’s use Cassette to build it in Julia with less than 1% of the code!
The first example in VCR’s README looks like this:

Walk on Spheres Method in Julia

By: Philip Zucker

Re-posted from: https://www.philipzucker.com/walk-on-spheres-method-in-julia/

I saw a cool tweet (and corresponding conference paper) by Keenan Crane

http://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/index.html

I was vaguely aware that one can use a Monte Carlo method to solve the boundary value Laplace equation \nabla^2 \phi = 0 , but I don’t think I had seen the walk on spheres variant of it before. I think Crane’s point is how similar all this is to stuff graphics people already do and do well. It’s a super cool paper. Check it out.

Conceptually, I think it is plausible that the Laplace equation and a monte carlo walk are related because the static diffusion equation \nabla^2 n = 0 from Fick’s law ultimately comes from the brownian motion of little guys wobbling about from a microscopic perspective.

Slightly more abstractly, both linear differential equations and random walks can be describe by matrices, a finite difference matrix (for concreteness) K and a transition matrix of jump probabilities T. The differential equation is discretized to Kx=b and the stationary probability distribution is Tp=b, where b are sources and sinks at the boundary.

The mean value property of the Laplace equation allows one to speed this process up. Instead of having a ton of little walks, you can just walk out randomly sampling on the surface of big spheres. en.wikipedia.org/wiki/Walk-on-spheres_method. Alternatively you can think of it as eventually every random walk exits a sphere, and it is at a random spot on it.

So here’s the procedure. Pick a point you want the value of \phi at. Make the biggest sphere you can that stays in the domain. Pick a random point on the sphere. If that point is on the boundary, record that boundary value, otherwise iterate. Do this many many times, then the average value of the boundaries you recorded it the value of \phi

This seems like a good example for Julia use. It would be somewhat difficult to code this up efficiently in python using vectorized numpy primitives. Maybe in the future we could try parallelize or do this on the GPU? Monte carlo methods like these are quite parallelizable.

The solution of the 1-d Laplace equation is absolutely trivial. If the second derivative is 0, then $\phi = a + b x $. This line is found by fitting it to the two endpoint values.

So we’re gonna get a line out

using LinearAlgebra
avg = 0
phi0 = 0
phi1 = 10
x_0 = 0.75
function monte_run(x)
    while true
            l = rand(Bool) # go left?
            if (l && x <= 0.5) # finish at left edge 0
                return phi0
            elseif (!l && x >= 0.5) # finish at right edge 1
                return phi1
            else
                if x <= 0.5 # move away from 0
                    x += x
                else
                    x -= 1 - x # move away from 1
                end
            end
    end
end

monte_runs = [monte_run(x) for run_num =1:100, x=0:0.05:1 ]
import Statistics
avgs = vec(Statistics.mean( monte_runs , dims=1))
stddevs = vec(Statistics.std(monte_runs, dims=1)) ./ sqrt(size(monte_runs)[1]) # something like this right?

plot(0:0.05:1, avgs, yerror=stddevs)
plot!(0:0.05:1,  (0:0.05:1) * 10 )

And indeed we do.

You can do a very similar thing in 2d. Here I use the boundary values on a disc corresponding to x^2 – y^2 (which is a simple exact solution of the Laplace equation).



function monte_run_2d(phi_b, x)
    while true
            r = norm(x)
            if r > 0.95 # good enough
                return phi_b(x)
            else
                dr = 1.0 - r #assuming big radius of 1
                θ = 2 * pi * rand(Float64) #
                x[1] += dr * cos(θ)
                x[2] += dr * sin(θ)
            end
    end
end


monte_run_2d( x -> x[1],  [0.0 0.0] )


monte_runs = [monte_run_2d(x -> x[1]^2 - x[2]^2 ,  [x 0.0] ) for run_num =1:1000, x=0:0.05:1 ]

import Statistics
avgs = vec(Statistics.mean( monte_runs , dims=1))
stddevs = vec(Statistics.std(monte_runs, dims=1)) ./ sqrt(size(monte_runs)[1]) # something like this right?
plot(0:0.05:1, avgs, yerror=stddevs)
plot!(0:0.05:1,  (0:0.05:1) .^2 )

There’s more notes and derivations in my notebook here https://github.com/philzook58/thoughtbooks/blob/master/monte_carlo_integrate.ipynb

Random-based indexing for arrays

By: Mosè Giordano

Re-posted from: http://giordano.github.io/blog/2020-05-23-random-array-indexing/

image

Image credit: “xkcd: Donald Knuth” (CC-BY-NC
2.5)

Let me introduce a package I recently wrote:
RandomBasedArrays.jl, a
hassle-free package in the Julia programming language
for dealing with arrays. Every time you access an element of an array, the
first index is random, so this package relieves you from having to remember
whether Julia uses 0- or 1-based indexing: you simply cannot ever know what the
initial element will be. As an additional benefit, you can use any Int to
index a RandomBasedArray.

Motivation

This package takes a new stance in the longstanding debate whether arrays should
have 0-based or 1-based
indexing.

The main motivation for this package is that I’m sick of reading about this
debate. It is incredibly hard to convince people that there is no “one size
fits all” indexing in programming, both alternatives have their merits:

  • 0-based indexing is natural every time you deal with offsets, e.g. when
    referencing memory addresses, or when doing modular arithmetic
  • 1-based indexing is natural when you are counting elements: the 1st element is
    “1”, the 2nd element is “2”, etc…

It is pointless to claim the superiority of one indexing over the other one, as
they’re useful in different situations.

As a matter of fact, many “math-oriented” languages (e.g., Fortran, Julia,
Mathematica, MATLAB, R), that are less likely to fiddle with pointers’
addresses, default to 1-based indexing, even though probably the majority of the
programming languages nowadays uses 0-based indexing.

Many people claim superiority of 0-based indexing over 1-based indexing because
of a popular note by Edsger
W. Dijkstra: “Why numbering
should start at
zero”.
However, I personally find this note unconvincing and partial, mainly based on
subjective arguments (like elegance and ugliness) rather than really compelling
reasons. That said, 0-based indexing is certainly useful in many situations.

A good programming language, whatever indexing convention it uses, should
provide an abstraction layer to let users forget which is the initial index.
For example, Fortran has lbound
to reference the first element of an array. Besides the
first
function to reference the first element of a collection, the Julia programming
language has different utilities to iterate over collections:

  • arrays are iterables, this means that you can write a for loop like

    for element in my_array
    # do things with the `element`...
end

    without using indices at all

  • the length
    function gives you the length of a collection, so that Dijkstra’s argument
    about the length of a sequence remains aesthetic rather than practical
  • eachindex is
    a function that returns an iterable object with all the indices of the array.

Some times you need to use the indices of an array and know which is the first
one. In this case, the abstraction layer above is not useful. Thus, I think it
is important for a programming language to provide also a way to easily switch
to the most appropriate indexing for the task at hand. In Fortran you can set a
different initial index for an array with the
dimension
statement. Julia allows you to define custom indices for you new array-like
type, as described in the
documentation. The
most notable application of custom indices is probably the
OffsetArrays.jl package.
Other use cases of custom indices are shown in this blog
post.

Usage

Let’s come to RandomBasedArrays.jl. You can install it with Julia built-in
package manager. In a Julia
session, after entering the package manager mode with ], run the command

pkg> add RandomBasedArrays

After that, just run using RandomBasedArrays in the REPL to load the package.
It defines a new AbstractArray type, RandomBasedArray, which is a thin
wrapper around any other AbstractArray.

julia> using RandomBasedArrays

julia> A = RandomBasedArray(reshape(collect(1:12), 3, 4))
3×4 Array{Int64,2}:
 10  2  6   6
  8  9  1  11
  7  8  8   7

Every time you access an element, you conveniently get a random element of the
parent array, making any doubt about first index go away. This means that you
can use any Int as index, including negative numbers:

julia> A[-35]
6

julia> A[-35]
9

julia> A[-35]
4

You can also set elements of the array in the usual way, just remember that
you’ll set a random element of the parent array:

julia> A[28,-19] = 42
42

julia> A
5×5 Array{Int64,2}:
 13  16   3  25   9
 23  20  16  18   1
  5  17  21   6   8
  5   3  42  10  13
 25   6  23   4  11

julia> A
5×5 Array{Int64,2}:
  9  25  25   3   3
  4  14   9   7  18
 22  14  13  21   2
 11  12  19  14  19
 19   2  21   2  21

There are other Julia packages that play with different indexing of arrays, e.g.:

  • OffsetArrays.jl:
    Fortran-like arrays with arbitrary, zero or negative starting indices
  • TwoBasedIndexing.jl:
    Two-based indexing (note: this is currently out-of-date and doesn’t work in
    Julia 1.0)