Script to determine the drift of the PyXID box

gistfile1.py

# -*- coding: utf-8 -*-
"""Script to determine the drift of the PyXID box"""

import numpy as np
import pyxid
import time
import matplotlib.pyplot as plt

n_run = 50
dur = 60.  # seconds
dt = 0.01   # time between ticks

devices = pyxid.get_xid_devices()
dev = devices[0]  # get the first device to use


def get_fits(ts, vals):
    # Do a fit
    b, m = np.dot(np.linalg.pinv(np.array((np.ones_like(vals), vals)).T), ts)
    # Check for outliers beyond the expected jitter range and re-fit
    diffs = (m * vals + b) - ts
    goods = np.abs(diffs) < 0.00075
    if goods.mean() < 0.99:
        raise RuntimeError('bad run')
    if (~goods).sum() > 0:
        print('  Discarding %d/%d points in refit'
              % ((~goods).sum(), len(goods)))
    # Refit data
    b, m = np.dot(np.linalg.pinv(np.array((np.ones_like(vals[goods]),
                                           vals[goods])).T), ts[goods])
    return b, m


all_ts, all_vals = [], []
all_slopes = np.zeros(n_run)
print('running')
for ii in range(n_run):
    ts, vals = [], []
    t0 = time.time()
    while(time.time() - t0 < dur):
        ts.append(time.time())
        vals.append(dev.query_base_timer())
        time.sleep(dt)
    all_ts.extend(ts)
    all_vals.extend(vals)

    ts = np.array(ts)
    ts -= ts[0]
    vals = np.array(vals) / 1e3

    b, m = get_fits(ts, vals)
    plt.plot(ts, np.zeros_like(ts), 'k')
    plt.plot(ts, vals + b - ts, 'b')
    plt.plot(ts, m * vals + b - ts, 'r')
    plt.ylabel('relative diff (ms)')
    plt.xlabel('Time (s)')
    plt.axis('tight')
    print('  slope[%d] = %s' % (ii, m))
    all_slopes[ii] = m

print('slope1:  %s sec/sec (std: %s μs/sec)'
      % (np.median(all_slopes), 1e6 * np.std(all_slopes)))

all_ts = np.array(all_ts)
all_ts -= all_ts[0]
all_vals = np.array(all_vals)
b, m = get_fits(all_ts, all_vals)
print('slope2: %s sec/sec' % (1e3 * m))

plt.show()

# Empirical results on Ubuntu 14.10 run as::
#
#     $ sudo nice -n -20 su -c "python -O pyxid_problem.py" $USER
#
# Note that these slopes are in addition to the
# expected 1e-3 ms->sec scale factor.
#
# n_run  dur   slope
#    20    1   1.00066084174 sec/sec (std: 128.293726741 μs/sec)
#    20    5   1.00065111157 sec/sec (std: 10.3710617819 μs/sec)
#    10   10   1.00064919427 sec/sec (std: 1.81778368813 μs/sec)
#    20   20   1.00064933238 sec/sec (std: 1.05694509903 μs/sec)
#    50   60   1.00064206973 sec/sec (std: 3.09985183877 μs/sec)