29 from __future__
import print_function, division
40 from scipy.integrate
import cumtrapz
41 from scipy.interpolate
import interp1d
44 from scipy.special
import logsumexp
46 from scipy.misc
import logsumexp
48 from six
import string_types
51 ARCSECTORAD = float(
"4.8481368110953599358991410235794797595635330237270e-6")
52 RADTOARCSEC = float(
"206264.80624709635515647335733077861319665970087963")
53 SECTORAD = float(
"7.2722052166430399038487115353692196393452995355905e-5")
54 RADTOSEC = float(
"13750.987083139757010431557155385240879777313391975")
55 RADTODEG = float(
"57.295779513082320876798154814105170332405472466564")
56 DEGTORAD = float(
"1.7453292519943295769236907684886127134428718885417e-2")
57 RADTOHRS = float(
"3.8197186342054880584532103209403446888270314977710")
58 HRSTORAD = float(
"2.6179938779914943653855361527329190701643078328126e-1")
59 PI = float(
"3.1415926535897932384626433832795028841971693993751")
60 TWOPI = float(
"6.2831853071795864769252867665590057683943387987502")
61 PIBYTWO = float(
"1.5707963267948966192313216916397514420985846996876")
62 SECPERDAY = float(
"86400.0")
63 SECPERJULYR = float(
"31557600.0")
64 KMPERPC = float(
"3.0856776e13")
65 KMPERKPC = float(
"3.0856776e16")
66 Tsun = float(
"4.925490947e-6")
67 Msun = float(
"1.9891e30")
68 Mjup = float(
"1.8987e27")
69 Rsun = float(
"6.9551e8")
70 Rearth = float(
"6.378e6")
71 SOL = float(
"299792458.0")
72 MSUN = float(
"1.989e+30")
73 G = float(
"6.673e-11")
75 KPC = float(
"3.0856776e19")
83 Convert radians to degrees, minutes, and seconds of arc.
89 arc = RADTODEG * np.fmod(np.fabs(rad), math.pi)
91 arc = (arc - d) * 60.0
94 if sign == -1
and d == 0:
95 return (sign * d, sign * m, sign * s)
97 return (sign * d, m, s)
102 dms_to_rad(deg, min, sec):
103 Convert degrees, minutes, and seconds of arc to radians.
107 elif deg == 0.0
and (mins < 0.0
or sec < 0.0):
114 * (60.0 * (60.0 * np.fabs(deg) + np.fabs(mins)) + np.fabs(sec))
120 dms_to_deg(deg, min, sec):
121 Convert degrees, minutes, and seconds of arc to degrees.
129 Convert radians to hours, minutes, and seconds of arc.
131 rad = np.fmod(rad, 2.0 * math.pi)
133 rad = rad + 2.0 * math.pi
136 arc = (arc - h) * 60.0
144 hms_to_rad(hour, min, sec):
145 Convert hours, minutes, and seconds of arc to radians
152 sign * SECTORAD * (60.0 * (60.0 * np.fabs(hour) + np.fabs(mins)) + np.fabs(sec))
158 coord_to_string(h_or_d, m, s):
159 Return a formatted string of RA or DEC values as
160 'hh:mm:ss.ssss' if RA, or 'dd:mm:ss.ssss' if DEC.
165 elif abs(h_or_d) == 0:
166 if (m < 0.0)
or (s < 0.0):
169 h_or_d, m, s = abs(h_or_d), abs(m), abs(s)
171 return retstr +
"%.2d:%.2d:%.4f" % (h_or_d, m, s)
173 return retstr +
"%.2d:%.2d:0%.4f" % (h_or_d, m, s)
178 rad_to_string(rad, ra_or_dec):
179 Convert an angle in radians to hours/degrees, minutes seconds and output
180 it as a string in the format 'hh:mm:ss.ssss' if RA, or 'dd:mm:ss.ssss' if DEC.
181 Whether to use hours or degrees is set by whether ra_or_dec is 'RA' or 'DEC'
183 if ra_or_dec.upper() ==
"RA":
185 elif ra_or_dec.upper() ==
"DEC":
189 "Unrecognised option: Expected 'ra_or_dec' to be 'RA' or 'DEC'"
197 ra_to_rad(ar_string):
198 Given a string containing RA information as
199 'hh:mm:ss.ssss', return the equivalent decimal
200 radians. Also deal with cases where input
201 string is just hh:mm, or hh.
203 hms = ra_string.split(
":")
211 print(
"Problem parsing RA string %s" % ra_string, file=sys.stderr)
217 dec_to_rad(dec_string):
218 Given a string containing DEC information as
219 'dd:mm:ss.ssss', return the equivalent decimal
220 radians. Also deal with cases where input string
223 dms = dec_string.split(
":")
224 if "-" in dms[0]
and float(dms[0]) == 0.0:
236 print(
"Problem parsing DEC string %s" % dec_string, file=sys.stderr)
240 def p_to_f(p, pd, pdd=None):
242 p_to_f(p, pd, pdd=None):
243 Convert period, period derivative and period second
244 derivative to the equivalent frequency counterparts.
245 Will also convert from f to p.
255 fdd = 2.0 * pd * pd / (p**3.0) - pdd / (p * p)
260 def pferrs(porf, porferr, pdorfd=None, pdorfderr=None):
262 pferrs(porf, porferr, pdorfd=None, pdorfderr=None):
263 Calculate the period or frequency errors and
264 the pdot or fdot errors from the opposite one.
267 return [1.0 / porf, porferr / porf**2.0]
269 forperr = porferr / porf**2.0
271 (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0
272 + pdorfderr**2.0 / porf**4.0
274 [forp, fdorpd] =
p_to_f(porf, pdorfd)
276 return [forp, forperr, fdorpd, fdorpderr]
392 This class parses a TEMPO(2)-style pulsar parameter file. If possible all parameters
393 are converted into SI units and angles are in radians. Epochs will be converted from
394 MJD values into GPS times.
398 for line
in pf.readlines():
404 line = line.replace(
"D-",
"E-")
405 line = line.replace(
"D+",
"E+")
407 line = line.replace(
"d-",
"e-")
408 line = line.replace(
"d+",
"e+")
409 splitline = line.split()
412 key = splitline[0].upper()
415 setattr(self, key, splitline[1])
416 elif key
in float_keys:
418 setattr(self, key, float(splitline[1]))
425 if splitline[2]
not in [
"0",
"1"]:
426 setattr(self, key +
"_ERR", float(splitline[2]))
427 setattr(self, key +
"_FIT", 0)
429 if len(splitline) == 4:
430 if splitline[2] ==
"1":
431 setattr(self, key +
"_FIT", 1)
433 setattr(self, key +
"_FIT", 0)
434 setattr(self, key +
"_ERR", float(splitline[3]))
437 if hasattr(self,
"RAJ"):
438 setattr(self,
"RA_RAD",
ra_to_rad(self.RAJ))
441 if hasattr(self,
"RAJ_ERR"):
442 setattr(self,
"RA_RAD_ERR",
hms_to_rad(0, 0, self.RAJ_ERR))
444 if hasattr(self,
"RAJ_FIT"):
445 setattr(self,
"RA_RAD_FIT", self[
"RAJ_FIT"])
446 if hasattr(self,
"DECJ"):
447 setattr(self,
"DEC_RAD",
dec_to_rad(self.DECJ))
450 if hasattr(self,
"DECJ_ERR"):
451 setattr(self,
"DEC_RAD_ERR",
dms_to_rad(0, 0, self.DECJ_ERR))
453 if hasattr(self,
"DECJ_FIT"):
454 setattr(self,
"DEC_RAD_FIT", self[
"DECJ_FIT"])
457 for pv
in [
"RA",
"DEC"]:
458 if hasattr(self,
"PM" + pv):
459 pmv = self[
"PM" + pv]
460 setattr(self,
"PM" + pv +
"_ORIGINAL", pmv)
462 self,
"PM" + pv, pmv * np.pi / (180.0 * 3600.0e3 * 365.25 * 86400.0)
465 if hasattr(self,
"PM" + pv +
"_ERR"):
466 pmve = self[
"PM" + pv +
"_ERR"]
468 self,
"PM" + pv +
"_ERR_ORIGINAL", pmv
473 pmve * np.pi / (180.0 * 3600.0e3 * 365.25 * 86400.0),
477 if hasattr(self,
"P"):
478 setattr(self,
"P0", self.P)
479 if hasattr(self,
"P0"):
480 setattr(self,
"F0", 1.0 / self.P0)
481 if hasattr(self,
"F0"):
482 setattr(self,
"P0", 1.0 / self.F0)
483 if hasattr(self,
"FB0"):
484 setattr(self,
"PB", (1.0 / self.FB0) / 86400.0)
485 if hasattr(self,
"P0_ERR"):
486 if hasattr(self,
"P1_ERR"):
487 f, ferr, fd, fderr =
pferrs(self.P0, self.P0_ERR, self.P1, self.P1_ERR)
488 setattr(self,
"F0_ERR", ferr)
489 setattr(self,
"F1", fd)
490 setattr(self,
"F1_ERR", fderr)
492 if hasattr(self,
"P1"):
496 ) =
p_to_f(self.P0, self.P1)
497 setattr(self,
"F0_ERR", self.P0_ERR / (self.P0 * self.P0))
498 setattr(self,
"F1", fd)
499 if hasattr(self,
"F0_ERR"):
500 if hasattr(self,
"F1_ERR"):
501 p, perr, pd, pderr =
pferrs(self.F0, self.F0_ERR, self.F1, self.F1_ERR)
502 setattr(self,
"P0_ERR", perr)
503 setattr(self,
"P1", pd)
504 setattr(self,
"P1_ERR", pderr)
506 if hasattr(self,
"F1"):
510 ) =
p_to_f(self.F0, self.F1)
511 setattr(self,
"P0_ERR", self.F0_ERR / (self.F0 * self.F0))
512 setattr(self,
"P1", pd)
516 from lalpulsar
import TTMJDtoGPS
527 if hasattr(self, epoch):
529 self, epoch +
"_ORIGINAL", self[epoch]
531 setattr(self, epoch, TTMJDtoGPS(self[epoch]))
537 self, epoch +
"_ERR_ORIGINAL", self[epoch +
"_ERR"]
539 setattr(self, epoch +
"_ERR", self[epoch +
"_ERR"] * SECPERDAY)
542 "Could not convert epochs to GPS times. They are all still MJD values.",
547 convfacs = {
"DIST": KPC,
"PX": 1e-3 * ARCSECTORAD}
548 for item
in convfacs:
549 if hasattr(self, item):
550 setattr(self, item +
"_ORIGINAL", self[item])
551 setattr(self, item, self[item] * convfacs[item])
553 if hasattr(self, item +
"_ERR"):
555 self, item +
"_ERR_ORIGINAL", self[item +
"_ERR"]
557 setattr(self, item +
"_ERR", self[item +
"_ERR"] * convfacs[item])
561 for om
in [
"OM",
"OM_2",
"OM_3",
"KIN",
"KOM"]:
562 if hasattr(self, om):
563 setattr(self, om, self[om] / RADTODEG)
565 if hasattr(self, om +
"_ERR"):
566 setattr(self, om +
"_ERR", self[om +
"_ERR"] / RADTODEG)
569 for pb
in [
"PB",
"PB_2",
"PB_3"]:
570 if hasattr(self, pb):
571 setattr(self, pb +
"_ORIGINAL", self[pb])
572 setattr(self, pb, self[pb] * SECPERDAY)
574 if hasattr(self, pb +
"_ERR"):
576 self, pb +
"_ERR_ORIGINAL", self[pb +
"_ERR"]
578 setattr(self, pb +
"_ERR", self[pb +
"_ERR"] * SECPERDAY)
581 if hasattr(self,
"OMDOT"):
582 setattr(self,
"OMDOT_ORIGINAL", self[
"OMDOT"])
583 setattr(self,
"OMDOT", self[
"OMDOT"] / (RADTODEG * 365.25 * SECPERDAY))
585 if hasattr(self,
"OMDOT_ERR"):
587 self,
"OMDOT_ERR_ORIGINAL", self[
"OMDOT_ERR"]
592 self[
"OMDOT_ERR"] / (RADTODEG * 365.25 * SECPERDAY),
595 if hasattr(self,
"EPS1")
and hasattr(self,
"EPS2"):
596 ecc = math.sqrt(self.EPS1 * self.EPS1 + self.EPS2 * self.EPS2)
597 omega = math.atan2(self.EPS1, self.EPS2)
598 setattr(self,
"E", ecc)
599 setattr(self,
"OM", omega)
600 setattr(self,
"T0", self.TASC + self.PB * omega / TWOPI)
603 and hasattr(self,
"A1")
604 and not (hasattr(self,
"E")
or hasattr(self,
"ECC"))
606 setattr(self,
"E", 0.0)
609 and not hasattr(self,
"TASC")
610 and hasattr(self,
"OM")
611 and hasattr(self,
"PB")
613 setattr(self,
"TASC", self.T0 - self.PB * self.OM / TWOPI)
616 for binu
in [
"XDOT",
"PBDOT",
"EDOT",
"EPS1DOT",
"EPS2DOT",
"XPBDOT"]:
617 if hasattr(self, binu):
618 setattr(self, binu +
"_ORIGINAL", self[binu])
619 if np.abs(self[binu]) > 1e-7:
620 setattr(self, binu, self[binu] * 1.0e-12)
623 if self[binu] > 10000.0:
624 setattr(self, binu, 0.0)
627 for mass
in [
"M2",
"MTOT"]:
628 if hasattr(self, mass):
629 setattr(self, mass +
"_ORIGINAL", self[mass])
630 setattr(self, mass, self[mass] * MSUN)
632 if hasattr(self, mass +
"_ERR"):
634 self, mass +
"_ERR_ORIGINAL", self[mass +
"_ERR"]
636 setattr(self, mass +
"_ERR", self[mass +
"_ERR"] * MSUN)
639 if hasattr(self,
"D_AOP"):
640 setattr(self,
"D_AOP_ORIGINAL", self[
"D_AOP"])
641 setattr(self,
"D_AOP", self[
"D_AOP"] * RADTODEG * 3600.0)
647 par = getattr(self, key)
655 for k, v
in self.__dict__.items():
657 if isinstance(self.__dict__[k], string_type):
658 out +=
"%10s = '%s'\n" % (k, v)
660 out +=
"%10s = %-20.15g\n" % (k, v)
669 pf = open(priorfilenm)
670 for line
in pf.readlines():
671 splitline = line.split()
674 key = splitline[0].upper()
678 setattr(self, key, [float(splitline[2]), float(splitline[3])])
679 elif key
in float_keys:
680 setattr(self, key, [float(splitline[2]), float(splitline[3])])
683 if hasattr(self,
"RA"):
688 setattr(self,
"RA_STR", [rastrl, rastru])
690 if hasattr(self,
"DEC"):
695 setattr(self,
"DEC_STR", [decstrl, decstru])
701 atr = getattr(self, key)
709 for k, v
in self.__dict__.items():
711 out +=
"%10s = %-20.15g, %-20.15g\n" % (k, float(v[0]), float(v[1]))
720 hsd = np.sqrt((5.0 / 2.0) * (G / C**3) * I38 * np.fabs(fdot) / freq) / (
730 ell = h0 * C**4.0 * dist * KPC / (16.0 * np.pi**2 * G * I38 * freq**2)
738 np.sqrt(15.0 / (8.0 * np.pi))
743 / (16.0 * np.pi**2 * G * freq**2)
753 * np.sqrt((8.0 * np.pi) / 15.0)
758 / (C**4.0 * dist * KPC)
767 chainlen = len(phipchain)
772 theta = math.atan2(1, 2)
776 for i
in range(0, chainlen):
777 phi0 = (1 / (2 * st)) * phipchain[i] - (1 / (2 * st)) * psipchain[i]
778 psi = (1 / (2 * ct)) * phipchain[i] + (1 / (2 * ct)) * psipchain[i]
781 if math.fabs(psi) > math.pi / 4.0:
783 phi0 = phi0 + math.pi
786 if psi > math.pi / 4.0:
787 psi = -(math.pi / 4.0) + math.fmod(psi + (math.pi / 4.0), math.pi / 2.0)
789 psi = (math.pi / 4.0) - math.fmod((math.pi / 4.0) - psi, math.pi / 2.0)
792 if phi0 > 2.0 * math.pi:
793 phi0 = math.fmod(phi0, 2.0 * math.pi)
795 phi0 = 2.0 * math.pi - math.fmod(2.0 * math.pi - phi0, 2.0 * math.pi)
797 phichain.append(phi0)
800 return phichain, psichain
810 parambounds=[float(
"-inf"), float(
"inf")],
818 from matplotlib
import pyplot
as plt
832 "axes.linewidth": 0.5,
834 "grid.linewidth": 0.5,
835 "font.family":
"serif",
839 matplotlib.rcParams.update(mplparams)
842 coldict = {
"H1":
"r",
"H2":
"c",
"L1":
"g",
"V1":
"b",
"G1":
"m",
"Joint":
"k"}
846 paraxis = paramlatexdict[param.upper()]
856 parval = parfile[param.upper()]
863 for idx, ifo
in enumerate(ifos):
865 if overplot
and idx == 0:
866 myfig = plt.figure(figsize=(4, 4), dpi=200)
869 myfig = plt.figure(figsize=(4, 4), dpi=200)
873 pos_samps = pos[param].samples
877 pos_samps,
int(nbins), parambounds[0], parambounds[1]
881 plt.step(bins, n, color=coldict[ifo])
883 if "h0" not in param:
884 plt.xlim(parambounds[0], parambounds[1])
886 plt.xlabel(
r"" + paraxis, fontsize=14, fontweight=100)
887 plt.ylabel(
r"Probability Density", fontsize=14, fontweight=100)
888 myfig.subplots_adjust(left=0.18, bottom=0.15)
891 plt.ylim(0, n.max() + 0.1 * n.max())
895 ax.set_axis_bgcolor(
"#F2F1F0")
898 ymax.append(n.max() + 0.1 * n.max())
902 ct = cumtrapz(n, bins)
905 ct = np.insert(ct, 0, 0)
908 ctu, ui = np.unique(ct, return_index=
True)
909 intf = interp1d(ctu, bins[ui], kind=
"linear")
910 ulvals.append(intf(float(upperlimit)))
916 plt.axvline(parval, color=
"k", ls=
"--", linewidth=1.5)
918 plt.axvline(parval, color=
"k", ls=
"--", linewidth=1.5)
921 plt.ylim(0,
max(ymax))
924 ax.set_axis_bgcolor(
"#F2F1F0")
928 return myfigs, ulvals
934 pos, upperlimit=0.95, parambounds=[float(
"-inf"), float(
"inf")], nbins=50
943 ct = cumtrapz(n, bins)
946 ct = np.insert(ct, 0, 0)
949 ctu, ui = np.unique(ct, return_index=
True)
950 intf = interp1d(ctu, bins[ui], kind=
"linear")
951 ulval = intf(float(upperlimit))
957 n, binedges = np.histogram(pos, bins=nbins)
958 dbins = binedges[1] - binedges[0]
964 frac += float(nv) / len(pos)
966 if frac > upperlimit:
970 ul = binedges[j - 1] + (upperlimit - prevfrac) * (dbins / (frac - prevfrac))
983 from matplotlib
import pyplot
as plt
987 from matplotlib
import gridspec
995 "axes.linewidth": 0.5,
997 "grid.linewidth": 0.5,
998 "font.family":
"sans-serif",
999 "font.sans-serif":
"Avant Garde, Helvetica, Computer Modern Sans serif",
1003 matplotlib.rcParams.update(mplparams)
1005 coldict = {
"H1":
"r",
"H2":
"c",
"L1":
"g",
"V1":
"b",
"G1":
"m",
"Joint":
"k"}
1014 paryaxis = paramlatexdict[p.upper()]
1023 minsamp = float(
"inf")
1024 maxsamp = -float(
"inf")
1026 for idx, ifo
in enumerate(ifos):
1028 myfig = plt.figure(figsize=(12, 4), dpi=200)
1029 myfig.subplots_adjust(bottom=0.15)
1032 gs = gridspec.GridSpec(1, 4, wspace=0)
1033 ax1 = plt.subplot(gs[:-1])
1034 ax2 = plt.subplot(gs[-1])
1036 pos = list(poslist)[idx]
1040 pos_samps = np.cos(pos[
"iota"].samples)
1042 pos_samps = pos[param].samples
1044 if np.min(pos_samps) < minsamp:
1045 minsamp = np.min(pos_samps)
1046 if np.max(pos_samps) > maxsamp:
1047 maxsamp = np.max(pos_samps)
1050 ax1.plot(pos_samps,
".", color=coldict[ifo], markersize=1)
1052 n, binedges = np.histogram(pos_samps, withhist, density=
True)
1054 ax2.step(n, binedges, color=coldict[ifo])
1056 if np.max(n) > maxn:
1059 plt.plot(pos_samps,
".", color=coldict[ifo], markersize=1)
1063 legendvals.append(
r"$R = %.2f$" % grr[idx][param])
1069 if len(pos_samps) > maxiter:
1070 maxiter = len(pos_samps)
1075 bounds = [minsamp, maxsamp]
1077 ax1.set_ylabel(
r"" + paryaxis, fontsize=16, fontweight=100)
1078 ax1.set_xlabel(
r"Iterations", fontsize=16, fontweight=100)
1080 ax1.set_xlim(0, maxiter)
1081 ax1.set_ylim(bounds[0], bounds[1])
1084 ax2.set_ylim(bounds[0], bounds[1])
1085 ax2.set_xlim(0, maxn + 0.1 * maxn)
1086 ax2.set_yticklabels([])
1087 ax2.set_axis_bgcolor(
"#F2F1F0")
1091 ax1.legend(legendvals, title=
"Gelman-Rubins test")
1111 fp = open(histfile,
"rb")
1113 print(
"Could not open prior file %s" % histfile, file=sys.stderr)
1114 return None,
None,
None
1119 print(
"Could not read in data from prior file %s" % histfile, file=sys.stderr)
1120 return None,
None,
None
1126 header = struct.unpack(
"d" * 6, pd[: 6 * 8])
1133 grid = struct.unpack(
"d" *
int(header[2]) *
int(header[5]), pd[6 * 8 :])
1138 header = list(header)
1142 histarr = np.array([g[:
int(header[5])]], ndmin=2)
1143 for i
in range(
int(header[2]) - 1):
1144 histarr = np.append(
1145 histarr, [g[(i + 1) *
int(header[5]) : (i + 2) *
int(header[5])]], axis=0
1148 xbins = np.linspace(
1149 header[0], header[0] + header[1] * (header[2] - 1),
int(header[2])
1151 ybins = np.linspace(
1152 header[3], header[3] + header[4] * (header[5] - 1),
int(header[5])
1155 return xbins, ybins, histarr
1163 histfile, ndimlabel, mdimlabel, margpars=True, mplparams=False
1166 from matplotlib
import pyplot
as plt
1173 print(
"Could not read binary histogram file", file=sys.stderr)
1182 "text.usetex":
True,
1183 "axes.linewidth": 0.5,
1185 "grid.linewidth": 0.5,
1186 "font.family":
"serif",
1190 matplotlib.rcParams.update(mplparams)
1192 fig = plt.figure(figsize=(4, 4), dpi=200)
1196 parxaxis = paramlatexdict[ndimlabel.upper()]
1198 parxaxis = ndimlabel
1201 paryaxis = paramlatexdict[mdimlabel.upper()]
1203 paryaxis = mdimlabel
1205 plt.xlabel(
r"" + parxaxis, fontsize=14, fontweight=100)
1206 plt.ylabel(
r"" + paryaxis, fontsize=14, fontweight=100, rotation=270)
1208 dx = xbins[1] - xbins[0]
1209 dy = ybins[1] - ybins[0]
1212 xbins[0] - dx / 2.0,
1213 xbins[-1] + dx / 2.0,
1214 ybins[-1] + dy / 2.0,
1215 ybins[0] - dy / 2.0,
1219 np.transpose(histarr),
1222 interpolation=
"bicubic",
1226 fig.subplots_adjust(left=0.18, bottom=0.15)
1227 fax = (fig.get_axes())[0].
axis()
1234 for i
in range(len(xbins)):
1235 xmarg.append(np.trapz(histarr[:][i], x=ybins))
1238 xarea = np.trapz(xmarg, x=xbins)
1239 xmarg = map(
lambda x: x / xarea, xmarg)
1242 for i
in range(len(ybins)):
1243 ymarg.append(np.trapz(np.transpose(histarr)[:][i], x=xbins))
1246 yarea = np.trapz(ymarg, x=ybins)
1247 ymarg = map(
lambda x: x / yarea, ymarg)
1250 figx = plt.figure(figsize=(4, 4), dpi=200)
1251 plt.step(xbins, xmarg, color=
"k")
1252 plt.ylim(0,
max(xmarg) + 0.1 *
max(xmarg))
1253 plt.xlim(fax[0], fax[1])
1254 plt.xlabel(
r"" + parxaxis, fontsize=14, fontweight=100)
1255 plt.ylabel(
r"Probability Density", fontsize=14, fontweight=100)
1256 figx.subplots_adjust(left=0.18, bottom=0.15)
1258 ax.set_axis_bgcolor(
"#F2F1F0")
1261 figy = plt.figure(figsize=(4, 4), dpi=200)
1262 plt.step(ybins, ymarg, color=
"k")
1263 plt.ylim(0,
max(ymarg) + 0.1 *
max(ymarg))
1264 plt.xlim(fax[3], fax[2])
1265 plt.xlabel(
r"" + paryaxis, fontsize=14, fontweight=100)
1266 plt.ylabel(
r"Probability Density", fontsize=14, fontweight=100)
1267 figy.subplots_adjust(left=0.18, bottom=0.15)
1269 ax.set_axis_bgcolor(
"#F2F1F0")
1283 if not h0bins.any():
1284 print(
"Could not read binary histogram file", file=sys.stderr)
1289 for i
in range(len(h0bins)):
1290 h0marg.append(np.trapz(histarr[:][i], x=cibins))
1293 h0bins = h0bins - (h0bins[1] - h0bins[0]) / 2
1294 h0area = np.trapz(h0marg, x=h0bins)
1295 h0margnorm = map(
lambda x: x / h0area, h0marg)
1298 ct = cumtrapz(h0margnorm, h0bins)
1301 ct = np.insert(ct, 0, 0)
1304 ctu, ui = np.unique(ct, return_index=
True)
1305 intf = interp1d(ctu, h0bins[ui], kind=
"linear")
1306 return intf(float(ulval))
1321 from matplotlib
import pyplot
as plt
1324 if len(params) != 2:
1325 print(
"Require 2 parameters", file=sys.stderr)
1332 "text.usetex":
True,
1333 "axes.linewidth": 0.5,
1335 "grid.linewidth": 0.5,
1336 "font.family":
"serif",
1340 matplotlib.rcParams.update(mplparams)
1356 parxaxis = paramlatexdict[params[0].upper()]
1358 parxaxis = params[0]
1361 paryaxis = paramlatexdict[params[1].upper()]
1363 paryaxis = params[1]
1369 parval1 = parfile[params[0].upper()]
1370 parval2 = parfile[params[1].upper()]
1375 for idx, ifo
in enumerate(ifos):
1379 myfig = plt.figure(figsize=(4, 4), dpi=200)
1383 cmap = plt.cm.get_cmap(
"Greys")
1386 cmap = plt.cm.get_cmap(coldict[ifo])
1389 myfig = plt.figure(figsize=(4, 4), dpi=200)
1391 cmap = plt.cm.get_cmap(
"Greys")
1393 posterior = poslist[idx]
1395 a = np.squeeze(posterior[params[0]].samples)
1396 b = np.squeeze(posterior[params[1]].samples)
1399 par1pos_min = a.min()
1400 par2pos_min = b.min()
1402 par1pos_max = a.max()
1403 par2pos_max = b.max()
1405 plt.xlabel(
r"" + parxaxis, fontsize=14, fontweight=100)
1406 plt.ylabel(
r"" + paryaxis, fontsize=14, fontweight=100, rotation=270)
1408 H, xedges, yedges = np.histogram2d(a, b, nbins, normed=
True)
1415 Hm = np.transpose(H)
1417 extent = [xedges[0], xedges[-1], yedges[-1], yedges[0]]
1422 interpolation=
"bilinear",
1429 plt.xlim(bounds[0][0], bounds[0][1])
1430 plt.ylim(bounds[1][0], bounds[1][1])
1433 if parval1
and parval2:
1434 if overplot
and idx == len(ifos) - 1:
1437 plt.plot(parval1, parval2,
"rx", markersize=8, mew=2)
1440 plt.plot(parval1, parval2,
"rx", markersize=8, mew=2)
1445 myfig.subplots_adjust(left=0.18, bottom=0.15)
1446 myfigs.append(myfig)
1449 myfig.subplots_adjust(left=0.18, bottom=0.15)
1450 myfigs.append(myfig)
1460 n, binedges = np.histogram(samples, nbins)
1463 binwidth = binedges[1] - binedges[0]
1466 bincentres = np.array([])
1467 for i
in range(0, len(binedges) - 1):
1468 bincentres = np.append(bincentres, binedges[i] + binwidth / 2)
1472 if bincentres[0] - binwidth > low:
1474 n = np.insert(n, 0, 0)
1477 bincentres = np.insert(bincentres, 0, bincentres[0] - binwidth)
1482 dx = bincentres[0] - low
1485 bincentres = np.insert(bincentres, 0, low)
1489 nbound = n[0] - (dn / binwidth) * dx
1494 n = np.insert(n, 0, nbound)
1497 if bincentres[-1] + binwidth < high:
1502 bincentres = np.append(bincentres, bincentres[-1] + binwidth)
1504 dx = high - bincentres[-1]
1507 bincentres = np.append(bincentres, high)
1511 nbound = n[-1] + (dn / binwidth) * dx
1516 n = np.append(n, nbound)
1519 area = np.trapz(n, x=bincentres)
1522 for i
in range(0, len(bincentres)):
1523 ns = np.append(ns, float(n[i]) / area)
1525 return ns, bincentres
1532 return np.hanning(N)
1535 x = np.linspace(0, 1, N)
1538 win = np.ones(x.shape)
1542 win[lhs] = 0.5 * (1 + np.cos(2 * np.pi / alpha * (x[lhs] - alpha / 2)))
1545 rhs = x >= (1 - alpha / 2)
1546 win[rhs] = 0.5 * (1 + np.cos(2 * np.pi / alpha * (x[rhs] - 1 + alpha / 2)))
1564 from matplotlib.mlab
import specgram
1565 from matplotlib
import colors
1566 from matplotlib
import pyplot
as plt
1579 "text.usetex":
True,
1580 "axes.linewidth": 0.5,
1582 "grid.linewidth": 0.5,
1583 "font.family":
"sans-serif",
1584 "font.sans-serif":
"Avant Garde, Helvetica, Computer Modern Sans serif",
1588 matplotlib.rcParams.update(mplparams)
1589 matplotlib.rcParams[
"text.latex.preamble"] =
r"\usepackage{xfrac}"
1592 coldict = {
"H1":
"r",
"H2":
"c",
"L1":
"g",
"V1":
"b",
"G1":
"m"}
1605 Bk = np.loadtxt(Bkdata[ifo], comments=[
"#",
"%"])
1607 print(
"Could not open file %s" % Bkdata[ifo], file=sys.stderr)
1611 Bk = Bk[Bk[:, 0].argsort()]
1614 uargs = np.unique(Bk[:, 0], return_index=
True)[1]
1623 n, binedges = np.histogram(
1624 np.log(np.fabs(np.concatenate((Bk[:, 1], Bk[:, 2])))), 50
1629 stdest = math.exp((binedges[j] + binedges[j + 1]) / 2)
1633 vals = np.where(np.fabs(Bk[:, 1]) < removeoutlier * stdest)
1634 Bknew = Bk[vals, :][-1]
1636 vals = np.where(np.fabs(Bknew[:, 2]) < removeoutlier * stdest)
1637 Bk = Bknew[vals, :][-1]
1642 mindt =
min(np.diff(gpstime))
1643 sampledt.append(mindt)
1645 Bkabs = np.sqrt(Bk[:, 1] ** 2 + Bk[:, 2] ** 2)
1648 Bkfig = plt.figure(figsize=(11, 3.5), dpi=200)
1649 Bkfig.subplots_adjust(bottom=0.15, left=0.09, right=0.94)
1651 tms = list(map(
lambda x: x - gpstime[0], gpstime))
1653 plt.plot(tms, Bkabs,
".", color=coldict[ifo], markersize=1)
1654 plt.xlabel(
r"GPS - %d" %
int(gpstime[0]), fontsize=14, fontweight=100)
1655 plt.ylabel(
r"$|B_k|$", fontsize=14, fontweight=100)
1657 plt.xlim(tms[0], tms[-1])
1659 Bkfigs.append(Bkfig)
1661 if plotpsds
or plotfscan:
1664 totlen = gpstime[-1] - gpstime[0]
1667 if math.fmod(mindt, 1) != 0.0
or mindt < 1:
1669 "Error... time steps between data points must be integers",
1677 datazeropad = np.zeros(
int(math.ceil(totlen / mindt)) + 1, dtype=complex)
1679 idx = list(map(
lambda x:
int(math.floor((x / mindt) + 0.5)), tms))
1680 for i
in range(0, len(idx)):
1681 datazeropad[idx[i]] = complex(Bk[i, 1], Bk[i, 2])
1683 win =
tukey_window(math.floor(delt / mindt), alpha=0.1)
1687 fscan, freqs, t = specgram(
1689 NFFT=
int(math.floor(delt / mindt)),
1692 noverlap=
int(math.floor(delt / (2.0 * mindt))),
1696 fshape = fscan.shape
1698 totalasd = np.zeros(fshape[0])
1700 for i
in range(0, fshape[0]):
1701 scanasd = np.sqrt(fscan[i, :])
1702 nonzasd = np.nonzero(scanasd)
1703 scanasd = scanasd[nonzasd]
1707 totalasd[i] = np.median(scanasd)
1710 asdlist.append(totalasd)
1713 psdfig = plt.figure(figsize=(4, 3.5), dpi=200)
1714 psdfig.subplots_adjust(left=0.18, bottom=0.15)
1716 plt.plot(freqs, totalasd, color=coldict[ifo])
1717 plt.xlim(freqs[0], freqs[-1])
1718 plt.xlabel(
r"Frequency (Hz)", fontsize=14, fontweight=100)
1719 plt.ylabel(
r"$h/\sqrt{\rm Hz}$", fontsize=14, fontweight=100)
1723 xt = [-Fs / 2.0, -Fs / 4.0, 0.0, Fs / 2.0, Fs / 4.0]
1731 xl.append(
r"$-\sfrac{1}{%d}$" % (-1.0 / item))
1733 xl.append(
r"$\sfrac{1}{%d}$" % (1.0 / item))
1734 ax.set_xticklabels(xl)
1736 plt.tick_params(axis=
"x", which=
"major", labelsize=14)
1738 psdfigs.append(psdfig)
1741 fscanfig = plt.figure(figsize=(11, 3.5), dpi=200)
1742 fscanfig.subplots_adjust(bottom=0.15, left=0.09, right=0.94)
1744 extent = [tms[0], tms[-1], freqs[0], freqs[-1]]
1746 np.sqrt(np.flipud(fscan)),
1750 cmap=colmapdic[ifo],
1751 norm=colors.Normalize(),
1753 plt.ylabel(
r"Frequency (Hz)", fontsize=14, fontweight=100)
1754 plt.xlabel(
r"GPS - %d" %
int(gpstime[0]), fontsize=14, fontweight=100)
1758 yt = [-Fs / 2.0, -Fs / 4.0, 0.0, Fs / 2.0, Fs / 4.0]
1766 yl.append(
r"$-\sfrac{1}{%d}$" % (-1.0 / item))
1768 yl.append(
r"$\sfrac{1}{%d}$" % (1.0 / item))
1769 ax.set_yticklabels(yl)
1770 plt.tick_params(axis=
"y", which=
"major", labelsize=14)
1772 fscanfigs.append(fscanfig)
1774 return Bkfigs, psdfigs, fscanfigs, asdlist, sampledt
1781 lims, param, ifos, prevlims=None, bins=20, overplot=False, mplparams=False
1784 from matplotlib
import pyplot
as plt
1790 "text.usetex":
True,
1791 "axes.linewidth": 0.5,
1793 "grid.linewidth": 0.5,
1794 "font.family":
"serif",
1798 matplotlib.rcParams.update(mplparams)
1803 coldict = {
"H1":
"r",
"H2":
"c",
"L1":
"g",
"V1":
"b",
"G1":
"m",
"Joint":
"k"}
1806 parxaxis = paramlatexdict[param.upper()]
1810 paryaxis =
"Number of Pulsars"
1814 if prevlims
is not None:
1815 logprevlims = np.log10(prevlims)
1828 for j, ifo
in enumerate(ifos):
1829 theselims = lims[ifo]
1832 for i, val
in enumerate(theselims):
1836 loglims = np.log10(theselims)
1838 stackdata.append(loglims)
1840 highbins.append(np.max(loglims))
1841 lowbins.append(np.min(loglims))
1843 if logprevlims
is not None:
1844 highbins.append(
max(logprevlims))
1845 lowbins.append(
min(logprevlims))
1847 hrange = (
min(lowbins),
max(highbins))
1849 for j, ifo
in enumerate(ifos):
1851 theselims = lims[ifo]
1854 for i, val
in enumerate(theselims):
1858 loglims = np.log10(theselims)
1864 myfig = plt.figure(figsize=(4, 4), dpi=200)
1871 for ifoname
in ifos:
1872 edgecolor.append(coldict[ifoname])
1873 facecolor.append(coldict[ifoname])
1875 edgecolor = [coldict[ifo]]
1876 facecolor = [coldict[ifo]]
1879 n, bins, patches = plt.hist(
1889 for i, patch
in enumerate(patches):
1890 plt.setp(patch,
"facecolor", facecolor[i])
1893 maxlims.append(np.max(loglims))
1894 minlims.append(np.min(loglims))
1896 if logprevlims
is not None:
1898 maxlims.append(
max(logprevlims))
1899 minlims.append(
min(logprevlims))
1901 maxlim =
max(maxlims)
1902 minlim =
min(minlims)
1920 maxlim =
max(maxlims)
1921 minlim =
min(minlims)
1930 tickvals =
range(
int(math.ceil(maxlim) - math.floor(minlim)))
1931 tickvals = map(
lambda x: x +
int(math.floor(minlim)), tickvals)
1932 ax.set_xticks(tickvals)
1933 tls = map(
lambda x:
"$10^{%d}$" % x, tickvals)
1934 ax.set_xticklabels(tls)
1936 plt.ylabel(
r"" + paryaxis, fontsize=14, fontweight=100)
1937 plt.xlabel(
r"" + parxaxis, fontsize=14, fontweight=100)
1939 myfig.subplots_adjust(left=0.18, bottom=0.15)
1941 myfigs.append(myfig)
1972 from matplotlib
import pyplot
as plt
1977 "text.usetex":
True,
1978 "axes.linewidth": 0.5,
1980 "grid.linewidth": 0.5,
1981 "font.family":
"serif",
1985 matplotlib.rcParams.update(mplparams)
1990 coldict = {
"H1":
"r",
"H2":
"c",
"L1":
"g",
"V1":
"b",
"G1":
"m",
"Joint":
"k"}
1992 parxaxis =
"Frequency (Hz)"
1995 for j, ifo
in enumerate(ifos):
2003 if len(h0lim) != len(f0s):
2006 if not overplot
or (overplot
and j == 0):
2007 myfig = plt.figure(figsize=(7, 5.5), dpi=200)
2010 if ulesttop
is not None and ulestbot
is not None:
2014 if len(ult) == len(ulb)
and len(ult) == len(f0s):
2015 for i
in range(len(ult)):
2017 [f0s[i], f0s[i]], [ulb[i], ult[i]], ls=
"-", lw=3, c=
"lightgrey"
2022 and prevlimf0gw
is not None
2023 and len(prevlim) == len(prevlimf0gw)
2025 if not overplot
or (overplot
and j == len(ifos) - 1):
2039 f0s, h0lim, marker=
"*", ms=10, mfc=coldict[ifo], mec=coldict[ifo], ls=
"None"
2042 plt.ylabel(
r"" + paryaxis, fontsize=14, fontweight=100)
2043 plt.xlabel(
r"" + parxaxis, fontsize=14, fontweight=100)
2045 plt.xlim(xlims[0], xlims[1])
2047 if not overplot
or (overplot
and j == len(ifos) - 1):
2049 myfig.subplots_adjust(left=0.12, bottom=0.10)
2050 myfigs.append(myfig)
2067 s = np.array([], dtype=complex)
2069 sphi = np.sin(pardict[
"phi0"])
2070 cphi = np.cos(pardict[
"phi0"])
2072 Xplus = 0.25 * (1.0 + pardict[
"cosiota"] * pardict[
"cosiota"]) * pardict[
"h0"]
2073 Xcross = 0.5 * pardict[
"cosiota"] * pardict[
"h0"]
2074 Xpsinphi = Xplus * sphi
2075 Xcsinphi = Xcross * sphi
2076 Xpcosphi = Xplus * cphi
2077 Xccosphi = Xcross * cphi
2080 while tmpts < starttime + duration:
2081 ts.append(starttime + (dt * i))
2085 ts[i], pardict[
"ra"], pardict[
"dec"], pardict[
"psi"], detector
2090 s, (fp * Xpcosphi + fc * Xcsinphi) + 1j * (fp * Xpsinphi - fc * Xccosphi)
2112 pardict, detector, starttime=900000000.0, duration=86400.0, dt=60.0, datatimes=None
2114 if "cosiota" in pardict:
2115 cosiota = pardict[
"cosiota"]
2116 iota = math.acos(cosiota)
2117 siniota = math.sin(iota)
2119 raise KeyError(
"cos(iota) not defined!")
2121 if "psi" in pardict:
2122 psi = pardict[
"psi"]
2124 raise KeyError(
"psi not defined!")
2126 if "C22" in pardict:
2127 C22 = pardict[
"C22"]
2130 C22 = -pardict[
"h0"] / 2.0
2134 if "phi22" in pardict:
2135 phi22 = pardict[
"phi22"]
2137 if "phi0" in pardict:
2138 phi22 = 2.0 * pardict[
"phi0"]
2142 ePhi22 = cmath.exp(phi22 * 1j)
2144 if "C21" in pardict:
2145 C21 = pardict[
"C21"]
2149 if "phi21" in pardict:
2150 phi21 = pardict[
"phi21"]
2154 ePhi21 = cmath.exp(phi21 * 1j)
2159 if "C21" in pardict
and "C22" in pardict
and "h0" not in pardict:
2166 if datatimes
is None:
2167 datatimes = np.arange(starttime, starttime + duration, dt)
2171 datatimes, pardict[
"ra"], pardict[
"dec"], pardict[
"psi"], detector
2176 -(C21 / 4.0) * ePhi21 * siniota * cosiota * fp
2177 + 1j * (C21 / 4.0) * ePhi21 * siniota * fc
2181 -(C22 / 2.0) * ePhi22 * (1.0 + cosiota**2.0) * fp
2182 + 1j * (C22) * ePhi22 * cosiota * fc
2185 ts.append(datatimes)
2194 costheta = pardict[
"costheta"]
2195 theta = np.arccos(costheta)
2196 sintheta = math.sin(theta)
2197 sin2theta = math.sin(2.0 * theta)
2198 sinlambda = math.sin(pardict[
"lambdapin"])
2199 coslambda = math.cos(pardict[
"lambdapin"])
2200 sin2lambda = math.sin(2.0 * pardict[
"lambdapin"])
2202 phi0 = pardict[
"phi0"]
2204 I21 = pardict[
"I21"]
2205 I31 = pardict[
"I31"]
2207 A22 = I21 * (sinlambda**2 - (coslambda * costheta) ** 2) - I31 * sintheta**2
2208 B22 = I21 * sin2lambda * costheta
2210 C22 = 2.0 * math.sqrt(A22**2 + B22**2)
2212 A21 = I21 * sin2lambda * sintheta
2213 B21 = (I21 * coslambda**2 - I31) * sin2theta
2215 C21 = 2.0 * math.sqrt(A21**2 + B21**2)
2217 phi22 = np.fmod(2.0 * phi0 - math.atan2(B22, A22), 2.0 * np.pi)
2218 phi21 = np.fmod(phi0 - math.atan2(B21, A21), 2.0 * np.pi)
2220 outvals = {
"C22": C22,
"C21": C21,
"phi22": phi22,
"phi21": phi21}
2233 iota = np.arccos(pardict[
"cosiota"])
2234 theta = np.arccos(pardict[
"costheta"])
2235 siniota = math.sin(iota)
2236 sintheta = math.sin(theta)
2237 sin2theta = math.sin(2.0 * theta)
2238 sinlambda = math.sin(pardict[
"lambdapin"])
2239 coslambda = math.cos(pardict[
"lambdapin"])
2240 sin2lambda = math.sin(2.0 * pardict[
"lambdapin"])
2242 ePhi = cmath.exp(0.5 * pardict[
"phi0"] * 1j)
2243 e2Phi = cmath.exp(pardict[
"phi0"] * 1j)
2245 f2_r = pardict[
"f0"] ** 2 / pardict[
"dist"]
2248 This model is a complex heterodyned time series for a pinned superfluid
2249 neutron star emitting at its roation frequency and twice its rotation
2250 frequency (as defined in Eqns 35-38 of Jones 2012 LIGO DCC T1200265-v3)
2253 Xplusf = -(f2_r / 2.0) * siniota * pardict[
"cosiota"]
2254 Xcrossf = (f2_r / 2.0) * siniota
2256 Xplus2f = -f2_r * (1.0 + pardict[
"cosiota"] ** 2)
2257 Xcross2f = 2.0 * f2_r * pardict[
"cosiota"]
2259 A21 = pardict[
"I21"] * sin2lambda * sintheta
2260 B21 = (pardict[
"I21"] * coslambda**2 - pardict[
"I31"]) * sin2theta
2263 pardict[
"I21"] * (sinlambda**2 - coslambda**2 * pardict[
"costheta"] ** 2)
2264 - pardict[
"I31"] * sintheta**2
2266 B22 = pardict[
"I21"] * sin2lambda * pardict[
"costheta"]
2274 s1 = np.array([], dtype=complex)
2275 s2 = np.array([], dtype=complex)
2278 while tmpts < starttime + duration:
2279 ts1.append(starttime + (dt * i))
2280 ts2.append(starttime + (dt * i))
2284 ts1[i], pardict[
"ra"], pardict[
"dec"], pardict[
"psi"], detector
2290 (fp * Xplusf * ePhi * (A21 - 1j * B21))
2291 + (fc * Xcrossf * ePhi * (B21 + 1j * A21)),
2297 (fp * Xplus2f * e2Phi * (A22 - 1j * B22))
2298 + (fc * Xcross2f * e2Phi * (B22 + 1j * A22)),
2306 ts = np.vstack([ts1, ts2])
2307 s = np.vstack([s1, s2])
2320 "H1": lal.LALDetectorIndexLHODIFF,
2321 "H2": lal.LALDetectorIndexLHODIFF,
2322 "L1": lal.LALDetectorIndexLLODIFF,
2323 "G1": lal.LALDetectorIndexGEO600DIFF,
2324 "V1": lal.LALDetectorIndexVIRGODIFF,
2325 "T1": lal.LALDetectorIndexTAMA300DIFF,
2326 "AL1": lal.LALDetectorIndexLLODIFF,
2327 "AH1": lal.LALDetectorIndexLHODIFF,
2328 "AV1": lal.LALDetectorIndexVIRGODIFF,
2329 "E1": lal.LALDetectorIndexE1DIFF,
2330 "E2": lal.LALDetectorIndexE2DIFF,
2331 "E3": lal.LALDetectorIndexE3DIFF,
2335 detector = detMap[det]
2337 raise KeyError(
"ERROR. Key {} is not a valid detector name.".format(det))
2340 detval = lal.CachedDetectors[detector]
2341 response = detval.response
2344 if not isinstance(ra, float):
2345 if isinstance(ra, string_types):
2350 "Right ascension string '{}' not formatted properly".format(ra)
2354 "Right ascension must be a 'float' in radians or a string of the format 'hh:mm:ss.s'"
2357 if not isinstance(dec, float):
2358 if isinstance(dec, string_types):
2363 "Declination string '{}' not formatted properly".format(ra)
2367 "Declination must be a 'float' in radians or a string of the format 'dd:mm:ss.s'"
2371 if isinstance(gpsTime, float)
or isinstance(gpsTime, int):
2372 gpsTime = np.array([gpsTime])
2374 gpsTime = np.copy(gpsTime)
2377 gpsTime = gpsTime.astype(
"float64")
2380 if len(gpsTime) == 1
or np.unique(np.diff(gpsTime)).size == 1:
2381 gpsStart = lal.LIGOTimeGPS(gpsTime[0])
2383 if len(gpsTime) > 1:
2384 dt = gpsTime[1] - gpsTime[0]
2385 fp, fc = lal.ComputeDetAMResponseSeries(
2386 response, ra, dec, psi, gpsStart, dt, len(gpsTime)
2390 return fp.data.data, fc.data.data
2392 fp = np.zeros(len(gpsTime))
2393 fc = np.zeros(len(gpsTime))
2395 for i
in range(len(gpsTime)):
2396 gps = lal.LIGOTimeGPS(gpsTime[i])
2397 gmst_rad = lal.GreenwichMeanSiderealTime(gps)
2400 fp[i], fc[i] = lal.ComputeDetAMResponse(response, ra, dec, psi, gmst_rad)
2420 if isinstance(detectors, string_types):
2421 detectors = [detectors]
2428 for det
in detectors:
2434 ns = np.sqrt((psd / 2.0) / (2.0 * dt))
2442 for j, det
in enumerate(detectors):
2445 tmpnpsds.append(np.sqrt((npsds[0] / 2.0) / (2.0 * dt)))
2447 tmpnpsds.append(np.sqrt((npsds[count] / 2.0) / (2.0 * dt)))
2453 if len(freqfac) == 1
and len(detectors) != len(npsds):
2455 "Number of detectors {} not the same as number of noises {}".format(
2456 len(detectors), len(npsds)
2460 if len(freqfac) == 2
and 2 * len(detectors) != len(npsds):
2462 "Number of detectors {} not half the number of noises {}".format(
2463 len(detectors), len(npsds)
2471 for j, det
in enumerate(detectors):
2473 if len(freqfac) == 1
and freqfac[0] == 2.0:
2477 tss = np.append(tss, ts)
2478 ss = np.append(ss, s)
2480 tss = np.vstack([tss, ts])
2481 ss = np.vstack([ss, s])
2485 elif len(freqfac) == 2:
2489 for k, frf
in enumerate(freqfac):
2490 if j == 0
and k == 0:
2491 tss = np.append(tss, ts[k][:])
2492 ss = np.append(ss, s[k][:])
2494 tss = np.vstack([tss, ts[k][:]])
2495 ss = np.vstack([ss, s[k][:]])
2498 snrtmp = snrtmp + snrtmp2 * snrtmp2
2500 snrtmp = np.sqrt(snrtmp)
2502 snrtot = snrtot + snrtmp * snrtmp
2505 snrtot = np.sqrt(snrtot)
2508 if snrscale
is not None:
2510 snrscale = snrscale / snrtot
2515 signalonly = ss.copy()
2518 for det
in detectors:
2520 if len(freqfac) == 1:
2522 rs = np.random.randn(len(ts[0]), 2)
2524 for j, t
in enumerate(ts[0]):
2525 if len(tss.shape) == 1:
2526 signalonly[j] = (snrscale * ss[j].real) + 1j * (
2527 snrscale * ss[j].imag
2529 ss[j] = (snrscale * ss[j].real + npsds[i] * rs[j][0]) + 1j * (
2530 snrscale * ss[j].imag + npsds[i] * rs[j][1]
2533 signalonly[i][j] = (snrscale * ss[i][j].real) + 1j * (
2534 snrscale * ss[i][j].imag
2536 ss[i][j] = (snrscale * ss[i][j].real + npsds[i] * rs[j][0]) + 1j * (
2537 snrscale * ss[i][j].imag + npsds[i] * rs[j][1]
2541 elif len(freqfac) == 2:
2543 rs = np.random.randn(len(ts[0][:]), 4)
2545 for j, t
in enumerate(ts[0][:]):
2546 signalonly[i][j] = (snrscale * ss[i][j].real) + 1j * (
2547 snrscale * ss[i][j].imag
2549 ss[i][j] = (snrscale * ss[i][j].real + npsds[i] * rs[j][0]) + 1j * (
2550 snrscale * ss[i][j].imag + npsds[i] * rs[j][1]
2552 signalonly[i + 1][j] = (snrscale * ss[i + 1][j].real) + 1j * (
2553 snrscale * ss[i + 1][j].imag
2556 snrscale * ss[i + 1][j].real + npsds[i + 1] * rs[j][2]
2557 ) + 1j * (snrscale * ss[i + 1][j].imag + npsds[i + 1] * rs[j][3])
2561 print(
"Something wrong with injection", file=sys.stderr)
2564 snrtot = snrtot * snrscale
2566 return tss, ss, signalonly, snrtot, snrscale, npsds
2575 import lalsimulation
2578 return lalsimulation.SimNoisePSDAdvVirgo(f)
2579 elif det ==
"H1" or det ==
"L1":
2580 return lalsimulation.SimNoisePSDiLIGOSRD(f)
2582 return lalsimulation.SimNoisePSDiLIGOSRD(f) * 2.0
2584 return lalsimulation.SimNoisePSDGEO(f)
2586 return lalsimulation.SimNoisePSDVirgo(f)
2588 return lalsimulation.SimNoisePSDTAMA(f)
2590 return lalsimulation.SimNoisePSDKAGRA(f)
2591 elif det ==
"AL1" or det ==
"AH1":
2592 return lalsimulation.SimNoisePSDaLIGOZeroDetHighPower(f)
2594 raise ValueError(
"{} is not a recognised detector".format(det))
2603 ss = ss + val.real**2 + val.imag**2
2605 return np.sqrt(ss / sig**2)
2612 chainMeans = [np.mean(data)
for data
in chains]
2613 chainVars = [np.var(data)
for data
in chains]
2615 BoverN = np.var(chainMeans)
2616 W = np.mean(chainVars)
2618 sigmaHat2 = W + BoverN
2622 VHat = sigmaHat2 + BoverN / m
2636 from lalinference
import bayespputils
as bppu
2644 for cfile
in chainfiles:
2645 if os.path.isfile(cfile):
2650 if mcmcChain
is None:
2651 return None,
None,
None,
None
2655 for j
in range(1, mcmcChain.shape[1]):
2656 neff, acl, acf = bppu.effectiveSampleSize(mcmcChain[:, j])
2657 neffstmp.append(neff)
2662 neffs.append(math.floor(np.mean(neffstmp)))
2665 nskip =
int(math.ceil(mcmcChain.shape[0] / np.mean(neffstmp)))
2668 mcmc.append(mcmcChain[::nskip, :])
2669 cl.append(mcmcChain.shape[0])
2671 print(
"File %s does not exist!" % cfile, file=sys.stderr)
2672 return None,
None,
None,
None
2676 cf = open(chainfiles[0],
"r")
2677 headers = cf.readline()
2678 headers = cf.readline()
2680 headers = re.sub(
"%",
"", headers)
2682 headers = re.sub(
"rads",
"", headers)
2684 headers = re.sub(
"[()]",
"", headers)
2688 for idx, parv
in enumerate(headers.split()):
2691 for j
in range(0, len(mcmc)):
2693 singlechain = achain[:, idx]
2694 lgr.append(singlechain)
2699 headers = headers.replace(
"\n",
"\tpost\n")
2702 comfile = chainfiles[0] +
"_common_tmp.dat"
2704 cf = open(comfile,
"w")
2706 print(
"Can't open common posterior file!", file=sys.stderr)
2711 for j
in range(0, narr.shape[0]):
2714 mline = np.append(mline, np.exp(mline[0]))
2716 strmline =
" ".join(
str(x)
for x
in mline) +
"\n"
2721 peparser = bppu.PEOutputParser(
"common")
2722 cf = open(comfile,
"r")
2723 commonResultsObj = peparser.parse(cf)
2730 pos = bppu.Posterior(commonResultsObj, SimInspiralTableEntry=
None, votfile=
None)
2735 cipos = bppu.PosteriorOneDPDF(
"cosiota", np.cos(pos[
"iota"].samples))
2743 return pos, neffs, grr, cl
2755 cfdata = np.loadtxt(cf, comments=
"%")
2761 headers = fc.readline()
2762 headers = fc.readline()
2765 headers = re.sub(
"%",
"", headers)
2767 headers = re.sub(
"rads",
"", headers)
2769 headers = re.sub(
"[()]",
"", headers)
2771 lh = len(headers.split())
2773 cfdata = np.array([])
2775 lines = cf.readlines()
2777 for i, line
in enumerate(lines):
2781 lvals = line.split()
2784 if len(lvals) != lh:
2789 lvalsf = map(float, lvals)
2795 cfdata = np.array(lvalsf)
2797 cfdata = np.vstack((cfdata, lvalsf))
2799 if cfdata.size == 0:
2809 This function will import a posterior sample file created by `lalapps_nest2pos` (or a nested
2810 sample file). It will be returned as a Posterior class object from `bayespputils`. The
2811 signal evidence and noise evidence are also returned.
2813 If there are samples in 'Q22' and not 'H0', and also a distance, or distance samples, then
2814 the Q22 samples will be converted into equivalent H0 values.
2818 postfile : str, required
2819 The file name of the posterior or nested sample file. In general this should be a HDF5 file with the extension
2820 '.hdf' or '.h5', although older format ascii files are still supported at the moment.
2821 nestedsamples : bool, default: False
2822 If the file being input contains nested samples, rather than posterior samples, then this flag must be set to
2824 removeuntrig : bool, default: True
2825 If this is True then any parameters that are sines or cosines of a value will have the value removed if present
2826 e.g. if cosiota and iota exist then iota will be removed.
2829 from lalinference
import bayespputils
as bppu
2832 fe = os.path.splitext(postfile)[-1].lower()
2835 if fe ==
".h5" or fe ==
".hdf":
2840 from lalinference
import LALInferenceHDF5NestedSamplesDatasetName
2842 samples = read_samples(
2843 postfile, tablename=LALInferenceHDF5NestedSamplesDatasetName
2845 params = samples.colnames
2848 for param
in params:
2849 replace_column(samples, param, samples[param].astype(float))
2853 samples.as_array().view(float).reshape(-1, len(params)),
2856 raise IOError(
"Could not import nested samples")
2858 peparser = bppu.PEOutputParser(
"hdf5")
2859 nsResultsObject = peparser.parse(postfile)
2863 peparser = bppu.PEOutputParser(
"common")
2864 nsResultsObject = peparser.parse(gzip.open(postfile,
"r"))
2866 peparser = bppu.PEOutputParser(
"common")
2867 nsResultsObject = peparser.parse(open(postfile,
"r"))
2869 pos = bppu.Posterior(nsResultsObject, SimInspiralTableEntry=
None)
2873 nsamps = len(pos[pnames[0]].samples)
2874 permarr = np.arange(nsamps)
2875 np.random.shuffle(permarr)
2878 for pname
in pnames:
2880 if pos[pname].samples.tolist().
count(pos[pname].samples[0]) == len(
2883 if pname ==
"f0_fixed":
2885 posfreqs = pos[pname].samples[0]
2886 elif pname ==
"dist":
2888 posdist = pos[pname].samples[0] / KPC
2894 shufpos = bppu.PosteriorOneDPDF(pname, pos[pname].samples[permarr])
2900 if "iota" in pos.names:
2901 posIota = pos[
"iota"].samples
2903 if posIota
is not None:
2905 cipos = bppu.PosteriorOneDPDF(
"cosiota", np.cos(posIota))
2912 if "i" in pos.names:
2913 posI = pos[
"i"].samples
2915 if posI
is not None:
2917 sinipos = bppu.PosteriorOneDPDF(
"sini", np.sin(posI))
2924 if "c21" in pos.names:
2925 posC21 = pos[
"c21"].samples
2928 if "c22" in pos.names:
2929 posC22 = pos[
"c22"].samples
2932 if "phi22" in pos.names:
2933 posphi22 = pos[
"phi22"].samples
2936 if posC22
is not None and posC21
is None:
2938 h0pos = bppu.PosteriorOneDPDF(
"h0", 2.0 * pos[
"c22"].samples)
2943 if posphi22
is not None:
2945 phi0pos = bppu.PosteriorOneDPDF(
2946 "phi0", np.fmod(pos[
"phi22"].samples + math.pi, 2.0 * math.pi)
2954 if "q22" in pos.names:
2955 posQ22 = pos[
"q22"].samples
2957 if "dist" in pos.names:
2959 pos[
"dist"].samples / KPC
2964 distpos = bppu.PosteriorOneDPDF(
"dist", posdist)
2967 if "f0" in pos.names:
2968 posfreqs = pos[
"f0"].samples
2971 if posQ22
is not None and posdist
is not None and "h0" not in pos.names:
2973 h0pos = bppu.PosteriorOneDPDF(
"h0",
quadrupole_to_h0(posQ22, posfreqs, posdist))
2978 if fe ==
".h5" or fe ==
".hdf":
2980 hdf = h5py.File(postfile,
"r")
2981 a = hdf[
"lalinference"][
"lalinference_nest"]
2982 sigev = a.attrs[
"log_evidence"]
2983 noiseev = a.attrs[
"log_noise_evidence"]
2986 B = np.loadtxt(postfile.replace(
".gz",
"") +
"_B.txt")
2991 return pos, sigev, noiseev
2996 return np.logaddexp(x, y)
3001 Perform trapezium rule integration for the logarithm of a function on a regular grid.
3005 lnf - a numpy array of values that are the natural logarithm of a function
3006 dx - a float giving the step size between values in the function
3010 The natural logarithm of the area under the function.
3013 return np.log(dx / 2.0) + logsumexp([logsumexp(lnf[:-1]), logsumexp(lnf[1:])])
3023 places = ranges[pname]
3024 axis = pnames.index(pname)
3028 x = np.apply_along_axis(logtrapz, axis, like, places[1] - places[0])
3029 elif len(places) == 1:
3032 q = np.arange(0, len(z)).astype(int) != axis
3033 newshape = tuple((np.array(list(z)))[q])
3034 x = np.reshape(like, newshape)
3038 pnames.remove(pname)
3044 def marginal(lnlike, pname, pnames, ranges, multflatprior=False):
3045 from copy
import deepcopy
3048 lnliketmp = deepcopy(lnlike)
3049 pnamestmp = deepcopy(pnames)
3050 rangestmp = deepcopy(ranges)
3057 + np.log(1.0 / (rangestmp[name][-1] - rangestmp[name][0])),
3063 lnliketmp =
marginalise(lnliketmp, name, pnamestmp, rangestmp)
3076 dt = np.min(np.diff(ts))
3085 chunkLengths.append(count)
3094 if t2 - t1 > 2.0 * dt
or count == chunkMax:
3095 chunkLengths.append(count)
3105 A function to estimate the signal-to-noise ratio of a pulsar signal (for a triaxial neutron
3106 star emitting at twice the rotation frequency from the l=m=2 quadrupole mode) in a given
3107 detector over a given time range, for a given one-sided power spectral density.
3111 source - a dictionary containing 'h0', 'cosiota', 'psi', 'ra' (rads or string), 'dec' (rads or string)
3112 det - the detector, e.g. 'H1'
3113 tstart - a GPS start time
3114 duration - a signal duration (seconds)
3115 Sn - a one-sided power spectral density
3116 dt - timestep for antenna response calculation [default: 1800s]
3120 rho - an estimate of the optimal matched filter SNR
3125 Ndays = duration / sidday
3126 Nsamps = np.floor(sidday / dt)
3127 tend = tstart + duration
3134 tts = np.arange(tstart, tstart + sidday, dt)
3137 Fp, Fc =
antenna_response(tts, source[
"ra"], source[
"dec"], source[
"psi"], det)
3139 sFp = np.floor(Ndays) * np.sum(Fp**2)
3140 sFc = np.floor(Ndays) * np.sum(Fc**2)
3143 if np.floor(Ndays) * sidday > dt:
3144 tts = np.arange(tstart + sidday * np.floor(Ndays), tend, dt)
3146 Fp, Fc =
antenna_response(tts, source[
"ra"], source[
"dec"], source[
"psi"], det)
3148 sFp += np.sum(Fp**2)
3149 sFc += np.sum(Fc**2)
3153 (source[
"h0"] ** 2 / Sn)
3156 sFp * (0.5 * (1 + source[
"cosiota"] ** 2)) ** 2
3157 + sFc * source[
"cosiota"] ** 2
3161 return np.sqrt(snr2)
3166 A function to estimate the signal amplitude of a pulsar signal (for a triaxial neutron star emitting at twice
3167 the rotation frequency from the l=m=2 quadrupole mode) for a given SNR in a particular detector over
3168 a given time range, for a given one-sided power spectral density.
3172 snr - the optimal matched filter signal-to-noise ratio
3173 source - a dictionary containing 'cosiota', 'psi', 'ra' (rads or string), 'dec' (rads or string)
3174 det - the detector, e.g. 'H1'
3175 tstart - a GPS start time for the signal
3176 duration - a signal duration (seconds)
3177 Sn - a one-sided power spectral density
3178 dt - timestep for antenna response calculation [default: 600s]
3182 h0 - an estimate of signal amplitude required to give the input SNR
3187 Ndays = duration / sidday
3188 Nsamps = np.floor(sidday / dt)
3189 tend = tstart + duration
3196 tts = np.arange(tstart, tstart + sidday, dt)
3199 Fp, Fc =
antenna_response(tts, source[
"ra"], source[
"dec"], source[
"psi"], det)
3201 sFp = np.floor(Ndays) * np.sum(Fp**2)
3202 sFc = np.floor(Ndays) * np.sum(Fc**2)
3205 if np.floor(Ndays) * sidday > dt:
3206 tts = np.arange(tstart + sidday * np.floor(Ndays), tend, dt)
3208 Fp, Fc =
antenna_response(tts, source[
"ra"], source[
"dec"], source[
"psi"], det)
3210 sFp += np.sum(Fp**2)
3211 sFc += np.sum(Fc**2)
3214 h02 = (snr**2 * Sn) / (
3217 sFp * (0.5 * (1 + source[
"cosiota"] ** 2)) ** 2
3218 + sFc * source[
"cosiota"] ** 2
3239 A function to calculate the 4-parameter posterior probability density for a continuous wave signal
3240 given a set of processed data from a set of detectors.
3244 dets - a list of strings containing the detectors being used in the likelihood calculation
3245 ts - a dictionary of 1d numpy arrays containing the GPS times of the data for each detector
3246 data - a dictionary of 1d numpy arrays containing the complex data values for each detector
3247 ra - the right ascension (in rads) of the source
3248 dec - the declination (in rads) of the source
3249 sigmas - a dictionary of 1d numpy arrays containing the times series of standard deviation
3250 values for each detector. If this is not given a Student's t likelihood will be used,
3251 but if it is given a Gaussian likelihood will be used (default: None)
3252 paramranges - a dictionary of tuples for each parameter ('h0', 'phi0', 'psi' and 'cosiota') giving
3253 the lower and upper ranges of the parameter grid and the number of grid points. If not
3254 given then defaults will be used.
3255 datachunks - in the calculation split the data into stationary chunks with a maximum length given
3256 by this value. If set to 0 or inf then the data will not be split. (default: 30)
3257 chunkmin - this is the shortest chunk length allowed to be included in the likelihood calculation
3259 ngrid - the number of grid points to use for each dimension of the likelihood calculation. This
3260 is used if the values are not specified in the paramranges argument (default: 50)
3261 outputlike - output the log likelihoods rather than posteriors (default: False)
3265 L - The 4d posterior (or likelihood) over all parameters
3266 h0pdf - The 1d marginal posterior for h0
3267 phi0pdf - The 1d marginal posterior for phi0 (the rotation frequency, not GW frequency)
3268 psipdf - The 1d marginal posterior for psi
3269 cosiotapdf - The 1d marginal posterior for cosiota
3270 lingrids - A dictionary of the grid points for each parameter
3271 sigev - The log evidence for the signal model
3272 noiseev - The log evidence for the noise model
3274 An example would be:
3278 # set the time series and data
3282 ts[det] = np.arange(900000000., 921000843., 60.)
3283 data[det] = np.random.randn(len(ts[det])) + 1j*np.random.randn(len(ts[det]))
3285 # set the parameter ranges
3289 paramranges['h0'] = (0., 2., 50)
3290 paramranges['psi'] = (0., np.pi/2., 50)
3291 paramranges['phi0'] = (0., np.pi, 50)
3292 paramranges['cosiota'] = (-1., 1., 50)
3294 L, h0pdf, phi0pdf, psipdf, cosiotapdf, grid, sigev, noiseev = pulsar_posterior_grid(dets, ts, data, ra, dec,
3295 paramranges=paramranges)
3301 from scipy.special
import gammaln
3305 liketype =
"studentst"
3307 liketype =
"gaussian"
3310 if not isinstance(dets, list):
3311 if isinstance(dets, string_types):
3314 print(
"Detector not, or incorrectly, set", file=sys.stderr)
3318 alloweddets = [
"H1",
"L1",
"V1",
"G1",
"T1",
"H2"]
3322 if det
not in alloweddets:
3324 "Detector not in list of allowed detectors ("
3325 +
",".join(alloweddets)
3333 print(
"No time stamps given for detector %s" % det, file=sys.stderr)
3338 print(
"No data time series given for detector %s" % det, file=sys.stderr)
3343 if det
not in sigmas:
3345 "No sigma time series given for detector %s" % det, file=sys.stderr
3350 if len(ts[det]) != len(data[det]):
3352 "Length of times stamps array and data array are inconsistent for %s"
3358 if len(ts[det]) != len(sigmas[det]):
3360 "Length of times stamps array and sigma array are inconsistent for %s"
3366 params = [
"h0",
"phi0",
"psi",
"cosiota"]
3368 defaultranges[
"phi0"] = (0.0, np.pi, ngrid)
3369 defaultranges[
"psi"] = (0.0, np.pi / 2.0, ngrid)
3370 defaultranges[
"cosiota"] = (-1.0, 1.0, ngrid)
3372 defaultranges[
"h0"] = (
3375 * np.max([np.std(data[dv])
for dv
in data])
3376 * np.sqrt(1440.0 / np.max([len(ts[dtx])
for dtx
in ts])),
3381 for param
in params:
3382 if param
in paramranges:
3383 if len(paramranges[param]) != 3:
3384 paramranges[param] = defaultranges[param]
3387 paramranges[param][1] < paramranges[param][0]
3388 or paramranges[param][2] < 1
3391 "Parameter ranges wrong for %s, reverting to defaults" % param,
3394 paramranges[param] = defaultranges[param]
3396 paramranges[param] = defaultranges[param]
3398 lingrids[param] = np.linspace(
3399 paramranges[param][0], paramranges[param][1], paramranges[param][2]
3401 shapetuple += (paramranges[param][2],)
3404 [H0S, PHI0S, PSIS, COSIS] = np.meshgrid(
3408 lingrids[
"cosiota"],
3412 APLUS = 0.5 * (1.0 + COSIS**2)
3414 SINPHI = np.sin(2.0 * PHI0S)
3415 COSPHI = np.cos(2.0 * PHI0S)
3416 SIN2PSI = np.sin(2.0 * PSIS)
3417 COS2PSI = np.cos(2.0 * PSIS)
3419 Apsinphi_2 = 0.5 * APLUS * SINPHI
3420 Acsinphi_2 = 0.5 * ACROSS * SINPHI
3421 Apcosphi_2 = 0.5 * APLUS * COSPHI
3422 Accosphi_2 = 0.5 * ACROSS * COSPHI
3423 Acpsphicphi_2 = 0.5 * APLUS * ACROSS * COSPHI * SINPHI
3427 TWOAPACSPCP = 2.0 * APLUS * ACROSS * SINPHI * COSPHI
3431 C2PS2P = COS2PSI * SIN2PSI
3438 like = np.zeros(shapetuple)
3444 logprior -= np.log(paramranges[p][1] - paramranges[p][0])
3448 if liketype ==
"gaussian":
3451 nstd = np.ones(len(ts[det]))
3457 if np.isfinite(datachunks)
and datachunks > 0:
3460 chunklengths = [len(ts[det])]
3463 for cl
in chunklengths:
3464 endidx = startidx + cl
3470 thisdata = data[det][startidx:endidx] / nstd[startidx:endidx]
3471 ast = As[startidx:endidx] / nstd[startidx:endidx]
3472 bst = Bs[startidx:endidx] / nstd[startidx:endidx]
3476 AB = np.dot(ast, bst)
3478 dA1real = np.dot(thisdata.real, ast)
3479 dA1imag = np.dot(thisdata.imag, ast)
3481 dB1real = np.dot(thisdata.real, bst)
3482 dB1imag = np.dot(thisdata.imag, bst)
3484 dd1real = np.sum(thisdata.real**2)
3485 dd1imag = np.sum(thisdata.imag**2)
3487 P1 = A * C2P2 + B * S2P2 + 2.0 * AB * C2PS2P
3488 P2 = B * C2P2 + A * S2P2 - 2.0 * AB * C2PS2P
3489 P3 = AB * (C2P2 - S2P2) - A * C2PS2P + B * C2PS2P
3491 hr2 = (Apcosphi_2**2) * P1 + (Acsinphi_2**2) * P2 + Acpsphicphi_2 * P3
3492 hi2 = (Apsinphi_2**2) * P1 + (Accosphi_2**2) * P2 - Acpsphicphi_2 * P3
3494 d1hr = Apcosphi_2 * (dA1real * COS2PSI + dB1real * SIN2PSI) + Acsinphi_2 * (
3495 dB1real * COS2PSI - dA1real * SIN2PSI
3497 d1hi = Apsinphi_2 * (dA1imag * COS2PSI + dB1imag * SIN2PSI) - Accosphi_2 * (
3498 dB1imag * COS2PSI - dA1imag * SIN2PSI
3501 chiSq = dd1real + dd1imag + H02 * (hr2 + hi2) - 2.0 * H0S * (d1hr + d1hi)
3503 if liketype ==
"gaussian":
3504 like -= 0.5 * (chiSq)
3505 like -= (cl * np.log(2.0 * np.pi)) + 2.0 * np.sum(
3506 np.log(nstd[startidx:endidx])
3508 noiselike -= 0.5 * (dd1real + dd1imag)
3509 noiselike -= (cl * np.log(2.0 * np.pi)) + 2.0 * np.sum(
3510 np.log(nstd[startidx:endidx])
3513 like -= float(cl) * np.log(chiSq)
3514 like += gammaln(cl) - np.log(2.0) - cl * np.log(np.pi)
3515 noiselike -= float(cl) * np.log(dd1real + dd1imag)
3516 noiselike += gammaln(cl) - np.log(2.0) - cl * np.log(np.pi)
3524 sigev =
marginal(like,
"all", params, lingrids)
3536 posts[p] = np.exp(
marginal(like, p, params, lingrids))
3539 like, p, params, lingrids, multflatprior=
True
3543 evrat = sigev - noiselike
3558 ATNF_VERSION =
"1.58"
3563 Get the pulsar (psr) distance (DIST in kpc), proper motion corrected period derivative (P1_I) and any association
3564 (ASSOC e.g. GC) from the ATNF catalogue.
3569 psrname = re.sub(
r"\+",
"%2B", psr)
3572 "http://www.atnf.csiro.au/people/pulsar/psrcat/proc_form.php?version="
3575 atnfurl +=
"&Dist=Dist&Assoc=Assoc&P1_i=P1_i"
3577 "&startUserDefined=true&c1_val=&c2_val=&c3_val=&c4_val=&sort_attr=jname&sort_order=asc&condition=&pulsar_names="
3580 atnfurl +=
"&ephemeris=selected&submit_ephemeris=Get+Ephemeris&coords_unit=raj%2Fdecj&radius=&coords_1=&coords_2="
3581 atnfurl +=
"&style=Long+with+last+digit+error&no_value=*&fsize=3&x_axis=&x_scale=linear&y_axis=&y_scale=linear&state=query"
3584 urldat = requests.get(atnfurl)
3586 r"<pre[^>]*>([^<]+)</pre>",
str(urldat.content)
3589 predat.group(1).strip().split(
r"\n")
3593 "Warning... could not get information from ATNF pulsar catalogue.",
3602 for line
in [p
for p
in pdat
if len(p) != 0]:
3603 if "WARNING" in line
or "not in catalogue" in line:
3606 if "DIST" in vals[0]:
3607 dist = float(vals[1])
3608 if "P1_I" in vals[0]:
3609 age = float(vals[1])
3610 if "ASSOC" in vals[0]:
3613 return (dist, p1_I, assoc, atnfurl)
3618 Convert a covariance matrix to a correlation coefficient matrix.
3619 Return the correlation coefficient matrix and the standard deviations
3620 from the covariance matrix.
3624 sigmas = np.sqrt(np.diag(cov))
3627 D = sigmas * np.identity(cov.shape[0])
3630 Dinv = np.linalg.inv(D)
3633 Ccor = np.dot(np.dot(Dinv, cov), Dinv)
def __init__(self, parfilenm)
This class parses a TEMPO(2)-style pulsar parameter file.
def __getitem__(self, key)
def __init__(self, priorfilenm)
def __getitem__(self, key)
def heterodyned_pulsar_signal(pardict, detector, starttime=900000000.0, duration=86400.0, dt=60.0, datatimes=None)
def quadrupole_to_h0(q22, freq, dist)
def rad_to_hms(rad)
rad_to_hms(rad): Convert radians to hours, minutes, and seconds of arc.
def dec_to_rad(dec_string)
dec_to_rad(dec_string): Given a string containing DEC information as 'dd:mm:ss.ssss',...
def hms_to_rad(hour, mins, sec)
hms_to_rad(hour, min, sec): Convert hours, minutes, and seconds of arc to radians
def cov_to_cor(cov)
Convert a covariance matrix to a correlation coefficient matrix.
def pulsar_mcmc_to_posterior(chainfiles)
def heterodyned_pinsf_pulsar(starttime, duration, dt, detector, pardict)
def plot_posterior_hist(poslist, param, ifos, parambounds=[float("-inf"), float("inf")], nbins=50, upperlimit=0, overplot=False, parfile=None, mplparams=False)
def dms_to_rad(deg, mins, sec)
dms_to_rad(deg, min, sec): Convert degrees, minutes, and seconds of arc to radians.
def pulsar_posterior_grid(dets, ts, data, ra, dec, sigmas=None, paramranges={}, datachunks=30, chunkmin=5, ngrid=50, outputlike=False)
A function to calculate the 4-parameter posterior probability density for a continuous wave signal gi...
def rad_to_string(rad, ra_or_dec)
rad_to_string(rad, ra_or_dec): Convert an angle in radians to hours/degrees, minutes seconds and outp...
def ra_to_rad(ra_string)
ra_to_rad(ar_string): Given a string containing RA information as 'hh:mm:ss.ssss',...
def pferrs(porf, porferr, pdorfd=None, pdorfderr=None)
pferrs(porf, porferr, pdorfd=None, pdorfderr=None): Calculate the period or frequency errors and the ...
def plot_limits_hist(lims, param, ifos, prevlims=None, bins=20, overplot=False, mplparams=False)
def marginalise(like, pname, pnames, ranges)
def pulsar_estimate_snr(source, det, tstart, duration, Sn, dt=1800)
A function to estimate the signal-to-noise ratio of a pulsar signal (for a triaxial neutron star emit...
def read_pulsar_mcmc_file(cf)
def heterodyned_triaxial_pulsar(starttime, duration, dt, detector, pardict)
def convert_model_parameters(pardict)
def coord_to_string(h_or_d, m, s)
coord_to_string(h_or_d, m, s): Return a formatted string of RA or DEC values as 'hh:mm:ss....
def logtrapz(lnf, dx)
Perform trapezium rule integration for the logarithm of a function on a regular grid.
def gelman_rubins(chains)
def pulsar_nest_to_posterior(postfile, nestedsamples=False, removeuntrig=True)
This function will import a posterior sample file created by lalapps_nest2pos (or a nested sample fil...
def phipsiconvert(phipchain, psipchain)
def upper_limit_greedy(pos, upperlimit=0.95, nbins=100)
def marginal(lnlike, pname, pnames, ranges, multflatprior=False)
def pulsar_estimate_h0_from_snr(snr, source, det, tstart, duration, Sn, dt=600)
A function to estimate the signal amplitude of a pulsar signal (for a triaxial neutron star emitting ...
def plot_posterior_hist2D(poslist, params, ifos, bounds=None, nbins=[50, 50], parfile=None, overplot=False, mplparams=False)
def get_chunk_lengths(ts, chunkMax)
def tukey_window(N, alpha=0.5)
def read_hist_from_file(histfile)
def h0_to_quadrupole(h0, freq, dist)
def get_optimal_snr(s, sig)
def hist_norm_bounds(samples, nbins, low=float("-inf"), high=float("inf"))
def plot_2Dhist_from_file(histfile, ndimlabel, mdimlabel, margpars=True, mplparams=False)
def spin_down_limit(freq, fdot, dist)
def inject_pulsar_signal(starttime, duration, dt, detectors, pardict, freqfac=[2.0], npsds=None, snrscale=None)
def antenna_response(gpsTime, ra, dec, psi, det)
def upper_limit(pos, upperlimit=0.95, parambounds=[float("-inf"), float("inf")], nbins=50)
def detector_noise(det, f)
def plot_posterior_chain(poslist, param, ifos, grr=None, withhist=0, mplparams=False)
def get_atnf_info(psr)
Get the pulsar (psr) distance (DIST in kpc), proper motion corrected period derivative (P1_I) and any...
def plot_h0_lims(h0lims, f0gw, ifos, xlims=[10, 1500], ulesttop=None, ulestbot=None, prevlim=None, prevlimf0gw=None, overplot=False, mplparams=False)
def dms_to_deg(deg, mins, sec)
dms_to_deg(deg, min, sec): Convert degrees, minutes, and seconds of arc to degrees.
def rad_to_dms(rad)
rad_to_dms(rad): Convert radians to degrees, minutes, and seconds of arc.
def plot_Bks_ASDs(Bkdata, delt=86400, plotpsds=True, plotfscan=False, removeoutlier=None, mplparams=False)
def p_to_f(p, pd, pdd=None)
p_to_f(p, pd, pdd=None): Convert period, period derivative and period second derivative to the equiva...
def h0ul_from_prior_file(priorfile, ulval=0.95)
def h0_to_ellipticity(h0, freq, dist)