Data¶
We need to set up our working directory with two subfolders, datafiles and datalogs.
datafiles will contain our RV catalogues. For each target or system we create a subfolder with the system name. In this case, 51Peg. Inside, we create a second subfolder, named RV, which will contain the data to be read.
We copy-paste the file downloaded from GitHub into /datafiles/51Peg/RV/.
📂working_directory
┣ 📜mini_test.py
┣ 📂datafiles
┃ ┣ 📂51Peg
┃ ┃ ┗ 📂RV
┃ ┃ ┗ 📜51peg.vels
┣ 📂datalogs
┃ ┣ 📂51Peg
┃ ┃ ┗ 📂run_1
Setting up EMPEROR¶
Under our working directory, we create a python file named mini_test.
First, we import the library and start our simulation:
In [1]:
Copied!
import astroemperor as emp
import numpy as np
np.random.seed(1234)
sim = emp.Simulation()
sim.load_data('51Peg') # folder read from /datafiles/
import astroemperor as emp
import numpy as np
np.random.seed(1234)
sim = emp.Simulation()
sim.load_data('51Peg') # folder read from /datafiles/
I couldnt grab the terminal size. Trying with pandas... Terminal size with pandas successful! ~~ Simulation Successfully Initialized ~~ Reading data from 51peg.vels
Setting the engine¶
For this example, we will use reddemcee, with 10 temperatures, 256 walkers, 2048 sweeps each of 1 step:
In [2]:
Copied!
sim.set_engine('reddemcee')
sim.engine_config['setup'] = [10, 256, 2048, 1]
sim.set_engine('reddemcee')
sim.engine_config['setup'] = [10, 256, 2048, 1]
Setting the model¶
We feed the name of the instrument (optional), as well as the starmass for calculating the minimum-mass and semi-major axis. We will use the Keplerian parameterisation ($P$, $K$, $\phi$, $e_{s}$, $e_{c}$). We add some boundaries to speed up the process, and add some initial positions:
In [3]:
Copied!
sim.instrument_names_RV = ['LICK']
sim.starmass = 1.12 # in solar masses
sim.keplerian_parameterisation = 1 # see above description
# Initial Positions
sim.add_condition(['Period 1', 'limits', [3, 5]])
sim.add_condition(['Amplitude 1', 'limits', [45, 60]])
sim.add_condition(['Offset 1', 'limits', [-10., 10.]])
sim.add_condition(['Period 1', 'init_pos', [4.1, 4.3]])
sim.add_condition(['Amplitude 1', 'init_pos', [50, 60]])
sim.instrument_names_RV = ['LICK']
sim.starmass = 1.12 # in solar masses
sim.keplerian_parameterisation = 1 # see above description
# Initial Positions
sim.add_condition(['Period 1', 'limits', [3, 5]])
sim.add_condition(['Amplitude 1', 'limits', [45, 60]])
sim.add_condition(['Offset 1', 'limits', [-10., 10.]])
sim.add_condition(['Period 1', 'init_pos', [4.1, 4.3]])
sim.add_condition(['Amplitude 1', 'init_pos', [50, 60]])
Plotting Options¶
We add some plotting options to speed up this test a little. We will plot the posteriors for the cold chain (0), and two intermediate chains (2, and 6). Also, we won't use the arviz optional plots.
In [4]:
Copied!
sim.plot_posteriors['temps'] = [0, 2, 6]
sim.plot_trace['plot'] = False # arviz plots
sim.plot_posteriors['temps'] = [0, 2, 6]
sim.plot_trace['plot'] = False # arviz plots
Finally, we run our simulation (it will take some minutes):
In [5]:
Copied!
sim.autorun(0, 1)
sim.autorun(0, 1)
Offset block added, OffsetBlock Jitter block added, JitterBlock Condition applied: Parameter Offset 1 attribute limits set to [-10.0, 10.0] ~~ Setup Info ~~ Current Engine is reddemcee 0.9.12 Number of cores is 24 Save location is datalogs/51Peg/run_5/k0 Dynamical Criteria is None Posterior fit method is Gaussian Mixtures Limits constrain method is range Model Selection method is BIC ~~ Automatically Saving ~~ Logger : ✔ Samples : ✘ Posteriors : ✔ Likelihoods : ✔ Plots: Posteriors : ✔ Plots: Keplerian Model : ✔ Plots: Gaussian Mixture : ✔ Plots: Parameter Histograms : ✔ Plots: Corner : ✔ ~~ Pre-Run Info ~~ Parameter Prior Limits ----------- ------------------ --------------------------------- Offset 1 ~𝓤 (-10.0, 10.0) [-10. 10.] ---------- ------------------ --------------------------------- Jitter 1 ~𝓝 (5, 5) [ 0. 75.852] Math for OffsetBlock: γ₀|ᵢ Math for JitterBlock: 𝝈ᵢ Generating Samples
100%|██████████| 20480/20480 [03:25<00:00, 99.51it/s]
temp_script.py took 206.036 seconds Autocorrelation tolerance=50 fails. Setting to 0.
/IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:294: RuntimeWarning: overflow encountered in exp C_series = np.exp(logC_series) /IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:54: RuntimeWarning: invalid value encountered in subtract bm = (csum[b:] - csum[:-b]) / b /IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:55: RuntimeWarning: invalid value encountered in subtract diff = bm - mu_hat /IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:56: RuntimeWarning: overflow encountered in matmul Sigma_hat = (b / K) * diff.T @ diff
Are you running a short chain? Evidence with hybrid failed! Trying fallback to TI+.
Calculating Gaussian Mixtures
100%|██████████| 2/2 [00:01<00:00, 1.02it/s]
~~ Best Fit ~~ Parameter Value(max) Range(±σ) Prior Limits ----------- ------------ --------------- ---------------- ----------- γ₀₁ -0.156 [-1.121 1.024] ~𝓤 (-10.0, 10.0) [-10. 10.] ----------- ------------ --------------- --------- --------------- J₁ 36.067 [-0.549 0.789] ~𝓝 (5, 5) [ 0. 75.852] ~~ Run Info ~~ Info Value ----------------------------------- -------------------------------------------------------------------------------- Star Name : 51Peg The sample sizes are : [524288, 524288, 524288, 524288, 524288, 524288, 524288, 524288, 524288, 524288] Temps, Walkers, Sweeps, Steps : [10, 256, 2048, 1] Model used is : ['OffsetBlock', 'JitterBlock'](2) N data : 256 t0 epoch is : 2450002.665695 Number of Dimensions : 2 Degrees of Freedom : 254 ---------------------------------------- Decay Timescale, Rate, Scheme : 1000, 1, SAR Beta Detail : [1.0, 0.4907, 0.2626, 0.1394, 0.07151, 0.03387, 0.01583, 0.005525, 0.001317, 2.478e-08] Mean Logl Detail : [-1309.564, -1313.947, -1323.203, -1341.768, -1377.180, -1439.486, -1550.105, -1766.061, -2302.283, -3832.927] Mean Acceptance Fraction : [0.713, 0.712, 0.711, 0.711, 0.707, 0.701, 0.692, 0.682, 0.676, 0.657] Autocorrelation Time : [16.363, 19.295] Temperature Swap Rate : [0.338, 0.339, 0.341, 0.343, 0.346, 0.349, 0.354, 0.358, 0.362] Mean Swap Distance : [1.165, 1.512, 1.811, 2.083, 2.269, 2.461, 2.644, 2.809, 2.996] ~~ Statistical Details ~~ Statistic Value ----------------------------- ------------------ The evidence is : -1330.864 +- 3.157 The maximum posterior is : -1333.444 The maximum likelihood is : -1308.790 The BIC is : 2628.671 The AIC is : 2621.581 The DIC is : 2680.599 The HQIC is : 2624.433 The Bayes Factor is : 22.074 The chi2 is : 302.629 The reduced chi2 is : 1.191 The RMSE is : 39.917 The RMSi is : [39.917] The Weights are : [1.] Plotting Posterior Scatter Plot
100%|██████████| 4/4 [00:03<00:00, 1.27it/s] 100%|██████████| 2/2 [00:02<00:00, 1.50s/it]
Plotting Histograms Plot
100%|██████████| 20/20 [00:05<00:00, 3.51it/s]
Plotting Keplerian Models Plotting E[log L](beta) Plot
100%|██████████| 1/1 [00:00<00:00, 7.36it/s]
Plotting Beta Density
100%|██████████| 1/1 [00:00<00:00, 5.45it/s]
Plotting Temperature Rates
100%|██████████| 9/9 [00:00<00:00, 10.24it/s]
Plotting Gaussian Mixtures
100%|██████████| 2/2 [00:00<00:00, 4.76it/s]
Time Table Time RUN : 00:03:28 Time POSTPROCESS : 00:00:04 Time CALCULATE GM : 00:00:01 Time plot_posteriors : 00:00:06 Time plot_histograms : 00:00:05 Time plot_keplerian_model : 00:00:00 Time plot_betas : 00:00:00 Time plot_beta_density : 00:00:00 Time plot_rates : 00:00:00 Time plot_trace : 00:00:00 past BIC - present BIC > 5 BIC condition met!! inf - 2628.671 > 5 ~~ Proceeding with the next run ! ~~ Keplerian block added, KeplerianBlock 1 Condition applied: Parameter Period 1 attribute limits set to [3, 5] Condition applied: Parameter Period 1 attribute init_pos set to [4.1, 4.3] Condition applied: Parameter Amplitude 1 attribute limits set to [45, 60] Condition applied: Parameter Amplitude 1 attribute init_pos set to [50, 60] Condition applied: Parameter Jitter 1 attribute limits set to [1e-05, 39.126244622733076] ~~ Pre-Run Info ~~ Parameter Prior Limits ------------- ------------------ --------------------------------- Period 1 ~𝓤 (3, 5) [3. 5.] Amplitude 1 ~𝓤 (45, 60) [45. 60.] Phase 1 ~𝓤 (0.0, 6.283) [0. 6.283] Ecc_sin 1 ~𝓤 (-1, 1) [-1. 1.] Ecc_cos 1 ~𝓤 (-1, 1) [-1. 1.] ------------- ------------------ --------------------------------- Offset 1 ~𝓤 (-10.0, 10.0) [-10. 10.] ------------- ------------------ --------------------------------- Jitter 1 ~𝓝 (5.0, 5.0) [ 0. 39.126] Math for KeplerianBlock 1: K⋅(cos(ν(t,P,𝜙,e)+𝜔 )+e⋅cos(𝜔 ))|₁ Math for OffsetBlock: γ₀|ᵢ Math for JitterBlock: 𝝈ᵢ Generating Samples
100%|██████████| 20480/20480 [03:36<00:00, 94.56it/s]
temp_script.py took 216.75 seconds Autocorrelation tolerance=50 fails. Setting to 0.
/IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:294: RuntimeWarning: overflow encountered in exp C_series = np.exp(logC_series) /IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:54: RuntimeWarning: invalid value encountered in subtract bm = (csum[b:] - csum[:-b]) / b /IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:55: RuntimeWarning: invalid value encountered in subtract diff = bm - mu_hat
Are you running a short chain? Evidence with hybrid failed! Trying fallback to TI+.
Calculating Gaussian Mixtures
100%|██████████| 11/11 [00:10<00:00, 1.07it/s]
~~ Best Fit ~~ Parameter Value(max) Range(±σ) Prior Limits ----------- ------------ --------------- --------------- ------------- P₁ 4.231 [-0. 0.] ~𝓤 (3, 5) [3. 5.] K₁ 55.95 [-0.29 0.331] ~𝓤 (45, 60) [45. 60.] M₀₁ 1.681 [-0.855 1.125] ~𝓤 (0.0, 6.283) [0. 6.283] eₛᵢₙ₁ 0.101 [-0.062 0.021] ~𝓤 (-1, 1) [-1. 1.] eₖₒₛ₁ 0.058 [-0.072 0.061] ~𝓤 (-1, 1) [-1. 1.] e₁ 0.014 [-0.005 0.004] ~ [0. 1.] ω₁ 1.05 [-1.036 2.4 ] ~ [0. 6.283] a₁ 0.053 [-0. 0.] ~ [ 0. 1000.] Mₘᵢₙ₁ 0.48 [-0.002 0.003] ~ [ 0. 1000.] ----------- ------------ --------------- ---------------- ----------- γ₀₁ 3.807 [-0.316 0.136] ~𝓤 (-10.0, 10.0) [-10. 10.] ----------- ------------ --------------- ------------- --------------- J₁ 3.004 [-0.181 0.571] ~𝓝 (5.0, 5.0) [ 0. 39.126] ~~ Run Info ~~ Info Value ----------------------------------- -------------------------------------------------------------------------------- Star Name : 51Peg The sample sizes are : [524288, 524288, 524288, 524288, 524288, 524288, 524288, 524288, 524288, 524288] Temps, Walkers, Sweeps, Steps : [10, 256, 2048, 1] Model used is : ['KeplerianBlock 1', 'OffsetBlock', 'JitterBlock'](7) N data : 256 t0 epoch is : 2450002.665695 Number of Dimensions : 7 Degrees of Freedom : 249 ---------------------------------------- Decay Timescale, Rate, Scheme : 1000, 1, SAR Beta Detail : [1.0, 0.4322, 0.1784, 0.08029, 0.03487, 0.01552, 0.009404, 0.004159, 0.001062, 2.478e-08] Mean Logl Detail : [-884.153, -894.830, -911.685, -941.604, -997.206, -1138.482, -1427.055, -1823.948, -2688.347, -5583.611] Mean Acceptance Fraction : [0.210, 0.202, 0.201, 0.203, 0.212, 0.190, 0.156, 0.183, 0.216, 0.234] Autocorrelation Time : [11.932, 71.329, 152.548, 127.035, 172.767, 75.182, 109.901] Temperature Swap Rate : [0.284, 0.284, 0.283, 0.282, 0.282, 0.281, 0.281, 0.281, 0.281] Mean Swap Distance : [0.600, 0.856, 1.198, 1.710, 2.346, 3.062, 3.028, 2.950, 2.992] ~~ Statistical Details ~~ Statistic Value ----------------------------- ----------------- The evidence is : -925.677 +- 5.897 The maximum posterior is : -881.635 The maximum likelihood is : -869.579 The BIC is : 1777.974 The AIC is : 1753.158 The DIC is : 14718.131 The HQIC is : 1763.139 The Bayes Factor is : 56.098 The chi2 is : 267.427 The reduced chi2 is : 1.074 The RMSE is : 7.614 The RMSi is : [7.614] The Weights are : [1.] Plotting Posterior Scatter Plot
100%|██████████| 4/4 [00:11<00:00, 2.98s/it] 100%|██████████| 4/4 [00:12<00:00, 3.17s/it] 100%|██████████| 1/1 [00:02<00:00, 2.71s/it]
Plotting Histograms Plot
40%|████ | 12/30 [00:07<00:09, 1.95it/s]
Failed to plot the Period 1 histogram
100%|██████████| 30/30 [00:19<00:00, 1.57it/s]
Plotting Keplerian Models Plotting E[log L](beta) Plot
100%|██████████| 1/1 [00:00<00:00, 7.68it/s]
Plotting Beta Density
100%|██████████| 1/1 [00:00<00:00, 4.89it/s]
Plotting Temperature Rates
100%|██████████| 9/9 [00:00<00:00, 9.75it/s]
Plotting Gaussian Mixtures
100%|██████████| 11/11 [00:02<00:00, 4.43it/s]
Time Table Time RUN : 00:03:39 Time POSTPROCESS : 00:00:15 Time CALCULATE GM : 00:00:10 Time plot_posteriors : 00:00:28 Time plot_histograms : 00:00:19 Time plot_keplerian_model : 00:00:00 Time plot_betas : 00:00:00 Time plot_beta_density : 00:00:00 Time plot_rates : 00:00:00 Time plot_trace : 00:00:00 past BIC - present BIC > 5 BIC condition met!! 2628.671 - 1777.974 > 5 ~~ End of the Run ~~
Now we can explore the results in the /datalogs/51Peg/run_1/ folder:
In [ ]:
Copied!
from IPython.display import Image
loc = sim.saveplace
Image(filename=f'{loc}/plots/models/KeplerianBlock_1.png')
from IPython.display import Image
loc = sim.saveplace
Image(filename=f'{loc}/plots/models/KeplerianBlock_1.png')
Out[ ]:
