BlackHolePerturbationToolkit / FastEMRIWaveforms

Blazingly fast EMRI waveforms
GNU General Public License v3.0
38 stars 26 forks source link

How does FEW calculate the length of generated GW time series? #75

Closed AminBoumerdassi closed 12 months ago

AminBoumerdassi commented 12 months ago

I'm using FEW to train an LSTM auto encoder on generated EMRI waveforms for the purpose of signal detection. This requires me to specify the input dimensionality of my ML model before actually generating EMRIs.

My first thought was to do the standard T/dt to get the length of the strain array, however this calculation doesn't exactly match FEW's array length. I then realised that T is in years, and the solar year is not exactly 365 days so I tried the more accurate 365.242 days, and the day length of 86400.002 seconds. So my new calculation for the length of a EMRI time series with T=1 and dt=10 was:

T=1#years dt=10#seconds array_length=round(86400.002*365.242*T/dt)+1

However, this gives 3155692 time steps when it should be 3155815. How exactly does FEW calculate the length of the strain array?

Thanks

mikekatz04 commented 12 months ago

Hello,

It calculates it like you say but by default it chops off anywhere after the signal plunges. If you add “pad_output=True” to the sum kwargs it will give you a fixed length array no matter what always returning the same length given by int(T/dt). Check the documentation and tutorial to make sure the call is right.

On Tue, Jul 4, 2023 at 8:28 PM AminBoumerdassi @.***> wrote:

I'm using FEW to train an LSTM auto encoder on generated EMRI waveforms for the purpose of signal detection. This requires me to specify the input dimensionality of my ML model before actually generating EMRIs.

My first thought was to do the standard T/dt to get the length of the strain array, however this calculation doesn't exactly match FEW's array length. I then realised that T is in years, and the solar year is not exactly 365 days so I tried the more accurate 365.242 days, and the day length of 86400.002 seconds. So my new calculation for the length of a EMRI time series with T=1 and dt=10 was:

T=1#years dt=10#seconds array_length=round(86400.002365.242T/dt)+1

However, this gives 3155692 time steps when it should be 3155815. How exactly does FEW calculate the length of the strain array?

Thanks

— Reply to this email directly, view it on GitHub https://github.com/BlackHolePerturbationToolkit/FastEMRIWaveforms/issues/75, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFTXYR6B4VU4BYLRBXE7BR3XOS7M7ANCNFSM6AAAAAAZ6JBHFE . You are receiving this because you are subscribed to this thread.Message ID: @.***>

-- Michael Katz Postdoctoral Researcher Astrophysical and Cosmological Relativity Group Max Planck Institute for Gravitational Physics (Albert Einstein Institute) @. @.

AminBoumerdassi commented 12 months ago

Hi Michael,

I've changed the pad_output argument to True however this hasn't changed the length of the output array. Here's what I tried based off of the tutorial:

import cupy as xp
from numpy.random import default_rng

#FEW imports
import sys
import os
from numpy.random import default_rng
from few.trajectory.inspiral import EMRIInspiral
from few.amplitude.romannet import RomanAmplitude
from few.amplitude.interp2dcubicspline import Interp2DAmplitude
from few.waveform import FastSchwarzschildEccentricFlux, GenerateEMRIWaveform

use_gpu = True

# keyword arguments for inspiral generator (RunSchwarzEccFluxInspiral)
inspiral_kwargs={
        "DENSE_STEPPING": 0,  # we want a sparsely sampled trajectory
        "max_init_len": int(1e3),  # all of the trajectories will be well under len = 1000
    }

# keyword arguments for inspiral generator (RomanAmplitude)
amplitude_kwargs = {
    "max_init_len": int(1e3),  # all of the trajectories will be well under len = 1000
    "use_gpu": use_gpu  # GPU is available in this class
}

# keyword arguments for Ylm generator (GetYlms)
Ylm_kwargs = {
    "assume_positive_m": False  # if we assume positive m, it will generate negative m for all m>0
}

# keyword arguments for summation generator (InterpolatedModeSum)
sum_kwargs = {
    "use_gpu": use_gpu,  # GPU is availabel for this type of summation
    "pad_output": True,
}

#Initialise the fast Scharwzchild waveform object
few = FastSchwarzschildEccentricFlux(
    inspiral_kwargs=inspiral_kwargs,
    amplitude_kwargs=amplitude_kwargs,
    Ylm_kwargs=Ylm_kwargs,
    sum_kwargs=sum_kwargs,
    use_gpu=use_gpu,
    # num_threads=num_threads,  # 2nd way for specific classes
)

##generating EMRIs with randomised parameters
no_of_EMRIs= 1
dt = 10.0 # time step (s)
T = 1.0#/365 # obs. time (years)
dim= round(86400*365.242*T/dt)#should be the output length of array

EMRIs= xp.empty((no_of_EMRIs, dim, 2))

rng = default_rng(seed=2023)

set_of_eta= rng.uniform(1e-6,1e-4, size= no_of_EMRIs)#Mass ratio for calculating mu
set_of_M= rng.uniform(1e4, 1e7, size = no_of_EMRIs)
set_of_mu= set_of_eta*set_of_M
set_of_e0= rng.uniform(0, .7, size = no_of_EMRIs)#Initial eccentricity
set_of_p0= rng.uniform(10, 16+set_of_e0)#Needs to be based on set_of_e0
set_of_theta= rng.uniform(-np.pi/2, +np.pi/2, size = no_of_EMRIs)#polar viewing angle
set_of_phi= rng.uniform(0, 2*np.pi, size = no_of_EMRIs)
set_of_phi_phi0= rng.uniform(0, 2*np.pi, size = no_of_EMRIs)
set_of_phi_r0= rng.uniform(0, 2*np.pi, size = no_of_EMRIs)
dist = None #effectively in source frame

#Iterating through parameter sets
for i, M, mu, p0, e0, theta, phi, phi_phi0, phi_r0  in zip(range(no_of_EMRIs), set_of_M, set_of_mu, set_of_p0,
                                                                set_of_e0, set_of_theta, set_of_phi, set_of_phi_phi0, set_of_phi_r0):
    wave= few(M, mu, p0, e0, theta, phi, dist=dist, Phi_phi0 = phi_phi0, Phi_r0=phi_r0, dt=dt, T=T)
    EMRIs[i,:,0]= wave.real
    EMRIs[i,:,1]= wave.imag

And I get the following error:

ValueError                                Traceback (most recent call last)
/dev/shm/jobs/37546032/ipykernel_161678/96694761.py in <module>
     25                                                                 set_of_e0, set_of_theta, set_of_phi, set_of_phi_phi0, set_of_phi_r0):
     26     wave= few(M, mu, p0, e0, theta, phi, dist=dist, Phi_phi0 = phi_phi0, Phi_r0=phi_r0, dt=dt, T=T)
---> 27     EMRIs[i,:,0]= wave.real
     28     EMRIs[i,:,1]= wave.imag

cupy/_core/core.pyx in cupy._core.core._ndarray_base.__setitem__()

cupy/_core/_routines_indexing.pyx in cupy._core._routines_indexing._ndarray_setitem()

cupy/_core/_routines_indexing.pyx in cupy._core._routines_indexing._scatter_op()

cupy/_core/_kernel.pyx in cupy._core._kernel.ufunc.__call__()

cupy/_core/internal.pyx in cupy._core.internal._broadcast_core()

ValueError: operands could not be broadcast together with shapes (3155815,) (3155691,)

Checking the array length of the generated EMRI with wave.shape returns (3155815,). Checking the expected array length with int(86400*365.242*T/dt) returns 3155690

CChapmanbird commented 12 months ago

FEW uses an internal constant that defines the length of a year in seconds, YRSID_SI found in few.utils.constants - you should use this for all datastream length-related calculations.

AminBoumerdassi commented 12 months ago

Thanks Christian, that did the trick!