Deploying to gh-pages from @ a6d36d6 🚀

latest/examples/catalog/data.py

@@ -0,0 +1,25 @@
+import matplotlib.pyplot as pp
+import pycbc.catalog
+
+
+m = pycbc.catalog.Merger("GW170817", source='gwtc-1')
+
+fig, axs = pp.subplots(2, 1, sharex=True, sharey=True)
+for ifo, ax in zip(["L1", "H1"], axs):
+    pp.sca(ax)
+    pp.title(ifo)
+    # Retreive data around the BNS merger
+    ts = m.strain(ifo).time_slice(m.time - 15, m.time + 6)
+
+    # Whiten the data with a 4s filter
+    white = ts.whiten(4, 4)
+
+    times, freqs, power = white.qtransform(.01, logfsteps=200,
+                                        qrange=(110, 110),
+                                        frange=(20, 512))
+    pp.pcolormesh(times, freqs, power**0.5, vmax=5)
+
+pp.yscale('log')
+pp.ylabel("Frequency (Hz)")
+pp.xlabel("Time (s)")
+pp.show()

latest/examples/catalog/stat.py

@@ -0,0 +1,13 @@
+import matplotlib.pyplot as pp
+import pycbc.catalog
+
+
+c = pycbc.catalog.Catalog(source='gwtc-2')
+mchirp, elow, ehigh = c.median1d('mchirp', return_errors=True)
+spin = c.median1d('chi_eff')
+
+pp.errorbar(mchirp, spin, xerr=[-elow, ehigh], fmt='o', markersize=7)
+pp.xlabel('Chirp Mass')
+pp.xscale('log')
+pp.ylabel('Effective Spin')
+pp.show()

latest/examples/dataquality/hwinj.py

@@ -0,0 +1,22 @@
+"""This example shows how to determine when a CBC hardware injection is present
+in the data from a detector.
+"""
+
+import matplotlib.pyplot as pp
+from pycbc import dq
+
+
+start_time = 1126051217
+end_time = start_time + 10000000
+
+# Get times that the Livingston detector has CBC injections into the data
+segs = dq.query_flag('L1', 'CBC_HW_INJ', start_time, end_time)
+
+pp.figure(figsize=[10, 2])
+for seg in segs:
+    start, end = seg
+    pp.axvspan(start, end, color='blue')
+
+pp.xlabel('Time (s)')
+pp.show()
+

latest/examples/dataquality/on.py

@@ -0,0 +1,27 @@
+"""This example shows how to determine when a detector is active."""
+
+import matplotlib.pyplot as pp
+from pycbc import dq
+from pycbc.results import ifo_color
+
+
+start_time = 1126051217
+end_time = start_time + 100000
+
+# Get times that the Hanford detector has data
+hsegs = dq.query_flag('H1', 'DATA', start_time, end_time)
+
+# Get times that the Livingston detector has data
+lsegs = dq.query_flag('L1', 'DATA', start_time, end_time)
+
+pp.figure(figsize=[10,2])
+for seg in lsegs:
+    start, end = seg
+    pp.axvspan(start, end, color=ifo_color('L1'), ymin=0.1, ymax=0.4)
+
+for seg in hsegs:
+    start, end = seg
+    pp.axvspan(start, end, color=ifo_color('H1'), ymin=0.6, ymax=0.9)
+
+pp.xlabel('Time (s)')
+pp.show()

latest/examples/distributions/mass_examples.py

@@ -0,0 +1,44 @@
+import matplotlib.pyplot as plt
+from pycbc import distributions
+
+# Create a mass distribution object that is uniform between 0.5 and 1.5
+# solar masses.
+mass1_distribution = distributions.Uniform(mass1=(0.5, 1.5))
+# Take 100000 random variable samples from this uniform mass distribution.
+mass1_samples = mass1_distribution.rvs(size=1000000)
+
+# Draw another distribution that is Gaussian between 0.5 and 1.5 solar masses
+# with a mean of 1.2 solar masses and a standard deviation of 0.15 solar
+# masses. Gaussian takes the variance as an input so square the standard
+# deviation.
+variance = 0.15*0.15
+mass2_gaussian = distributions.Gaussian(mass2=(0.5, 1.5), mass2_mean=1.2,
+                                        mass2_var=variance)
+
+# Take 100000 random variable samples from this gaussian mass distribution.
+mass2_samples = mass2_gaussian.rvs(size=1000000)
+
+# We can make pairs of distributions together, instead of apart.
+two_mass_distributions = distributions.Uniform(mass3=(1.6, 3.0),
+                                               mass4=(1.6, 3.0))
+two_mass_samples = two_mass_distributions.rvs(size=1000000)
+
+# Choose 50 bins for the histogram subplots.
+n_bins = 50
+
+# Plot histograms of samples in subplots
+fig, axes = plt.subplots(nrows=2, ncols=2)
+ax0, ax1, ax2, ax3, = axes.flat
+
+ax0.hist(mass1_samples['mass1'], bins = n_bins)
+ax1.hist(mass2_samples['mass2'], bins = n_bins)
+ax2.hist(two_mass_samples['mass3'], bins = n_bins)
+ax3.hist(two_mass_samples['mass4'], bins = n_bins)
+
+ax0.set_title('Mass 1 samples')
+ax1.set_title('Mass 2 samples')
+ax2.set_title('Mass 3 samples')
+ax3.set_title('Mass 4 samples')
+
+plt.tight_layout()
+plt.show()

latest/examples/distributions/mchirp_q_from_uniform_m1m2_example.py

@@ -0,0 +1,64 @@
+import matplotlib.pyplot as plt
+from pycbc import distributions
+from pycbc import conversions
+import numpy as np
+
+# Create chirp mass and mass ratio distribution object that is uniform
+# in mass1 and mass2
+minmc = 5
+maxmc = 60
+mc_distribution = distributions.MchirpfromUniformMass1Mass2(mc=(minmc,maxmc))
+# generate q in a symmetric range [min, 1/min] to make mass1 and mass2
+# symmetric
+minq = 1/4
+maxq = 1/minq
+q_distribution = distributions.QfromUniformMass1Mass2(q=(minq,maxq))
+
+# Take 100000 random variable samples from this chirp mass and mass ratio
+# distribution.
+n_size = 100000
+mc_samples = mc_distribution.rvs(size=n_size)
+q_samples = q_distribution.rvs(size=n_size)
+
+# Convert chirp mass and mass ratio to mass1 and mass2
+m1 = conversions.mass1_from_mchirp_q(mc_samples['mc'],q_samples['q'])
+m2 = conversions.mass2_from_mchirp_q(mc_samples['mc'],q_samples['q'])
+
+# Check the 1D marginalization of mchirp and q is consistent with the 
+# expected analytical formula
+n_bins = 200
+xq = np.linspace(minq,maxq,100)
+yq = ((1+xq)/(xq**3))**(2/5)
+xmc = np.linspace(minmc,maxmc,100)
+ymc = xmc
+
+plt.figure(figsize=(10,10))
+# Plot histograms of samples in subplots
+plt.subplot(221)
+plt.hist2d(mc_samples['mc'], q_samples['q'], bins=n_bins, cmap='Blues')
+plt.xlabel('chirp mass')
+plt.ylabel('mass ratio')
+plt.colorbar(fraction=.05, pad=0.05,label='number of samples')
+
+plt.subplot(222)
+plt.hist2d(m1, m2, bins=n_bins, cmap='Blues')
+plt.xlabel('mass1')
+plt.ylabel('mass2')
+plt.colorbar(fraction=.05, pad=0.05,label='number of samples')
+
+plt.subplot(223)
+plt.hist(mc_samples['mc'],density=True,bins=100,label='samples')
+plt.plot(xmc,ymc*mc_distribution.norm,label='$P(M_c)\propto M_c$')
+plt.xlabel('chirp mass')
+plt.ylabel('PDF')
+plt.legend()
+
+plt.subplot(224)
+plt.hist(q_samples['q'],density=True,bins=n_bins,label='samples')
+plt.plot(xq,yq*q_distribution.norm,label='$P(q)\propto((1+q)/q^3)^{2/5}$')
+plt.xlabel('mass ratio')
+plt.ylabel('PDF')
+plt.legend()
+
+plt.tight_layout()
+plt.show()

latest/examples/distributions/sampling_from_config_example.py

@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from pycbc.distributions.utils import draw_samples_from_config
+
+
+# A path to the .ini file.
+CONFIG_PATH = "./pycbc_bbh_prior.ini"
+random_seed = np.random.randint(low=0, high=2**32-1)
+
+# Draw a single sample.
+sample = draw_samples_from_config(
+            path=CONFIG_PATH, num=1, seed=random_seed)
+
+# Print all parameters.
+print(sample.fieldnames)
+print(sample)
+# Print a certain parameter, for example 'mass1'.
+print(sample[0]['mass1'])
+
+# Draw 1000000 samples, and select all values of a certain parameter.
+n_bins = 50
+samples = draw_samples_from_config(
+            path=CONFIG_PATH, num=1000000, seed=random_seed)
+
+fig, axes = plt.subplots(nrows=3, ncols=2)
+ax1, ax2, ax3, ax4, ax5, ax6 = axes.flat
+
+ax1.hist(samples[:]['srcmass1'], bins=n_bins)
+ax2.hist(samples[:]['srcmass2'], bins=n_bins)
+ax3.hist(samples[:]['comoving_volume'], bins=n_bins)
+ax4.hist(samples[:]['redshift'], bins=n_bins)
+ax5.hist(samples[:]['distance'], bins=n_bins)
+ax6.hist(samples[:]['mass1'], bins=n_bins)
+
+ax1.set_title('srcmass1')
+ax2.set_title('srcmass2')
+ax3.set_title('comoving_volume')
+ax4.set_title('redshift')
+ax5.set_title('distance')
+ax6.set_title('mass1 or mass2')
+
+plt.tight_layout()
+plt.show()

latest/examples/distributions/spin_examples.py

@@ -0,0 +1,56 @@
+import matplotlib.pyplot as plt
+import numpy
+import pycbc.coordinates as co
+from pycbc import distributions
+
+# We can choose any bounds between 0 and pi for this distribution but in
+# units of pi so we use between 0 and 1
+theta_low = 0.
+theta_high = 1.
+
+# Units of pi for the bounds of the azimuthal angle which goes from 0 to 2 pi
+phi_low = 0.
+phi_high = 2.
+
+# Create a distribution object from distributions.py. Here we are using the
+# Uniform Solid Angle function which takes
+# theta = polar_bounds(theta_lower_bound to a theta_upper_bound), and then
+# phi = azimuthal_ bound(phi_lower_bound to a phi_upper_bound).
+uniform_solid_angle_distribution = distributions.UniformSolidAngle(
+                                          polar_bounds=(theta_low,theta_high),
+                                          azimuthal_bounds=(phi_low,phi_high))
+
+# Now we can take a random variable sample from that distribution. In this
+# case we want 50000 samples.
+solid_angle_samples = uniform_solid_angle_distribution.rvs(size=500000)
+
+# Make spins with unit length for coordinate transformation below.
+spin_mag = numpy.ndarray(shape=(500000), dtype=float)
+
+for i in range(0,500000):
+    spin_mag[i] = 1.
+
+# Use the pycbc.coordinates as co spherical_to_cartesian function to convert
+# from spherical polar coordinates to cartesian coordinates
+spinx, spiny, spinz = co.spherical_to_cartesian(spin_mag,
+                                                solid_angle_samples['phi'],
+                                                solid_angle_samples['theta'])
+
+# Choose 50 bins for the histograms.
+n_bins = 50
+
+plt.figure(figsize=(10,10))
+plt.subplot(2, 2, 1)
+plt.hist(spinx, bins = n_bins)
+plt.title('Spin x samples')
+
+plt.subplot(2, 2, 2)
+plt.hist(spiny, bins = n_bins)
+plt.title('Spin y samples')
+
+plt.subplot(2, 2, 3)
+plt.hist(spinz, bins = n_bins)
+plt.title('Spin z samples')
+
+plt.tight_layout()
+plt.show()

latest/examples/distributions/spin_spatial_distr_example.py

@@ -0,0 +1,50 @@
+import numpy
+import matplotlib.pyplot as plt
+import pycbc.coordinates as co
+from pycbc import distributions
+
+# We can choose any bounds between 0 and pi for this distribution but in units
+# of pi so we use between 0 and 1.
+theta_low = 0.
+theta_high = 1.
+
+# Units of pi for the bounds of the azimuthal angle which goes from 0 to 2 pi.
+phi_low = 0.
+phi_high = 2.
+
+# Create a distribution object from distributions.py
+# Here we are using the Uniform Solid Angle function which takes
+# theta = polar_bounds(theta_lower_bound to a theta_upper_bound), and then
+# phi = azimuthal_bound(phi_lower_bound to a phi_upper_bound).
+uniform_solid_angle_distribution = distributions.UniformSolidAngle(
+                                          polar_bounds=(theta_low,theta_high),
+                                          azimuthal_bounds=(phi_low,phi_high))
+
+# Now we can take a random variable sample from that distribution.
+# In this case we want 50000 samples.
+solid_angle_samples = uniform_solid_angle_distribution.rvs(size=10000)
+
+# Make a spin 1 magnitude since solid angle is only 2 dimensions and we need a
+# 3rd dimension for a 3D plot that we make later on.
+spin_mag = numpy.ndarray(shape=(10000), dtype=float)
+
+for i in range(0,10000):
+    spin_mag[i] = 1.
+
+# Use pycbc.coordinates as co. Use  spherical_to_cartesian function to
+# convert from spherical polar coordinates to cartesian coordinates.
+spinx, spiny, spinz = co.spherical_to_cartesian(spin_mag,
+                                                solid_angle_samples['phi'],
+                                                solid_angle_samples['theta'])
+
+# Plot the spherical distribution of spins to make sure that we
+# distributed across  the surface of a sphere.
+
+fig = plt.figure(figsize=(10,10))
+ax = fig.add_subplot(111, projection='3d')
+ax.scatter(spinx, spiny, spinz, s=1)
+
+ax.set_xlabel('Spin X Axis')
+ax.set_ylabel('Spin Y Axis')
+ax.set_zlabel('Spin Z Axis')
+plt.show()

latest/examples/filter/chisq.py

@@ -0,0 +1,41 @@
+"""This example shows how to calculate the chi^2 discriminator described in
+https://arxiv.org/abs/gr-qc/0405045, also known as the "power chi^2" or "Allen
+chi^2" discriminator.
+"""
+
+import matplotlib.pyplot as pp
+import pycbc.noise
+import pycbc.psd
+import pycbc.waveform
+import pycbc.vetoes
+
+
+# Generate some noise with an advanced ligo psd
+flow = 30.0
+delta_f = 1.0 / 16
+flen = int(2048 / delta_f) + 1
+psd = pycbc.psd.aLIGOZeroDetHighPower(flen, delta_f, flow)
+
+# Generate 16 seconds of noise at 4096 Hz
+delta_t = 1.0 / 4096
+tsamples = int(16 / delta_t)
+strain = pycbc.noise.noise_from_psd(tsamples, delta_t, psd, seed=127)
+stilde = strain.to_frequencyseries()
+
+# Calculate the power chisq time series
+hp, hc = pycbc.waveform.get_fd_waveform(approximant='IMRPhenomD',
+                             mass1=25, mass2=25,
+                             f_lower=flow, delta_f=stilde.delta_f)
+
+hp.resize(len(stilde))
+num_bins = 16
+chisq = pycbc.vetoes.power_chisq(hp, stilde, num_bins, psd,
+                                      low_frequency_cutoff=flow)
+
+# convert to a reduced chisq
+chisq /= (num_bins * 2) - 2
+
+pp.plot(chisq.sample_times, chisq)
+pp.ylabel('$\chi^2_r$')
+pp.xlabel('time (s)')
+pp.show()