Track mergers and disruptions

[katiebreivik]

Function that selects out bcm and bpp arrays for merged and disrupted systems. See bpp.evol_type and here: https://cosmic-popsynth.github.io/examples/index.html#single-binary to get the key for the different evol_type (e.g. contact or disrupted).

This would ideally be set as a command line flag in the runFixedPop call (e.g. --merger_rates)

[scottcoughlin2014]

@katiebreivik I am trying to tackle this issue this morning. I think I have identified the impact of this feature change to be here for specifying whether or not to keep disrupted binaries

And for adding evol_type options something like

 74     parser.add_argument("--merger_rates", help="filter for evolutionary change "
 75                         "type: default is to keep all",
 76                         default=[1,2,3,4,5,6,7,8,9,10,11,12,13,14],
 77                         type=int, nargs='+')

I am looking for feedback on how this filtering should be structured, and in fact if we should dedicate a whole part of the code to some really heirachal and clearly laid out filtering.

Right now looks like the user experiences unknown hardcoded filtering of their bcm output

255         # Filter out any binaries that have has mass transfer from a WD onto a compact object
256         idx_MT = bpp.loc[(bpp.kstar_1.isin([10,11,12,13,14])) &(bpp.kstar_2.isin([10,11,12])) &
257                          (bpp.evol_type == 3.0)].bin_num
258         bcm = bcm.loc[~bcm.bin_num.isin(idx_MT)]
259 
260         # Select the state of the binary today and filter out 
261         # any disrupted binaries
262         bcm_filtered = bcm.loc[bcm.tphys > 1.0]
263         bcm_filtered = bcm_filtered.loc[bcm_filtered.sep > 0.0]

[katiebreivik]

@scottcoughlin2014 I think we should definitely have an end goal of hierarchical filters. All filtering right now is hard coded and removes a lot of flexibility we could have (e.g. tracking how things are disrupted/merged).

I think it would be helpful to also solicit input from @chasebk and @michaelzevin since they have requested for different filtering schemes (e.g. binary parameters just before SN explosions).

As a start though, since the user is already specifying the final kstars that they are interested in, I suspect that a good default is to filter for mergers or disruptions that occur for the star types listed as final kstars.

[scottcoughlin2014]

In fat, @katiebreivik I wonder if much like your filter flag

266         if args.LISA_convergence:

if we should have a whole utility function detected to sepecific string flags that represent some integer filtering in the bpp and bcm array. I fear that with all the obviously desireable custom filtering things will get too confusing if there is not a separate well documented utility function describing what say

bpp.filter(method=['lisa_convergence', 'mt_wd_on_mc', 'common envelope'])

[scottcoughlin2014]

Now that I wrote that our I wonder if a method on the bpp/bcm mtable is perferabelt o a method that would go like

from cosmic.utils import filter
bpp = filter(bpp, method = ['lisa_convergence', 'mt_wd_on_mc', 'common envelope'])

both seem good actually now that I write it out

[katiebreivik]

@scottcoughlin2014: excellent point. I really like the idea of having a separate filter module (method?). Is this the type of thing that they could then just specify in the inifile and be able to track exactly what filters are placed? If so, then I'm all for it since it will also allow for easier future upgrades and more filtering schemes.

[scottcoughlin2014]

Ya I think putting filtering methods in the inifile is definitely a good idea. I will work on structuring this ina initial PR and we can work from there

[scottcoughlin2014]

A niave attempt to demonstrate what this may look like logically

[filters]
mass_transfer_white_dwarf_to_co = False
lisa_convergence = True
select_binary_state = today
disrupted_binaries = False

so you have a descriptive phrase and in general it should have a boolean state yes I want those kind of things or no I do not. So right now we do not have disrupted binaries or mass_transfer_white_dwarf_to_co so you would sya False in the ini file

[katiebreivik]

I think this is maybe still confusing. I know this runs the risk of over complicating things, but I wonder if we split this into two structures:

[filters]
; True == retain binaries
; False == throw away binaries
mass_transfer_white_dwarf_to_co = False
select_final_state = True
disrupted_binaries = False
porb_above_1e4 = False

[convergence]
; perform convergence over selected region
lisa_convergence=True

The reason I think this might be helpful is that the LISA convergence should, in principle, not filter out sources. Rather, it should compute the convergence over a certain subset of data that fulfills a range. I'm open to suggestions though-

[scottcoughlin2014]

Ya that is totally fine.

[scottcoughlin2014]

63 Is a WIP PR to start the ball rolling on handling filtering bpp and bcm through the ini file. The convergence logic is not in. I am confused about

275             if args.LISA_convergence:
276                 bcm_match = bcm_save.loc[bcm_save.porb < 3.5]
277                 bcm_match_filtered = bcm_save_filtered.loc[bcm_save_filtered.porb < 3.5]
278             else:
279                 bcm_match = bcm_save
280                 bcm_match_filtered = bcm_save_filtered

and

253         if args.LISA_convergence:
254             bcm_save_filtered = bcm_save_filtered.loc[bcm_save_filtered.porb < 4]

why are these different?

[scottcoughlin2014]

@katiebreivik is it the case that if we were tracking disrupted binaries (either they merged before t=today or whatever) would we want to make sure we were tracking the convergence of the parameters of disrupted systems? would anyone want to retain both disrupted and those that still exists today? I ask this because I am realizing our True False methodology for discrupted binaries assumed you either want to track disrupted binaries or those alive today but not both.

[scottcoughlin2014]

Now that I think about it probably if the method is commented out then we simply do not use it. and keep both disrupted and not disrupted binaries.

Scotty

[katiebreivik]

63 Is a WIP PR to start the ball rolling on handling filtering bpp and bcm through the ini file. The convergence logic is not in. I am confused about

275             if args.LISA_convergence:
276                 bcm_match = bcm_save.loc[bcm_save.porb < 3.5]
277                 bcm_match_filtered = bcm_save_filtered.loc[bcm_save_filtered.porb < 3.5]
278             else:
279                 bcm_match = bcm_save
280                 bcm_match_filtered = bcm_save_filtered

and

253         if args.LISA_convergence:
254             bcm_save_filtered = bcm_save_filtered.loc[bcm_save_filtered.porb < 4]

why are these different?

Lines 253-254 filter out all systems with orbital periods greater than 1e4 seconds which you should do for a 'LISA' Galaxy since LISA won't observe them.

The other code sets the systems that the 'match' or convergence is computed for.

[katiebreivik]

Now that I think about it probably if the method is commented out then we simply do not use it. and keep both disrupted and not disrupted binaries.

Yeah, I had thought that if you do a False for retain disrupted binaries, then you don't retain them. This same logic then also goes for the 'alive today' binaries.

So if you do a False for disrupted and a False for alive today, you get nothing. If you do a True for both, you get both.

[scottcoughlin2014]

^^^ gotcha

[scottcoughlin2014]

@katiebreivik we have hit a place where #71 made it now easier to track mergers and disruptions in the bcm array but now the logic of setting the appropriate filters in the ini file is a little tricky. The way I see it for bin state we now have to smartly constrcut a way to ask for just one of merger distrupt still alive, a combo of two, and of course all of those states. Thoughts?

[scottcoughlin2014]

I think I may have a solution.

In the inifile

; 0 alive today, 1 merged, 2 disrupted
binary_state = [0,1,2]

then in the code it is super simple

        elif (meth == 'binary_state'):
            bcm = bcm.loc[bcm.bin_state.isin(use)]

[scottcoughlin2014]

@katiebreivik I am running into a weird initCthat leads to a bin_state=1 but a sytem that does not seem to have actually merged or the sep column is set incorrectly. The code below can reproduce the odd bcm array I am looking at. Thoughts?

from cosmic.sample.initialbinarytable import InitialBinaryTable
from cosmic.evolve import Evolve
import pandas as pd
import numpy as np
from cosmic import _evolvebin

BSEDict = {'neta': 0.5, 'bwind': 0.0, 'hewind': 1.0, 'alpha1': 1.0, 'lambdaf': 1.0, 'ceflag': 0, 
           'tflag': 1, 'ifflag': 0, 'wdflag': 0, 'bhflag': 1, 'nsflag': 3, 'mxns': 3.0, 'pts1': 0.05, 
           'pts2': 0.01, 'pts3': 0.02, 'sigma': 265.0, 'beta': -1.0, 'xi': 0.5, 'acc2': 1.5, 'epsnov': 0.001, 
           'eddfac': 1.0, 'gamma': -2.0, 'bconst': -3000, 'CK': -1000, 'merger': 0, 'windflag': 3, 
           'B_0': [0.0, 0.0], 'bacc': [0.0, 0.0], 'tacc': [0.0, 0.0],
           'bkick': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 'massc': [0.0, 0.0], 
           'opsin': [0.0, 0.0], 'epoch': [0.0, 0.0], 'tphys': 0.0, 'ppsn': 1}
single_binary = InitialBinaryTable.SingleBinary(m1=36.019611, m2=29.477283, porb=5693.803735, ecc=0.792748,
                                                tphysf=13700.0,
                                                kstar1=1.0, kstar2=1.0, metallicity=0.03)
bpp, bcm, initC = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)

initC['randomseed'] = 707971

bpp_columns = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2' , 'sep', 'ecc', 'RROL_1', 'RROL_2', 'evol_type']
bcm_columns = ['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lumin_1', 'rad_1', 'teff_1', 'massc_1',
'radc_1', 'menv_1', 'renv_1', 'epoch_1', 'ospin_1', 'deltam_1', 'RROL_1', 'kstar_2', 'mass0_2', 'mass_2',
'lumin_2', 'rad_2', 'teff_2', 'massc_2', 'radc_2', 'menv_2', 'renv_2', 'epoch_2', 'ospin_2', 'deltam_2', 'RROL_2',
'porb', 'sep', 'ecc', 'B_0_1', 'B_0_2', 'formation_1', 'formation_2', 'bin_state']
f = np.array(initC)
f =  f[0]
[tmp1, tmp2] = _evolvebin.evolv2(f[0], f[1], f[2], f[3], f[4], f[5], f[6], f[7], f[8], f[9],
                                        f[10], f[11], f[12], f[13], f[14], f[15], f[16], f[17], f[18], f[19],
                                        f[20], f[21], f[22], f[23], f[24], f[25], f[26], f[27], f[28], f[29],
                                        f[30], f[31], f[32], f[33], f[34], f[35], f[36])

bpp_tmp = tmp1[np.argwhere(tmp1[:,0]>0),:].squeeze(1)
bcm_tmp = tmp2[np.argwhere(tmp2[:,0]>0),:].squeeze(1)

bpp_tmp = pd.DataFrame(bpp_tmp, columns=bpp_columns, index=[int(f[37])] * len(bpp_tmp))
bpp_tmp['bin_num'] = int(f[37])

bcm_tmp = pd.DataFrame(bcm_tmp, columns=bcm_columns, index=[int(f[37])] * len(bcm_tmp))
bcm_tmp['bin_num'] = int(f[37])

[katiebreivik]

Ah- good find. I need to do some digging in evolv2 to find where this is happening. I think the problem is that evolv2 sets conditions for contact, then overrides the merger if the system somehow survives the common envelope. In this case, I don't have bin_state resetting to zero if this override happens. I'll look at this right now-

NOTE: this is now fixed in this PR: https://github.com/COSMIC-PopSynth/COSMIC/pull/73

[katiebreivik]

I think I may have a solution.

In the inifile

; 0 alive today, 1 merged, 2 disrupted
binary_state = [0,1,2]

then in the code it is super simple

        elif (meth == 'binary_state'):
            bcm = bcm.loc[bcm.bin_state.isin(use)]

I like this solution. We'll need to update it, per the PR I just submitted, to include bin_state=[0,1,2,3]

[scottcoughlin2014]

@katiebreivik I will get working on this. However, I feel like this bring me to another issue we should address shortly. #75

[scottcoughlin2014]

@katiebreivik copied here from Skype for posterity

2:52 PM Hey Katie, I am trying to do a run with the following settings along with the track-mergers branch to see if I can produce a merged BNS from runFixedPop

runFixedPop --final_kstar1 13 --final_kstar2 13 --inifile /projects/b1011/scoughlin/COSMIC/Emporium2/Params.ini --galaxy_component FIRE --metallicity 0.02 --convergence-params mass_1 mass_2 porb ecc --initial_samp multidim --Nstep 100000 --Niter 100000000 -n 20

I have gotten

The total mass sampled so far is: 142701409.77271318

and nothing so far. Does this fit your expectation? in the ini

[filters]
; True == retain binaries
; False == throw away binaries
mass_transfer_white_dwarf_to_co = False
select_final_state = True
; 0 alive today, 1 merged, 2 disrupted, 3 no remnant left after SN
binary_state = [1]
;porb_above_1e4 = False
;LISA_sources = False

[convergence]
; perform convergence over selected region
LISA_convergence=True

[scottcoughlin2014]

@katiebreivik What did we ever say about why the bcm array sometimes has only one entry in it which is both first state and last state for some reason see

(Pdb) tmp.loc[tmp.tphys == 1]
     tphys  kstar_1   mass0_1    mass_1   lumin_1     rad_1    teff_1  massc_1   ...           sep       ecc  B_0_1  B_0_2  formation_1  formation_2  bin_state  bin_num
135    1.0      1.0  3.518897  3.518897  2.166811  0.335699  4.137123      0.0   ...     19.662191  0.020919    0.0    0.0          0.0          0.0        0.0      135
198    1.0      1.0  4.260725  4.260725  2.474677  0.382867  4.190505      0.0   ...     13.938964  0.003453    0.0    0.0          0.0          0.0        0.0      198
385    1.0      1.0  3.518897  3.518897  2.166811  0.335699  4.137123      0.0   ...     19.662191  0.020919    0.0    0.0          0.0          0.0        0.0      385
448    1.0      1.0  4.260725  4.260725  2.474677  0.382867  4.190505      0.0   ...     13.938964  0.003453    0.0    0.0          0.0          0.0        0.0      448

More specifically after applying the final_state filter you go from

(Pdb) len(bcm)
996
(Pdb) len(filter_bpp_bcm(bcm, bpp, filters))
500

[scottcoughlin2014]

@katiebreivik For first example

from cosmic.sample.initialbinarytable import InitialBinaryTable
from cosmic.evolve import Evolve
import pandas as pd
import numpy as np
from cosmic import _evolvebin

BSEDict = {'neta': 0.5, 'bwind': 0.0, 'hewind': 1.0, 'alpha1': 1.0, 'lambdaf': 1.0, 'ceflag': 0, 
           'tflag': 1, 'ifflag': 0, 'wdflag': 0, 'bhflag': 1, 'nsflag': 3, 'mxns': 3.0, 'pts1': 0.05, 
           'pts2': 0.01, 'pts3': 0.02, 'sigma': 265.0, 'beta': -1.0, 'xi': 0.5, 'acc2': 1.5, 'epsnov': 0.001, 
           'eddfac': 1.0, 'gamma': -2.0, 'bconst': -3000, 'CK': -1000, 'merger': 0, 'windflag': 3, 
           'B_0': [0.0, 0.0], 'bacc': [0.0, 0.0], 'tacc': [0.0, 0.0],
           'bkick': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 'massc': [0.0, 0.0], 
           'opsin': [0.0, 0.0], 'epoch': [0.0, 0.0], 'tphys': 0.0, 'ppsn': 1}
single_binary = InitialBinaryTable.SingleBinary(m1=3.518899, m2=3.508906, porb=3.811531, ecc=0.020918,
                                                tphysf=5016.353066,
                                                kstar1=1.0, kstar2=1.0, metallicity=0.02)
bpp, bcm, initC = Evolve.evolve(initialbinarytable=single_binary, BSEDict=BSEDict)

initC['randomseed'] = 243848

bpp_columns = ['tphys', 'mass_1', 'mass_2', 'kstar_1', 'kstar_2' , 'sep', 'ecc', 'RROL_1', 'RROL_2', 'evol_type']
bcm_columns = ['tphys', 'kstar_1', 'mass0_1', 'mass_1', 'lumin_1', 'rad_1', 'teff_1', 'massc_1',
'radc_1', 'menv_1', 'renv_1', 'epoch_1', 'ospin_1', 'deltam_1', 'RROL_1', 'kstar_2', 'mass0_2', 'mass_2',
'lumin_2', 'rad_2', 'teff_2', 'massc_2', 'radc_2', 'menv_2', 'renv_2', 'epoch_2', 'ospin_2', 'deltam_2', 'RROL_2',
'porb', 'sep', 'ecc', 'B_0_1', 'B_0_2', 'formation_1', 'formation_2', 'bin_state']
f = np.array(initC)
f =  f[0]
[tmp1, tmp2] = _evolvebin.evolv2(f[0], f[1], f[2], f[3], f[4], f[5], f[6], f[7], f[8], f[9],
                                        f[10], f[11], f[12], f[13], f[14], f[15], f[16], f[17], f[18], f[19],
                                        f[20], f[21], f[22], f[23], f[24], f[25], f[26], f[27], f[28], f[29],
                                        f[30], f[31], f[32], f[33], f[34], f[35], f[36])

bpp_tmp = tmp1[np.argwhere(tmp1[:,0]>0),:].squeeze(1)
bcm_tmp = tmp2[np.argwhere(tmp2[:,0]>0),:].squeeze(1)

bpp_tmp = pd.DataFrame(bpp_tmp, columns=bpp_columns, index=[int(f[37])] * len(bpp_tmp))
bpp_tmp['bin_num'] = int(f[37])

bcm_tmp = pd.DataFrame(bcm_tmp, columns=bcm_columns, index=[int(f[37])] * len(bcm_tmp))
bcm_tmp['bin_num'] = int(f[37])

Output for

In [7]: bcm_tmp
Out[7]: 
   tphys  kstar_1   mass0_1    mass_1   lumin_1     rad_1    teff_1  massc_1   ...           sep       ecc  B_0_1  B_0_2  formation_1  formation_2  bin_state  bin_num
0    1.0      1.0  3.518896  3.518896  2.166811  0.335699  4.137123      0.0   ...     19.662193  0.020919    0.0    0.0          0.0          0.0        0.0        0

[1 rows x 38 columns]

In [8]: bpp_tmp
Out[8]: 
        tphys    mass_1    mass_2  kstar_1  kstar_2         sep       ecc    RROL_1    RROL_2  evol_type  bin_num
0  248.428925  3.518372  3.508392      2.0      1.0   20.148928  0.017226  0.659561  0.795973        2.0        0
0  248.810394  3.518325  3.508400      2.0      1.0   20.120745  0.016865  1.000270  0.791744        3.0        0
0  249.009445  3.510758  3.515946      2.0      1.0   20.054161  0.016164  1.240760  0.799469       14.0        0
0  249.764359  1.336033  5.690557      3.0      1.0   47.185654  0.000000  2.212173  0.173002        2.0        0
0  250.997009  0.551731  6.474408      4.0      1.0  205.458649  0.000000  1.702142  0.035024        2.0        0
0  251.051926  0.547529  6.478584      4.0      1.0  207.748657  0.000000  0.984973  0.034632        4.0        0
0  251.714752  0.544148  6.480565      7.0      1.0  207.643723  0.000000  0.002984  0.034952        2.0        0
0  276.313171  0.543800  6.462473      7.0      2.0  208.190201  0.000000  0.003274  0.057946        2.0        0
0  276.512970  0.543798  6.461960      7.0      3.0  208.637085  0.000000  0.003269  0.664089        2.0        0
0  276.528839  0.543798  6.461858      7.0      3.0  155.714661  0.000000  0.004380  1.001624        3.0        0
0  276.528839  0.543798  1.217554      7.0      7.0    1.065858  0.000000  0.004380  1.001624        7.0        0
0  276.528839  0.543798  1.217554      7.0      7.0    1.065858  0.000000  0.403579  0.493768        4.0        0
0  290.090546  0.544579  1.193009      7.0      8.0    1.034981  0.000000  0.434602  0.503179        2.0        0
0  291.067932  0.544946  1.185558      7.0      8.0    1.034092  0.000000  0.435963  1.000834        3.0        0
0  291.067932  0.544946  1.668191     15.0      8.0    0.059324  0.000000  0.435963  1.000834        7.0        0
0  291.513367  0.000000  1.639101     15.0      9.0    0.000000 -1.000000 -1.000000  0.000100        2.0        0
0  291.655060  0.000000  1.193009     15.0     15.0    0.000000  0.000000  0.000000 -2.000000        9.0        0

[scottcoughlin2014]

89 #63 Allows us to filter for only disrupted, alive today, or merged systems as well as filtering on type of merged system.