Calculating derivatives using numdifftools

[harshalc03]

Hi, I'm trying to use waveforms present in LAL but facing an issue when the Fisher Matrix is calculated(to be specific, when using finite differentiation method through numdifftools). The code is as follows:

dets1 = {det:alldets[det] for det in ['H1', 'L1']}

dets1['L1']['asd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', 'LIGO_L_ASD_GW150914.txt')
dets1['H1']['asd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', 'LIGO_H_ASD_GW150914.txt')

signals_L1H1_XPHM = {}

for det in dets1.keys():
    signals_L1H1_XPHM[det] = GWSignal(wave.LAL_WF('IMRPhenomXPHM', 
                                                  is_HigherModes=True,
                                                  is_Precessing=True),
                             psd_path= dets1[det]['asd_path'],
                             detector_shape= dets1[det]['shape'],
                             det_lat= dets1[det]['lat'],
                             det_long= dets1[det]['long'],
                             det_xax= dets1[det]['xax'],
                             verbose=False,
                             useEarthMotion=True,
                             fmin=20.,
                             fmax=1000.,
                             IntTablePath=None)

detnet_L1H1_XPHM = DetNet(signals_L1H1_XPHM)

fish_mat_L1H1_XPHM = detnet_L1H1_XPHM.FisherMatr(GW150914)

The error I get is as follows:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[23], line 1
----> 1 fish_mat_L1H1_XPHM = detnet_L1H1_XPHM.FisherMatr(GW150914)

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/gwfast/network.py:103, in DetNet.FisherMatr(self, evParams, return_all, **kwargs)
    101 if self.verbose:
    102     print('Computing Fisher for %s...' %d)
--> 103 F_ = self.signals[d].FisherMatr(evParams, return_all=return_all, **kwargs) 
    104 #totF +=  self.signals[d].FisherMatr(evParams, **kwargs) 
    105 if self.signals[d].detector_shape=='T' and return_all:

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/gwfast/signal.py:909, in GWSignal.FisherMatr(self, evParams, res, df, spacing, use_m1m2, use_chi1chi2, use_prec_ang, computeDerivFinDiff, computeAnalyticalDeriv, return_all, **kwargs)
    905 allFishers=[]
    907 if self.detector_shape=='L': 
    908     # Compute derivatives
--> 909     FisherDerivs = self._SignalDerivatives_use(fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=0., use_m1m2=use_m1m2, use_chi1chi2=use_chi1chi2, use_prec_ang=use_prec_ang, computeAnalyticalDeriv=computeAnalyticalDeriv, computeDerivFinDiff=computeDerivFinDiff, **kwargs)
    910     # Change the units of the tcoal derivative from days to seconds (this improves conditioning)
    911     FisherDerivs = onp.array(FisherDerivs)

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/gwfast/signal.py:1209, in GWSignal._SignalDerivatives(self, fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot, use_m1m2, use_chi1chi2, use_prec_ang, computeDerivFinDiff, computeAnalyticalDeriv, stepNDT, methodNDT, **kwargs)
   1206                     evpars = [Mc, eta, iota, chiS, chiA, ecc]
   1208 dh = ndt.Jacobian(GWstrainUse, step=stepNDT, method=methodNDT, order=2, n=1)
-> 1209 FisherDerivs = np.asarray(dh(evpars))
   1210 if len(FisherDerivs.shape) == 2: #len(Mc) == 1:
   1211     FisherDerivs = FisherDerivs[:,:,np.newaxis]

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/numdifftools/core.py:431, in Jacobian.__call__(self, x, *args, **kwds)
    430 def __call__(self, x, *args, **kwds):
--> 431     return super(Jacobian, self).__call__(np.atleast_1d(x), *args, **kwds)

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/numdifftools/core.py:289, in Derivative.__call__(self, x, *args, **kwds)
    287 with np.errstate(divide='ignore', invalid='ignore'):
    288     results, f_xi = self._derivative(x_i, args, kwds)
--> 289     derivative, info = self._extrapolate(*results)
    290 if self.full_output:
    291     return derivative, self.info(f_xi, *info)

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/numdifftools/limits.py:206, in _Limit._extrapolate(self, results, steps, shape)
    204 if len(der1) > 2:
    205     der1, errors1, steps = self._wynn_extrapolate(der1, steps)
--> 206 der, info = self._get_best_estimate(der1, errors1, steps, shape)
    207 return der, info

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/numdifftools/limits.py:186, in _Limit._get_best_estimate(der, errors, steps, shape)
    184 @staticmethod
    185 def _get_best_estimate(der, errors, steps, shape):
--> 186     errors += _Limit._add_error_to_outliers(der)
    187     idx = _Limit._get_arg_min(errors)
    188     final_step = steps.flat[idx].reshape(shape)

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/numdifftools/limits.py:170, in _Limit._add_error_to_outliers(der, trim_fact)
    168         p25, median, p75 = np.nanpercentile(der, [25,50, 75], axis=0) 
    169     else:
--> 170         p25, median, p75 = np.percentile(der, [25,50, 75], axis=0)
    172     iqr = np.abs(p75 - p25)
    173 except ValueError as msg:

File ~/miniconda3/envs/ProjSB/lib/python3.11/site-packages/numpy/lib/function_base.py:4277, in percentile(a, q, axis, out, overwrite_input, method, keepdims, interpolation)
   4275 a = np.asanyarray(a)
   4276 if a.dtype.kind == "c":
-> 4277     raise TypeError("a must be an array of real numbers")
   4279 q = np.true_divide(q, 100)
   4280 q = asanyarray(q)  # undo any decay that the ufunc performed (see gh-13105)

TypeError: a must be an array of real numbers

I get the same error when I initialise a signal with python IMRPhenomHM or IMRPhenomD and set computeDerivFinDiff=True in the FisherMatr function. So it is an issue with numdifftools.

I'm not able to understand how to resolve the issue. Any help will be appreciated.

[harshalc03]

Ok, so I was able to resolve the issue by just commenting out these 2 lines in the source code:

if a.dtype.kind == "c":
     raise TypeError("a must be an array of real numbers")

But please note that I'm not aware about the other consequences that I'll face after doing this so this might not be the best solution.