Commit c9ca9111 authored by Richard O'Shaughnessy's avatar Richard O'Shaughnessy

Merge branch 'temp-RIT-Tides-port_master-GPUIntegration' into temp-RIT-Tides

parents 43c5ae68 b588e0df
setup:
# Retrieve Pankow 'blind injection challenge'
gsiscp -r ldas-grid.ligo.caltech.edu:/home/pankow/richard_mdc/signal_hoft .
gsiscp -r ldas-grid.ligo.caltech.edu:/home/pankow/richard_mdc/signal_hoft_no_noise .
gsiscp ldas-grid.ligo.caltech.edu:/home/pankow/richard_mdc/mdc.xml.gz .
find -L signal_hoft -name '*.gwf' -print | lalapps_path2cache > test1.cache
find -L signal_hoft_no_noise -name '*.gwf' -print | lalapps_path2cache > test1-0noise.cache
build_gpu:
# Default configuration
# - cuda10
# - tesla P40/P100 [head nodes] or K10 https://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/index.html#virtual-architecture-feature-list
# -arch=sm_60 # ldas-pcdev13 architecture (GTX 10 series also pascal-powered, should work for all)
# Check configuraito with
# nvidia-smi # prints available machines
nvcc -arch=sm_60 -cubin -o cuda_Q_inner_product.cubin cuda_Q_inner_product.cu
# iLIGO noise
test-fiducial-data:
python test_like_and_samp.py --inj-xml mdc.xml.gz --cache-file test1-0noise.cache --channel-name H1=FAKE-STRAIN --channel-name L1=FAKE-STRAIN --channel-name V1=FAKE-STRAIN --force-read-full-frames --show-psd --show-likelihood-versus-time --show-sampler-inputs --show-sampler-results --Niter 1000000 --save-sampler-file fiducial-0noise --save-threshold-fraction 0.95
# iLIGO noise
test-fiducial-data-noisy:
python test_like_and_samp.py --inj-xml mdc.xml.gz --cache-file test1.cache --channel-name H1=FAKE-STRAIN --channel-name L1=FAKE-STRAIN --channel-name V1=FAKE-STRAIN --force-read-full-frames --show-psd --show-likelihood-versus-time --show-sampler-inputs --show-sampler-results --Niter 1000000 --save-sampler-file fiducial-0noise --save-threshold-fraction 0.95
# AVAILABLE GPU MACHINES (2019-01 edition)
# CIT:
# 560 nodes with single GTX 1050 Ti
# 32 nodes with dual K10 PCI cards (shows up as 4 GPU devices) and dual GTX 750 Ti
# 16 nodes with dual GTX 750 Ti (and dual Xeon Phi co-processor)
# LHO and LLO each:
# 40 nodes with GTX 1050 Ti
# single DGX-1 system with 8 tightly coupled V100
# And all 3 sites advertise the login/head node GPU devices in their MOTD.
# In general, you can use condor_status to query all the machine in a local Condor pool to see what resources are being advertised as available to run jobs. This includes information about GPUs when available. For example, the CIT cluster is currently advertising 872 Condor slots with CUDA 9.2,
# ldas-grid:~> condor_status -af CUDADeviceName -constraint "DynamicSlot =!= True" | sort | uniq -c
# 549 GeForce GTX 1050 Ti
# 16 GeForce GTX 750 Ti
# 1911 undefined
# ldas-grid:~> condor_status -af CUDA0DeviceName -constraint "DynamicSlot =!= True" | sort | uniq -c
# 30 GeForce GTX 750 Ti
# 2 Tesla K10.G2.8GB
# 1 Tesla K80
# 2442 undefined
test:
# test hlm generation, overlaps
......@@ -34,10 +53,3 @@ test:
injection-info:
# Parameters
ligolw_print mdc.xml.gz -t sim_inspiral -c mass1 -c mass2 -c longitude -c latitude
# Frame files
FrDump -i signal_hoft/H-ER_signal_hoft-10000/H-ER_signal_hoft-1000000000-64.gwf
FrDump -i signal_hoft/L-ER_signal_hoft-10000/L-ER_signal_hoft-1000000000-64.gwf
FrDump -i signal_hoft/V-ER_signal_hoft-10000/V-ER_signal_hoft-1000000000-64.gwf
from __future__ import division
import cupy
def Q_inner_product_cupy(Q, A, start_indices, window_size):
num_time_points, num_lms = Q.shape
num_extrinsic_samples, _ = A.shape
assert not cupy.isfortran(Q)
assert not cupy.isfortran(A)
out = cupy.empty(
(num_extrinsic_samples, window_size),
dtype=cupy.complex128,
order="C",
)
mod = cupy.cuda.function.Module()
mod.load_file("cuda_Q_inner_product.cubin")
Q_prod_fn = mod.get_function("Q_inner")
float_prec = 16
num_threads_x = 4
num_threads_y = 1024 // 4
block_size = num_threads_x, num_threads_y, 0
grid_size = (
(num_extrinsic_samples+num_threads_x-1)//num_threads_x,
0,
0,
)
args = (
Q, A, start_indices, window_size,
num_time_points, num_extrinsic_samples, num_lms,
out,
)
Q_prod_fn(
grid_size, block_size, args,
shared_mem=cupy.int32(num_threads_x*num_lms*float_prec),
)
return out
#include <cuComplex.h>
extern "C" {
__global__ void Q_inner(
cuDoubleComplex * Q, cuDoubleComplex * A,
int * index_start,
int window_size,
int num_time_points,
int num_extrinsic_samples,
int num_lms,
cuDoubleComplex * out
){
extern __shared__ cuDoubleComplex A_sample[];
/* Figure out which extrinsic sample number we're on. */
size_t sample_idx = threadIdx.x + blockDim.x*blockIdx.x;
// time index in the window for each sample
size_t t_idx = threadIdx.y + blockDim.y * blockIdx.y;
/* Only do something if we're not out of bounds. */
if (sample_idx < num_extrinsic_samples) {
for (size_t i = 0; i<num_lms; ++i) {
A_sample[threadIdx.x*num_lms+i] = A[sample_idx*num_lms+i];
}
__syncthreads();
/* Determine the time index we need to use. */
size_t i_first_time = index_start[sample_idx];
/* Iterate over the time window. */
for (size_t i_time = t_idx; i_time < window_size; i_time+=blockDim.y) {
/* Determine the index we're going to output to. */
size_t i_output = sample_idx*window_size + i_time;
/* Ensure the output is initialized to zero. */
cuDoubleComplex out_tmp = make_cuDoubleComplex(0.0, 0.0);
/* Take the outer product over the lm axis. */
for (size_t i_lm = 0; i_lm < num_lms; ++i_lm) {
out_tmp = cuCadd(out_tmp, cuCmul(
A_sample[threadIdx.x*num_lms + i_lm],
Q[(i_first_time+i_time)*num_lms + i_lm]
));
}
out[i_output] = out_tmp;
}
}
}
}
......@@ -12,12 +12,26 @@ import functools
from optparse import OptionParser, OptionGroup
import numpy
import numpy as np
try:
import cupy
xpy_default=cupy
identity_convert = cupy.asnumpy
except:
print ' no cupy'
import numpy as cupy # will automatically replace cupy calls with numpy!
optimized_gpu_tools=None
Q_inner_product=None
xpy_default=numpy # just in case, to make replacement clear and to enable override
identity_convert = lambda x: x # trivial return itself
import lal
from glue.ligolw import utils, lsctables, table, ligolw
from glue.ligolw.utils import process
import glue.lal
import pylal
#import pylal
import lalsimutils
import factored_likelihood
......@@ -91,11 +105,14 @@ optp.add_option("-t", "--event-time", type=float, help="GPS time of the event --
optp.add_option("-i", "--inv-spec-trunc-time", type=float, default=8., help="Timescale of inverse spectrum truncation in seconds (Default is 8 - give 0 for no truncation)")
optp.add_option("-w", "--window-shape", type=float, default=0, help="Shape of Tukey window to apply to data (default is no windowing)")
optp.add_option("-m", "--time-marginalization", action="store_true", help="Perform marginalization over time via direct numerical integration. Default is false.")
optp.add_option("--vectorized", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, not LAL data structures. (Combine with --gpu to enable GPU use, where available)")
optp.add_option("--gpu", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, CONVERTING TO GPU when available. You MUST use this option with --vectorized (otherwise it is a no-op). You MUST have a suitable version of cupy installed, your cuda operational, etc")
optp.add_option("-o", "--output-file", help="Save result to this file.")
optp.add_option("-O", "--output-format", default='xml', help="[xml|hdf5]")
optp.add_option("-S", "--save-samples", action="store_true", help="Save sample points to output-file. Requires --output-file to be defined.")
optp.add_option("-L", "--save-deltalnL", type=float, default=float("Inf"), help="Threshold on deltalnL for points preserved in output file. Requires --output-file to be defined")
optp.add_option("-P", "--save-P", type=float,default=0, help="Threshold on cumulative probability for points preserved in output file. Requires --output-file to be defined")
optp.add_option("--verbose",action='store_true')
#
# Add the integration options
......@@ -395,8 +412,9 @@ for inst, chan in map(lambda c: c.split("="), opts.channel_name):
print "Reading channel %s from cache %s" % (inst+":"+chan, opts.cache_file)
data_dict[inst] = lalsimutils.frame_data_to_non_herm_hoff(opts.cache_file,
inst+":"+chan, start=start_time, stop=end_time,
window_shape=opts.window_shape)
print "Frequency binning: %f, length %d" % (data_dict[inst].deltaF,
window_shape=opts.window_shape,verbose=opts.verbose)
if opts.verbose:
print "Frequency binning: %f, length %d" % (data_dict[inst].deltaF,
data_dict[inst].data.length)
flow_ifo_dict = {}
......@@ -419,18 +437,21 @@ for inst, psdf in map(lambda c: c.split("="), opts.psd_file):
if isinstance(psd_dict[inst], lal.REAL8FrequencySeries):
psd_fvals = psd_dict[inst].f0 + deltaF*numpy.arange(psd_dict[inst].data.length)
psd_dict[inst].data.data[ psd_fvals < flow_ifo_dict[inst]] = 0 #
elif isinstance(psd_dict[inst], pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries): # for backward compatibility
psd_fvals = psd_dict[inst].f0 + deltaF*numpy.arange(len(psd_dict[inst].data))
psd_dict[inst].data[psd_fvals < ifo_dict[inst]] =0
else:
print 'FAIL'
sys.exit(0)
# elif isinstance(psd_dict[inst], pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries): # for backward compatibility
# psd_fvals = psd_dict[inst].f0 + deltaF*numpy.arange(len(psd_dict[inst].data))
# psd_dict[inst].data[psd_fvals < ifo_dict[inst]] =0
assert psd_dict[inst].deltaF == deltaF
# Highest freq. at which PSD is defined
if isinstance(psd_dict[inst],
pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries):
fmax = psd_dict[inst].f0 + deltaF * (len(psd_dict[inst].data) - 1)
elif isinstance(psd_dict[inst], lal.REAL8FrequencySeries):
#if isinstance(psd_dict[inst],
# pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries):
# fmax = psd_dict[inst].f0 + deltaF * (len(psd_dict[inst].data) - 1)
if isinstance(psd_dict[inst], lal.REAL8FrequencySeries):
fmax = psd_dict[inst].f0 + deltaF * (psd_dict[inst].data.length - 1)
# Assert upper limit of IP integral does not go past where PSD defined
......@@ -881,7 +902,7 @@ if not opts.time_marginalization:
res, var, neff, dict_return = sampler.integrate(likelihood_function, *unpinned_params, **pinned_params)
else: # Sum over time for every point in other extrinsic params
if not opts.rom_integrate_intrinsic:
if not (opts.rom_integrate_intrinsic or opts.vectorized) :
def likelihood_function(right_ascension, declination, phi_orb, inclination,
psi, distance):
dec = numpy.copy(declination).astype(numpy.float64) # get rid of 'object', and allocate space
......@@ -914,6 +935,67 @@ else: # Sum over time for every point in other extrinsic params
return numpy.exp(lnL -manual_avoid_overflow_logarithm)
args = likelihood_function.func_code.co_varnames[:likelihood_function.func_code.co_argcount]
print " --------> Arguments ", args
res, var, neff, dict_return = sampler.integrate(likelihood_function, *unpinned_params, **pinned_params)
elif opts.vectorized: # use array-based multiplications, fewer for loops
lookupNKDict = {}
lookupKNDict={}
lookupKNconjDict={}
ctUArrayDict = {}
ctVArrayDict={}
rholmArrayDict={}
rholms_intpArrayDict={}
epochDict={}
nEvals=0
for det in rholms_intp.keys():
print " Packing ", det
lookupNKDict[det],lookupKNDict[det], lookupKNconjDict[det], ctUArrayDict[det], ctVArrayDict[det], rholmArrayDict[det], rholms_intpArrayDict[det], epochDict[det] = factored_likelihood.PackLikelihoodDataStructuresAsArrays( rholms[det].keys(), rholms_intp[det], rholms[det], cross_terms[det],cross_terms_V[det])
if opts.gpu and xpy_default==cupy:
lookupNKDict[det] = cupy.asarray(lookupNKDict[det])
rholmArrayDict[det] = cupy.asarray(rholmArrayDict[det])
ctUArrayDict[det] = cupy.asarray(ctUArrayDict[det])
ctVArrayDict[det] = cupy.asarray(ctVArrayDict[det])
epochDict[det] = cupy.asarray(epochDict[det])
if (not opts.gpu):
def likelihood_function(right_ascension, declination, phi_orb, inclination,
psi, distance):
global nEvals
tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally
# use EXTREMELY many bits
lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128)
P.phi = right_ascension.astype(float) # cast to float
P.theta = declination.astype(float)
P.tref = float(fiducial_epoch)
P.phiref = phi_orb.astype(float)
P.incl = inclination.astype(float)
P.psi = psi.astype(float)
P.dist = (distance* 1.e6 * lalsimutils.lsu_PC).astype(float) # luminosity distance
lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVector(tvals,
P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max)
nEvals +=len(right_ascension)
return numpy.exp(lnL-manual_avoid_overflow_logarithm)
else:
print " Using CUDA GPU likelihood, if cupy available "
def likelihood_function(right_ascension, declination, phi_orb, inclination,
psi, distance):
global nEvals
tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally
P.phi = xpy_default.asarray(right_ascension, dtype=np.float64) # cast to float
P.theta = xpy_default.asarray(declination.astype,dtype=np.float64)
P.tref = float(fiducial_epoch)
P.phiref = xpy_default.asarray(phi_orb, dtype=np.float64)
P.incl = xpy_default.asarray(inclination, dtype=np.float64)
P.psi = xpy_default.asarray(psi, dtype=np.float64)
P.dist = xpy_default.asarray(distance* 1.e6 * lalsimutils.lsu_PC,dtype=np.float64) # luminosity distance
lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoops(tvals,
P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max)
nEvals +=len(right_ascension)
return identity_convert(xpy_default.exp(lnL-manual_avoid_overflow_logarithm))
args = likelihood_function.func_code.co_varnames[:likelihood_function.func_code.co_argcount]
print " --------> Arguments ", args
res, var, neff, dict_return = sampler.integrate(likelihood_function, *unpinned_params, **pinned_params)
......
import cupy
def tupleset(t, i, value):
l = list(t)
l[i] = value
return tuple(l)
def _basic_simps(y, start, stop, x, dx, axis):
import cupy
nd = len(y.shape)
if start is None:
start = 0
step = 2
slice_all = (slice(None),)*nd
slice0 = tupleset(slice_all, axis, slice(start, stop, step))
slice1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
slice2 = tupleset(slice_all, axis, slice(start+2, stop+2, step))
if x is None: # Even spaced Simpson's rule.
result = cupy.sum(dx/3.0 * (y[slice0]+4*y[slice1]+y[slice2]),
axis=axis)
else:
# Account for possibly different spacings.
# Simpson's rule changes a bit.
h = cupy.diff(x, axis=axis)
sl0 = tupleset(slice_all, axis, slice(start, stop, step))
sl1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
h0 = h[sl0]
h1 = h[sl1]
hsum = h0 + h1
hprod = h0 * h1
h0divh1 = h0 / h1
tmp = hsum/6.0 * (y[slice0]*(2-1.0/h0divh1) +
y[slice1]*hsum*hsum/hprod +
y[slice2]*(2-h0divh1))
result = cupy.sum(tmp, axis=axis)
return result
def simps(y, x=None, dx=1, axis=-1, even='avg'):
"""
Integrate y(x) using samples along the given axis and the composite
Simpson's rule. If x is None, spacing of dx is assumed.
If there are an even number of samples, N, then there are an odd
number of intervals (N-1), but Simpson's rule requires an even number
of intervals. The parameter 'even' controls how this is handled.
Parameters
----------
y : array_like
Array to be integrated.
x : array_like, optional
If given, the points at which `y` is sampled.
dx : int, optional
Spacing of integration points along axis of `y`. Only used when
`x` is None. Default is 1.
axis : int, optional
Axis along which to integrate. Default is the last axis.
even : {'avg', 'first', 'str'}, optional
'avg' : Average two results:1) use the first N-2 intervals with
a trapezoidal rule on the last interval and 2) use the last
N-2 intervals with a trapezoidal rule on the first interval.
'first' : Use Simpson's rule for the first N-2 intervals with
a trapezoidal rule on the last interval.
'last' : Use Simpson's rule for the last N-2 intervals with a
trapezoidal rule on the first interval.
See Also
--------
quad: adaptive quadrature using QUADPACK
romberg: adaptive Romberg quadrature
quadrature: adaptive Gaussian quadrature
fixed_quad: fixed-order Gaussian quadrature
dblquad: double integrals
tplquad: triple integrals
romb: integrators for sampled data
cumtrapz: cumulative integration for sampled data
ode: ODE integrators
odeint: ODE integrators
Notes
-----
For an odd number of samples that are equally spaced the result is
exact if the function is a polynomial of order 3 or less. If
the samples are not equally spaced, then the result is exact only
if the function is a polynomial of order 2 or less.
"""
import cupy
y = cupy.asarray(y)
nd = len(y.shape)
N = y.shape[axis]
last_dx = dx
first_dx = dx
returnshape = 0
if x is not None:
x = cupy.asarray(x)
if len(x.shape) == 1:
shapex = [1] * nd
shapex[axis] = x.shape[0]
saveshape = x.shape
returnshape = 1
x = x.reshape(tuple(shapex))
elif len(x.shape) != len(y.shape):
raise ValueError("If given, shape of x must be 1-d or the "
"same as y.")
if x.shape[axis] != N:
raise ValueError("If given, length of x along axis must be the "
"same as y.")
if N % 2 == 0:
val = 0.0
result = 0.0
slice1 = (slice(None),)*nd
slice2 = (slice(None),)*nd
if even not in ['avg', 'last', 'first']:
raise ValueError("Parameter 'even' must be "
"'avg', 'last', or 'first'.")
# Compute using Simpson's rule on first intervals
if even in ['avg', 'first']:
slice1 = tupleset(slice1, axis, -1)
slice2 = tupleset(slice2, axis, -2)
if x is not None:
last_dx = x[slice1] - x[slice2]
val += 0.5*last_dx*(y[slice1]+y[slice2])
result = _basic_simps(y, 0, N-3, x, dx, axis)
# Compute using Simpson's rule on last set of intervals
if even in ['avg', 'last']:
slice1 = tupleset(slice1, axis, 0)
slice2 = tupleset(slice2, axis, 1)
if x is not None:
first_dx = x[tuple(slice2)] - x[tuple(slice1)]
val += 0.5*first_dx*(y[slice2]+y[slice1])
result += _basic_simps(y, 1, N-2, x, dx, axis)
if even == 'avg':
val /= 2.0
result /= 2.0
result = result + val
else:
result = _basic_simps(y, 0, N-2, x, dx, axis)
if returnshape:
x = x.reshape(saveshape)
return result
......@@ -178,6 +178,7 @@ def ParseStandardArguments():
# Unused: To suck up parameters by letting me do search/replace on text files
parser.add_argument("--unused-argument", default=None, help="Used to sop up arguments I want to replace")
parser.add_argument("--save-no-samples",action='store_true')
args = parser.parse_args()
# Connect debugging options to argument parsing
......
......@@ -72,6 +72,8 @@ P.print_params()
T_est = lalsimutils.estimateWaveformDuration(P)
T_est = P.deltaT*lalsimutils.nextPow2(T_est/P.deltaT)
if T_est < opts.seglen:
T_est =opts.seglen
P.deltaF = 1./T_est
print " Duration ", T_est
if T_est < opts.seglen:
......
import numpy
try:
import cupy
xpy_default=cupy
except:
xpy_default=numpy
def TimeDelayFromEarthCenter(
detector_earthfixed_xyz_metres,
source_right_ascension_radians,
source_declination_radians,
greenwich_mean_sidereal_time,
xpy=xpy_default, dtype=numpy.float64,
):
"""
Parameters
----------
detector_earthfixed_xyz_metres : array_like, shape = det_shape + (3,)
Location of detector(s) relative to Earth's center in meters. May provide
multiple detectors, last axis must be (x,y,z) but other axes can take
whatever form is desired.
source_right_ascension_radians : array_like, shape = sample_shape
Right ascension of source in radians, can be an arbitrary dimensional
array.
source_declination_radians : array_like, shape = sample_shape
Declination of source in radians, can be an arbitrary dimensional array.
greenwich_mean_sidereal_time : float
Should be equivalent to XLALGreenwichMeanSiderealTime(gpstime).
Returns
-------
time_delay_from_earth_center : array_like, shape = det_shape + sample_shape
"""
negative_speed_of_light = xpy.asarray(-299792458.0)
det_shape = detector_earthfixed_xyz_metres.shape[:-1]
sample_shape = source_right_ascension_radians.shape
cos_dec = xpy.cos(source_declination_radians)
greenwich_hour_angle = (
greenwich_mean_sidereal_time - source_right_ascension_radians
)
ehat_src = xpy.empty(sample_shape + (3,), dtype=dtype)
ehat_src[...,0] = cos_dec * xpy.cos(greenwich_hour_angle)
ehat_src[...,1] = -cos_dec * xpy.sin(greenwich_hour_angle)
ehat_src[...,2] = xpy.sin(source_declination_radians)
neg_separation = xpy.inner(detector_earthfixed_xyz_metres, ehat_src)
return xpy.divide(
neg_separation, negative_speed_of_light,
out=neg_separation,
)
def ComputeDetAMResponse(
detector_response_matrix,
source_right_ascension_radians,
source_declination_radians,
source_polarization_radians,
greenwich_mean_sidereal_time,
xpy=xpy_default, dtype_real=numpy.float64, dtype_complex=numpy.complex128,
):
"""
Parameters
----------
detector_response_matrix : array_like, shape = det_shape + (3, 3)
Detector response matrix, or matrices for multiple detectors. Last two
axes must be 3-by-3 response matrix, and may include arbitrary axes before
that for various detectors.
source_right_ascension_radians : array_like, shape = sample_shape
Right ascension of source in radians, can be an arbitrary dimensional
array.
source_declination_radians : array_like, shape = sample_shape
Declination of source in radians, can be an arbitrary dimensional array.
source_polarization_radians : array_like, shape = sample_shape
Polarization angle of source in radians, can be an arbitrary dimensional
array.
greenwich_mean_sidereal_time : float
Should be equivalent to XLALGreenwichMeanSiderealTime(gpstime).
Returns
-------
F : array_like, shape = det_shape + sample_shape
"""
det_shape = detector_response_matrix.shape[:-1]
sample_shape = source_right_ascension_radians.shape
matrix_shape = 3, 3
# Initialize trig matrices.
X = xpy.empty(sample_shape+(3,), dtype=dtype_real)
Y = xpy.empty(sample_shape+(3,), dtype=dtype_real)
# Greenwich hour angle of source in radians.
source_greenwich_radians = (
greenwich_mean_sidereal_time - source_right_ascension_radians
)
# Pre-compute trig functions
cos_gha = xpy.cos(source_greenwich_radians)
sin_gha = xpy.sin(source_greenwich_radians)
cos_dec = xpy.cos(source_declination_radians)
sin_dec = xpy.sin(source_declination_radians)
cos_psi = xpy.cos(source_polarization_radians)
sin_psi = xpy.sin(source_polarization_radians)
# Populate trig matrices.
X[...,0] = -cos_psi*sin_gha - sin_psi*cos_gha*sin_dec
X[...,1] = -cos_psi*cos_gha + sin_psi*sin_gha*sin_dec
X[...,2] = sin_psi*cos_dec
Y[...,0] = sin_psi*sin_gha - cos_psi*cos_gha*sin_dec
Y[...,1] = sin_psi*cos_gha + cos_psi*sin_gha*sin_dec
Y[...,2] = cos_psi*cos_dec
# Compute F for each polarization state.
F_plus = (
X*xpy.inner(X, detector_response_matrix) -
Y*xpy.inner(Y, detector_response_matrix)
).sum(axis=-1)
F_cross = (
X*xpy.inner(Y, detector_response_matrix) +
Y*xpy.inner(X, detector_response_matrix)
).sum(axis=-1)
return F_plus + 1.0j*F_cross
......@@ -4,18 +4,16 @@ import lalsimulation
from glue.ligolw import utils
from pylal.series import make_psd_xmldoc
from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS
from pylal.xlal.datatypes.real8frequencyseries import REAL8FrequencySeries
from pylal.xlal.datatypes.lalunit import LALUnit
from lal.series import make_psd_xmldoc
import lal
from optparse import OptionParser
optp = OptionParser()
optp.add_option("-n", "--nyquist", type=float, default=2048.0, help="Set the nyquist frequency of the PSD. Default is 2048 Hz.")
optp.add_option("-d", "--delta-f", type=float, default=1.0/16, help="Set frequency bin size in Hz. Default is 1/16 Hz.")
optp.add_option("-f", "--low-frequency", type=float, default=0.0, help="Set the low frequency in Hz of the output PSD. Default is zero.")
#optp.add_option( "--ifo", type=str, default="H1", help="Set the low frequency in Hz of the output PSD. Default is zero.")
opts, args = optp.parse_args()
# Generate iLIGO PSD
......@@ -26,10 +24,12 @@ psd = numpy.array(map(lalsimulation.SimNoisePSDaLIGOZeroDetHighPower, f))
psd[0] = 0
# Generate the frequency series
epoch = LIGOTimeGPS(0.0)
epoch = lal.LIGOTimeGPS(0.0)
psddict = {}
d=opts.ifo
for d in ["H1", "L1", "V1"]:
psd_s = REAL8FrequencySeries(name="aLIGO PSD", epoch=epoch, f0=f0, deltaF=df, sampleUnits=LALUnit(""), data=psd)
psd_s = lal.CreateREAL8FrequencySeries(name=d, epoch=epoch, f0=f0, deltaF=df, sampleUnits="s", length=len(psd))
psd_s.data.data = psd
psddict[d] = psd_s
xmldoc = make_psd_xmldoc(psddict)
......
......@@ -4,11 +4,8 @@ import lalsimulation
from glue.ligolw import utils
from pylal.series import make_psd_xmldoc