Debugging time decoding performance problems

[rabernat]

Over in #1385 (comment), @jtomfarrar posted an example of a NetCDF file which was really slow to open unless one specifies decode_times=False. I spent some time debugging this, and I thought I would post my findings here.

The file in question is here: https://drive.google.com/file/d/1-05bG2kF8wbvldYtDpZ3LYLyqXnvZyw1/view?usp=sharing

When we open it with decode_times=False it is very fast.

import xarray as xr
fname = '/Users/rpa/tmp/SMODE_PFC_Wavegliders_WHOI43.nc'
%time ds = xr.open_dataset(fname, decode_cf=False)
# CPU times: user 49 ms, sys: 78 ms, total: 127 ms
# Wall time: 127 ms

If we don't specify that option, it takes forever and also produces the following warning

/opt/miniconda3/envs/pangeo-forge-recipes/lib/python3.9/site-packages/xarray/coding/times.py:527: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range

This was a 🚩 for me. Why would recent data from an instrument need anything but a totally standard datetime. So I looked at time variable values and metadata.

print(ds.time_20Hz.values[:4])
# array([738457.04166667, 738457.04166725, 738457.04166782, 738457.0416684 ])
print(ds.time_20Hz.units)
# 'days since 1950-01-01T00:00:00Z'

This did not make sense. If you believe these units, that would place us in the year 3971. This is out of the range of what pandas time indexes can hold, which was triggering the usage of cftime.

So next I tried removing an offset from the time coordinate

with xr.set_options(keep_attrs=True):
    dsa = ds.assign_coords(
        time_20Hz=(ds.time_20Hz - 738457),
        time_1Hz=(ds.time_1Hz - 738457)
    )

and then decoding again

%time dsb = xr.decode_cf(dsa)
print(dsb)

CPU times: user 185 ms, sys: 55.1 ms, total: 240 ms
Wall time: 239 ms
<xarray.Dataset>
Dimensions:                      (time_1Hz: 771961, time_20Hz: 15439201)
Coordinates:
  * time_1Hz                     (time_1Hz) float64 0.04167 0.04168 ... 8.976
    latitude_1Hz                 (time_1Hz) float64 ...
    longitude_1Hz                (time_1Hz) float64 ...
  * time_20Hz                    (time_20Hz) datetime64[ns] 1950-01-01T00:59:...
    latitude_20Hz                (time_20Hz) float64 ...
    longitude_20Hz               (time_20Hz) float64 ...
    altitude_20Hz                (time_20Hz) float32 ...

Things are looking a lot better. We now have a datetime index for time_20Hz. However, I noticed that time_1Hz was not getting decoded at all. Looking deeper at the metadata

print(dsb.time_1Hz.attrs)

{'long_name': 'Time recorded by Campbell CR6 datalogger at 1Hz',
 'standard_name': 'time',
 'unit': 'days since 1950-01-01T00:00:00Z',
 'axis': 'T',
 'coverage_content_type': 'coordinate'}

I noticed that the wrong attribute was used for the units (unit instead of units). Fixing this allows both variables to be decoded quickly and correctly.

dsa.time_1Hz.attrs['units'] = dsa.time_1Hz.attrs['unit']
del dsa.time_1Hz.attrs['unit']
dsb = xr.decode_cf(dsa)

I think it is useful to step through this because it shows how a metadata problem can translate to a performance problem.

[kmuehlbauer]

@rabernat Great sleuthing. 👍 💯 This is another example of how much needed metadata sanity checks are.

[dcherian]

These are matlab datetimes... (hi tom!)

Ryan, does it look like creating the cftimeindex is really slow?

[rabernat]

Creating the cftimeindex is insanely slow. I did not have the patience to even let it finish and interrupted after a few minutes.

[jtomfarrar]

My colleague @bpgreenwood is fixing the time variable. We were wondering, is there an optimal or 'native' time format that would give the best performance, or is it just that we need to not go outside the Timestamp-valid range?

[dcherian]

It's incredibly annoying because matlab references to year 0 which doesn't exist in the default calendar

dnum = np.array([738457.04166667, 738457.04166725, 738457.04166782, 738457.0416684])

da = xr.DataArray(
    dnum-366, name="time", dims="time", attrs={"units": "days since 0001-01-01 00:00"}
)
xr.decode_cf(da.to_dataset())

This decodes to

array(['2021-10-27T01:00:00.000288000', '2021-10-27T01:00:00.050399000',
       '2021-10-27T01:00:00.099644000', '2021-10-27T01:00:00.149764000'],
      dtype='datetime64[ns]')

If that looks right then you could

1. subtracting 366 from your matlab datenum
2. setting the units to "days since 0001-01-01 00:00"

Though with 20Hz data, you might want to convert the units to "seconds since XXXX" to preserve precision.

[rabernat]

is there an optimal or 'native' time format that would give the best performance

The optimal format would probably be to encode the times as an integer, rather than a float. So for 20Hz data you might use ms, [0, 50, 100, 150, ...] together with units of "ms since 0001-01-01 00:00". If there is jitter in your clock and you need a floating point type, you will get the best precision by having the numbers be O(1).

[spencerkclark]

+1 for integers where possible. I would also recommend using a reference date that can be represented using a np.datetime64[ns] value ("0001-01-01 00:00:00" cannot). This way cftime will not need to be used at all in the decoding process.

I may dig into the performance issues with cftime at some point later.

[bpgreenwood]

Great suggestions, thanks for the tips everyone.

[spencerkclark]

For posterity, here are some notes regarding the performance when using cftime (though I think we've established that ideally cftime would not be used at all for this example).

So the cftime performance is not xarray-related. It has to do with two things:

1. The array of times is very long (15+ million elements) -- this is much longer than arrays we typically encounter (and so is not necessarily something that was prioritized in recent design decisions in cftime).
2. The times are encoded using large values (thousands of years from their reference date).

A while back I made a change to cftime.num2date such that it would use exact timedelta arithmetic to decode times rather than an approximate algorithm that used floating point numbers (the inexactness of time decoding of cftime dates was a long-standing issue). One of the performance implications of this I noted at the time was that times encoded with a distant reference date took longer to decode than times encoded with a more proximal reference date. The new approach was chosen despite this slowdown on the basis that exact was better than approximate, and that it was actually faster than the old algorithm for typical ranges of times (see more discussion in this PR: Unidata/cftime#176).

The upshot is that it further emphasizes the benefits of choosing a proximal reference date, at least when using cftime. In this particular example, choosing a proximal reference date affords a nearly 300x speedup:

In [6]: times = ds.time_20Hz.values[:10000]

In [7]: proximal_times = times - times[0]

In [8]: %timeit cftime.num2date(times, units="days since 0001-01-01", calendar="proleptic_gregorian")
3.16 s ± 31.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [9]: %timeit cftime.num2date(proximal_times, units="days since 2021-10-27T01:00:00.000288", calendar="proleptic_gregorian")
10.6 ms ± 96.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Here I chose to limit the size of the problem, because decoding 15 million distant times with this approach seems like it would take over an hour. In contrast, with a proximal reference date it might take on the order of fifteen seconds (still long, but bearable).

[spencerkclark]

See Unidata/cftime#269 for a thorough analysis of the possible optimization approaches -- computing each time relative to the previous seems like the best approach overall.

[spencerkclark]

With Unidata/cftime#270 merged, the performance of cftime.num2date in these circumstances should be much improved in the next cftime release!

[jtomfarrar]

Thank you, Spencer! I appreciate the effort you and others are putting into developing these tools.

[rabernat]

Just noting what success story this is! 🎉 In about one week, we went from a user-reported problem to an order-of-magnitude speedup in the upstream library.

Also worth observing that xarray was (as usual) the reporting point for a performance issue that was not in fact resolved in xarray.

[dcherian]

Time for a blog post?

[jtomfarrar]

I looked a little more on size of int. Numpy's docs explain:

The behaviour of NumPy and Python integer types differs significantly for integer overflows and may confuse users expecting NumPy integers to behave similar to Python’s int. Unlike NumPy, the size of Python’s int is flexible. This means Python integers may expand to accommodate any integer and will not overflow.

Indeed, I think this is what was confusing me. (I'm still confused-- I guess integers in a netcdf file would be interpreted by xarray as Numpy int's? Maybe I should go read the xarray docs.)

https://numpy.org/doc/stable/user/basics.types.html

[rabernat]

Yes, all data in xarray (and netCDF) in python is represented using numpy arrays.

Now it's my turn to learn something! I definitely did not know that about python's builtin ints.

[dopplershift]

Also important to keep in mind that by the time you're talking about what's being written to a file, you're talking about file format specifications dictated by netCDF4 itself, and really HDF5. At that point Python/Numpy are irrelevant, and it's the raw primitive C datatypes.

[bpgreenwood]

Thank you all for your help. I'll share a Matlab function that I wrote that converts Matlab datenum array to 'ms since [user defined epoch]'. Tom, our dataset contains ~9 days of 20Hz data. We can either use:

• uint32 datatype "ms since 2021-10-01T00:00:00.000Z" (non-standard epoch but less space)
• uint64 datatype "ms since 2000-01-01T00:00:00.000Z" (more standard epoch but 2x the space)

% convert Matlab datenum timestamp array to 'ms since epoch' for use with NetCDF
% refer to https://github.com/pydata/xarray/discussions/6284#discussioncomment-2205538
% epoch    - datenum scalar representing start of time period
% time     - datenum array of timestamps
% int_type - 'uint32' or 'uint64' string defines time datatype that is returned
% returns uint32 or uint64 array of timestamps in 'ms since epoch' format
function [time_ms] = convert_time_ms(epoch,time,int_type)

    fprintf(' * converting timestamps from Matlab datenum to NetCDF "ms since %s"\n',datestr(epoch,'YYYY-mm-ddTHH:MM:SS.fffZ'));
    
    if any(time < epoch)
        fprintf(' * warning: timestamps found prior to epoch; see convert_time_ms()\n');
    end
    
    % unit conversion using user specified integer data type
    if strcmp(int_type,'uint32')
        time_ms = uint32((time - epoch)*24*60*60*1000);
    elseif strcmp(int_type,'uint64')
        time_ms = uint64((time-epoch)*24*60*60*1000);
    else
        % handle unknown data type exception
        fprintf(' * error: convert_time_ms() unknown data type\n');
        keyboard;
    end
    
    % detect overflow
    if (time_ms(end) == intmax(int_type))
        fprintf(' * warning: timestamps exceed %s max range: [0,%u]\n',int_type,intmax(int_type));
    end
end

[dopplershift]

Just note that, in terms of file size compression can make that much less than 2x due to so many of those extra bytes being 0:

import os

from netCDF4 import Dataset
import numpy as np

num = 100_000_000

for fname, dtype in [('small.nc', np.uint32), ('big.nc', np.uint64)]:
    nc = Dataset(fname, 'w')
    nc.createDimension('time', num)
    time = small.createVariable('time', dtype, ('time',), zlib=True, complevel=4)
    time[:] = (np.arange(num) * 5).astype(time.dtype)
    nc.close()
    with open(fname, 'rb') as readfile:
        print(f'{fname}: {len(readfile.read())} bytes')
        os.remove(fname)

gives for me:

small.nc: 5793654 bytes
big.nc: 6310239 bytes