From: andrew cooke <andrew@...>

Date: Tue, 27 Aug 2013 20:52:40 -0400

I wanted to have a closer look at Julia so I decided to revisit a program I worked on years ago - generating an image of crumpled paper using a "functional images" approach (the original was in Haskell with Pancito, based on Conal Elliot's Pan). I'll include the code below, but here I want to focus more on the "higher level" details of what it was like using the language. Getting Started There don't seem to be pre-built binaries for OpenSuse, but building from source works fine (although it takes a long time - an hour or so I guess). Simply clone https://github.com/JuliaLang/julia and "make". Positives and Negatives In general, it was a very pleasant, productive experience - quite similar to Python (at least, the functional part). I didn't use anything fancy (no multimethods), but I did use generators and I did try "automatic" parallelism (via pmap). Note that it is still very much in development (apparently about to release 0.2) so you're building and working with a moving target. One surprising detail was that indices / ranges are 1-based and inclusive - I guess for consistency with Fortran and Matlab. Another, similarly small but annoying wrinkle is the use of "elseif" (sic) to solve the nested if/else parsing problem (why can't we standardise on "elif"?). Back to a positive: there are already many packages, and the packaging system worked fine. I installed the Cairo (graphics) wrapper with little fuss. The package itself had no documentation, but it's just a simple wrapper around the standard C API, so it was easy enough to translate existing C examples. Unfortunately it's not a complete wrapper - details like line end caps aren't exposed. I had a problem understanding what is automatically broadcast when code runs in parallel. It turns out that it's similar to MP or Python's multiprocessing, in that you have to do some work to ensure that the environment in each process is initialised correctly. In this case by using "require()" to load the same source file on each. The way to use pmap correctly is in the documentation (which is pretty good), but I somehow missed it. Instead, I asked on the google group / mailing list and got support there. The list is not terribly busy (took a while to get an answer, and I have another question waiting), but is friendly enough and seems to be read by the developers. I also tried IRC, but no-one responded. I can't really comment on speed as I don't have alternative implementations to compare, and the problem is one where it's easy to increase parameters until you get bored waiting for results (my current best takes 500m to run). However, adding type annotations didn't change anything - looking at the docs more closely it seems like the speed comes mainly from JIT generation of specialised code, with annotations only helping you find places where types jump around (and so cause prediction issues). Conclusions In summary: positive. Way, way better than using C or Fortran. So, if it has similar speed (and I assume it does), then it's freaking awesome. BUT it's still very much in development, and with a small community for support. I imagine it will get better with time, and I hope the developers stick with it until they get the success they deserve. It has replaced Clojure as my best candidate for a "better Python". Code And here's the code, in two files. paper.jl: # 0 is black, 1 is white # the paper extends from (0,0) to (1,1) # the underlying image is a function from position and incident light to # brightness. the need for the light is something of a hack - it allows # us to calculate the effect of folds on brightness even when we're # still unsure whether we're evaluating the brightness of the paper (or # are outside, in the background). const border_shadow = 0.075 const border = 0.05 # less than shadow as entire image shrinks from folds const plain_shade = 0.65 function to_bottom_corner(x) if 0 < x < 1 return 0 elseif x > 1 return 1 - x else return x end end function undistorted(xy, light) x, y = viewport(xy) if 0 < x < 1 && 0 < y < 1 return min(1, max(0, light)) else x = to_bottom_corner(x) y = to_bottom_corner(y) r = sqrt(x^2 + y^2) return min(1, r / border_shadow) end end # this transform shifts the paper (originally also in the unit square) # and border (originally outside) into the unit square (it just makes # things neater if we can use rand() at the top level). function viewpoint(x) return x * (1 + 2 * border) - border end function viewport(xy) x, y = xy return (viewpoint(x), viewpoint(y)) end # now the meat. given a fold (position and angle), transform the x, y # coordinates and the incident light. these are just nice-looking # approximations - there's no real physics here. function distance(xy, fold) x, y = xy (p, q), theta = fold return (x - p) * sin(theta) - (y - q) * cos(theta) end const n_folds = 200 const fold_width = 0.04 const fold_shift = 0.15 / n_folds const fold_grey = 0.15 function normalize(d) d = d / fold_width if d > 1 return 1 elseif d < -1 return -1 else return sign(d) * abs(d) ^ 1.5 end end function shift(xy, theta, d) x, y = xy return (x + d * fold_shift * sin(theta), y - d * fold_shift * cos(theta)) end function shade(g, theta, d) return g * (1 + (1 - abs(d)) * sign(d) * cos(theta + pi/3) * fold_grey) end # a stream of random folds. function random_folds() while true produce(((rand(), rand()), pi * rand())) end end # take 'n' folds from the stream, and an underlying function (the paper) # and create a new function that accumulates the effect of all the folds. # an alternative would be to make an explicit iteration over the folds # (i have no idea which is faster, but this seems neater). function combine_folds(n, folds, base) if n == 0 return base else fold = consume(folds) pq, theta = fold function wrapped(xy, light) d = normalize(distance(xy, fold)) xy = shift(xy, theta, d) return base(xy, shade(light, theta, d)) end return combine_folds(n-1, folds, wrapped) end end # to render the image we evaluate it and random points. function random_points() while true produce((rand(), rand())) end end function take(n, seq) function inner() for i = 1:n produce(consume(seq)) end end return Task(inner) end function driver(n, image, output) points = take(n, Task(random_points)) for result in map(point -> (point, image(point, plain_shade)), points) output(result...) end end typealias Point (Float64, Float64) function pdriver(n, image, output) points = take(n, Task(random_points)) for result::(Point, Float64) in pmap(point::Point -> (point, image(point, plain_shade)), points) output(result...) end end # testing function print_output(xy, g) x, y = xy println("$x $y $g") end #driver(20, undistorted, print_output) # output to cairo pdf using Cairo function cairo_pdf(file, size) s = CairoPDFSurface(file, size, size) c = CairoContext(s) set_line_width(c, 1) set_source_rgb(c, 0, 0, 0) function output(xy, g) if g < rand() # set_source_rgb(c, g, g, g) x, y = xy move_to(c, x * size, y * size) rel_line_to(c, 1/sqrt(2), 1/sqrt(2)) end end for x = 0:1 for y = 0:1 output((x,y), 0) end end function close() stroke(c) finish(s) end return (output, close) end paper-run.jl: require("/home/andrew/project/paper/git/paper.jl") # here we go... const width = 4000 const n_points = 10000000 output, close = cairo_pdf("paper.pdf", width) paper = combine_folds(n_folds, Task(random_folds), undistorted) pdriver(n_points, paper, output) close() Above I am using the parallel code. Change pdriver to driver to run the serial version. The "const" aren't necessary - I just wondered if they would help seed things up... Andrew