· 5 years ago · Oct 11, 2020, 08:52 PM
1# -*- coding: utf-8 -*-
2"""
3Created on Wed Aug 12 12:24:13 2020
4
5@author: norma
6"""
7import sqlite3
8import numpy as np
9import pandas as pd
10
11import populations
12
13import matplotlib.pyplot as plt
14from matplotlib.transforms import Affine2D
15from subprocess import run
16
17from uuid import uuid4
18
19class ulx:
20 def __init__(self):
21 self.is_ulx = None
22 self.is_transient = None
23 self.P_cycle_1_ulx = None
24 self.P_transient = [None, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
25 self.transient_cycle = None
26
27 @classmethod
28 def persistent_alive_system(cls):
29 ulx = cls()
30 ulx.P_transient = [None, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
31 ulx.P_cycle_1_ulx = 1.0
32 return ulx
33
34 @classmethod
35 def persistent_dead_system(cls):
36 ulx = cls()
37 ulx.P_transient = [None, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
38 ulx.P_cycle_1_ulx = 0.0
39 return ulx
40
41 @classmethod
42 def from_df_classifications_row(cls, Series, period):
43 ulx = cls()
44
45 ulx.P_cycle_1_ulx = Series.erass_1_ulx_prob
46 ulx.P_transient = [None, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
47 if period == 'P_wind':
48 ulx.P_transient[1] = Series.erass_2_P_wind_transient_prob
49 ulx.P_transient[2] = Series.erass_3_P_wind_transient_prob
50 ulx.P_transient[3] = Series.erass_4_P_wind_transient_prob
51 ulx.P_transient[4] = Series.erass_5_P_wind_transient_prob
52 ulx.P_transient[5] = Series.erass_6_P_wind_transient_prob
53 ulx.P_transient[6] = Series.erass_7_P_wind_transient_prob
54 ulx.P_transient[7] = Series.erass_8_P_wind_transient_prob
55 elif period == 'P_sup':
56 ulx.P_transient[1] = Series.erass_2_P_sup_transient_prob
57 ulx.P_transient[2] = Series.erass_3_P_sup_transient_prob
58 ulx.P_transient[3] = Series.erass_4_P_sup_transient_prob
59 ulx.P_transient[4] = Series.erass_5_P_sup_transient_prob
60 ulx.P_transient[5] = Series.erass_6_P_sup_transient_prob
61 ulx.P_transient[6] = Series.erass_7_P_sup_transient_prob
62 ulx.P_transient[7] = Series.erass_8_P_sup_transient_prob
63 return ulx
64
65
66 def is_ulx_first_cycle(self):
67 self.is_ulx = np.random.random() < self.P_cycle_1_ulx
68 return self.is_ulx
69
70 def observe(self, cycle):
71 if cycle == 0:
72 return [self.is_ulx_first_cycle(), self.is_transient]
73 else:
74 self.is_transient = np.random.random() < self.P_transient[cycle]
75 if self.is_transient:
76 self.is_ulx = not self.is_ulx
77 if self.transient_cycle == None:
78 self.transient_cycle = cycle
79 return [self.is_ulx, self.is_transient]
80
81
82class Sim_eRASS:
83 def __init__(self, systems):
84 self.systems = systems
85 self.size = len(self.systems)
86
87 self.N_new_systems = np.array([0,0,0,0,0,0,0,0])
88 self.N_old_system_become_transient = np.array([0,0,0,0,0,0,0,0])
89
90 def calc_secondary_quantities(self):
91 self.N_delta_obs_ulx = self.N_new_systems - self.N_old_system_become_transient
92 self.N_observed_ulxs = np.cumsum(self.N_delta_obs_ulx)
93 self.N_transients = self.N_new_systems[1:] + self.N_old_system_become_transient[1:]
94 self.N_transients = np.insert(self.N_transients, 0, 0)
95 self.N_transients_cum = np.cumsum(self.N_transients)
96 self.N_total_systems = np.cumsum(self.N_new_systems)
97 self.N_persistent_ulx_systems = self.N_new_systems[0] - np.cumsum(self.N_old_system_become_transient)
98
99 @classmethod
100 def from_run_id(cls, db, run_id, period):
101 sql1 = f"""SELECT * FROM ERASS_MC_SAMPLED_SYSTEMS
102 WHERE run_id='{run_id}'"""
103 sql2 = f"""SELECT * FROM CLASSIFICATIONS
104 WHERE run_id='{run_id}'"""
105 sql3 = f"""SELECT * FROM TRANSIENT
106 WHERE run_id='{run_id}'"""
107
108
109 conn = sqlite3.connect(db)
110 df_sampled_systems = pd.read_sql_query(sql1, conn)
111 df_classifications = pd.read_sql_query(sql2, conn)
112 df_transient = pd.read_sql_query(sql3, conn)
113 conn.close()
114
115 class_count = df_classifications['lc_classification'].value_counts()
116
117 size = len(df_sampled_systems)
118 N_systems_not_simulated = size - len(df_classifications) # These systems are persistent
119 try:
120 N_alive_systems = class_count[1]
121 except KeyError:
122 N_alive_systems = 0
123 try:
124 N_dead_systems = class_count[0]
125 except KeyError:
126 N_dead_systems = 0
127 N_transient_not_sampled = N_alive_systems - len(df_transient) # P_wind & P_sup > 4 year systems
128 N_persistent_systems = N_systems_not_simulated + N_alive_systems + N_transient_not_sampled
129
130
131 persistent_systems = [ulx.persistent_alive_system() for i in range(N_persistent_systems)]
132 dead_systems = [ulx.persistent_dead_system() for i in range(N_dead_systems)]
133 transient_systems = [ulx.from_df_classifications_row(row, period) for i, row in df_transient.iterrows()]
134
135 systems = persistent_systems + dead_systems + transient_systems
136 Sim_eRASS = cls(systems)
137 return Sim_eRASS
138
139
140 def run(self):
141 for cycle in range(8):
142 new_systems = 0
143 old_system_become_transient = 0
144
145 for s in self.systems:
146 is_ulx, is_transient = s.observe(cycle)
147 if is_ulx and is_transient == None:
148 new_systems += 1
149 continue
150 elif is_ulx and is_transient and s.transient_cycle == cycle:
151 new_systems += 1
152 continue
153 elif not is_ulx and is_transient and s.transient_cycle == cycle:
154 old_system_become_transient += 1
155 continue
156
157 self.N_new_systems[cycle] = new_systems
158 if cycle>0:
159 self.N_old_system_become_transient[cycle] = old_system_become_transient
160 self.calc_secondary_quantities()
161
162 def collect_results(self):
163 res = pd.DataFrame()
164 res['erass_cycle'] = np.arange(1, 9)
165 res['N_new_systems'] = self.N_new_systems
166 res['N_old_system_become_transient'] = self.N_old_system_become_transient
167 res['N_observed_ulxs'] = self.N_observed_ulxs
168 res['N_delta_obs_ulx'] = self.N_delta_obs_ulx
169 res['N_transients'] = self.N_transients
170 res['N_transients_cum'] = self.N_transients_cum
171 res['N_total_systems'] = self.N_total_systems
172 res['N_persistent_ulx_systems'] = self.N_persistent_ulx_systems
173 self.res = res
174 return res
175
176 def sim_write_to_sql(self):
177 pass
178
179
180class ResultsProcessor:
181 def __init__(self):
182 self.plot_set_latex_font()
183
184 # matplotlib colors
185 self.color_dead = 'C3'
186 self.color_trans = 'C1'
187 self.color_alive = 'C2'
188
189 self.linestyle_dead = 'dotted'
190 self.linestyle_trans = '--'
191 self.linestyle_alive = '-'
192
193
194 #Code for these figures are in ../reports/investigations.ipynb
195 self.PERCENT_ALIVE_EARNSHAW = 0.8271604938271605 * 100
196 self.PERCENT_ALIVE_EARNSHAW_ERROR = 0.12256472421344072 * 100
197
198 self.PERCENT_ALIVE_EARNSHAW_UPPER = self.PERCENT_ALIVE_EARNSHAW + self.PERCENT_ALIVE_EARNSHAW_ERROR
199 self.PERCENT_ALIVE_EARNSHAW_LOWER = self.PERCENT_ALIVE_EARNSHAW - self.PERCENT_ALIVE_EARNSHAW_ERROR
200
201 self.PERCENT_TRANS_EARNSHAW = 0.1728395061728395 * 100
202 self.PERCENT_TRANS_EARNSHAW_ERROR = 0.03744750536124969 * 100
203 self.PERCENT_TRANS_EARNSHAW_UPPER = self.PERCENT_TRANS_EARNSHAW + self.PERCENT_TRANS_EARNSHAW_ERROR
204 self.PERCENT_TRANS_EARNSHAW_LOWER = self.PERCENT_TRANS_EARNSHAW - self.PERCENT_TRANS_EARNSHAW_ERROR
205
206 def set_active_db(self, path):
207 self.db = path
208
209 def set_parent_population(self, Population):
210 self.pop = Population
211
212 def table_create_erass_mc_info(self):
213 conn = sqlite3.connect(self.db)
214 sql = """CREATE TABLE IF NOT EXISTS ERASS_MC_INFO(
215 run_id CHAR,
216 size INT,
217 bh_ratio REAL,
218 dincl_cutoff INT,
219 Z CHAR,
220 erass_system_period_cutoff INT,
221 erass_lmxrb_duty_cycle REAL);"""
222 conn.execute(sql)
223 conn.close()
224
225 def table_create_erass_mc_results(self):
226 conn = sqlite3.connect(self.db)
227 sql = """CREATE TABLE IF NOT EXISTS ERASS_MC_RESULTS(
228 erass_cycle INT,
229 N_new_systems INT,
230 N_old_system_become_transient INT,
231 N_observed_ulxs INT,
232 N_delta_obs_ulx INT,
233 N_transients INT,
234 N_transients_cum INT,
235 N_total_systems INT,
236 N_persistent_ulx_systems INT,
237 period CHAR,
238 run_id CHAR);"""
239 conn.execute(sql)
240 conn.close()
241
242 def table_create_erass_mc_sampled_systems(self):
243 conn = sqlite3.connect(self.db)
244 sql = """CREATE TABLE IF NOT EXISTS ERASS_MC_SAMPLED_SYSTEMS(
245 row_id INT,
246 theta_half_deg REAL,
247 Lx1 REAL,
248 P_wind_days REAL,
249 P_sup_days REAL,
250 lmxrb INT,
251 inclination INT,
252 dincl INT,
253 idum_run INT,
254 iidd_old INT,
255 run_id CHAR);"""
256 conn.execute(sql)
257 sql = """CREATE INDEX IF NOT EXISTS idx_sampled_systems_run_id
258 ON ERASS_MC_SAMPLED_SYSTEMS (run_id);"""
259 conn.execute(sql)
260 conn.close()
261
262 def table_load(self, table_name):
263 conn = sqlite3.connect(self.db)
264 df = pd.read_sql_query(f"SELECT * from {table_name}", conn)
265 conn.close()
266 return df
267
268 def table_load_transient(self):
269 self.df_transient = self.table_load('TRANSIENT')
270
271 def table_load_classifications(self):
272 self.df_classifications = self.table_load('CLASSIFICATIONS')
273
274 def table_load_erass_mc_info(self):
275 self.df_erass_mc_info = self.table_load('ERASS_MC_INFO')
276
277 def table_load_erass_mc_results(self):
278 self.df_erass_mc_results = self.table_load('ERASS_MC_RESULTS')
279
280 def table_load_erass_mc_sampled_systems(self):
281 self.df_erass_mc_sampled_systems = self.table_load('ERASS_MC_SAMPLED_SYSTEMS')
282
283 def table_load_all(self, load_sampled=False):
284 self.table_load_transient()
285 self.table_load_classifications()
286 self.table_load_erass_mc_info()
287 self.table_load_erass_mc_results()
288 if load_sampled == True:
289 self.table_load_erass_mc_sampled_systems()
290
291
292 def table_erass_mc_results_map_info(self):
293 try:
294 self.df_erass_mc_info
295 except:
296 self.table_load_erass_mc_info()
297 try:
298 self.df_erass_mc_results
299 except:
300 self.table_load_erass_mc_results()
301
302 info = self.df_erass_mc_info.set_index('run_id')
303 self.df_erass_mc_results['bh_ratio'] = self.df_erass_mc_results['run_id'].map(info['bh_ratio'])
304 self.df_erass_mc_results['size'] = self.df_erass_mc_results['run_id'].map(info['size'])
305 self.df_erass_mc_results['dincl_cutoff'] = self.df_erass_mc_results['run_id'].map(info['dincl_cutoff'])
306 self.df_erass_mc_results['Z'] = self.df_erass_mc_results['run_id'].map(info['Z'])
307 self.df_erass_mc_results['erass_system_period_cutoff'] = self.df_erass_mc_results['run_id'].map(info['erass_system_period_cutoff'])
308
309
310 def table_classifications_map_systems(self):
311 try:
312 self.df_classifications
313 except:
314 self.table_load_classifications()
315
316 self.df_classifications['is_bh'] = self.df_classifications['system_row_id'].map(self.pop.df['is_bh'])
317 self.df_classifications['P_wind_days'] = self.df_classifications['system_row_id'].map(self.pop.df['P_wind_days'])
318 self.df_classifications['a*'] = self.df_classifications['system_row_id'].map(self.pop.df['a*'])
319 self.df_classifications['Z'] = self.df_classifications['system_row_id'].map(self.pop.df['Z'])
320 self.df_classifications['lmxrb'] = self.df_classifications['system_row_id'].map(self.pop.df['lmxrb'])
321
322
323 def table_classifications_map_info(self):
324 info = self.df_erass_mc_info.set_index('run_id')
325 self.df_classifications['bh_ratio'] = self.df_erass_mc_results['run_id'].map(info['bh_ratio'])
326 self.df_classifications['size'] = self.df_erass_mc_results['run_id'].map(info['size'])
327 self.df_classifications['dincl_cutoff'] = self.df_erass_mc_results['run_id'].map(info['dincl_cutoff'])
328 self.df_classifications['Z'] = self.df_erass_mc_results['run_id'].map(info['Z'])
329
330
331 def table_classifications_pivot(self, margins=True, split_Z=True):
332 try:
333 self.df_classifications['Z']
334 except KeyError:
335 self.table_classifications_map_systems()
336 if split_Z==True:
337 piv = pd.pivot_table(self.df_classifications,
338 columns=['is_bh'],
339 index=['Z', 'lc_classification'],
340 aggfunc='count',
341 margins=margins,
342 margins_name='total').run_id
343
344
345 if split_Z==False:
346 piv = pd.pivot_table(self.df_classifications,
347 columns=['is_bh'],
348 index=['lc_classification'],
349 aggfunc='count',
350 margins=margins,
351 margins_name='total').run_id
352
353 if margins==True:
354 piv['%_NS'] = piv[0]/piv['total']*100
355 piv['%_BH'] = piv[1]/piv['total']*100
356
357 self.df_pivot_lc_classifcations = piv
358 return self.df_pivot_lc_classifcations
359
360 def table_classifications_split_by_metallicity(self):
361 try:
362 self.df_classifications['Z']
363 except:
364 self.table_classifications_map_systems()
365
366 self.df_c_02 = self.df_classifications[self.df_classifications['Z']==0.02]
367 self.df_c_002 = self.df_classifications[self.df_classifications['Z']==0.002]
368 self.df_c_0002 = self.df_classifications[self.df_classifications['Z']==0.0002]
369
370 def table_classifications_calc_intermediate(self):
371 self.df_d = self.df_classifications[self.df_classifications['lc_classification']==0]
372 self.df_t = self.df_classifications[self.df_classifications['lc_classification']==1]
373 self.df_a = self.df_classifications[self.df_classifications['lc_classification']==2]
374
375 systems_per_dincl_bin = self.df_classifications['system_dincl'].value_counts().unique()[0]
376 systems_per_i_bin = self.df_classifications['system_inclination'].value_counts().unique()[0]
377
378 # Number of systems for each dincl
379 self.a_dincl_N = self.df_a['system_dincl'].value_counts().sort_index()
380 self.t_dincl_N = self.df_t['system_dincl'].value_counts().sort_index()
381 self.d_dincl_N = self.df_d['system_dincl'].value_counts().sort_index()
382
383 # Number of systems for each inclination
384 self.a_i_N = self.df_a['system_inclination'].value_counts().sort_index()
385 self.t_i_N = self.df_t['system_inclination'].value_counts().sort_index()
386 self.d_i_N = self.df_d['system_inclination'].value_counts().sort_index()
387
388 # Percentages
389 self.a_dincl_percent = (self.a_dincl_N/systems_per_dincl_bin * 100)
390 self.t_dincl_percent = (self.t_dincl_N/systems_per_dincl_bin * 100)
391 self.d_dincl_percent = (self.d_dincl_N/systems_per_dincl_bin * 100)
392
393 self.a_i_percent = (self.a_i_N/systems_per_i_bin * 100).sort_index()
394 self.t_i_percent = (self.t_i_N/systems_per_i_bin * 100).sort_index()
395 self.d_i_percent = (self.d_i_N/systems_per_i_bin * 100).sort_index()
396
397
398 def MC_bh_ratio_classifications_sampler(self, N=10000, size=500, dincl_cutoff=46, Z='all'):
399 def create_classification_dict():
400 df_classifications_reindexed = self.df_classifications.set_index(['system_row_id', 'system_dincl', 'system_inclination'])
401 class_dict = dict(df_classifications_reindexed['lc_classification'])
402 return class_dict
403
404 try:
405 self.class_dict
406 except:
407 self.class_dict = create_classification_dict()
408
409 if Z!='all':
410 self.pop.df_ulx[self.pop.df_ulx['Z'] == Z]
411
412 res = []
413
414 for i in range(N):
415 print(i)
416 for bh_ratio in np.arange(0, 1.05, 0.05):
417
418 selected_systems = self.pop.sample_ulxs(bh_ratio, size=size)
419 selected_dincls = np.random.randint(0, dincl_cutoff, size=size)
420 selected_inclinations = np.random.randint(0, 91, size=size)
421 selected_keys = list(zip(selected_systems, selected_dincls, selected_inclinations))
422
423 res_classications = [self.class_dict.get(key) for key in selected_keys]
424
425 # None systems correspond to opening angles > 45 and are considered alive
426 N_alive = res_classications.count(2) + res_classications.count(None)
427 N_transient = res_classications.count(1)
428 N_dead = res_classications.count(0)
429
430 res.append([bh_ratio, N_alive, N_transient, N_dead, i])
431
432 df_res = pd.DataFrame(res, columns=['bh_ratio', 'alive', 'transient', 'dead', 'iter'])
433 self.df_bh_ratio = df_res
434
435 #savepath = '../data/interim/bh_ns_sampling_results/'
436 #filename = f'Z = {Z}.csv'
437
438 #df_res.to_csv(savepath+filename, index=False)
439 # df_res.to_csv(savepath+filename, mode='a', header=False, index=False)
440
441
442 def MC_ERASS_simulation_run(self,
443 size=500,
444 dincl_cutoff=46,
445 Z='all',
446 bh_ratio=0.5,
447 erass_system_period_cutoff=1460,
448 erass_lmxrb_duty_cycle=0.1):
449 """
450 eRASS Monte-Carlo Simulation.
451
452 Parameters
453 ----------
454 size : int, optional
455 Sampled population size.
456 The default is 500.
457 dincl_cutoff : int, optional
458 Maximum precessional angle.
459 The default is 46.
460 Z : float
461 Metallicity to use, one of 'all', 0.02, 0.002, 0.0002.
462 The default is 'all'.
463 bh_ratio : float
464 Ratio of black holes to neutron stars.
465 The default is 0.5.
466 erass_system_period_cutoff : int, optional
467 Maximum period in days after which systems will be considered as persistent.
468 The default is 1460.
469 erass_lmxrb_duty_cycle : float, optional
470 Duty cycle to use for LMXRB systems.
471 The default is 0.1.
472
473 Returns
474 -------
475 None.
476
477 """
478 print('Running erass MC...')
479
480 # Create SQLite tables
481 self.table_create_erass_mc_results()
482 self.table_create_erass_mc_info()
483 self.table_create_erass_mc_sampled_systems()
484
485 # Get columns needed for simulation
486 sim_cols = ['theta_half_deg', 'Lx1', 'P_wind_days', 'P_sup_days', 'idum_run', 'iidd_old', 'lmxrb']
487
488 # Delete large dataframes to save memory
489 # del(self.pop.df)
490
491 # Check if we have the correct Z population
492 if self.pop.df_ulx_Z_subset != float(Z):
493 self.pop.filter_df_ulx_by_Z(Z)
494 self.pop.calc_bh_ns_ulx_sub_populations()
495
496 # Filter population df to only have columns we are interested in to save memory
497 self.pop.df_ulx = self.pop.df_ulx[sim_cols]
498
499 # Create run info dataframe for storing id, size etc..
500 df_info = pd.DataFrame()
501
502 # Random
503 run_id = str(uuid4())
504
505 df_info = pd.DataFrame()
506 df_info['run_id'] = [run_id]
507 df_info['size'] = [size]
508 df_info['bh_ratio'] = [bh_ratio]
509 df_info['dincl_cutoff'] = [dincl_cutoff]
510 df_info['Z'] = [Z]
511 df_info['erass_system_period_cutoff'] = [erass_system_period_cutoff]
512 df_info['erass_lmxrb_duty_cycle'] = [erass_lmxrb_duty_cycle]
513
514 sampled_indexs = self.pop.sample_ulxs(bh_ratio, size=size)
515 selected_inclinations = np.random.randint(0, 91, size=size)
516 selected_dincls = np.random.randint(0, dincl_cutoff, size=size)
517
518 df_sampled = self.pop.df_ulx.loc[sampled_indexs]
519 df_sampled['Lx1'] = df_sampled['Lx1'] / 1e39
520 df_sampled['inclination'] = selected_inclinations
521 df_sampled['dincl'] = selected_dincls
522 df_sampled['row_id'] = sampled_indexs
523 df_sampled['run_id'] = run_id
524
525 conn = sqlite3.connect(self.db)
526 df_sampled.to_sql('ERASS_MC_SAMPLED_SYSTEMS', conn, if_exists='append', index=False)
527 df_info.to_sql('ERASS_MC_INFO', conn, if_exists='append', index=False)
528 conn.close()
529
530 run(["ulxlc/ulxlc_code_v0.1/a.out", str(run_id), str(size), str(erass_system_period_cutoff), str(erass_lmxrb_duty_cycle)])
531
532 for period in ['P_wind', 'P_sup']:
533 sim = Sim_eRASS.from_run_id(self.db, run_id, period)
534 sim.run()
535 res = sim.collect_results()
536 res['period'] = period
537 res['run_id'] = run_id
538 conn = sqlite3.connect(self.db)
539 res.to_sql('ERASS_MC_RESULTS', conn, if_exists='append', index=False)
540 conn.close()
541
542
543 def MC_calc_N_persistent(self, size=500):
544 """Calculate the number of systems that had half opening angles > 45
545 and thus were not put through the lightcurve simulation"""
546 self.df_N_persistent = size - self.df_classifications['run_id'].value_counts()
547
548
549 def MC_get_run_ids(self, group_precession_cuttoff=True):
550 try:
551 self.df_erass_mc_info
552 except AttributeError:
553 self.table_load_erass_mc_info()
554
555 info = self.df_erass_mc_info.copy()
556 info = info.set_index('run_id')
557 if group_precession_cuttoff:
558 self.dict_MC_run_ids = info.groupby(by=['size', 'bh_ratio', 'dincl_cutoff', 'Z']).groups
559 else:
560 self.dict_MC_run_ids = info.groupby(by=['size', 'bh_ratio', 'dincl_cutoff', 'Z', 'erass_system_period_cutoff']).groups
561
562 def MC_get_run_counts(self):
563 count_dict = {}
564 for k, v in self.dict_MC_run_ids.items():
565 count_dict[k] = len(v)
566 return count_dict
567
568 def MC_get_classifications(self, key):
569 """
570 Parameters
571 ----------
572 key : tuple from key in self.dict_MC_run_ids
573 """
574 try:
575 self.dict_MC_run_ids
576 except:
577 self.MC_get_run_ids()
578
579 ids = self.dict_MC_run_ids[key]
580 df_classifications = self.df_classifications.set_index('run_id')
581 # See https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#deprecate-loc-reindex-listlike
582 df_classifications = df_classifications.loc[df_classifications.index.intersection(ids).unique()]
583 return df_classifications
584
585 def MC_get_erass_mc_results(self, key):
586 """
587 Parameters
588 ----------
589 key : tuple from key in self.dict_MC_run_ids
590 """
591 try:
592 self.dict_MC_run_ids
593 except:
594 self.MC_get_run_ids(group_precession_cuttoff=False)
595
596 ids = self.dict_MC_run_ids[key]
597 df_erass_mc_results = self.df_erass_mc_results.set_index('run_id')
598 df_erass_mc_results = df_erass_mc_results.loc[df_erass_mc_results.index.intersection(ids).unique()]
599 return df_erass_mc_results
600
601 def MC_ADT_result_stats(self):
602 try:
603 self.df_classifications['Z']
604 except:
605 self.table_classifications_map_info()
606 self.df_classifications_stats = self.df_classifications.groupby(['dincl_cutoff', 'Z', 'bh_ratio']).agg(['min', 'mean', 'max', 'std', 'count'])
607
608 def MC_ERASS_result_stats(self):
609 try:
610 self.df_erass_mc_results['dincl_cutoff']
611 except:
612 self.table_erass_mc_results_map_info()
613
614 self.df_erass_mc_stats = self.df_erass_mc_results.groupby(['dincl_cutoff', 'Z', 'bh_ratio', 'period', 'erass_cycle']).agg(['min', 'mean', 'max', 'std', 'count'])
615
616
617 def MC_calc_classification_counts(self, key):
618 try:
619 self.df_N_persistent
620 except AttributeError:
621 self.MC_calc_N_persistent()
622
623 if len(key) != 4:
624 print('key length !=4, probably seperated by precessional angle cut.')
625
626 df_c = self.MC_get_classifications(key)
627
628 df_d = df_c[df_c['lc_classification']==0]
629 df_t = df_c[df_c['lc_classification']==1]
630 df_a = df_c[df_c['lc_classification']==2]
631
632 df_N_d = df_d.index.value_counts().rename('N_dead')
633 df_N_t = df_t.index.value_counts().rename('N_trans')
634 df_N_a = df_a.index.value_counts().rename('N_alive')
635 df_N_p = self.df_N_persistent.rename('N_persistent')
636
637 df_classification_counts = pd.concat([df_N_d, df_N_t, df_N_a, df_N_p], axis=1).dropna()
638
639 # Add in those systems with opening angles > 45 that weren't simulated
640 df_classification_counts['N_alive_persistent'] = df_classification_counts['N_alive'] + df_classification_counts['N_persistent']
641
642 df_classification_counts['frac_alive_visible'] = 100 * df_classification_counts['N_alive_persistent'] / (df_classification_counts['N_alive_persistent'] + df_classification_counts['N_trans'])
643 df_classification_counts['frac_trans_visible'] = 100 * df_classification_counts['N_trans'] / (df_classification_counts['N_alive_persistent'] + df_classification_counts['N_trans'])
644 return df_classification_counts
645
646
647 def MC_calc_all_classifications_count_stats(self):
648 res = pd.DataFrame()
649 for k in self.dict_MC_run_ids.keys():
650 print(k)
651 df = self.MC_calc_classification_counts(k)
652 df['size'] = k[0]
653 df['bh_ratio'] = k[1]
654 df['dincl_cutoff'] = k[2]
655 df['Z'] = k[3]
656 res = res.append(df)
657 gb = res.groupby(['Z', 'dincl_cutoff', 'bh_ratio']).agg(['min', 'mean', 'max', 'std', 'count'])
658 gb = gb.drop('size', axis=1)
659 self.df_classifications_counts_stats = gb
660
661
662 def plot_set_latex_font(self):
663 import matplotlib
664 matplotlib.rcParams['mathtext.fontset'] = 'stix'
665 matplotlib.rcParams['font.family'] = 'STIXGeneral'
666
667 def plot_classifications_hist(self, Z, dincl_cut, frac_visible=False, save=False):
668
669 fig, ax = plt.subplots(5,1, sharex=True, sharey=True, figsize=(6, 4.5))
670 ax[0].set_title(fr'Z = {Z} | $\Delta i_{{max}} =$ {dincl_cut}$^{{\circ}}$')
671
672
673 for i, bh in enumerate([0.0, 0.25, 0.5, 0.75, 1.0]):
674 print(i, bh)
675 key = (500, bh, dincl_cut, Z)
676 df_classification_counts = self.MC_calc_classification_counts(key)
677 if frac_visible==False:
678 df_classification_counts['N_dead'].hist(bins=np.arange(0,500,5), label='Dead', alpha=0.8, edgecolor='black', linestyle=self.linestyle_dead, histtype='step', ax=ax[i], grid=False)
679 df_classification_counts['N_trans'].hist(bins=np.arange(0,500,5), label='Transient', alpha=0.8, edgecolor='black', linestyle=self.linestyle_trans, histtype='step', ax=ax[i], grid=False)
680 df_classification_counts['N_alive_persistent'].hist(bins=np.arange(0,500,5), label='Alive', alpha=0.8, edgecolor='black', linestyle=self.linestyle_alive, histtype='step', ax=ax[i], grid=False)
681 ax[-1].set_xlabel('Number of observed ULXs')
682 else:
683 ax[i].axvspan(self.PERCENT_ALIVE_EARNSHAW_LOWER, self.PERCENT_ALIVE_EARNSHAW_UPPER, alpha=0.3, color='grey')
684 ax[i].axvspan(self.PERCENT_TRANS_EARNSHAW_LOWER, self.PERCENT_TRANS_EARNSHAW_UPPER, alpha=0.7, color='grey')
685
686 df_classification_counts['frac_alive_visible'].hist(bins=np.arange(0,100,1), label='% Alive', alpha=0.8, edgecolor='black', linestyle=self.linestyle_alive, histtype='step', ax=ax[i], grid=False)
687 df_classification_counts['frac_trans_visible'].hist(bins=np.arange(0,100,1), label='% Transient', alpha=0.8, edgecolor='black', linestyle=self.linestyle_trans, histtype='step', ax=ax[i], grid=False)
688 ax[-1].set_xlabel('% of observed ULXs')
689
690 ax[i].text(x=350, y=ax[i].get_ylim()[1]*(1/2), s=r'$\%_{{BH}} = $ {}'.format((bh)*100), fontsize=9)
691 ax[i].legend(fontsize=7, loc='right')
692
693 if frac_visible:
694 ax[0].text(x=(self.PERCENT_ALIVE_EARNSHAW_LOWER)+0.25, y=ax[0].get_ylim()[1]+5, s='Alive $1 \sigma$ XMM', fontsize=7)
695 ax[0].text(x=(self.PERCENT_TRANS_EARNSHAW_LOWER)+0.25, y=ax[0].get_ylim()[1]+5, s='Transient $1 \sigma$ XMM', fontsize=7)
696
697
698
699 plt.tight_layout()
700 plt.subplots_adjust(hspace=0.0)
701
702 if save:
703 if frac_visible==False:
704 plt.savefig(f'../reports/figures/ADT_{Z}_{dincl_cut}.png', dpi=500)
705 plt.savefig(f'../reports/figures/ADT_{Z}_{dincl_cut}.eps')
706 plt.savefig(f'../reports/figures/ADT_{Z}_{dincl_cut}.pdf')
707 else:
708 plt.savefig(f'../reports/figures/ADT_frac_{Z}_{dincl_cut}.png', dpi=500)
709 plt.savefig(f'../reports/figures/ADT_frac_{Z}_{dincl_cut}.eps')
710 plt.savefig(f'../reports/figures/ADT_frac_{Z}_{dincl_cut}.pdf')
711
712
713
714 def plot_classifications_dincl_i(self, percent=False, savefig=False):
715 try:
716 self.df_d
717 except:
718 self.table_classifications_calc_intermediate()
719
720 fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(8,3))
721 plt.tight_layout()
722 ax[0].set_xlim(0,max(self.df_d['system_dincl']))
723 ax[0].set_xlabel(r'Precessional angle $\Delta i$')
724
725 ax[0].minorticks_on()
726
727
728 ax[1].set_xlim(0,max(self.df_d['system_inclination']))
729 ax[1].set_xlabel(r'Inclination $i$')
730 ax[1].minorticks_on()
731
732
733
734
735 if percent == False:
736 ax[0].set_ylabel(r'Number of systems')
737 ax[1].set_ylabel(r'Number of systems')
738
739 # dincl vs Number systems
740 self.a_dincl_N.plot(label='Alive', linestyle=self.linestyle_alive, color='black', linewidth=0.8, ax=ax[0])
741 self.t_dincl_N.plot(label='Transient', linestyle=self.linestyle_trans, color='black', linewidth=0.8, ax=ax[0])
742 self.d_dincl_N.plot(label='Dead', linestyle=self.linestyle_dead, color='black', linewidth=0.8, ax=ax[0])
743
744 # inclination vs Number systems
745 self.a_i_N.plot(label='Alive', linestyle=self.linestyle_alive, color='black', linewidth=0.8, ax=ax[1])
746 self.t_i_N.plot(label='Transient', linestyle=self.linestyle_trans, color='black', linewidth=0.8, ax=ax[1])
747 self.d_i_N.sort_index().plot(label='Dead', linestyle=self.linestyle_dead, color='black', linewidth=0.8, ax=ax[1])
748
749 ax[1].set_ylabel(r'Number of systems')
750
751 if savefig:
752 plt.savefig('../reports/figures/dincl_i_classifications.png', dpi=1000)
753 plt.savefig('../reports/figures/dincl_i_classifications.eps')
754
755 else:
756 ax[0].set_ylabel(r'% of systems')
757 ax[1].set_ylabel(r'% of systems')
758 ax[0].set_ylim(0,100)
759 ax[1].set_ylim(0,100)
760
761 # dincl vs percentages
762 self.a_dincl_percent.plot(label='Alive', linestyle=self.linestyle_alive, color='black', linewidth=0.8, ax=ax[0])
763 self.t_dincl_percent.plot(label='Transient', linestyle=self.linestyle_trans, color='black', linewidth=0.8, ax=ax[0])
764 self.d_dincl_percent.plot(label='Dead', linestyle=self.linestyle_dead, color='black', linewidth=0.8, ax=ax[0])
765
766 # inclination vs percentages
767 self.a_i_percent.plot(label='Alive', linestyle='-', color='black', linewidth=0.8, ax=ax[1])
768 self.t_i_percent.plot(label='Transient', linestyle='--', color='black', linewidth=0.8, ax=ax[1])
769 self.d_i_percent.plot(label='Dead', linestyle='dotted', color='black', linewidth=0.8, ax=ax[1])
770
771
772 if savefig:
773 plt.savefig('../reports/figures/dincl_i_classifications_percent.png', dpi=1000)
774 plt.savefig('../reports/figures/dincl_i_classifications_percent.eps')
775 ax[0].legend()
776 ax[1].legend()
777 plt.show()
778
779
780 def plot_erass_mc_results_hist(self, key, by='cycle', save=False):
781 try:
782 self.df_erass_mc_results['Z']
783 except:
784 self.table_erass_mc_results_map_info()
785
786 params = ['N_new_systems',
787 'N_old_system_become_transient',
788 'N_observed_ulxs',
789 'N_delta_obs_ulx',
790 'N_transients',
791 'N_transients_cum',
792 'N_total_systems']
793
794 res = self.MC_get_erass_mc_results(key)
795
796 if by == 'cycle':
797 for p in params:
798 print(f'Plotting {p}')
799 fig, axes = plt.subplots(2, 4, figsize=(12,10))
800 for i, ax in enumerate(axes.flat):
801 cycle = i+1
802 ax.set_title(f'eRASS_{cycle}_{p}')
803 sub = res[res['erass_cycle'] == cycle]
804 for b in unique_bhs:
805 sub1 = sub[sub['bh_ratio'] == b]
806 ax.hist(sub1[p], bins=np.arange(sub1[p].min(), sub1[p].max()+5, 1), label=f'%bh={b}', edgecolor='black', histtype='stepfilled', alpha=0.5)
807 ax.set_xlabel('N')
808 ax.legend()
809
810 if by == 'bh_ratio':
811 for p in params:
812 print(f'Plotting {p}')
813 fig, axes = plt.subplots(len(unique_bhs), 2, figsize=(12,10), sharex=True)
814 for i, b in enumerate(unique_bhs):
815 sub = res[res['bh_ratio'] == b]
816 for c in range(1, 9):
817 sub1 = sub[sub['erass_cycle'] == c]
818 sub_wind = sub1[sub1['period'] == 'P_wind']
819 sub_sup = sub1[sub1['period'] == 'P_sup']
820
821 axes[i][0].tick_params(axis='both', labelsize=8, labelbottom=True)
822 axes[i][1].tick_params(axis='both', labelsize=8, labelbottom=True)
823
824 axes[i][0].hist(sub_wind[p], bins=np.arange(sub_wind[p].min(), sub_wind[p].max()+1, 1), label=f'cycle={c}', edgecolor='black', histtype='stepfilled', alpha=0.5)
825 axes[i][1].hist(sub_sup[p], bins=np.arange(sub_sup[p].min(), sub_sup[p].max()+1, 1), label=f'cycle={c}', edgecolor='black', histtype='stepfilled', alpha=0.5)
826
827 axes[i][0].set_title(f'{p} | bh_ratio = {b} |P_wind')
828 # axes[i][0].set_xlabel('N')
829 axes[i][0].legend(fontsize=7)
830
831 axes[i][1].set_title(f'{p} | bh_ratio = {b} |P_sup')
832 # axes[i][1].set_xlabel('N')
833 axes[i][1].legend(fontsize=7)
834 plt.tight_layout()
835 if save==True:
836 plt.savefig(f'../reports/figures/e_{p}_{by}_{erass_system_period_cutoff}_{dincl_cutoff}_{Z}.png',dpi=500)
837 plt.savefig(f'../reports/figures/e_{p}_{by}_{erass_system_period_cutoff}_{dincl_cutoff}_{Z}.eps')
838 plt.savefig(f'../reports/figures/e_{p}_{by}_{erass_system_period_cutoff}_{dincl_cutoff}_{Z}.pdf')
839
840 def plot_bar_classifications_Z(self):
841 piv = self.table_classifications_pivot(margins=False)
842 piv.plot(kind='bar')
843
844 def plot_erass_transients(self, erass_system_period_cutoff=1460, Z='all', save=False):
845
846 cycles = ['1', '2', '3', '4', '5', '6', '7', '8']
847 dincls = [46, 20]
848 bh_percents = [0, 0.25, 0.5, 0.75, 1.0]
849
850 nbars = 8
851 spacing = 0.1
852 linewidth = 1.0
853
854 fig, ax = plt.subplots(ncols=2, nrows=len(dincls), figsize=(6,6), sharey='row')
855
856 ax[-1][0].set_xlabel('eRASS cycle')
857 ax[-1][1].set_xlabel('eRASS cycle')
858
859 clist = ["#ff595e", "#ffca3a", "#8ac926", "#1982c4", "#6a4c93"]
860
861 for i, dincl in enumerate(dincls):
862 for j, bh in enumerate(bh_percents):
863 key = (500, bh, dincl, Z, erass_system_period_cutoff)
864 df = self.MC_get_erass_mc_results(key)
865
866 sub_wind = df[df['period'] == 'P_wind']
867 sub_sup = df[df['period'] == 'P_sup']
868
869
870 stats_wind = sub_wind.groupby(['erass_cycle']).agg(['mean', 'std'])['N_transients']
871 stats_sup = sub_sup.groupby(['erass_cycle']).agg(['mean', 'std'])['N_transients']
872
873 bh_label = f'$\%_{{BH}}$ = {bh*100}'
874
875 trans1 = Affine2D().translate(-spacing*(nbars)/4 + spacing*j, 0.0) + ax[0][i].transData
876 trans2 = Affine2D().translate(-spacing*(nbars)/4 + spacing*j, 0.0) + ax[1][i].transData
877
878 ax[0][i].errorbar(cycles, stats_wind['mean'], yerr=stats_wind['std'], linestyle="none",
879 linewidth=linewidth, capsize=1.0, transform=trans1, label=bh_label, c=clist[j])
880
881 ax[1][i].errorbar(cycles, stats_sup['mean'], yerr=stats_sup['std'], linestyle="none",
882 linewidth=linewidth, capsize=1.0, transform=trans2, label=bh_label, c=clist[j])
883
884
885 ax[0][i].set_title(f'$\Delta i_{{max}} = {dincl}^{{\circ}}$ | $P_{{wind}}$ | Z = {Z}')
886 ax[1][i].set_title(f'$\Delta i_{{max}} = {dincl}^{{\circ}}$ | $P_{{sup}}$ | Z = {Z}')
887
888
889 ax[i][0].grid(axis='y')
890 ax[i][1].grid(axis='y')
891 #ax[1][i].grid(axis='y')
892 #ax[0][i].tick_params(labelrotation=90)
893 #ax[1][i].tick_params(axis='x', labelrotation=90)
894
895 ax[0][0].set_ylabel('N Transients')
896 ax[1][0].set_ylabel('N Transients')
897
898 ax[0][0].set_ylim(0,35)
899 ax[1][0].set_ylim(0,95)
900
901 ax[0][0].legend()
902 #ax[0][0].legend()
903
904 plt.subplots_adjust(wspace=0)
905
906 plt.tight_layout()
907 if save:
908 plt.savefig(f'../reports/figures/erass_N_transients_Z={Z}.eps', bbox_inches='tight')
909 plt.savefig(f'../reports/figures/erass_N_transients_Z={Z}.png', bbox_inches='tight', dpi=1000)
910
911
912
913
914if __name__ == "__main__":
915 # Load population
916 df = populations.startrack_v2_mt_1_all()
917 pop = populations.Population(df)
918
919 # Create sim object
920 rp = ResultsProcessor()
921
922 # Set Database
923 db_path = 'ulxlc2.db'
924 rp.set_active_db(db_path)
925
926 # Insert population
927 rp.set_parent_population(pop)
928
929 # # Run Simulations
930 # for Z in ['0.02', '0.002', '0.0002']:
931 # pop = populations.Population(df)
932 # rp.set_parent_population(pop)
933 # for i in range(10):
934 # for dincl in [21, 46]:
935 # for dc in [0.1, 0.2, 0.3, 1.0]:
936 # for bh_ratio in [0.0, 0.25, 0.5, 0.75, 1.0]:
937 # print(i, bh_ratio, dincl, Z, dc)
938 # # import pdb; pdb.set_trace()
939 # rp.MC_ERASS_simulation_run(size=500, dincl_cutoff=dincl, Z=Z, erass_lmxrb_duty_cycle=dc, erass_system_period_cutoff=999999)
940
941
942 # rp.MC_bh_ratio_classifications_sampler(N=10)
943
944 # Load tables
945 # rp.table_load_all()
946 rp.table_load_classifications()
947 # rp.table_load_transient()
948 rp.table_load_erass_mc_info()
949 # rp.table_load_erass_mc_results()
950
951 # Process tables
952 # rp.table_erass_mc_results_map_info()
953 # rp.table_classifications_map_systems()
954 # rp.table_classifications_pivot()
955
956 # rp.table_classifications_calc_intermediate()
957 # rp.table_classifications_map_systems()
958 # rp.table_classifications_pivot()
959
960 # MC related functions
961 # rp.MC_get_run_ids(group_precession_cuttoff=True)
962 # rp.MC_get_run_counts()
963
964 # Plot results
965
966 # rp.plot_set_latex_font()
967 # rp.plot_classifications_hist(Z='all', dincl_cut=20, frac_visible=True, save=True)
968 # rp.plot_erass_transients()
969
970
971 # for key in rp.dict_MC_run_ids.keys():
972 # print(key)
973
974
975 # rp.plot_erass_mc_results_hist(key)
976 # rp.plot_classifications_hist('0.002', 46, frac_visible=True, save=False)
977 # rp.plot_classifications_dincl_i()
978 # rp.plot_classifications_dincl_i(percent=True)
979 # rp.plot_classifications_hist('all', 20, frac_visible=True, save=True)
980
981 # rp.plot_erass_mc_results_hist(period='P_sup',
982 # by='bh_ratio',
983 # erass_system_period_cutoff=1460,
984 # dincl_cutoff=46,
985 # Z='0.002',
986 # save=True)
987
988 # rp.MC_calc_N_persistent()
989