Tag Archives: Programming

What’s Wrong with Statistics in Julia?

By: John Myles White

Re-posted from: http://www.johnmyleswhite.com/notebook/2014/11/29/whats-wrong-with-statistics-in-julia/

Introduction

Several months ago, I promised to write an updated version of my old post, “The State of Statistics in Julia”, that would describe how Julia’s support for statistical computing has evolved since December 2012.

I’ve kept putting off writing that post for several reasons, but the most important reason is that all of my attention for the last few months has been focused on what’s wrong with how Julia handles statistical computing. As such, the post I’ve decided to write isn’t a review of what’s already been done in Julia, but a summary of what’s being done right now to improve Julia’s support for statistical computing.

In particular, this post focuses on several big changes to the core data structures that are used in Julia to represent statistical data. These changes should all ship when Julia 0.4 is released.

What’s Wrong with Statistics in Julia Today?

The primary problem with statistical computing in Julia is that the current tools were all designed to emulate R. Unfortunately, R’s approach to statistical computing isn’t amenable to the kinds of static analysis techniques that Julia uses to produce efficient machine code.

In particular, the following differences between R and Julia have repeatedly created problems for developers:

  • In Julia, computations involving scalars are at least as important as computations involving vectors. In particular, iterative computations are first-class citizens in Julia. This implies that statistical libraries must allow developers to write efficient code that iterates over the elements of a vector in pure Julia. Because Julia’s compiler can only produce efficient machine code for computations that are type-stable, the representations of missing values, categorical values and ordinal values in Julia programs must all be type-stable. Whether a value is missing or not, its type must remain the same.
  • In Julia, almost all end-users will end up creating their own types. As such, any tools for statistical computing must be generic enough that they can be extended to arbitrary types with little to no effort. In contrast to R, which can heavily optimize its algorithms for a very small number of primitive types, Julia developers must ensure that their libraries are both highly performant and highly abstract.
  • Julia, like most mainstream languages, eagerly evaluates the arguments passed to functions. This implies that idioms from R which depend upon non-standard evaluation are not appropriate for Julia, although it is possible to emulate some forms of non-standard evaluation using macros. In addition, Julia doesn’t allow programmers to reify scope. This implies that idioms from R that require access to the caller’s scope are not appropriate for Julia.

The most important way in which these issues came up in the first generation of statistical libraries was in the representation of a single scalar missing value. In Julia 0.3, this concept is represented by the value NA, but that representation will be replaced when 0.4 is released. Most of this post will focus on the problems created by NA.

In addition to problems involving NA, there were also problems with how expressions were being passed to some functions. These problems have been resolved by removing the function signatures for statistical functions that involved passing expressions as arguments to those functions. A prototype package called DataFramesMeta, which uses macros to emulate some kinds of non-standard evaluation, is being developed by Tom Short.

Representing Missing Values

In Julia 0.3, missing values are represented by a singleton object, NA, of type NAtype. Thus, a variable x, which might be either a Float64 value or a missing value encoded as NA, will end up with type Union(Float64, NAtype). This Union type is a source of performance problems because it defeats Julia’s compiler’s attempts to assign a unique concrete type to every variable.

We could remove this type-instability by ensuring that every type has a specific value, such as NaN, that signals missingness. This is the approach that both R and pandas take. It offers acceptable performance, but does so at the expense of generic handling of non-primitive types. Given Julia’s rampant usage of custom types, the sentinel values approach is not viable.

As such, we’re going to represent missing values in Julia 0.4 by borrowing some ideas from functional languages. In particular, we’ll be replacing the singleton object NA with a new parametric type Nullable{T}. Unlike NA, a Nullable object isn’t a direct scalar value. Rather, a Nullable object is a specialized container type that either contains one value or zero values. An empty Nullable container is taken to represent a missing value.

The Nullable approach to representing a missing scalar value offers two distinct improvements:

  • Nullable{T} provides radically better performance than Union(T, NA). In some benchmarks, I find that iterative constructs can be as much as 100x faster when using Nullable{Float64} instead of Union(Float64, NA). Alternatively, I’ve found that Nullable{Float64} is about 60% slower than using NaN to represent missing values, but involves a generic approach that trivially extends to arbitrary new types, including integers, dates, complex numbers, quaternions, etc…
  • Nullable{T} provides more type safety by requiring that all attempts to interact with potentially missing values explicitly indicate how missing values should be treated.

In a future blog post, I’ll describe how Nullable works in greater detail.

Categorical Values

In addition to revising the representation of missing values, I’ve also been working on revising our representation of categorical values. Working with categorical data in Julia has always been a little strange, because the main tool for representing categorical data, the PooledDataArray, has always occupied an awkward intermediate position between two incompatible objectives:

  • A container that keeps track of the unique values present in the container and uses this information to efficiently represent values as pointers to a pool of unique values.
  • A container that contains values of a categorical variable drawn from a well-defined universe of possible values. The universe can include values that are not currently present in the container.

These two goals come into severe tension when considering subsets of a PooledDataArray. The uniqueness constraint suggests that the pool should shrink, whereas the categorical variable definition suggests that the pool should be maintained without change. In Julia 0.4, we’re going to commit completely to the latter behavior and leave the problem of efficiently representing highly compressible data for another data structure.

We’ll also begin representing scalar values of categorical variables using custom types. The new CategoricalVariable and OrdinalVariable types that will ship with Julia 0.4 will further the efforts to put scalar computations on an equal footing with vector computations. This will be particularly notable for dealing with ordinal variables, which are not supported at all in Julia 0.3.

Metaprogramming

Many R functions employ non-standard evaluation as a mechanism for augmenting the current scope with the column names of a data.frame. In Julia, it’s often possible to emulate this behavior using macros. The in-progress DataFramesMeta package explores this alternative to non-standard evaluation. We will also be exploring other alternatives to non-standard evaluation in the future.

What’s Next

In the long-term future, I’m hoping to improve several other parts of Julia’s core statistical infrastructure. In particular, I’d like to replace DataFrames with a new type that no longer occupies a strange intermediate position between matrices and relational tables. I’ll write another post about those issues later.

A bioinformatics walk-through: Accessing protein-protein interaction interfaces for all known protein structures with PDBe PISA

If this summer’s posting became a little infrequent, part of the blame lies with computational research I’ve been working on, regarding the systems biology of chromosomal translocations and the ensuing chimeric proteins at the Medical Research Council Laboratory of Molecular Biology in Cambridge.

A sizeable part of bioinformatics ‘dry lab’ work falls into what has been described in the NYT as ‘data wrangling’ (or the work of a ‘data janitor’). This post is about accessing the data held in the Protein Databank in Europe’s repository of Proteins, Interfaces, Structures and Assemblies (PDBe PISA).

Sent out onto the web to find a source of structural protein-protein interaction data with amino acid-level resolution, my first port of call was the Nucleic Acids Research Molecular Biology Online Database Collection (which I’d read of in the opening chapters of Arthur Lesk’s Introduction to Bioinformatics) where I found a sizeable list of PPI databases.

Not wanting to click through each, I chose to browse this programmatically, using Javascript-automated AJAX requests (effectively asking the website to give me web pages but without displaying them) and just ‘scrape’ what I wanted (full workings here), as follows:

From these results, here’s a little background info on PDBe:

  {
    "name": "PDBe",
    "url": "http://www.ebi.ac.uk/pdbe/",
    "entryurl": "http://www.oxfordjournals.org/nar/database/summary/456",
    "desc": "EMBL-EBI's Protein Data Bank in Europe (PDBe) is the European resource for the collection, organization and dissemination of data about biological macromolecular structures. PDBe is one of four partners in the worldwide Protein Data Bank (wwPDB), the consortium entrusted with the collation, maintenance and distribution of the global repository of macromolecular structure data. PDBe uses a relational database that presents the data derived from the Protein Data Bank (PDB) in a consistent way and allows users to retrieve meaningful data using complex and sophisticated searches including simple textual queries or more complex 3D structure-based queries. PDBe has also developed a number of advanced tools for analysis of macromolecules. The \"Structure Integration with Function, Taxonomy and Sequence\" (SIFTS) initiative integrates data from a number of bioinformatics resources that is used by major global sequence, structure and protein-family resources. Furthermore, PDBe works actively with the X-ray crystallography, Nuclear Magnetic Resonance (NMR) spectroscopy and cryo-Electron Microscopy (EM) communities and is a partner in the Electron Microscopy Data Bank (EMDB). The active involvement with the scientific communities has resulted in improved tools for structure deposition and analysis.",
    "ref": null,
    "absurl": "http://nar.oxfordjournals.org/cgi/content/abstract/42/D1/D285",
    "email": "pdbe@ebi.ac.uk"
  },

Web scraping can feel quite kludgy, and there are doubtless better ways to do the above. Having said that, it’s great for prototyping: you can use Javascript within a web browser console, i.e. without littering your computer with temporary files. What’s more, dedicated communities like the ScraperWiki forum are around to support and develop the associated tools, and in its more elaborate incarnations ‘scraping’ features in journals like Briefings in Bioinformatics (“Web scraping technologies in an API World” was published there just this week).

After having decided on PDBe PISA thanks to my scraped-together report, and finding no guidance on how to tackle the task, I turned to the bioinformatician’s equivalent of [computing/programming Q&A site] Stack Overflow known as Biostars. My question got a grand total of 0 answers(!), so what follows is my approach — which may either be of interest as a peek into the work going under the banner of ‘bioinformatics’ or as a guide to other scientists seeking to access the same information.

First off, a Python script parcelled up a list of every PDB code (the unique identifier to an author-deposited structure from X-ray crystallography, NMR etc.) in PDB into comma-separated chunks of 50, which were stuck onto the end of a web-service query as recommended. The server would process these queries, understood through its “API”: the CGI of cgi-bin in the URL means it’s invoking a script on the server, which in turn expects interfaces.pisa? to be followed by comma-separated PDB codes. Given these expectations, the API will respond in a regular manner each time, enabling reliable scripting.

With over 2000 such queries for interface data (each of them requesting 50 PDB-code-identified structures), this isn’t something you want to be doing manually. It wasn’t clear exactly which PDB entries were needed at the time, so the full complement was downloaded.

This download script just works for one query, putting the received XML in one file – to handle all 2029 queries, a bit of lateral thinking was required. 50 queries (each containing 50 PDB codes) were executed to make up a single interfacesi–j.xml file, where i is an integer 1 to 4, and likewise j from 1 to 10 (plus a bonus 4-11 to get those final 29). Download scripts (named similarly as getxmli–j.py) were written individually by another script — code writing code…

With download scripts written, the task of running each of them consecutively fell to yet another Python script, playing the sound of Super Mario picking up a coin when each file finished downloading, or the Mario pause-game sound upon encountering an error, because I could because clear feedback becomes necessary on something taking days across multiple computers.

Inevitably a minority of the queries failed, and had to be obtained separately.

Once downloaded, various pattern matching text-processing programs were run on the XML from within a shell script — readers unfamiliar with programming may have heard of these this week thanks to the 22 year old security bug(s) being referred to as shellshock. Shell scripts make looping through files in this manner a simple task, and are becoming essential for everyday file manipulation now that I’m a reformed Windows user. For the 41 XML files, a function runprocessor was called, with instructions to:

  1. Split each file successively at every <pisa_interfaces> tag through to the closing </pisa_interfaces> tag, the line numbers of which were stored together in an ordered list (an “array variable”) pisapairs
  2. Write each of these sections to a cache file xmlcache.xml, of suitable size for parsing by a Python XML parser.
  3. Reduce the time spent by the parser by in turn splitting this cache into just the PDB entries in the shortlist of interest with a function extractsubsets
  4. Initiate a Python script to read the entire cachesubset.xml file into memory, and write the pertinent structural data into a report formatted as tab-separated values (TSV). This file is a mere few hundred megabytes compared to the 120 GB grand total for the XML.

Clicking Details for an interface on the list of all interfaces for a given protein structure, e.g. for the only one in spider silk precursor protein spidroin, shows the interfacial residues in yellow:
image

The output threads together all interfacial residues and the associated statistical figures for each on a single line for every interface, but it’s simple enough to separate out each according to commas (then colons) to get a longform residue-per-line output once all XML is processed.

Progress is indicated in terminal output, where the current i and j values are printed followed by the pisapair (i.e. which of the 50 pisa_interfaces tags) is being worked through:

image

As shown in the logfile, there are inevitable errors, such as Entry not found: it’s simple enough to find the difference between the output report file’s list of PDB codes and the input ‘shortlist’, which can be mapped back to the constituent files for any follow-up investigation (the “wrangling” facet of computational science alluded to earlier) since the order of the original 2029 queries is known:

I’m putting these together in a code repository on GitHub, with a disclaimer that it’s not fit for all purposes (for instance if you’re interested in H-bonds, in brown on the PISA website residue table, above).

A lot of this was painfully slow — there’s nothing to be done about the speed of downloading the files, given that its rate is limited by the server. Yes there was a lot of data to get through, but Python’s sluggishness at the final step makes me wonder if I could implement some leaner algorithm, parallelise, etc., but with term recommenced code optimisation on a successfully completed task isn’t top priority. Advice on improvements would be appreciated if you have any.

I’m currently reading Jones & Pevzner’s An Introduction to Bioinformatics Algorithms which gives insight into how you can analyse and improve these types of operations (the book is core reading for a Coursera.org lecture series which kicks off next month), and have been recommended Goldwasser & Tamassia’s Data Structures and Algorithms in Python (a few online resources in a similar vein are available here).

I’ve also been fiddling with Julia, an R-like language with C-like speeds — in a 2012 blog post its creators say they “created Julia, in short, because we are greedy”. Fernando Perez is overseeing its incorporation into IPython Notebooks as ‘Project Jupyter’ and a port of R’s ggplot2 library has recently emerged for Julia under the name of Gadfly (a tutorial IPy NB is up here).

I’m starting a final year undergraduate project as of this week, on mapping small RNA-seq data to miRNAs, under the supervision of the founder of the database central to cataloguing this class of non-coding RNA ‒ super exciting stuff! :¬)

If you’ve got questions on PISA you think I could help with, feel free to ask here, or shoot me an email.

✣ PDBe PISA homepage

✣ Peter Briggs, a scientific programmer at STFC Daresbury Laboratory, has a nice little guide to the service here.

Values vs. Bindings: The Map is Not the Territory

By: John Myles White

Re-posted from: http://www.johnmyleswhite.com/notebook/2014/09/06/values-vs-bindings-the-map-is-not-the-territory/

Many newcomers to Julia are confused by the seemingly dissimilar behaviors of the following two functions:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

julia> a = [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3
 
julia> function foo!(a)
           a[1] = 10
           return
       end
foo! (generic function with 1 method)
 
julia> foo!(a)
 
julia> a
3-element Array{Int64,1}:
 10
  2
  3
 
julia> function bar!(a)
           a = [1, 2]
           return
       end
bar! (generic function with 1 method)
 
julia> bar!(a)
 
julia> a
3-element Array{Int64,1}:
 10
  2
  3

Why does the first function successfuly alter the global variable a, but the second function does not?

To answer that question, we need to explain the distinction between values and bindings. We’ll start with a particularly simple example of a value and a binding.

In Julia, the number 1 is a value:

1
2

julia> 1
1

In contrast to operating on a value, the Julia assignment operation shown below creates a binding:

1
2

julia> a = 1
1

This newly created binding is an association between the symbolic name a and the value 1. In general, a binding operation always associates a specific value with a specific name. In Julia, the valid names that can be used to create bindings are symbols, because it is important that the names be parseable without ambiguity. For example, the string "a = 1" is not an acceptable name for a binding, because it would be ambiguous with the code that binds the value 1 to the name a.

This first example of values vs. bindings might lead one to believe that values and bindings are very easy to both recognize and distinguish. Unfortunately, the values of many common objects are not obvious to many newcomers.

What, for example, is the value of the following array?

1
2
3
4
5

julia> [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3

To answer this question, note that the value of this array is not defined by the contents of the array. You can confirm this by checking whether Julia considers two objects to be exactly identical using the === operator:

1
2
3
4
5

julia> 1 === 1
true
 
julia> [1, 2, 3] === [1, 2, 3]
false

The general rule is simple, but potentially non-intuitive: two arrays with identical contents are not the same array. To motivate this, think of arrays as if they were cardboard boxes. If I have two cardboard boxes, each of which contains a single ream of paper, I would not claim that the two boxes are the exact same box just because they have the same contents. Our intuitive notion of object identity is rich enough to distinguish between two containers with the same contents, but it takes some time for newcomers to programming languages to extend this notion to their understanding of arrays.

Because every container is distinct regardless of what it contains, every array is distinct because every array is its own independent container. An array’s identity is not defined by what it contains. As such, its value is not equivalent to its contents. Instead, an array’s value is a unique identifier that allows one to reliably distinguish each array from every other array. Think of arrays like numbered cardboard boxes. The value of an array is its identifier: thus the value of [1, 2, 3] is something like the identifier “Box 1″. Right now, “Box 1″ happens to contain the values 1, 2 and 3, but it will continue to be “Box 1″ even after its contents have changed.

Hopefully that clarifies what the value of an array is. Starting from that understanding, we need to re-examine bindings because bindings themselves behave like containers.

A binding can be thought of as a named box that can contain either 0 or 1 values. Thus, when a new Julia session is launched, the name a has no value associated with it: it is an empty container. But after executing the line, a = 1, the name has a value: the container now has one element in it. Being a container, the name is distinct from its contents. As such, the name can be rebound by a later operation: the line a = 2 will change the contents of the box called a to refer to the value 2.

The fact that bindings behave like containers becomes a source of confusion when the value of a binding is itself a container:

1

a = [1, 2, 3]

In this case, the value associated with the name a is the identifier of an array that happens to have the values 1, 2, and 3 in it. But if the contents of that array are changed, the name a will still refer to the same array — because the value associated with a is not the contents of the array, but the identifier of the array.

As such, there is a very large difference between the following two operations:

1
2

a[1] = 10
a = [1, 2]
  • In the first case, we are changing the contents of the array that a refers to.
  • In the second case, we are changing which array a refers to.

In this second case, we are actually creating a brand new container as an intermediate step to changing the binding of a. This new container has, as its initial contents, the values 1 and 2. After creating this new container, the name a is changed to refer to the value that is the identifier of this new container.

This is why the two functions at the start of this post behave so differently: one mutates the contents of an array, while the other mutates which array a name refers to. Because variable names in functions are local, changing bindings inside of a function does not change the bindings outside of that function. Thus, the function bar! does not behave as some would hope. To change the contents of an array wholesale, you must not change bindings: you must change the contents of the array. To do that, bar! should be written as:

1
2
3
4

function bar!(a)
    a[:] = [1, 2]
    return
end

The notation a[:] allows one to talk about the contents of an array, rather than its identifier. In general, you should not expect that you can change the contents of any container without employing some indexing syntax that allows you to talk about the contents of the container, rather than the container itself.