## scipy.stats.chi2.ppf() without scipy

Allantools has a dependence on scipy because it uses scipy.stats.chi2.ppf() (inverse of the chi-squared cumulative distribution function) in the code for confidence intervals.

When testing scipy takes quite a long time to install and the whole package seems a bit overkill for just this one function.

So I tried implementing it using just numpy. It sort of works, but not for corner cases where p is very close to 1.0 or k is a big number. I get for example:

```    top = pow(x, s)*math.exp(-x)*pow(x,k)
OverflowError: (34, 'Numerical result out of range')
```

Maybe there's something trick when numerically evaluating lower_gamma() using the series expansion that I am missing...

```import scipy.stats import math import numpy   def bisection(function, k, p, a, b, tol): # http://code.activestate.com/recipes/578417-bisection-method-in-python/ assert (function(a,k)-p)*(function(b,k)-p) < 0 # a, b must bracket root c = (a+b)/2.0 while (b-a)/2.0 > tol: if (function(c, k)-p) == 0: return c elif (function(a,k)-p)*(function(c,k)-p) < 0: b = c else : a = c c = (a+b)/2.0 return c   def lower_gamma(s,x): # lower incomplete gamma function # https://en.wikipedia.org/wiki/Incomplete_gamma_function#Evaluation_formulae g = 0 last_g = 1.0 done = False tol = 1.0e-6 k=0 while not done: top = pow(x, s)*math.exp(-x)*pow(x,k) bot = numpy.prod( [float(s+j) for j in range(k+1) ] ) dg = float(top)/float(bot) if dg == float("Inf"): break g += dg k += 1 if k>100: # get at least 100 terms in the sum if g==0: break delta = abs(dg/g) if delta == float("Inf"): break if delta < tol: done = True last_g = g return g   def chi2_cdf(x, k): # chi-squared cumulative density function # cdf(x; k) = lower_gamma(k/2, x/2) / gamma(k/2) return lower_gamma(k/2.0, x/2.0) / math.gamma(k/2.0)   def chi2_ppf(p, k): # chi-squared Percent point function (inverse of cdf percentiles). # look for x such that # p = chi2_cdf( x=chi2_ppf(p, k), k) tol = 1e-8 lolim = 0 hilim = k while (chi2_cdf(lolim,k)-p)*(chi2_cdf(hilim,k)-p) > 0: hilim *= 1.5 return bisection( chi2_cdf, k, p, lolim, hilim, tol)   print "scipy cdf: ",scipy.stats.chi2.cdf(55, 33) print "own cdf: ",chi2_cdf(55, 33)     print "scipy ppf ", scipy.stats.chi2.ppf(0.4, 33) print " own ppf ", chi2_ppf(0.4, 33)   # test that we really found the inverse print scipy.stats.chi2.cdf(scipy.stats.chi2.ppf(0.4, 33), 33) print chi2_cdf( chi2_ppf(0.4, 33), 33 )   # try to check the scipy function against our own function # for some random input of (p, k) for n in range(100): k = numpy.random.randint(20, 200) p = numpy.random.random() print k, p, a=scipy.stats.chi2.ppf(p, k) b=chi2_ppf(p, k) ok = numpy.isclose(a, b) if ok: print ok else: print ok, a, b assert ok```

## Leap second time-stamps for leap-seconds.list

This function generates time-stamps in the format used in leap-seconds.list (on e.g. NTP or PTP servers). The format is integer seconds since 1900-01-01T00:00+00.

```import datetime import pytz     def generate_ntp_timestamp(year,month,day,hour,minute,second): t = datetime.datetime(year,month,day,hour,minute,second,tzinfo=pytz.utc) #NTP leap-seconds list wants seconds since 1900 epoch_start = datetime.datetime(1900,1,1,0,0,0,tzinfo=pytz.utc) delta = t-epoch_start delta_s = delta.total_seconds() return int(delta_s)     # test that this generates OK values. # from existing leap-seconds.list # 2918937600 27 # 1 Jul 1992 # 2950473600 28 # 1 Jul 1993 # 2982009600 29 # 1 Jul 1994 # 3029443200 30 # 1 Jan 1996 # 3076704000 31 # 1 Jul 1997 # 3124137600 32 # 1 Jan 1999 # 3345062400 33 # 1 Jan 2006 # 3439756800 34 # 1 Jan 2009 print "1 Jan 2009: ", generate_ntp_timestamp(2009,1,1,0,0,0) print "1 Jan 2006: ", generate_ntp_timestamp(2006,1,1,0,0,0) print "1 Jul 1997: ", generate_ntp_timestamp(1997,7,1,0,0,0) #output: #1 Jan 2009: 3439756800 #1 Jan 2006: 3345062400 #1 Jul 1997: 3076704000```

## Clock display with Python and Tkinter

On github: https://github.com/aewallin/digiclock

A simple clock display with local time, UTC, date (iso8601 format of course!), MJD, day-of-year, and week number. Useful as a permanent info-display in the lab - just make sure your machines are synced with NTP.

```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 # use Tkinter to show a digital clock # tested with Python24 vegaseat 10sep2006 # https://www.daniweb.com/software-development/python/code/216785/tkinter-digital-clock-python # http://en.sharejs.com/python/19617   # AW2015-04-24 # added: UTC, localtime, date, MJD, DOY, week from Tkinter import * import time import jdutil # https://gist.github.com/jiffyclub/1294443   root = Tk() root.attributes("-fullscreen", True) # this should make Esc exit fullscrreen, but doesn't seem to work.. #root.bind('',root.attributes("-fullscreen", False)) root.configure(background='black')   #root.geometry("1280x1024") # set explicitly window size time1 = '' clock_lt = Label(root, font=('arial', 230, 'bold'), fg='red',bg='black') clock_lt.pack()   date_iso = Label(root, font=('arial', 75, 'bold'), fg='red',bg='black') date_iso.pack()   date_etc = Label(root, font=('arial', 40, 'bold'), fg='red',bg='black') date_etc.pack()   clock_utc = Label(root, font=('arial', 230, 'bold'),fg='red', bg='black') clock_utc.pack()   def tick(): global time1 time2 = time.strftime('%H:%M:%S') # local time_utc = time.strftime('%H:%M:%S', time.gmtime()) # utc # MJD date_iso_txt = time.strftime('%Y-%m-%d') + " %.5f" % jdutil.mjd_now() # day, DOY, week date_etc_txt = "%s DOY: %s Week: %s" % (time.strftime('%A'), time.strftime('%j'), time.strftime('%W'))   if time2 != time1: # if time string has changed, update it time1 = time2 clock_lt.config(text=time2) clock_utc.config(text=time_utc) date_iso.config(text=date_iso_txt) date_etc.config(text=date_etc_txt) # calls itself every 200 milliseconds # to update the time display as needed # could use &gt;200 ms, but display gets jerky clock_lt.after(20, tick)   tick() root.mainloop()```

This uses two small additions to jdutil:

```1 2 3 4 5 6 7 def mjd_now(): t = dt.datetime.utcnow() return dt_to_mjd(t)   def dt_to_mjd(dt): jd = datetime_to_jd(dt) return jd_to_mjd(jd)```

## Faster allantools with numpy

Update: these slopes are both linear, or O(n), where n is the input data size.

However for the worst performing algorithm, MTIE, an O(n*log(n)) algorithm that uses a binary tree to store intermediate results is possible. This algorithm restricts the tau-values to 2**k (with k integer) times the data interval.

Danny Price has made a fantastic contribution to allantools  by taking my previous pure python code (hacked together in January 2014, shown in red) and numpyifying it all (shown in blue)!

For large datasets the speedups are 100-fold or more in most cases!
The speedups are calculated only for datasets larger than 1e5, since python's time.time() doesn't seem suitable for measuring 1 ms or shorter running times.

## ZeromMQ REQuest/REPly pattern

This pattern is for a reply&request situation where one machine asks for data from another machine. The pattern enforces an exchange of messages like this:

REQuester always first sends, then receives.

Requester

```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 import zmq import time   # request data from this server:port server = "192.168.3.38" port = 6000   context = zmq.Context() socket = context.socket(zmq.REQ) socket.connect ("tcp://%s:%s" % (server,port) )   while True: tstamp = time.time() msg = "%f" % tstamp # just send a timestamp # REQuester: TX first, then RX socket.send(msg) # TX first rep = socket.recv() # RX second done = time.time() oneway = float(rep)-tstamp twoway = done - tstamp print "REQuest ",msg print "REPly %s one-way=%.0f us two-way= %.0f us" % ( rep, 1e6*oneway,1e6*twoway ) time.sleep(1)   # typical output: # REQuest 1392043172.864768 # REPly 1392043172.865097 one-way=329 us two-way= 1053 us # REQuest 1392043173.866911 # REPly 1392043173.867331 one-way=420 us two-way= 1179 us # REQuest 1392043174.869177 # REPly 1392043174.869572 one-way=395 us two-way= 1144 us```

The one-way delay here assumes we have machines with synchronized clocks. This might be true to within 1ms or so if both machines are close to the same NTP server.

The two-way delay should be more accurate, since we subtract two time-stamps that both come from the REQuesters system clock. Much better performance can probably be achieved by switching to C or C++.

```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 import zmq import time   port = 6000 context = zmq.Context() socket = context.socket(zmq.REP) # reply to anyone who asks on this port socket.bind("tcp://*:%s" % port)   print "waiting.." while True: # Wait for next request from client message = socket.recv() # RX first print "Received REQuest: ", message tstamp = time.time() sendmsg = "%f" % tstamp # send back a time-stamp print "REPlying: ",sendmsg socket.send(sendmsg) # TX second   # typical output: # Received REQuest: 1392043244.206799 # REPlying: 1392043244.207214 # Received REQuest: 1392043245.209014 # REPlying: 1392043245.209458 # Received REQuest: 1392043246.211170 # REPlying: 1392043246.211567```

## ZeroMQ PUB/SUB pattern

I'm playing with zeromq for simple network communication/control. For a nice introduction with diagrams see Nicholas Piel's post "ZeroMQ an introduction"

Here is the PUBlish/SUBscribe patter, which usually consist of one Publisher and one or many Subscribers.

Publisher:

```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 import zmq from random import choice import time   context = zmq.Context() socket = context.socket(zmq.PUB) socket.bind("tcp://*:6000") # publish to anyone listening on port 6000   countries = ['netherlands','brazil','germany','portugal'] events = ['yellow card', 'red card', 'goal', 'corner', 'foul']   while True: msg = choice( countries ) +" "+ choice( events ) print "PUB-> " ,msg socket.send( msg ) # PUBlisher only sends messages time.sleep(2)   # output: # PUB-> netherlands yellow card # PUB-> netherlands red card # PUB-> portugal corner```

And here is the Subscriber:

```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 import zmq   context = zmq.Context() socket = context.socket(zmq.SUB) # specific server:port to subscribe to socket.connect("tcp://192.168.3.35:6000") # subscription filter, "" means subscribe to everything socket.setsockopt(zmq.SUBSCRIBE, "")   print "waiting.." while True: msg = socket.recv() # SUBscriber only receives messages print "SUB -> %s" % msg   # typical output: # waiting.. # SUB -> portugal corner # SUB -> portugal foul # SUB -> portugal foul```

## allantools

Python code for calculating Allan deviation and related statistics: https://github.com/aewallin/allantools

Sample output:

Cheat-sheet (see http://en.wikipedia.org/wiki/Colors_of_noise).

 Phase PSD Frequency PSD ADEV 1/f^4, "Black?" 1/f^2, Brownian, random walk sqrt(tau) 1/f^3 ?, "Black?" 1/f Pink constant 1/f^2, Brownian, random walk f^0, White 1/sqrt(tau) f^0, White f^2 Violet 1/tau

Phase PSD is frequency dependence of the the power-spectral-density of the phase noise.
Frequency PSD is the frequency dependence of the power-spectral-density of the frequency noise.
ADEV is the tau-dependence of the Allan deviation.

These can get confusing! Please comment below if I made an error!

## NTP failure detection

Update: Here's how the graph looks like when NTP traffic is enabled again:

A computer that doesn't receive NTP traffic for a while will have its system time drift quite a lot.

In the absence of NTP (UDP port 123) traffic, we can try to roughly ask for the current time over HTTP using wget from google with this shell script:

```#!/bin/sh /usr/bin/wget --no-cache -S -O /dev/null google.com 2&gt;&amp;1 | \ /bin/sed -n -e '/ *Date: */ {' -e s///p -e q -e '}'```

This outputs a string such as "Sat, 23 Nov 2013 09:14:45 GMT"

Now we can put together a python script that calls this shell script, converts the text-format time-stamp into UTC seconds, and compares against the system time. For plotting we then store the time-error in an RRDTool database. This script is called once per minute using cron.

```import subprocess import time import datetime import rrdtool import syslog   args = ['googletime.sh'] datetimestring = subprocess.check_output(args) syslog.syslog( "googletime {0}".format(datetimestring))   # input: Sat, 23 Nov 2013 09:18:02 GMT # output: 1385191082.0 (seconds since 1.1.1970) timestamp = time.mktime( time.strptime(datetimestring, '%a, %d %b %Y %H:%M:%S GMT\n'))   # system time, e.g.: 1385191082.0 stime = time.mktime( time.gmtime() )   # should be zero, if all is well terror = stime-timestamp   # store the measured error in a database datastring = 'N:{0}'.format(str(terror)) # 'N:1234' syslog.syslog( "rrd update: {0}".format(datastring) ) ret = rrdtool.updatev( "time_error.rrd" ,datastring) syslog.syslog( "rrd updatev: {0}".format(ret) )```

Once we have all the values in our time_error.rrd database we can plot them with rrdtool graph. This is what I get:

There is about -4 seconds of drift during 24 hours, or 46 us/s (46 ppm). If the drift is steady we can guess that the computer was on time 14/4 = ~4 days ago.

The script for creating the rrdtool database is this:

```import rrdtool import time   # DS:ds-name:GAUGE | COUNTER | DERIVE | ABSOLUTE:heartbeat:min:max data_sources=[ 'DS:TERROR:GAUGE:70:U:U'] # RRA:AVERAGE | MIN | MAX | LAST:xff:steps:rows   utcsecs = int( time.time() ) pts_day= 24*60 primary = 'RRA:AVERAGE:0.5:1:{0}'.format(pts_day) # 2016 points rrdtool.create( 'time_error.rrd', # filename '--start', str(utcsecs), # when to start '--step', '60', # step between datapoints data_sources, primary)```

And the graph is created by this script. I am using the simple python-rrdtool python bindings - the object-oriented python-pyrrd may have neater syntax and be more pythonic.

```import rrdtool import time   graphname = '/var/www/test.png' day = 24*3600 span = 1*day starttime = int(time.time()-span) endtime = int( time.time() + 0.1*span) updatestamp = time.strftime("%Y-%m-%d %H:%M:%S UTC", time.gmtime(time.time())) graphtitle = '192.168.1.55 System time - google.com time upated: '+updatestamp rrdtool.graph( graphname, '--start', str(starttime), '--end',str(endtime), '--title',graphtitle, '--width',str(1024), '--height',str(600), '--full-size-mode', '--upper-limit',str(20), '--lower-limit',str(-20), '--vertical-label','Error (s)', '--right-axis', '1:0', 'DEF:terror=time_error.rrd:TERROR:AVERAGE', 'LINE2:terror#FF0000')```

## GPIB on Linux

Update 2015 June: additions for Ubuntu 14.04LTS (where a symbolic link is required), and for the NI-USB-B dongle (which needs firmware loading). See also: Using the National Instruments GPIB-USB-B on linux

Notes on how to use GPIB from python on Ubuntu 12.04LTS. Tested with National Instruments USB-GPIB device and a recently bought NI GPIB PCI-E card.

• Get linux-gpib source from http://sourceforge.net/projects/linux-gpib/
• Building the python bindings requires libboost-python (and maybe libboost-python-dev)
• build linux-gpib by running: ./configure, make, sudo make install
• edit /etc/gpib.conf to fit your hardware. The default config board_type = "ni_pci" seems to work with the NI PCI-E card. For the USB-GPIB device we need board_type = "ni_usb_b"
• load the kernel module, for the PCI-E card this is sudo modprobe tnt4882, for the USB-dongle this is sudo modprobe ni_usb_gpib
• On Ubuntu 14.04LTS the library that gpib_config is linked against is installed in a location where the binary doesn't find it. The fix is to add a symbolic link in /lib/
/lib\$ sudo ln -s /usr/local/lib/libgpib.so.0 libgpib.so.0
• configure: sudo gpib_config --minor 0 this reads and parses /etc/gpib.conf so it needs to be re-done if something is changed in gpib.conf
• If gpib_config fails, and you have the older/slower NI-USB-B dongle, then you need to first load the firmware. This seems to not be required with the newer NI-USB-HS (to be confirmed?). Firmware loading requires fxload and the firmware from http://linux-gpib.sourceforge.net/firmware/. Find out where your USB-dongle is with lsusb and then load the firmware with something like:
fxload -D /proc/bus/usb/001/002 -I niusbb_firmware.hex -s niusbb_loader.hex where you replace /001/002 with the numbers for your USB device. Now gpib_config should work.
• By default only members of the gpib-group can use gpib, and only root has access to /etc/gpib0. If we want normal users to use gpib, add them to the gpib group (on Ubuntu the graphical user&group editor is gnome-system-tools)
• Test!
note: If "import gpib" fails in python you might have forgotten to install libboost-python before building linux-gpib. Python bindings will not be built&installed unless you have libboost-python.

I have two devices, so I've configured them like this in gpib.conf

```/* agilent mux/dmm */ device { minor = 0 /* first gpib-card/usb-dongle */ name = "dmm" pad = 8 /* primary address, configure on instrument front-panel*/ sad = 0 /* secondary address, always zero? */ } /* kikusui psu controller */ device { minor = 0 name = "kikusui" pad = 1 sad = 0 }```

And now we can try some python:

```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 import gpib import time dmm= gpib.find("dmm") # define in /etc/gpib.conf kik= gpib.find("kikusui")   def query(handle, command, numbytes=100): gpib.write(handle,command) time.sleep(0.1) response = gpib.read(handle,numbytes) return response   print query(dmm,"*IDN?") print query(kik,"*IDN?")   gpib.close(dmm) gpib.close(kik)   # output: # HEWLETT-PACKARD,34970A,0,13-2-2 # KIKUSUI ELECTRONICS CORP.,PIA4830,0,2.20```

## Jämsä-Jukola statistics

Here's the distribution of total time for teams in the Jämsä-Jukola relay:

Here are the distributions of individual times for the seven legs of the relay.

The histogram from the first leg shows bunching of runners into groups - this is different from all the remaining legs.

The 2nd and 3rd leg look quite similar. Note how the histograms lean to the left - there are many amateurs of varying ability which contribute to the long tail to the right.

The 4th and 5th legs are shorter and faster.

The 6th leg is again longer with the final 7th leg being the longest of them all.

Scraping scripts (requires minor modifications to downloaded results html file): jukola_scrape