COSMIC-PopSynth · katiebreivik · Aug 8, 2024 · Jun 21, 2024 · Jul 15, 2024 · Aug 7, 2024
diff --git a/.github/workflows/build_wheels_and_publish.yml b/.github/workflows/build_wheels_and_publish.yml
@@ -7,10 +7,11 @@ on:
 
 
 env:
-  CIBW_BUILD: "cp37-* cp38-* cp39-* cp310-*"
+  CIBW_BUILD: "cp39-* cp310-*"
   CIBW_SKIP: "*-win32 *musllinux*"
   CIBW_MANYLINUX_X86_64_IMAGE: manylinux2014
   CIBW_BEFORE_BUILD: pip install numpy --config-settings=setup-args="-Dallow-noblas=true"
+  CIBW_BUILD_VERBOSITY: "1"
 
 jobs:
   build-wheels:
@@ -27,15 +28,12 @@ jobs:
     - uses: actions/setup-python@v4
       name: Install Python
       with:
-        python-version: "3.7"
+        python-version: "3.10"
 
     - name: Setup Mac
       if: runner.os == 'macOS'
       run: |
         python -m pip install numpy h5py versioneer
-        sudo ln -s /usr/local/bin/gfortran-11 /usr/local/bin/gfortran
-        sudo mkdir /usr/local/gfortran
-        sudo ln -s /usr/local/Cellar/gcc@9/9.3.0_1/lib/gcc/9 /usr/local/gfortran/lib
 
     - name: Setup Linux
       if: runner.os == 'Linux'
@@ -44,8 +42,11 @@ jobs:
 
     - name: Build wheels
       run: |
-        python -m pip install cibuildwheel
-        python -m cibuildwheel --output-dir wheelhouse
+        sudo ln -s /opt/homebrew/bin/gfortran-11 /opt/homebrew/bin/gfortran || true
+        python -m pip install meson ninja build
+        meson setup build
+        meson compile -C build
+        python -m build --wheel --outdir wheelhouse 
     - uses: actions/upload-artifact@v3
       with:
         path: ./wheelhouse/*.whl
@@ -60,16 +61,17 @@ jobs:
       - uses: actions/setup-python@v4
         name: Install Python
         with:
-          python-version: "3.7" 
+          python-version: "3.10" 
 
       - name: Build dist
         run: |
           python -m pip install --upgrade pip
-          pip install setuptools wheel twine numpy
-          python setup.py sdist
+          pip install setuptools wheel twine numpy build meson ninja
+          pip install .
+          python -m build
       - uses: actions/upload-artifact@v3
         with:
-          path: dist/*.tar.gz
+          path: ./dist/*.tar.gz
 
   publish:
     needs: [build-wheels, build-source-dist]
@@ -78,7 +80,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        python-version: [3.8]
+        python-version: ["3.10"]
 
     steps:
     - uses: actions/checkout@v2

diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
@@ -12,7 +12,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ['3.7', '3.8', '3.9', '3.10']
+        python-version: ['3.9', '3.10']
 
     steps:
     - uses: actions/checkout@v2

diff --git a/bin/cosmic-pop b/bin/cosmic-pop
@@ -101,7 +101,7 @@ def parse_commandline():
     parser.add_argument("--binary_state", nargs='+', type=int)
     parser.add_argument("--sampling_method")
     parser.add_argument("--primary_model", help="Chooses the initial primary mass function from: salpeter55, kroupa93, kroupa01", type=str)
-    parser.add_argument("--binfrac_model", help="Chooses the binary fraction model from: a float between [0,1] and vanHaaften", type=float)
+    parser.add_argument("--binfrac_model", help="Chooses the binary fraction model from: a float between [0,1], vanHaaften, and offner22", type=float)
     parser.add_argument("--ecc_model", help="Chooses the initial eccentricity distribution model from: thermal, uniform, and sana12", type=str)
     parser.add_argument("--porb_model", help="Chooses the initial orbital period distribution model from: log_uniform and sana12", type=str)
     parser.add_argument("--SF_start", help="Sets the time in the past when star formation initiates in Myr", type=float)

diff --git a/meson.build b/meson.build
@@ -1,7 +1,7 @@
 project('cosmic', 
     'c',
     'fortran',
-    version : '3.4.10',
+    version : '3.4.13',
     default_options: ['warning_level=0', 'optimization=3'],
 )
 

diff --git a/pyproject.toml b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
-name = "cosmic"
-distname = "cosmic-popsynth"
+name = "cosmic-popsynth"
+distname = "cosmic"
 description = "a Python-interfaced binary population synthesis suite"
 authors = [
     { name="Katelyn Breivik", email="katie.breivik@gmail.com" },
@@ -11,7 +11,7 @@ authors = [
 readme = "README.md"
 
 
-requires-python = ">=3.7"
+requires-python = ">=3.9"
 license = { text = "MIT License" }
 classifiers = [
     "Development Status :: 5 - Production/Stable",

diff --git a/requirements.txt b/requirements.txt
@@ -1,11 +1,11 @@
 scipy >= 0.12.1
-numpy >= 1.16
+numpy == 1.23.5
 astropy >= 1.1.1
 configparser
 tqdm >= 4.0
 pandas >= 0.24
 tables > 3.5.0
-h5py >= 1.3
+h5py >= 3.7.0
 schwimmbad >= 0.3.1
 matplotlib >= 2.0.0
 sphinx >= 1.6.1

diff --git a/src/cosmic/_version.py b/src/cosmic/_version.py
@@ -1 +1 @@
-__version__ = "3.4.10"
+__version__ = "3.4.13"
diff --git a/src/cosmic/sample/sampler/cmc.py b/src/cosmic/sample/sampler/cmc.py
@@ -83,7 +83,7 @@ def get_cmc_sampler(
         Sets the pairing of stars M>msort only with stars with M>msort
 
     binfrac_model : `str or float`
-        Model for binary fraction; choices include: vanHaaften or a fraction where 1.0 is 100% binaries
+        Model for binary fraction; choices include: vanHaaften, offner22, or a fraction where 1.0 is 100% binaries
 
     binfrac_model_msort : `str or float`
         Same as binfrac_model for M>msort

diff --git a/src/cosmic/sample/sampler/independent.py b/src/cosmic/sample/sampler/independent.py
@@ -114,7 +114,7 @@ def get_independent_sampler(
         Duration of constant star formation beginning from SF_Start in Myr
 
     binfrac_model : `str or float`
-        Model for binary fraction; choices include: vanHaaften or a fraction where 1.0 is 100% binaries
+        Model for binary fraction; choices include: vanHaaften, offner22, or a fraction where 1.0 is 100% binaries
 
     binfrac_model_msort : `str or float`
         Same as binfrac_model for M>msort
@@ -601,6 +601,7 @@ def binary_select(self, primary_mass, binfrac_model=0.5, **kwargs):
         either a binary fraction specified by a float or a
         primary-mass dependent binary fraction following
         `van Haaften et al.(2009) <http://adsabs.harvard.edu/abs/2013A%26A...552A..69V>`_ in appdx
+        or `Offner et al.(2022) <https://arxiv.org/abs/2203.10066>`_ in fig 1
 
         Parameters
         ----------
@@ -609,6 +610,9 @@ def binary_select(self, primary_mass, binfrac_model=0.5, **kwargs):
 
             binfrac_model : str or float
                 vanHaaften - primary mass dependent and ONLY VALID up to 100 Msun
+
+                offner22 - primary mass dependent
+
                 float - fraction of binaries; 0.5 means 2 in 3 stars are a binary pair while 1
                 means every star is in a binary pair
 
@@ -650,13 +654,35 @@ def binary_select(self, primary_mass, binfrac_model=0.5, **kwargs):
                 binary_choose_low = np.random.uniform(
                     0, 1.0, primary_mass[lowmassIdx].size)
 
+                (singleIdx_low,) = np.where(
+                    binary_fraction_low < binary_choose_low)
+                (binaryIdx_low,) = np.where(
+                    binary_fraction_low >= binary_choose_low)
+            elif binfrac_model == "offner22":
+                from scipy.interpolate import BSpline
+                t = [0.0331963853, 0.0331963853, 0.0331963853, 0.0331963853, 0.106066017,
+                     0.212132034, 0.424264069, 0.866025404, 1.03077641, 1.11803399,
+                     1.95959179, 3.87298335, 6.32455532, 11.6619038, 29.1547595,
+                     40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 150, 150, 150, 150]
+                c = [0.08, 0.15812003, 0.20314101, 0.23842953, 0.33154153, 0.39131739,
+                     0.46020725, 0.59009569, 0.75306454, 0.81652502, 0.93518422, 0.92030594,
+                     0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96]
+                k = 3
+                def offner_curve(x):
+                    a = -0.16465041
+                    b = -0.11616329
+                    return np.piecewise(x, [x < 6.4, x >= 6.4], [BSpline(t,c,k), lambda x : a * np.exp(b * x) + 0.97])
+                binary_fraction_low = offner_curve(primary_mass[lowmassIdx])
+                binary_choose_low = np.random.uniform(
+                    0, 1.0, primary_mass[lowmassIdx].size)
+
                 (singleIdx_low,) = np.where(
                     binary_fraction_low < binary_choose_low)
                 (binaryIdx_low,) = np.where(
                     binary_fraction_low >= binary_choose_low)
             else:
                 raise ValueError(
-                    "You have supplied a non-supported binary fraction model. Please choose vanHaaften or a float"
+                    "You have supplied a non-supported binary fraction model. Please choose vanHaaften, offner22, or a float"
                 )
         elif type(binfrac_model) == float:
             if (binfrac_model <= 1.0) & (binfrac_model >= 0.0):
@@ -675,7 +701,7 @@ def binary_select(self, primary_mass, binfrac_model=0.5, **kwargs):
                 )
         else:
             raise ValueError(
-                "You have not supplied a model or a fraction. Please choose either vanHaaften or a float"
+                "You have not supplied a model or a fraction. Please choose either vanHaaften, offner22, or a float"
             )
 
         # --- if using a different binary fraction for high-mass systems
@@ -686,13 +712,35 @@ def binary_select(self, primary_mass, binfrac_model=0.5, **kwargs):
                 binary_choose_high = np.random.uniform(
                     0, 1.0, primary_mass[highmassIdx].size)
 
+                (singleIdx_high,) = np.where(
+                    binary_fraction_high < binary_choose_high)
+                (binaryIdx_high,) = np.where(
+                    binary_fraction_high >= binary_choose_high)
+            elif binfrac_model_msort == "offner22":
+                from scipy.interpolate import BSpline
+                t = [0.0331963853, 0.0331963853, 0.0331963853, 0.0331963853, 0.106066017,
+                     0.212132034, 0.424264069, 0.866025404, 1.03077641, 1.11803399,
+                     1.95959179, 3.87298335, 6.32455532, 11.6619038, 29.1547595,
+                     40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 150, 150, 150, 150]
+                c = [0.08, 0.15812003, 0.20314101, 0.23842953, 0.33154153, 0.39131739,
+                     0.46020725, 0.59009569, 0.75306454, 0.81652502, 0.93518422, 0.92030594,
+                     0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96]
+                k = 3
+                def offner_curve(x):
+                    a = -0.16465041
+                    b = -0.11616329
+                    return np.piecewise(x, [x < 6.4, x >= 6.4], [BSpline(t,c,k), lambda x : a * np.exp(b * x) + 0.97])
+                binary_fraction_high = offner_curve(primary_mass[highmassIdx])
+                binary_choose_high = np.random.uniform(
+                    0, 1.0, primary_mass[highmassIdx].size)
+
                 (singleIdx_high,) = np.where(
                     binary_fraction_high < binary_choose_high)
                 (binaryIdx_high,) = np.where(
                     binary_fraction_high >= binary_choose_high)
             else:
                 raise ValueError(
-                    "You have supplied a non-supported binary fraction model. Please choose vanHaaften or a float"
+                    "You have supplied a non-supported binary fraction model. Please choose vanHaaften, offner22, or a float"
                 )
         elif (binfrac_model_msort is not None) and (type(binfrac_model_msort) == float):
             if (binfrac_model_msort <= 1.0) & (binfrac_model_msort >= 0.0):
@@ -711,7 +759,7 @@ def binary_select(self, primary_mass, binfrac_model=0.5, **kwargs):
                 )
         elif (binfrac_model_msort is not None):
             raise ValueError(
-                "You have not supplied a model or a fraction. Please choose either vanHaaften or a float"
+                "You have not supplied a model or a fraction. Please choose either vanHaaften, offner22, or a float"
             )
 
 
@@ -903,7 +951,7 @@ def sample_porb(self, mass1, mass2, rad1, rad2, porb_model, porb_max=None, size=
         elif porb_model == "moe19":
             from scipy.interpolate import interp1d
             from scipy.stats import norm
-            from scipy.integrate import trapz
+            from scipy.integrate import trapezoid
             try:
                 met = kwargs.pop('met')
             except:

diff --git a/src/cosmic/tests/test_sample.py b/src/cosmic/tests/test_sample.py
@@ -42,6 +42,9 @@
 N_BINARY_SELECT = 85
 VANHAAFTEN_BINFRAC_MAX = 0.9989087986493874
 VANHAAFTEN_BINFRAC_MIN = 0.6192803136799157
+OFFNER_MASS_RANGES = [(0.075,0.15), (0.15,0.30), (0.3,0.6), (0.75,1.00), (0.85,1.25), (1.00,1.25), (1.6,2.4), (3,5), (5,8), (8,17), (17,50)]
+OFFNER_DATA = [0.19, 0.23, 0.30, 0.42, 0.47, 0.50, 0.68, 0.81, 0.89, 0.93, 0.96]
+OFFNER_ERRORS = [0.03, 0.02, 0.02, 0.03, 0.03, 0.04, 0.07, 0.06, 0.05, 0.04, 0.04]
 MULTIDIM_BINFRAC_MAX = 0.6146916774140262
 MULTIDIM_BINFRAC_MIN = 0.13786300908773025
 CONST_SFR_SUM = 460028.2453521937
@@ -207,6 +210,20 @@ def test_binary_fraction(self):
         m1_b, m1_s, binfrac, bin_index = SAMPLECLASS.binary_select(primary_mass=np.arange(1,100), binfrac_model='vanHaaften')
         self.assertEqual(binfrac.max(), VANHAAFTEN_BINFRAC_MAX)
         self.assertEqual(binfrac.min(), VANHAAFTEN_BINFRAC_MIN)
+        test_fracs = []
+        test_errs = []
+        primary_mass = np.array([float(x) for x in np.logspace(np.log10(0.08), np.log10(150), num=100000)])
+        m1_b, m1_s, binfrac, bin_index = SAMPLECLASS.binary_select(primary_mass=primary_mass, binfrac_model='offner22')
+        for i in range(len(OFFNER_MASS_RANGES)):
+            low, high = OFFNER_MASS_RANGES[i][0], OFFNER_MASS_RANGES[i][1]
+            offner_value = OFFNER_DATA[i]
+            offner_error = OFFNER_ERRORS[i]
+            bins_count = len(m1_b[(m1_b >= low) & (m1_b <= high)])
+            singles_count
CEB7
 = len(m1_s[(m1_s >= low) & (m1_s <= high)])
+            bin_frac = bins_count / (bins_count + singles_count)
+            error = abs(offner_value - bin_frac)
+            self.assertLess(error, offner_error)
+
 
     def test_msort(self):
         np.random.seed(2)

diff --git a/src/cosmic/utils.py b/src/cosmic/utils.py
@@ -855,7 +855,7 @@ def get_porb_norm(Z, close_logP=4.0, wide_logP=6.0, binfrac_tot_solar=0.66, Z_su
         normalization factor for kde for wide binaries
     '''
     from scipy.stats import norm
-    from scipy.integrate import trapz
+    from scipy.integrate import trapezoid
     from scipy.interpolate import interp1d
 
     # fix to values used in Moe+19
@@ -866,7 +866,7 @@ def get_porb_norm(Z, close_logP=4.0, wide_logP=6.0, binfrac_tot_solar=0.66, Z_su
     logP_pdf = norm.pdf(log_P, loc=4.9, scale=2.3)
 
     # set up the wide binary fraction inflection point
-    norm_wide = binfrac_tot_solar/trapz(logP_pdf, log_P)
+    norm_wide = binfrac_tot_solar/trapezoid(logP_pdf, log_P)
 
     # set up the close binary fraction inflection point
     FeHclose = np.linspace(-3.0, 0.5, 100)
@@ -877,7 +877,7 @@ def get_porb_norm(Z, close_logP=4.0, wide_logP=6.0, binfrac_tot_solar=0.66, Z_su
     fclose_interp = interp1d(Zclose, fclose)
 
     fclose_Z = fclose_interp(Z)
-    norm_close = fclose_Z/trapz(logP_pdf[log_P < close_logP], log_P[log_P < close_logP])
+    norm_close = fclose_Z/trapezoid(logP_pdf[log_P < close_logP], log_P[log_P < close_logP])
 
     return norm_wide, norm_close
 
@@ -904,7 +904,7 @@ def get_met_dep_binfrac(met):
     neval = 5000
 
     from scipy.interpolate import interp1d
-    from scipy.integrate import trapz
+    from scipy.integrate import trapezoid
     from scipy.stats import norm
 
     norm_wide, norm_close = get_porb_norm(met)
@@ -919,7 +919,7 @@ def get_met_dep_binfrac(met):
                        np.linspace(wide_logP, logP_hi_lim, neval),])
     y_dat = np.hstack([prob_close, prob_interp_int(np.linspace(close_logP, wide_logP, neval)), prob_wide])
 
-    binfrac = trapz(y_dat, x_dat)/0.66 * 0.5
+    binfrac = trapezoid(y_dat, x_dat)/0.66 * 0.5
 
     return float(np.round(binfrac, 2))