Vaex tutorial

Dataset

Central to vaex is the dataset (similar, but not to be confused with a pandas dataframe, hence a different name), and we often use the variables ds to represent it. A dataset is an efficient representation for large tabular data, and has:

• A bunch of columns, say x, y and z
• Backed by a numpy array, e.g. ds.data.x (but you shouldn’t work with this)
• Wrapped by an expression system, e.g. ds.x, ds['x'] or ds.col.x is an expression
• Columns/expression can perform lazy computations, e.g. ds.x * np.sin(ds.y) does nothing, until the result is needed
• A set of virtual columns, columns that are backed by a (lazy) computation, e.g. ds['r'] = ds.x/ds.y
• A set of selection, that can be used to explore the dataset, e.g. ds.select(ds.x < 0)
• Filtered datasets, that does not copy the data, ds_negative = ds[ds.x < 0]

Lets start with an example dataset, included in vaex

import vaex
ds = vaex.example()
ds  # begin the last statement it will print out the tabular data

#	x	y	z	vx	vy	vz	E	L	Lz	FeH
0	-0.77747076699999995	2.1062629199999998	1.93743467	53.276721999999999	288.38604700000002	-95.264907800000003	-121238.171875	831.0799560546875	-336.426513671875	-2.3092276091645179
1	3.7742731599999999	2.2338719399999998	3.76209331	252.81079099999999	-69.949844400000003	-56.312103299999997	-100819.9140625	1435.1839599609375	-828.75677490234375	-1.788735491591229
2	1.3757626999999999	-6.3283844	2.6325001700000001	96.276473999999993	226.440201	-34.752716100000001	-100559.9609375	1039.2989501953125	920.802490234375	-0.76181090224787984
3	-7.0673780400000004	1.31737781	-6.1054353700000004	204.968842	-205.67901599999999	-58.977703099999999	-70174.8515625	2441.724853515625	1183.5899658203125	-1.5208778422936413
4	0.243441463	-0.82278168200000001	-0.20659387100000001	-311.74237099999999	-238.41217	186.824127	-144138.75	374.81643676757812	-314.53530883789062	-2.6553413584273611
...	...	...	...	...	...	...	...	...	...	...
329,995	3.7688379300000001	4.6625165900000001	-4.4290413900000001	107.432999	-2.13771296	17.5130272	-119687.3203125	746.88336181640625	-508.96484375	-1.6499842518381402
329,996	9.1740932500000003	-8.8709135099999994	-8.61707687	32.0	108.089264	179.06063800000001	-68933.8046875	2395.633056640625	1275.490234375	-1.4336036247720836
329,997	-1.1404100699999999	-8.4957694999999998	2.2574982600000002	8.4671134899999991	-38.276523599999997	-127.541473	-112580.359375	1182.436279296875	115.58557891845703	-1.9306227597361942
329,998	-14.298593500000001	-5.5175042200000002	-8.6547231700000005	110.221558	-31.3925591	86.272682200000006	-74862.90625	1324.5926513671875	1057.017333984375	-1.2250198188385679
329,999	10.5450506	-8.86106777	-4.6583542800000002	-2.1054141500000001	-27.6108856	3.80799961	-95361.765625	351.09555053710938	-309.81439208984375	-2.5689636894079477

Columns

The above preview shows this dataset contains \(> 300,000\) rows, and columns named x,y,z (positions), vx, vy, vz (velocities), E (energy), L (angular momentum). Printing out a column, shows it is not a numpy array, but an expression

ds.x  # ds.col.x or ds['x'] are equivalent, but may be preferred because it is more tab completion friend or programmatics friendly respectively

<vaex.expression.Expression(expressions='x')> instance at 0x10f3d82e8 [-0.777470767, 3.77427316, 1.3757627, -7.06737804, 0.243441463 ... (total 330000 values) ... 3.76883793, 9.17409325, -1.14041007, -14.2985935, 10.5450506]

The underlying data is often accessible using ds.data.x, but should not be used, since selections and filtering are not reflected in this. However sometimes it is useful to access the raw numpy array.

ds.data.x

array([ -0.77747077,   3.77427316,   1.3757627 , ...,  -1.14041007,
       -14.2985935 ,  10.5450506 ])

A better way, if you need a numpy array (for instance for plotting, or passing to a different library) it to use evalulate, which will also work with virtual columns, selections and filtered datasets (more on that below).

ds.evaluate(ds.x)

array([ -0.77747077,   3.77427316,   1.3757627 , ...,  -1.14041007,
       -14.2985935 ,  10.5450506 ])

Most numpy function (ufuncs) can be performed on expressions, and will not result in a direct result, but in a new expression.

import numpy as np
np.sqrt(ds.x**2 + ds.y**2 + ds.z**2)

<vaex.expression.Expression(expressions='sqrt((((x ** 2) + (y ** 2)) + (z ** 2)))')> instance at 0x10f3d8a58 [2.96554503966, 5.77829281049, 6.9907960395, 9.43184275271, 0.882561312135 ... (total 330000 values) ... 7.45383176151, 15.3984124911, 8.86425027393, 17.601047186, 14.540181525]

Virtual functions

Sometimes it is convenient to store an expression as a column, or virtual column, a column that does not take up memory, but will be computed on the fly. A virtual column can be treated as a normal column.

ds['r'] = np.sqrt(ds.x**2 + ds.y**2 + ds.z**2)
ds[['x', 'y', 'z', 'r']]

#	x	y	z	r
0	-0.77747076699999995	2.1062629199999998	1.93743467	2.9655450396553587
1	3.7742731599999999	2.2338719399999998	3.76209331	5.7782928104901803
2	1.3757626999999999	-6.3283844	2.6325001700000001	6.9907960395025599
3	-7.0673780400000004	1.31737781	-6.1054353700000004	9.4318427527075368
4	0.243441463	-0.82278168200000001	-0.20659387100000001	0.88256131213479672
...	...	...	...	...
329,995	3.7688379300000001	4.6625165900000001	-4.4290413900000001	7.4538317615146807
329,996	9.1740932500000003	-8.8709135099999994	-8.61707687	15.398412491068198
329,997	-1.1404100699999999	-8.4957694999999998	2.2574982600000002	8.8642502739256326
329,998	-14.298593500000001	-5.5175042200000002	-8.6547231700000005	17.601047186042507
329,999	10.5450506	-8.86106777	-4.6583542800000002	14.540181524970293

Selections and filtering

Vaex can be efficient when exploring subsets of the data, for instance to remove outlier or to inspect only a part of the data. Instead of making copies, internally vaex keeps track which rows is selected.

ds.select(ds.x < 0)
ds.evaluate(ds.x, selection=True)

array([ -0.77747077,  -7.06737804,  -5.17174435, ...,  -1.87310386,
        -1.14041007, -14.2985935 ])

Selections can be useful if you want to change what you select frequently, as in visualization, or when you want to compute statistics on several selections efficiently. Instead, you can also create a filtered dataset, and is similar in use to pandas, except that it does not copy the data.

ds_negative = ds[ds.x < 0]
ds_negative[['x', 'y', 'z', 'r']]

#	x	y	z	r
0	-0.77747076699999995	2.1062629199999998	1.93743467	2.9655450396553587
1	-7.0673780400000004	1.31737781	-6.1054353700000004	9.4318427527075368
2	-5.17174435	7.8291530600000003	1.8266882900000001	9.5592555864715436
3	-15.953885100000001	5.7712588299999998	-9.0247230500000004	19.21664654397474
4	-12.3994961	13.9181805	-5.4348230400000004	19.416502090763164
...	...	...	...	...
165,935	-9.8855323800000008	-6.5925359700000001	6.5374202700000001	13.561826747838182
165,936	-2.38018084	4.7354030600000003	0.14176586299999999	5.3018299229296861
165,937	-1.8731038600000001	-0.50309121599999995	-0.95197701499999998	2.1605275001840565
165,938	-1.1404100699999999	-8.4957694999999998	2.2574982600000002	8.8642502739256326
165,939	-14.298593500000001	-5.5175042200000002	-8.6547231700000005	17.601047186042507

Statistics on N-d grids

A core feature of vaex, and used for visualization, is calculation of statistics on N dimensional grids.

ds.count(), ds.mean(ds.x), ds.mean(ds.x, selection=True)

(330000.0, -0.067131491264005971, -5.2110379721119671)

Similar to SQL’s groupby, vaex uses the binby concept, which tells vaex that a statistic should be calculated on a regular grid (for performance reasons)

xcounts = ds.count(binby=ds.x, limits=[-10, 10], shape=64)
xcounts

array([ 1310.,  1416.,  1452.,  1519.,  1599.,  1810.,  1956.,  2005.,
        2157.,  2357.,  2653.,  2786.,  3012.,  3215.,  3619.,  3890.,
        3973.,  4400.,  4782.,  5126.,  5302.,  5729.,  6042.,  6562.,
        6852.,  7167.,  7456.,  7633.,  7910.,  8415.,  8619.,  8246.,
        8358.,  8769.,  8294.,  7870.,  7749.,  7389.,  7174.,  6901.,
        6557.,  6173.,  5721.,  5367.,  4963.,  4655.,  4246.,  4110.,
        3939.,  3611.,  3289.,  3018.,  2811.,  2570.,  2505.,  2267.,
        2013.,  1803.,  1687.,  1563.,  1384.,  1326.,  1257.,  1189.])

This results in a numpy array with the number counts in 64 bins distributed between x = -10, and x = 10. We can quickly visualize this using matplotlib.

import matplotlib.pylab as plt
plt.plot(np.linspace(-10, 10, 64), xcounts)
plt.show()

We can instead of doing 1d binning, do it in 2d as well (N-d actually), and visualize it using imshow.

xycounts = ds.count(binby=[ds.x, ds.y], limits=[[-10, 10], [-10, 20]], shape=(64, 128))
xycounts

array([[  9.,   3.,   3., ...,   3.,   2.,   1.],
       [  5.,   3.,   1., ...,   1.,   3.,   3.],
       [ 11.,   3.,   2., ...,   1.,   1.,   4.],
       ...,
       [ 12.,   6.,   8., ...,   0.,   1.,   0.],
       [  7.,   6.,  12., ...,   3.,   0.,   0.],
       [ 11.,  10.,   7., ...,   1.,   1.,   1.]])

plt.imshow(xycounts.T, origin='lower', extent=[-10, 10, -10, 20])
plt.show()

v = np.sqrt(ds.vx**2 + ds.vy**2 + ds.vz**2)
xy_mean_v = ds.mean(v, binby=[ds.x, ds.y], limits=[[-10, 10], [-10, 20]], shape=(64, 128))
xy_mean_v

array([[ 144.38495511,  183.45775869,  187.78325557, ...,  138.99392387,
         168.66141282,  142.55018784],
       [ 143.72427758,  152.14679337,  107.90949865, ...,  119.65318885,
          94.00098292,  104.35109636],
       [ 172.08240652,  137.47896886,   72.51331138, ...,  179.85933835,
          33.36968912,  111.81826254],
       ...,
       [ 186.56949934,  161.3747346 ,  174.27411865, ...,           nan,
         105.96746091,           nan],
       [ 179.55997022,  137.48979882,  113.82121826, ...,  104.90205692,
                  nan,           nan],
       [ 151.94323763,  135.44083212,   84.81787495, ...,  175.79289144,
         129.63799565,  108.19069385]])

plt.imshow(xy_mean_v.T, origin='lower', extent=[-10, 10, -10, 20])
plt.show()

Other statistics can be computed, such as:

• Dataset.count
• Dataset.mean
• Dataset.std
• Dataset.var
• Dataset.median_approx
• Dataset.percentile_approx
• Dataset.mode
• Dataset.min
• Dataset.max
• Dataset.minmax
• Dataset.mutual_information
• Dataset.correlation

Or see the full list at the API docs

Getting your data in

Before continuing, you may want to read in your own data. Ultimately, a vaex Dataset just wraps a set of numpy arrays. If you can access your data as a set of numpy arrays, you can therefore make dataset using from_arrays

import vaex
import numpy as np
x = np.arange(5)
y = x**2
ds = vaex.from_arrays(x=x, y=y)
ds

#	y	x
0	0	0
1	1	1
2	4	2
3	9	3
4	16	4

Other quick ways to get your data in are:

• from_csv: Comma separated files
• from_ascii: Space/tab separated files
• from_pandas: Converts a pandas DataFrame
• from_astropy_table: Converts a astropy table

Plotting

1d and 2d

Most visualization can be done in 1 and 2d, and vaex wraps matplotlib to provide most use cases.

import vaex
import numpy as np
ds = vaex.example()
%matplotlib inline

The simpelest visualization is a 1d plot using Dataset.plot1d. When only given one arguments, it will show a histogram showing 99.8% of the data.

ds.plot1d(ds.x)

[<matplotlib.lines.Line2D at 0x11375f630>]

A slighly more complication visualization, is to not plot the counts, but a different statistic for that bin. In most cases, passing the what='<statistic>(<expression>) argument will do, where <statistic> is any of the statistics mentioned in the list above, or in the API docs

ds.plot1d(ds.x, what='mean(E)')

[<matplotlib.lines.Line2D at 0x1134a90f0>]

An equivalent method, is to use the vaex.stat.<statistic> functions, e.g. vaex.stat.mean

ds.plot1d(ds.x, what=vaex.stat.mean(ds.E))

[<matplotlib.lines.Line2D at 0x1134d7128>]

These objects are very similar to vaex’ expression, in that they represent an underlying calculation, while normal arithmetic and numpy functions can be applied to it. However, these object represent a statistics computation, and not a column.

np.log(vaex.stat.mean(ds.x)/vaex.stat.std(ds.x))

log((mean(x) / std(x)))

These statistical objects can be passed to the what argument. The advantage being that the data will only have to be passed over once.

ds.plot1d(ds.x, what=np.clip(np.log(-vaex.stat.mean(ds.E)), 11, 11.4))

[<matplotlib.lines.Line2D at 0x1134fce10>]

A similar result can be obtained by calculating the statistic ourselves, and passing it to plot1d’s grid argument. Care has to be taken that the limits used for calculating the statistics and the plot are the same, otherwise the x axis may not correspond to the real data.

limits = [-30, 30]
shape  = 64
meanE  = ds.mean(ds.E, binby=ds.x, limits=limits, shape=shape)
grid   = np.clip(np.log(-meanE), 11, 11.4)
ds.plot1d(ds.x, grid=grid, limits=limits, ylabel='clipped E')

[<matplotlib.lines.Line2D at 0x11822bc88>]

The same applies for 2d plotting.

ds.plot(ds.x, ds.y, what=vaex.stat.mean(ds.E)**2)

<matplotlib.image.AxesImage at 0x116034278>

Selections for plotting

While filtering is useful for narrowing down a selection (e.g. ds_negative = ds[ds.x < 0]) there are a few downsides to this. First, a practical issue is that when you filter 4 different ways, you will need to have 4 different objects, polluting your namespace. However, more importantly, when vaex executes a bunch of statistical computations, it will do that per dataset, meaning for 4 different datasets (although pointing to the same underlying data) it will do a total of 4 passes over the data. If instead, we have 4 (named) selections in our dataset, it can calculate statistics in one single pass over the data, which can speed up especially when you dataset is larger than your memory.

In the plot below, we show three selection, which by default are blended together, requiring just one pass over the data.

ds.plot(ds.x, ds.y, what=np.log(vaex.stat.count()+1),
       selection=[None, ds.x < ds.y, ds.x < -10])

<matplotlib.image.AxesImage at 0x112648198>

Advanced Plotting

Lets say we would like to see two plots next to eachother, we can pass a list of expression pairs.

ds.plot([["x", "y"], ["x", "z"]],
        title="Face on and edge on", figsize=(10,4));

By default, if you have multiple plots, they are shows as columns, multiple selections are overplotted, and multiple ‘whats’ (statistics) are shows as rows.

ds.plot([["x", "y"], ["x", "z"]],
        what=[np.log(vaex.stat.count()+1), vaex.stat.mean(ds.E)],
        selection=[None, ds.x < ds.y],
        title="Face on and edge on", figsize=(10,10));

(Note that the selection has no effect in the bottom rows)

However, this behaviour can be changed using the visual argument.

ds.plot([["x", "y"], ["x", "z"]],
        what=vaex.stat.mean(ds.E),
        selection=[None, ds.Lz < 0],
        visual=dict(column='selection'),
        title="Face on and edge on", figsize=(10,10));

Slices in a 3rd dimension

If a 3rd axis (z) is given, you can ‘slice’ through the data, displaying the z slices as rows. Note that here the rows are wrapped, which can be changed using the wrap_columns argument.

ds.plot("Lz", "E", z="FeH:-3,-1,10", show=True, visual=dict(row="z"),
        figsize=(12,8), f="log", wrap_columns=3);

Smaller datasets / scatter plot

Although vaex focusses on large datasets, sometimes you end up with a fraction of the data (due to a selection) and you want to make a scatter plot. You could try the following approach:

import vaex
ds = vaex.example()
%matplotlib inline

import matplotlib.pylab as plt
x = ds.evaluate("x", selection=ds.Lz < -2500)
y = ds.evaluate("y", selection=ds.Lz < -2500)
plt.scatter(x, y, c="red", alpha=0.5, s=4);

ds.scatter(ds.x, ds.y, selection=ds.Lz < -2500, c="red", alpha=0.5, s=4)
ds.scatter(ds.x, ds.y, selection=ds.Lz > 1500, c="green", alpha=0.5, s=4);

In control

While vaex provides a wrapper for matplotlib, there are situations where you want to use the Dataset.plot method, but want to be in control of the plot. Vaex simply uses the current figure and axes, so that is easy to do

import numpy as np

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14,7))
plt.sca(ax1)
selection = ds.Lz < -2500
x = ds[selection].x.evaluate()#selection=selection)
y = ds[selection].y.evaluate()#selection=selection)
ds.plot(ds.x, ds.y)
plt.scatter(x, y)
plt.xlabel('my own label $\gamma$')
plt.xlim(-20, 20)
plt.ylim(-20, 20)

plt.sca(ax2)
ds.plot1d(ds.x, label='counts', n=True)
x = np.linspace(-30, 30, 100)
std = ds.std(ds.x.expression)
y = np.exp(-(x**2/std**2/2)) / np.sqrt(2*np.pi) / std
plt.plot(x, y, label='gaussian fit')
plt.legend()

<matplotlib.legend.Legend at 0x1193bf438>

Healpix (Plotting)

Using healpix is made available by the vaex-healpix package using the healpy package. Vaex does not need special support for healpix, only for plotting, but some helper functions are introduced to make working with healpix easier. By diving the source_id by 34359738368 you get a healpix index level 12, and diving it further will take you to lower levels.

To understand healpix better, we will start from the beginning. If we want to make a density sky plot, we would like to pass healpy a 1d numpy array where each value represents the density at a location of the sphere, where the location is determined by the array size (the healpix level) and the offset (the location). Since the Gaia data includes the healpix index encoded in the source_id. By diving the source_id by 34359738368 you get a healpix index level 12, and diving it further will take you to lower levels.

import vaex
import healpy as hp
%matplotlib inline
tgas = vaex.datasets.tgas.fetch()

We will start showing how you could manually do statistics on healpix bins using vaex.count. We will do a really course healpix scheme (level 2).

level = 2
factor = 34359738368 * (4**(12-level))
nmax = hp.nside2npix(2**level)
epsilon = 1e-16
counts = tgas.count(binby=tgas.source_id/factor, limits=[-epsilon, nmax-epsilon], shape=nmax)
counts

array([  4021.,   6171.,   5318.,   7114.,   5755.,  13420.,  12711.,
        10193.,   7782.,  14187.,  12578.,  22038.,  17313.,  13064.,
        17298.,  11887.,   3859.,   3488.,   9036.,   5533.,   4007.,
         3899.,   4884.,   5664.,  10741.,   7678.,  12092.,  10182.,
         6652.,   6793.,  10117.,   9614.,   3727.,   5849.,   4028.,
         5505.,   8462.,  10059.,   6581.,   8282.,   4757.,   5116.,
         4578.,   5452.,   6023.,   8340.,   6440.,   8623.,   7308.,
         6197.,  21271.,  23176.,  12975.,  17138.,  26783.,  30575.,
        31931.,  29697.,  17986.,  16987.,  19802.,  15632.,  14273.,
        10594.,   4807.,   4551.,   4028.,   4357.,   4067.,   4206.,
         3505.,   4137.,   3311.,   3582.,   3586.,   4218.,   4529.,
         4360.,   6767.,   7579.,  14462.,  24291.,  10638.,  11250.,
        29619.,   9678.,  23322.,  18205.,   7625.,   9891.,   5423.,
         5808.,  14438.,  17251.,   7833.,  15226.,   7123.,   3708.,
         6135.,   4110.,   3587.,   3222.,   3074.,   3941.,   3846.,
         3402.,   3564.,   3425.,   4125.,   4026.,   3689.,   4084.,
        16617.,  13577.,   6911.,   4837.,  13553.,  10074.,   9534.,
        20824.,   4976.,   6707.,   5396.,   8366.,  13494.,  19766.,
        11012.,  16130.,   8521.,   8245.,   6871.,   5977.,   8789.,
        10016.,   6517.,   8019.,   6122.,   5465.,   5414.,   4934.,
         5788.,   6139.,   4310.,   4144.,  11437.,  30731.,  13741.,
        27285.,  40227.,  16320.,  23039.,  10812.,  14686.,  27690.,
        15155.,  32701.,  18780.,   5895.,  23348.,   6081.,  17050.,
        28498.,  35232.,  26223.,  22341.,  15867.,  17688.,   8580.,
        24895.,  13027.,  11223.,   7880.,   8386.,   6988.,   5815.,
         4717.,   9088.,   8283.,  12059.,   9161.,   6952.,   4914.,
         6652.,   4666.,  12014.,  10703.,  16518.,  10270.,   6724.,
         4553.,   9282.,   4981.])

And using healpy’s mollview we can visualize this.

hp.mollview(counts, nest=True)

To simplify life, vaex includes Dataset.healpix_count to take care of this.

counts = tgas.healpix_count(healpix_level=6)
hp.mollview(counts, nest=True)

Or even simpler, use Dataset.healpix_plot

tgas.healpix_plot(f="log1p", healpix_level=6, figsize=(10,8),
                  healpix_output="ecliptic")

Parallel computations

As mentioned in the sections on selections, vaex can do computations on a dataset in parallel. Often, this is taken care of, when for instance passing multiple selections, or multiple arguments to one of the statistical functions. However, sometimes it is difficult or impossible to express a computation in one expression, and we need to resort to doing so called ‘delayed’ computationed, similar as in joblib and dask.

import vaex
ds = vaex.example()
limits = [-10, 10]
delayed_count = ds.count(ds.E, binby=ds.x, limits=limits,
                         shape=4, delay=True)
delayed_count

<vaex.promise.Promise at 0x1164989b0>

Note that now the returned value is not a promise (TODO: a more Pythonic way would be to return a Future). This may be subject to change, and the best way to work with this is to use the delayed decorator. And call Dataset.execute when the result is needed.

In addition to the above delayed computation, we schedule another computation, such that both the count and mean are execute in parallel such that we only do a single pass over the data. We schedule the execution of two extra functions using the vaex.delayed decorator, and run the whole pipeline using ds.execute().

delayed_sum = ds.sum(ds.E, binby=ds.x, limits=limits,
                         shape=4, delay=True)

@vaex.delayed
def calculate_mean(sums, counts):
    print('calculating mean')
    return sums/counts

print('before calling mean')
# since calculate_mean is decorator with vaex.delayed
# this now also returns a 'delayed' object (a promise)
delayed_mean = calculate_mean(delayed_sum, delayed_count)

# if we'd like to perform operations on that, we can again
# use the same decorator
@vaex.delayed
def print_mean(means):
    print('means', means)
print_mean(delayed_mean)

print('before calling execute')
ds.execute()

# Using the .get on the promise will also return the resut
# However, this will only work after execute, and may be
# subject to change
means = delayed_mean.get()
print('same means', means)

before calling mean
before calling execute
calculating mean
means [ -94415.16581227 -118856.63989386 -118919.86423543  -95000.5998913 ]
same means [ -94415.16581227 -118856.63989386 -118919.86423543  -95000.5998913 ]

Interactive widgets

Note: The interactive widgets require a running Python kernel, if you are viewing this documentation online you mean get a feeling for what the widgets can do, but computation will not be possible!

Using the vaex-jupyter package, we get access to interactive widgets.

import vaex
import vaex.jupyter
import numpy as np
import pylab as plt
%matplotlib inline
ds = vaex.example()

The simplest way to get a more interactive visualization (or even print out statistics) is to the the vaex.jupyter.interactive_selection decorator, which will execute the decorated function each time the selection is changed.

ds.select(ds.x > 0)
@vaex.jupyter.interactive_selection(ds)
def plot():
    print("Mean x for the selection is:", ds.mean(ds.x, selection=True))
    ds.plot(ds.x, ds.y, what=np.log(vaex.stat.count()+1), selection=[None, True])
    plt.show()

After changing the selection programmatically, the visualization will update, as well as the print output.

ds.select(ds.x > ds.y)

However, to get truly interactive visualization, we need to use widgets, such as the bqplot library. Again, if we make a selection here, the above visualization will also update, so lets select a square region. One issue is that if you have installed ipywidget, bqplot, ipyvolume etc, it may not be enabled if you installed them from pip (from conda-forge will enabled it automagically). To enable it, run the next cell, and refresh the notebook if there were not enabled already. (Note that these commands will execute in the environment where the notebook is running, not where the kernel is running)

import sys
!jupyter nbextension enable --sys-prefix --py widgetsnbextension
!jupyter nbextension enable --sys-prefix --py bqplot
!jupyter nbextension enable --sys-prefix --py ipyvolume
!jupyter nbextension enable --sys-prefix --py ipympl
!jupyter nbextension enable --sys-prefix --py ipyleaflet

# the default backend is bqplot, but we pass it here explicity
ds.plot_widget(ds.x, ds.y, f='log', backend='bqplot')

Joining

Joining in vaex is similar to pandas, except the data will no be copied. Internally an index array is kept for each row on the left dataset, pointing to the right dataset, requiring about 8GB for a billion row \(10^9\) dataset. Lets start with 2 small dataset, ds1 and ds2:

a = np.array(['a', 'b', 'c'])
x = np.arange(1,4)
ds1 = vaex.from_arrays(a=a, x=x)
ds1

#	a	x
0	'a'	1
1	'b'	2
2	'c'	3

b = np.array(['a', 'b', 'd'])
y = x**2
ds2 = vaex.from_arrays(b=b, y=y)
ds2

#	y	b
0	1	'a'
1	4	'b'
2	9	'd'

The default join, is a ‘left’ join, where all rows for the left dataset (ds1) are kept, and matching rows of the right dataset (ds2) are added. We see for for the columns b and y, some values are missing, as expected.

ds1.join(ds2, left_on='a', right_on='b')

#	a	x	y	b
0	'a'	1	1	'a'
1	'b'	2	4	'b'
2	'c'	3	masked	masked

A ‘right’ join, is basically the same, but now the roles of the left and right label swapped, so now we have some values from columns x and a missing.

ds1.join(ds2, left_on='a', right_on='b', how='right')

#	y	b	a	x
0	1	'a'	'a'	1
1	4	'b'	'b'	2
2	9	'd'	masked	masked

Other joins (inner and outer) aren’t supported, feel free open an issue on github for this.

Just-In-Time compilation

Lets start with a function that converts from two angles, to an angular distance. The function assumes as input, 2 pairs on angular coordinates, in radians.

import vaex
import numpy as np
# From http://pythonhosted.org/pythran/MANUAL.html
def arc_distance(theta_1, phi_1, theta_2, phi_2):
    """
    Calculates the pairwise arc distance
    between all points in vector a and b.
    """
    temp = (np.sin((theta_2-2-theta_1)/2)**2
           + np.cos(theta_1)*np.cos(theta_2) * np.sin((phi_2-phi_1)/2)**2)
    distance_matrix = 2 * np.arctan2(np.sqrt(temp), np.sqrt(1-temp))
    return distance_matrix

Let us use the New York Taxi dataset of 2015, as can be downloaded in hdf5 format

nytaxi = vaex.open("/Users/maartenbreddels/vaex/data/nytaxi/nyc_taxi2015.hdf5")
# lets use just 20% of the data, since we want to make sure it fits
# into memory (so we don't measure just hdd/ssd speed)
nytaxi.set_active_fraction(0.2)

Although the function above expected numpy arrays, vaex can pass in columns or expression, which will delay execution till needed, and add the resulting expression as a virtual column.

nytaxi['arc_distance'] = arc_distance(nytaxi.pickup_longitude * np.pi/180,
                                      nytaxi.pickup_latitude * np.pi/180,
                                      nytaxi.dropoff_longitude * np.pi/180,
                                      nytaxi.dropoff_latitude * np.pi/180)

When we calculate the mean angular distance of a taxi trip, we encounter some invalid data, that will give warnings, which we can savely ignore for this demonstration.

%%time
nytaxi.mean(nytaxi.arc_distance)

CPU times: user 10.1 s, sys: 3.77 s, total: 13.8 s
Wall time: 5.1 s

1.9999877196037044

This computation uses quite some heavy mathematical operation, and since it’s (internally) using numpy arrays, also uses quite some temporary arrays. We can optimize this calculation by doing a Just-In-Time compilation, based on numba or pythran. Choose whichever gives the best performance or is easiest to install.

nytaxi['arc_distance_jit'] = nytaxi.arc_distance.jit_numba()
# nytaxi['arc_distance_jit'] = nytaxi.arc_distance.jit_pythran()

%%time
nytaxi.mean(nytaxi.arc_distance_jit)

CPU times: user 4.4 s, sys: 15.6 ms, total: 4.42 s
Wall time: 1.19 s

1.9999877196037033

We can that we can get a significant speedup (\(\approx 3 x\)) in this case.