A suggestion about the EM part of ICAROGW

I discovered an error while using your program, ICAROGW. I suspect it may be due to a recent update because the previous version did not encounter any issues when independently using the GW170817 gravitational wave event to infer cosmological parameters. However, the latest version appears to have introduced this error. (Specifically, I ran the program using the provided example in the documentation, which aims to constrain the Hubble constant using the GW170817 event.)
After reviewing the code, I believe the error originates from the CBC_vanilla_EM_counterpart class in wrappers.py. This class is designed to support multiple bright sirens but encounters an issue when dealing with a single event. In the case of a single event, the input for the posterior samples from GW170817, kwargs['mass_1'], is a list [m1_1, m1_2, m1_3, ..., m1_N] where len(kwargs['mass_1'].shape) = 1. However, when there are multiple events, kwargs['mass_1'] becomes a list of lists representing individual event samples, resulting in len(kwargs['mass_1'].shape) = 2. The error occurs because the code has a condition (from line 316 to line 317 in wrappers.py file) that assumes len(kwargs['mass_1'].shape) must be 2, which is not met in the single event case.

'''python if len(kwargs['mass_1'].shape) != 2: raise ValueError('The EM counterpart rate wants N_ev x N_samples arrays') ''' This error can be avoided by optimizing the data format, or simplely replace the log_rate_PE function of CBC_vanilla_EM_counterpart class with the following code.

'''python def log_rate_PE(self,prior,**kwargs): xp = get_module_array(prior) sx = get_module_array_scipy(prior)

    ms1, ms2, z = detector2source(kwargs['mass_1'],kwargs['mass_2'],kwargs['luminosity_distance'],self.cw.cosmology) 
    log_dVc_dz=xp.log(self.cw.cosmology.dVc_by_dzdOmega_at_z(z)*4*xp.pi)

    # Sum over posterior samples in Eq. 1.1 on the icarogw2.0 document
    log_weights=self.mw.log_pdf(ms1,ms2)+self.rw.rate.log_evaluate(z)+log_dVc_dz \
    -xp.log(prior)-xp.log(detector2source_jacobian(z,self.cw.cosmology))-xp.log1p(z)

    if self.sw is not None:
        log_weights+=self.sw.log_pdf(**{key:kwargs[key] for key in self.sw.event_parameters})

    if len(kwargs['mass_1'].shape) == 2:
        n_ev = kwargs['mass_1'].shape[0]
        lwtot = xp.empty(kwargs['z_EM'].shape)
        for i in range(n_ev): 
            ww = xp.exp(log_weights[i,:])
            kde_fit = gaussian_kde(z[i,:],weights=ww/ww.sum())   
            lwtot[i,:] = sx.special.logsumexp(log_weights[i,:])-xp.log(kwargs['mass_1'].shape[1])+np2cp(kde_fit.logpdf(cp2np(kwargs['z_EM'][i,:])))
    else:
        lwtot = xp.empty(kwargs['z_EM'].shape)
        ww = xp.exp(log_weights)
        kde_fit = gaussian_kde(z,weights=ww/ww.sum())
        lwtot = sx.special.logsumexp(log_weights)-xp.log(len(kwargs['mass_1']))+np2cp(kde_fit.logpdf(cp2np(kwargs['z_EM'])))

    if not self.scale_free:
        log_out = lwtot + xp.log(self.R0)
    else:
        log_out = lwtot

    return log_out

'''

simone-mastrogiovanni / icarogw

A suggestion about the EM part of ICAROGW #6