diff --git a/.github/workflows/build_wheels_and_publish.yml b/.github/workflows/build_wheels_and_publish.yml index c4efa6b99..7813d006d 100644 --- 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 index b1b843a9d..f3bd3af3b 100755 --- 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 index 6e46c2373..561e029b9 100644 --- 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` \ No newline at end of file + - 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 index 4e7f210cc..c169885e4 100644 --- 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 index b0b273599..883697170 100644 --- 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`` diff --git a/docs/pages/sample/independent.rst b/docs/pages/sample/independent.rst index dd354bd10..0d5994289 100644 --- a/docs/pages/sample/independent.rst +++ b/docs/pages/sample/independent.rst @@ -47,9 +47,9 @@ If you don't want to filter the binaries, you can supply final kstars as In [6]: final_kstars = np.linspace(0, 14, 15) - In [7]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = 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=10000) + In [7]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstars, final_kstars, binfrac_model=0.5, primary_model='kroupa01', ecc_model='sana12', porb_model='sana12', qmin=-1, SF_start=13700.0, SF_duration=0.0, met=0.02, size=10000) -Additionally if you are interested in single stars then you can specify ``keep_singles=True``. +Additionally if you are interested in single stars then you can specify ``keep_singles=True``. In this case, the singles will be added onto the end of the InitialBinaryTable where ``kstar_1`` will host the singles, ``kstar_2`` will be filled with 15s only, and all orbital properties (e.g. ``porb`` or ``ecc``) will be indicated with -1. Understanding parameter sampling models ======================================= @@ -76,7 +76,7 @@ Using the final kstar inputs we mentioned above, the initial binary population c .. ipython:: - In [9]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstar1, final_kstar2, 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=10000) + In [9]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstar1, final_kstar2, binfrac_model=0.5, primary_model='kroupa01', ecc_model='sana12', porb_model='sana12', qmin=-1, SF_start=13700.0, SF_duration=0.0, met=0.02, size=10000) In [10]: print(InitialBinaries) @@ -95,7 +95,7 @@ Alternatively, we could do the same thing but now instead set our ``sampling_tar .. ipython:: - In [10]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstar1, final_kstar2, 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, sampling_target="total_mass", total_mass=15000) + In [10]: InitialBinaries, mass_singles, mass_binaries, n_singles, n_binaries = InitialBinaryTable.sampler('independent', final_kstar1, final_kstar2, binfrac_model=0.5, primary_model='kroupa01', ecc_model='sana12', porb_model='sana12', qmin=-1, SF_start=13700.0, SF_duration=0.0, met=0.02, sampling_target="total_mass", total_mass=15000) In [11]: print(InitialBinaries) diff --git a/docs/pages/sample/multidim.rst b/docs/pages/sample/multidim.rst index b1cdebe98..dd53e993a 100644 --- a/docs/pages/sample/multidim.rst +++ b/docs/pages/sample/multidim.rst @@ -33,3 +33,7 @@ The multidimensional sample is generated as follows: NOTE that in the multidimensional case, the binary fraction is a parameter in the sample. This results in the size of the initial binary data matching the size provided to the sampler. As in the independent sampling case, we keep track of the total sampled mass of singles and binaries as well as the total number of single and binary stars to scale the simulated population to astrophysical populations. +.. note:: + + NOTE that you can also keep singles for the multidim sampelr as well. As with the independent sampler, the singles will be added onto the end of the InitialBinaryTable where ``kstar_1`` will host the singles, ``kstar_2`` will be filled with 15s only, and all orbital properties (e.g. ``porb`` or ``ecc``) will be indicated with -1. + diff --git a/meson.build b/meson.build index 7e0a669c9..a9c6c98c4 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,7 @@ project('cosmic', 'c', 'fortran', - version : '3.5.0', + version : '3.6.1', default_options: ['warning_level=0', 'optimization=3'], ) diff --git a/pyproject.toml b/pyproject.toml index 9be1bf3df..7b747791b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,7 +10,7 @@ authors = [ { name="Tom Wagg" }, ] readme = "README.md" -version = "3.6.0" +version = "3.6.1" requires-python = ">=3.9" license = { text = "MIT License" } classifiers = [ diff --git a/src/cosmic/evolve.py b/src/cosmic/evolve.py index 0d46e0254..f09efcf40 100644 --- a/src/cosmic/evolve.py +++ b/src/cosmic/evolve.py @@ -46,7 +46,7 @@ __all__ = ['Evolve'] - +# Make this match the ordering of all_cols in bpp_array.f ALL_COLUMNS = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2', 'sep', 'porb', 'ecc', 'RRLO_1', 'RRLO_2', 'evol_type', 'aj_1', 'aj_2', 'tms_1', 'tms_2', 'massc_1', 'massc_2', 'rad_1', 'rad_2', 'mass0_1', diff --git a/src/cosmic/sample/initialcmctable.py b/src/cosmic/sample/initialcmctable.py index c346ed4ab..a17042826 100644 --- a/src/cosmic/sample/initialcmctable.py +++ b/src/cosmic/sample/initialcmctable.py @@ -332,9 +332,10 @@ def write(cls, Singles, Binaries, filename="input.hdf5", **kwargs): ) singles = pd.concat([singles, Singles]) singles = pd.concat([singles, singles_bottom]) - singles["r"].iloc[-1] = 1e40 - singles["r"].iloc[0] = 2.2250738585072014e-308 - singles["m"].iloc[0] = Singles.central_bh + + singles.iloc[-1, singles.columns.get_loc("r")] = 1e40 + singles.iloc[0, singles.columns.get_loc("r")] = 2.2250738585072014e-308 + singles.iloc[0, singles.columns.get_loc("m")] = Singles.central_bh # Add a special row to the end of Bianries table binaries = pd.DataFrame( diff --git a/src/cosmic/sample/sampler/independent.py b/src/cosmic/sample/sampler/independent.py index 010401c4f..fb6a052c5 100644 --- a/src/cosmic/sample/sampler/independent.py +++ b/src/cosmic/sample/sampler/independent.py @@ -362,8 +362,7 @@ def get_independent_sampler( m_sampled_singles, m_sampled_binaries, n_singles, - n_binaries, - len(mass1_singles) + n_binaries ) @@ -506,7 +505,6 @@ def sample_secondary(self, primary_mass, q_power_law=0, **kwargs): sampled secondary masses with array size matching size of primary_mass """ - qmin = kwargs["qmin"] if "qmin" in kwargs.keys() else 0.0 m1_min = kwargs["m1_min"] if "m1_min" in kwargs.keys() else 0.08 m2_min = kwargs["m2_min"] if "m2_min" in kwargs.keys() else None @@ -516,6 +514,9 @@ def sample_secondary(self, primary_mass, q_power_law=0, **kwargs): raise ValueError("The m2_min you specified is above the minimum" " primary mass of the IMF, either lower m2_min or" " raise the lower value of your sampled primaries") + + if (m2_min is not None) & (qmin != 0): + raise ValueError("You cannot specify both m2_min and qmin, please choose one or the other") # --- `msort` kwarg can be set to have different qmin above `msort` msort = kwargs["msort"] if "msort" in kwargs.keys() else None diff --git a/src/cosmic/sample/sampler/multidim.py b/src/cosmic/sample/sampler/multidim.py index f658b27a2..219bcbae1 100644 --- a/src/cosmic/sample/sampler/multidim.py +++ b/src/cosmic/sample/sampler/multidim.py @@ -176,7 +176,7 @@ def get_multidim_sampler( metallicity[metallicity < 1e-4] = 1e-4 metallicity[metallicity > 0.03] = 0.03 - if kwargs.pop("keep_singles", False): + if kwargs.pop("keep_singles", True): binary_table = InitialBinaryTable.InitialBinaries( mass1_binary, mass2_binary, @@ -201,7 +201,7 @@ def get_multidim_sampler( np.ones_like(single_mass_list)*-1, tphysf, kstar1, - np.ones_like(single_mass_list)*0, + np.ones_like(single_mass_list)*15, # # kstar2 is not used for singles metallicity, ) binary_table = pd.concat([binary_table, singles_table], ignore_index=True) diff --git a/src/cosmic/tests/test_sample.py b/src/cosmic/tests/test_sample.py index 2f642bd9d..5b6ac9321 100644 --- a/src/cosmic/tests/test_sample.py +++ b/src/cosmic/tests/test_sample.py @@ -186,11 +186,12 @@ def test_sample_secondary(self): slope = linear_fit(q) self.assertEqual(np.round(slope, 1), FLAT_SLOPE) - np.random.seed(2) - mass2 = SAMPLECLASS.sample_secondary(primary_mass=mass1, qmin=0.1, m2_min=0.08) - q = mass2[ind_massive] / mass1[ind_massive] - slope = linear_fit(q) - self.assertEqual(np.round(slope, 1), FLAT_SLOPE) + # This is now redundant since you should only sample with qmin or m2_min + #np.random.seed(2) + #mass2 = SAMPLECLASS.sample_secondary(primary_mass=mass1, qmin=0.1, m2_min=0.08) + #q = mass2[ind_massive] / mass1[ind_massive] + #slope = linear_fit(q) + #self.assertEqual(np.round(slope, 1), FLAT_SLOPE) def test_sample_q(self): """Test you can sample different mass ratio distributions""" @@ -239,15 +240,15 @@ def test_msort(self): np.random.seed(2) mass1, total_mass = SAMPLECLASS.sample_primary(primary_model='kroupa01', size=1000000) # Check that qmin_msort and m2_min_msort are workings as expected - mass2 = SAMPLECLASS.sample_secondary(primary_mass = mass1, qmin=0.1, m2_min=0.08, msort=15, qmin_msort=0.7, m2_min_msort=12) + mass2 = SAMPLECLASS.sample_secondary(primary_mass = mass1, qmin=0.1, msort=15, qmin_msort=0.7, m2_min_msort=12) ind_light, = np.where(mass1 < 15.0) ind_massive, = np.where(mass1 >= 15.0) m2_light = mass2[ind_light] m2_massive = mass2[ind_massive] q_light = mass2[ind_light]/mass1[ind_light] q_massive = mass2[ind_massive]/mass1[ind_massive] - assert m2_light.min() > M2MIN_LOWMASS - assert m2_massive.min() > M2MIN_HIGHMASS + assert m2_light.min() > np.min(mass1[ind_light]) * 0.1 + assert m2_massive.min() > np.min(mass1[ind_massive]) * 0.7 assert q_light.min() > QMIN_LOWMASS assert q_massive.min() > QMIN_HIGHMASS # Check that the binary fraction tracking is correct when using msort @@ -324,7 +325,6 @@ def test_sample_porb(self): metallicity = 0.001 # this is a metallicity dependent population: binfrac = get_met_dep_binfrac(metallicity) - print(binfrac) mass1, total_mass = SAMPLECLASS.sample_primary(primary_model='kroupa01', size=100000) (mass1_binaries, mass_single, binfrac_binaries, binary_index,) = SAMPLECLASS.binary_select( mass1, binfrac_model=binfrac, diff --git a/src/cosmic/utils.py b/src/cosmic/utils.py index 124a1e7e7..2d97ca898 100644 --- a/src/cosmic/utils.py +++ b/src/cosmic/utils.py @@ -415,6 +415,12 @@ def pop_write( bpp_singles : `pandas.DataFrame` kwargs bpp_singles array to write + initC_singles : `pandas.DataFrame` + kwargs initC_singles array to write + + kick_info_singles : `pandas.DataFrame` + kwargs kick_info_singles array to write + Returns ------- Nothing! @@ -436,7 +442,7 @@ def pop_write( # Save the initial binaries # ensure that the index corresponds to bin_num - dat_store.append("initCond", initC.set_index("bin_num", drop=False)) + dat_store.append("initC", initC.set_index("bin_num", drop=False)) # Save the converging dataframe dat_store.append("conv", conv) @@ -455,15 +461,21 @@ def pop_write( if "conv_singles" in kwargs.keys(): - # Save the singles dataframe + # Save the singles conv dataframe dat_store.append("conv_singles", kwargs["conv_singles"]) - # Save the singles dataframe + # Save the singles bcm dataframe dat_store.append("bcm_singles", kwargs["bcm_singles"]) - # Save the singles dataframe + # Save the singles bpp dataframe dat_store.append("bpp_singles", kwargs["bpp_singles"]) + # save the singles initCond dataframe + dat_store.append("initC_singles", kwargs["initC_singles"]) + + # save the singles kick_info dataframe + dat_store.append("kick_info_singles", kwargs["kick_info_singles"]) + return @@ -577,7 +589,7 @@ def mass_min_max_select(kstar_1, kstar_2, **kwargs): if ((primary_min < 0.08) | (secondary_min < 0.08)): warnings.warn("Tread carefully, BSE is not equipped to handle stellar masses less than 0.08 Msun!") if primary_max > 150: - warnings.warn("Tread carefully, BSE is not equipped to handle stellar masses greater than 150 Msun!") + warnings.warn("Tread carefully, BSE is not equipped to handle stellar masses greater than 150 Msun! And to be honest, we are extrapolating beyond 50 Msun :-/") min_mass = [primary_min, secondary_min] max_mass = [primary_max, secondary_max]