From 1241c64674aa45a43c25b043f1bd16778bade737 Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Mon, 4 May 2026 20:24:00 -0300 Subject: [PATCH 01/10] bump requirements Namely, SSPtools and limepy have had important updates that we are now going to explicitly require. --- pyproject.toml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 0231d81..1adc7e8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "GCfit" -version = "3.1.2" +version = "3.2.0" description = "Fitting of multimass equilibrium models of globular clusters" authors = [ {name = "Nolan Dickson", email = "nolan.dickson@smu.ca"}, @@ -8,15 +8,15 @@ authors = [ {name = "Vincent Henault-Brunet", email= "Vincent.Henault@smu.ca"} ] readme = "README.md" -requires-python = ">=3.9" +requires-python = ">=3.10" keywords = ["GCfit", "Globular Cluster", "Star Cluster"] license = {file = "LICENSE"} dependencies = [ "numpy", "scipy", "astropy", - "astro-limepy", - "astro-ssptools", + "astro-limepy>=1.3.0", + "astro-ssptools>=3.0.0", "emcee", "dynesty>=3.0.0", "h5py", From 4b902b32b36ac521556facf8955a0d1171d247cb Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Mon, 4 May 2026 20:26:27 -0300 Subject: [PATCH 02/10] remove all references to `BHret` for `BH_ret_dyn` With ssptools>3.0.0, the `BH_ret_int` parameter is removed entirely and the `BH_ret_dyn` parameter now means something different, as it does not include the natal kicks. To make this breaking change explicit, I'm removing `BHret` entirely in favour of the actual `BH_ret_dyn`. Nothing functionally changes about how it is used, but its meaning has changed. --- gcfit/analysis/runs.py | 2 ++ gcfit/core/data.py | 31 ++++++++++++++-------------- gcfit/probabilities/priors.py | 10 ++++++--- gcfit/probabilities/probabilities.py | 2 +- gcfit/scripts/GCfitter.py | 2 +- gcfit/util/probabilities.py | 2 ++ 6 files changed, 28 insertions(+), 21 deletions(-) diff --git a/gcfit/analysis/runs.py b/gcfit/analysis/runs.py index 8656e95..f5db82d 100644 --- a/gcfit/analysis/runs.py +++ b/gcfit/analysis/runs.py @@ -41,6 +41,7 @@ 'a2': r'\alpha_2', 'a3': r'\alpha_3', 'BHret': r'\mathrm{BH}_{ret}', + 'BH_ret_dyn': r'\mathrm{BH}_{\mathrm{ret},\mathrm{dyn}}', 'd': r'd', # Evolved Model Parameters 'M0': r'M_{0}', @@ -94,6 +95,7 @@ 'rhoh0': r'M_\odot\ \mathrm{pc^{-3}}', 's2': r'\mathrm{arcmin^{-4}}', 'BHret': r'\%', + 'BH_ret_dyn': r'\%', 'd': r'\mathrm{kpc}', 'Ndot': r'\dot{N}', 'RA': r'\deg', diff --git a/gcfit/core/data.py b/gcfit/core/data.py index 63e716d..1a0335e 100644 --- a/gcfit/core/data.py +++ b/gcfit/core/data.py @@ -21,7 +21,7 @@ DEFAULT_FREE_PARAMS = ( 'W0', 'M', 'rh', 'ra', 'g', 'delta', - 's2', 'F', 'a1', 'a2', 'a3', 'BHret', 'd', + 's2', 'F', 'a1', 'a2', 'a3', 'BH_ret_dyn', 'd', ) DEFAULT_FREE_EV_PARAMS = ( @@ -800,10 +800,10 @@ class Model(lp.limepy): The high-mass IMF exponent (representing masses between `m_breaks[2:4]`). Defaults to 2.3, matching Kroupa (2001). - BHret : float, optional - The black hole retention fraction, representing the percentage (between - 0 and 100) of black holes retained after dynamical ejections and natal - kicks. + BH_ret_dyn : float, optional + The dynamical black hole retention fraction, representing the percentage + (between 0 and 100) of black holes retained after dynamical ejections. + Note this does *not* include natal kicks. See `ssptools` for details. d : float or astropy.Quantity, optional Distance to the cluster, from Earth, in kiloparsecs. Mainly used for any @@ -992,7 +992,7 @@ def __str__(self): return "Model" def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, - NS_ret, BH_ret_int, BHret, natal_kicks, vesc, + NS_ret, BH_ret_dyn, natal_kicks, vesc, kick_method, f_kick, SNe_method, kick_vdisp, kick_slope, kick_scale, **kwargs): '''Compute an evolved mass function using `ssptools.EvolvedMF`''' @@ -1010,8 +1010,7 @@ def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, esc_rate=esc_rate, tcc=tcc, NS_ret=NS_ret, - BH_ret_int=BH_ret_int, - BH_ret_dyn=BHret / 100., + BH_ret_dyn=BH_ret_dyn / 100., natal_kicks=natal_kicks, vesc=vesc.value, kick_method=kick_method, @@ -1102,7 +1101,7 @@ def _extract_indiv_attrs(self, mask): rhoj=rhoj, Sigmaj=Sigmaj, f=f, rh=rh) def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, - a1=1.3, a2=2.3, a3=2.3, BHret=5.0, d=5, + a1=1.3, a2=2.3, a3=2.3, BH_ret_dyn=5.0, d=5, s2=0., F=1., J=1., *, observations=None, age=None, FeH=None, m_breaks=[0.1, 0.5, 1.0, 150], nbins=[5, 5, 20], tracer_masses=None, tcc=0.0, NS_ret=0.1, BH_ret_int=1.0, @@ -1131,7 +1130,7 @@ def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, self.theta = dict(W0=W0.value, M=M.to_value('1e6 Msun'), rh=rh.value, ra=np.log10(ra.value), g=g, delta=delta, - a1=a1, a2=a2, a3=a3, BHret=BHret, + a1=a1, a2=a2, a3=a3, BH_ret_dyn=BH_ret_dyn, s2=s2, F=F, J=J, d=d.value) self.d = d @@ -1181,7 +1180,7 @@ def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, self._mf = self._evolve_mf(m_breaks, a1, a2, a3, nbins, self.FeH, self.age, esc_rate, tcc, - NS_ret, BH_ret_int, BHret, + NS_ret, BH_ret_dyn, natal_kicks, self.vesc0, kick_method, f_kick, SNe_method, kick_vdisp, kick_slope, kick_scale, **MF_kwargs) @@ -1265,6 +1264,7 @@ def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, mssg = (f"Model extent is not finite (rt>{self.rt:.2f}). " "Model parameters must be adjusted") + raise ValueError(mssg) from err elif "ode not successful" in cause: @@ -1756,7 +1756,7 @@ class EvolvedModel(Model): ''' def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, - NS_ret, BH_ret_int, BHret, natal_kicks, vesc, + NS_ret, BH_ret_dyn, natal_kicks, vesc, kick_method, f_kick, SNe_method, kick_vdisp, kick_slope, kick_scale, **kwargs): '''Alternative MF init using prior-computed IMF and clusterBH outputs''' @@ -1772,7 +1772,6 @@ def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, N0=self._clusterbh.N, # N is N0 tcc=tcc, NS_ret=NS_ret, - BH_ret_int=BH_ret_int, natal_kicks=natal_kicks, vesc=vesc.value, esc_norm='M', @@ -1976,11 +1975,11 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, mssg = f'Too few clusterBH timesteps created: t={self._clusterbh.t}' raise ValueError(mssg) - BHret = -1 # Spoof unneeded BH retention fraction for `Model` + BH_ret_dyn = -1 # Spoof unneeded BH retention fraction for `Model` # Explicitly specify everything so we can get the correct Signature super().__init__(W0, M, rh, g=g, delta=delta, ra=ra, - a1=a1, a2=a2, a3=a3, BHret=BHret, d=d, + a1=a1, a2=a2, a3=a3, BH_ret_dyn=BH_ret_dyn, d=d, meq=meq, eta=eta, zeta=zeta, s2=s2, F=F, J=J, observations=observations, age=age, FeH=FeH, m_breaks=m_breaks, nbins=nbins, @@ -1998,7 +1997,7 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, # reset theta to use initial values self.theta = dict(W0=W0, M0=M0.to_value('1e6 Msun'), rh0=rh0.value, ra=np.log10(ra), g=g, delta=delta, - a1=a1, a2=a2, a3=a3, BHret=BHret, + a1=a1, a2=a2, a3=a3, BH_ret_dyn=BH_ret_dyn, s2=s2, F=F, J=J, d=d.value) def get_visualizer(self): diff --git a/gcfit/probabilities/priors.py b/gcfit/probabilities/priors.py index c9d3780..9fe0cca 100644 --- a/gcfit/probabilities/priors.py +++ b/gcfit/probabilities/priors.py @@ -653,6 +653,10 @@ class FunctionalUniformPrior(UniformPrior): *not* be applied. This may change how the bounding functions must be written, to account for this. + For example, a prior on `rh0` which sets bounds on the initial + density of 1e2 Date: Mon, 4 May 2026 20:34:26 -0300 Subject: [PATCH 03/10] switch to the new relative stopping conditions Introduced in limepy 1.3.0, these new conditions will mean that we no longer have to guess at a much lower diffcrit than the default just to avoid having the BH bins way off of the inputs. The default 0.1% in each bin should be enough now. --- gcfit/core/data.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/gcfit/core/data.py b/gcfit/core/data.py index 1a0335e..36de951 100644 --- a/gcfit/core/data.py +++ b/gcfit/core/data.py @@ -1110,7 +1110,8 @@ def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, f_kick=None, SNe_method='rapid', vesc=90, kick_vdisp=265., kick_slope=1, kick_scale=20, MF_kwargs=None, meanmassdef='global', ode_maxstep=1e10, ode_rtol=1e-7, - diffcrit=1e-8, max_mf_iter=100, mf_iter_index=0.5): + diffcrit=1e-3, max_mf_iter=100, mf_iter_index=0.5, + diffdef='rel'): # ------------------------------------------------------------------ # Add/convert units of some quantities. Supports quantities as inputs @@ -1252,7 +1253,8 @@ def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, ode_rtol=ode_rtol, diffcrit=diffcrit, max_mf_iter=max_mf_iter, - mf_iter_index=mf_iter_index + mf_iter_index=mf_iter_index, + diffdef=diffdef ) try: @@ -1797,8 +1799,8 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, f_kick=None, SNe_method='rapid', kick_vdisp=265., kick_slope=1, kick_scale=20, cbh_kwargs=None, MF_kwargs=None, meanmassdef='global', - ode_maxstep=1e10, ode_rtol=1e-7, diffcrit=1e-8, - max_mf_iter=100, mf_iter_index=0.5): + ode_maxstep=1e10, ode_rtol=1e-7, diffcrit=1e-3, + max_mf_iter=100, mf_iter_index=0.5, diffdef='rel'): import clusterbh M0 <<= u.Msun @@ -1992,7 +1994,8 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, kick_scale=kick_scale, meanmassdef=meanmassdef, ode_maxstep=ode_maxstep, ode_rtol=ode_rtol, diffcrit=diffcrit, max_mf_iter=max_mf_iter, - mf_iter_index=mf_iter_index, MF_kwargs=MF_kwargs) + mf_iter_index=mf_iter_index, MF_kwargs=MF_kwargs, + diffdef=diffdef) # reset theta to use initial values self.theta = dict(W0=W0, M0=M0.to_value('1e6 Msun'), rh0=rh0.value, From f6de62aa039fa15e23dfd22f655574b259a4a971 Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Wed, 6 May 2026 13:32:17 -0300 Subject: [PATCH 04/10] update kick plotting functions for new ssptools --- gcfit/analysis/models.py | 42 +++++++++++++++++++++++++++------------- 1 file changed, 29 insertions(+), 13 deletions(-) diff --git a/gcfit/analysis/models.py b/gcfit/analysis/models.py index ee64d89..8b1fefa 100644 --- a/gcfit/analysis/models.py +++ b/gcfit/analysis/models.py @@ -3274,7 +3274,7 @@ def plot_vesc(self, fig=None, ax=None, *, x_unit='pc', @_support_units def plot_BH_kick_fret(self, fig=None, ax=None, *, x_unit='Msun', label_position='left', verbose_label=True, - stairplot=False, **kwargs): + stairplot=True, **kwargs): r'''Plot model BH natal kick retention fraction. Plots the retention fraction of BHs caused by natal kicks in this @@ -3395,7 +3395,7 @@ def plot_BH_mass_func(self, fig=None, ax=None, *, initial=False, else: label = r'$\frac{\mathrm{d}\,N}{\mathrm{d}\,m}_{BH}$' - self._set_ylabel(ax, label, self.BH_kick_ret.unit, label_position) + self._set_ylabel(ax, label, None, label_position) self._set_xlabel(ax, r'$m_{\mathrm{BH}}$', unit=x_unit) return fig @@ -3718,6 +3718,7 @@ def __init__(self, model, observations=None): self.r = model.r self.t = [model.age] << u.Gyr + # TODO uppermost BH bin edge will often be inf, causing plot issues self._mbh_edges = np.r_[model._mf.massbins.bins.BH.lower, model._mf.massbins.bins.BH.upper[-1]] << u.Msun @@ -3993,6 +3994,9 @@ def _init_BH_dNdm(self, model): b = np.r_[BH_bins.lower, BH_bins.upper[-1]] << u.Msun bw = (BH_bins.upper - BH_bins.lower) << u.Msun + # Spline must have no inf, so just make it large (?) + bw[~np.isfinite(bw)] = 1000 << u.Msun + model_dN0dm = model._mf.Nr.BH / bw bhmf_interp = util.QuantitySpline(b[:-1] + (bw / 2), model_dN0dm, k=1) @@ -4002,14 +4006,19 @@ def _init_BH_dNdm(self, model): return bhmf_interp(mbh) def _init_kicks(self, model): - from ssptools import kicks - ks = model._mf._kick_stats - fret = kicks._get_kick_method(model._mf_kwargs['kick_method']) + fret = model._mf._kick_stats.retention << u.dimensionless_unscaled - mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) + return fret + + # from ssptools import kicks + + # ks = model._mf._kick_stats + # fret = kicks._get_kick_method(model._mf_kwargs['kick_method']) - return fret(mbh.value, **ks.parameters) << u.dimensionless_unscaled + # mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) + + # return fret(mbh.value, **ks.parameters) << u.dimensionless_unscaled class CIModelVisualizer(_ClusterVisualizer): @@ -5082,6 +5091,8 @@ def _init_BH_dNdm(self, model): b = np.r_[BH_bins.lower, BH_bins.upper[-1]] << u.Msun bw = (BH_bins.upper - BH_bins.lower) << u.Msun + bw[~np.isfinite(bw)] = 1000 << u.Msun + model_dN0dm = model._mf.Nr.BH / bw bhmf_interp = util.QuantitySpline(b[:-1] + (bw / 2), model_dN0dm, k=1) @@ -5091,18 +5102,20 @@ def _init_BH_dNdm(self, model): return bhmf_interp(mbh) def _init_kicks(self, model): - from ssptools import kicks + # from ssptools import kicks # This holds nans wherever kicks are not actually done (e.g. 0 BH bins) - # model_ret = model._mf._kick_stats['retention'] + fret = model._mf._kick_stats.retention << u.dimensionless_unscaled # So instead, recompute the kicks (which are really fast) - ks = model._mf._kick_stats - fret = kicks._get_kick_method(model._mf_kwargs['kick_method']) + # ks = model._mf._kick_stats + # fret = kicks._get_kick_method(model._mf_kwargs['kick_method']) - mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) + # mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) - return fret(mbh.value, **ks.parameters) << u.dimensionless_unscaled + # return fret(mbh.value, **ks.parameters) << u.dimensionless_unscaled + + return fret # ---------------------------------------------------------------------- # Save and load confidence intervals to a file @@ -5712,10 +5725,13 @@ def from_theta(cls, theta, observations, model_params): def _init_BH_dN0dm(self, model): + # clusterbh ibh should match mf ibh, but with natal kicks applied BH_bins = model._clusterbh.ibh.bins b = np.r_[BH_bins.lower, BH_bins.upper[-1]] << u.Msun bw = (BH_bins.upper - BH_bins.lower) << u.Msun + bw[~np.isfinite(bw)] = 1000 << u.Msun + model_dN0dm = model._clusterbh.ibh.N / bw bhmf_interp = util.QuantitySpline(b[:-1] + (bw / 2), model_dN0dm, k=1) From 17dd2841afd2ec523e64a751e7803fcff3b60f2a Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Wed, 6 May 2026 13:33:41 -0300 Subject: [PATCH 05/10] add BH IFMR options back in These must be explicitly arguments, not just in `MF_kwargs` because they must be passed to clusterBH as well --- gcfit/core/data.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/gcfit/core/data.py b/gcfit/core/data.py index 36de951..2745d2f 100644 --- a/gcfit/core/data.py +++ b/gcfit/core/data.py @@ -1108,8 +1108,10 @@ def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, meq=0.0, eta=0.0, zeta=1.0, esc_rate=0.0, natal_kicks=True, kick_method='maxwellian', f_kick=None, SNe_method='rapid', vesc=90, kick_vdisp=265., - kick_slope=1, kick_scale=20, MF_kwargs=None, - meanmassdef='global', ode_maxstep=1e10, ode_rtol=1e-7, + kick_slope=1, kick_scale=20, + BH_IFMR_method='banerjee20', BH_IFMR_kwargs=None, + MF_kwargs=None, meanmassdef='global', + ode_maxstep=1e10, ode_rtol=1e-7, diffcrit=1e-3, max_mf_iter=100, mf_iter_index=0.5, diffdef='rel'): @@ -1184,7 +1186,8 @@ def __init__(self, W0, M, rh, g=1.5, delta=0.45, ra=1e8, NS_ret, BH_ret_dyn, natal_kicks, self.vesc0, kick_method, f_kick, SNe_method, kick_vdisp, - kick_slope, kick_scale, **MF_kwargs) + kick_slope, kick_scale, + BH_IFMR_method, BH_IFMR_kwargs, **MF_kwargs) if not self._mf.converged: mssg = ("Mass function evolution ODE failed to converge" @@ -1760,7 +1763,8 @@ class EvolvedModel(Model): def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, NS_ret, BH_ret_dyn, natal_kicks, vesc, kick_method, f_kick, SNe_method, kick_vdisp, - kick_slope, kick_scale, **kwargs): + kick_slope, kick_scale, BH_IFMR_method, BH_IFMR_kwargs, + **kwargs): '''Alternative MF init using prior-computed IMF and clusterBH outputs''' from ssptools import EvolvedMFWithBH @@ -1784,6 +1788,8 @@ def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, kick_vdisp=kick_vdisp, kick_slope=kick_slope, kick_scale=kick_scale, + BH_IFMR_method=BH_IFMR_method, + BH_IFMR_kwargs=BH_IFMR_kwargs, **kwargs # will error here if MF_kwargs included any of above args ) @@ -1798,6 +1804,7 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, natal_kicks=True, kick_method='maxwellian', f_kick=None, SNe_method='rapid', kick_vdisp=265., kick_slope=1, kick_scale=20, + BH_IFMR_method='banerjee20', BH_IFMR_kwargs=None, cbh_kwargs=None, MF_kwargs=None, meanmassdef='global', ode_maxstep=1e10, ode_rtol=1e-7, diffcrit=1e-3, max_mf_iter=100, mf_iter_index=0.5, diffdef='rel'): @@ -1821,9 +1828,12 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, # Make sure the relevant kickparams are passed to clusterBH by default # but still respect any explicitly passed in `ibh_kwargs` too - ibh_kwargs = dict(kick_method=kick_method, f_kick=f_kick, - SNe_method=SNe_method, kick_vdisp=kick_vdisp, - kick_slope=kick_slope, kick_scale=kick_scale) + ibh_kwargs = dict( + kick_method=kick_method, f_kick=f_kick, + SNe_method=SNe_method, kick_vdisp=kick_vdisp, + kick_slope=kick_slope, kick_scale=kick_scale, + BH_IFMR_method=BH_IFMR_method, BH_IFMR_kwargs=BH_IFMR_kwargs + ) ibh_kwargs |= cbh_kwargs.get('ibh_kwargs', {}).copy() From 119fdd1cf106804e7e446a79fe471450a6ffb9fc Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Wed, 6 May 2026 13:36:04 -0300 Subject: [PATCH 06/10] switch from clusterBH to cBHBd library The underlying models should be the exact same, for the most part, but the cBHBd package will be the stable package. --- gcfit/core/data.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/gcfit/core/data.py b/gcfit/core/data.py index 2745d2f..3775c8c 100644 --- a/gcfit/core/data.py +++ b/gcfit/core/data.py @@ -1808,7 +1808,8 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, cbh_kwargs=None, MF_kwargs=None, meanmassdef='global', ode_maxstep=1e10, ode_rtol=1e-7, diffcrit=1e-3, max_mf_iter=100, mf_iter_index=0.5, diffdef='rel'): - import clusterbh + + import cbhbd M0 <<= u.Msun rh0 <<= u.pc @@ -1924,6 +1925,7 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, # clusterBH fit parameters should use defaults, or given in cbh_kwargs # ------------------------------------------------------------------ + cbh_kwargs.setdefault('dtout', 2.0) cbh_kwargs.setdefault('ssp', True) cbh_kwargs.setdefault('kick', True) cbh_kwargs.setdefault('tidal', True) @@ -1944,8 +1946,9 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, self.cbh_kwargs = cbh_kwargs - self._clusterbh = clusterbh.clusterBH(N0, self.rhoh0.value, - **self.cbh_kwargs) + self._clusterbh = cbhbd.cbhbd.CBHBD(N=N0, rhoh0=self.rhoh0.value, + compute_mergers=False, + **self.cbh_kwargs).cluster # Make sure no negative f_BH values are allowed self._clusterbh.fbh[self._clusterbh.fbh < 0] = 0. @@ -1971,7 +1974,8 @@ def __init__(self, W0, M0, rh0, g=1.5, delta=0.45, ra=1e8, / self._clusterbh.tev) # Ejection - bf = self._clusterbh.balance_function(self._clusterbh.t) + bf = self._clusterbh.balance_function((self._clusterbh.t * 1e3) - tcc) + alpha_c = (self._clusterbh.alpha_ci * bf) alpha_c += ((self._clusterbh.alpha_cf * bf - alpha_c) * (1 - self._clusterbh.beta_function(self._clusterbh.S))) From 0d3541d2bdea756030f3f714734289aeb768e87c Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Wed, 6 May 2026 13:39:14 -0300 Subject: [PATCH 07/10] remove some commented code --- gcfit/analysis/models.py | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/gcfit/analysis/models.py b/gcfit/analysis/models.py index 8b1fefa..604c4d4 100644 --- a/gcfit/analysis/models.py +++ b/gcfit/analysis/models.py @@ -4011,15 +4011,6 @@ def _init_kicks(self, model): return fret - # from ssptools import kicks - - # ks = model._mf._kick_stats - # fret = kicks._get_kick_method(model._mf_kwargs['kick_method']) - - # mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) - - # return fret(mbh.value, **ks.parameters) << u.dimensionless_unscaled - class CIModelVisualizer(_ClusterVisualizer): '''Analysis and visualization of a model, with confidence intervals. @@ -5102,19 +5093,9 @@ def _init_BH_dNdm(self, model): return bhmf_interp(mbh) def _init_kicks(self, model): - # from ssptools import kicks - # This holds nans wherever kicks are not actually done (e.g. 0 BH bins) fret = model._mf._kick_stats.retention << u.dimensionless_unscaled - # So instead, recompute the kicks (which are really fast) - # ks = model._mf._kick_stats - # fret = kicks._get_kick_method(model._mf_kwargs['kick_method']) - - # mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) - - # return fret(mbh.value, **ks.parameters) << u.dimensionless_unscaled - return fret # ---------------------------------------------------------------------- From 36c6ec707fe227adb1db442d058c6a8294874da1 Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Wed, 6 May 2026 14:01:57 -0300 Subject: [PATCH 08/10] make sure IFMR args are passed to base class --- gcfit/core/data.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gcfit/core/data.py b/gcfit/core/data.py index 3775c8c..45ee1d2 100644 --- a/gcfit/core/data.py +++ b/gcfit/core/data.py @@ -994,7 +994,8 @@ def __str__(self): def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, NS_ret, BH_ret_dyn, natal_kicks, vesc, kick_method, f_kick, SNe_method, kick_vdisp, - kick_slope, kick_scale, **kwargs): + kick_slope, kick_scale, BH_IFMR_method, BH_IFMR_kwargs, + **kwargs): '''Compute an evolved mass function using `ssptools.EvolvedMF`''' # Total mass of this will be wrong due to N0 but Mj is scaled in limepy @@ -1019,6 +1020,8 @@ def _evolve_mf(self, m_breaks, a1, a2, a3, nbins, FeH, age, esc_rate, tcc, kick_vdisp=kick_vdisp, kick_slope=kick_slope, kick_scale=kick_scale, + BH_IFMR_method=BH_IFMR_method, + BH_IFMR_kwargs=BH_IFMR_kwargs, **kwargs # will error here if MF_kwargs included any of above args ) From cae5bf9674273f8e3a395c641fe5a19522c52cf6 Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Wed, 6 May 2026 14:02:59 -0300 Subject: [PATCH 09/10] fix bugs around M scaling in static model plots The non-evolved model does not know N0, so the `EvolvedMF` it uses is not normalized to the correct mass, which only happens in limepy. So whenever `model._mf` is used in static models, the overall mass needs to have this scale factor applied again. --- gcfit/analysis/models.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/gcfit/analysis/models.py b/gcfit/analysis/models.py index 604c4d4..ef68b1e 100644 --- a/gcfit/analysis/models.py +++ b/gcfit/analysis/models.py @@ -3780,7 +3780,8 @@ def __init__(self, model, observations=None): self.N_WD = model.WD.Nj.sum() self.BH_massfunc = self.BH0_massfunc = self._init_BH_dNdm(model)[bh_slc] self.BH_kick_ret = self._init_kicks(model)[bh_slc] - self.M_kicked = model._mf._kick_stats.total_kicked << u.Msun + Mscale = model._MS / model._mf.M.sum() # non-ev models need mf scaled + self.M_kicked = (model._mf._kick_stats.total_kicked * Mscale) << u.Msun self.Ms_t = model.nonBH.Mj.sum()[t_slc] self.mmean_t = model.mmean[t_slc] self.rt_t = model.rt[t_slc] @@ -3997,9 +3998,12 @@ def _init_BH_dNdm(self, model): # Spline must have no inf, so just make it large (?) bw[~np.isfinite(bw)] = 1000 << u.Msun - model_dN0dm = model._mf.Nr.BH / bw + model_dNdm = model._mf.Nr.BH / bw - bhmf_interp = util.QuantitySpline(b[:-1] + (bw / 2), model_dN0dm, k=1) + # Sometimes nans sneak around + model_dNdm[~np.isfinite(model_dNdm)] = 0. << model_dNdm.unit + + bhmf_interp = util.QuantitySpline(b[:-1] + (bw / 2), model_dNdm, k=1) mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) @@ -4698,7 +4702,9 @@ def from_chain(cls, chain, observations, model_params, N=100, *, M_BH0[model_ind] = M_BH[model_ind] N_BH0[model_ind] = N_BH[model_ind] - M_kicked[model_ind] = model._mf._kick_stats.total_kicked << u.Msun + Mscale = model._MS / model._mf.M.sum() # scaling for non-ev eMFs + ks = model._mf._kick_stats + M_kicked[model_ind] = (ks.total_kicked * Mscale) << u.Msun bhslc = (slice(None), model_ind, 0) BH_massfunc[bhslc] = BH0_massfunc[bhslc] = viz._init_BH_dNdm(model) @@ -5084,9 +5090,12 @@ def _init_BH_dNdm(self, model): bw[~np.isfinite(bw)] = 1000 << u.Msun - model_dN0dm = model._mf.Nr.BH / bw + model_dNdm = model._mf.Nr.BH / bw - bhmf_interp = util.QuantitySpline(b[:-1] + (bw / 2), model_dN0dm, k=1) + # Sometimes nans sneak around + model_dNdm[~np.isfinite(model_dNdm)] = 0. << model_dNdm.unit + + bhmf_interp = util.QuantitySpline(b[:-1] + (bw / 2), model_dNdm, k=1) mbh = 0.5 * (self._mbh_edges[1:] + self._mbh_edges[:-1]) From 64722223f72e235945ad00ef51a91fd0905c569f Mon Sep 17 00:00:00 2001 From: Nolan Dickson Date: Sun, 7 Jun 2026 12:31:22 -0300 Subject: [PATCH 10/10] update rtfd build versions --- .readthedocs.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 9607f0a..be1954a 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -1,9 +1,9 @@ version: 2 build: - os: "ubuntu-22.04" + os: "ubuntu-24.04" tools: - python: "3.9" + python: "3.14" # Build from the docs/ directory with Sphinx sphinx: