vectorized orbital separations in inspirals, and prevent resampling spin initial conditions

[mdmould]

Two quick things I noticed while trying to do precession-averaged inspirals from finite orbital separation to infinity:

1. According to the docstring, when doing vectorized inspirals for multiple binaries at once, meaning the source parameters q, chi1, chi2, etc. are arrays with shape (N,), I should be able to specify an output array for r or u with shape (M,N), where M is the number of output points. However, I instead get the error AssertionError: u must be monotonic; minimal working example below.

import numpy as np
import scipy
import precession

np.random.seed(0)

n = 1_000
fGW = 20.0
M_msun = np.random.randn(n) + 30
q = scipy.stats.truncnorm.rvs(
    a = -1 / 0.1, b = 0, loc = 1, scale = 0.1, size = n,
)
chi1 = np.random.rand(n)
chi2 = np.random.rand(n)
theta1 = np.arccos(np.random.rand(n) * 2 - 1)
theta2 = np.arccos(np.random.rand(n) * 2 - 1)
deltaphi = np.random.rand(n) * 2 * np.pi
chieff = (chi1 * np.cos(theta1) + q * chi2 * np.cos(theta2)) / (1 + q)

r = precession.gwfrequency_to_pnseparation(
    theta1, theta2, deltaphi, fGW, q, chi1, chi2, M_msun,
)
rs = np.array([r, np.ones_like(r) * np.inf])
print(rs.shape)

outputs = precession.inspiral_precav(
    theta1 = theta1,
    theta2 = theta2,
    deltaphi = deltaphi,
    r = rs,
    q = q,
    chi1 = chi1,
    chi2 = chi2,
    requested_outputs = ['theta1', 'theta2', 'kappa', 'chieff'],
)

AssertionError: r must be monotonic

This is because of the map on this line https://github.com/dgerosa/precession/blob/master/precession/precession.py#L5012, where the first axis of r or u should have the same size as that of the source parameters. This is a minor fix that just requires r or u have shape (N,M) and only requires changing the docstring. I.e., the following works:

outputs = precession.inspiral_precav(
    theta1 = theta1,
    theta2 = theta2,
    deltaphi = deltaphi,
    r = rs.T,
    q = q,
    chi1 = chi1,
    chi2 = chi2,
    requested_outputs = ['theta1', 'theta2', 'kappa', 'chieff'],
)

2. I would expect that if I give initial conditions in terms of the spin angles theta1, theta2, deltaphi, that the code does not resample those values in the output array. Any resampled values are consistent with the same backwards-in-time evolution, but that said, given those initial conditions the evolution is deterministic and the initial value problem is overconstrained.

outputs['theta1'][:, 0] == theta1

array([False, False, False, False, False, False, False, False, False,
       False])

Would it be preferable to check if the evolution is backwards in time and, if so, not resample the initial condition? Perhaps this isn't an issue in practice as the user will have the initial values in memory anyway and can just ignore the output in the initial condition position.

I'll also note that the docstrings say:

These need to be arrays with lenght >=1, where e.g. r[0] corresponds to the initial condition and r[1:] corresponds to the location where outputs are returned.

This suggests the number of outputs for each binary should be one less than the number of points in the r or u of each binary. But in the above example, where I specify two values of r for each binary (initial condition and infinity), the outputs contain values for both of those r values, contrary to the documentation.

[mdmould]

I also just found the opposite issue happens in precession.deltachioft which states it accepts a time array of shape (N,M) when the source parameters have have shape (N,), but instead it requires (M,N).