From 598da84055df9d274f4ac2c492f9323bbc0df26f Mon Sep 17 00:00:00 2001 From: Sam Moore Date: Mon, 15 Oct 2012 10:45:38 +0800 Subject: [PATCH] TCS - Commit process.py before adding numpy BALRGSHA --- research/TCS/process.py | 550 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 503 insertions(+), 47 deletions(-) diff --git a/research/TCS/process.py b/research/TCS/process.py index a5f21f0c..88a536fb 100755 --- a/research/TCS/process.py +++ b/research/TCS/process.py @@ -16,6 +16,11 @@ import copy import Gnuplot, Gnuplot.funcutils import string +import time +import math +import cmath +import random +import numpy gnuplot = Gnuplot.Gnuplot() @@ -46,15 +51,28 @@ def DirectoryName(f, start=0,back=1): a = f.split("/") return string.join(a[start:(len(a)-back)], "/") -def GetData(filename): +def GetData(filename, key=1): input_file = open(filename, "r") - data = [] + data = {} for line in input_file: line = re.sub("#.*", "", line).strip("\r\n ") if len(line) == 0: continue - data.append(map(lambda e : float(e), line.split("\t"))) - return data + + line = map(lambda e : float(e), line.split("\t")) + + if line[key] in data: + for i in range(len(line)): + data[line[key]][0][i] += line[i] + data[line[key]][1] += 1 + else: + data.update({line[key] : [line, 1]}) + + d = map(lambda e : map(lambda f : float(f) / float(e[1][1]), e[1][0]), data.items()) + d.sort(key = lambda e : e[key]) + #for l in d: + # print str(l) + return d def DoNothing(data): return data @@ -125,12 +143,15 @@ def Derivative(data, a=1, b=2, sigma=None,step=1): deltaI += dI[1] count += 1 - deltaI /= float(count) - deltaE /= float(count) + if (count > 0): + deltaI /= float(count) + deltaE /= float(count) - if (deltaE != 0): - result[len(result)-1][b] = (deltaI / deltaE) + if (deltaE != 0): + result[len(result)-1][b] = (deltaI / deltaE) + else: + result[len(result)-1][b] = 0.0 else: result[len(result)-1][b] = 0.0 result.append([]) @@ -195,9 +216,9 @@ def SaveData(filename, data): out.write("\t") out.write("\n") -def AverageAllData(directory, save=None) +def AverageAllData(directory, save=None): data_sets = [] - if save == None: save = directory+"/average.dat": + if save == None: save = directory+"/average.dat" for d in FindDataFiles(directory): data_sets.append(GetData(d)) @@ -213,7 +234,12 @@ def CalibrateData(original, ammeter_scale=1e-6): 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,with_="lp", step=1, output=None, title=""): +def ShowTCS(filename, raw=True,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,with_="lp", step=1, output=None, title="", master_title="", smooth=0, show_peak=True, inflection=1): + + if raw == False: + calibrate = False + normalise = False + if type(filename) == type(""): data = GetData(filename) else: @@ -226,6 +252,14 @@ def ShowTCS(filename, calibrate=True, normalise=False, show_error=False, plot=gn if (len(data) <= 0): return data + if (smooth > 0): + if type(smooth) == type([]): + for i in range(smooth[0]): + data = Smooth(data, m=smooth[1]) + else: + data = Smooth(data, m=smooth) + + if calibrate: data = CalibrateData(data) units = ["V", "uA / V"] @@ -242,24 +276,54 @@ def ShowTCS(filename, calibrate=True, normalise=False, show_error=False, plot=gn gnuplot("set term png size 640,480") gnuplot("set output \""+str(output)+"\"") - gnuplot("set title \"Total Current Spectrum S(E)\"") + if master_title == "": + master_title = "Total Current Spectrum S(E)" + if type(filename) == type("") and plot == gnuplot.plot: + if filename != "tcs data": + p = ReadParameters(filename) + if "Sample" in p: + master_title += "\\nSample: "+p["Sample"] + + gnuplot("set title \""+str(master_title)+"\"") gnuplot("set xlabel \"U ("+str(units[0])+")\"") - d = Derivative(data, 1, 2, step=step) + if raw: + d = Derivative(data, 1, 2, step=step) + else: + d = data + + ymax = 0.01 + 1.2 * max(d, key=lambda e : e[2])[2] + ymin = -0.01 + 1.2 * min(d, key=lambda e : e[2])[2] + gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]") + + plotList = [] + plotList.append(Gnuplot.Data(d, using="2:3", with_=with_,title=title)) - plot(Gnuplot.Data(d, using="2:3", with_=with_,title=title)) 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="-sigma/2")) - gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="+sigma/2")) + plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="-sigma/2")) + plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="+sigma/2")) + if (show_peak): + peak = SmoothPeakFind(d, ap=DoNothing, stop=1, inflection=inflection) + plotList += PlotPeaks(peak,with_="l lt -1", plot=None) + + + + if (plot != None): + plot(*plotList) + time.sleep(0.2) + if (output != None and type(output) == type("")): gnuplot("set term wxt") + + if (plot == None): + return plotList return data -def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,with_="lp", step=1, output=None, title=""): +def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,with_="lp", step=1, output=None, title="", master_title="Sample Current I(E)", smooth=0): if type(filename) == type(""): data = GetData(filename) else: @@ -269,8 +333,14 @@ def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gn if (title == ""): title = BaseName(filename) + if len(data) <= 0: return data + + + if (smooth > 0): + data = Smooth(data) + if calibrate: data = CalibrateData(data) units = ["V", "uA"] @@ -287,13 +357,16 @@ def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gn gnuplot("set term png size 640,480") gnuplot("set output \""+str(output)+"\"") - gnuplot("set title \"Sample Current I(E)\"") + gnuplot("set title \""+str(master_title)+"\"") gnuplot("set xlabel \"U ("+str(units[0])+")\"") #d = Derivative(data, 1, 2, step=step) + + plotList = [] - plot(Gnuplot.Data(data, using="2:3", with_=with_,title=title)) + plotList.append(Gnuplot.Data(data, using="2:3", with_=with_,title=title)) + time.sleep(0.1) if (show_error): error1 = copy.deepcopy(data) error2 = copy.deepcopy(data) @@ -301,47 +374,430 @@ def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gn #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")) + plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound")) + plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound")) + + if plot != None: + + plot(*plotList) + if (output != None and type(output) == type("")): + gnuplot("set term wxt") + return data + else: + return plotList + +def ReadParameters(filename): + parameters = odict.odict() + input_file = open(filename, "r") + for line in input_file: + k = line.split("=") + item = None + #print str(k) + if (len(k) >= 2): + item = k[0].strip("# \r\n") + value = k[1].strip("# \r\n") + if (item in parameters): + parameters[item] = value + else: + parameters.update({str(item) : value}) + input_file.close() + return parameters + +def PlotParameters(filename): + ReadParameters(filename) + +def Smooth(data, m, k=2): + smooth = copy.deepcopy(data) + for i in range(len(smooth)): + count = 0 + smooth[i][k] = 0.0 + for j in range(i-m,i+m): + if j >= 0 and j < len(smooth): + count += 1 + smooth[i][k] += data[j][k] + if count > 0: + smooth[i][k] = smooth[i][k] / float(count) + else: + smooth[i][k] = data[i][k] + + return smooth + +def PeakFind(data, k=2,threshold=0.00, inflection=0): + results = [] + for i in range(len(data)): + if i == 0 or i == len(data)-1: + continue + #if abs(data[i][k]) < threshold * abs(max(data, key = lambda e : abs(e[k]))[k]): + # continue + + left = data[i-1][k] - data[i][k] + right = data[i+1][k] - data[i][k] + if abs(left) < threshold*abs(data[i][k]): + continue + if abs(right) < threshold*abs(data[i][k]): + continue + if left*right > 0: + results.append(data[i] + [inflection]) + + if inflection > 0: + results += PeakFind(Derivative(data), k=k, threshold=threshold, inflection=inflection-1) + + return results + +def SmoothPeakFind(data, a=1, k=2, ap=DoNothing, stop=10,smooth=5, inflection=0): + s = data + #results = [] + + peakList = [] + + m = 0 + while m < stop: + #results.append([]) + peaks = PeakFind(ap(s),k=k, inflection=inflection) + #print "m = " +str(m) + for p in peaks: + add = [m] + [add.append(f) for f in p] + + if m == 0: + #print "*New peak at " + str(p) + peakList.append([add]) + else: + score = [] + for i in range(len(peakList)): + p2 = peakList[i][len(peakList[i])-1] + if m - p2[0] > 1: + continue + score.append([i, abs(p[a] - p2[1+a])]) + + score.sort(key = lambda e : e[1]) + if len(score) == 0 or score[0][1] > 100: + #print "New peak at " + str(p) + peakList.append([add]) + else: + #print "Peak exists near " + str(p) + " ("+str(score[0][1])+") " + str(peakList[score[0][0]][len(peakList[score[0][0]])-1]) + peakList[score[0][0]].append(add) + + + + #results.append([m, []]) + #[results[len(results)-1].append(f) for f in p] + m += 1 + s = Smooth(s, m=smooth,k=k) + + #results.sort(key = lambda e : e[2]) + + #peaks = [] + return peakList + + + +def PlotPeaks(peaks, calibrate=True, with_="lp", plot=gnuplot.replot): + + plotList = [] + for p in peaks: + p.append(copy.deepcopy(p[len(p)-1])) + + p[len(p)-1][0] += 1 + + + #print "Adding " + str(p) + " to list" + if len(p) >= 0: + l = p[len(p)-1] + if l[len(l)-1] < 1: + with_ = with_.split(" lt")[0] + " lt 9" + plotList.append(Gnuplot.Data(p, using="3:1", with_=with_)) + - if (output != None and type(output) == type("")): - gnuplot("set term wxt") - return data + if len(plotList) > 0 and plot != None: + plot(*plotList) + time.sleep(0.2) + + #print str(plotList) + #for p in peaks: + # p = p[0:len(p)-1] + return plotList + def main(): - return 0 + if (len(sys.argv) < 2): sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n") return 1 - tcs = [] - gnuplot("set style data lp") - gnuplot("set key outside right") - #gnuplot("set title \"Au on Si (50min 3.5A 3-6 e-8mbar)\"") - #gnuplot("set xlabel \"E (DAC Counts)\"") - #gnuplot("set ylabel \"S(E) (ADC/DAC Counts)\"") - #gnuplot("set term postscript colour") - #gnuplot("set output \"test.eps\"") - for i in range(1, len(sys.argv)): - if (len(tcs[i-1]) > 0): - gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="lp")) - - # Now average the data + i = 1 + plotFunc = ShowTCS + normalise = False + title = "" + master_title = "" + while i < len(sys.argv): + if sys.argv[i] == "--raw": + plotFunc = ShowData + elif sys.argv[i] == "--tcs": + plotFunc = ShowTCS + elif sys.argv[i] == "--output": + if i+1 >= len(sys.argv): + sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n") + sys.exit(1) + gnuplot("set term postscript colour") + gnuplot("set output \""+sys.argv[i+1]+"\"") + i += 1 + elif sys.argv[i] == "--wxt": + gnuplot("set term wxt") + elif sys.argv[i] == "--normalise": + normalise = True + elif sys.argv[i] == "--unnormalise": + normalise = False + elif sys.argv[i] == "--title": + if i+1 >= len(sys.argv): + sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n") + sys.exit(1) + title = sys.argv[i+1] + i += 1 + elif sys.argv[i] == "--master_title": + if i+1 >= len(sys.argv): + sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n") + sys.exit(1) + master_title = sys.argv[i+1] + i += 1 + else: + plotFunc(sys.argv[i], plot=gnuplot.replot, normalise=normalise, title=title, master_title=master_title) + + i += 1 + + print "Done. Press enter to exit, or type name of file to save as." + out = sys.stdin.readline().strip("\t\r\n #") + if out != "": + gnuplot("set term postscript colour") + gnuplot("set output \""+out+"\"") + gnuplot.replot() + + +def ModelTCS(f, sigma, Emin, Emax, dE): + data = [] + E = Emin + while E < Emax: + S = (1 - sigma(0))*f(-E) + FuncIntegrate(lambda e : f(e - E) * FuncDerivative(sigma, E, dE), Emin, Emax, dE) + data.append([0.00, E, S,0.00]) + E += dE + return data + +def IntegrateTCS(data, imin, imax=0, di=1): + i = imin + if imax == 0: + imax = len(data)-1 + total = 0.0 + while i < imax: + total += data[i][2] * (data[i+1][1] - data[i][1]) + i += di + return total + +def FuncIntegrate(f, xmin, xmax, dx): + x = xmin + total = 0.0 + while x <= xmax: + total += f(x) * dx + x += dx + return total + +def FuncDerivative(f, x, dx): + return 0.50*(f(x+dx) - f(x-dx))/dx + +def FitTCS(data, min_mse=1e-4, max_fail=100, max_adjust=4,divide=10, plot=gnuplot.plot,smooth=0): + if type(data) == type(""): + d = GetData(data) + d = CalibrateData(d) + d = MaxNormalise(d) + d = Derivative(d) + else: + d = data + + if smooth != 0: + if type(smooth) == type([]): + for _ in range(smooth[0]): + d = Smooth(d, m=smooth[1]) + else: + d = Smooth(d, m=smooth) - avg = Average(tcs) - 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")) + + plotItems = ShowTCS(d, raw=False,smooth=smooth,plot=None) + plotItems.append(None) - sys.stdout.write("Save averaged data as (blank for no save): ") - filename = sys.stdin.readline().strip(" \r\n\t") - if (filename != ""): - SaveData(filename, avg) + peaks = SmoothPeakFind(d, smooth=5, stop=1, inflection=0) + peaks.sort(key = lambda e : e[len(e)-1][1]) + + fits = [] + + for i in range(0,len(peaks)): - return 0 + p = peaks[i] + l = p[len(p)-1] + if l[len(l)-1] == 0: + fits.append([l[3], l[2], 1.0]) + else: + if i-2 >= 0: + l = peaks[i-2][len(peaks[i-2])] + if l[len(l)-1] == 0: + fits.append([l[3], l[2], 1.0]) + if i+2 <= len(peaks)-1: + l = peaks[i+2][len(peaks[i+2])] + if l[len(l)-1] == 0: + fits.append([l[3], l[2], 1.0]) + + for i in range(len(fits)): + left = 2.0 + right = 2.0 + if i > 0: + left = fits[i-1][1] - fits[i][1] + if i < len(fits)-1: + right = fits[i+1][1] - fits[i][1] + + fits[i][2] = min([abs(0.5*left), abs(0.5*right)]) + + + #print "Fits are " + str(fits) + #stdin.readline() + + + + def tcs(E): + total = 0.0 + for f in fits: + dt = f[0] * gaussian(E - f[1], f[2]) + #print " Increase total by " + str(dt) + total += dt + #print "tcs returns " + str(total) + return total + + mse = 1 + old_mse = 1 + cycle = 0 + failcount = 0 + + + adjust = 1.0 + + iterations = 0 + + while failcount < max_fail and mse > min_mse: + i = random.randint(0, len(fits)-1) + j = random.randint(0, len(fits[i])-1) + #while j == 1: + # j = random.randint(0, len(fits[i])-1) + + #print "Adjust " + str(i) + ","+str(j) + ": Iteration " + str(iterations) + " mse: " + str(mse) + + old = fits[i][j] + old_mse = mse + + fits[i][j] += adjust * (random.random() - 0.50) + if i == 2: + while fits[i][j] <= 0.0005: + fits[i][j] = adjust * (random.random() - 0.50) + + + model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, divide*d[len(d)-1][1]/(len(d))) + mse = MeanSquareError(model, d[0::divide]) + if mse >= old_mse: + fits[i][j] = old + failcount += 1 + if failcount > max_fail / 2: + if adjust > 1.0/(2.0**max_adjust): + adjust /= 2 + mse = old_mse + else: + #adjust /= 2.0 + failcount = 0 + + + iterations += 1 + + + #model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, 16.8/len(d)) + plotItems[len(plotItems)-1] = Gnuplot.Data(model, using="2:3", with_="l lt 3", title="model") + time.sleep(0.1) + + fits.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True) + + if plot != None: + gnuplot("set title \"MSE = "+str(mse)+"\\nfailcount = "+str(failcount)+"\\nadjust = "+str(adjust)+"\"") + gnuplot.plot(*plotItems) + + return [fits, model] + else: + return [fits, model,plotItems] + + + #return model + +def SaveFit(filename, fit): + out = open(filename, "w", 0) + out.write("# TCS Fit\n") + + for f in fit: + out.write(str(f[0]) + "\t" + str(f[1]) + "\t" + str(f[2]) + "\n") + + out.close() + + +def LoadFit(filename): + infile = open(filename, "r") + if (infile.readline() != "# TCS Fit\n"): + sys.stderr.write("Error loading fit from file " + str(filename) + "\n") + sys.exit(0) + + fit = [] + while True: + f = infile.readline().strip("# \r\n\t").split("\t") + if len(f) != 3: + break + fit.append(map(float, f)) + + infile.close() + fit.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True) + #def model(e): + # total = 0.0 + # for f in fit: + # total += f[0] * gaussian(e - f[1], f[2]) + # return total + + #return table(lambda e : [0.0, e, model(e), 0.0], 0.0, 16.8, 16.8/400) + return fit + +def MeanSquareError(model, real, k = 2): + mse = 0.0 + for i in range(len(real)): + + mse += (model[i][k] - real[i][k])**2 + + mse /= len(model) + return mse + +def delta(x): + if (x == 0): + return 1.0 + else: + return 0.0 + +def table(f, xmin, xmax, dx): + result = [] + x = xmin + while (x <= xmax): + result.append(f(x)) + x += dx + return result + +def gaussian(x, sigma): + if (sigma == 0.0): + return 0.0 + return math.exp(- (x**2.0)/(2.0 * sigma**2.0)) / (sigma * (2.0 * math.pi)**0.50) +def step(x, sigma, T): + if T == 0: + return 1.0 + return 1.0 / (math.exp((x - sigma)/T) + 1.0) if __name__ == "__main__": sys.exit(main()) -- 2.20.1