Uniformly shuffle 5 gigabytes of numpy data
Question: I'm training a neural network with about five gigabytes of data stored as
`numpy` arrays. The data are split into chunks of 100000 rows, and I've done
six cycles of training over all the chunks in a random order. Unfortunately,
the network has begun to overfit. I think it still has capacity to fit the
data more closely; my suspicion is that internal regularities within each
chunk are starting to contradict one another, and I need to shuffle the data
more thoroughly so that it can train on different combinations. I want to try
this before going to the trouble of getting more training data.

Does anyone know a good way to generate a new permutation of 3.6 million (very
long) rows of numpy data? I thought about using [one of
these](http://stackoverflow.com/questions/19304279/how-to-shuffle-a-text-file-
on-disk-in-python) techniques, but writing these arrays using `numpy.savetxt`
produces _unbelievably_ huge files, and I can't tell how to manipulate
individual rows from a standard `npy` file in a way that helps to solve this
problem.

Right now, my best idea is to create a permutation of paired indices `(c, r)`
into the data, where `c` choses a chunk and `r` choses a row from that chunk.
I could store each row in a new preallocated array, and then save it. But I
wonder if there's a less horribly I/O-bound solution. Is there some principled
way to shuffle random pairs of chunks together until you get a permutation
that's statistically independent from the starting permutation?


Answer: Among the things I've tried so far, a PyTables solution is currently the best,
followed by a solution that uses `numpy`'s support for memmapped arrays. The
PyTables solution is not straightforward though. If you use a shuffled array
of integers to directly index a PyTables array, it's very slow. Much faster is
the following two-step process:

  1. Select a random subset of the array using a boolean index array. _This must be done in a chunkwise fashion_. If you pass the index array directly to the PyTables array, it's slow. 
     * Preallocate a numpy array and create a list of slices that split the PyTables array into chunks. 
     * Read each chunk _entirely_ into memory, and then use the corresponding chunk of the index array to select the correct values for that chunk.
     * Store the selected values in the preallocated array. 
  2. _Then_ shuffle the preallocated array.

This process produces a permutation as random as a normal shuffling process
would. If that doesn't seem obvious, consider this: `(n choose x) * x! = x! *
n! / (x! * (n - x)!) = n! / (n - x)!`. This method is fast enough to do a
shuffle-on-load for every training cycle. It's also able to compress the data
down to ~650M -- nearly a 90% deflation.

Here's my current implementation; this is called once for every training chunk
in the corpus. (The returned arrays are shuffled elsewhere.)

    
    
    def _h5_fast_bool_ix(self, h5_array, ix, read_chunksize=100000):
        '''Iterate over an h5 array chunkwise to select a random subset
        of the array. `h5_array` should be the array itself; `ix` should
        be a boolean index array with as many values as `h5_array` has
        rows; and you can optionally set the number of rows to read per
        chunk with `read_chunksize` (default is 100000). For some reason
        this is much faster than using `ix` to index the array directly.'''
    
        n_chunks = h5_array.shape[0] / read_chunksize
        slices = [slice(i * read_chunksize, (i + 1) * read_chunksize)
                  for i in range(n_chunks)]
    
        a = numpy.empty((ix.sum(), h5_array.shape[1]), dtype=float)
        a_start = 0
        for sl in slices:
            chunk = h5_array[sl][ix[sl]]
            a_end = a_start + chunk.shape[0]
            a[a_start:a_end] = chunk
            a_start = a_end
    
        return a
    

It's somewhat crazy to me that an O(n^2) approach (iterating over the entire
PyTables array for every chunk) is faster in this case than an O(n) approach
(randomly selecting each row in one pass). But hey, it works. With a bit more
indirection, this could be adapted for loading arbitrary non-random
permutations, but that adds more complexity than it's worth here.

The `mmap` solution is here for reference, for those people who need a pure
numpy solution for whatever reason. It shuffles all the data in about 25
minutes, while the above solution manages the same in less than half that
time. This should scale linearly too, because `mmap` allows (relatively)
efficient random access.

    
    
    import numpy
    import os
    import random
    
    X = []
    Y = []
    
    for filename in os.listdir('input'):
        X.append(numpy.load(os.path.join('input', filename), mmap_mode='r'))
    
    for filename in os.listdir('output'):
        Y.append(numpy.load(os.path.join('output', filename), mmap_mode='r'))
    
    indices = [(chunk, row) for chunk, rows in enumerate(X) 
                            for row in range(rows.shape[0])]
    random.shuffle(indices)
    
    newchunks = 50
    newchunksize = len(indices) / newchunks
    
    for i in range(0, len(indices), newchunksize):
        print i
        rows = [X[chunk][row] for chunk, row in indices[i:i + newchunksize]]
        numpy.save('X_shuffled_' + str(i), numpy.array(rows))
        rows = [Y[chunk][row] for chunk, row in indices[i:i + newchunksize]]
        numpy.save('Y_shuffled_' + str(i), numpy.array(rows))
    

