TCS - Improve process.py script

Going to use this for everything to do with data analysis
Abandoning the terrible plot.sh scripts due to fits of rage induced by bash special character escaping

diff --git a/research/TCS/process.py b/research/TCS/process.py
index ff186fb..31c98d2 100755 (executable)
--- a/research/TCS/process.py
+++ b/research/TCS/process.py
@@ -31,41 +31,85 @@ def GetData(filename):
 def DoNothing(data):
        return data
 
-def GetDataSets(directory=".", function=DoNothing):
+def AverageAllDataSets(directory=".", function=DoNothing):
        dirs = {}
        for f in os.listdir(directory):
                if os.path.isdir(directory+"/"+str(f)) == True:
                        data_set = []
                        for datafile in os.listdir(directory+"/"+str(f)):
                                if datafile.split(".")[1] == "dat":
-                                       data_set.append(function(map(lambda e : [e[1], e[2]], GetData("./"+str(f)+"/"+str(datafile)))))
+                                       data_set.append(GetData(f))
 
                        avg = Average(data_set)
                        dirs.update({f : avg})
        return dirs
 
+def GetDataSets(directory="."):
+       data_sets = []
+       for f in os.listdir(directory):
+               if os.path.isdir(directory+"/"+str(f)) == False:
+                       if (len(f.split(".")) > 1 and f.split(".")[1] == "dat"):
+                               d = GetData(directory+"/"+str(f))
+                               if len(d) > 0:
+                                       data_sets.append(d)
+       return data_sets        
+
 
 
-def Derivative(data, a=0, b=1):
-       result = []
+def Derivative(data, a=1, b=2, sigma=None,step=1):
+       result = [[]]
        n = 0
-       dI = 0
-       dE = 0
-       for i in range(1, len(data)-1):
-               dE = data[i+1][a] - data[i][a]
-               if (dE != 0):
-                       n = 0
-                       dI = 0
-               
-               n += 1
-               dI += data[i+1][b] - data[i][b]                 
-               if (dE != 0):                   
-                       result.append([data[i][a], (dI / (n * dE)) ] ) #/ data[i][2]])
-       return result
+       dI = [0,0]
+       dE = [0,0]
+       
+       for i in range(0, len(data),step):
+               result[len(result)-1] = [d for d in data[i]]
+               if (i >= step):
+                       dE[0] = data[i][a] - data[i-step][a]
+                       dI[0] = data[i][b] - data[i-step][b]
+               else:
+                       dI[0] = None
 
-def MaxNormalise(data, u=1):
+               if (i < len(data)-step):
+                       dE[1] = data[i+step][a] - data[i][a]
+                       dI[1] = data[i+step][b] - data[i][b]
+               else:
+                       dI[1] = None            
 
-       
+               if sigma != None:
+                       #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))])
+                       if sigma < 0:
+                               if dI[0] != None: dI[0] -= 0.5*data[i][int(abs(sigma))]
+                               if dI[1] != None: dI[1] -= 0.5*data[i][int(abs(sigma))]
+                       else:
+                               if dI[0] != None: dI[0] += 0.5*data[i][int(abs(sigma))]
+                               if dI[1] != None: dI[1] += 0.5*data[i][int(abs(sigma))]
+
+               deltaE = 0.0
+               deltaI = 0.0
+               count = 0
+               if dI[0] != None:
+                       deltaE += dE[0]
+                       deltaI += dI[0]
+                       count += 1
+               if dI[1] != None:
+                       deltaE += dE[1]
+                       deltaI += dI[1]
+                       count += 1
+
+               deltaI /= float(count)
+               deltaE /= float(count)
+
+
+               if (deltaE != 0):       
+                       result[len(result)-1][b] = (deltaI / deltaE)
+               else:
+                       result[len(result)-1][b] = 0.0
+               result.append([])
+                       
+       return result[0:len(result)-1]
+
+def MaxNormalise(data, u=1):   
        result = copy.deepcopy(data)
        if (len(data) <= 0):
                return result
@@ -80,21 +124,19 @@ def Average(data_sets, u=1):
        avg = odict.odict()
        for t in data_sets:
                for p in t:
-                       if p[0] in avg:
-                               avg[p[0]][0] += p[u]
-                               avg[p[0]][1] += 1
+                       if p[u] in avg:
+                               #print "Already have " + str(p[u])
+                               avg[p[u]][1] += 1
+                               for i in range(0, len(p)):
+                                       avg[p[u]][0][i] += p[i]
                        else:
-                               avg.update({p[0] : [p[u], 1]})
+                               #print "Create key for " + str(p[u])
+                               avg.update({p[u] : [p, 1]})
 
        for a in avg.keys():
-               avg[a] = float(avg[a][0]) / float(avg[a][1])
-       return sorted(avg.items(), key = lambda e : e[0])
-
-def Plot(*args):
-       gnuplot.plot(args)
-
-def FitTCS(data):
-       pass
+               for i in range(0, len(avg[a][0])):
+                       avg[a][0][i] /= float(avg[a][1])
+       return map(lambda e : e[1][0], sorted(avg.items(), key = lambda e : e[0]))
 
 def FullWidthAtHalfMax(data, u=1):
        maxval = max(data, key = lambda e : e[u])
@@ -119,8 +161,92 @@ def FullWidthAtHalfMax(data, u=1):
 def SaveData(filename, data):
        out = open(filename, "w", 0)
        for a in data:
-               out.write(str(a[0]) + "\t" + str(a[1]) + "\n")
+               for i in range(0, len(a)):
+                       out.write(str(a[i]))
+                       if (i < len(a) - 1):
+                               out.write("\t")
+                       out.write("\n")
+
+def CalibrateData(data, ammeter_scale=1e-6):
+       for i in range(0, len(data)):
+               data[i][1] = 16.8 * float(data[i][1]) / 4000.0
+               data[i][2] = ammeter_scale * 0.170 * float(data[i][2]) / 268.0
+               data[i][3] = ammeter_scale * 0.170 * float(data[i][3]) / 268.0
+       return data
+
+def ShowTCS(filename, calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,w="lp", step=1):
+       if type(filename) == type(""): 
+               data = GetData(filename)
+       else:
+               data = filename
+               filename = "data"
+
+       if calibrate: 
+               data = CalibrateData(data)
+               units = ["V", "uA / V"]
+       else:
+               units = ["DAC counts", "ADC counts / DAC counts"]
 
+       if not normalise:
+               gnuplot("set ylabel \"dI(E)/dE ("+str(units[1])+")\"")
+       else:
+               data = MaxNormalise(data)
+               gnuplot("set ylabel \"dI(E)/dE (normalised)\"")
+
+       gnuplot("set title \"S(E)\"")
+       gnuplot("set xlabel \"U ("+str(units[0])+")\"")
+
+
+       d = Derivative(data, 1, 2, step=step)
+       
+       plot(Gnuplot.Data(d, using="2:3", with_=w,title="S(E) : " + str(filename)))
+       if (show_error):
+               error1 = Derivative(data, 1, 2, -3,step=step)
+               error2 = Derivative(data, 1, 2, +3,step=step)
+               gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
+               gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
+
+       return data
+
+def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,w="lp", step=1):
+       if type(filename) == type(""): 
+               data = GetData(filename)
+       else:
+               data = filename
+               filename = "data"
+
+       if len(data) <= 0:
+               return
+       if calibrate: 
+               data = CalibrateData(data)
+               units = ["V", "uA"]
+       else:
+               units = ["DAC counts", "ADC counts"]
+
+       if not normalise:
+               gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
+       else:
+               data = MaxNormalise(data)
+               gnuplot("set ylabel \"I(E) (normalised)\"")
+
+       gnuplot("set title \"S(E)\"")
+       gnuplot("set xlabel \"U ("+str(units[0])+")\"")
+
+
+       #d = Derivative(data, 1, 2, step=step)
+       
+       plot(Gnuplot.Data(data, using="2:3", with_=w,title="S(E) : " + str(filename)))
+       if (show_error):
+               error1 = copy.deepcopy(data)
+               error2 = copy.deepcopy(data)
+               for i in range(len(data)):
+                       #print str(data[i])
+                       error1[i][2] -= 0.50*float(data[i][3])
+                       error2[i][2] += 0.50*float(data[i][3])
+               gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
+               gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
+               
+       return data
 
 def main():    
        if (len(sys.argv) < 2):
@@ -136,7 +262,7 @@ def main():
        #gnuplot("set term postscript colour")
        #gnuplot("set output \"test.eps\"")
        for i in range(1, len(sys.argv)):
-               tcs.append(DoNothing(map(lambda e : [e[1], e[2]], GetData(sys.argv[i]))))
+               tcs.append(Derivative(GetData(sys.argv[i]), 1, 2))
                #tcs.append(GetTCS(GetData(sys.argv[i])))
                if (len(tcs[i-1]) > 0):
                        gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="lp"))
@@ -146,7 +272,9 @@ def main():
 
        
        avg = Average(tcs)
-       #gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2"))
+       for a in avg:
+               sys.stdout.write(str(a[0]) + "\t" + str(a[1]) + "\t" + str(a[1]) + "\n")
+       gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2"))
        
        sys.stdout.write("Save averaged data as (blank for no save): ")
        filename = sys.stdin.readline().strip(" \r\n\t")