23"""Search for long-transient gravitational waves using a piecewise frequency model"""
48start_time = time.time()
51signal.signal(signal.SIGINT, signal.SIG_DFL)
53parser = ap.ArgumentParser()
63 help=
"Specifies the tbank object to use for a search",
71 help=
"The job array number on OzSTAR, to ensure unique random seeds, j < 10,000",
76parser.add_argument(
"--tstart", type=int, help=
"Start time of SFTs", default=
None)
80 help=
"What percentage tref occurs on each segment. E.g. if 1/2, then tref will be halfway through each segment",
86 help=
"Noise level to be added to data",
93 help=
"The number of templates with the greatest fstats to record",
98 help=
"Whether to return the 2F stat averaged over all segments",
103 help=
"Whether to write SFT files or have them stored in memory",
108 help=
"If selected all files will be written to local storage and deleted once job finished except return values of pwf.semifstatcatalogue",
114 help=
"BaseSeed * 10e7 + job number * 10e3 + iterations is the seed for generating random signal params",
120 "--Tsft", type=int, help=
"The length of each SFT to be built", default=10
122parser.add_argument(
"--h0", type=float, help=
"Strain of the injected signal", default=0)
124 "--cosi", type=float, help=
"Default injection cos-iota", default=
None
126parser.add_argument(
"--psi", type=float, help=
"Default injection psi", default=
None)
127parser.add_argument(
"--phi0", type=float, help=
"Default injection phi0", default=
None)
129 "--dt_wf", type=float, help=
"Sampling time for waveform generation", default=10
134 help=
"Right Ascension to sky position of source",
138 "--Delta", type=float, help=
"Declination to sky position of source", default=
None
143 help=
"Set detector being used. Options H1, V1, etc.",
148 "--dfreq", type=float, help=
"Frequency spacing for ComputeFstat()", default=0
153 help=
"sourceDeltaT option for ComputeFstat()",
159 "--s", type=int, help=
"Set the s parameter", nargs=
"?", const=1, default=
None
164 help=
"Set the fmin parameter",
172 help=
"Set the fmax parameter",
180 help=
"Set the nmin parameter",
188 help=
"Set the nmax parameter",
196 help=
"Inertia of the NS. Used to calculate default k values",
204 help=
"Ellipticity of NS. Used to calculate default k values",
212 help=
"Radius of the NS. Used to calculate default k values",
220 help=
"Set the kmin parameter. Overrides Izz, ellip and radius values",
228 help=
"Set the kmax parameter. Overrides Izz, ellip and radius values",
236 help=
"The maximum duration the knot algorithm to calculate up to",
242 "--knots", type=float, help=
"Use user defined knots", nargs=
"+", default=[0]
247 help=
"Set the mismatch parameter",
255 help=
"Use a knotnum number of knots generated by algorithm",
263 help=
"Uses knots generated by the algorithm extending to this duration",
271 help=
"Fractional bounding box padding to use on each dimension",
278 help=
"Integer points to use for padding on each dimension",
286 help=
"If selected, does an iterations number of searches, with signal h0 varying from min_h0 to max_h0. Also includes h0 = 0 for computing 2F* simultaneously",
291 help=
"Does an iterations number of searches with no injected signal. Useful for calculating 2F* without calculating det_prob",
297 help=
"Minimum h0 for calculating detection probabilities, det_prob must be true",
303 help=
"Maximum h0 for calculating detection probabilities, det_prob must be true",
308 help=
"Returns the template count for the parameter space. Does not compute anything else",
313 help=
"If selected with Fstat_mismatch, only returns the lowest mis_match",
318 help=
"If selected, carries out a search on fake data with an injected signal",
323 help=
"If selected, injects a randomly generated signal into the SFTs located at the path set by . Can be selected alongside other modes",
330 help=
"If selected, determines mismatches using the 2F definition, otherwise uses the metric",
334 "--signal_as_template",
335 help=
"If selected, the random injected signal will coincide with a template",
341 help=
"Number of loops for the code to perform. For unique random seeds, iterations < 1000",
347 help=
"Defines frequency bands for search to cover. Frequency bands chosen as fmin=freq_bands[j], fmax=freq_bands[j+1]",
353 help=
"If selected, returns data on the distribution of 2Fs as a histogram",
358 help=
"Use simulated data with a noise level equal to the PSD of data",
363parser.add_argument(
"--outbasedir", type=str, help=
"Output base directory", default=
".")
366 help=
"If true, log to file; if false, log to standard output",
370 "--loglevel", help=
"Level at which to log messages", default=logging.INFO
372parser.add_argument(
"--profile", help=
"Profile the code", action=
"store_true")
376 help=
"Path to noise data.",
382 help=
"Path to psd data.",
388 help=
"Path to sfts.",
396args = parser.parse_args()
397outbasedirectory = args.outbasedir
399ek.setknotarchivepath(outbasedirectory)
402logging_format =
"%(asctime)s : %(levelname)s : %(message)s"
404 logging_filename = os.path.join(outbasedirectory, f
"PiecewiseSearchLog_{j}.log")
406 filename=logging_filename,
409 format=logging_format,
411 logging.info(f
"Logging to {logging_filename}")
413 logging.basicConfig(level=args.loglevel, format=logging_format)
414 logging.info(f
"Logging to standard output")
415logging.captureWarnings(
True)
416if logging.getLogger().isEnabledFor(logging.INFO):
417 lp.globalvar.LatticeTilingProgressLogLevel = lal.LOG_NORMAL
419job_num = os.getenv(
"SLURM_ARRAY_TASK_ID")
420logging.info(
"Job_array_num_is_" +
str(job_num))
424 logging.info(
"Making the profiler")
425 pr = cProfile.Profile()
428logging.debug(
"Doing the arg stuff")
434tbankcode = args.tbankcode
441trefsegfrac = args.trefsegfrac
443tempsperfile = args.tempsperfile
444SFTFiles = args.SFTFiles
445UseLocal = args.UseLocal
446BaseSeed = args.BaseSeed
447iterations = args.iterations
449Random_Seed = BaseSeed * 10000000 + j * 1000
450Random_Seed_Str =
"{:0>9d}".
format(Random_Seed)
459elif args.threshold_2F:
460 mode =
"threshold_2F"
461elif args.template_count:
462 mode =
"template_count"
468logging.info(f
"Search mode: {mode}")
475if tbankcode ==
"GW170817":
476 logging.info(f
"Starting with default parameter space: {tbankcode}")
477 tbank.SetDefaultGW170817()
479elif tbankcode ==
"GW190425":
480 logging.info(f
"Starting with default parameter space: {tbankcode}")
481 tbank.SetDefaultGW190425()
483elif tbankcode ==
"1987A":
484 logging.info(f
"Starting with default parameter space: {tbankcode}")
485 tbank.SetDefault1987A()
487elif tbankcode ==
"CasA":
488 logging.info(f
"Starting with default parameter space: {tbankcode}")
489 tbank.SetDefaultCasA()
491 logging.info(
"No tbank selected! Put an error here!")
495 args.fmin = args.freq_bands[j]
496 args.fmax = args.freq_bands[j + 1]
499tbank.SetTBankParams(args)
505 logging.info(
"Setting knots from command line")
506 bf.knotslist = args.knots
507elif args.knot_alg_dur:
508 logging.info(
"Setting knots from algorithm")
509 tbank.dur = args.knot_alg_dur
510 tbank.SetKnotsByAlg()
512 logging.info(
"Using default knots")
513 bf.knotslist = tbank.knots
515tbank.knots = bf.knotslist
516tbank.dur = tbank.knots[-1] - tbank.knots[0]
519if bf.knotslist[0] == 0:
520 for i, knot
in enumerate(bf.knotslist):
521 bf.knotslist[i] = knot + tbank.tstart
523logging.info(
"Knots: %s",
str(bf.knotslist))
528if args.flags_bbox == []:
529 tbank.flags_bbox = [-1] * tbank.s * len(bf.knotslist)
531if args.flags_intbox == []:
532 tbank.flags_intbox = [-1] * tbank.s * len(bf.knotslist)
538detectors = tbank.detectors
540logging.info(f
"Detectors: {detectors}")
541logging.info(f
"Arg Detectors: {args.detectors}")
542logging.info(f
"TBank Detectors: {tbank.detectors}")
544num_det = len(detectors)
546detector_str_list = lal.CreateStringVector(detectors[0])
549for detector
in detectors[1:]:
550 lal.AppendString2Vector(detector_str_list, detector)
552noise_path = tbank.noise_path
553psd_path = tbank.psd_path
557if args.noise_sqrt_Sh == []
and mode !=
"search":
559 for detector
in detectors:
560 path_to_data = noise_path + detector +
"asd.txt"
563 path_to_data = psd_path + detector +
".txt"
568 if mode ==
"template_count":
571 with open(path_to_data)
as noise_curve:
572 for line
in noise_curve:
574 this_line = line.split()
578 if min_freq_strain == -1
and frequency > args.fmin:
579 min_freq_strain = strain
581 if max_freq_strain == -1
and frequency > args.fmax:
582 max_freq_strain = strain
585 this_det_noise = (min_freq_strain + max_freq_strain) / 2
586 noise_curve_data.append(this_det_noise)
588logging.info(
"Noise curve data for given frequency band: %s",
str(noise_curve_data))
589logging.info(
"User specified noise used (if given): %s",
str(args.noise_sqrt_Sh))
597finputdata = cd.FInput()
598finputdata.Tsft = args.Tsft
599finputdata.tstart = tbank.tstart
600finputdata.trefsegfrac = args.trefsegfrac
601finputdata.h0 = args.h0
602finputdata.cosi = args.cosi
603finputdata.psi = args.psi
604finputdata.phi0 = args.phi0
605finputdata.dt_wf = args.dt_wf
606finputdata.Alpha = tbank.Alpha
607finputdata.Delta = tbank.Delta
608finputdata.detectors = detector_str_list
609finputdata.noise_sqrt_Sh = (
610 args.noise_sqrt_Sh
if args.noise_sqrt_Sh != []
else noise_curve_data
612finputdata.dfreq = args.dfreq
613finputdata.sourceDeltaT = args.sourceDeltaT
614finputdata.inject_data = args.inject_data
622append_string =
"_BS-" +
str(args.BaseSeed) +
"_j-" +
str(j)
625 local_disk_directory = os.getenv(
"JOBFS")
627 os.path.join(local_disk_directory, tbank.toString()) + append_string
630 tbankdirectory = os.path.join(outbasedirectory, tbank.toString()) + append_string
632if not os.path.isdir(tbankdirectory):
633 os.mkdir(tbankdirectory)
635logging.info(f
"TBank directory is {tbankdirectory}")
640minimum_h0 = args.min_h0
641maximum_h0 = args.max_h0
643power_list = np.linspace(np.log10(minimum_h0), np.log10(maximum_h0), iterations)
645h0_list = [0] + [10**i
for i
in power_list]
654for i
in range(iterations + plus_one):
655 logging.info(
"Performing iteration: %s",
str(i))
659 finputdata.h0 = h0_list[i]
661 if args.threshold_2F:
664 logging.debug(
"Signal strength: %s",
str(finputdata.h0))
667 this_seed = Random_Seed + i
674 iteration_string =
"_i-" +
str(i)
676 tempsfilename =
"Temps_For_" + tbank.toString() + iteration_string + sumstr +
".txt"
678 "Mismatch_Temps_For_" + tbank.toString() + iteration_string + sumstr +
".txt"
682 if mode !=
"search" or args.inject_data:
684 seed = random.seed(this_seed)
686 finputdata.phi0 = random.uniform(0, 2 * np.pi)
687 finputdata.psi = random.uniform(0, 2 * np.pi)
688 finputdata.cosi = random.uniform(-1, 1)
690 logging.info(
"phi0: %s",
str(finputdata.phi0))
691 logging.info(
"psi: %s",
str(finputdata.psi))
692 logging.info(
"cosi: %s",
str(finputdata.cosi))
703 signalparams, SFTdirectory = pwsim.buildSFTs(
707 trefsegfrac=trefsegfrac,
708 parentdirectory=tbankdirectory,
712 elif mode !=
"search" or (mode ==
"search" and args.inject_data):
713 logging.debug(
"Creating tiling lattice")
714 tiling = lp.CreateLatticeTiling(tbank.s * len(bf.knotslist))
716 logging.debug(
"Building metric")
717 metric = scmm.PreCompMetric(tbank.s)
719 logging.debug(
"Setting Bounds")
721 tbe.setbounds(tiling, tbank)
723 logging.debug(
"Setting tiling lattice and metric")
724 lp.SetTilingLatticeAndMetric(
725 tiling, lp.TILING_LATTICE_ANSTAR, metric, tbank.maxmismatch
728 logging.debug(
"Creating random signal params")
730 randparams = lal.CreateRandomParams(this_seed)
731 signalparams = lal.gsl_matrix(tbank.s * len(bf.knotslist), 1)
732 lp.RandomLatticeTilingPoints(tiling, 0, randparams, signalparams)
738 if args.signal_as_template:
739 locator = lp.CreateLatticeTilingLocator(tiling)
741 signal_params_vec = lal.gsl_vector(tbank.s * len(bf.knotslist))
744 for i, elem
in enumerate(signalparams.data[0]):
745 signal_params_vec.data[i] = elem
747 nearest_point = lal.gsl_vector(tbank.s * len(bf.knotslist))
748 nearest_index = lal.CreateUINT8Vector(tbank.s * len(bf.knotslist))
750 lp.NearestLatticeTilingPoint(
751 locator, signal_params_vec, nearest_point, nearest_index
754 signalparams = np.transpose(nearest_point.data)
757 signalparams = np.transpose(signalparams.data)[0]
759 logging.info(
"Random Signal Params are: %s",
str(signalparams))
767 SFTdirectory = os.path.join(
769 f
"SFTs_h0-{finputdata.h0:.2e}_f0-{f0:.3f}_f1-{f1:.2e}_f2-{f2:.2e}_dur-{tbank.dur}_tstart-{tbank.tstart}/",
772 if finputdata.inject_data:
774 SFTdirectory = tbank.sft_path
776 elif mode ==
"search" and not args.inject_data:
778 SFTdirectory = tbank.sft_path
780 signalparams = [-1] * (tbank.s * len(bf.knotslist))
782 logging.info(
"SFT's for this iteration have path: %s", SFTdirectory)
786 with open(os.path.join(tbankdirectory, tempsfilename),
"w")
as reader:
788 commentline =
"{:20s}".
format(
"#FStat") +
" "
790 for detector_string
in detectors:
791 commentline +=
"{:20s}".
format(detector_string) +
" "
793 commentline +=
str(
"{:20s}".
format(
"Mismatch")) +
" "
797 for i
in range(len(signalparams)):
802 derivnum = i % tbank.s
807 columntitle =
"f_" +
str(knotnum) +
"_" +
str(derivnum)
808 commentline +=
" " +
"{:^20s}".
format(columntitle)
810 reader.write(commentline +
"\n")
815 with open(os.path.join(tbankdirectory, mismatchtempname),
"w")
as reader:
818 "{:20s}".
format(
"#Mismatch")
824 for detector_string
in detectors:
825 commentline +=
"{:20s}".
format(detector_string) +
" "
829 for i
in range(len(signalparams)):
831 derivnum = i % tbank.s
836 columntitle =
"f_" +
str(knotnum) +
"_" +
str(derivnum)
837 commentline +=
" " +
"{:^20s}".
format(columntitle)
839 reader.write(commentline +
"\n")
842 logging.debug(
"Setting Antenna Pattern")
843 lp.SetAntennaPatternMaxCond(10**5)
847 SFTfmin = gom.gte(tbank.dur, tbank.fmin, tbank.nmax, tbank.kmax) - SFTfpad - 10
848 SFTfmax = gom.gte(0, tbank.fmax, tbank.nmin, tbank.kmin) + SFTfpad
850 logging.info(f
"SFTfmin/fmax: [{SFTfmin}, {SFTfmax}]")
851 logging.info(f
"SFT files is: {SFTFiles}")
853 Fstat_mismatch = args.Fstat_mismatch
854 logging.info(f
"2F mismatch is: {Fstat_mismatch}")
857 if mode ==
"search" and not args.fstat_hist:
859 lowest_mismatch_metric,
860 lowest_mismatch_Fstat,
863 ) = pwf.pw_fstat_search_catalogue(
869 tbankdirectory=tbankdirectory,
870 SFTdirectory=SFTdirectory,
871 trefsegfrac=trefsegfrac,
873 tempsperfile=tempsperfile,
877 lowest_mismatch_metric,
878 lowest_mismatch_Fstat,
881 ) = pwf.semifstatcatalogue(
888 out_directory=tbankdirectory,
889 SFTdirectory=SFTdirectory,
890 trefsegfrac=trefsegfrac,
893 tempsperfile=tempsperfile,
894 Fstat_mismatch=args.Fstat_mismatch,
895 fstat_hist=args.fstat_hist,
903 if template_count != 0:
904 temp_count_directory =
"Template_Count_And_Timing_For_" + tbank.toString()
905 temp_count_file =
"Template_Count_And_Timing_" + Random_Seed_Str +
".txt"
907 if not os.path.isdir(temp_count_directory):
908 os.mkdir(temp_count_directory)
911 time_taken = time.time() - start_time
913 with open(temp_count_directory +
"/" + temp_count_file,
"a+")
as reader:
914 line =
str(time_taken) +
" " +
str(template_count) +
"\n"
924 Two_F_directory =
"Largest_2Fs_For_" + tbank.toString()
925 Two_F_file =
"Largest_2F_" + Random_Seed_Str +
".txt"
927 if not os.path.isdir(Two_F_directory):
928 os.mkdir(Two_F_directory)
930 with open(Two_F_directory +
"/" + Two_F_file,
"a+")
as reader:
931 line =
str(finputdata.h0) +
" " +
str(fstat_max) +
"\n"
936 if lowest_mismatch_metric != 0:
937 mismatch_directory =
"Mismatches_For_" + tbank.toString()
938 mismatch_file =
"Mismatch_" + Random_Seed_Str +
".txt"
940 if not os.path.isdir(mismatch_directory):
941 os.mkdir(mismatch_directory)
943 with open(mismatch_directory +
"/" + mismatch_file,
"a+")
as reader:
945 str(lowest_mismatch_metric)
947 +
str(lowest_mismatch_Fstat)
958 ps = pstats.Stats(pr, stream=s).sort_stats(
"cumtime")
961 os.path.join(basedirectory, f
"PiecewiseSearchProfile_{j}.txt"),
"w+"
963 f.write(s.getvalue())