Commit 0ed8a149 authored by Sebastian Khan's avatar Sebastian Khan
Browse files

Merge branch 'add-peak-check' into 'master'

Add peak check

See merge request !10
parents 118f6a28 c255ceee
Pipeline #25134 passed with stage
in 2 minutes and 6 seconds
......@@ -61,7 +61,8 @@ format1 = OrderedDict((
("mean_anomaly", MeanAnomaly))))))
format1Interfield = OrderedDict((
("mass-ordering", MassOrdering),))
("mass-ordering", MassOrdering),
("peak-near-zero", PeakNearZero)))
format2 = OrderedDict((
('mass1-vs-time', Mass1VsTime),
......
......@@ -17,6 +17,7 @@ import h5py as h5
from . import errors as err
import numpy as np
import re
from scipy.interpolate import insert, splrep, PPoly
class Spec(object):
......@@ -215,6 +216,7 @@ def isclose(a, b, rel_tol=1e-09, abs_tol=0.0):
"""
return abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
def ismonotonic(x):
"""A Function for testing is a list/array is monotonic
......@@ -234,6 +236,79 @@ def ismonotonic(x):
dx = np.diff(x)
return np.all(dx <= 0) or np.all(dx >= 0)
def get_waveform_peak(sim):
"""Determine the peak of an NR waveform
Determine the maximum point and maximum value of the specified waveform as
defined by (4) in https://arxiv.org/abs/1703.01076.
The peak is determined analytically as in general the maximum point will
not be an element of the `X` dataset of any of the amplitude groups.
To analytically determine the waveform peak, the amplitude splines are
constructed from their respective groups. Knots are inserted such that
all splines share common knots. The splines are then expressed as
peicewise polynomials which can then be summed to form a single piecewise
polynomial.
Parameters
----------
sim: h5.File
An open HDF5 file like interface that contains LVCNR waveform data.
Returns
-------
tPeak: float
Maximum point of the specified waveform.
hPeak:
Maximum value of the specified waveform.
"""
# Determine internal knots for waveform piecewise polynomial
ampKeys = [k for k in sim.keys() if k.startswith('amp_')]
interiorKnots = np.array([])
for k in ampKeys:
interiorKnots = np.union1d(interiorKnots, sim[k]['X'][3:-3])
# Prepare container for piecewise polynomial
amp22 = sim['amp_l2_m2']
tck22 = splrep(amp22['X'][:], amp22['Y'][:]**2, s=0, k=5)
interiorKnots22 = np.setdiff1d(
interiorKnots, amp22['X'][3:-3])
tck22 = reduce(
lambda tck, knot: insert(knot, tck),
interiorKnots22,
tck22)
pAll = PPoly.from_spline(tck22)
# Add remaining multipole moments to piecewise polynomial
ampKeys = [k for k in ampKeys if k != 'amp_l2_m2']
for k in ampKeys:
amplm = sim[k]
tcklm = splrep(amplm['X'][:], amplm['Y'][:]**2, s=0, k=5)
interiorKnotslm = np.setdiff1d(
interiorKnots, amplm['X'][3:-3])
tcklm = reduce(
lambda tck, knot: insert(knot, tck),
interiorKnotslm,
tcklm)
plm = PPoly.from_spline(tcklm)
pAll.c += plm.c
# Find the global maximum point of the piecewise polynomial in domain
dpAll = pAll.derivative()
extremumPoints = dpAll.roots()
t0, t1 = sim['amp_l2_m2/X'][:][np.r_[0, -1]]
mask = (extremumPoints >= t0) * (extremumPoints <= t1)
extremumPoints = extremumPoints[mask]
extremum = pAll(extremumPoints)
idx = np.argmax(extremum)
maximumPoint = extremumPoints[idx]
maximum = extremum[idx]
return (maximumPoint, maximum)
# General Fields
class Type(Spec):
"""Specification for the `type` field"""
......@@ -800,3 +875,35 @@ class Amplm(ROMSplineSpec):
lm = [Amplm(name=k) for k in keys]
return lm
class PeakNearZero(InterfieldSpec):
"""Specification for the location of the waveform peak"""
name = "peak-near-zero"
validMsg = "waveform peak is near zero"
invalidMsg = "waveform peak is not near zero"
def valid(self, sim, tol=10.0):
# Validate multipole moment amplitudes
ampModes = Amplm.getlm(sim)
for ampMode in ampModes:
ampModeValid = ampMode.valid(sim)
if not isinstance(ampModeValid, err.Valid):
self.invalidMsg = "amp_l{0:d}_m{1:d} is invalid".format(
ampMode.l, ampMode.m)
return err.InvalidInterfields(self)
maximumPoint, maximum = get_waveform_peak(sim)
if np.abs(maximumPoint) < tol:
self.validMsg = (
"waveform peak is at {0:.2f}M"
" which is less than {1:.2f}M from zero").format(
maximumPoint, tol)
return err.ValidInterfield(self)
else:
self.invalidMsg = (
"waveform peak is at {0:.2f}M"
" which is greater than {1:.2f}M from zero").format(
maximumPoint, tol)
return err.InvalidInterfield(self)
......@@ -70,6 +70,7 @@ templateOutput = """# Format 1
# Format 1 (Interfield)
- [=] mass-ordering (mass1 >= mass2)
- [=] peak-near-zero (waveform peak is at 0.05M which is less than 10.00M from zero)
# Format 2
......@@ -1824,6 +1825,10 @@ class TestAmplm(TestROMSpline):
("- [WRONG TYPE] {0:s} "
"(<class \'h5py._hl.dataset.Dataset\'>) "
"(Type must be <class \'h5py._hl.group.Group\'>)").format(name))
self.setNamedOutput(
'peak-near-zero',
('- [INVALID FIELDS] peak-near-zero '
'(Field dependencies are invalid)'))
(output, returncode) = helper.lvcnrcheck(['-f', '3', self.f.name], returncode=True)
assert output.strip() == self.output
......@@ -1836,6 +1841,30 @@ class TestAmplm(TestROMSpline):
"(<class \'h5py._hl.group.Group\'>) "
"(Field has subfields [invalid-1, invalid-2] "
"but should have [X, Y, deg, errors, tol])").format(name))
self.setNamedOutput(
'peak-near-zero',
('- [INVALID FIELDS] peak-near-zero '
'(Field dependencies are invalid)'))
(output, returncode) = helper.lvcnrcheck(['-f', '3', self.f.name], returncode=True)
assert output.strip() == self.output
assert returncode == 1
class TestPeakNearZero(TestInterfield):
name = 'peak-near-zero'
def test_invalid_peak_location(self):
self.setOutput(
('- [INVALID INTERFIELD] peak-near-zero (waveform peak is at '
'-99.95M which is greater than 10.00M from zero)'))
nr = h5.File(self.f.name)
keys = [k for k in nr.keys() if k.startswith('amp_l')]
for k in keys:
nr[k]['X'][:] -= 100.0
nr.close()
(output, returncode) = helper.lvcnrcheck(['-f', '3', self.f.name], returncode=True)
assert output.strip() == self.output
......
......@@ -84,6 +84,14 @@ def createValidSim():
Phaselm(l, m).createField(sim)
Amplm(l, m).createField(sim)
# Populate simulation with real fixture data
fixture = h5.File('{0:s}/fixture/SXS_BBH_0002_Res5.min.h5'.format(
os.path.dirname(os.path.abspath(__file__))), 'r')
for k in fixture.keys():
del(sim[k])
fixture.copy(fixture[k], sim)
fixture.close()
sim.close()
return f
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment