Fitting️
MRExo
can be used to fit samples with two to four different observables (along with their measurement uncertainties).
An example of how to load a sample and perform a fit in three dimensions is included in sample_fit.py , where the same script can be adapted to fit two or four dimensions.
Following this, inference can be performed following the examples here to predict one observable from another.
1import os, sys
2from astropy.table import Table
3import numpy as np
4import numpy as np
5import pandas as pd
6import shutil
7
8from mrexo.mle_utils_nd import InputData, MLE_fit
9from mrexo.fit_nd import fit_relation
10from mrexo.plotting_nd import Plot2DJointDistribution, Plot2DWeights, Plot1DInputDataHistogram
11import matplotlib.pyplot as plt
12
13Platform = sys.platform
14
15if Platform == 'win32':
16 HomeDir = 'C:\\Users\\skanodia\\Documents\\\\GitHub\\'
17else:
18 HomeDir = r"/storage/home/szk381/work/"
19 #HomeDir = r"/home/skanodia/work/"
20
21
22try :
23 pwd = os.path.dirname(__file__)
24except NameError:
25 pwd = os.path.join(HomeDir, 'mrexo', 'sample_scripts')
26 print('Could not find pwd')
27
28# Directory with dataset to be fit
29DataDirectory = os.path.join(HomeDir, 'Mdwarf-Exploration', 'Data', 'MdwarfPlanets')
30print(DataDirectory)
31
32t = pd.read_csv(os.path.join(DataDirectory, 'Teff_7000_ExcUpperLimits_20230306RpLt20.csv'))
33
34# Mask NaNs
35t = t[~np.isnan(t['pl_insolerr1'])]
36t = t[~np.isnan(t['pl_masse'])]
37
38# Define bounds in different dimensions
39RadiusBounds = [0, 10]# None# [0, 100]
40MassBounds = None# [0, 6000]
41InsolationBounds = None# [0.01, 5000]
42StellarMassBounds = None# [0.2, 1.2]
43
44# Measurements with very small errors are set to NaN to avoid integration errors
45t['st_masserr1'][t['st_masserr1'] < 0.005] = np.nan
46t['st_masserr2'][t['st_masserr2'] < 0.005] = np.nan
47
48if RadiusBounds is not None:
49 t = t[(t.pl_rade > RadiusBounds[0]) & (t.pl_rade < RadiusBounds[1])]
50
51if MassBounds is not None:
52 t = t[(t.pl_masse > MassBounds[0]) & (t.pl_masse < MassBounds[1])]
53
54if InsolationBounds is not None:
55 t = t[(t.pl_insol > InsolationBounds[0]) & (t.pl_insol < InsolationBounds[1])]
56
57if StellarMassBounds is not None:
58 t = t[(t.st_mass > StellarMassBounds[0]) & (t.st_mass < StellarMassBounds[1])]
59
60# Remove particular planets
61RemovePlanets = ['Kepler-54 b', 'Kepler-54 c']
62t = t[~np.isin(t.pl_name, RemovePlanets)]
63
64print(len(t))
65
66# In Earth units
67Mass = np.array(t['pl_masse'])
68# Asymmetrical errorbars
69MassUSigma = np.array(abs(t['pl_masseerr1']))
70MassLSigma = np.array(abs(t['pl_masseerr2']))
71
72Radius = np.array(t['pl_rade'])
73# Asymmetrical errorbars
74RadiusUSigma = np.array(abs(t['pl_radeerr1']))
75RadiusLSigma = np.array(abs(t['pl_radeerr2']))
76
77StellarMass = np.array(t['st_mass'])
78StellarMassUSigma = np.array(t['st_masserr1'])
79StellarMassLSigma = np.array(t['st_masserr2'])
80
81Insolation = np.array(t['pl_insol'])
82InsolationUSigma = np.array(t['pl_insolerr1'])
83InsolationLSigma = np.array(t['pl_insolerr2'])
84
85# Let the script pick the max and min bounds, or can hard code those in. Note that the dataset must fall within the bounds if they are specified.
86Max, Min = np.nan, np.nan
87
88# Define input dictionary for each dimension
89RadiusDict = {'Data': Radius, 'LSigma': RadiusLSigma, "USigma":RadiusUSigma, 'Max':Max, 'Min':Min, 'Label':'Radius ($R_{\oplus}$)', 'Char':'r'}
90MassDict = {'Data': Mass, 'LSigma': MassLSigma, "USigma":MassUSigma, 'Max':Max, 'Min':Min, 'Label':'Mass ($M_{\oplus}$)', 'Char':'m'}
91# PeriodDict = {'Data': Period, 'LSigma': PeriodSigma, "USigma":PeriodSigma, 'Max':Max, 'Min':Min, 'Label':'Period (d)', 'Char':'p'}
92StellarMassDict = {'Data': StellarMass, 'LSigma': StellarMassLSigma, "USigma":StellarMassUSigma, 'Max':Max, 'Min':Min, 'Label':'Stellar Mass (M$_{\odot}$)', 'Char':'stm'}
93InsolationDict = {'Data': Insolation, 'LSigma': InsolationLSigma, "USigma":InsolationUSigma, 'Max':Max, 'Min':Min, 'Label':'Pl Insol ($S_{\oplus}$)', 'Char':'insol'}
94
95# 3D fit with planetary radius, mass and insolation
96InputDictionaries = [RadiusDict, MassDict, InsolationDict]
97DataDict = InputData(InputDictionaries)
98ndim = len(InputDictionaries)
99
100RunName = 'AllPlanet_RpLt20_MRS_test'
101save_path = os.path.join(pwd, 'TestRuns', RunName)
102
103# Harcode the number of degrees per dimension. Read the `fit_relation()` documentation for alternatives
104deg_per_dim = [25, 25, 25]
105
106if __name__ == '__main__':
107
108 outputs= fit_relation(DataDict, select_deg=select_deg, save_path=save_path, degree_max=30, cores=2,SymmetricDegreePerDimension=True, NumMonteCarlo=0, NumBootstrap=0)
109 shutil.copy(os.path.join(pwd, 'sample_fit.py'), os.path.join(save_path, 'sample_fit_{}.py'.format(RunName)))
110
111 _ = Plot1DInputDataHistogram(save_path)
112
113 if ndim==2:
114
115 _ = Plot2DJointDistribution(save_path)
116 _ = Plot2DWeights(save_path)
117