+ 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)