COSMIC-PopSynth · katiebreivik · Jul 10, 2025 · Jun 4, 2025 · Jun 4, 2025 · Jun 4, 2025
diff --git a/.github/workflows/build_wheels_and_publish.yml b/.github/workflows/build_wheels_and_publish.yml
@@ -6,7 +6,7 @@ on:
   workflow_dispatch:
 
 env:
-  CIBW_BUILD: "cp39-* cp310-*"
+  CIBW_BUILD: "cp39-* cp310-* cp311-* cp312-*"
   CIBW_ARCHS_LINUX: "x86_64"
   CIBW_SKIP: "*-win32 *musllinux*"
   CIBW_MANYLINUX_X86_64_IMAGE: manylinux2014
@@ -20,7 +20,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [3.9, "3.10"]
+        python-version: [3.9, "3.10", "3.11", "3.12"]
 
     steps:
     - uses: actions/checkout@v4

diff --git a/bin/cosmic-pop b/bin/cosmic-pop
@@ -11,29 +11,24 @@
 ##############################################################################
 #  IMPORT ALL NECESSARY PYTHON PACKAGES
 ##############################################################################
-from collections import OrderedDict
-import warnings
 import argparse
 import schwimmbad
 
-import math
-import random
 import time
 from time import sleep
-import string
-import os.path
 import json
 
 import numpy as np
-import scipy.special as ss
 import pandas as pd
+from pandas.errors import PerformanceWarning
 import warnings
 
 from cosmic.sample.initialbinarytable import InitialBinaryTable
 from cosmic import Match, utils
 from cosmic.evolve import Evolve
 
-from schwimmbad import MultiPool, MPIPool
+from schwimmbad import MPIPool
+from os import sys
 
 def str2bool(v):
     if isinstance(v, bool):
@@ -243,10 +238,12 @@ if __name__ == '__main__':
         kstar2_range = args.final_kstar2
         kstar2_range_string = str(int(args.final_kstar2[0]))
 
+    dat_store_fname = 'dat_kstar1_{0}_kstar2_{1}_SFstart_{2}_SFduration_{3}_metallicity_{4}.h5'.format(kstar1_range_string, kstar2_range_string, sampling['SF_start'], sampling['SF_duration'], sampling['metallicity'])
     # Open the hdf5 file to store the fixed population data
     try:
-        dat_store = pd.HDFStore('dat_kstar1_{0}_kstar2_{1}_SFstart_{2}_SFduration_{3}_metallicity_{4}.h5'.format(kstar1_range_string, kstar2_range_string, sampling['SF_start'], sampling['SF_duration'], sampling['metallicity']),complib=args.complib,complevel=args.complevel)
-        conv_save = pd.read_hdf(dat_store, 'conv')
+        with pd.HDFStore(dat_store_fname,complib=args.complib,complevel=args.complevel) as dat_store:
+            # If the file exists, we will read it and continue from where we left off
+            conv_save = pd.read_hdf(dat_store, 'conv')
         log_file = open('log_kstar1_{0}_kstar2_{1}_SFstart_{2}_SFduration_{3}_metallicity_{4}.txt'.format(kstar1_range_string, kstar2_range_string, sampling['SF_start'], sampling['SF_duration'], sampling['metallicity']), 'a')
         log_file.write('There are already: '+str(conv_save.shape[0])+' '+kstar1_range_string+'_'+kstar2_range_string+' binaries evolved\n')
         log_file.write('\n')
@@ -258,8 +255,8 @@ if __name__ == '__main__':
         total_n_stars = np.max(pd.read_hdf(dat_store, 'n_stars'))[0]
         idx = int(np.max(pd.read_hdf(dat_store, 'idx'))[0])
     except:
+        #dat_store = pd.HDFStore('dat_kstar1_{0}_kstar2_{1}_SFstart_{2}_SFduration_{3}_metallicity_{4}.h5'.format(kstar1_range_string, kstar2_range_string, sampling['SF_start'], sampling['SF_duration'], sampling['metallicity']),complib=args.complib,complevel=args.complevel)
         conv_save = pd.DataFrame()
-        dat_store = pd.HDFStore('dat_kstar1_{0}_kstar2_{1}_SFstart_{2}_SFduration_{3}_metallicity_{4}.h5'.format(kstar1_range_string, kstar2_range_string, sampling['SF_start'], sampling['SF_duration'], sampling['metallicity']),complib=args.complib,complevel=args.complevel)
         total_mass_singles = 0  
         total_mass_binaries = 0
         total_mass_stars = 0
@@ -273,10 +270,14 @@ if __name__ == '__main__':
     configuration_settings = {'BSEDict' : BSEDict, 'filters' : filters,
                               'convergence' : convergence, 'sampling' : sampling}
 
-    for k, v in configuration_settings.items():
-        for k1, v1 in v.items():
-            dat_store.put('config/{0}/{1}/'.format(k, k1), pd.Series(v1))
-    dat_store.put('config/rand_seed/', pd.Series(seed_int))
+    with warnings.catch_warnings():
+        warnings.simplefilter(action="ignore", category=PerformanceWarning)
+
+        with pd.HDFStore(dat_store_fname,complib=args.complib,complevel=args.complevel) as dat_store:
+            for k, v in configuration_settings.items():
+                for k1, v1 in v.items():
+                    dat_store.put('config/{0}/{1}/'.format(k, k1), pd.Series(v1))
+            dat_store.put('config/rand_seed/', pd.Series(seed_int))
 
     # Initialize the step counter and convergence array/list
     Nstep = idx - np.mod(idx, args.Nstep)
@@ -304,6 +305,9 @@ if __name__ == '__main__':
 
         # Select the initial binary sample method from user input
         if sampling['sampling_method'] == 'independent':
+            if hasattr(args,'qmin') and hasattr(args,'m2_min'):
+                raise ValueError("You cannot specify both qmin and m2_min in the inifile if you are using the independent sampler. Please choose one or the other.")
+            # If qmin is specified, use it to sample the initial binary table
             if hasattr(args,'qmin'):
                 init_samp_list = InitialBinaryTable.sampler(format_ = sampling['sampling_method'],
                                                             final_kstar1 = kstar1_range,
@@ -319,6 +323,7 @@ if __name__ == '__main__':
                                                             size = args.Nstep,
                                                             qmin = args.qmin,
                                                             params = args.inifile)
+            # if m2_min is specified, use it to sample the initial binary table
             elif hasattr(args,'m2_min'):
                 init_samp_list = InitialBinaryTable.sampler(format_ = sampling['sampling_method'],
                                                             final_kstar1 = kstar1_range,
@@ -337,20 +342,21 @@ if __name__ == '__main__':
             else:
                 raise ValueError("You must specify either qmin or m2_min in the",
                                  " inifile if you are using the independent sampler")
-            IBT, mass_singles, mass_binaries, n_singles, n_binaries, n_singles_table = init_samp_list
+            IBT, mass_singles, mass_binaries, n_singles, n_binaries = init_samp_list
 
         if sampling['sampling_method'] == 'multidim':
             init_samp_list = InitialBinaryTable.sampler(format_ = sampling['sampling_method'],
                                                         final_kstar1 = kstar1_range,
                                                         final_kstar2 = kstar2_range,
-                                                        keep_singles = args.keep_singles,
                                                         rand_seed = rand_seed,
                                                         nproc = args.nproc,
                                                         SF_start = sampling['SF_start'],
                                                         SF_duration = sampling['SF_duration'],
                                                         met = sampling['metallicity'],
                                                         size = args.Nstep,
-                                                        pool=pool)
+                                                        pool=pool,
+                                                        keep_singles = args.keep_singles
+                                                        )
             IBT, mass_singles, mass_binaries, n_singles, n_binaries = init_samp_list
 
         # Log the total sampled mass from the initial binary sample
@@ -390,30 +396,30 @@ if __name__ == '__main__':
         # extract single stars
 
         if (args.keep_singles==True):
-            bcm_init = bcm.loc[bcm.tphys==0]
-            if sampling['sampling_method'] == 'multidim':
-                bcm_init_singles = bcm_init.iloc[-n_singles:]
-            if sampling['sampling_method'] == 'independent':
-                bcm_init_singles = bcm_init.iloc[-n_singles_table:]
-
-            bcm_singles = bcm.loc[bcm.bin_num.isin(bcm_init_singles.bin_num)]
-            bpp_singles = bpp.loc[bpp.bin_num.isin(bcm_init_singles.bin_num)]
-
-            bpp = bpp.loc[~bpp.bin_num.isin(bcm_init_singles.bin_num)]
-            bcm = bcm.loc[~bcm.bin_num.isin(bcm_init_singles.bin_num)]
-            initCond = initCond.loc[~initCond.bin_num.isin(bcm_init_singles.bin_num)]
-            kick_info = kick_info.loc[~kick_info.bin_num.isin(bcm_init_singles.bin_num)]
+            singles_bin_num = initCond.loc[initCond.kstar_2 == 15].bin_num.unique()
+            # get the singles from the bcm and bpp arrays
+            bcm_singles = bcm.loc[bcm.bin_num.isin(singles_bin_num)]
+            bpp_singles = bpp.loc[bpp.bin_num.isin(singles_bin_num)]
+            initCond_singles = initCond.loc[initCond.bin_num.isin(singles_bin_num)]
+            kick_info_singles = kick_info.loc[kick_info.bin_num.isin(singles_bin_num)]
+
+            bpp = bpp.loc[~bpp.bin_num.isin(singles_bin_num)]
+            bcm = bcm.loc[~bcm.bin_num.isin(singles_bin_num)]
+            initCond = initCond.loc[~initCond.bin_num.isin(singles_bin_num)]
+            kick_info = kick_info.loc[~kick_info.bin_num.isin(singles_bin_num)]
         # get any nans and pull them out for now
         nans = np.isnan(bpp.sep)
         if nans.any():
             nan_bin_nums = np.unique(bpp[nans]["bin_num"].values)
             initCond_nan = initCond.loc[initCond.bin_num.isin(nan_bin_nums)]
-            if pd.__version__<="2.0.0":
-                dat_store.append("nan_initC", initCond_nan)
-            else:
-                dat_store["nan_initC"] = initCond_nan
-            log_file.write(f"There are {len(nan_bin_nums)} NaNs stored in the datfile with key: 'nan_initC'")
-            log_file.write(f"These NaNs likely arise because you have pts1 = 0.001, try running with pts1 = 0.01")
+            with pd.HDFStore(dat_store_fname,complib=args.complib,complevel=args.complevel) as dat_store:
+                if pd.__version__<="2.0.0":
+                    dat_store.append("nan_initC", initCond_nan)
+                else:
+                    dat_store["nan_initC"] = initCond_nan
+            log_file.write(f"There are {len(nan_bin_nums)} NaNs stored in the datfile with key: 'nan_initC'\n")
+            log_file.write(f"You might want to check them out carefully to see if there is something that impacts your results\n")
+            #log_file.write(f"These NaNs likely arise because you have pts1 = 0.001, try running with pts1 = 0.01")
 
             bcm = bcm.loc[~bcm.bin_num.isin(nan_bin_nums)]
             bpp = bpp.loc[~bpp.bin_num.isin(nan_bin_nums)]
@@ -444,11 +450,14 @@ if __name__ == '__main__':
 
         bcm_filter = bcm.loc[bcm.bin_num.isin(conv_filter.bin_num)]
         bpp_filter = bpp.loc[bpp.bin_num.isin(conv_filter.bin_num)]
+        initC_filter = initCond.loc[initCond.bin_num.isin(conv_filter.bin_num)]
+        kick_info_filter = kick_info.loc[kick_info.bin_num.isin(conv_filter.bin_num)]
+
         if (args.keep_singles==True):
             bpp_singles_filter = bpp_singles.loc[bpp_singles.bin_num.isin(conv_singles_filter.bin_num)]
             bcm_singles_filter = bcm_singles.loc[bcm_singles.bin_num.isin(conv_singles_filter.bin_num)]
-        initC_filter = initCond.loc[initCond.bin_num.isin(conv_filter.bin_num)]
-        kick_info_filter = kick_info.loc[kick_info.bin_num.isin(conv_filter.bin_num)]
+            initC_singles_filter = initCond_singles.loc[initCond_singles.bin_num.isin(conv_singles_filter.bin_num)]
+            kick_info_singles_filter = kick_info_singles.loc[kick_info_singles.bin_num.isin(conv_singles_filter.bin_num)]
 
         # Filter the bin_state based on user specified filters
         bcm_filter, bin_state_nums = utils.filter_bin_state(bcm_filter, bpp_filter, filters, kstar1_range, kstar2_range)
@@ -476,7 +485,9 @@ if __name__ == '__main__':
         if (args.keep_singles==True):
             conv_singles_filter_match = conv_singles_filter.copy()
             bpp_singles_filter_match = bpp_singles_filter.copy()
-            bcm_singles_filter_match = bcm_singles_filter.copy()          
+            bcm_singles_filter_match = bcm_singles_filter.copy()
+            initC_filter_singles_match = initC_singles_filter.copy()
+            kick_info_singles_filter_match = kick_info_singles_filter.copy()          
 
         if len(conv_filter_match) >= np.min([50, args.Niter]):
             conv_save = pd.concat([conv_save, pd.DataFrame(conv_filter_match)], ignore_index=True)
@@ -499,14 +510,19 @@ if __name__ == '__main__':
             mass_list = [total_mass_singles, total_mass_binaries, total_mass_stars]
             n_list = [total_n_singles, total_n_binaries, total_n_stars]
 
-            if (args.keep_singles==True):
-                utils.pop_write(dat_store, log_file, mass_list, n_list, bcm_filter_match,
-                            bpp_filter_match, initC_filter_match, conv_filter_match, kick_info_filter_match,
-                            bin_state_nums, match_save, idx, conv_singles=conv_singles_filter_match, bcm_singles=bcm_singles_filter_match, bpp_singles=bpp_singles_filter_match)
-            else:
-                utils.pop_write(dat_store, log_file, mass_list, n_list, bcm_filter_match,
-                            bpp_filter_match, initC_filter_match, conv_filter_match, kick_info_filter_match,
-                            bin_state_nums, match_save, idx)
+            # write the data to the dat_store
+            with pd.HDFStore(dat_store_fname,complib=args.complib,complevel=args.complevel) as dat_store:
+                if (args.keep_singles==True):
+                    utils.pop_write(dat_store, log_file, mass_list, n_list, bcm_filter_match,
+                                bpp_filter_match, initC_filter_match, conv_filter_match, kick_info_filter_match,
+                                bin_state_nums, match_save, idx, 
+                                conv_singles=conv_singles_filter_match, bcm_singles=bcm_singles_filter_match, 
+                                bpp_singles=bpp_singles_filter_match, initC_singles=initC_filter_singles_match,
+                                kick_info_singles=kick_info_singles_filter_match)
+                else:
+                    utils.pop_write(dat_store, log_file, mass_list, n_list, bcm_filter_match,
+                                bpp_filter_match, initC_filter_match, conv_filter_match, kick_info_filter_match,
+                                bin_state_nums, match_save, idx)
 
             # reset the bcm_filter DataFrame
             bcm_filter_match = []
@@ -518,12 +534,16 @@ if __name__ == '__main__':
                 conv_singles_filter_match = []
                 bpp_singles_filter_match = []
                 bcm_singles_filter_match = []
+                initC_filter_singles_match = []
+                kick_info_singles_filter_match = []
             log_file.write('\n')
         Nstep += args.Nstep
         log_file.flush()
-    # Close the data storage file
-    dat_store.close()
-
+
+    # close the log file and print the final message
     log_file.write('All done friend!')
     log_file.close()
 
+    pool.close()
+    pool.join()
+
diff --git a/changelog.md b/changelog.md
@@ -50,6 +50,8 @@ See the discussed changes in our previous releases here: https://github.com/COSM
 
 ## 3.6.0
  - Overhaul documentation and added debugging environment
+ - Feature: Added Disberg+2025 kick prescription as a new choice of `kickflag` (`kickflag=5`). Applies log-normal distribution to regular CCSN, ECSN/USSN still use `sigmadiv` Maxwellian and BH fallback scaling is still applied via `bhflag` and `bhsigmafrac` as with `kickflag=1`
 
 ## 3.6.1
- - Feature: Added Disberg+2025 kick prescription as a new choice of `kickflag` (`kickflag=5`). Applies log-normal distribution to regular CCSN, ECSN/USSN still use `sigmadiv` Maxwellian and BH fallback scaling is still applied via `bhflag` and `bhsigmafrac` as with `kickflag=1`
+ - Add support for single stars in both independent and multidim sampling
+ - update documentation
diff --git a/docs/pages/evolve/evolve_sample.rst b/docs/pages/evolve/evolve_sample.rst
@@ -27,8 +27,8 @@ about the independent sampler in the :ref:`independent` page.
        ...:       InitialBinaryTable.sampler('independent', final_kstars, final_kstars,
        ...:                                  binfrac_model=0.5, primary_model='kroupa01',
        ...:                                  ecc_model='sana12', porb_model='sana12',
-       ...:                                  qmin=-1, m2_min=0.08, SF_start=13700.0,
-       ...:                                  SF_duration=0.0, met=0.02, size=10)
+       ...:                                  qmin=-1, SF_start=13700.0,
+       ...:                                  SF_duration=0.0, met=0.02, size=10, keep_singles=True)
 
 
 And finally, we can evolve the initial binary population using the Evolve class as we've done in the previous

diff --git a/docs/pages/fixedpop.rst b/docs/pages/fixedpop.rst
@@ -92,7 +92,15 @@ The fixed population contains several pandas DataFrames accessed by the followin
 
 * ``kick_info`` : The magnitude and direction of natal kicks, three dimensional systemic velocity changes, total tilt of orbital plane, and azimuthal angle of orbital angular momentum axis with respect to spins
 
-* ``initCond`` : The initial conditions for each binary which satisfies the user-specified final kstars and filter in the ``convergence`` subsection
+* ``initC`` : The initial conditions for each binary which satisfies the user-specified final kstars and filter in the ``convergence`` subsection
+
+* ``bpp_singles`` : The evolutionary history of single stars which satisfy the user-specified final kstars and filter in the ``convergence`` subsection
+
+* ``bcm_singles`` : The final state of single stars in the bcm array which satisfy the user-specified final kstars and filter in the ``convergence`` subsection
+
+* ``kick_info_singles`` : The magnitude and direction of natal kicks, three dimensional systemic velocity changes, total tilt of orbital plane, and azimuthal angle of orbital angular momentum axis with respect to spins
+
+* ``initC_singles`` : The initial conditions for each single star which satisfies the user-specified final kstars and filter in the ``convergence`` subsection
 
 * ``idx`` : An integer that keeps track of the total number of simulated binaries to maintain proper indexing across several runs of ``cosmic-pop``