In [216]:

    

%pylab inline 

pylab.rcParams['figure.figsize'] = (10.0, 8.0)
import numpy as np
import cv 
import cv2
import pylab as pl
from scipy.stats import linregress


    

    




Populating the interactive namespace from numpy and matplotlib

In [217]:

    

data = np.genfromtxt("data.csv", delimiter=',')
data = (data[1:,1:]).astype(np.float32)
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
_, labs_masked, centre = cv2.kmeans(data, 2, criteria, 10, cv2.KMEANS_PP_CENTERS)
day_idx = np.argmax(centre[[0,1], [0,0]])


    

In [218]:

    

def my_hist(data, col):
    fig = pl.figure()
    pl.hist(data, bins=50, color=col, normed=True)
    m, q1, q3 =  np.percentile(data, [50, 25, 75])
    pl.xlim(0, 255)
    pl.ylim(0, 0.06)
    pl.xlabel("pixel value")
    pl.ylabel("frequency")
    bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=0.9)   
    pl.text(150, 0.05,  "Median = %i\nQ1 = %i\nQ3=%i" % (m, q1, q3), ha="left", va="center", size=18, bbox=bbox_props)
    return fig


    

In [219]:

    

day_img = cv2.imread("day_shot.png")
fig = my_hist(day_img.flatten(), 'r')
pl.savefig("day_hist.svg")


    

In [220]:

    

night_img = cv2.imread("night_shot.png")
fig = my_hist(night_img.flatten(), 'k')
pl.savefig("night_hist.svg")


    

In [220]:

    




    

In [221]:

    

pl.figure()
pl.plot(data[:,0], data[:,2] -data[:,1] , "o",alpha=0.1, color=(0.5,0.5,0.5))
pl.ylim((50, 80))
pl.xlim((60, 130))
pl.ylabel("Interquartile range")
pl.xlabel("Median")
pl.savefig("before_clust.svg")


    

In [221]:

    




    

In [221]:

    




    

In [222]:

    

#colours = ['r' if c == day_idx else 'b' for c in labs_masked.flatten()]
colours = np.where(labs_masked.flatten() == day_idx, 'r', 'b')
pl.figure()
for i in (0,1):
    filt = labs_masked.flatten() == i
    x = data[filt,0]
    y = data[filt,2] - data[filt,1]
    if i == day_idx:            
        col = 'r'
    else:
        col= 'k'
        
    pl.plot(x, y , "o",alpha=0.1, color=col)
    
pl.ylim((50, 80))
pl.xlim((60, 130))
pl.ylabel("Interquartile range")
pl.xlabel("Median")
pl.savefig("after_clust.svg")


    

In [225]:

    

x = labs_masked.flatten()
result = np.correlate(x, x, mode = 'full')
#pl.acorr()


maxcorr = np.argmax(result)
# print 'maximum = ', result[maxcorr]
result = result / result[maxcorr]
result[maxcorr]
observed_time = np.where(np.diff(x == 1))[0]

expected_time = np.arange(np.min(observed_time), np.max(observed_time) + 720, 720)
pl.figure()
pl.plot(expected_time/60, observed_time/60, 'x')
pl.plot(expected_time/60,expected_time/60, '-')
pl.ylabel("Expected transition time(h)")
pl.xlabel("Observed transition time(h)")
slope = linregress(expected_time, observed_time)[0]
min_per_day = (24 - slope  * 24) * 60

bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=0.9)   
pl.text(10, 120,  "slope = %.03f\nSo we drift %.01f min/day!" % (slope, min_per_day), ha="left", va="center", size=18, bbox=bbox_props)

pl.savefig("time_drift.svg")


pl.show()


    

In [223]:

    




    

In [ ]: