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