Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-3a66518
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
lalpulsar_PiecewiseSearch.py
Go to the documentation of this file.
1##python
2# Copyright (C) 2019--2023 Benjamin Grace
3#
4# This program is free software; you can redistribute it and/or modify it
5# under the terms of the GNU General Public License as published by the
6# Free Software Foundation; either version 2 of the License, or (at your
7# option) any later version.
8#
9# This program is distributed in the hope that it will be useful, but
10# WITHOUT ANY WARRANTY; without even the implied warranty of
11# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
12# Public License for more details.
13#
14# You should have received a copy of the GNU General Public License along
15# with this program; if not, write to the Free Software Foundation, Inc.,
16# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17
18## \defgroup lalpulsar_bin_PiecewiseSearch Piecewise Search Applications
19## \ingroup lalpulsar_bin_Apps
20
21## \file
22## \ingroup lalpulsar_bin_PiecewiseSearch
23"""Search for long-transient gravitational waves using a piecewise frequency model"""
24
25import argparse as ap
26import cProfile
27import io
28import logging
29import os
30import pstats
31import random
32import signal
33import time
34
35import numpy as np
36
37import lal
38import lalpulsar as lp
47
48start_time = time.time()
49
50# Allow liblalpulsar C code to be interrupted with Ctrl+C
51signal.signal(signal.SIGINT, signal.SIG_DFL)
52
53parser = ap.ArgumentParser()
54
55# Thoughts on this module
56# It is beginning to get quite convoluted. Is it worth writing a new python module where a lot of this set up is written as functions? It may make this module more succinct
57# and readable
58
59# Required parameters
60parser.add_argument(
61 "--tbankcode",
62 type=str,
63 help="Specifies the tbank object to use for a search",
64 required=True,
65)
66
67# Random seed, typically chosen as the job array number of OzSTAR. To ensure unique random seeds, j < 10,000
68parser.add_argument(
69 "--j",
70 type=int,
71 help="The job array number on OzSTAR, to ensure unique random seeds, j < 10,000",
72 default=0,
73)
74
75# Parameters used for specifying signals to inject in SFTs as well as other parameters for F-Statistic calculations
76parser.add_argument("--tstart", type=int, help="Start time of SFTs", default=None)
77parser.add_argument(
78 "--trefsegfrac",
79 type=float,
80 help="What percentage tref occurs on each segment. E.g. if 1/2, then tref will be halfway through each segment",
81 default=0.0,
82)
83parser.add_argument(
84 "--noise_sqrt_Sh",
85 type=float,
86 help="Noise level to be added to data",
87 nargs="+",
88 default=[],
89)
90parser.add_argument(
91 "--tempsperfile",
92 type=int,
93 help="The number of templates with the greatest fstats to record",
94 default=1000,
95)
96parser.add_argument(
97 "--rtnsum",
98 help="Whether to return the 2F stat averaged over all segments",
99 action="store_true",
100)
101parser.add_argument(
102 "--SFTFiles",
103 help="Whether to write SFT files or have them stored in memory",
104 action="store_true",
105)
106parser.add_argument(
107 "--UseLocal",
108 help="If selected all files will be written to local storage and deleted once job finished except return values of pwf.semifstatcatalogue",
109 action="store_true",
110)
111parser.add_argument(
112 "--BaseSeed",
113 type=int,
114 help="BaseSeed * 10e7 + job number * 10e3 + iterations is the seed for generating random signal params",
115 default=0,
116)
117
118# Parameters to be used when creating the necessary elements to use lp.ComputeFStat
119parser.add_argument(
120 "--Tsft", type=int, help="The length of each SFT to be built", default=10
121)
122parser.add_argument("--h0", type=float, help="Strain of the injected signal", default=0)
123parser.add_argument(
124 "--cosi", type=float, help="Default injection cos-iota", default=None
125)
126parser.add_argument("--psi", type=float, help="Default injection psi", default=None)
127parser.add_argument("--phi0", type=float, help="Default injection phi0", default=None)
128parser.add_argument(
129 "--dt_wf", type=float, help="Sampling time for waveform generation", default=10
130)
131parser.add_argument(
132 "--Alpha",
133 type=float,
134 help="Right Ascension to sky position of source",
135 default=None,
136)
137parser.add_argument(
138 "--Delta", type=float, help="Declination to sky position of source", default=None
139)
140parser.add_argument(
141 "--detectors",
142 type=str,
143 help="Set detector being used. Options H1, V1, etc.",
144 nargs="+",
145 default=[],
146)
147parser.add_argument(
148 "--dfreq", type=float, help="Frequency spacing for ComputeFstat()", default=0
149)
150parser.add_argument(
151 "--sourceDeltaT",
152 type=float,
153 help="sourceDeltaT option for ComputeFstat()",
154 default=0.1,
155)
156
157# Parameters used in defining the parameter space
158parser.add_argument(
159 "--s", type=int, help="Set the s parameter", nargs="?", const=1, default=None
160)
161parser.add_argument(
162 "--fmin",
163 type=float,
164 help="Set the fmin parameter",
165 nargs="?",
166 const=1,
167 default=None,
168)
169parser.add_argument(
170 "--fmax",
171 type=float,
172 help="Set the fmax parameter",
173 nargs="?",
174 const=1,
175 default=None,
176)
177parser.add_argument(
178 "--nmin",
179 type=float,
180 help="Set the nmin parameter",
181 nargs="?",
182 const=1,
183 default=None,
184)
185parser.add_argument(
186 "--nmax",
187 type=float,
188 help="Set the nmax parameter",
189 nargs="?",
190 const=1,
191 default=None,
192)
193parser.add_argument(
194 "--Izz",
195 type=float,
196 help="Inertia of the NS. Used to calculate default k values",
197 nargs="?",
198 const=1,
199 default=None,
200)
201parser.add_argument(
202 "--ellip",
203 type=float,
204 help="Ellipticity of NS. Used to calculate default k values",
205 nargs="?",
206 const=1,
207 default=None,
208)
209parser.add_argument(
210 "--radius",
211 type=float,
212 help="Radius of the NS. Used to calculate default k values",
213 nargs="?",
214 const=1,
215 default=None,
216)
217parser.add_argument(
218 "--kmin",
219 type=float,
220 help="Set the kmin parameter. Overrides Izz, ellip and radius values",
221 nargs="?",
222 const=1,
223 default=None,
224)
225parser.add_argument(
226 "--kmax",
227 type=float,
228 help="Set the kmax parameter. Overrides Izz, ellip and radius values",
229 nargs="?",
230 const=1,
231 default=None,
232)
233parser.add_argument(
234 "--dur",
235 type=float,
236 help="The maximum duration the knot algorithm to calculate up to",
237 nargs="?",
238 const=1,
239 default=None,
240)
241parser.add_argument(
242 "--knots", type=float, help="Use user defined knots", nargs="+", default=[0]
243)
244parser.add_argument(
245 "--maxmismatch",
246 type=float,
247 help="Set the mismatch parameter",
248 nargs="?",
249 const=1,
250 default=None,
251)
252parser.add_argument(
253 "--knotnum",
254 type=int,
255 help="Use a knotnum number of knots generated by algorithm",
256 nargs="?",
257 const=1,
258 default=None,
259)
260parser.add_argument(
261 "--knot_alg_dur",
262 type=float,
263 help="Uses knots generated by the algorithm extending to this duration",
264 nargs="?",
265 const=1,
266 default=None,
267)
268parser.add_argument(
269 "--flags_bbox",
270 type=float,
271 help="Fractional bounding box padding to use on each dimension",
272 nargs="+",
273 default=[],
274)
275parser.add_argument(
276 "--flags_intbox",
277 type=int,
278 help="Integer points to use for padding on each dimension",
279 nargs="+",
280 default=[],
281)
282
283# Different modes for the code to run
284parser.add_argument(
285 "--det_prob",
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",
287 action="store_true",
288)
289parser.add_argument(
290 "--threshold_2F",
291 help="Does an iterations number of searches with no injected signal. Useful for calculating 2F* without calculating det_prob",
292 action="store_true",
293)
294parser.add_argument(
295 "--min_h0",
296 type=float,
297 help="Minimum h0 for calculating detection probabilities, det_prob must be true",
298 default=1e-30,
299)
300parser.add_argument(
301 "--max_h0",
302 type=float,
303 help="Maximum h0 for calculating detection probabilities, det_prob must be true",
304 default=1e-10,
305)
306parser.add_argument(
307 "--template_count",
308 help="Returns the template count for the parameter space. Does not compute anything else",
309 action="store_true",
310)
311parser.add_argument(
312 "--mis_hist",
313 help="If selected with Fstat_mismatch, only returns the lowest mis_match",
314 action="store_true",
315)
316parser.add_argument(
317 "--fake_data",
318 help="If selected, carries out a search on fake data with an injected signal",
319 action="store_true",
320)
321parser.add_argument(
322 "--inject_data",
323 help="If selected, injects a randomly generated signal into the SFTs located at the path set by . Can be selected alongside other modes",
324 action="store_true",
325)
326
327# Other options available
328parser.add_argument(
329 "--Fstat_mismatch",
330 help="If selected, determines mismatches using the 2F definition, otherwise uses the metric",
331 action="store_true",
332)
333parser.add_argument(
334 "--signal_as_template",
335 help="If selected, the random injected signal will coincide with a template",
336 action="store_true",
337)
338parser.add_argument(
339 "--iterations",
340 type=int,
341 help="Number of loops for the code to perform. For unique random seeds, iterations < 1000",
342 default=1,
343)
344parser.add_argument(
345 "--freq_bands",
346 type=float,
347 help="Defines frequency bands for search to cover. Frequency bands chosen as fmin=freq_bands[j], fmax=freq_bands[j+1]",
348 nargs="+",
349 default=[],
350)
351parser.add_argument(
352 "--fstat_hist",
353 help="If selected, returns data on the distribution of 2Fs as a histogram",
354 action="store_true",
355)
356parser.add_argument(
357 "--data_psd",
358 help="Use simulated data with a noise level equal to the PSD of data",
359 action="store_true",
360)
361
362# Parameters controlling outputs
363parser.add_argument("--outbasedir", type=str, help="Output base directory", default=".")
364parser.add_argument(
365 "--logfile",
366 help="If true, log to file; if false, log to standard output",
367 action="store_true",
368)
369parser.add_argument(
370 "--loglevel", help="Level at which to log messages", default=logging.INFO
371)
372parser.add_argument("--profile", help="Profile the code", action="store_true")
373parser.add_argument(
374 "--noise_path",
375 type=str,
376 help="Path to noise data.",
377 default=None,
378)
379parser.add_argument(
380 "--psd_path",
381 type=str,
382 help="Path to psd data.",
383 default=None,
384)
385parser.add_argument(
386 "--sft_path",
387 type=str,
388 help="Path to sfts.",
389 default=None,
390)
391
392# -------------------------------------------------------------------------------------------------------------------------------------------------------
393# Logging, profiling and setting the directory for the knots list
394# -------------------------------------------------------------------------------------------------------------------------------------------------------
395
396args = parser.parse_args()
397outbasedirectory = args.outbasedir
398
399ek.setknotarchivepath(outbasedirectory)
400
401# For logging
402logging_format = "%(asctime)s : %(levelname)s : %(message)s"
403if args.logfile:
404 logging_filename = os.path.join(outbasedirectory, f"PiecewiseSearchLog_{j}.log")
405 logging.basicConfig(
406 filename=logging_filename,
407 filemode="w",
408 level=args.loglevel,
409 format=logging_format,
410 )
411 logging.info(f"Logging to {logging_filename}")
412else:
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
418
419job_num = os.getenv("SLURM_ARRAY_TASK_ID")
420logging.info("Job_array_num_is_" + str(job_num))
421
422# For profiling
423if args.profile:
424 logging.info("Making the profiler")
425 pr = cProfile.Profile()
426 pr.enable()
427
428logging.debug("Doing the arg stuff")
429
430# -------------------------------------------------------------------------------------------------------------------------------------------------------
431# Required arguments
432# -------------------------------------------------------------------------------------------------------------------------------------------------------
433
434tbankcode = args.tbankcode
435
436# -------------------------------------------------------------------------------------------------------------------------------------------------------
437# Optional arguments and Random Seed
438# -------------------------------------------------------------------------------------------------------------------------------------------------------
439
440j = args.j
441trefsegfrac = args.trefsegfrac
442rtnsum = args.rtnsum
443tempsperfile = args.tempsperfile
444SFTFiles = args.SFTFiles
445UseLocal = args.UseLocal
446BaseSeed = args.BaseSeed
447iterations = args.iterations
448
449Random_Seed = BaseSeed * 10000000 + j * 1000
450Random_Seed_Str = "{:0>9d}".format(Random_Seed)
451
452# -------------------------------------------------------------------------------------------------------------------------------------------------------
453# Setting the mode the code is run in, influencing its output
454# -------------------------------------------------------------------------------------------------------------------------------------------------------
455mode = "search"
456
457if args.det_prob:
458 mode = "det_prob"
459elif args.threshold_2F:
460 mode = "threshold_2F"
461elif args.template_count:
462 mode = "template_count"
463elif args.mis_hist:
464 mode = "mis_hist"
465elif args.fake_data:
466 mode = "fake_data"
467
468logging.info(f"Search mode: {mode}")
469
470# -------------------------------------------------------------------------------------------------------------------------------------------------------
471# Setting the template bank object and changing user specified parameters
472# -------------------------------------------------------------------------------------------------------------------------------------------------------
473tbank = cd.TBank()
474
475if tbankcode == "GW170817":
476 logging.info(f"Starting with default parameter space: {tbankcode}")
477 tbank.SetDefaultGW170817()
478
479elif tbankcode == "GW190425":
480 logging.info(f"Starting with default parameter space: {tbankcode}")
481 tbank.SetDefaultGW190425()
482
483elif tbankcode == "1987A":
484 logging.info(f"Starting with default parameter space: {tbankcode}")
485 tbank.SetDefault1987A()
486
487elif tbankcode == "CasA":
488 logging.info(f"Starting with default parameter space: {tbankcode}")
489 tbank.SetDefaultCasA()
490else:
491 logging.info("No tbank selected! Put an error here!")
492
493# If using custom frequency bands
494if args.freq_bands:
495 args.fmin = args.freq_bands[j]
496 args.fmax = args.freq_bands[j + 1]
497
498# Override any tbank arguments set on the command line
499tbank.SetTBankParams(args)
500
501# -------------------------------------------------------------------------------------------------------------------------------------------------------
502# Checking for user defined knots, or if knots need to be recalculated based on any changed default parameters
503# -------------------------------------------------------------------------------------------------------------------------------------------------------
504if args.knots != [0]:
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()
511else:
512 logging.info("Using default knots")
513 bf.knotslist = tbank.knots
514
515tbank.knots = bf.knotslist
516tbank.dur = tbank.knots[-1] - tbank.knots[0]
517
518# Adjusting knot start time and duration
519if bf.knotslist[0] == 0:
520 for i, knot in enumerate(bf.knotslist):
521 bf.knotslist[i] = knot + tbank.tstart
522
523logging.info("Knots: %s", str(bf.knotslist))
524
525# -------------------------------------------------------------------------------------------------------------------------------------------------------
526# Padding flag customisation
527# -------------------------------------------------------------------------------------------------------------------------------------------------------
528if args.flags_bbox == []:
529 tbank.flags_bbox = [-1] * tbank.s * len(bf.knotslist)
530
531if args.flags_intbox == []:
532 tbank.flags_intbox = [-1] * tbank.s * len(bf.knotslist)
533
534# -------------------------------------------------------------------------------------------------------------------------------------------------------
535# Setting detector specifics for noise curve data
536# -------------------------------------------------------------------------------------------------------------------------------------------------------
537
538detectors = tbank.detectors
539
540logging.info(f"Detectors: {detectors}")
541logging.info(f"Arg Detectors: {args.detectors}")
542logging.info(f"TBank Detectors: {tbank.detectors}")
543
544num_det = len(detectors)
545
546detector_str_list = lal.CreateStringVector(detectors[0])
547
548# I'm sure there's a better way to initialise the detector LALStringVector
549for detector in detectors[1:]:
550 lal.AppendString2Vector(detector_str_list, detector)
551
552noise_path = tbank.noise_path
553psd_path = tbank.psd_path
554
555noise_curve_data = []
556
557if args.noise_sqrt_Sh == [] and mode != "search":
558
559 for detector in detectors:
560 path_to_data = noise_path + detector + "asd.txt"
561
562 if args.data_psd:
563 path_to_data = psd_path + detector + ".txt"
564
565 min_freq_strain = -1
566 max_freq_strain = -1
567
568 if mode == "template_count":
569 break
570
571 with open(path_to_data) as noise_curve:
572 for line in noise_curve:
573
574 this_line = line.split()
575 frequency = float(this_line[0])
576 strain = float(this_line[1])
577
578 if min_freq_strain == -1 and frequency > args.fmin:
579 min_freq_strain = strain
580
581 if max_freq_strain == -1 and frequency > args.fmax:
582 max_freq_strain = strain
583 break
584
585 this_det_noise = (min_freq_strain + max_freq_strain) / 2
586 noise_curve_data.append(this_det_noise)
587
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))
590
591# -------------------------------------------------------------------------------------------------------------------------------------------------------
592# Setting the FInput object
593# -------------------------------------------------------------------------------------------------------------------------------------------------------
594
595# Note that the parameters cosi, psi and phi0 are randomised in a for loop below.
596
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
611)
612finputdata.dfreq = args.dfreq
613finputdata.sourceDeltaT = args.sourceDeltaT
614finputdata.inject_data = args.inject_data
615
616# -------------------------------------------------------------------------------------------------------------------------------------------------------
617# Building the name of the directory where the SFTs will be saved. Directory name is based off the tbank object. In this way, any searches completed
618# using the same tbank will have all SFTs in the same directory for neatness. Within the tbankdirectory, sub folders are used to contain the SFTs
619# and results of each different search. If UseLocal, then all files are written to OzSTAR local disks, which are deleted when jobs are finished.
620# -------------------------------------------------------------------------------------------------------------------------------------------------------
621
622append_string = "_BS-" + str(args.BaseSeed) + "_j-" + str(j)
623
624if UseLocal:
625 local_disk_directory = os.getenv("JOBFS")
626 tbankdirectory = (
627 os.path.join(local_disk_directory, tbank.toString()) + append_string
628 )
629else:
630 tbankdirectory = os.path.join(outbasedirectory, tbank.toString()) + append_string
631
632if not os.path.isdir(tbankdirectory):
633 os.mkdir(tbankdirectory)
634
635logging.info(f"TBank directory is {tbankdirectory}")
636
637# -------------------------------------------------------------------------------------------------------------------------------------------------------
638# Building strain list for injected signals, only necessary if --det_prob option is chosen
639# -------------------------------------------------------------------------------------------------------------------------------------------------------
640minimum_h0 = args.min_h0
641maximum_h0 = args.max_h0
642
643power_list = np.linspace(np.log10(minimum_h0), np.log10(maximum_h0), iterations)
644
645h0_list = [0] + [10**i for i in power_list]
646
647plus_one = 0
648if args.det_prob:
649 plus_one = 1
650
651# -------------------------------------------------------------------------------------------------------------------------------------------------------
652# Begin search set up. Some options may benefit from running through multiple independent searches, hence all search setup is done in a for loop
653# -------------------------------------------------------------------------------------------------------------------------------------------------------
654for i in range(iterations + plus_one):
655 logging.info("Performing iteration: %s", str(i))
656
657 # Change strain
658 if args.det_prob:
659 finputdata.h0 = h0_list[i]
660
661 if args.threshold_2F:
662 finputdata.h0 = 0
663
664 logging.debug("Signal strength: %s", str(finputdata.h0))
665
666 # Random seed for this iteration
667 this_seed = Random_Seed + i
668
669 # The name of the file that 2F and templates are stored to is only changed by whether the 2F values are averaged or only summed
670 sumstr = ""
671 if rtnsum:
672 sumstr = "_Sum"
673
674 iteration_string = "_i-" + str(i)
675
676 tempsfilename = "Temps_For_" + tbank.toString() + iteration_string + sumstr + ".txt"
677 mismatchtempname = (
678 "Mismatch_Temps_For_" + tbank.toString() + iteration_string + sumstr + ".txt"
679 )
680
681 # Randomise the cosi, phi0 and psi parameters if we are not doing a targetted search
682 if mode != "search" or args.inject_data:
683
684 seed = random.seed(this_seed)
685
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)
689
690 logging.info("phi0: %s", str(finputdata.phi0))
691 logging.info("psi: %s", str(finputdata.psi))
692 logging.info("cosi: %s", str(finputdata.cosi))
693
694 # This chain of else-if statements determines the path to the SFTs being used and the signal parameters that will be/have been injected into them. The first option
695 # is if we are generating our own fake data with gaussian noise by building SFT files. Using the SFT files is advantageous as we can inject any type of signal we wish,
696 # depending on what signal model is defined in pwsim. Signal parameters for the first option are randomly generated inside pwsim. The second option randomly generates
697 # signal parameters inside the parameter space. This option is used if we do not want to write SFT files, but rather store SFTs in memory. This option should be used for
698 # all configurations which are not a search on experimental data OR we are doing a search with experimental data but we want to inject a signal into that data. If the case
699 # is the latter, the random signal parameters are generated, but an enclosed if statement resets the SFT directory to the appropriate path for the experimental data.
700 # The final option is for if we are doing a search, in which case we use the SFT directory to the data SFTs and set the injected signal to all zeros.
701
702 if SFTFiles:
703 signalparams, SFTdirectory = pwsim.buildSFTs(
704 tbank.dur,
705 tbank,
706 finputdata,
707 trefsegfrac=trefsegfrac,
708 parentdirectory=tbankdirectory,
709 rand_seed=this_seed,
710 )
711
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))
715
716 logging.debug("Building metric")
717 metric = scmm.PreCompMetric(tbank.s)
718
719 logging.debug("Setting Bounds")
720
721 tbe.setbounds(tiling, tbank)
722
723 logging.debug("Setting tiling lattice and metric")
724 lp.SetTilingLatticeAndMetric(
725 tiling, lp.TILING_LATTICE_ANSTAR, metric, tbank.maxmismatch
726 )
727
728 logging.debug("Creating random signal params")
729
730 randparams = lal.CreateRandomParams(this_seed)
731 signalparams = lal.gsl_matrix(tbank.s * len(bf.knotslist), 1)
732 lp.RandomLatticeTilingPoints(tiling, 0, randparams, signalparams)
733
734 # -------------------------------------------------------------------------------------------------------------------------------------------------------
735 # If we want a signal to coincide exactly with a template from the template bank. A random point in the parameter space is generated, then the nearest
736 # template found to that point is chosen to be the injected signal
737 # -------------------------------------------------------------------------------------------------------------------------------------------------------
738 if args.signal_as_template:
739 locator = lp.CreateLatticeTilingLocator(tiling)
740
741 signal_params_vec = lal.gsl_vector(tbank.s * len(bf.knotslist))
742
743 # Done as a for loop because setting signal_params_vec.data = signalparams.data[0] doesn't seem to work
744 for i, elem in enumerate(signalparams.data[0]):
745 signal_params_vec.data[i] = elem
746
747 nearest_point = lal.gsl_vector(tbank.s * len(bf.knotslist))
748 nearest_index = lal.CreateUINT8Vector(tbank.s * len(bf.knotslist))
749
750 lp.NearestLatticeTilingPoint(
751 locator, signal_params_vec, nearest_point, nearest_index
752 )
753
754 signalparams = np.transpose(nearest_point.data)
755
756 else:
757 signalparams = np.transpose(signalparams.data)[0]
758
759 logging.info("Random Signal Params are: %s", str(signalparams))
760
761 # Details for the SFT directory
762
763 f0 = signalparams[0]
764 f1 = signalparams[1]
765 f2 = signalparams[2]
766
767 SFTdirectory = os.path.join(
768 tbankdirectory,
769 f"SFTs_h0-{finputdata.h0:.2e}_f0-{f0:.3f}_f1-{f1:.2e}_f2-{f2:.2e}_dur-{tbank.dur}_tstart-{tbank.tstart}/",
770 )
771
772 if finputdata.inject_data:
773
774 SFTdirectory = tbank.sft_path
775
776 elif mode == "search" and not args.inject_data:
777
778 SFTdirectory = tbank.sft_path
779
780 signalparams = [-1] * (tbank.s * len(bf.knotslist))
781
782 logging.info("SFT's for this iteration have path: %s", SFTdirectory)
783
784 # Create the text file where 2F and template data will be stored. The first line is added to the file which contains the column titles for all data. We also write
785 # the signal parameters as the first data entry to this file, however we skip that step if we are carrying out a search
786 with open(os.path.join(tbankdirectory, tempsfilename), "w") as reader:
787
788 commentline = "{:20s}".format("#FStat") + " "
789
790 for detector_string in detectors:
791 commentline += "{:20s}".format(detector_string) + " "
792
793 commentline += str("{:20s}".format("Mismatch")) + " "
794
795 knotnum = -1
796
797 for i in range(len(signalparams)):
798
799 if mode == "search":
800 break
801
802 derivnum = i % tbank.s
803
804 if derivnum == 0:
805 knotnum += 1
806
807 columntitle = "f_" + str(knotnum) + "_" + str(derivnum)
808 commentline += " " + "{:^20s}".format(columntitle)
809
810 reader.write(commentline + "\n")
811
812 # Create the text file where 2F and template data will be stored but ordered by mismatch. The first line is added which contains the column titles for all data. We skip
813 # this if we are doing a search as we can't calculate the mismatch
814 if mode != "search":
815 with open(os.path.join(tbankdirectory, mismatchtempname), "w") as reader:
816
817 commentline = (
818 "{:20s}".format("#Mismatch")
819 + " "
820 + str("{:20s}".format("FStat"))
821 + " "
822 )
823
824 for detector_string in detectors:
825 commentline += "{:20s}".format(detector_string) + " "
826
827 knotnum = -1
828
829 for i in range(len(signalparams)):
830
831 derivnum = i % tbank.s
832
833 if derivnum == 0:
834 knotnum += 1
835
836 columntitle = "f_" + str(knotnum) + "_" + str(derivnum)
837 commentline += " " + "{:^20s}".format(columntitle)
838
839 reader.write(commentline + "\n")
840
841 # Set the maximum condition number for the antenna pattern above its default (get a warning using the piecewise model if not manually changed to be larger)
842 logging.debug("Setting Antenna Pattern")
843 lp.SetAntennaPatternMaxCond(10**5)
844
845 # The minimum and maximum frequencies needed to load in from SFTs to cover all frequencies any template may cover
846 SFTfpad = 10 # 1800 / tbank.Tsft + tbank.dur / 86400 + 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
849
850 logging.info(f"SFTfmin/fmax: [{SFTfmin}, {SFTfmax}]")
851 logging.info(f"SFT files is: {SFTFiles}")
852
853 Fstat_mismatch = args.Fstat_mismatch
854 logging.info(f"2F mismatch is: {Fstat_mismatch}")
855
856 # Output from the search, which is carried out in pw_fstat.py. Output changes depending on run mode. All search results are written to files within this method
857 if mode == "search" and not args.fstat_hist:
858 (
859 lowest_mismatch_metric,
860 lowest_mismatch_Fstat,
861 fstat_max,
862 template_count,
863 ) = pwf.pw_fstat_search_catalogue(
864 SFTfmin,
865 SFTfmax,
866 tbank,
867 finputdata,
868 tempsfilename,
869 tbankdirectory=tbankdirectory,
870 SFTdirectory=SFTdirectory,
871 trefsegfrac=trefsegfrac,
872 rtnsum=rtnsum,
873 tempsperfile=tempsperfile,
874 )
875 else:
876 (
877 lowest_mismatch_metric,
878 lowest_mismatch_Fstat,
879 fstat_max,
880 template_count,
881 ) = pwf.semifstatcatalogue(
882 SFTfmin,
883 SFTfmax,
884 tbank,
885 finputdata,
886 signalparams,
887 tempsfilename,
888 out_directory=tbankdirectory,
889 SFTdirectory=SFTdirectory,
890 trefsegfrac=trefsegfrac,
891 rtnsum=rtnsum,
892 SFTFiles=SFTFiles,
893 tempsperfile=tempsperfile,
894 Fstat_mismatch=args.Fstat_mismatch,
895 fstat_hist=args.fstat_hist,
896 mode=mode,
897 )
898
899 # For each non-zero return value above, we create a directory and write the data into a file based on Random_Seed (not the same random seed used inside the four loop)
900
901 # Save the template count of the search and the time taken to conduct the search (if template_count option selected, the time taken is not indicative of
902 # how long the search takes)
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"
906
907 if not os.path.isdir(temp_count_directory):
908 os.mkdir(temp_count_directory)
909
910 # Writing the total time taken and total templates counted
911 time_taken = time.time() - start_time
912
913 with open(temp_count_directory + "/" + temp_count_file, "a+") as reader:
914 line = str(time_taken) + " " + str(template_count) + "\n"
915
916 reader.write(line)
917
918 # If we are in search mode, the lowest mismatches aren't calculated, so we can skip this step. Mismatches written in two columns. First column is metric mismatches, second is Fstat mismatches.
919 # We also don't save the largest 2Fs in search mode, as they should be saved with our templates separately.
920 if mode != "search":
921
922 # Save the highest found fstat
923 if fstat_max != 0:
924 Two_F_directory = "Largest_2Fs_For_" + tbank.toString()
925 Two_F_file = "Largest_2F_" + Random_Seed_Str + ".txt"
926
927 if not os.path.isdir(Two_F_directory):
928 os.mkdir(Two_F_directory)
929
930 with open(Two_F_directory + "/" + Two_F_file, "a+") as reader:
931 line = str(finputdata.h0) + " " + str(fstat_max) + "\n"
932
933 reader.write(line)
934
935 # Save the lowest mismatches
936 if lowest_mismatch_metric != 0:
937 mismatch_directory = "Mismatches_For_" + tbank.toString()
938 mismatch_file = "Mismatch_" + Random_Seed_Str + ".txt"
939
940 if not os.path.isdir(mismatch_directory):
941 os.mkdir(mismatch_directory)
942
943 with open(mismatch_directory + "/" + mismatch_file, "a+") as reader:
944 line = (
945 str(lowest_mismatch_metric)
946 + " "
947 + str(lowest_mismatch_Fstat)
948 + "\n"
949 )
950 reader.write(line)
951
952logging.info("Done")
953
954# For logging and profiling
955if args.profile:
956 pr.disable()
957 s = io.StringIO()
958 ps = pstats.Stats(pr, stream=s).sort_stats("cumtime")
959 ps.print_stats()
960 with open(
961 os.path.join(basedirectory, f"PiecewiseSearchProfile_{j}.txt"), "w+"
962 ) as f:
963 f.write(s.getvalue())