Several minor changes to comb tracking, added comb finding user script, added known comb plotting

chplot.py

@@ -17,6 +17,7 @@ def makeparser():
 	parser.add_argument("--verbose",help="Verbose output for errors (default=False)",type='bool', default=False)
 	parser.add_argument("--truefcutoff",help="Truncate the frequency range (not just limiting the x-axis of initial plotting). Good for high resolution. (default=False)",type='bool',default=False)
 	parser.add_argument("--knownLinesFile",help="A file containing known lines and their descriptions, which should be added to the plot",default=None)
+	parser.add_argument("--knowncombs",nargs="+",help="Combs which are already known")
 	parser.add_argument("--zeroKnownLines",help="Zero out the known lines (requires knownLinesFile)",type='bool',default=False)
 
 	return parser
@@ -59,6 +60,12 @@ def main():
 						else:
 							print(fname)
 							freq,asd=h5.load_spectrum(foldername+fname)
+
+						if "fscan" in fname:
+							yaxlabel="Normalized average power"
+						else:
+							yaxlabel="ASD"
+
 						if args.scalerange!=0:
 							asd=cf.scale_data(asd,args.scalerange)
 
@@ -77,7 +84,7 @@ def main():
 							asd=keepasd[:]
 							freq=freq[(freq>args.fmin)&(freq<args.fmax)]
 						p=figure(webgl=True,tools="pan,box_zoom,wheel_zoom,resize,reset,save",y_axis_type="log",
-							x_axis_label="Frequency",y_axis_label="ASD",x_range=[args.fmin,args.fmax],y_range=[np.percentile(keepasd,0),np.percentile(keepasd,100)],
+							x_axis_label="Frequency",y_axis_label=yaxlabel,x_range=[args.fmin,args.fmax],y_range=[np.percentile(keepasd,0),np.percentile(keepasd,100)],
 							plot_width=1000,plot_height=500)
 						l1=p.line(freq,asd)
 						p.add_tools(HoverTool(renderers=[l1], tooltips='@x{1.1111} Hz',line_policy='nearest'))
@@ -99,11 +106,34 @@ def main():
 							for tagcomb in args.tagcombs:
 								combspacing,comboffset=np.array(tagcomb.split(","),dtype=np.float)
 								comb_inds=cf.predict_comb(freq,comboffset,combspacing)
+								#bg,stddev=cf.get_background(asd,np.arange(0,len(asd),1),30,calcStd=True) # TEMP
+								#mark_inds=cf.mark_comb(asd,comb_inds,args.checkrange,args.consecutive,spectBg=bg,spectStd=stddev,threshold=0) #TEMP
 								mark_inds=cf.mark_comb(asd,comb_inds,args.checkrange,args.consecutive)
 								if len(mark_inds)>0:
 									p.circle(freq[mark_inds],asd[mark_inds],size=10,color=colors[iColors],fill_alpha=.5,legend="Spacing {}, offset {}".format(combspacing,comboffset))
 								iColors+=1
+						if args.knowncombs:
+							knownCombInds=np.zeros(0)
+							knownCombTags=np.zeros(0)
+							for knowncomb in args.knowncombs:
+								combspacing,comboffset=np.array(knowncomb.split(","),dtype=np.float)
+								comb_inds=cf.predict_comb(freq,comboffset,combspacing)
+								mark_inds=cf.mark_comb(asd,comb_inds,args.checkrange,args.consecutive)
+								names=np.repeat(["Comb {},{}".format(combspacing,comboffset)],len(mark_inds))
+								knownCombInds=np.append(knownCombInds,mark_inds)
+								knownCombTags=np.append(knownCombTags,names)
+							knownCombInds=np.hstack(knownCombInds).astype(int)
+							knownCombTags=np.hstack(knownCombTags)
+							knownSource=ColumnDataSource(data=dict(xK=freq[knownCombInds],yK=asd[knownCombInds],descK=knownCombTags))
+							known=p.square('xK','yK',source=knownSource,color='black',size=3)
+							p.add_tools(HoverTool(renderers=[known],tooltips="@xK{1.1111} Hz<br>@descK",line_policy='nearest'))
+						try:
 							save(p)
+						except:
+							print("... Plotting failed; data may not be valid.")
+							if args.verbose==True:
+								print(traceback.format_exc())
+							continue
 						print("... Finished successfully.")
 					else:
 						print("... Output file already exists; skipping.")

findComb.py

+import generalio as g
+import combfinder as cf
+import argparse
+import numpy as np
+
+def makeparser():
+	parser=argparse.ArgumentParser(description="Early test of comb finder.")
+	g.register_custom_types(parser)
+	parser.add_argument("--inputfile",help="File to process")
+	parser.add_argument("--fmin",type=float,default=0,help="Minimum frequency to consider")
+	parser.add_argument("--fmax",type=float,help="Maximim frequency to consider")
+	parser.add_argument("--knownlinesfile",default=None,help="Text file containing known lines with text tag")
+	parser.add_argument("--knowncombs",nargs="+",default=None,help="Combs which are already known")
+	parser.add_argument("--printwithstat",type='bool',default=False,help="Print combs sequentially with the associated strength stat (default false)")
+	parser.add_argument("--printforplotting",type='bool',default=True,help="Print combs in the correct format for plotting (default true)")
+	parser.add_argument("--showbeforezero",type='bool',default=False,help="Print combs before iterative zeroing process")
+	parser.add_argument("--spacingmin",type=float,default=None,help="Minimum comb spacing to accept")
+	parser.add_argument("--spacingmax",type=float,default=None,help="Maximum comb spacing to accept")
+	parser.add_argument("--scalerange",type=int,default=5,help="Number of nearest bins to use for background calculation (default 5)")
+	parser.add_argument("--ntoppeaks",type=int,default=600,help="Top peaks to consider (default 600)")
+	parser.add_argument("--nampneighbors",type=int,default=30,help="Nearest amplitude neighbors to consider (default 30)")
+	parser.add_argument("--nmostcommoncombs",type=int,default=1000,help="Number of most common combs to check (default 300)")
+	parser.add_argument("--consecutive",type=int,default=3,help="Number of consecutive peaks required to define a comb (default 3)")
+	parser.add_argument("--checkrange",type=int,default=2,help="Number of nearest neighbors which must be lower to define a peak (default 2)")
+	parser.add_argument("--threshold",type=float,default=0,help="Number of standard deviations by which a peak must exceed the background (default 0)")
+	parser.add_argument("--verbose",type='bool',default=False)	
+	return parser
+
+def main():
+	''' '''
+	import hdf5io as h5
+
+	parser=makeparser()
+	args=parser.parse_args()
+
+	print("Loading data...")
+	vFreq,vData=h5.load_spectrum(args.inputfile)
+
+	print("Clipping spectrum to selected frequency range...")
+	vFreq,vData=cf.clipspect(vFreq,vData,args.fmin,args.fmax)
+	
+	if args.knownlinesfile:
+		print("Computing indices for known lines...")
+		cat=np.genfromtxt(args.knownlinesfile,delimiter=",",dtype=np.str)
+		knownLines=cat[:,0].astype(np.float)
+		knownLineTags=cat[:,1]
+		knownLineTags=knownLineTags[(knownLines>args.fmin)&(knownLines<args.fmax)]
+		knownLines=knownLines[(knownLines>args.fmin)&(knownLines<args.fmax)]
+		knownLineInds=cf.arrNearest(knownLines,vFreq)
+
+	else:
+		knownLineInds=None
+
+	if args.knowncombs:
+		for knowncomb in args.knowncombs:
+			combspacing,comboffset=np.array(knowncomb.split(","),dtype=np.float)
+			comb_inds=cf.predict_comb(vFreq,comboffset,combspacing)
+			if knownLineInds is None:
+				knownLineInds=comb_inds
+			else:
+				knownLineInds=np.append(knownLineInds,comb_inds)
+
+	knownLineInds=np.hstack(knownLineInds).astype(int)
+	print("Running combs algorithm...")
+	combs=cf.auto_find_comb(vFreq,vData,scaleRange=args.scalerange,spacingmin=args.spacingmin,spacingmax=args.spacingmax,
+		nTopPeaks=args.ntoppeaks,nNearestAmpNeighbors=args.nampneighbors,nMostCommonCombs=args.nmostcommoncombs,
+		nConsecReq=args.consecutive,nCheckAdj=args.checkrange,verbose=args.verbose,showBeforeZeroing=args.showbeforezero,
+		zeroInds=knownLineInds,threshold=args.threshold)
+
+	if args.printwithstat:
+		for i in combs[:]:
+			print("\"{:.6f},{:.6f}\" {:.2f}".format(i['sp'],i['off'],i['str']))
+	if args.printforplotting:
+		l=""
+		for i in combs[:]:
+			l+="\"{:.6f},{:.6f}\" ".format(i['sp'],i['off'])
+		print(l)
+
+
+if __name__=='__main__':
+	main()

hdf5io.py

@@ -2,6 +2,7 @@ from __future__ import division
 import h5py
 import numpy as np
 import combfinder as cf
+import sys
 
 def truncate_file(fname,fmax):
 	freqs,values=load_spectrum(fname)
@@ -140,7 +141,14 @@ def load_spectrum(fname):
 	dset=f['spec_data']
 	values=dset[:]
 	freqs=np.arange(dset.attrs['fmin'],dset.attrs['fmax']+dset.attrs['fstep'],dset.attrs['fstep'])
-	assert len(values)==len(freqs)
+	if not len(values)==len(freqs):
+		print("ERROR in hdfio.load_spectrum: `values` has length {} while `freqs` has length {}".format(len(values),len(freqs)))
+		print(freqs[-1],dset.attrs['fmax'],dset.attrs['fstep'])
+		if len(freqs)-1==len(values) and freqs[-1]>dset.attrs['fstep']:
+			print("Looks like rounding caused a minor error in np.arange . Correcting.")
+			freqs=freqs[:-1]
+		else:
+			sys.exit()
 	f.close()
 	return freqs,values
 

trackComb.py

@@ -23,12 +23,15 @@ def makeparser():
 	parser.add_argument("--consecutive",help="Require that marked peaks be part of a chain of N consecutive comb peaks (default=5)",type=int,default=5)
 	parser.add_argument("--parentfolder",type='folder',help="Parent folder (e.g. /home/aneunzert/public_html/)")
 	parser.add_argument("--parenturl",type=str,help="URL corresponding to parent folder (e.g. ldas-jobs.ligo-wa.caltech.edu/~aneunzert/)")
+	parser.add_argument("--outputfolder",type=str,help="URL corresponding to output folder, if different from parent folder)")
 	parser.add_argument("--overwrite",type='bool',help="Overwrite existing files (default false)",default=False)
 	return parser
 
 def main():
 	''' '''
 	args=makeparser().parse_args()
+	if not args.outputfolder:
+		args.outputfolder=args.parentfolder
 	combstr=np.dtype([('unmarked','>f4'),('marked','>f4'),('pervasive','>f4'),('lowerb','>f4'),('upperb','>f4'),('file','|S10000'),('date','|S10')])
 	chlist = [line[:-1] for line in open(args.chlist,'r')]
 	if args.channels and args.chlist:
@@ -47,7 +50,7 @@ def main():
 
 		for tagcomb in args.tagcomb:
 			combspacing,comboffset=np.array(tagcomb.split(","),dtype=np.float)
-			fnamehtml=args.parentfolder+"combtracking/sp"+g.num_to_pstr(combspacing)+"_off"+g.num_to_pstr(comboffset)+"_"+channel+".html"
+			fnamehtml=args.outputfolder+"combtracking/sp"+g.num_to_pstr(combspacing)+"_off"+g.num_to_pstr(comboffset)+"_"+channel+".html"
 			if os.path.isfile(fnamehtml) and args.overwrite==False:
 				print("File exists, skipping")
 				continue
@@ -77,7 +80,10 @@ def main():
 
 			rejectinds=cf.bootstrap_reject(data['unmarked'],sigma=1000,nmax=5)
 			print("Rejected {} outliers with bad data.".format(len(rejectinds)))
+			if len(rejectinds)<.3*len(data['unmarked']):
 				data=np.delete(data,rejectinds)
+			else:
+				print("Rejected points would comprise more than 30% of data; rejecting nothing.")
 
 			print(combspacing,comboffset)
 			print(np.percentile(data['pervasive']*data['marked']/data['upperb'],99))