Subtracting two interleaved, differently based time-series arrays with numpy?
Question: I have two datasets, `a[ts1]` and `b[ts2]`, where `ts1` and `ts2` are
timestamps taken at different times (in different bases?). I wanted to plot
`b[ts2]-a[ts1]`, but I think I made a mistake, in that the plotting software
understood that I want `b[i]-a[i]` instead, where `i` is the index order of
the value.

So I wanted to make a small example of this with `numpy`, and I realized I
have no idea whether, and how, `numpy` can perform this operation - but using
vectors, and avoiding `for` loops. I have made an example (below), that
defines `a[ts1]` and `b[ts2]` as `numpy` structured arrays titled `a_np` and
`b_np`:

    
    
    array([(0.0, 0.0), (0.8865606188774109, 0.30000001192092896),
           (1.6939274072647095, 0.6000000238418579),
           (2.3499808311462402, 0.8999999761581421)], ...
          dtype=[('a', '<f4'), ('ts1', '<f4')])
    array([(0.3973386585712433, 0.10000000149011612),
           (0.7788366675376892, 0.20000000298023224),
           (1.4347121715545654, 0.4000000059604645), (1.6829419136047363, 0.5)], ...
          dtype=[('b', '<f4'), ('ts2', '<f4')])
    

So my questions here are:

  * What is this class of arrays/problems called? Is it just "time-series" arrays? Basically, they describe a 1-D signal, but since timestamps have to be kept, it's a 2D array; and since the "time" column can have any meaning, I guess it can be generalized to interpolation of values of arrays in (value) column, over the "indexing" values in the (time/index) column.
  * Can `numpy` do a vectorized operation where the arrays are properly interpolated in "time", before being subtracted?

Looking for info on this, I found [pandas: Python Data Analysis
Library](http://pandas.pydata.org/); I guess I should use this instead, given
it has a "time-series" functionality - but in this case, I don't need any
fancy interpolation of sample values - just a "step" or "hold" one (basically,
no interpolation); which is why I was wandering if `numpy` can do that in a
vectorized fashion. Otherwise, the example below employs `for` loops.

* * *

The example will result with an image like this:

![array-timeseries](http://i.stack.imgur.com/tmH1O.png)

Arrays `a` and `b` represent values taken at different times, which is
indicated by their respective `impulses`; `a` is plotted with `lines` (so,
linearly interpolated before plotting), but `b` with `steps` (to indicate the
actual values that exist)

The array `d1` represents the "original" difference `b[t]-a[t]` taken when
constructing the array - obviously, I do not have access to this data in
reality, so I have to work from the sampled values. In that case, the
difference `b[ts2]-a[ts1]` is shown as the array/signal `d2`, again as `steps`
to emphasize the errors being made in respect to the "original". This `d2` is
what I'd like to calculate with `numpy` (but below, it is calculated in the
same `for` loop).

The error I made with my plotting software is getting an index-by-index
difference of `b` and `a`, or `b[i]-a[i]`; this is shown as array/signal `e`
\- and as shown, it is _way_ off from what it's otherwise supposed to
represent. This is only the case if sampling intervals in the two signals are
uneven; try `modulowith = 2` in the code, then `e` is actually not that bad -
however, my real life case has uneven timestamps, so `b[i]-a[i]` doesn't help
me at all.

Here is the code, which also calls `gnuplot` (tested on Python 2.7, `numpy`
1.5 I think):

    
    
    import subprocess
    import math, random
    import numpy as np
    from pprint import pprint
    from numpy.lib.recfunctions import append_fields
    
    step = 0.1
    modulowith = 3
    # must init all arrays separately;
    # a=b=[] makes a==b by reference!
    ts1 = []; ts2 = [] ; tsd = []
    valsa = []; valsb = []; valsd1 = []; valsd2 = []
    stra = strb = strd1 = strd2 = "" ; kval1 = kval2 = 0
    for ix in range(0, 100, 1):
      ts = ix*step
      val1 = 3.0*math.sin(ts) #+random.random()
      val2 = 2.0*math.sin(2.0*ts)
      if ( ix%modulowith == 0):
        ts1.append(ts) ; valsa.append(val1)
        stra += "%.03f %.06f\n" % (ts, val1)
        kval1 = val1
      else:
        ts2.append(ts) ; valsb.append(val2)
        strb += "%.03f %.06f\n" % (ts, val2)
        kval2 = val2
      tsd.append(ts)
      valb = val2 - val1 ; valsd1.append(valb)
      strd1 += "%.03f %.06f\n" % (ts, valb)
      valc = kval2 - kval1 ; valsd2.append(valc)
      strd2 += "%.03f %.06f\n" % (ts, valc)
    
    a_np = np.array(
      [(_valsa,) for _valsa in valsa],
      dtype=[('a','f4')]
    )
    b_np = np.array(
      [(_valsb,) for _valsb in valsb],
      dtype=[('b','f4')]
    )
    a_np = append_fields(a_np, names='ts1', data=ts1, dtypes='f4', usemask=False)
    b_np = append_fields(b_np, names='ts2', data=ts2, dtypes='f4', usemask=False)
    
    pprint(a_np[:4])
    pprint(b_np[:4])
    
    # e_np = np.subtract(b_np['b'],a_np['a'])
    # (via field reference) is same as doing:
    # e_np = np.subtract(np.array(valsa, dtype="f4"), np.array(valsb, dtype="f4"))
    # but for different sized arrays, must do:
    e_np = b_np['b'] - np.resize(a_np, b_np.shape)['a']
    pprint(e_np[:4])
    e_str = ""
    for ts, ie in zip(ts2, e_np):
      e_str += "%.03f %.06f\n" % (ts, ie)
    
    
    gpscript = """
    plot "-" using 1:2 with lines lc rgb "green" t"a", \\
    "" using 1:2 with impulses lc rgb "green" t"", \\
    "" using 1:2 with steps lc rgb "blue" t"b", \\
    "" using 1:2 with impulses lc rgb "blue" t"", \\
    "" using 1:2 with lines lc rgb "red" t"d1", \\
    "" using 1:2 with steps lc rgb "orange" t"d2", \\
    "" using 1:2 with steps lc rgb "brown" t"e"
    -
    {0}
    e
    {0}
    e
    {1}
    e
    {1}
    e
    {2}
    e
    {3}
    e
    {4}
    e
    """.format(stra, strb, strd1, strd2, e_str)
    
    proc = subprocess.Popen(
      ['gnuplot','--persist'],
      shell=False,
      stdin=subprocess.PIPE,
    )
    proc.communicate(gpscript)
    

* * *

Thanks to the answer by
[@runnerup](http://stackoverflow.com/a/20449087/277826), here is a slightly
verbose (for syntax example purposes) `numpy`-only solution:

    
    
    # create union of both timestamp arrays as tsz
    ntz = np.union1d(b_np['ts2'], a_np['ts1'])
    # interpolate `a` values over tsz
    a_z = np.interp(ntz, a_np['ts1'], a_np['a'])
    # interpolate `b` values over tsz
    b_z = np.interp(ntz, b_np['ts2'], b_np['b'])
    # create structured arrays for resampled `a` and `b`,
    # indexed against tsz timestamps
    a_npz = np.array( [ (tz,az) for tz,az in zip(ntz,a_z) ],
      dtype=[('tsz', 'f4'), ('a', 'f4')] )
    b_npz = np.array( [ (tz,bz) for tz,bz in zip(ntz,b_z) ],
      dtype=[('tsz', 'f4'), ('b', 'f4')] )
    # subtract resized array
    e_npz = np.subtract(b_npz['b'], a_npz['a'])
    e_str = ""
    
    # check:
    pprint(e_npz[:4])
    # gnuplot string:
    for ts, ie in zip(ntz, e_npz):
      e_str += "%.03f %.06f\n" % (ts, ie)
    

This is linearly interpolated, so it will be different from `d2` above, but
still nicely fitting.

If there wasn't for those `for` loops creating the arrays, then it is
vectorized - and in principle I don't even have to create those arrays - just
wanted to see how they will look like as structured ones.. In all, I guess I
was hoping for a one-liner that would do this, with structured arrays (that
is, handling field names as well).


Answer: this is an attempt to sell you on switching to `pandas` : )

    
    
    import numpy as np
    import pandas as pd
    import datetime as dt
    import matplotlib.pyplot as plt
    
    # one minute interval
    start = dt.datetime.now( )
    end = start + dt.timedelta( minutes=1 )
    
    # sin curve at seconds frequancy
    idx1 = pd.date_range( start, end, freq='S' )
    ts1 = pd.Series( np.sin( np.linspace( 0, 4 * np.pi, len( idx1 ) ) ), index=idx1 )
    
    # cosine curve at milisecond frequency
    idx2 = pd.date_range( start, end, freq='L' )
    ts2 = pd.Series( np.cos( np.linspace( 0, 4 * np.pi, len( idx2 ) ) ), index=idx2 )
    

now `len( ts1 ) = 61` and `len( ts2 ) = 6001`, with different frequencies

    
    
    fig = plt.figure( figsize=(8, 6) )
    ax = fig.add_axes( [.05, .05, .9, .9] )
    
    ts1.plot( ax, color='DarkBlue' )
    ts2.plot( ax, color='DarkRed' )
    
    # reindex ts2 like ts1
    ts2 = ts2.reindex_like( ts1 )
    (ts1 - ts2).plot( ax, color='DarkGreen' )
    

and you get:

![time series](http://i.stack.imgur.com/T184l.png)

**edit** : for interpolation purposes you may use nonparametric methods in
statsmodels, so basically you may interpolate one series at the frequencies of
the other one and then subtract the two:

    
    
    import statsmodels.api as sm
    n = 1000
    x = np.linspace( 0, 1, n )
    y = np.random.randn( n ).cumsum( )
    z = sm.nonparametric.lowess( y, x, return_sorted=False, frac=.05)
    
    ax.plot( x, y, 'Blue', linestyle='--' )
    ax.plot( x, z, color='DarkRed' )
    

![bm](http://i.stack.imgur.com/DFyox.png)

