search

j-faria / kima

Exoplanet detection in RVs with DNest4 and GPs
http://www.kima.science
MIT License
15 stars 10 forks source link

Patches for GP beta sampling and KimaResults column counting #72

Closed accameron closed 2 years ago

accameron commented 2 years ago

Thank you for sending a pull request! :tada:

Proposed changes

Describe your changes here and tell the maintainers why we should accept this pull request:

  1. When the GP option is used with activity decorrelation, sampling from the beta prior was not being carried out. To remedy this I have replaced lines 546-586 of kima/src/RVmodel.cpp with a copy of lines 621-681.
  2. When the GP option is used with activity decorrelation, columns are counted incorrectly in pykima/results.py with the result that results.etas contains values read from the columns of posterior_sample.txt that actually contain betas. To remedy this I have added n_act_ind to the calculation of start_hyperpars at line 276.

Types of changes

What types of changes does your code introduce in kima?

Tests

Did you test the new code locally and/or added new tests to ensure everything is working correctly?
I tested the code locally with the attached data file and kima_setup.cpp containing:

RVmodel::RVmodel():fix(false),npmax(5) { auto data = get_data(); //kernel = celerite;

// Priors for orbital elements
auto conditional = planets.get_conditional_prior();
conditional->Pprior = make_prior<LogUniform>(0.5, 25); // days
conditional->Kprior = make_prior<ModifiedLogUniform>(0.5, 15.); // m/s
conditional->eprior = make_prior<Kumaraswamy>(0.867, 3.03);
conditional->phiprior = make_shared<Uniform>(0, 2*PI);
conditional->wprior = make_shared<Uniform>(0, 2*PI);

//Prior on systemic velocity in m/s
Cprior = make_prior<Uniform>(data.get_RV_min(), data.get_RV_max());

// use the default priors
// priors for first known_object
KO_Pprior[0] = make_prior<Gaussian>(0.85359159, 5.7e-7);
KO_Kprior[0] = make_prior<LogUniform>(1, 10);
KO_eprior[0] = make_prior<Kumaraswamy>(0.867, 3.03);
// KO_eprior[0] = make_prior<Uniform>(0, .1);
KO_wprior[0] = make_prior<Uniform>(-PI, PI);
// KO_phiprior[0] = make_prior<Uniform>(0, 2*PI);
KO_phiprior[0] = make_prior<Gaussian_from_Tc>(54398.07756, 0.0007, 0.85359159, 5.7e-7);

//Priors for additional white noise
Jprior = make_prior<ModifiedLogUniform>(0.5, data.get_max_RV_span());

// Priors for GP hyperparameters (Tailored for CoRoT-7)
log_eta1_prior = make_prior<Uniform>(1.6, 2.7);
eta2_prior = make_prior<LogUniform>(15, 60);
eta3_prior = make_prior<Uniform>(19, 25);
log_eta4_prior = make_prior<Uniform>(-1, -.35);

//Priors for correlation coefficients with activity indicators
betaprior = make_prior<Gaussian>(0,1);

}

int main(int argc, char** argv) { datafile = "corot7activ.txt"; indicators = {"U0", "U1", "U2", "U3", "U4"}; load(datafile, "ms", 4, indicators);

Sampler<RVmodel> sampler = setup<RVmodel>(argc, argv);
sampler.run(50);

return 0;

}

Did you add/change the documentation? No data_activ.txt o

j-faria commented 2 years ago

Thank you very much, Andrew!