Tutorial: Running populations with binary_c-python

This notebook will show you how to evolve a population of stars

Much of the code in the binarycpython package is written to evolve a population of stars through the Population object, rather than running a single system. Let’s go through the functionality of this object step by step and set up some example populations.

At the bottom of this notebook there are some complete example scripts

[1]:
import os

from binarycpython.utils.custom_logging_functions import temp_dir
from binarycpython import Population

TMP_DIR = temp_dir("notebooks", "notebook_population")

# help(Population) # Uncomment to see the public functions of this object

Setting up the Population object

To set up and configure the population object we need to make an object instance of the Population object, and add configuration via the .set() function.

There are three categories of options that the Population object can set: - BSE options: these options will be used for the binary_c calls, and are recognized by comparing the arguments to a known list of available arguments of binary_c. To see which options are available, see section `binary_c parameters in the documentation <https://ri0005.pages.surrey.ac.uk/binary_c-python/binary_c_parameters.html>`__. You can access these through population.bse_options['<bse option name>'] after you have set them.

  • Grid options: these options will be used to configure the behaviour of the Population object. To see which options are available, see section `Population grid code options in the documentation <https://ri0005.pages.surrey.ac.uk/binary_c-python/grid_options_descriptions.html>`__. They can be accessed via population.grid_options['<grid option name>'] after you have set them.

  • Custom options: these options are not recognized as either of the above, so they will be stored in the custom_options, and can be accessed via population.custom_options['<custom option name>']

[2]:
# Create population object
example_pop = Population()

# If you want verbosity, set this before other things
example_pop.set(verbosity=1)

# Setting values can be done via .set(<parameter_name>=<value>)
# Values that are known to be binary_c_parameters are loaded into bse_options.
# Those that are present in the default grid_options are set in grid_options
# All other values that you set are put in a custom_options dict
example_pop.set(
    # binary_c physics options
    M_1=10,  # bse_options
    orbital_period=45000000080,  # bse_options
    max_evolution_time=15000,  # bse_options
    eccentricity=0.02,  # bse_options


    # grid_options
    num_cores=2,
    tmp_dir=TMP_DIR,

    # Custom options # TODO: need to be set in grid_options probably
    data_dir=os.path.join(
        TMP_DIR, "example_python_population_result"
    ),  # custom_options
    base_filename="example_pop.dat",  # custom_options
)

# We can use the options through
print(example_pop.grid_options['verbosity'])
print(example_pop.custom_options['base_filename'])
print(example_pop.bse_options['M_1'])
adding: num_cores=2 to grid_options
adding: tmp_dir=/tmp/binary_c_python-david/notebooks/notebook_population to grid_options
<<<< Warning: Key does not match previously known parameter:                     adding: data_dir=/tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result to custom_options >>>>
<<<< Warning: Key does not match previously known parameter:                     adding: base_filename=example_pop.dat to custom_options >>>>
1
example_pop.dat
10

After configuring the population, but before running the actual population, its usually a good idea to export the full configuration (including version info of binary_c and all the parameters) to a file. To do this we use example_pop.export_all_info().

On default this exports everything, each of the sections can be disabled:

  • population settings (bse_options, grid_options, custom_options), turn off with include_population settings=False

  • binary_c_defaults (all the commandline arguments that binary c accepts, and their defaults). turn off with include_binary_c_defaults=False

  • include_binary_c_version_info (all the compilation info, and information about the compiled parameters), turn off with include_binary_c_version_info=False

  • include_binary_c_help_all (all the help information for all the binary_c parameters), turn off with include_binary_c_help_all=Fase

On default it will write this to the custom_options[‘data_dir’], but that can be overriden by setting use_datadir=False and providing an outfile=<>

[3]:
example_pop.export_all_info()
ok
File at /tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/example_pop_settings.json.gz already exists: cannot write to it
ok
ok pre
ok ret
[3]:
'/tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/example_pop_settings.json.gz'

Adding grid variables

The main purpose of the Population object is to handle the population synthesis side of running a set of stars. The main method to do this with binarycpython, as is the case with Perl binarygrid, is to use grid variables. These are loops over a predefined range of values, where a probability will be assigned to the systems based on the chosen probability distributions.

Usually we use either 1 mass grid variable, or a trio of mass, mass ratio and period (See below for full examples of all of these). We can, however, also add grid sampling for e.g. eccentricity, metallicity or other parameters.

In some cases it could be easier to set up a for loop that sets that parameter and calls the evolve function several times, e.g. when you want to vary a prescription (usually a discrete, unweighted parameter)

A notable special type of grid variable is that of the Moe & di Stefano 2017 dataset (see further down in the notebook).

To add a grid variable to the population object we use example_pop.add_grid_variable (see next cell)

[4]:
help(example_pop.add_grid_variable)
Help on method add_grid_variable in module binarycpython.utils.population_extensions.gridcode:

add_grid_variable(name: str, parameter_name: str, longname: str, valuerange: Union[list, str], samplerfunc: str, probdist: str, dphasevol: Union[str, int] = -1, gridtype: str = 'centred', branchpoint: int = 0, branchcode: Optional[str] = None, precode: Optional[str] = None, postcode: Optional[str] = None, topcode: Optional[str] = None, bottomcode: Optional[str] = None, condition: Optional[str] = None, index: Optional[int] = None, dry_parallel: Optional[bool] = False) -> None method of binarycpython.utils.grid.Population instance
    Function to add grid variables to the grid_options.

    The execution of the grid generation will be through a nested for loop.
    Each of the grid variables will get create a deeper for loop.

    The real function that generates the numbers will get written to a new file in the TMP_DIR,
    and then loaded imported and evaluated.
    beware that if you insert some destructive piece of code, it will be executed anyway.
    Use at own risk.

    Args:
        name:
            name of parameter used in the grid Python code.
            This is evaluated as a parameter and you can use it throughout
            the rest of the function

            Examples::

                name = 'lnM_1'

        parameter_name:
            name of the parameter in binary_c

            This name must correspond to a Python variable of the same name,
            which is automatic if parameter_name == name.

            Note: if parameter_name != name, you must set a
                  variable in "precode" or "postcode" to define a Python variable
                  called parameter_name

        longname:
            Long name of parameter

            Examples::

                longname = 'Primary mass'

        range:
            Range of values to take. Does not get used really, the samplerfunc is used to
            get the values from

            Examples::

                range = [math.log(m_min), math.log(m_max)]

        samplerfunc:
            Function returning a list or numpy array of samples spaced appropriately.
            You can either use a real function, or a string representation of a function call.

            Examples::

                samplerfunc = "self.const_linear(math.log(m_min), math.log(m_max), {})".format(resolution['M_1'])

        precode:
            Extra room for some code. This code will be evaluated within the loop of the
            sampling function (i.e. a value for lnM_1 is chosen already)

            Examples::

                precode = 'M_1=math.exp(lnM_1);'

        postcode:
            Code executed after the probability is calculated.

        probdist:
            Function determining the probability that gets assigned to the sampled parameter

            Examples::

                probdist = 'self.Kroupa2001(M_1)*M_1'

        dphasevol:
            part of the parameter space that the total probability is calculated with. Put to -1
            if you want to ignore any dphasevol calculations and set the value to 1

            Examples::

                dphasevol = 'dlnM_1'

        condition:
            condition that has to be met in order for the grid generation to continue

            Examples::

                condition = "self.grid_options['binary']==1"

        gridtype:
            Method on how the value range is sampled. Can be either 'edge' (steps starting at
            the lower edge of the value range) or 'centred'
            (steps starting at ``lower edge + 0.5 * stepsize``).

        dry_parallel:
            If True, try to parallelize this variable in dry runs.

        topcode:
            Code added at the very top of the block.

        bottomcode:
            Code added at the very bottom of the block.

All the distribution functions that we can use are stored in the binarycpython.utils.distribution_functions or binarycpython/utils/distribution_functions.py on git. If you uncomment the help statement below you can see which functions are available now:

[5]:
# import binarycpython.utils.distribution_functions
# help(binarycpython.utils.distribution_functions)

The next cell contains an example of adding the mass grid variable, but sampling in log mass. The commented grid variables are examples of the mass ratio sampling and the period sampling.

[6]:
# Add grid variables
resolution = {"M_1": 20}

# Mass
example_pop.add_grid_variable(
    name="lnm1",
    longname="Primary mass",
    valuerange=[2, 150],
    samplerfunc="self.const_linear(math.log(2), math.log(150), {})".format(resolution["M_1"]),
    precode="M_1=math.exp(lnm1)",
    probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
    dphasevol="dlnm1",
    parameter_name="M_1",
    condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
)

# # Mass ratio
# test_pop.add_grid_variable(
#     name="q",
#     longname="Mass ratio",
#     valuerange=["0.1/M_1", 1],
#     samplerfunc="self.const_linear(0.1/M_1, 1, {})".format(resolution['q']),
#     probdist="self.flatsections(q, [{'min': 0.1/M_1, 'max': 1.0, 'height': 1}])",
#     dphasevol="dq",
#     precode="M_2 = q * M_1",
#     parameter_name="M_2",
#     condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
# )

# #
# test_pop.add_grid_variable(
#    name="log10per", # in days
#    longname="log10(Orbital_Period)",
#    valuerange=[0.15, 5.5],
#    samplerfunc="self.const_linear(0.15, 5.5, {})".format(resolution["per"]),
#    precode="""orbital_period = 10** log10per
# sep = calc_sep_from_period(M_1, M_2, orbital_period)
# sep_min = calc_sep_from_period(M_1, M_2, 10**0.15)
# sep_max = calc_sep_from_period(M_1, M_2, 10**5.5)""",
#    probdist="self.sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**5.5), -0.55)",
#    parameter_name="orbital_period",
#    dphasevol="dlog10per",
# )

Setting logging and handling the output

On default, binary_c will not output anything (except for ‘SINGLE STAR LIFETIME’). It is up to us to determine what will be printed. We can either do that by hardcoding the print statements into binary_c (see documentation binary_c). Or, we can use the custom logging functionality of binarycpython (see notebook notebook_custom_logging.ipynb), which is faster to set up and requires no recompilation of binary_c, but is somewhat more limited in its functionality.

After configuring what will be printed, we need to make a function to parse the output. This can be done by setting the parse_function parameter in the population object (see also notebook notebook_individual_systems.ipynb).

In the code below we will set up both the custom logging, and a parse function to handle that output

[7]:
# Create custom logging statement: in this case we will log when the star turns into a compact object, and then terminate the evolution.
custom_logging_code = """
if(stardata->star[0].stellar_type >= 13)
{
    if (stardata->model.time < stardata->model.max_evolution_time)
    {
        Printf("EXAMPLE_COMPACT_OBJECT %30.12e %g %g %g %d\\n",
            //
            stardata->model.time, // 1
            stardata->star[0].mass, // 2
            stardata->common.zero_age.mass[0], // 3
            stardata->model.probability, // 4
            stardata->star[0].stellar_type // 5
      );
    };
    /* Kill the simulation to save time */
    stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;
};
"""

example_pop.set(
    C_logging_code=custom_logging_code
)
adding: C_logging_code=
if(stardata->star[0].stellar_type >= 13)
{
    if (stardata->model.time < stardata->model.max_evolution_time)
    {
        Printf("EXAMPLE_COMPACT_OBJECT %30.12e %g %g %g %d\n",
            //
            stardata->model.time, // 1
            stardata->star[0].mass, // 2
            stardata->common.zero_age.mass[0], // 3
            stardata->model.probability, // 4
            stardata->star[0].stellar_type // 5
      );
    };
    /* Kill the simulation to save time */
    stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;
};
 to grid_options

The parse function must now catch lines that start with “EXAMPLE_COMPACT_OBJECT”, and write that line to a file

[8]:
def parse_function(self, output):
    """
    Example parse function
    """

    # get info from the population instance
    data_dir = self.custom_options["data_dir"]
    base_filename = self.custom_options["base_filename"]

    # Check directory, make if necessary
    os.makedirs(data_dir, exist_ok=True)

    seperator = " "

    # Create filename
    outfilename = os.path.join(data_dir, base_filename)

    parameters = ["time", "mass", "zams_mass", "probability", "stellar_type"]

    # Go over the output.
    for line in output.splitlines():
        headerline = line.split()[0]

        # CHeck the header and act accordingly
        if headerline == "EXAMPLE_COMPACT_OBJECT":
            values = line.split()[1:]
            print(line)

            if not len(parameters) == len(values):
                print("Number of column names isnt equal to number of columns")
                raise ValueError

            if not os.path.exists(outfilename):
                with open(outfilename, "w") as f:
                    f.write(seperator.join(parameters) + "\n")

            with open(outfilename, "a") as f:
                f.write(seperator.join(values) + "\n")

# Add the parsing function
example_pop.set(
    parse_function=parse_function,
)
adding: parse_function=<function parse_function at 0x7f2b6ca163a0> to grid_options

Evolving the grid

Now that we configured all the main parts of the population object, we can actually run the population! Doing this is straightforward: example_pop.evolve()

This will start up the processing of all the systems. We can control how many cores are used by settings num_cores. By setting the verbosity of the population object to a higher value we can get a lot of verbose information about the run, but for now we will set it to 0.

There are many grid_options that can lead to different behaviour of the evolution of the grid. Please do have a look at those: grid options docs, and try

[9]:
# change verbosity
example_pop.set(verbosity=0)

## Executing a population
## This uses the values generated by the grid_variables
analytics = example_pop.evolve()  # TODO: update this function call
adding: verbosity=0 to grid_options
Do dry run? True
Doing dry run to calculate total starcount and probability
Grid has handled 19 stars with a total probability of 0.0443872
**********************************
*             Dry run            *
*      Total starcount is 19     *
* Total probability is 0.0443872 *
**********************************

EXAMPLE_COMPACT_OBJECT             3.598268106227e+01 1.30592 8.75988 0.00193614 13
EXAMPLE_COMPACT_OBJECT             2.436983545111e+01 1.35842 10.9948 0.00144093 13
EXAMPLE_COMPACT_OBJECT             1.690157944401e+01 1.43124 13.7998 0.00107238 13
EXAMPLE_COMPACT_OBJECT             1.242397939068e+01 1.52416 17.3205 0.000798096 13
EXAMPLE_COMPACT_OBJECT             9.756794139032e+00 1.66914 21.7394 0.000593966 13
EXAMPLE_COMPACT_OBJECT             8.401414766976e+00 1.73729 27.2857 0.000442046 13
EXAMPLE_COMPACT_OBJECT             7.536373523810e+00 1.80677 34.247 0.000328983 13
EXAMPLE_COMPACT_OBJECT             7.393982410080e+00 1.82164 42.9844 0.000244839 13
EXAMPLE_COMPACT_OBJECT             7.396470605248e+00 1.82129 53.9508 0.000182216 13
Do join of subprocesses ...
EXAMPLE_COMPACT_OBJECT             7.399005684057e+00 1.82041 67.7151 0.00013561 13
EXAMPLE_COMPACT_OBJECT             7.443375325717e+00 1.81645 84.9909 0.000100925 13
EXAMPLE_COMPACT_OBJECT             7.451195752942e+00 1.81559 106.674 7.51114e-05 13
EXAMPLE_COMPACT_OBJECT             7.452661646076e+00 1.81543 133.89 5.59e-05 13
Joined subprocesses.
**********************************************************
*  Population-ce756bb317f64099a459bf8b55a746ac finished! *
*           The total probability is 0.0443872.          *
*  It took a total of 0.73s to run 19 systems on 2 cores *
*                   = 1.46s of CPU time.                 *
*              Maximum memory use 293.406 MB             *
**********************************************************

No failed systems were found in this run.
Do analytics
Added analytics to metadata

After the run is complete, some technical report on the run is returned. I stored that in analytics. As we can see below, this dictionary is like a status report of the evolution. Useful for e.g. debugging.

[10]:
print(analytics)
{'population_id': 'ce756bb317f64099a459bf8b55a746ac', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 0.044387171445641534, 'total_count': 19, 'start_timestamp': 1646563001.7193637, 'end_timestamp': 1646563002.4480088, 'time_elapsed': 0.7286450862884521, 'total_mass_run': 649.905447944397, 'total_probability_weighted_mass_run': 0.28133908148630704, 'zero_prob_stars_skipped': 0}

Noteworthy functionality

Some extra features that are available from via the population object are: - write_binary_c_calls_to_file: Function to write the calls that would be passed to binary_c to a file

[11]:
help(example_pop.write_binary_c_calls_to_file)
Help on method write_binary_c_calls_to_file in module binarycpython.utils.population_extensions.dataIO:

write_binary_c_calls_to_file(output_dir: Optional[str] = None, output_filename: Optional[str] = None, include_defaults: bool = False, encoding='utf-8') -> None method of binarycpython.utils.grid.Population instance
    Function that loops over the grid code and writes the generated parameters to a file.
    In the form of a command line call

    Only useful when you have a variable grid as system_generator. MC wouldn't be that useful

    Also, make sure that in this export there are the basic parameters
    like m1,m2,sep, orb-per, ecc, probability etc.

    On default this will write to the datadir, if it exists

    Args:
        output_dir: (optional, default = None) directory where to write the file to. If custom_options['data_dir'] is present, then that one will be used first, and then the output_dir
        output_filename: (optional, default = None) filename of the output. If not set it will be called "binary_c_calls.txt"
        include_defaults: (optional, default = None) whether to include the defaults of binary_c in the lines that are written. Beware that this will result in very long lines, and it might be better to just export the binary_c defaults and keep them in a separate file.

    Returns:
        filename: filename that was used to write the calls to

[12]:
example_pop.set(verbosity=1)
calls_filename = example_pop.write_binary_c_calls_to_file()
print(calls_filename)

with open(calls_filename, 'r') as f:
    print('\n'.join(f.read().splitlines()[:4]))
print("(abridged)")
Generating grid code
Save grid code to grid_options
Write grid code to /tmp/binary_c_python-david/notebooks/notebook_population/binary_c_grid_ce756bb317f64099a459bf8b55a746ac.py [dry_run = False]
Symlinked grid code to /tmp/binary_c_python-david/notebooks/notebook_population/binary_c_grid-latest2
Load grid code function from /tmp/binary_c_python-david/notebooks/notebook_population/binary_c_grid_ce756bb317f64099a459bf8b55a746ac.py
Grid code loaded
Writing binary_c calls to /tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/binary_c_calls.txt
Grid has handled 19 stars with a total probability of 0.0443872
/tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/binary_c_calls.txt
binary_c M_1 2.2406484012210224 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191106 probability 0.011394572976608001
binary_c M_1 2.812296769855663 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191117 probability 0.008480166685456411
binary_c M_1 3.5297876799548944 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191106 probability 0.006311182276049824
binary_c M_1 4.430329401616038 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191106 probability 0.004696962123378559
(abridged)

Full examples of population scripts

Single star population

Below is a full setup for a population of single stars

[13]:
def parse_function(self, output):
    """
    Example parsing function
    """

    # extract info from the population instance

    # Get some information from the
    data_dir = self.custom_options["data_dir"]
    base_filename = self.custom_options["base_filename"]

    # Check directory, make if necessary
    os.makedirs(data_dir, exist_ok=True)

    #
    seperator = " "

    # Create filename
    outfilename = os.path.join(data_dir, base_filename)

    # The header columns matching this
    parameters = ["time", "mass", "zams_mass", "probability", "radius", "stellar_type"]

    # Go over the output.
    for el in output.splitlines():
        headerline = el.split()[0]

        # CHeck the header and act accordingly
        if headerline == "MY_STELLAR_DATA":
            values = el.split()[1:]

            if not len(parameters) == len(values):
                print("Number of column names isnt equal to number of columns")
                raise ValueError

            if not os.path.exists(outfilename):
                with open(outfilename, "w") as f:
                    f.write(seperator.join(parameters) + "\n")

            with open(outfilename, "a") as f:
                f.write(seperator.join(values) + "\n")


# Create population object
example_pop = Population()

# If you want verbosity, set this before other things
example_pop.set(verbosity=0)

# Setting values can be done via .set(<parameter_name>=<value>)
example_pop.set(
    # binary_c physics options
    M_1=10,  # bse_options
    separation=0,  # bse_options
    orbital_period=45000000080,  # bse_options
    max_evolution_time=15000,  # bse_options
    eccentricity=0.02,  # bse_options

    # grid_options
    num_cores=2,
    tmp_dir=TMP_DIR,

    # Custom options: the data directory and the output filename
    data_dir=os.path.join(
        TMP_DIR, "example_python_population_result"
    ),  # custom_options
    base_filename="example_pop.dat",  # custom_options
)

# Creating a parsing function
example_pop.set(
    parse_function=parse_function,  # Setting the parse function thats used in the evolve_population
)

### Custom logging
# Log the moment when the star turns into neutron
example_pop.set(
    C_logging_code="""
if(stardata->star[0].stellar_type >= 13)
{
    if (stardata->model.time < stardata->model.max_evolution_time)
    {
        Printf("MY_STELLAR_DATA %30.12e %g %g %g %g %d\\n",
            //
            stardata->model.time, // 1
            stardata->star[0].mass, // 2
            stardata->common.zero_age.mass[0], // 4
            stardata->model.probability, // 5
            stardata->star[0].radius, // 6
            stardata->star[0].stellar_type // 7
      );
    };
    /* Kill the simulation to save time */
    stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;
};
"""
)

# Add grid variables
resolution = {"M_1": 20}

# Mass
example_pop.add_grid_variable(
    name="lnm1",
    longname="Primary mass",
    valuerange=[2, 150],
    samplerfunc="self.const_linear(math.log(2), math.log(150), {})".format(resolution["M_1"]),
    precode="M_1=math.exp(lnm1)",
    probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
    dphasevol="dlnm1",
    parameter_name="M_1",
    condition="",
)

# Exporting of all the settings can be done with .export_all_info()
example_pop.export_all_info()

# remove the result file if it exists
if os.path.isfile(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat")):
    os.remove(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"))


# Evolve the population
example_pop.evolve()

#
with open(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"), 'r') as f:
    output = f.read()
print("\n")
print(output)
<<<< Warning: Key does not match previously known parameter:                     adding: data_dir=/tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result to custom_options >>>>
<<<< Warning: Key does not match previously known parameter:                     adding: base_filename=example_pop.dat to custom_options >>>>
ok
File at /tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/example_pop_settings.json.gz already exists: cannot write to it
ok
ok pre
ok ret
Do dry run? True
Doing dry run to calculate total starcount and probability
Grid has handled 19 stars with a total probability of 0.0443872
**********************************
*             Dry run            *
*      Total starcount is 19     *
* Total probability is 0.0443872 *
**********************************

Do join of subprocesses ...
Joined subprocesses.
**********************************************************
*  Population-0fa4c2b8707741a5ab41d209ef95a3a4 finished! *
*           The total probability is 0.0443872.          *
*  It took a total of 0.61s to run 19 systems on 2 cores *
*                   = 1.23s of CPU time.                 *
*              Maximum memory use 299.531 MB             *
**********************************************************

No failed systems were found in this run.
Do analytics
Added analytics to metadata


time mass zams_mass probability radius stellar_type
3.598268106227e+01 1.30592 8.75988 0.00193614 1.72498e-05 13
2.436983545111e+01 1.35842 10.9948 0.00144093 1.72498e-05 13
1.690157944401e+01 1.43124 13.7998 0.00107238 1.72498e-05 13
1.242397939068e+01 1.52416 17.3205 0.000798096 1.72498e-05 13
9.756794139032e+00 1.66914 21.7394 0.000593966 1.72498e-05 13
8.401414766976e+00 1.73729 27.2857 0.000442046 1.72498e-05 13
7.536373523810e+00 1.80677 34.247 0.000328983 1.72498e-05 13
7.393982410080e+00 1.82164 42.9844 0.000244839 1.72498e-05 13
7.396470605248e+00 1.82129 53.9508 0.000182216 1.72498e-05 13
7.399005684057e+00 1.82041 67.7151 0.00013561 1.72498e-05 13
7.451195752942e+00 1.81559 106.674 7.51114e-05 1.72498e-05 13
7.443375325717e+00 1.81645 84.9909 0.000100925 1.72498e-05 13
7.452661646076e+00 1.81543 133.89 5.59e-05 1.72498e-05 13

Binary star population

We can also set up a population that samples binary systems, by adding extra grid variables. Below is an example of a full script that runs a binary population and registers when a double compact object is formed. The logging is rather compact and should be expanded to be more useful. Also note that we run very little systems in the following example, as its just intended to show how the code works.

[14]:
def parse_function(self, output):
    """
    Example parsing function
    """

    # extract info from the population instance

    # Get some information from the
    data_dir = self.custom_options["data_dir"]
    base_filename = self.custom_options["base_filename"]

    # Check directory, make if necessary
    os.makedirs(data_dir, exist_ok=True)

    #
    seperator = " "

    # Create filename
    outfilename = os.path.join(data_dir, base_filename)

    # The header columns matching this
    parameters = [
        "time",
        "mass_1", "zams_mass_1", "mass_2", "zams_mass_2",
        "stellar_type_1", "prev_stellar_type_1", "stellar_type_2", "prev_stellar_type_2",
        "metallicity", "probability"
    ]

    # Go over the output.
    for el in output.splitlines():
        headerline = el.split()[0]

        # CHeck the header and act accordingly
        if headerline == "EXAMPLE_DCO":
            values = el.split()[1:]

            if not len(parameters) == len(values):
                print("Number of column names isnt equal to number of columns")
                raise ValueError

            if not os.path.exists(outfilename):
                with open(outfilename, "w") as f:
                    f.write(seperator.join(parameters) + "\n")

            with open(outfilename, "a") as f:
                f.write(seperator.join(values) + "\n")


# Create population object
example_pop = Population()

# If you want verbosity, set this before other things
example_pop.set(verbosity=0)

# Setting values can be done via .set(<parameter_name>=<value>)
example_pop.set(
    # binary_c physics options
    M_1=10,  # bse_options
    separation=0,  # bse_options
    orbital_period=45000000080,  # bse_options
    max_evolution_time=15000,  # bse_options
    eccentricity=0.02,  # bse_options

    # grid_options
    num_cores=2,  # grid_options
    tmp_dir=TMP_DIR,

    # Custom options: the data directory and the output filename
    data_dir=os.path.join(
        TMP_DIR, "example_python_population_result"
    ),  # custom_options
    base_filename="example_pop.dat",  # custom_options
)

# Creating a parsing function
example_pop.set(
    parse_function=parse_function,  # Setting the parse function thats used in the evolve_population
)

### Custom logging
# Log the moment when the star turns into neutron
example_pop.set(
    C_logging_code="""
// logger to find gravitational wave progenitors
if(stardata->star[0].stellar_type>=NS && stardata->star[1].stellar_type>=NS)
{
    if (stardata->model.time < stardata->model.max_evolution_time)
    {
        Printf("EXAMPLE_DCO %30.12e " // 1
            "%g %g %g %g " // 2-5
            "%d %d %d %d " // 6-9
            "%g %g\\n", // 10-11

            stardata->model.time, // 1

            stardata->star[0].mass, //2
            stardata->common.zero_age.mass[0], //3
            stardata->star[1].mass, //4
            stardata->common.zero_age.mass[1], //5

            stardata->star[0].stellar_type, //6
            stardata->previous_stardata->star[0].stellar_type, //7
            stardata->star[1].stellar_type, //8
            stardata->previous_stardata->star[1].stellar_type, //9

            // model stuff
            stardata->common.metallicity, //10
            stardata->model.probability //11
        );
    }
    /* Kill the simulation to safe time */
    stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;
}
"""
)

# Add grid variables
resolution = {"M_1": 3, "q": 3, "per": 3}

# Mass
example_pop.add_grid_variable(
    name="lnm1",
    longname="Primary mass",
    valuerange=[2, 150],
    samplerfunc="self.const_linear(math.log(2), math.log(150), {})".format(resolution["M_1"]),
    precode="M_1=math.exp(lnm1)",
    probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
    dphasevol="dlnm1",
    parameter_name="M_1",
    condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
)

# Mass ratio
example_pop.add_grid_variable(
    name="q",
    longname="Mass ratio",
    valuerange=["0.1/M_1", 1],
    samplerfunc="self.const_linear(0.1/M_1, 1, {})".format(resolution['q']),
    probdist="self.flatsections(q, [{'min': 0.1/M_1, 'max': 1.0, 'height': 1}])",
    dphasevol="dq",
    precode="M_2 = q * M_1",
    parameter_name="M_2",
    condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
)

#
example_pop.add_grid_variable(
   name="log10per", # in days
   longname="log10(Orbital_Period)",
   valuerange=[0.15, 5.5],
   samplerfunc="self.const_linear(0.15, 5.5, {})".format(resolution["per"]),
   precode="""orbital_period = 10** log10per
sep = calc_sep_from_period(M_1, M_2, orbital_period)
sep_min = calc_sep_from_period(M_1, M_2, 10**0.15)
sep_max = calc_sep_from_period(M_1, M_2, 10**5.5)""",
   probdist="self.sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**5.5), -0.55)",
   parameter_name="orbital_period",
   dphasevol="dlog10per",
)

# Exporting of all the settings can be done with .export_all_info()
example_pop.export_all_info()

# remove the result file if it exists
if os.path.isfile(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat")):
    os.remove(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"))

# Evolve the population
example_pop.evolve()

#
with open(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"), 'r') as f:
    output = f.read()
print("\n")
print(output)
<<<< Warning: Key does not match previously known parameter:                     adding: data_dir=/tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result to custom_options >>>>
<<<< Warning: Key does not match previously known parameter:                     adding: base_filename=example_pop.dat to custom_options >>>>
ok
File at /tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/example_pop_settings.json.gz already exists: cannot write to it
ok
ok pre
ok ret
Do dry run? True
Doing dry run to calculate total starcount and probability
Grid has handled 8 stars with a total probability of 0.0211592
**********************************
*             Dry run            *
*      Total starcount is 8      *
* Total probability is 0.0211592 *
**********************************

Do join of subprocesses ...
Joined subprocesses.
**********************************************************
*  Population-0eb5c0c9abd34607a6ee060b26a7e32f finished! *
*           The total probability is 0.0211592.          *
*  It took a total of 0.84s to run 8 systems on 2 cores  *
*                   = 1.68s of CPU time.                 *
*              Maximum memory use 300.125 MB             *
**********************************************************

No failed systems were found in this run.
Do analytics
Added analytics to metadata


time mass_1 zams_mass_1 mass_2 zams_mass_2 stellar_type_1 prev_stellar_type_1 stellar_type_2 prev_stellar_type_2 metallicity probability
1.378266748188e+01 1.66293 50.9713 1.78767 12.8178 13 13 13 8 0.02 0.000339963
1.817608462595e+01 1.82104 50.9713 1.41436 12.8178 13 13 13 5 0.02 0.000193036
7.422997711686e+00 1.82479 50.9713 1.82171 38.2535 13 13 13 8 0.02 0.000193036
1.205711924468e+01 1.73765 50.9713 0 38.2535 13 13 15 8 0.02 0.000339963