Skip to content
Snippets Groups Projects
notebook_population.ipynb 50.4 KiB
Newer Older
{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bbbaafbb-fd7d-4b73-a970-93506ba35d71",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Tutorial: Running populations with binary_c-python\n",
    "This notebook will show you how to evolve a population of stars\n",
    "\n",
    "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. \n",
    "\n",
    "At the bottom of this notebook there are some complete example scripts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bf6b8673-a2b5-4b50-ad1b-e90671f57470",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
David Hendriks's avatar
David Hendriks committed
    "\n",
    "from binarycpython.utils.custom_logging_functions import temp_dir\n",
    "from binarycpython.utils.grid import Population\n",
David Hendriks's avatar
David Hendriks committed
    "TMP_DIR = temp_dir(\"notebooks\", \"notebook_population\")\n",
    "# help(Population) # Uncomment to see the public functions of this object"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f268eff3-4e08-4f6b-8b59-f22dba4d2074",
   "metadata": {},
   "source": [
    "## Setting up the Population object\n",
    "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.\n",
    "\n",
    "There are three categories of options that the Population object can set:\n",
    "- 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. \n",
    "- 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. \n",
    "\n",
    "- 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>']`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "79ab50b7-591f-4883-af09-116d1835a751",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adding: M_1=10 to BSE_options\n",
      "adding: orbital_period=45000000080 to BSE_options\n",
      "adding: max_evolution_time=15000 to BSE_options\n",
      "adding: eccentricity=0.02 to BSE_options\n",
      "adding: num_cores=2 to grid_options\n",
      "adding: tmp_dir=/tmp/binary_c_python-izzard/notebooks/notebook_population to grid_options\n",
      "<<<< Warning: Key does not match previously known parameter:                     adding: data_dir=/tmp/binary_c_python-izzard/notebooks/notebook_population/example_python_population_result to custom_options >>>>\n",
      "<<<< Warning: Key does not match previously known parameter:                     adding: base_filename=example_pop.dat to custom_options >>>>\n",
      "1\n",
      "example_pop.dat\n",
      "10\n"
     ]
    }
   ],
   "source": [
    "# Create population object\n",
    "example_pop = Population()\n",
    "\n",
    "# If you want verbosity, set this before other things\n",
    "example_pop.set(verbosity=1)\n",
    "\n",
    "# Setting values can be done via .set(<parameter_name>=<value>)\n",
    "# Values that are known to be binary_c_parameters are loaded into bse_options.\n",
    "# Those that are present in the default grid_options are set in grid_options\n",
    "# All other values that you set are put in a custom_options dict\n",
    "example_pop.set(\n",
    "    # binary_c physics options\n",
    "    M_1=10,  # bse_options\n",
    "    orbital_period=45000000080,  # bse_options\n",
    "    max_evolution_time=15000,  # bse_options\n",
    "    eccentricity=0.02,  # bse_options\n",
    "\n",
    "\n",
    "    # grid_options\n",
David Hendriks's avatar
David Hendriks committed
    "    tmp_dir=TMP_DIR,\n",
    "    \n",
    "    # Custom options # TODO: need to be set in grid_options probably\n",
    "    data_dir=os.path.join(\n",
David Hendriks's avatar
David Hendriks committed
    "        TMP_DIR, \"example_python_population_result\"\n",
    "    ),  # custom_options\n",
    "    base_filename=\"example_pop.dat\",  # custom_options\n",
    ")\n",
    "\n",
    "# We can use the options through\n",
    "print(example_pop.grid_options['verbosity'])\n",
    "print(example_pop.custom_options['base_filename'])\n",
    "print(example_pop.bse_options['M_1'])"
   "cell_type": "markdown",
   "id": "f8d46d19-633d-4911-821d-a59daed31816",
   "metadata": {},
   "source": [
    "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()`.\n",
    "On default this exports everything, each of the sections can be disabled:\n",
    "  - population settings (bse_options, grid_options, custom_options), turn off with include_population\n",
    "      settings=False\n",
    "  - binary_c_defaults (all the commandline arguments that binary c accepts, and their defaults).\n",
    "      turn off with include_binary_c_defaults=False\n",
    "  - include_binary_c_version_info (all the compilation info, and information about the compiled\n",
    "      parameters), turn off with include_binary_c_version_info=False\n",
    "  - include_binary_c_help_all (all the help information for all the binary_c parameters),\n",
    "      turn off with include_binary_c_help_all=Fase\n",
    "      \n",
    "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=<>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b9c2471a-a5b0-48b7-a50b-2f0d22100926",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing settings to /tmp/binary_c_python-izzard/notebooks/notebook_population/example_python_population_result/example_pop_settings.json\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'/tmp/binary_c_python-izzard/notebooks/notebook_population/example_python_population_result/example_pop_settings.json'"
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "example_pop.export_all_info()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9a65554-36ab-4a04-96ca-9f1422c307fd",
   "metadata": {},
   "source": [
    "## Adding grid variables\n",
    "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.\n",
    "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. \n",
    "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) \n",
    "A notable special type of grid variable is that of the Moe & di Stefano 2017 dataset (see further down in the notebook).\n",
    "To add a grid variable to the population object we use `example_pop.add_grid_variable` (see next cell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "68c84521-9ae8-4020-af7a-5334173db969",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method add_grid_variable in module binarycpython.utils.grid:\n",
      "\n",
      "add_grid_variable(name: str, parameter_name: str, longname: str, valuerange: Union[list, str], samplerfunc: str, probdist: str, dphasevol: Union[str, int], 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) -> None method of binarycpython.utils.grid.Population instance\n",
      "    Function to add grid variables to the grid_options.\n",
      "    \n",
      "    The execution of the grid generation will be through a nested for loop.\n",
      "    Each of the grid variables will get create a deeper for loop.\n",
      "    \n",
      "    The real function that generates the numbers will get written to a new file in the TMP_DIR,\n",
      "    and then loaded imported and evaluated.\n",
      "    beware that if you insert some destructive piece of code, it will be executed anyway.\n",
      "    Use at own risk.\n",
      "    \n",
      "    Tasks:\n",
      "        - TODO: Fix this complex function.\n",
      "    \n",
      "    Args:\n",
      "        name:\n",
      "            name of parameter used in the grid Python code.\n",
      "            This is evaluated as a parameter and you can use it throughout\n",
      "            the rest of the function\n",
David Hendriks's avatar
David Hendriks committed
      "    \n",
      "            Examples:\n",
      "                name = 'lnm1'\n",
      "    \n",
      "        parameter_name:\n",
      "            name of the parameter in binary_c\n",
      "    \n",
      "            This name must correspond to a Python variable of the same name,\n",
      "            which is automatic if parameter_name == name.\n",
      "    \n",
      "            Note: if parameter_name != name, you must set a\n",
      "                  variable in \"precode\" or \"postcode\" to define a Python variable\n",
      "                  called parameter_name\n",
      "    \n",
      "        longname:\n",
      "            Long name of parameter\n",
David Hendriks's avatar
David Hendriks committed
      "    \n",
      "            Examples:\n",
      "                longname = 'Primary mass'\n",
      "        range:\n",
      "            Range of values to take. Does not get used really, the samplerfunc is used to\n",
      "            get the values from\n",
David Hendriks's avatar
David Hendriks committed
      "    \n",
      "            Examples:\n",
      "                range = [math.log(m_min), math.log(m_max)]\n",
      "        samplerfunc:\n",
      "            Function returning a list or numpy array of samples spaced appropriately.\n",
      "            You can either use a real function, or a string representation of a function call.\n",
      "    \n",
      "            Examples:\n",
      "                samplerfunc = \"const(math.log(m_min), math.log(m_max), {})\".format(resolution['M_1'])\n",
      "    \n",
      "        precode:\n",
      "            Extra room for some code. This code will be evaluated within the loop of the\n",
      "            sampling function (i.e. a value for lnm1 is chosen already)\n",
David Hendriks's avatar
David Hendriks committed
      "    \n",
      "            Examples:\n",
      "                precode = 'M_1=math.exp(lnm1);'\n",
      "        postcode:\n",
      "            Code executed after the probability is calculated.\n",
      "        probdist:\n",
      "            Function determining the probability that gets assigned to the sampled parameter\n",
David Hendriks's avatar
David Hendriks committed
      "    \n",
      "            Examples:\n",
      "                probdist = 'Kroupa2001(M_1)*M_1'\n",
      "        dphasevol:\n",
      "            part of the parameter space that the total probability is calculated with. Put to -1\n",
      "            if you want to ignore any dphasevol calculations and set the value to 1\n",
      "            Examples:\n",
      "                dphasevol = 'dlnm1'\n",
      "        condition:\n",
      "            condition that has to be met in order for the grid generation to continue\n",
      "            Examples:\n",
      "                condition = 'self.grid_options['binary']==1'\n",
      "        gridtype:\n",
      "            Method on how the value range is sampled. Can be either 'edge' (steps starting at\n",
David Hendriks's avatar
David Hendriks committed
      "            the lower edge of the value range) or 'centred'\n",
      "            (steps starting at lower edge + 0.5 * stepsize).\n",
      "    \n",
      "        topcode:\n",
      "            Code added at the very top of the block.\n",
      "    \n",
      "        bottomcode:\n",
      "            Code added at the very bottom of the block.\n",
    "help(example_pop.add_grid_variable)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd75cebe-2152-4025-b680-dc020b80889b",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "048db541-3e92-4c5d-a25c-9c5a34b9c857",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import binarycpython.utils.distribution_functions\n",
    "# help(binarycpython.utils.distribution_functions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b3a007b-5c17-42a7-a981-7e268e6f545c",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "47979841-2c26-4b26-8945-603d013dc93a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Added grid variable: {\n",
      "    \"name\": \"lnm1\",\n",
      "    \"parameter_name\": \"M_1\",\n",
      "    \"longname\": \"Primary mass\",\n",
      "    \"valuerange\": [\n",
      "        2,\n",
      "        150\n",
      "    ],\n",
      "    \"samplerfunc\": \"const(math.log(2), math.log(150), 20)\",\n",
      "    \"precode\": \"M_1=math.exp(lnm1)\",\n",
      "    \"postcode\": null,\n",
      "    \"probdist\": \"three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1\",\n",
      "    \"dphasevol\": \"dlnm1\",\n",
      "    \"condition\": \"\",\n",
David Hendriks's avatar
David Hendriks committed
      "    \"gridtype\": \"centred\",\n",
      "    \"branchpoint\": 0,\n",
      "    \"branchcode\": null,\n",
      "    \"topcode\": null,\n",
      "    \"bottomcode\": null,\n",
      "    \"grid_variable_number\": 0\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# Add grid variables\n",
    "resolution = {\"M_1\": 20}\n",
    "\n",
    "# Mass\n",
    "example_pop.add_grid_variable(\n",
    "    name=\"lnm1\",\n",
    "    longname=\"Primary mass\",\n",
    "    valuerange=[2, 150],\n",
    "    samplerfunc=\"const(math.log(2), math.log(150), {})\".format(resolution[\"M_1\"]),\n",
    "    precode=\"M_1=math.exp(lnm1)\",\n",
    "    probdist=\"three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1\",\n",
    "    dphasevol=\"dlnm1\",\n",
    "    parameter_name=\"M_1\",\n",
    "    condition=\"\",  # Impose a condition on this grid variable. Mostly for a check for yourself\n",
    ")\n",
    "\n",
    "# # Mass ratio\n",
    "# test_pop.add_grid_variable(\n",
    "#     name=\"q\",\n",
    "#     longname=\"Mass ratio\",\n",
    "#     valuerange=[\"0.1/M_1\", 1],\n",
    "#     samplerfunc=\"const(0.1/M_1, 1, {})\".format(resolution['q']),\n",
    "#     probdist=\"flatsections(q, [{'min': 0.1/M_1, 'max': 1.0, 'height': 1}])\",\n",
    "#     dphasevol=\"dq\",\n",
    "#     precode=\"M_2 = q * M_1\",\n",
    "#     parameter_name=\"M_2\",\n",
    "#     condition=\"\",  # Impose a condition on this grid variable. Mostly for a check for yourself\n",
    "# )\n",
    "\n",
    "# #\n",
    "# test_pop.add_grid_variable(\n",
    "#    name=\"log10per\", # in days\n",
    "#    longname=\"log10(Orbital_Period)\",\n",
    "#    valuerange=[0.15, 5.5],\n",
    "#    samplerfunc=\"const(0.15, 5.5, {})\".format(resolution[\"per\"]),\n",
    "#    precode=\"\"\"orbital_period = 10** log10per\n",
    "# sep = calc_sep_from_period(M_1, M_2, orbital_period)\n",
    "# sep_min = calc_sep_from_period(M_1, M_2, 10**0.15)\n",
    "# sep_max = calc_sep_from_period(M_1, M_2, 10**5.5)\"\"\",\n",
    "#    probdist=\"sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**5.5), -0.55)\",\n",
    "#    parameter_name=\"orbital_period\",\n",
    "#    dphasevol=\"dlog10per\",\n",
    "# )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "163f13ae-fec1-4ee8-b9d4-c1b75c19ff39",
   "metadata": {},
   "source": [
    "## Setting logging and handling the output\n",
    "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. \n",
    "\n",
    "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`). \n",
    "\n",
    "In the code below we will set up both the custom logging, and a parse function to handle that output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0c986215-93b1-4e30-ad79-f7c397e9ff7d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adding: C_logging_code=\n",
      "if(stardata->star[0].stellar_type >= 13)    \n",
      "{\n",
      "    if (stardata->model.time < stardata->model.max_evolution_time)\n",
      "    {\n",
      "        Printf(\"EXAMPLE_COMPACT_OBJECT %30.12e %g %g %g %d\\n\",\n",
      "            // \n",
      "            stardata->model.time, // 1\n",
      "            stardata->star[0].mass, // 2\n",
      "            stardata->common.zero_age.mass[0], // 3\n",
      "            stardata->model.probability, // 4\n",
      "            stardata->star[0].stellar_type // 5\n",
      "      );\n",
      "    };\n",
      "    /* Kill the simulation to save time */\n",
      "    stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;\n",
      "};\n",
      " to grid_options\n"
     ]
    "# Create custom logging statement: in this case we will log when the star turns into a compact object, and then terminate the evolution.\n",
    "custom_logging_statement = \"\"\"\n",
    "if(stardata->star[0].stellar_type >= 13)    \n",
    "{\n",
    "    if (stardata->model.time < stardata->model.max_evolution_time)\n",
    "    {\n",
    "        Printf(\"EXAMPLE_COMPACT_OBJECT %30.12e %g %g %g %d\\\\n\",\n",
    "            // \n",
    "            stardata->model.time, // 1\n",
    "            stardata->star[0].mass, // 2\n",
    "            stardata->common.zero_age.mass[0], // 3\n",
    "            stardata->model.probability, // 4\n",
    "            stardata->star[0].stellar_type // 5\n",
    "      );\n",
    "    };\n",
    "    /* Kill the simulation to save time */\n",
    "    stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;\n",
    "};\n",
    "\"\"\"\n",
    "\n",
    "example_pop.set(\n",
    "    C_logging_code=custom_logging_statement\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae1f1f0c-1f8b-42d8-b051-cbf8c6b51514",
   "metadata": {},
   "source": [
    "The parse function must now catch lines that start with \"EXAMPLE_COMPACT_OBJECT\", and write that line to a file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fd197154-a8ce-4865-8929-008d3483101a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adding: parse_function=<function parse_function at 0x1528ac7290d0> to grid_options\n"
    "def parse_function(self, output):\n",
    "    \"\"\"\n",
    "    Example parse function\n",
    "    \"\"\"\n",
    "    \n",
    "    # get info from the population instance\n",
    "    data_dir = self.custom_options[\"data_dir\"]\n",
    "    base_filename = self.custom_options[\"base_filename\"]\n",
    "\n",
    "    # Check directory, make if necessary\n",
    "    os.makedirs(data_dir, exist_ok=True)\n",
    "\n",
    "    seperator = \" \"\n",
    "\n",
    "    # Create filename\n",
    "    outfilename = os.path.join(data_dir, base_filename)\n",
    "\n",
    "    parameters = [\"time\", \"mass\", \"zams_mass\", \"probability\", \"stellar_type\"]\n",
    "\n",
    "    # Go over the output.\n",
    "    for line in output.splitlines():\n",
    "        headerline = line.split()[0]\n",
    "\n",
    "        # CHeck the header and act accordingly\n",
    "        if headerline == \"EXAMPLE_COMPACT_OBJECT\":\n",
    "            values = line.split()[1:]\n",
    "            print(line)\n",
    "            \n",
    "            if not len(parameters) == len(values):\n",
    "                print(\"Number of column names isnt equal to number of columns\")\n",
    "                raise ValueError\n",
    "\n",
    "            if not os.path.exists(outfilename):\n",
    "                with open(outfilename, \"w\") as f:\n",
    "                    f.write(seperator.join(parameters) + \"\\n\")\n",
    "\n",
    "            with open(outfilename, \"a\") as f:\n",
    "                f.write(seperator.join(values) + \"\\n\")\n",
    "\n",
    "# Add the parsing function\n",
    "example_pop.set(\n",
    "    parse_function=parse_function,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91509ce5-ffe7-4937-aa87-6d7baac9ac04",
   "metadata": {},
   "source": [
    "## Evolving the grid\n",
    "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()`\n",
    "\n",
    "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.\n",
    "\n",
    "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](https://ri0005.pages.surrey.ac.uk/binary_c-python/grid_options_descriptions.html), and try  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8ea376c1-1e92-45af-8cab-9d7fdca564eb",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adding: verbosity=0 to grid_options\n",
      "Doing dry run to calculate total starcount and probability\n",
      "Generating grid code\n",
      "Grid has handled 20 stars with a total probability of 0.0444029\n",
      "**************************************\n",
      "* Total starcount for this run is 20 *\n",
      "*   Total probability is 0.0444029   *\n",
      "**************************************\n",
      "\n",
      "Generating grid code\n",
David Hendriks's avatar
David Hendriks committed
      "EXAMPLE_COMPACT_OBJECT             4.139293101586e+01 1.29427 8.13626 0.00202467 13\n",
      "EXAMPLE_COMPACT_OBJECT             2.802986496151e+01 1.33699 10.0967 0.00152924 13\n",
      "EXAMPLE_COMPACT_OBJECT             1.963621764679e+01 1.39754 12.5294 0.00115504 13\n",
      "EXAMPLE_COMPACT_OBJECT             1.427601421985e+01 1.47745 15.5483 0.000872405 13\n",
      "EXAMPLE_COMPACT_OBJECT             1.094409257247e+01 1.57571 19.2947 0.00065893 13\n",
      "EXAMPLE_COMPACT_OBJECT             9.181971798545e+00 1.68748 23.9436 0.000497691 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.905335716621e+00 1.77287 29.7128 0.000375908 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.451192744924e+00 1.81495 36.872 0.000283924 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.396133472739e+00 1.82088 45.7561 0.000214449 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.396675941641e+00 1.82123 56.7809 0.000161974 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.404641347602e+00 1.82074 70.4621 0.000122339 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.444217227690e+00 1.81636 87.4397 9.2403e-05 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.453317880232e+00 1.81536 108.508 6.97923e-05 13\n",
      "EXAMPLE_COMPACT_OBJECT             7.450828476487e+00 1.81563 134.653 5.27143e-05 13\n",
      "**********************************************************\n",
      "*  Population-50fb66cc659c46c8bbc29fe0c8651c2f finished! *\n",
      "*           The total probability is 0.0444029.          *\n",
      "*  It took a total of 3.30s to run 20 systems on 2 cores *\n",
      "*                   = 6.60s of CPU time.                 *\n",
      "*              Maximum memory use 433.070 MB             *\n",
      "**********************************************************\n",
      "\n",
      "There were no errors found in this run.\n"
     ]
    }
   ],
    "# change verbosity\n",
    "example_pop.set(verbosity=0)\n",
    "\n",
    "## Executing a population\n",
    "## This uses the values generated by the grid_variables\n",
    "analytics = example_pop.evolve()  # TODO: update this function call"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91ab45c7-7d31-4543-aee4-127ab58e891f",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e1f0464b-0424-4022-b34b-5b744bc2c59d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'population_name': '50fb66cc659c46c8bbc29fe0c8651c2f', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 0.04440288843805411, 'total_count': 20, 'start_timestamp': 1635760967.3245144, 'end_timestamp': 1635760970.6249793, 'total_mass_run': 684.2544031669784, 'total_probability_weighted_mass_run': 0.28134439269236855, 'zero_prob_stars_skipped': 0}\n"
     ]
    }
   ],
   "source": [
    "print(analytics)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6460df56-9fba-4817-9a1e-593ef15d98c1",
   "metadata": {},
    "## Noteworthy functionality\n",
    "Some extra features that are available from via the population object are:\n",
    "- write_binary_c_calls_to_file: Function to write the calls that would be passed to binary_c to a file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "83f8e519-4f7c-474a-ad95-f175a34fae81",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method write_binary_c_calls_to_file in module binarycpython.utils.grid:\n",
      "\n",
      "write_binary_c_calls_to_file(output_dir: Optional[str] = None, output_filename: Optional[str] = None, include_defaults: bool = False) -> None method of binarycpython.utils.grid.Population instance\n",
      "    Function that loops over the grid code and writes the generated parameters to a file.\n",
      "    In the form of a command line call\n",
      "    \n",
      "    Only useful when you have a variable grid as system_generator. MC wouldn't be that useful\n",
      "    \n",
      "    Also, make sure that in this export there are the basic parameters\n",
      "    like m1,m2,sep, orb-per, ecc, probability etc.\n",
      "    \n",
      "    On default this will write to the datadir, if it exists\n",
      "    \n",
      "    Tasks:\n",
      "        - TODO: test this function\n",
      "        - TODO: make sure the binary_c_python .. output file has a unique name\n",
      "    \n",
      "    Args:\n",
      "        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\n",
      "        output_filename: (optional, default = None) filename of the output. If not set it will be called \"binary_c_calls.txt\"\n",
      "        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.\n",
      "    \n",
      "    Returns:\n",
      "        filename: filename that was used to write the calls to\n",
      "\n"
    "help(example_pop.write_binary_c_calls_to_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "dacfed75-dfe3-4afd-a0ff-a4be17746021",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating grid code\n",
      "Generating grid code\n",
      "Saving grid code to grid_options\n",
      "Writing grid code to /tmp/binary_c_python-izzard/notebooks/notebook_population/binary_c_grid_50fb66cc659c46c8bbc29fe0c8651c2f.py [dry_run = False]\n",
      "Symlinked grid code to /tmp/binary_c_python-izzard/notebooks/notebook_population/binary_c_grid-latest2 \n",
      "Loading grid code function from /tmp/binary_c_python-izzard/notebooks/notebook_population/binary_c_grid_50fb66cc659c46c8bbc29fe0c8651c2f.py\n",
      "Grid code loaded\n",
      "Writing binary_c calls to /tmp/binary_c_python-izzard/notebooks/notebook_population/example_python_population_result/binary_c_calls.txt\n",
      "Grid has handled 20 stars with a total probability of 0.0444029\n",
      "/tmp/binary_c_python-izzard/notebooks/notebook_population/example_python_population_result/binary_c_calls.txt\n",
David Hendriks's avatar
David Hendriks committed
      "binary_c M_1 2.227955577093495 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.21587440567681548 probability 0.010905083645619543\n",
      "binary_c M_1 2.7647737053496777 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.2158744056768156 probability 0.00823663875514986\n",
      "binary_c M_1 3.430936289925951 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.21587440567681537 probability 0.0062211552141636295\n",
      "binary_c M_1 4.2576084265970895 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.2158744056768156 probability 0.004698855121516281\n"
David Hendriks's avatar
David Hendriks committed
    "example_pop.set(verbosity=1)\n",
    "calls_filename = example_pop.write_binary_c_calls_to_file()\n",
    "print(calls_filename)\n",
    "with open(calls_filename, 'r') as f:\n",
    "    print('\\n'.join(f.read().splitlines()[:4]))"
  },
  {
   "cell_type": "markdown",
   "id": "60359eb1-4d0c-4d2d-8265-ec5171b944a2",
   "metadata": {},
   "source": [
    "## Full examples of population scripts\n",
    "Below is a full setup for a population of single stars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7212b6be-9787-4122-a7f1-86538cf38179",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<<<< Warning: Key does not match previously known parameter:                     adding: data_dir=/tmp/binary_c_python-izzard/notebooks/notebook_population/example_python_population_result to custom_options >>>>\n",
      "<<<< Warning: Key does not match previously known parameter:                     adding: base_filename=example_pop.dat to custom_options >>>>\n",
      "Doing dry run to calculate total starcount and probability\n",
      "Generating grid code\n",
      "Grid has handled 20 stars with a total probability of 0.0444029\n",
      "**************************************\n",
      "* Total starcount for this run is 20 *\n",
      "*   Total probability is 0.0444029   *\n",
      "**************************************\n",
      "\n",
      "Generating grid code\n",
      "**********************************************************\n",
      "*  Population-05e51ba114934b37bab48f1db40b7333 finished! *\n",
      "*           The total probability is 0.0444029.          *\n",
      "*  It took a total of 3.46s to run 20 systems on 2 cores *\n",
      "*                   = 6.93s of CPU time.                 *\n",
      "*              Maximum memory use 437.047 MB             *\n",
      "**********************************************************\n",
      "\n",
      "There were no errors found in this run.\n",
      "\n",
      "\n",
      "time mass zams_mass probability radius stellar_type\n",
David Hendriks's avatar
David Hendriks committed
      "4.139293101586e+01 1.29427 8.13626 0.00202467 1.72498e-05 13\n",
      "2.802986496151e+01 1.33699 10.0967 0.00152924 1.72498e-05 13\n",
      "1.963621764679e+01 1.39754 12.5294 0.00115504 1.72498e-05 13\n",
      "1.427601421985e+01 1.47745 15.5483 0.000872405 1.72498e-05 13\n",
      "1.094409257247e+01 1.57571 19.2947 0.00065893 1.72498e-05 13\n",
      "9.181971798545e+00 1.68748 23.9436 0.000497691 1.72498e-05 13\n",
      "7.905335716621e+00 1.77287 29.7128 0.000375908 1.72498e-05 13\n",
      "7.451192744924e+00 1.81495 36.872 0.000283924 1.72498e-05 13\n",
      "7.396133472739e+00 1.82088 45.7561 0.000214449 1.72498e-05 13\n",
      "7.396675941641e+00 1.82123 56.7809 0.000161974 1.72498e-05 13\n",
      "7.404641347602e+00 1.82074 70.4621 0.000122339 1.72498e-05 13\n",
      "7.444217227690e+00 1.81636 87.4397 9.2403e-05 1.72498e-05 13\n",
      "7.453317880232e+00 1.81536 108.508 6.97923e-05 1.72498e-05 13\n",
      "7.450828476487e+00 1.81563 134.653 5.27143e-05 1.72498e-05 13\n",
   "source": [
    "def parse_function(self, output):\n",
    "    \"\"\"\n",
    "    Example parsing function\n",
    "    \"\"\"\n",
    "    \n",
    "    # extract info from the population instance\n",
    "\n",
    "    # Get some information from the\n",
    "    data_dir = self.custom_options[\"data_dir\"]\n",
    "    base_filename = self.custom_options[\"base_filename\"]\n",
    "\n",
    "    # Check directory, make if necessary\n",
    "    os.makedirs(data_dir, exist_ok=True)\n",
    "\n",
    "    #\n",
    "    seperator = \" \"\n",
    "\n",
    "    # Create filename\n",
    "    outfilename = os.path.join(data_dir, base_filename)\n",
    "\n",
    "    # The header columns matching this \n",
    "    parameters = [\"time\", \"mass\", \"zams_mass\", \"probability\", \"radius\", \"stellar_type\"]\n",
    "\n",
    "    # Go over the output.\n",
    "    for el in output.splitlines():\n",
    "        headerline = el.split()[0]\n",
    "\n",
    "        # CHeck the header and act accordingly\n",
    "        if headerline == \"MY_STELLAR_DATA\":\n",
    "            values = el.split()[1:]\n",
    "\n",
    "            if not len(parameters) == len(values):\n",
    "                print(\"Number of column names isnt equal to number of columns\")\n",
    "                raise ValueError\n",
    "\n",
    "            if not os.path.exists(outfilename):\n",
    "                with open(outfilename, \"w\") as f:\n",
    "                    f.write(seperator.join(parameters) + \"\\n\")\n",
    "\n",
    "            with open(outfilename, \"a\") as f:\n",
    "                f.write(seperator.join(values) + \"\\n\")\n",
    "\n",
    "\n",
    "# Create population object\n",
    "example_pop = Population()\n",
    "\n",
    "# If you want verbosity, set this before other things\n",
    "example_pop.set(verbosity=0)\n",
    "\n",
    "# Setting values can be done via .set(<parameter_name>=<value>)\n",
    "example_pop.set(\n",
    "    # binary_c physics options\n",
    "    M_1=10,  # bse_options\n",
    "    separation=0,  # bse_options\n",
    "    orbital_period=45000000080,  # bse_options\n",
    "    max_evolution_time=15000,  # bse_options\n",
    "    eccentricity=0.02,  # bse_options\n",
    "    \n",
    "    # grid_options\n",
David Hendriks's avatar
David Hendriks committed
    "    tmp_dir=TMP_DIR,\n",
    "\n",
    "    # Custom options: the data directory and the output filename\n",
    "    data_dir=os.path.join(\n",
David Hendriks's avatar
David Hendriks committed
    "        TMP_DIR, \"example_python_population_result\"\n",
    "    ),  # custom_options\n",
    "    base_filename=\"example_pop.dat\",  # custom_options\n",
    ")\n",
    "\n",
    "# Creating a parsing function\n",
    "example_pop.set(\n",
    "    parse_function=parse_function,  # Setting the parse function thats used in the evolve_population\n",
    ")\n",
    "\n",
    "### Custom logging\n",
    "# Log the moment when the star turns into neutron\n",
    "example_pop.set(\n",
    "    C_logging_code=\"\"\"\n",
    "if(stardata->star[0].stellar_type >= 13)    \n",
    "{\n",
    "    if (stardata->model.time < stardata->model.max_evolution_time)\n",
    "    {\n",
    "        Printf(\"MY_STELLAR_DATA %30.12e %g %g %g %g %d\\\\n\",\n",
    "            // \n",
    "            stardata->model.time, // 1\n",
    "            stardata->star[0].mass, // 2\n",
    "            stardata->common.zero_age.mass[0], // 4\n",
    "            stardata->model.probability, // 5\n",
    "            stardata->star[0].radius, // 6\n",
    "            stardata->star[0].stellar_type // 7\n",
    "      );\n",
    "    };\n",
    "    /* Kill the simulation to save time */\n",
    "    stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;\n",
    "};\n",
    "\"\"\"\n",
    ")\n",
    "\n",
    "# Add grid variables\n",
    "resolution = {\"M_1\": 20}\n",
    "\n",
    "# Mass\n",
    "example_pop.add_grid_variable(\n",
    "    name=\"lnm1\",\n",
    "    longname=\"Primary mass\",\n",
    "    valuerange=[2, 150],\n",
    "    samplerfunc=\"const(math.log(2), math.log(150), {})\".format(resolution[\"M_1\"]),\n",
    "    precode=\"M_1=math.exp(lnm1)\",\n",
    "    probdist=\"three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1\",\n",
    "    dphasevol=\"dlnm1\",\n",
    "    parameter_name=\"M_1\",\n",
    "    condition=\"\",\n",
    ")\n",
    "\n",
    "# Exporting of all the settings can be done with .export_all_info()\n",
    "example_pop.export_all_info()\n",
    "\n",
    "# remove the result file if it exists\n",
David Hendriks's avatar
David Hendriks committed
    "if os.path.isfile(os.path.join(TMP_DIR, \"example_python_population_result\", \"example_pop.dat\")):\n",
    "    os.remove(os.path.join(TMP_DIR, \"example_python_population_result\", \"example_pop.dat\"))\n",
    "\n",
    "\n",
    "# Evolve the population\n",
    "example_pop.evolve()\n",
    "\n",
    "# \n",
David Hendriks's avatar
David Hendriks committed
    "with open(os.path.join(TMP_DIR, \"example_python_population_result\", \"example_pop.dat\"), 'r') as f:\n",
    "    output = f.read()\n",
    "print(\"\\n\")\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2ab0979-6575-481d-9c1c-ca98517b2437",
   "metadata": {},
   "source": [
    "We can also set up a population that samples biinary 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 top be more useful"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "79acdbb2-7dd6-45c4-9609-80994f03619a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<<<< Warning: Key does not match previously known parameter:                     adding: data_dir=/tmp/binary_c_python-izzard/notebooks/notebook_population/example_python_population_result to custom_options >>>>\n",
      "<<<< Warning: Key does not match previously known parameter:                     adding: base_filename=example_pop.dat to custom_options >>>>\n",
      "Doing dry run to calculate total starcount and probability\n",
      "Generating grid code\n",
      "Grid has handled 27 stars with a total probability of 0.0248684\n",
      "**************************************\n",
      "* Total starcount for this run is 27 *\n",
      "*   Total probability is 0.0248684   *\n",
      "**************************************\n",
      "\n",
      "Generating grid code\n",
      "**********************************************************\n",
      "*  Population-8bc1eafea1c34b05894c1618639d8c37 finished! *\n",
      "*           The total probability is 0.0248684.          *\n",
      "* It took a total of 16.10s to run 27 systems on 2 cores *\n",
      "*                  = 32.20s of CPU time.                 *\n",
      "*              Maximum memory use 437.695 MB             *\n",
      "**********************************************************\n",
      "\n",
      "There were no errors found in this run.\n",
      "\n",
      "\n",
      "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\n",
David Hendriks's avatar
David Hendriks committed
      "1.219029061236e+01 1.60007 17.3205 0 2.97008 13 5 15 15 0.02 0.000498487\n",
      "1.935920339886e+01 1.29448 17.3205 0 8.71025 13 13 15 2 0.02 0.000498487\n",
      "2.123794969278e+01 1.30902 17.3205 1.58518 8.71025 13 13 13 5 0.02 0.000287968\n",
David Hendriks's avatar
David Hendriks committed
      "3.579099761269e+01 1.52414 17.3205 1.30642 8.71025 13 13 13 5 0.02 0.000220016\n",
      "1.674063083432e+01 1.29457 17.3205 0 14.4504 13 13 15 2 0.02 0.000498487\n",
      "1.548740826516e+01 1.52415 17.3205 1.45407 14.4504 13 13 13 5 0.02 0.000220016\n",
      "1.779197348711e+01 1.3228 17.3205 1.71196 14.4504 13 13 13 8 0.02 0.000287968\n",
      "1.367065497322e+01 1.66003 73.0434 1.79487 12.2572 13 13 13 8 0.02 7.67586e-05\n",
      "1.772169325355e+01 1.81957 73.0434 1.46573 12.2572 13 13 13 5 0.02 4.43422e-05\n",
David Hendriks's avatar
David Hendriks committed
      "2.021960493499e+01 1.82061 73.0434 1.39205 12.2572 13 13 13 5 0.02 3.38788e-05\n",
      "9.012246630357e+00 1.81529 73.0434 0 36.5717 13 8 15 15 0.02 7.67586e-05\n",
      "7.462779538274e+00 1.82255 73.0434 1.81499 36.5717 13 13 13 8 0.02 3.38788e-05\n",
      "1.030499912298e+01 1.80592 73.0434 1.81066 36.5717 13 13 13 8 0.02 4.43422e-05\n",
      "9.823059079115e+00 2.43711 73.0434 1.81689 60.8862 14 14 13 8 0.02 7.67586e-05\n",
      "7.394722435913e+00 1.79092 73.0434 1.79092 60.8862 13 8 13 8 0.02 4.43422e-05\n",
David Hendriks's avatar
David Hendriks committed
      "7.396288708628e+00 1.8216 73.0434 1.8216 60.8862 13 8 13 8 0.02 3.38788e-05\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def parse_function(self, output):\n",
    "    \"\"\"\n",
    "    Example parsing function\n",
    "    \"\"\"\n",
    "    \n",
    "    # extract info from the population instance\n",
    "\n",
    "    # Get some information from the\n",
    "    data_dir = self.custom_options[\"data_dir\"]\n",
    "    base_filename = self.custom_options[\"base_filename\"]\n",
    "\n",
    "    # Check directory, make if necessary\n",
    "    os.makedirs(data_dir, exist_ok=True)\n",
    "\n",
    "    #\n",
    "    seperator = \" \"\n",
    "\n",
    "    # Create filename\n",
    "    outfilename = os.path.join(data_dir, base_filename)\n",
    "\n",
    "    # The header columns matching this \n",
    "    parameters = [\n",
    "        \"time\", \n",
    "        \"mass_1\", \"zams_mass_1\", \"mass_2\", \"zams_mass_2\",\n",
    "        \"stellar_type_1\", \"prev_stellar_type_1\", \"stellar_type_2\", \"prev_stellar_type_2\", \n",
    "        \"metallicity\", \"probability\"\n",
    "    ]\n",
    "    \n",
    "    # Go over the output.\n",
    "    for el in output.splitlines():\n",
    "        headerline = el.split()[0]\n",
    "\n",
    "        # CHeck the header and act accordingly\n",