#Fabio Eboli 2012

from numpy import arange,array,ones
import pylab
from pylab import *
from scipy import stats
import time

#pylab.ion()
fig=pylab.figure(1)
sp1=fig.add_subplot(221)
sp2=fig.add_subplot(222)
sp3=fig.add_subplot(223)
sp4=fig.add_subplot(224)

manager=pylab.get_current_fig_manager()

data=[0]
xdata=[0]
datapos=[0]
correctionvalues=[0]
dataunstep=[0]
dataundrift=[0]
correctionvalues=[0]
driftslope=0
driftintercept=0
numsamples=0
outliers=0

p1=sp1.plot(xdata,data,"r.")
sp1.set_title("Raw Data, outliers removed")
p2=sp2.plot(xdata,correctionvalues,"b")
sp2.set_title("Step correction values")
p3=sp3.plot(xdata,dataunstep,"g")
p31=sp3.plot([0,len(dataunstep)],[driftintercept,driftintercept+driftslope*len(dataunstep)],"g")
sp3.set_title("Steps removed and lin. drift est")
p4=sp4.plot(xdata,dataundrift,"y")
sp4.set_title("Drift corrected")

def Calculations(arg):
    global data
    global datapos
    global xdata
    global dataunstep
    global dataundrift
    global driftslope
    global driftintercept
    global correctionvalues
    global numsamples
    global outliers
    data=[]
    xdata=[]
    datapos=[]
    correctionsteps=[]
    correctionvalues=[]
    dataundrift=[]
    dataunstep=[]
    stepvalue=1e-7
    threshold=5e-8
    correctionstep=0
    outliers=0
    logfile=open("..\\TI.txt")
    header=logfile.readline()
    for lines in logfile:
        data.append(float(lines))
    logfile.close()
    numsamples=len(data)
    xdata=arange(0,numsamples)
    print(str(numsamples)+" Samples")
    print("Removing negative values...")
    for i in data:
        if i<0 :
            datapos.append(i+1e-7)
        else :
            datapos.append(i)
    print("Estimating steps correction and removing ouliers...")
    xi=arange(1,len(datapos))
    correctionsteps.append(0)
    correctionvalues.append(0)
    dataunstep.append(datapos[0])
    for i in xi:
        olddata=datapos[i-1]
        newdata=datapos[i]
        diff=newdata-olddata
        correctionstep=0
        if diff>0:
            if diff> 3*threshold:
                datapos[i]=datapos[i-1] # filetring out values
                outliers=outliers+1
            else:
                while diff>threshold:
                    correctionstep=correctionstep-1
                    diff=newdata+correctionstep*stepvalue-olddata
        else :
            if diff< -3*threshold:
                datapos[i]=datapos[i-1] # filetring out values
                outliers=outliers+1
            else:
                while diff<-threshold:
                    correctionstep=correctionstep+1
                    diff=newdata+correctionstep*stepvalue-olddata
        correctionsteps.append(correctionstep)
        correctionvalues.append(correctionvalues[i-1]+correctionstep*stepvalue)
        dataunstep.append(datapos[i]+correctionvalues[i])
    xi=arange(0,len(dataunstep))
    print("Drift estimation:")
    driftline=stats.linregress(xi,dataunstep)
    driftslope=driftline[0]
    driftintercept=driftline[1]
    print("Slope="+str(driftslope)+"; Intercept:"+str(driftintercept))
    print("Drift correction...")
    for i in xi:
        dataundrift.append(dataunstep[i]-driftslope*i)

def Plotter(arg):
    global datapos
    global dataunstep
    global dataundrift
    global driftslope
    global driftintercept
    global correctionvalues
    global numsamples
    global outliers
    p1[0].set_data(arange(0,len(datapos)),datapos)
    sp1.axis([0,len(datapos),min(datapos),max(datapos)])
    sp1.set_title("Raw Data: "+str(numsamples)+"samples, "+str(outliers)+" outliers removed")
    p2[0].set_data(arange(0,len(correctionvalues)),correctionvalues)
    sp2.axis([0,len(correctionvalues),min(correctionvalues),max(correctionvalues)])
    sp2.set_title("Step correction values")
    p3[0].set_data(arange(0,len(dataunstep)),dataunstep)
    p31[0].set_data([0,len(dataunstep)],[driftintercept,driftintercept+driftslope*len(dataunstep)])
    sp3.axis([0,len(dataunstep),min(dataunstep),max(dataunstep)])
    sp3.set_title("Steps removed and lin. drift est. : "+str(driftslope))
    p4[0].set_data(arange(0,len(dataundrift)),dataundrift)
    sp4.axis([0,len(dataundrift),min(dataundrift),max(dataundrift)])
    sp4.set_title("Drift Removed")
    manager.canvas.draw()

Calculations(0)
Plotter(0)

#time.sleep(20)
timer1=fig.canvas.new_timer(interval=50000)
timer1.add_callback(Calculations,())
timer=fig.canvas.new_timer(interval=50000)
timer.add_callback(Plotter,())
timer.start()
timer1.start()
pylab.show()