X-Git-Url: https://git.ucc.asn.au/?a=blobdiff_plain;f=research%2FTCS%2Fprocess.py;h=afb38c67591a75fdbca7ec7961496866701f6078;hb=8c40cf6d18b192436651eb671b2eacd761d98ae9;hp=a5f21f0cbbd8f3425b7f3887643c006ec19f1e91;hpb=4ac4dae6c58dda0530d237f412241c0ce3f7c1a9;p=matches%2Fhonours.git diff --git a/research/TCS/process.py b/research/TCS/process.py index a5f21f0c..afb38c67 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,31 +51,46 @@ 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): + + if type(filename) != type(""): + if type(filename) == type([]): + return filename + else: + return [[0,0,0,0]] + + + if os.path.isdir(filename): + if os.path.exists(filename.strip("/")+"/average.dat"): + os.remove(filename.strip("/")+"/average.dat") + AverageAllData(filename) + return GetData(filename.strip("/")+"/average.dat") + 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 -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(GetData(f)) - - avg = Average(data_set) - dirs.update({f : avg}) - return dirs def GetDataSets(directory="."): data_sets = [] @@ -125,12 +145,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([]) @@ -142,6 +165,9 @@ def MaxNormalise(data, u=2): if (len(data) <= 0): return result maxval = max(data, key = lambda e : e[u])[u] + + if maxval == 0: + return result for d in result: d[u] = d[u] / maxval @@ -195,11 +221,14 @@ def SaveData(filename, data): out.write("\t") out.write("\n") -def AverageAllData(directory, save=None) +def AverageAllData(directory, save=None, normalise=True): data_sets = [] - if save == None: save = directory+"/average.dat": - for d in FindDataFiles(directory): - data_sets.append(GetData(d)) + if save == None: save = directory+"/average.dat" + for f in FindDataFiles(directory): + d = GetData(f) + if normalise: + d = MaxNormalise(d) + data_sets.append(d) a = Average(data_sets) SaveData(save, a) @@ -213,7 +242,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=False, inflection=1): + + if raw == False: + calibrate = False + normalise = False + if type(filename) == type(""): data = GetData(filename) else: @@ -226,6 +260,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 +284,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 +341,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 +365,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 +382,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 (output != None and type(output) == type("")): - gnuplot("set term wxt") - return data + 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 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 #lambda e : ShowTCS(e, show_peak=False) + 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())