Add Offner et al. 2023 multiplicity models (unresolved + resolved) by jluastro · Pull Request #160 · MovingUniverseLab/SPISEA

jluastro · 2026-08-13T16:36:08Z

Adds opt-in Offner et al. 2023 (PPVII) multiplicity for SPISEA. Default cluster generation is unchanged (MultiplicityUnresolved / MultiplicityResolvedDK). There is no MultiplicityOffner2023 alias: pick unresolved or resolved explicitly.

Models

MultiplicityUnresolvedOffner2023: Table 1 logistic MF/CSF and error-weighted logistic mass-ratio index γ(M).
MultiplicityResolvedOffner2023: same plus Table 2 μ(a) (smooth broken power law) and σ(log10 a) (2-param logistic). Does not subclass MultiplicityResolvedDK.
Generic bases: MultiplicityLogistic (MF/CSF and γ(M) logistics) and MultiplicityPiecewisePowerLaw.

API

Companion masses and “is multiple?” are drawn in IMF.generate_cluster / IMF.calc_multi.
random_q(x, mass=None) is mass-aware.
Resolved orbits are duck-typed (log_semimajoraxis / random_is_hard).
Both Offner class __init__s list the same model-shape kwargs (MF_*, CSF_*, q_*, a_*, sig_*) with the numbers hardcoded in that function. Resolved also has sep_sig / sep_sig_mass and passes the shape kwargs to super().__init__(...) explicitly. No OFFNER2023_* module globals.
H_BURNING_MASS (0.08 Msun) remains only for the Fontanive bd_q_power γ step on MultiplicityUnresolved.q_power_at_mass. It is not a companion-count policy.
log_semimajoraxis(mass, rng=None) and random_keplarian_parameters(x, y, z, rng=None) take a numpy.random.Generator. If rng is omitted they use a new numpy.random.default_rng(). synthetic.py passes the cluster self.rng into both. Keplerian sign draws use rng.choice, not stdlib random.choice. Semimajor-axis draws pass random_state=rng to truncnorm.rvs.

Companion counts

Higher-order multiples (triples+) are allowed at all masses, including brown dwarfs.
There is no binary_only_mass_max and no CSF=MF / n_comp≤1 mass cut.
companion_max=True is the only companion-count cap: it clips at CSF_max across the entire mass range, in every model.

Docs / figures

docs/multiplicity.rst is a short landing page (wiring + autoclass + figure list). Functional forms live in class/method docstrings.
Comparison figures vs SPISEA v2.5, generated from multiplicity objects (not hardcoded formulas), with Offner Table 1/2 points:
- docs/figures/plot_mf_offner_vs_spisea2.5.py → mf_offner_vs_spisea2.5.png
- docs/figures/plot_csf_offner_vs_spisea2.5.py → csf_offner_vs_spisea2.5.png
- docs/figures/plot_q_sep_offner_vs_spisea2.5.py regenerates the last four PNGs (q, sep, sig_loga, meanq).
Duchêne & Kraus public methods on MultiplicityResolvedDK: log_a_mean, a_mean, sigma_log_a.

Tests

Same-seed default_rng reproducibility for random_companion_count, draw_n_companions, draw_companion_masses (companion assignment), log_semimajoraxis, and random_keplarian_parameters.
pytest --noconftest spisea/tests/test_multiplicity.py -k "not test_resolvedmult" .
test_resolvedmult still needs the full photometry stack (synphot).

Introduce MultiplicityPiecewisePowerLaw plus unresolved/resolved Offner 2023 classes fitted to Table 1 MF/CF (including brown dwarfs). Companion mass and separation draws now live on the multiplicity object; Lu+2013 defaults are unchanged. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

calc_multi now only delegates to the multiplicity object. random_q is mass-aware so Offner γ_trunc applies to BD primaries instead of the hardcoded Fontanive 6.1. random_companion_count owns the BD binaries-only cap. synthetic.py duck-types resolved orbits on log_semimajoraxis, random_e, and random_keplarian_parameters. Docs and IMF tests updated. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Replace the discontinuous 8-segment two-point Table 1 fit with a 3-segment broken power law continuous at 0.08 and 1.5 Msun. Add a two-panel MF vs mass figure (Lu+2013 vs Offner vs Table 1) to the docs and PR. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Offner et al. 2023 now subclasses MultiplicityLogistic instead of the piecewise power law. MF/CSF use equal-weight Table 1 logistic coefficients; the generic piecewise class remains for tests and other surveys. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Two-panel BD-zoom plus full-range plots of Table 1 gamma_trunc and characteristic a, plus a mean-q companion panel. Curves come from the multiplicity objects so they cannot drift from the code. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Two-panel plot of Table 2 separation scatter versus the Duchene-Kraus linear-in-log-M fit. Offner holds sigma=0.7 for brown dwarfs; Lu DK shows the 0.08 Msun blend dip. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Replace Table 1/2 log-mass interpolations with an error-weighted logistic for gamma, a logcosh smooth broken power law for mu(a), and a 2-parameter logistic for sigma(log10 a). Resolved draws use those as loc and scale. MF/CSF logistic and Lu+2013 are unchanged. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Expand docs/multiplicity.rst with formulas for MF/CSF, gamma, mu(a), and sigma(log a), plus the five Lu+2013 comparison figures. Point imf.rst and the changelog at the opt-in Offner model. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

The comparison plots and docs describe Offner 2023 vs the default MultiplicityUnresolved / MultiplicityResolvedDK that shipped in SPISEA v2.5. Keep Lu et al. 2013 citations on the original classes. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Filenames, RST figure paths, and plot-script helpers now match the SPISEA v2.5 comparison baseline. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Updated the documentation for the multiplicity object, clarifying its functions and usage in the IMF. Revised sections on companion evolution and recommended multiplicity classes.

Add Parameters and Returns (types and units) to every new helper, class, and method in multiplicity.py. Expand modified methods that already had a Parameters block but omitted units. No math changes. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Private helpers sit after the public classes and MultiplicityOffner2023 alias so Sphinx/source order documents the classes first. Table arrays that call _offner2023_table1_geom_mass follow that helper. No behavior change. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Evaluate the logistic and smooth broken power law on positive masses only. Keep M>0 math unchanged. Tests and docs no longer treat M<=0 as mapping to the low-mass asymptote. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

hustonm · 2026-08-19T17:59:54Z

same comment as the mf plot

hustonm · 2026-08-19T18:08:22Z

+SPISEA v2.5 uses :math:`\mathrm{MF} = 0.44\,M^{0.51}` (clipped
+to 1) for arrays, plus a scalar-only brown-dwarf staircase
+(0 / 8% / 16%). Cluster generation on that class still uses the
+stellar power law for brown-dwarf primaries.


Not clear what "Cluster generation on that class still uses the stellar power law for brown-dwarf primaries." means directly after saying the mass function differs for BDs.

hustonm · 2026-08-19T18:09:58Z

+
+with :math:`(A, B, M_0, k) = (6.6, -1.77, 0.0651, 0.629)`. Call
+``q_power_at_mass(mass)`` or ``random_q(x, mass=...)``. Without
+``mass``, ``random_q(x)`` keeps the historical stellar-only power law.


Do we want that to be allowed? I don't know why we would want the Offner model to be able to use a different q distribution.

hustonm · 2026-08-19T18:11:05Z

+The err-weighted fit undershoots Fontanive et al. (2018)
+:math:`8\pm 6\%` MF and :math:`\gamma = 4.8\pm 2.2`: at
+:math:`0.033\,M_\odot`, :math:`\gamma \approx 3.3`. That is the
+fit, not a bug. SPISEA v2.5 is a step: :math:`\gamma = 6.1` for


It might be confusing that we say SPISEA v2.5 throughout rather than specifying the specific older Multiplicity classes that do this.

hustonm · 2026-08-19T18:13:25Z

+The comparison figures are generated from the multiplicity object
+methods (``multiplicity_fraction``, ``q_power_at_mass``,
+``a_mean``, ``sigma_log_a``) so they cannot drift from the code.
+From the repository root::


They don't actually use the older Multiplicity models for the comparison though

hustonm · 2026-08-19T18:14:41Z

+			 sigma_log_a
+
+
+Offner et al. 2023 multiplicity


How much of this overlaps with the documentation that already exists within the module? It's a bit weird that we have so much extra info here for this model but not for the other Multiplicity classes. Maybe we should just have more info in all of their docstrings so the auto-doc can handle it

hustonm · 2026-08-19T18:16:26Z

+# Equal-weight logistic-in-log-mass fit to Offner et al. 2023 Table 1
+# geom-mean (M, MF) and (M, CF) points:
+#   y(M) = A + (B - A) / (1 + (M / M0)**(-k))
+OFFNER2023_MF_A = 0.14


We had talked about potentially just putting the numbers in as defaults for the module instead of all these variables. That reads cleaner to me, and I think would be more consistent with the handling of the other modules.

hustonm · 2026-08-19T18:18:44Z

+        """
+        mass_arr = np.atleast_1d(np.asarray(mass, dtype=float))
+        q_pow = np.full(mass_arr.shape, self.q_pow, dtype=float)
+        q_pow[mass_arr <= self.binary_only_mass_max] = FONTANIVE2018_BD_Q_POWER


This could get weird if someone wants to use a cutoff other than the BD mass limit for the binary only cut.

hustonm · 2026-08-19T18:19:53Z

-        return  q
+        if mass is None:
+            return _q_from_powerlaw(x, self.q_pow, self.q_min)
+        return _q_from_powerlaw(x, self.q_power_at_mass(mass), self.q_min)


We should probably force it to use mass in this calculation? That's easy for me to do in my refactor if we want.

hustonm · 2026-08-19T18:21:20Z

+            [q_min, 1]. Shape (len(prim_subset), n_comp).
+        """
+        q_values = np.empty((len(prim_subset), n_comp))
+        bd_mask = prim_subset <= self.binary_only_mass_max


again, this is equating the binary_only_mass_max with the brown dwarf mass limit, which I thought we were not locking in as equivalent

hustonm · 2026-08-19T18:22:31Z

+            mass_arr, self.mass_limits, self.MF_amps, self.MF_powers,
+            clip_min=0.0, clip_max=1.0)
+        csf = np.maximum(csf, mf)
+        bd = mass_arr <= self.binary_only_mass_max


same concern about BD mass limit versus binary_only limit

hustonm · 2026-08-19T18:22:52Z

+        csf = np.maximum(csf, mf)
+
+        # Fix all brown dwarf binaries so they only have one companion
+        bd = mass_arr <= self.binary_only_mass_max


BD versus binary-only mass limit conflation

hustonm · 2026-08-19T18:24:35Z

+            CSF_B=OFFNER2023_CSF_B,
+            CSF_M0=OFFNER2023_CSF_M0,
+            CSF_k=OFFNER2023_CSF_K,
+            CSF_max=CSF_max, q_power=q_power, q_min=q_min,


why is there q_power and q_min here? shouldn't it be the the logistic parameters?

hustonm · 2026-08-19T18:24:55Z

+            Mass-ratio power-law index γ, dimensionless.
+            Python float if ``mass`` is scalar, ndarray otherwise.
+        """
+        gamma = _logistic_in_logm(mass, OFFNER2023_Q_A, OFFNER2023_Q_B, OFFNER2023_Q_M0, OFFNER2023_Q_K)


these should probably be optional kwargs with defaults, right, not hard-coded?

hustonm · 2026-08-19T18:25:28Z

+        """
+        # Calculate the characteristic log10(a / 1 AU) using a smooth broken power law 
+        log_a_mean = _smooth_broken_loglog(
+            mass, OFFNER2023_A_MUP, OFFNER2023_A_MP,


again, do we want these hard-coded or no? I think either is justifiable, but it's a bit inconsistent throughout

hustonm · 2026-08-19T18:26:31Z

-        return inclination, Omega, omega
+        return gamma
+
+    def log_a_mean(self, mass):


Why are there orbital parameters in the Unresolved class?

hustonm · 2026-08-19T18:26:56Z

+        return log_a
+
+# Convenience alias; unresolved Table 1 model is the usual opt-in object.
+MultiplicityOffner2023 = MultiplicityUnresolvedOffner2023


I think we should remove this, like we talked about

hustonm · 2026-08-19T18:31:58Z

+    return out
+
+
+def _offner2023_table1_geom_mass(m_lo, m_hi):


this is just a geometric mean, we might make the name clearer that it's generic

hustonm · 2026-08-19T18:33:37Z

+            comps = compMass[i].compressed()
+            if len(comps):
+                q_bd.extend(list(comps / mass[i]))
+        if len(q_bd) >= 5:


If we set a seed, we can hard-code in a case where these are guaranteed to be properly tested

hustonm · 2026-08-19T18:37:15Z

+    # Force multiples so we test the count draw, not the MF coin flip.
+    n_comp = multi.draw_n_companions(masses, csf, mf, rng)
+    assert np.all(n_comp <= 1)
+    assert np.all(n_comp >= 1)


Confusing asserts here. Shouldn't a BD allowed to have no companions?

hustonm · 2026-08-19T18:48:44Z

+    assert np.abs(500.0 - sys_mass.sum()) < 500.0 * 0.05
+
+
+def test_calc_multi_uses_multiplicity_q_and_counts():


Do we want all these tests to be active on the user side, or was a lot of this to keep the AI in check?

hustonm · 2026-08-19T18:49:58Z

-            companions.add_column(Column(self.imf._multi_props.log_semimajoraxis(star_systems['mass'][companions['system_idx']]), name='log_a'))
-            companions.add_column(Column(self.imf._multi_props.random_e(self.rng.random(N_comp_tot)), name='e'))
-            companions['i'], companions['Omega'], companions['omega'] = self.imf._multi_props.random_keplarian_parameters(
+        # Duck-type resolved multiplicity: any object with orbital methods


The comparison figure now evaluates the original object: the array path for the dashed power law and the scalar path for the BD staircase. Hardcoded 0.44 M**0.51 and 0/8%/16% helpers are gone. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Expose MultiplicityResolvedDK.log_a_mean, a_mean, and sigma_log_a as the characteristic mean and width already used by log_semimajoraxis, and generate each comparison PNG from a matching script that calls the multiplicity objects rather than hardcoded functional forms. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Users must construct MultiplicityUnresolvedOffner2023 or MultiplicityResolvedOffner2023 explicitly; there is no shorthand. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

docs/multiplicity.rst is now a short wiring intro plus autoclass and one-line comparison figures. Formulas, coefficients, BD policy, and Table 1 caveats live on the public multiplicity classes. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Keep docs/multiplicity.rst as a short landing page plus autoclass and the five comparison figures. List the plot_* scripts used to regenerate them. Align Piecewise/Logistic Notes with the other classes. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Drop leftover piecewise-fit helper and Table 1/2 arrays that were stored on the objects but never read. Keep sep_sig_mass / sep_sig for the σ comparison plot. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Delete OFFNER2023_* and FONTANIVE2018_BD_Q_POWER module globals. MultiplicityUnresolved takes bd_q_power=6.1; Offner classes take the Table 1/2 logistic and smooth-break values as kwargs, stored on self and used by q_power_at_mass / log_a_mean / sigma_log_a. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

List the Offner shape kwargs with hardcoded defaults on both unresolved and resolved __init__s. Drop binary_only_mass_max so companion_max is the only count cap. Add the CSF vs mass comparison figure from companion_star_fraction on the multiplicity objects. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

log_semimajoraxis and random_keplarian_parameters now take rng=None and default to a new numpy Generator. synthetic.py passes self.rng. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Five tests check that the same default_rng seed reproduces random_companion_count, draw_n_companions, draw_companion_masses, log_semimajoraxis, and random_keplarian_parameters. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

…ingUniverseLab/SPISEA into cursor/offner-2023-multiplicity-ad2c

jluastro marked this pull request as ready for review August 13, 2026 16:36

cursoragent and others added 15 commits August 13, 2026 16:38

Add Offner vs Lu+2013 sigma(log a) comparison figure

8ba5a50

Two-panel plot of Table 2 separation scatter versus the Duchene-Kraus linear-in-log-M fit. Offner holds sigma=0.7 for brown dwarfs; Lu DK shows the 0.08 Msun blend dip. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Rename Offner comparison figures from lu2013 to spisea2.5

b98efcf

Filenames, RST figure paths, and plot-script helpers now match the SPISEA v2.5 comparison baseline. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Revise multiplicity documentation and clarify usage

13d722f

Updated the documentation for the multiplicity object, clarifying its functions and usage in the IMF. Revised sections on companion evolution and recommended multiplicity classes.

Code cleanup (spacing, etc.).

48eafdc

Drop massless-primary (M<=0) support in multiplicity helpers

af8eb8e

Evaluate the logistic and smooth broken power law on positive masses only. Keep M>0 math unchanged. Tests and docs no longer treat M<=0 as mapping to the low-mass asymptote. Co-authored-by: Jessica Lu <jlu.astro@berkeley.edu>

Reviewed code, cleaned up formatting, added more documentation.

a66554b