Astropy.Tables performance issue on array arithmetic

I was working on a project where I was taking multiple spectra and then dividing the counts per pixel by the exposure time to get a count rate. The spectra were stored in astropy.Tables, the exposure time in the first column followed by 2048 columns of counts per pixel. I wanted to perform a simple operation to divide each column (containing the counts in a given pixel over various exposure times) by the exposure column. This is very simple to do in a numpy array, something like:

array_data[:,1:] /= data[:,0, np.newaxis]

but the only way I could figure out to do the equivalent using astropy.Tables was to iterate column by column in a loop, which looked something like:

# I previously defined 'speccols', a list of columns with spectral data 
for cols in speccols:
      orig_Table[cols] = orig_Table[cols] / orig_Table['exposure']

It turns out to be a huge performance hit. Working with a structured array as produced by simply calling np.array(orig_Table) was not useful here, since I could not perform the simple operation above.

It was roughly 50 times faster to convert the floating point data in the Table orig_Table into a numpy array array_data using:

# list of columns with floating point data
colnames = orig_Table.colnames
floatcols = [cols for cols in colnames if orig_Table[cols].dtype == np.dtype('float64')]
# Copy just the floating point data to numpy array
fp_data = Table.copy(orig_data[floatcols])
fp_rows = len(fp_data)   # Number of rows of floating point data in table
fp_cols = len(floatcols)  # Number of columns of floating point data in table

# Create array of zeros same size as table and copy data over column by column
array_data = np.zeros((fp_rows,fp_cols),dtype=np.float64)
for n, name in enumerate(floatcols):
    array_data[:,n] = fp_data[name].astype(np.float64)

Is there supposed to be a faster way to do this sort of operation within the astropy.Tables context? Is there a simple way to convert a table into an numpy array of floats, since np.array() doesn't seem to cut it (presumably because to keep all the information in the astropy.Table intact requires structured arrays).

I'm not 100% on where the performance hit comes from, since in-place should be fast for Column too, but it may be that you are stuck with creating a large float array. Note that you could have a bit of "best of both worlds" if you make the table columns as views of the large array of spectra (use copy=False in the Column constructor).

Maybe using a line profiler can help track down the bottleneck here. You can use it as %lprun magic if you use Jupyter Notebook.

Doing the following is worst-case slow:

# I previously defined 'speccols', a list of columns with spectral data 
for cols in speccols:
      orig_Table[cols] = orig_Table[cols] / orig_Table['exposure']

The problem is that this syntax is actually creating a brand new column and updating the table to effectively remove the original column and add the newly created column. I haven't checked, but using orig_table[cols] /= orig_table['exposure'] should be much faster as it should be doing an in-place update.

Beyond that, is there a reason you are not using a vector column where the entire 2048-pixel spectrum is stored in a single column 'spectrum'? This is perfectly allowed and much preferred over 2048 individual columns (as implied in your second sentence).

This is a perfect candidate for adding documentation as described in #5753 😄

I get Tom's point about the creation of a new column in memory as opposed to in situ operations (you will notice when I performed the operation using numpy arrays, I did the operation in situ). I don't know if that was the main performance bottleneck. I may try playing around with it a bit.

To answer Tom's question about the "reason" for this organization, there were actually two instruments whose data I was recording and I was trying to keep the records together. So the first few columns were data from another instrument, a timestamp, and then the spectral data. I did originally have exposure time, then the spectra in a column, but as I started putting more pertinent data in there, it seemed to make sense to use the ability to refer to columns by name and besides, I couldn't store some of the information (e.g. timestamp, device flags) in the same column since they were of a different dtype.

This prompts some questions:

Is there a way to store the entire 2048 spectra vector in a single entry in one column with a 'vector' or array datatype? I am thinking in terms of storing the entire spectrum in a column, but with each entry in the column being a vector. A 2-D column of sorts. This allows the storage with other metadata (with different dtypes) in a row that associates it with the spectra.
Clearly astropy.Tables are very much Column-centric in their assumed data organization. Some of the Column operations are not supported for Rows (e.g. 'astype()'), which makes sense since each column is locked to one dtype. Does the documentation to explicitly mention this in terms of suggestions for organizing data? If not, it may be worth documenting this organizational approach as suggested by astrofrog. 👍

In-place operations should be faster on columns, but at the moment it seems like they aren't:

In [21]: t['a'] = np.random.random(1e8)

In [22]: %timeit t['a'] /= 1.3
1.2 s ± 185 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [23]: %timeit t['a'] = t['a'] / 1.3
1.14 s ± 78.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

This would indicate a bug in astropy.table?

And I can confirm that:

# I previously defined 'speccols', a list of columns with spectral data 
for col in speccols:
    orig_Table[col] /= orig_Table['exposure']

is much, much slower than the conversion to numpy approach. So the in situ advantage just isn't there for astropy.Tables.

@astrofrog

@astrofrog - the right test is using a small column so you are timing the table overhead rather than the time for doing the operation. BUT sadly that doesn't help!

However, here is a workaround which gives back a factor of 100 that suggests we need to implement the in-place operators, e.g. __iadd__ or whatnot. Timing of the in-place ops suggests they are still going through the heavy Column init machinery (that surprises me, but numpy subclassing is never easy).

In [21]: t = simple_table()

In [22]: %timeit t['b'] += 1e-6
1000 loops, best of 3: 319 µs per loop

In [23]: t
Out[23]: 
<Table length=3>
  a      b      c  
int64 float64  str1
----- -------- ----
    1 1.004111    c
    2 2.004111    d
    3 3.004111    e

In [24]: t = simple_table()

In [25]: %timeit x = t['b'].data; x += 1e-6
The slowest run took 13.41 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 2.45 µs per loop

In [26]: t
Out[26]: 
<Table length=3>
  a         b        c  
int64    float64    str1
----- ------------- ----
    1 1.41111099997    c
    2 2.41111100006    d
    3 3.41111100006    e

Note that t['b'].data += 0.1 fails because it thinks you are trying to set the data attribute which is not allowed. Fun!

About vector columns, no problem:

In [28]: specs = [[1,2], [3,4], [5,6]]

In [29]: exptimes = [5, 6, 7]

In [30]: t = Table([exptimes, specs], names=['exptime', 'spectrum'])

In [31]: t
Out[31]: 
<Table length=3>
exptime spectrum [2]
 int64     int64    
------- ------------
      5       1 .. 2
      6       3 .. 4
      7       5 .. 6

In [32]: t['spectrum'][1]
Out[32]: array([3, 4])

About docs, 👍 on the performance docs section. I also think it would be worth funding somebody to work on table performance, but that is a separate discussion. 😄

There is the following that makes it pretty clear how tables are internally implemented:
http://astropy.readthedocs.io/en/latest/table/implementation_details.html
Perhaps a sentence up front talking about column-centric implementation and a link here would be worthwhile.

@taldcroft

@taldcroft:

@astrofrog - the right test is using a small column so you are timing the table overhead rather than the time for doing the operation. BUT sadly that doesn't help!

The same operation with Numpy arrays takes only 200ms so I am assuming that the overhead (5x) is from the memory copy not from the actual operation.

The problem does not seem to be with Column:

a = np.random.random(1000000)
%timeit b = np.array(a, copy=False); b/=2.
# 100 loops, best of 3: 4.03 ms per loop
%timeit c = Column(a, name='a', copy=False); c/=2.
# 100 loops, best of 3: 4.23 ms per loop

However, as part of a table it is not so good:

t = QTable()
t['a'] = c
%timeit t['a'] /= 2.
# 100 loops, best of 3: 7.35 ms per loop

What is worse, though, is that there must be some memory leak: after doing this, python uses about 1GB...

This does not happen if I do:

%timeit c = t['a']; c /= 2.
# 100 loops, best of 3: 4.05 ms per loop

(or c = t.columns[0])

@JuanCab

@JuanCab - from my tests, I think the work-around is clear: you'd need to do something like:

for cols in speccols:
    col = orig_Table[cols]
    col /= orig_Table['exposure']

Of course, ideally, this should not be necessary... (and I don't really understand why it is happening)

Note: this also happens with quantity, so it is not Column specific. Indeed,

for i in range(50): t['a'] = a

blows up memory usage, so the problem seems to be in Table.__setitem__.

Following the trail: t.replace_column also blows up memory... I'm actually a bit confused with what is done there. Is it meant to copy the data in a?

Yes, replace_column copies the data in a. This comes from (and is consistent with) the default of copy=True for table initialization.

About the memory leak...

    def replace_column(self, name, col):
        t = self.__class__([col], names=[name])
        cols = OrderedDict(self.columns)
        cols[name] = t[name]
        self._init_from_cols(cols.values())

So this is making a temporary table in order to use the generic Table init machinery to make sense of whatever the user has provided for col (which might be a column or an array or list ...). Then it makes a reference to that column and uses it in self._init_from_cols, which does NOT copy. So this should all be fine? Maybe _init_from_cols on an existing table somehow leaves around a ref to the original cols? But I don't quite see how...

docs: Reference issue astropy#6253

Navigation Menu

Repository navigation

Astropy.Tables performance issue on array arithmetic #6253

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Astropy.Tables performance issue on array arithmetic #6253

Description

Activity

mhvk commented on Jun 22, 2017

pllim commented on Jun 22, 2017

taldcroft commented on Jun 22, 2017

astrofrog commented on Jun 22, 2017

JuanCab commented on Jun 22, 2017

astrofrog commented on Jun 22, 2017

JuanCab commented on Jun 22, 2017

taldcroft commented on Jun 22, 2017

taldcroft commented on Jun 22, 2017

taldcroft commented on Jun 22, 2017

astrofrog commented on Jun 23, 2017

mhvk commented on Jun 23, 2017

mhvk commented on Jun 23, 2017

mhvk commented on Jun 23, 2017

mhvk commented on Jun 23, 2017

taldcroft commented on Jun 23, 2017

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions