| Andrew Cooke | Contents | Latest | RSS | Twitter | Previous | Next


Welcome to my blog, which was once a mailing list of the same name and is still generated by mail. Please reply via the "comment" links.

Always interested in offers/projects/new ideas. Eclectic experience in fields like: numerical computing; Python web; Java enterprise; functional languages; GPGPU; SQL databases; etc. Based in Santiago, Chile; telecommute worldwide. CV; email.

Personal Projects

Lepl parser for Python.

Colorless Green.

Photography around Santiago.

SVG experiment.

Professional Portfolio

Calibration of seismometers.

Data access via web services.

Cache rewrite.

Extending OpenSSH.

Last 100 entries

British Words; Chinese Govt Intercepts External Web To DDOS github; Numbering Permutations; Teenage Engineering - Low Price Synths; GCHQ Can Do Whatever It Wants; Dublinesque; A Cryptographic SAT Solver; Security Challenges; Word Lists for Crosswords; 3D Printing and Speaker Design; Searchable Snowden Archive; XCode Backdoored; Derived Apps Have Malware (CIA); Rowhammer - Hacking Software Via Hardware (DRAM) Bugs; Immutable SQL Database (Kinda); Tor GPS Tracker; That PyCon Dongle Mess...; ASCII Fluid Dynamics; Brandalism; Table of Shifter, Cassette and Derailleur Compatability; Lenovo Demonstrates How Bad HTTPS Is; Telegraph Owned by HSBC; Smaptop - Sunrise (Music); Equation Group (NSA); UK Torture in NI; And - A Natural Extension To Regexps; This Is The Future Of Religion; The Shazam (Music Matching) Algorithm; Tributes To Lesbian Community From AIDS Survivors; Nice Rust Summary; List of Good Fiction Books; Constructing JSON From Postgres (Part 2); Constructing JSON From Postgres (Part 1); Postgres in Docker; Why Poor Places Are More Diverse; Smart Writing on Graceland; Satire in France; Free Speech in France; MTB Cornering - Where Should We Point Our Thrusters?; Secure Secure Shell; Java Generics over Primitives; 2014 (Charlie Brooker); How I am 7; Neural Nets Applied to Go; Programming, Business, Social Contracts; Distributed Systems for Fun and Profit; XML and Scheme; Internet Radio Stations (Curated List); Solid Data About Placebos; Half of Americans Think Climate Change Is a Sign of the Apocalypse; Saturday Surf Sessions With Juvenile Delinquents; Ssh, tty, stdout and stderr; Feathers falling in a vacuum; Santiago 30m Bike Route; Mapa de Ciclovias en Santiago; How Unreliable is UDP?; SE Santiago 20m Bike Route; Cameron's Rap; Configuring libxml with Eclipse; Reducing Combinatorial Complexity With Occam - AI; Sentidos Comunes (Chilean Online Magazine); Hilary Mantel: The Assassination of Margaret Thatcher - August 6th 1983; NSA Interceptng Gmail During Delivery; General IIR Filters; What's happening with Scala?; Interesting (But Largely Illegible) Typeface; Retiring Essentialism; Poorest in UK, Poorest in N Europe; I Want To Be A Redneck!; Reverse Racism; The Lost Art Of Nomography; IBM Data Center (Photo); Interesting Account Of Gamma Hack; The Most Interesting Audiophile In The World; How did the first world war actually end?; Ky - Restaurant Santiago; The Black Dork Lives!; The UN Requires Unaninmous Decisions; LPIR - Steganography in Practice; How I Am 6; Clear Explanation of Verizon / Level 3 / Netflix; Teenage Girls; Formalising NSA Attacks; Switching Brakes (Tektro Hydraulic); Naim NAP 100 (Power Amp); AKG 550 First Impressions; Facebook manipulates emotions (no really); Map Reduce "No Longer Used" At Google; Removing RAID metadata; New Bike (Good Bike Shop, Santiago Chile); Removing APE Tags in Linux; Compiling Python 3.0 With GCC 4.8; Maven is Amazing; Generating Docs from a GitHub Wiki; Modular Shelves; Bash Best Practices; Good Emergency Gasfiter (Santiago, Chile); Readings in Recent Architecture; Roger Casement; Integrated Information Theory (Or Not); Possibly undefined macro AC_ENABLE_SHARED; Update on Charges

© 2006-2013 Andrew Cooke (site) / post authors (content).

First impressions using Julia

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

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).


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".


And here's the code, in two files.


# 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
        return x

function undistorted(xy, light)
    x, y = viewport(xy)
    if 0 < x < 1 && 0 < y < 1
        return min(1, max(0, light))
        x = to_bottom_corner(x)
        y = to_bottom_corner(y)
        r = sqrt(x^2 + y^2)
        return min(1, r / border_shadow)

# 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

function viewport(xy)
    x, y = xy
    return (viewpoint(x), viewpoint(y))

# 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)

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
        return sign(d) * abs(d) ^ 1.5

function shift(xy, theta, d)
    x, y = xy
    return (x + d * fold_shift * sin(theta), y - d * fold_shift * cos(theta))

function shade(g, theta, d)
    return g * (1 + (1 - abs(d)) * sign(d) * cos(theta + pi/3) * fold_grey)

# a stream of random folds.

function random_folds()
    while true
        produce(((rand(), rand()), pi * rand()))

# 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
        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))
        return combine_folds(n-1, folds, wrapped)

# to render the image we evaluate it and random points.

function random_points()
    while true
        produce((rand(), rand()))

function take(n, seq)
    function inner()
        for i = 1:n
    return Task(inner)

function driver(n, image, output)
    points = take(n, Task(random_points))
    for result in map(point -> (point, image(point, plain_shade)), points)

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)

# testing

function print_output(xy, g)
    x, y = xy
    println("$x $y $g")

#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))
    for x = 0:1
        for y = 0:1
            output((x,y), 0)
    function close()
    return (output, close)



# 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)

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...


Increasing Julia Parallel Throughput

From: andrew cooke <andrew@...>

Date: Thu, 29 Aug 2013 16:01:36 -0400

In the example above I included pmap, which runs each calculation of the map
on a remote mode.  However, in practice, it was slower than running everything
on a single core.

The problem, it seems, was that the amount of work done, for the network
traffic required, was too small.  The coordinating core was at 100% while the
saves stayed at 10-20%.

The fix was to make each call do more work.  I did this by modifying the code
to (1) broadcast a random seed (rather than a point) and (2) calculating 10000
points, instead of 1, returning the results in an array.

Note - turns out that Julia has typed arrays (as well as untyped ones) so you
get an efficient, minimal block of memory.

This was a huge win.  Workers saturated at 100% (in fact it made sense to
specify one more process than I have cores, since the master process was doing
very little).

The new code is basically:

  typealias Point (Float64, Float64)
  typealias Result (Point, Float64)

  const block_size = 10000

  function block_eval(seed, image)
      a = Array(Result, block_size)
      for i = 1:block_size
	  point = random_point()
	  a[i] = (point, image(point, plain_shade))
      return a

  function driver(n, image, output)
      blocks = div(n, block_size)
      for block in map(seed -> block_eval(seed, image), 1:blocks)
	  for result in block

  function pdriver(n, image, output)
      blocks = div(n, block_size)
      for block in pmap(seed -> block_eval(seed, image), 1:blocks)
	  for result in block

Which is still simple and understandable (imho).


Calling C from Julia

From: andrew cooke <andrew@...>

Date: Thu, 29 Aug 2013 23:32:41 -0400

Missing the ability to set Cairo line caps, I read the manual and tried to do
it myself.  Turns out that calling C is easy:

  function set_line_cap(ctx::CairoContext, style)
      ccall((:cairo_set_line_cap, "libcairo"), 
            Void, (Ptr{Void},Int), ctx.ptr, style)

  set_line_cap(ctx, 1)

works just fine.  Nice.


Comment on this post