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

In [1]:
import vaex
ds = vaex.example()
ds  # begin the last statement it will print out the tabular data
Out[1]:
#xyzvxvyvzELLzFeH
0-0.777470766999999952.10626291999999981.9374346753.276721999999999288.38604700000002-95.264907800000003-121238.171875831.0799560546875-336.426513671875-2.3092276091645179
13.77427315999999992.23387193999999983.76209331252.81079099999999-69.949844400000003-56.312103299999997-100819.91406251435.1839599609375-828.75677490234375-1.788735491591229
21.3757626999999999-6.32838442.632500170000000196.276473999999993226.440201-34.752716100000001-100559.96093751039.2989501953125920.802490234375-0.76181090224787984
3-7.06737804000000041.31737781-6.1054353700000004204.968842-205.67901599999999-58.977703099999999-70174.85156252441.7248535156251183.5899658203125-1.5208778422936413
40.243441463-0.82278168200000001-0.20659387100000001-311.74237099999999-238.41217186.824127-144138.75374.81643676757812-314.53530883789062-2.6553413584273611
.................................
329,9953.76883793000000014.6625165900000001-4.4290413900000001107.432999-2.1377129617.5130272-119687.3203125746.88336181640625-508.96484375-1.6499842518381402
329,9969.1740932500000003-8.8709135099999994-8.6170768732.0108.089264179.06063800000001-68933.80468752395.6330566406251275.490234375-1.4336036247720836
329,997-1.1404100699999999-8.49576949999999982.25749826000000028.4671134899999991-38.276523599999997-127.541473-112580.3593751182.436279296875115.58557891845703-1.9306227597361942
329,998-14.298593500000001-5.5175042200000002-8.6547231700000005110.221558-31.392559186.272682200000006-74862.906251324.59265136718751057.017333984375-1.2250198188385679
329,99910.5450506-8.86106777-4.6583542800000002-2.1054141500000001-27.61088563.80799961-95361.765625351.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

In [2]:
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
Out[2]:
<vaex.expression.Expression(expressions='x')> instance at 0x7f21c16d3ac8 values=[-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.

In [3]:
ds.data.x
Out[3]:
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).

In [4]:
ds.evaluate(ds.x)
Out[4]:
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.

In [5]:
import numpy as np
np.sqrt(ds.x**2 + ds.y**2 + ds.z**2)
Out[5]:
<vaex.expression.Expression(expressions='sqrt((((x ** 2) + (y ** 2)) + (z ** 2)))')> instance at 0x7f21978f42e8 values=[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.

In [6]:
ds['r'] = np.sqrt(ds.x**2 + ds.y**2 + ds.z**2)
ds[['x', 'y', 'z', 'r']]
Out[6]:
#xyzr
0-0.777470766999999952.10626291999999981.937434672.9655450396553587
13.77427315999999992.23387193999999983.762093315.7782928104901803
21.3757626999999999-6.32838442.63250017000000016.9907960395025599
3-7.06737804000000041.31737781-6.10543537000000049.4318427527075368
40.243441463-0.82278168200000001-0.206593871000000010.88256131213479672
...............
329,9953.76883793000000014.6625165900000001-4.42904139000000017.4538317615146807
329,9969.1740932500000003-8.8709135099999994-8.6170768715.398412491068198
329,997-1.1404100699999999-8.49576949999999982.25749826000000028.8642502739256326
329,998-14.298593500000001-5.5175042200000002-8.654723170000000517.601047186042507
329,99910.5450506-8.86106777-4.658354280000000214.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.

In [7]:
ds.select(ds.x < 0)
ds.evaluate(ds.x, selection=True)
Out[7]:
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.

In [8]:
ds_negative = ds[ds.x < 0]
ds_negative[['x', 'y', 'z', 'r']]
Out[8]:
#xyzr
0-0.777470766999999952.10626291999999981.937434672.9655450396553587
1-7.06737804000000041.31737781-6.10543537000000049.4318427527075368
2-5.171744357.82915306000000031.82668829000000019.5592555864715436
3-15.9538851000000015.7712588299999998-9.024723050000000419.21664654397474
4-12.399496113.9181805-5.434823040000000419.416502090763164
...............
165,935-9.8855323800000008-6.59253597000000016.537420270000000113.561826747838182
165,936-2.380180844.73540306000000030.141765862999999995.3018299229296861
165,937-1.8731038600000001-0.50309121599999995-0.951977014999999982.1605275001840565
165,938-1.1404100699999999-8.49576949999999982.25749826000000028.8642502739256326
165,939-14.298593500000001-5.5175042200000002-8.654723170000000517.601047186042507

Statistics on N-d grids

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

In [9]:
ds.count(), ds.mean(ds.x), ds.mean(ds.x, selection=True)
Out[9]:
(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)

In [10]:
xcounts = ds.count(binby=ds.x, limits=[-10, 10], shape=64)
xcounts
Out[10]:
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.

In [12]:
import matplotlib.pylab as plt
plt.plot(np.linspace(-10, 10, 64), xcounts)
plt.show()
_images/tutorial_22_0.png

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

In [13]:
xycounts = ds.count(binby=[ds.x, ds.y], limits=[[-10, 10], [-10, 20]], shape=(64, 128))
xycounts
Out[13]:
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.]])
In [14]:
plt.imshow(xycounts.T, origin='lower', extent=[-10, 10, -10, 20])
plt.show()
_images/tutorial_25_0.png
In [15]:
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
Out[15]:
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]])
In [16]:
plt.imshow(xy_mean_v.T, origin='lower', extent=[-10, 10, -10, 20])
plt.show()
_images/tutorial_27_0.png

Other statistics can be computed, such as:

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

In [17]:
import vaex
import numpy as np
x = np.arange(5)
y = x**2
ds = vaex.from_arrays(x=x, y=y)
ds
Out[17]:
#xy
000
111
224
339
4416

Other quick ways to get your data in are:

Plotting

1d and 2d

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

In [18]:
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.

In [19]:
ds.plot1d(ds.x)
Out[19]:
[<matplotlib.lines.Line2D at 0x7f2151713a58>]
_images/tutorial_35_1.png

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

In [20]:
ds.plot1d(ds.x, what='mean(E)')
Out[20]:
[<matplotlib.lines.Line2D at 0x7f215167e390>]
_images/tutorial_37_1.png

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

In [21]:
ds.plot1d(ds.x, what=vaex.stat.mean(ds.E))
Out[21]:
[<matplotlib.lines.Line2D at 0x7f215056fba8>]
_images/tutorial_39_1.png

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.

In [22]:
np.log(vaex.stat.mean(ds.x)/vaex.stat.std(ds.x))
Out[22]:
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.

In [23]:
ds.plot1d(ds.x, what=np.clip(np.log(-vaex.stat.mean(ds.E)), 11, 11.4))
Out[23]:
[<matplotlib.lines.Line2D at 0x7f2150559908>]
_images/tutorial_43_1.png

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.

In [24]:
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')
Out[24]:
[<matplotlib.lines.Line2D at 0x7f21504c3a58>]
_images/tutorial_45_1.png

The same applies for 2d plotting.

In [25]:
ds.plot(ds.x, ds.y, what=vaex.stat.mean(ds.E)**2)
Out[25]:
<matplotlib.image.AxesImage at 0x7f21515cec18>
_images/tutorial_47_1.png

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.

In [26]:
ds.plot(ds.x, ds.y, what=np.log(vaex.stat.count()+1),
       selection=[None, ds.x < ds.y, ds.x < -10])
Out[26]:
<matplotlib.image.AxesImage at 0x7f21515c7cf8>
_images/tutorial_49_1.png

Advanced Plotting

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

In [27]:
ds.plot([["x", "y"], ["x", "z"]],
        title="Face on and edge on", figsize=(10,4));
/net/gaia/data/users/breddels/python/anaconda3/envs/dev/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
_images/tutorial_51_1.png

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

In [28]:
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));
_images/tutorial_53_0.png

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

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

In [29]:
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));
/net/gaia/data/users/breddels/python/anaconda3/envs/dev/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
_images/tutorial_55_1.png

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.

In [30]:
ds.plot("Lz", "E", z="FeH:-3,-1,10", show=True, visual=dict(row="z"),
        figsize=(12,8), f="log", wrap_columns=3);
/net/gaia/data/users/breddels/python/anaconda3/envs/dev/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
_images/tutorial_57_1.png

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:

In [31]:
import vaex
ds = vaex.example()
%matplotlib inline
In [32]:
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);
_images/tutorial_60_0.png
In [33]:
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);
_images/tutorial_61_0.png

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

In [34]:
import numpy as np
In [35]:
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()
Out[35]:
<matplotlib.legend.Legend at 0x7f2150316208>
_images/tutorial_64_1.png

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.

In [36]:
import vaex
import healpy as hp
%matplotlib inline
tgas = vaex.datasets.tgas.fetch()
Downloading http://vaex.astro.rug.nl/data/tgas.hdf5 to /Users/users/breddels/.vaex/data/tgas.hdf5
wget failed, using urlretrieve

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

In [37]:
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
Out[37]:
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.

In [38]:
hp.mollview(counts, nest=True)
_images/tutorial_70_0.png

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

In [39]:
counts = tgas.healpix_count(healpix_level=6)
hp.mollview(counts, nest=True)
_images/tutorial_72_0.png

Or even simpler, use Dataset.healpix_plot

In [40]:
tgas.healpix_plot(f="log1p", healpix_level=6, figsize=(10,8),
                  healpix_output="ecliptic")
_images/tutorial_74_0.png

Propagation of uncertainties

In science we often deal with measurement uncertainties (sometimes refererred to as measurement errors). When transformations are made with quantities that have uncertainties associated with them, the uncertainties on these transformed quantities can be calculated automatically by vaex. Note that propagation of uncertainties requires derivatives and matrix multiplications of lengthy equations, which is not complex, but tedious. Vaex can automatically calculate all dependencies, derivatives and compute the full covariance matrix.

In [41]:
import vaex
import pylab as plt
%matplotlib inline
tgas = vaex.datasets.tgas_1percent.fetch()
Downloading http://vaex.astro.rug.nl/data/tgas_1percent.hdf5 to /Users/users/breddels/.vaex/data/tgas_1percent.hdf5
wget failed, using urlretrieve

Even though the TGAS dataset already contains galactic sky coordiantes (l and b), we add them again as virtual columns such that the transformation between RA. and Dec. and the galactic sky coordinates is know.

In [42]:
# convert parallas to distance
tgas.add_virtual_columns_distance_from_parallax(tgas.parallax)
# 'overwrite' the real columns 'l' and 'b' with virtual columns
tgas.add_virtual_columns_eq2gal('ra', 'dec', 'l', 'b')
# and combined with the galactic sky coordinates gives galactic cartesian coordinates of the stars
tgas.add_virtual_columns_spherical_to_cartesian(tgas.l, tgas.b, tgas.distance, 'x', 'y', 'z')
j2000

Since RA. and Dec. are in degrees, while ra_error and dec_error is in miliarcseconds, so we put them on the same scale

In [43]:
tgas['ra_error'] = tgas.ra_error / 1000 / 3600
tgas['dec_error'] = tgas.dec_error / 1000 / 3600

We now let vaex sort out what the covariance matrix is for the cartesian coordinates x, y, and z. And take 50 samples from the datasets for visualization.

In [44]:
tgas.propagate_uncertainties([tgas.x, tgas.y, tgas.z])
tgas_50 = tgas.sample(50, random_state=42)

For this small dataset we visualize the uncertainties, with and without the covariance.

In [45]:
tgas_50.scatter(tgas_50.x, tgas_50.y, xerr=tgas_50.x_uncertainty, yerr=tgas_50.y_uncertainty)
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.show()
tgas_50.scatter(tgas_50.x, tgas_50.y, xerr=tgas_50.x_uncertainty, yerr=tgas_50.y_uncertainty, cov=tgas_50.y_x_covariance)
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.show()
_images/tutorial_84_0.png
_images/tutorial_84_1.png

From the second plot, we see that showing error ellipses (so narrow that they appear as lines) instead of error bars reveal that the distance information dominates the uncertainty in this case.

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.

In [46]:
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
Out[46]:
<vaex.promise.Promise at 0x7f214809d2e8>

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

In [47]:
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.

In [48]:
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.

In [49]:
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.

In [50]:
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)

In [50]:
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

Enabling notebook extension jupyter-js-widgets/extension...
      - Validating: OK
Enabling notebook extension bqplot/extension...
      - Validating: OK
Enabling notebook extension ipyvolume/extension...
      - Validating: OK
Enabling notebook extension jupyter-matplotlib/extension...
      - Validating: OK
Enabling notebook extension jupyter-leaflet/extension...
      - Validating: OK
In [51]:
# the default backend is bqplot, but we pass it here explicity
ds.plot_widget(ds.x, ds.y, f='log1p', 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:

In [52]:
a = np.array(['a', 'b', 'c'])
x = np.arange(1,4)
ds1 = vaex.from_arrays(a=a, x=x)
ds1
Out[52]:
#ax
0'a'1
1'b'2
2'c'3
In [53]:
b = np.array(['a', 'b', 'd'])
y = x**2
ds2 = vaex.from_arrays(b=b, y=y)
ds2
Out[53]:
#by
0'a'1
1'b'4
2'd'9

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.

In [54]:
ds1.join(ds2, left_on='a', right_on='b')
Out[54]:
#axby
0'a'1'a'1
1'b'2'b'4
2'c'3maskedmasked

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.

In [55]:
ds1.join(ds2, left_on='a', right_on='b', how='right')
Out[55]:
#byax
0'a'1'a'1
1'b'4'b'2
2'd'9maskedmasked

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.

In [56]:
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

In [58]:
nytaxi = vaex.open("/Users/maartenbreddels/datasets/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.

In [59]:
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.

In [60]:
%%time
nytaxi.mean(nytaxi.arc_distance)
CPU times: user 7.64 s, sys: 2.8 s, total: 10.4 s
Wall time: 2.63 s
Out[60]:
1.9999877196037037

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.

In [61]:
nytaxi['arc_distance_jit'] = nytaxi.arc_distance.jit_numba()
# nytaxi['arc_distance_jit'] = nytaxi.arc_distance.jit_pythran()
In [62]:
%%time
nytaxi.mean(nytaxi.arc_distance_jit)
CPU times: user 3.22 s, sys: 23.3 ms, total: 3.24 s
Wall time: 511 ms
Out[62]:
1.9999877196036921

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

In [ ]: