51 Peg¶

This section is a hands-on tutorial on how to make a simple run. We use the 51 Peg RV data available on GitHub.

In [1]:

Copied!





import numpy as np
import astroemperor as emp

np.random.seed(1234)

target_folder = 'HD55693'
my_cores = 24
eng = 'reddemcee'
import numpy as np
import astroemperor as emp

np.random.seed(1234)

target_folder = 'HD55693'
my_cores = 24
eng = 'reddemcee'

I couldnt grab the terminal size. Trying with pandas...
Terminal size with pandas successful!

In [2]:

Copied!





sim = emp.Simulation()

sim.cores__ = my_cores
sim.acceleration = 0
sim.keplerian_parameterisation = 0

sim.starmass = 1.109 # HD55693 http://exoplanets.astro.yale.edu/science/analysis/spocs_evol_95.php
sim.instrument_names_RV = ['TERRA1',
                           'PFS',
                           'TERRA2']

sim.jitter_prargs = [1, 5]
sim.eccentricity_prargs = [0, 0.1]

if True:
    sim.add_condition(['Period 1', 'limits', [2000, 3000]])
    sim.add_condition(['Amplitude 1', 'limits', [2., 5.]])
    sim.add_condition(['Eccentricity 1', 'limits', [0., 0.2]])

    sim.add_condition(['Period 2', 'limits', [1.0, 300]])
    sim.add_condition(['Period 3', 'limits', [1.0, 300]])

    sim.add_condition(['Period 1', 'init_pos', [2500, 2700]])
    sim.add_condition(['Amplitude 1', 'init_pos', [2.5, 4.5]])
    sim.add_condition(['Eccentricity 1', 'init_pos', [0, 0.1]])
    #sim.add_condition(['Phase 1', 'init_pos', [4, 4.8]])
    #sim.add_condition(['Longitude 1', 'init_pos', [2.1, 3.]])

    sim.add_condition(['Period 2', 'init_pos', [10, 40]])
    sim.add_condition(['Amplitude 2', 'init_pos', [1.5, 2.5]])
    sim.add_condition(['Eccentricity 2', 'init_pos', [0, 0.3]])

    sim.add_condition(['Period 2', 'init_pos', [10, 150]])
    sim.add_condition(['Amplitude 2', 'init_pos', [1.5, 2.5]])
    sim.add_condition(['Eccentricity 2', 'init_pos', [0, 0.3]])

if True:
    sim.plot_posteriors['plot'] = True
    sim.plot_posteriors['modes'] = [0]
    sim.plot_posteriors['temps'] = [0]
    sim.plot_posteriors['format'] = 'png'

    sim.plot_rates['window'] = 10
    sim.plot_histograms['plot'] = False
    
    sim.plot_trace['plot'] = False

    sim.plot_betas['plot'] = True
    sim.plot_rates['plot'] = True
    sim.plot_rates['window'] = 10

    sim.plot_keplerian_model['uncertain'] = False
    sim.plot_all['format'] = 'png'
    sim.plot_all['paper_mode'] = True


sim.set_engine('reddemcee')
sim = emp.Simulation()

sim.cores__ = my_cores
sim.acceleration = 0
sim.keplerian_parameterisation = 0

sim.starmass = 1.109 # HD55693 http://exoplanets.astro.yale.edu/science/analysis/spocs_evol_95.php
sim.instrument_names_RV = ['TERRA1',
                           'PFS',
                           'TERRA2']

sim.jitter_prargs = [1, 5]
sim.eccentricity_prargs = [0, 0.1]

if True:
    sim.add_condition(['Period 1', 'limits', [2000, 3000]])
    sim.add_condition(['Amplitude 1', 'limits', [2., 5.]])
    sim.add_condition(['Eccentricity 1', 'limits', [0., 0.2]])

    sim.add_condition(['Period 2', 'limits', [1.0, 300]])
    sim.add_condition(['Period 3', 'limits', [1.0, 300]])

    sim.add_condition(['Period 1', 'init_pos', [2500, 2700]])
    sim.add_condition(['Amplitude 1', 'init_pos', [2.5, 4.5]])
    sim.add_condition(['Eccentricity 1', 'init_pos', [0, 0.1]])
    #sim.add_condition(['Phase 1', 'init_pos', [4, 4.8]])
    #sim.add_condition(['Longitude 1', 'init_pos', [2.1, 3.]])

    sim.add_condition(['Period 2', 'init_pos', [10, 40]])
    sim.add_condition(['Amplitude 2', 'init_pos', [1.5, 2.5]])
    sim.add_condition(['Eccentricity 2', 'init_pos', [0, 0.3]])

    sim.add_condition(['Period 2', 'init_pos', [10, 150]])
    sim.add_condition(['Amplitude 2', 'init_pos', [1.5, 2.5]])
    sim.add_condition(['Eccentricity 2', 'init_pos', [0, 0.3]])

if True:
    sim.plot_posteriors['plot'] = True
    sim.plot_posteriors['modes'] = [0]
    sim.plot_posteriors['temps'] = [0]
    sim.plot_posteriors['format'] = 'png'

    sim.plot_rates['window'] = 10
    sim.plot_histograms['plot'] = False
    
    sim.plot_trace['plot'] = False

    sim.plot_betas['plot'] = True
    sim.plot_rates['plot'] = True
    sim.plot_rates['window'] = 10

    sim.plot_keplerian_model['uncertain'] = False
    sim.plot_all['format'] = 'png'
    sim.plot_all['paper_mode'] = True


sim.set_engine('reddemcee')

                                                                                
                   ~~ Simulation Successfully Initialized ~~

In [3]:

Copied!

sim.engine_config['setup'] = [8, 64, 256, 1]  # ntemps, nwalkers, nsweeps, nsteps
sim.engine_config['setup'] = [8, 64, 256, 1]  # ntemps, nwalkers, nsweeps, nsteps

In [4]:

Copied!

sim.load_data('HD55693')  # folder read from /datafiles/
sim.autorun(2, 2)
sim.load_data('HD55693')  # folder read from /datafiles/
sim.autorun(2, 2)



Reading data from HD55693_1_TERRA_select.vels
Reading data from HD55693_2_PFS.vels
Reading data from HD55693_3_TERRA2_select.vels






                              Offset block added, OffsetBlock


                              Jitter block added, JitterBlock


                              Keplerian block added, KeplerianBlock 1


                              Keplerian block added, KeplerianBlock 2



Condition applied: Parameter Period 1 attribute limits set to [2000, 3000]

Condition applied: Parameter Period 1 attribute init_pos set to [2500, 2700]

Condition applied: Parameter Amplitude 1 attribute limits set to [2.0, 5.0]

Condition applied: Parameter Amplitude 1 attribute init_pos set to [2.5, 4.5]

Condition applied: Parameter Eccentricity 1 attribute limits set to [0.0, 0.2]

Condition applied: Parameter Eccentricity 1 attribute init_pos set to [0, 0.1]

Condition applied: Parameter Period 2 attribute limits set to [1.0, 300]

Condition applied: Parameter Period 2 attribute init_pos set to [10, 40]

Condition applied: Parameter Period 2 attribute init_pos set to [10, 150]

Condition applied: Parameter Amplitude 2 attribute init_pos set to [1.5, 2.5]

Condition applied: Parameter Eccentricity 2 attribute init_pos set to [0, 0.3]


                                ~~ Setup Info ~~                                


Current Engine is            reddemcee 0.9.12
Number of cores is           24
Save location is             datalogs/HD55693/run_2/k2
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
----------------  --------------------  -----------------
Period 1          ~𝓤 (2000, 3000)       [2000, 3000]
Amplitude 1       ~𝓤 (2.0, 5.0)         [2, 5]
Phase 1           ~𝓤 (0.0, 6.283)       [0, 6.283]
Eccentricity 1    ~𝓝 (0.0, 0.1)         [0, 0.2]
Longitude 1       ~𝓤 (0.0, 6.283)       [0, 6.283]

----------------  --------------------  -----------------
Period 2          ~𝓤 (1.0, 300.0)       [1, 300]
Amplitude 2       ~𝓤 (0.0, 7.07)        [0, 7.07]
Phase 2           ~𝓤 (0.0, 6.283)       [0, 6.283]
Eccentricity 2    ~𝓝 (0.0, 0.1)         [0, 1]
Longitude 2       ~𝓤 (0.0, 6.283)       [0, 6.283]

----------------  --------------------  -----------------
Offset 1          ~𝓤 (-8.488, 8.488)    [-8.488, 8.488]
Offset 2          ~𝓤 (-6.684, 6.684)    [-6.684, 6.684]
Offset 3          ~𝓤 (-5.74, 5.74)      [-5.74, 5.74]

----------------  --------------------  -----------------
Jitter 1          ~𝓝 (1, 5)             [0, 8.488]
Jitter 2          ~𝓝 (1, 5)             [0, 6.684]
Jitter 3          ~𝓝 (1, 5)             [0, 5.74]




Math for KeplerianBlock 1: 
                       K⋅(cos(ν(t,P,𝜙,e)+𝜔 )+e⋅cos(𝜔 ))|₁                       
Math for KeplerianBlock 2: 
                       K⋅(cos(ν(t,P,𝜙,e)+𝜔 )+e⋅cos(𝜔 ))|₂                       
Math for OffsetBlock: 
                                      γ₀|ᵢ                                      
Math for JitterBlock: 
                                       𝝈ᵢ                                       
                               Generating Samples

100%|██████████| 2048/2048 [00:11<00:00, 183.96it/s]

temp_script.py took 11.259 seconds
Autocorrelation tolerance=50 fails. Setting to 0.

/IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/utils.py:56: RuntimeWarning: overflow encountered in matmul
  Sigma_hat = (b / K) * diff.T @ diff

                         Calculating Gaussian Mixtures

100%|██████████| 20/20 [00:00<00:00, 26.16it/s]


                                 ~~ Best Fit ~~                                 


Parameter           Value(max)  Range(±σ)               Prior                 Limits
----------------  ------------  ----------------------  --------------------  -----------------
Period 1              2609.92   [2566.006, 2656.397]    ~𝓤 (2000, 3000)       [2000, 3000]
Amplitude 1              3.814  [3.343, 4.294]          ~𝓤 (2.0, 5.0)         [2, 5]
Phase 1                  3.406  [3.011, 3.939]          ~𝓤 (0.0, 6.283)       [0, 6.283]
Eccentricity 1           0.073  [0.041, 0.094]          ~𝓝 (0.0, 0.1)         [0, 0.2]
Longitude 1              3.805  [2.755, 4.96]           ~𝓤 (0.0, 6.283)       [0, 6.283]

----------------  ------------  ----------------------  --------------------  -----------------
Period 2                55.847  [27.53, 91.987]         ~𝓤 (1.0, 300.0)       [1, 300]
Amplitude 2              2.279  [1.865, 2.527]          ~𝓤 (0.0, 7.07)        [0, 7.07]
Phase 2                  3.407  [2.527, 3.711]          ~𝓤 (0.0, 6.283)       [0, 6.283]
Eccentricity 2           0.089  [0.029, 0.185]          ~𝓝 (0.0, 0.1)         [0, 1]
Longitude 2              3.904  [2.129, 5.122]          ~𝓤 (0.0, 6.283)       [0, 6.283]

----------------  ------------  ----------------------  --------------------  -----------------
Offset 1                -0.402  [-1.408, -0.363]        ~𝓤 (-8.488, 8.488)    [-8.488, 8.488]
Offset 2                 0.391  [-0.648, 1.257]         ~𝓤 (-6.684, 6.684)    [-6.684, 6.684]
Offset 3                -0.203  [-0.978, 0.227]         ~𝓤 (-5.74, 5.74)      [-5.74, 5.74]

----------------  ------------  ----------------------  --------------------  -----------------
Jitter 1                 2.459  [1.474, 4.102]          ~𝓝 (1, 5)             [0, 8.488]
Jitter 2                 2.808  [2.467, 3.083]          ~𝓝 (1, 5)             [0, 6.684]
Jitter 3                 2.345  [1.668, 3.309]          ~𝓝 (1, 5)             [0, 5.74]




                                 ~~ Run Info ~~                                 






Info                                 Value
-----------------------------------  --------------------------------------------------------------------------
Star Name                      :     HD55693
The sample sizes are           :     [16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384]
Temps, Walkers, Sweeps, Steps  :     [8, 64, 256, 1]
Model used is                  :     ['KeplerianBlock 1', 'KeplerianBlock 2', 'OffsetBlock', 'JitterBlock'](16)
N data                         :     84
t0 epoch is                    :     2452940.8583322
Number of Dimensions           :     16
Degrees of Freedom             :     68


----------------------------------------



Decay Timescale, Rate, Scheme   :   1000, 1, SAR

Beta Detail                     :   [1.0, 0.7759, 0.5287, 0.3583, 0.2062, 0.1189, 0.05119, 0.01627]

Mean Logl Detail                :   [-210.491, -216.626, -224.466, -233.610, -245.816, -264.855, -300.731, -416.933]

Mean Acceptance Fraction        :   [0.061, 0.083, 0.102, 0.121, 0.128, 0.138, 0.142, 0.141]

Autocorrelation Time            :   [18.624, 16.376, 17.371, 16.884, 17.285, 14.247, 17.946, 15.208, 15.330, 17.127, 13.427, 11.199, 14.609, 21.159, 15.203, 14.349]

Temperature Swap Rate           :   [0.397, 0.401, 0.405, 0.408, 0.412, 0.416, 0.417]

Mean Swap Distance              :   [29.834, 30.597, 36.936, 40.783, 42.391, 44.823, 48.487]




                           ~~ Statistical Details ~~                            




Statistic                     Value
----------------------------  -------------------
Evidence   ln(Z)        :     -590.940 +- inf
Posterior  ln(P) (max)  :     -225.854
Likelihood ln(L) (max)  :     -193.147
Bayes Factor            :     397.793
BIC                     :     457.186
AIC                     :     418.293
DIC                     :     3680.372
HQIC                    :     433.928
χ²                      :     65.245
Reduced χ²              :     0.959
RMSE                    :     2.326
RMSi                    :     [2.257 2.209 2.628]
Weights (RMS)           :     [0.359 0.375 0.265]






                        Plotting Posterior Scatter Plot

100%|██████████| 4/4 [00:01<00:00,  2.35it/s]
0it [00:00, ?it/s]


                           Plotting Keplerian Models                            
data1=            BJD        RV       eRV  Flag    eRV_og  residuals
0      0.000000  4.912826  2.541371     1  0.640157   1.624782
1     59.930089  4.434408  2.668277     1  1.034862   2.361263
2    108.761071  3.277795  2.481395     1  0.329479   0.254507
3    397.936897 -0.775525  2.511916     1  0.510838   0.506437
4    397.942036 -0.538872  2.511765     1  0.510097   0.744319
..          ...       ...       ...   ...       ...        ...
31  6662.795638 -6.743544  3.102287     2  1.319130  -3.395605
32  6662.800288 -6.294784  3.115435     2  1.349760  -2.946442
33  6711.647318  2.666146  3.075738     2  1.255420   4.299660
34  6711.649908 -1.291774  3.072120     2  1.246530   0.342292
35  6711.652458  1.802906  3.071507     2  1.245020   3.437516

[84 rows x 6 columns]
data_keps=            BJD        RV       eRV  Flag    eRV_og  residuals
0      0.000000  4.912826  2.541371     1  0.640157   1.624782
1     59.930089  4.434408  2.668277     1  1.034862   2.361263
2    108.761071  3.277795  2.481395     1  0.329479   0.254507
3    397.936897 -0.775525  2.511916     1  0.510838   0.506437
4    397.942036 -0.538872  2.511765     1  0.510097   0.744319
..          ...       ...       ...   ...       ...        ...
31  6662.795638 -6.743544  3.102287     2  1.319130  -3.395605
32  6662.800288 -6.294784  3.115435     2  1.349760  -2.946442
33  6711.647318  2.666146  3.075738     2  1.255420   4.299660
34  6711.649908 -1.291774  3.072120     2  1.246530   0.342292
35  6711.652458  1.802906  3.071507     2  1.245020   3.437516

[84 rows x 6 columns]
data1=            BJD        RV       eRV  Flag    eRV_og  residuals
0      0.000000  3.799448  2.541371     1  0.640157   1.624782
0     36.899961  2.051683  3.070883     2  1.243480   0.113863
1     59.930089  4.144167  2.668277     1  1.034862   2.361263
5    101.814282  3.149020  2.350560     3  0.154521   1.660798
6    104.779345  2.741009  2.436845     3  0.661031   1.274247
..          ...       ...       ...   ...       ...        ...
0   2512.025052 -0.313657  2.349396     3  0.135676  -3.041615
1   2512.989258 -0.308165  2.347816     3  0.104818  -3.031246
2   2514.090094 -0.057631  2.374180     3  0.368069  -2.775130
3   2516.015901  2.681919  2.350716     3  0.156883  -0.025775
4   2517.039496  3.816059  2.349484     3  0.137190   1.113596

[84 rows x 6 columns]
data1=          BJD        RV       eRV  Flag    eRV_og  residuals
0    0.000000  2.738161  2.541371     1  0.640157   1.624782
22   0.189767 -1.989284  2.482531     1  0.337927  -3.068220
9    1.175086 -0.448357  2.471904     1  0.248078  -1.341767
17   1.251406  6.332444  2.348890     3  0.126618   5.453852
18   2.228863  1.526952  2.350846     3  0.158821   0.843385
..        ...       ...       ...   ...       ...        ...
2   54.961181  0.706180  3.008409     2  1.080020  -0.562030
15  55.021128  1.349096  2.481958     1  0.333691   0.091052
8   55.038571  3.274913  2.467368     1  0.197846   2.019836
14  55.043484  3.054355  2.477114     1  0.295517   1.800115
16  55.072281  1.580842  2.350499     3  0.153597   0.331510

[84 rows x 6 columns]


                          Plotting E[log L](beta) Plot

100%|██████████| 1/1 [00:00<00:00,  7.79it/s]


                             Plotting Beta Density

100%|██████████| 1/1 [00:00<00:00,  4.77it/s]


                           Plotting Temperature Rates

100%|██████████| 7/7 [00:00<00:00, 11.37it/s]


                           Plotting Gaussian Mixtures

100%|██████████| 20/20 [00:04<00:00,  4.50it/s]




Time Table
Time RUN                   : 00:00:13
Time POSTPROCESS           : 00:00:01
Time CALCULATE GM          : 00:00:00
Time plot_posteriors       : 00:00:02
Time plot_histograms       : 00:00:00
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
Time plot_astrometric_model: 00:00:00

(Past BIC) - (Present BIC) > 5

 BIC condition met!!
inf - 457.186 > 5


                                                                                
                              ~~ End of the Run ~~