diff --git a/pyfstat/mcmc_based_searches.py b/pyfstat/mcmc_based_searches.py index 9bb1c223428ad7e382a93ab90029d31d00ec3325..e34dd836873a0035b295d2a44815e67e14640eb4 100644 --- a/pyfstat/mcmc_based_searches.py +++ b/pyfstat/mcmc_based_searches.py @@ -37,7 +37,7 @@ class MCMCSearch(core.BaseSearchClass): def __init__(self, label, outdir, theta_prior, tref, minStartTime, maxStartTime, sftfilepath=None, nsteps=[100, 100], nwalkers=100, ntemps=1, log10temperature_min=-5, - theta_initial=None, scatter_val=1e-10, + theta_initial=None, scatter_val=1e-10, rhohatmax=1000, binary=False, BSGL=False, minCoverFreq=None, maxCoverFreq=None, detectors=None, earth_ephem=None, sun_ephem=None, injectSources=None, assumeSqrtSX=None): @@ -70,6 +70,10 @@ class MCMCSearch(core.BaseSearchClass): log10temperature_min float < 0 The log_10(tmin) value, the set of betas passed to PTSampler are generated from np.logspace(0, log10temperature_min, ntemps). + rhohatmax: float + Upper bound for the SNR scale parameter (required to normalise the + Bayes factor) - this needs to be carefully set when using the + evidence. binary: Bool If true, search over binary parameters detectors: str @@ -107,6 +111,8 @@ class MCMCSearch(core.BaseSearchClass): if args.clean and os.path.isfile(self.pickle_path): os.rename(self.pickle_path, self.pickle_path+".old") + self.lnlikelihoodcoef = np.log(70./self.rhohatmax**4) + self._log_input() def _log_input(self): @@ -139,7 +145,7 @@ class MCMCSearch(core.BaseSearchClass): self.fixed_theta[theta_i] = theta[j] FS = search.compute_fullycoherent_det_stat_single_point( *self.fixed_theta) - return FS + return FS + self.lnlikelihoodcoef def _unpack_input_theta(self): full_theta_keys = ['F0', 'F1', 'F2', 'Alpha', 'Delta'] @@ -1161,7 +1167,7 @@ class MCMCSearch(core.BaseSearchClass): maxtwoF = self.logl(p, self.search) self.search.BSGL = self.BSGL else: - maxtwoF = maxlogl + maxtwoF = maxlogl - self.lnlikelihoodcoef repeats = [] for i, k in enumerate(self.theta_keys): @@ -1431,9 +1437,10 @@ class MCMCGlitchSearch(MCMCSearch): def __init__(self, label, outdir, sftfilepath, theta_prior, tref, minStartTime, maxStartTime, nglitch=1, nsteps=[100, 100], nwalkers=100, ntemps=1, log10temperature_min=-5, - theta_initial=None, scatter_val=1e-10, dtglitchmin=1*86400, - theta0_idx=0, detectors=None, BSGL=False, minCoverFreq=None, - maxCoverFreq=None, earth_ephem=None, sun_ephem=None): + theta_initial=None, scatter_val=1e-10, rhohatmax=1000, + dtglitchmin=1*86400, theta0_idx=0, detectors=None, + BSGL=False, minCoverFreq=None, maxCoverFreq=None, + earth_ephem=None, sun_ephem=None): """ Parameters ---------- @@ -1466,6 +1473,10 @@ class MCMCGlitchSearch(MCMCSearch): dtglitchmin: int The minimum duration (in seconds) of a segment between two glitches or a glitch and the start/end of the data + rhohatmax: float + Upper bound for the SNR scale parameter (required to normalise the + Bayes factor) - this needs to be carefully set when using the + evidence. nwalkers, ntemps: int, The number of walkers and temperates to use in the parallel tempered PTSampler. @@ -1513,6 +1524,8 @@ class MCMCGlitchSearch(MCMCSearch): self.old_data_is_okay_to_use = self._check_old_data_is_okay_to_use() self._log_input() + self.lnlikelihoodcoef = (self.nglitch+1)*np.log(70./self.rhohatmax**4) + def _initiate_search_object(self): logging.info('Setting up search object') self.search = core.SemiCoherentGlitchSearch( @@ -1546,7 +1559,7 @@ class MCMCGlitchSearch(MCMCSearch): for j, theta_i in enumerate(self.theta_idxs): self.fixed_theta[theta_i] = theta[j] FS = search.compute_nglitch_fstat(*self.fixed_theta) - return FS + return FS + self.lnlikelihoodcoef def _unpack_input_theta(self): glitch_keys = ['delta_F0', 'delta_F1', 'tglitch'] @@ -1682,10 +1695,11 @@ class MCMCSemiCoherentSearch(MCMCSearch): def __init__(self, label, outdir, theta_prior, tref, sftfilepath=None, nsegs=None, nsteps=[100, 100, 100], nwalkers=100, binary=False, ntemps=1, log10temperature_min=-5, - theta_initial=None, scatter_val=1e-10, detectors=None, - BSGL=False, minStartTime=None, maxStartTime=None, - minCoverFreq=None, maxCoverFreq=None, earth_ephem=None, - sun_ephem=None, injectSources=None, assumeSqrtSX=None): + theta_initial=None, scatter_val=1e-10, rhohatmax=1000, + detectors=None, BSGL=False, minStartTime=None, + maxStartTime=None, minCoverFreq=None, maxCoverFreq=None, + earth_ephem=None, sun_ephem=None, injectSources=None, + assumeSqrtSX=None): """ """ @@ -1713,6 +1727,8 @@ class MCMCSemiCoherentSearch(MCMCSearch): self._log_input() + self.lnlikelihoodcoef = self.nsegs * np.log(70./self.rhohatmax**4) + def _get_data_dictionary_to_save(self): d = dict(nsteps=self.nsteps, nwalkers=self.nwalkers, ntemps=self.ntemps, theta_keys=self.theta_keys, @@ -1742,7 +1758,7 @@ class MCMCSemiCoherentSearch(MCMCSearch): self.fixed_theta[theta_i] = theta[j] FS = search.run_semi_coherent_computefstatistic_single_point( *self.fixed_theta) - return FS + return FS + self.lnlikelihoodcoef class MCMCFollowUpSearch(MCMCSemiCoherentSearch): @@ -2097,7 +2113,7 @@ class MCMCTransientSearch(MCMCSearch): if in_theta[1] > self.maxStartTime: return -np.inf FS = search.run_computefstatistic_single_point(*in_theta) - return FS + return FS + self.lnlikelihoodcoef def _unpack_input_theta(self): full_theta_keys = ['transient_tstart',