Commit 1a91e3ba authored by Edward Fauchon-Jones's avatar Edward Fauchon-Jones
Browse files

Add function to find waveform peak

parent eea7de2a
......@@ -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
......@@ -233,13 +235,76 @@ def ismonotonic(x):
"""
dx = np.diff(x)
return np.all(dx <= 0) or np.all(dx >= 0)
def is_peak_near_zero():
"""
A function that checks that the waveform data (amp and phase)
have a peak that is near zero.
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.
"""
raise NotImplementedError
# 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'][:], 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'][:], 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
dpAll = pAll.derivative()
extremumPoints = dpAll.roots()
extremum = pAll(extremumPoints)
idx = np.argmax(extremum)
maximumPoint = extremumPoints[idx]
maximum = extremum[idx]
return (maximumPoint, maximum)
# General Fields
class Type(Spec):
......
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