diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a2ddc3..82317c0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -73,6 +73,10 @@ nanobind_add_module(SPLEAFmodel STABLE_ABI NB_STATIC src/kima/SPLEAFmodel.cpp src/kima/Data.cpp src/kima/ConditionalPrior.cpp src/kima/spleaf.cpp src/kima/kepler.cpp src/kima/AMDstability.cpp) +nanobind_add_module(BINARIESmodel STABLE_ABI NB_STATIC + src/kima/BINARIESmodel.cpp src/kima/Data.cpp src/kima/ConditionalPrior.cpp + src/kima/postkepler.cpp src/kima/kepler.cpp src/kima/AMDstability.cpp) + # compile the Sampler module to create the kima.run function nanobind_add_module(Sampler STABLE_ABI NB_STATIC @@ -83,8 +87,10 @@ nanobind_add_module(Sampler STABLE_ABI NB_STATIC src/kima/GPmodel.cpp src/kima/RVFWHMmodel.cpp src/kima/SPLEAFmodel.cpp + src/kima/BINARIESmodel.cpp src/kima/Data.cpp src/kima/ConditionalPrior.cpp - src/kima/spleaf.cpp src/kima/kepler.cpp src/kima/transits.cpp src/kima/AMDstability.cpp) + src/kima/spleaf.cpp src/kima/kepler.cpp + src/kima/postkepler.cpp src/kima/transits.cpp src/kima/AMDstability.cpp) nanobind_add_module(kepler STABLE_ABI NB_STATIC @@ -100,7 +106,7 @@ target_compile_features(transits PRIVATE cxx_std_17) set(INCLUDES src/vendor/cpp-loadtxt/src ${DNEST4_PATH}) # compilation settings for each module -foreach(targ distributions Data ConditionalPrior RVmodel OutlierRVmodel TRANSITmodel) +foreach(targ distributions Data ConditionalPrior RVmodel OutlierRVmodel TRANSITmodel BINARIESmodel) target_compile_features(${targ} PRIVATE cxx_std_17) target_include_directories(${targ} PRIVATE ${INCLUDES}) target_link_libraries(${targ} PRIVATE dnest4) @@ -152,7 +158,7 @@ foreach(module Data ConditionalPrior) install(TARGETS ${module} LIBRARY DESTINATION kima) endforeach() -foreach(module RVmodel OutlierRVmodel TRANSITmodel GPmodel RVFWHMmodel SPLEAFmodel) +foreach(module RVmodel OutlierRVmodel TRANSITmodel GPmodel RVFWHMmodel SPLEAFmodel BINARIESmodel) install(TARGETS ${module} LIBRARY DESTINATION kima) endforeach() diff --git a/src/kima/BINARIESmodel.cpp b/src/kima/BINARIESmodel.cpp index 724175c..0a6fb3a 100644 --- a/src/kima/BINARIESmodel.cpp +++ b/src/kima/BINARIESmodel.cpp @@ -1,7 +1,7 @@ -#include "RV_binaries_model.h" +#include "BINARIESmodel.h" using namespace std; -using namespace Eigen; +// using namespace Eigen; using namespace DNest4; using namespace nijenhuis; using namespace postKep; @@ -23,7 +23,6 @@ void BINARIESmodel::initialize_from_data(RVData& data) // resize RV model vector mu.resize(data.N()); - mu_2.resize(data.N()); // set default conditional priors that depend on data @@ -61,7 +60,7 @@ void BINARIESmodel::setPriors() // BUG: should be done by only one thread! } // if offsets_prior is not (re)defined, assume a default - if (multi_instrument && !offsets_prior) + if (data._multi && !offsets_prior) offsets_prior = make_prior( -data.get_RV_span(), data.get_RV_span() ); for (size_t j = 0; j < data.number_instruments - 1; j++) @@ -76,7 +75,7 @@ void BINARIESmodel::setPriors() // BUG: should be done by only one thread! // if (n_known_object == 0) cout << "Warning: `known_object` is true, but `n_known_object` is set to 0"; for (int i = 0; i < n_known_object; i++){ if (!KO_Pprior[i] || !KO_Kprior[i] || !KO_eprior[i] || !KO_phiprior[i] || !KO_wprior[i] || !KO_wdotprior[i]) - throw std::logic_error("When known_object=true, please set priors for each (KO_Pprior, KO_Kprior, KO_eprior, KO_phiprior, KO_wprior, KO_wdotprior)"); + throw std::logic_error("When using BINARIESmodel, please set priors on the binary parameters (KO_Pprior, KO_Kprior, KO_eprior, KO_phiprior, KO_wprior, KO_wdotprior)"); if (double_lined && !KO_qprior[i]) throw std::logic_error("When double_lined=true, please set prior for KO_qprior"); } @@ -122,11 +121,11 @@ void BINARIESmodel::from_prior(RNG& rng) } - if (data.indicator_correlations) - { - for (unsigned i=0; igenerate(rng); - } +// if (data.indicator_correlations) +// { +// for (unsigned i=0; igenerate(rng); +// } if (known_object) { // KO mode! KO_P.resize(n_known_object); @@ -238,29 +237,28 @@ void BINARIESmodel::calculate_mu() //else //P = components[j][0]; - P = components[j][0] + P = components[j][0]; K = components[j][1]; phi = components[j][2]; ecc = components[j][3]; omega = components[j][4]; - omegadot = components[j][5]; - auto v = brandt::keplerian(data.t, P, K, ecc, omega, phi, data.M0_epoch); - for(size_t i=0; iperturb(jitters[i], rng); @@ -566,17 +563,17 @@ double BINARIESmodel::perturb(RNG& rng) if(trend) { mu[i] -= slope*(data.t[i]-tmid) + quadr*pow(data.t[i]-tmid, 2) + cubic*pow(data.t[i]-tmid, 3); } - if(data._multi)) { + if(data._multi) { for(size_t j=0; jperturb(bkg2, rng); // propose new instrument offsets - if (data._multi)){ + if (data._multi){ for(unsigned j=0; jperturb(offsets[j], rng); } @@ -622,11 +619,11 @@ double BINARIESmodel::perturb(RNG& rng) } // propose new indicator correlations - if(data.indicator_correlations){ - for(size_t j = 0; j < data.number_indicators; j++){ - betaprior->perturb(betas[j], rng); - } - } +// if(data.indicator_correlations){ +// for(size_t j = 0; j < data.number_indicators; j++){ +// betaprior->perturb(betas[j], rng); +// } +// } for(size_t i=0; ieprior << endl; fout << "phiprior: " << *conditional->phiprior << endl; fout << "wprior: " << *conditional->wprior << endl; - fout << "wdotprior: " << *conditional->wdotprior << endl; } if (known_object) { @@ -1040,7 +1036,7 @@ auto BINARIESMODEL_DOC = R"D( Implements a sum-of-Keplerians model where the number of Keplerians can be free. This model assumes white, uncorrelated noise. This modules is tailored for the analysis of stellar binaries through the known object mode (without it this defaults to the RVmodel) -The binary can have one set of RVs or two (one on each star) and the model adds apsidal +The binary can have one set of RVs or two (one on each star). This model adds apsidal precession as a free parameter and accounts for GR and Tidal effects on the radial velocities. Args: @@ -1049,22 +1045,22 @@ precession as a free parameter and accounts for GR and Tidal effects on the radi data (RVData): the RV data )D"; -class RVmodel_publicist : public RVmodel +class BINARIESmodel_publicist : public BINARIESmodel { public: - using RVmodel::trend; - using RVmodel::degree; - using RVmodel::studentt; - using RVmodel::known_object; - using RVmodel::n_known_object; - using RVmodel::star_mass; - using RVmodel::binary_mass; - using RVmodel::star_radius; - using RVmodel::binary_radius; - using RVmodel::enforce_stability; - using RVmodel::relativistic_correction; - using RVmodel::tidal_correction; - using RVmodel::double_lined; + using BINARIESmodel::trend; + using BINARIESmodel::degree; + using BINARIESmodel::studentt; + using BINARIESmodel::known_object; + using BINARIESmodel::n_known_object; + using BINARIESmodel::star_mass; + using BINARIESmodel::binary_mass; + using BINARIESmodel::star_radius; + using BINARIESmodel::binary_radius; + using BINARIESmodel::enforce_stability; + using BINARIESmodel::relativistic_correction; + using BINARIESmodel::tidal_correction; + using BINARIESmodel::double_lined; }; @@ -1103,7 +1099,7 @@ NB_MODULE(BINARIESmodel, m) { // .def_rw("relativistic_correction", &BINARIESmodel_publicist::relativistic_correction, "whether to perform the GR correction") - .def_rw("tidal_correction", &BINARIESmodel_publicist::tidal_correection, + .def_rw("tidal_correction", &BINARIESmodel_publicist::tidal_correction, "whether to perform the tidal correction") // @@ -1139,6 +1135,33 @@ NB_MODULE(BINARIESmodel, m) { [](BINARIESmodel &m) { return m.nu_prior; }, [](BINARIESmodel &m, distribution &d) { m.nu_prior = d; }, "Prior for the degrees of freedom of the Student-t likelihood") + + // known object priors + // ? should these setters check if known_object/sb2 is true? + .def_prop_rw("KO_Pprior", + [](BINARIESmodel &m) { return m.KO_Pprior; }, + [](BINARIESmodel &m, std::vector& vd) { m.KO_Pprior = vd; }) + .def_prop_rw("KO_Kprior", + [](BINARIESmodel &m) { return m.KO_Kprior; }, + [](BINARIESmodel &m, std::vector& vd) { m.KO_Kprior = vd; }) + .def_prop_rw("KO_eprior", + [](BINARIESmodel &m) { return m.KO_eprior; }, + [](BINARIESmodel &m, std::vector& vd) { m.KO_eprior = vd; }) + .def_prop_rw("KO_wprior", + [](BINARIESmodel &m) { return m.KO_wprior; }, + [](BINARIESmodel &m, std::vector& vd) { m.KO_wprior = vd; }) + .def_prop_rw("KO_phiprior", + [](BINARIESmodel &m) { return m.KO_phiprior; }, + [](BINARIESmodel &m, std::vector& vd) { m.KO_phiprior = vd; }) + .def_prop_rw("KO_wdotprior", + [](BINARIESmodel &m) { return m.KO_wdotprior; }, + [](BINARIESmodel &m, std::vector& vd) { m.KO_wdotprior = vd; }) + .def_prop_rw("KO_qprior", + [](BINARIESmodel &m) { return m.KO_qprior; }, + [](BINARIESmodel &m, std::vector& vd) { m.KO_qprior = vd; }) + + + // conditional object .def_prop_rw("conditional", [](BINARIESmodel &m) { return m.get_conditional_prior(); }, diff --git a/src/kima/BINARIESmodel.h b/src/kima/BINARIESmodel.h index c9f691d..dcb061a 100644 --- a/src/kima/BINARIESmodel.h +++ b/src/kima/BINARIESmodel.h @@ -7,6 +7,7 @@ #include "ConditionalPrior.h" #include "utils.h" #include "kepler.h" +#include "postkepler.h" #include "AMDstability.h" using namespace std; @@ -20,7 +21,8 @@ using namespace nb::literals; #include "nb_shared.h" -class BINARIES_model + +class KIMA_API BINARIESmodel { protected: /// whether the model includes a polynomial trend @@ -99,9 +101,9 @@ class BINARIES_model std::vector KO_wdot; // The signal - std::vector mu = // the RV model + std::vector mu;// the RV model //std::vector(RVData::get_instance().N()); - std::vector mu_2 = // the RV model for secondary + std::vector mu_2;// the RV model for secondary //std::vector(RVData::get_instance().N()); // changed to imitate RVFWHM get_data replaced by get_instance void calculate_mu(); void calculate_mu_2(); @@ -139,7 +141,7 @@ class BINARIES_model distribution cubic_prior; /// (Common) prior for the between-instruments offsets. distribution offsets_prior; - std::vector individual_offset_prior + std::vector individual_offset_prior ; // priors for KO mode! diff --git a/src/kima/Data.cpp b/src/kima/Data.cpp index 3cf1ed2..7aa98f7 100644 --- a/src/kima/Data.cpp +++ b/src/kima/Data.cpp @@ -329,7 +329,6 @@ for (size_t n = 0; n < t.size(); n++) { y[n] = y[n] * factor; sig[n] = sig[n] * factor; - sig[n] = sig[n] * factor; if (sb2) { y2[n] = y2[n] * factor; diff --git a/src/kima/Data.h b/src/kima/Data.h index ae2e187..89b0fa5 100644 --- a/src/kima/Data.h +++ b/src/kima/Data.h @@ -46,7 +46,7 @@ vector sort_indexes(const vector &v) { } -class RVData { +class KIMA_API RVData { friend class RVmodel; friend class GPmodel; diff --git a/src/kima/__init__.py b/src/kima/__init__.py index 1b08e4e..9d17d8c 100644 --- a/src/kima/__init__.py +++ b/src/kima/__init__.py @@ -6,6 +6,7 @@ from .SPLEAFmodel import SPLEAFmodel from .TRANSITmodel import TRANSITmodel from .OutlierRVmodel import OutlierRVmodel +from .BINARIESmodel import BINARIESmodel from .Sampler import run diff --git a/src/kima/examples/Kepler16/Kepler16.rdb b/src/kima/examples/Kepler16/Kepler16.rdb new file mode 100755 index 0000000..17b9f08 --- /dev/null +++ b/src/kima/examples/Kepler16/Kepler16.rdb @@ -0,0 +1,142 @@ +bjd vrad svrad fwhm contrast bis_span berv sn6 sn26 sn35 +--- ---- ----- ---- ---- ---- ---- ---- ---- ---- +57578.42404 -22.7567 0.0126 10.4138 22.7274 -0.0233 3.8255 9.9 38.6 49.1 +57582.52743 -20.4322 0.0176 10.3201 27.222 -0.0199 3.1132 2.9 22.8 30.4 +57587.43166 -23.9495 0.0113 10.2935 25.8863 0.0192 2.5266 7.9 37.5 47.7 +57591.44122 -29.2376 0.0146 10.2883 26.1452 0.0538 1.9226 4.8 28.5 36.4 +57595.47979 -35.3428 0.0096 10.334 29.0424 0.0125 1.2659 6.4 39.5 50.2 +57604.43905 -46.8841 0.0145 10.4605 25.3376 -0.0054 -0.0596 5.6 29.9 38.3 +57607.38886 -47.6294 0.0145 10.4123 28.2874 0.0037 -0.4504 3.5 26.8 34.6 +57612.47728 -40.2406 0.011 10.3402 28.9602 0.0079 -1.3426 5.2 34.7 44.5 +57623.37736 -20.3825 0.01 10.2928 29.0664 -0.0093 -2.8515 6.2 37.8 48.6 +57626.38491 -21.8607 0.011 10.3395 28.928 0.0329 -3.3029 5 34.3 44.4 +57634.40704 -32.0631 0.0109 10.3925 28.9876 -0.067 -4.4706 5.3 35 45 +57719.26505 -36.2013 0.0095 10.2807 28.7526 0.0100 -8.5613 10.2 40.2 50.1 +57746.26923 -20.3670 0.0085 10.3421 28.9562 -0.0183 -6.0605 11.1 45.1 57.6 +57815.65363 -44.9323 0.0100 10.2523 26.6492 0.0043 4.6822 11.8 40.9 50.4 +57815.67731 -44.8881 0.0083 10.2621 28.9193 0.0079 4.6710 11.4 46.0 56.8 +57850.59947 -46.6705 0.0092 10.2708 26.6208 0.0055 8.2528 12.5 44.3 55.0 +57850.62308 -46.6764 0.0120 10.2166 27.8034 -0.0089 8.2335 8.2 32.7 40.7 +57858.58652 -41.1541 0.0150 10.1911 23.2800 -0.0157 8.6830 10.0 31.0 38.1 +57860.61247 -35.6321 0.0192 10.1749 22.7570 0.0001 8.7402 8.0 24.9 30.6 +57881.49587 -32.9924 0.0182 10.2268 21.4730 -0.0099 9.0574 9.7 27.8 34.5 +57881.51944 -33.0022 0.0160 10.2041 25.4235 -0.0059 9.0393 7.5 26.7 33.3 +57890.53416 -45.6721 0.0090 10.2166 26.9172 0.0145 8.7767 12.9 44.7 54.6 +57890.57816 -45.7175 0.0082 10.1824 27.3311 0.0022 8.7247 14.1 48.4 58.8 +57914.55102 -22.4007 0.0112 10.2007 24.8603 -0.0160 7.1242 12.5 38.8 47.8 +57914.57472 -22.4184 0.0085 10.2088 27.7574 0.0047 7.0911 13.4 46.0 56.6 +57924.49590 -35.9622 0.0101 10.2442 26.7931 0.0275 6.1723 11.4 40.2 50.0 +57924.52118 -35.9972 0.0103 10.2056 28.6725 0.0157 6.1373 9.3 36.7 46.0 +57936.48767 -47.4512 0.0108 10.2167 26.3998 0.0047 4.7140 11.0 38.0 46.9 +57936.51133 -47.4445 0.0113 10.2790 28.3570 0.0015 4.6802 8.5 33.9 42.1 +57954.37721 -21.3710 0.0137 10.2775 24.7272 -0.0031 2.3473 9.9 32.2 39.9 +57955.36873 -22.1571 0.0132 10.2851 24.6882 -0.0219 2.2090 10.4 33.5 42.4 +57970.44413 -43.0795 0.0086 10.1809 27.9862 0.0258 -0.1818 12.5 45.2 55.7 +57970.46806 -43.1248 0.0082 10.1928 28.7579 0.0006 -0.2158 12.4 46.1 56.6 +57989.42736 -22.4538 0.0099 10.2214 26.7263 0.0009 -3.0409 12.0 41.1 51.5 +57989.45189 -22.4256 0.0091 10.2462 28.9240 -0.0128 -3.0731 10.7 41.3 51.7 +57999.37737 -25.3649 0.0165 10.2100 23.2439 0.0027 -4.3921 9.4 28.3 36.6 +58003.42095 -31.0559 0.0146 10.1719 22.5221 -0.0485 -4.9856 11.5 33.0 41.0 +58007.42352 -37.1454 0.0142 10.2436 28.1567 -0.0170 -5.5047 6.5 27.2 34.3 +58026.32236 -31.5715 0.0102 10.2507 26.1549 0.0080 -7.4837 12.3 40.9 51.1 +58029.34747 -24.2833 0.0087 10.2234 30.6550 -0.0043 -7.7808 11.4 40.6 50.6 +58034.33565 -20.3785 0.0117 10.2511 30.4850 -0.0444 -8.1673 9.0 30.4 38.8 +58038.28757 -22.8841 0.0086 10.1955 28.8814 -0.0074 -8.3917 12.0 43.8 55.1 +58043.30107 -29.2523 0.0089 10.1540 27.1799 0.0023 -8.7098 13.6 44.8 56.0 +58049.33656 -38.4059 0.0195 10.2348 22.8639 0.0170 -9.0165 7.6 24.4 30.6 +58057.38946 -47.3980 0.0127 10.2675 27.4942 -0.0127 -9.2448 6.2 31.3 40.0 +58066.32858 -34.6719 0.0178 10.2258 23.2093 -0.0015 -9.2345 8.0 26.3 33.1 +58076.24290 -20.6152 0.0112 10.2370 25.8315 0.0194 -8.9394 11.4 37.8 47.6 +58084.24172 -29.0392 0.0424 10.2126 16.3483 0.0261 -8.5534 6.0 15.7 19.5 +58227.56803 -42.8022 0.0095 10.2644 29.0967 0.0023 8.8458 8.7 39.3 49.5 +58231.59284 -31.9117 0.0116 10.2791 28.0351 -0.0215 8.9311 8.4 33.5 42.3 +58255.55341 -39.6711 0.0113 10.2359 26.2808 -0.0078 8.7672 10.2 36.7 46.3 +58271.53928 -35.1047 0.0817 10.2827 13.3499 0.0391 7.8481 3.7 10.0 12.2 +58289.48662 -28.8532 0.0087 10.2042 27.2210 0.0031 6.2062 13.3 45.9 56.9 +58294.41403 -36.2925 0.0248 10.0754 18.8088 -0.0413 5.7173 8.7 23.1 27.8 +58301.52589 -45.8755 0.0095 10.1888 26.5342 -0.0116 4.6984 12.8 43.0 52.5 +58319.49236 -21.0031 0.0096 10.2279 26.9129 -0.0296 2.2166 12.3 41.9 52.4 +58329.35035 -27.1402 0.0124 10.1420 26.6103 -0.0095 0.9139 6.6 32.9 42.1 +58364.49871 -20.9600 0.0092 10.2883 30.7118 -0.0072 -4.4962 9.0 38.4 48.9 +58373.38714 -31.4314 0.0085 10.2044 28.8165 -0.0127 -5.5610 8.1 44.1 55.9 +58389.28389 -46.7867 0.0085 10.2039 28.3455 -0.0011 -7.2120 8.6 45.0 56.5 +58410.34703 -25.6062 0.0108 10.2338 30.3847 -0.0018 -8.8412 4.7 33.1 43.4 +58414.31942 -31.2263 0.0136 10.1724 26.2082 -0.0249 -8.9835 5.8 30.5 39.2 +58438.23468 -28.4894 0.0083 10.2378 30.2861 -0.0021 -9.0298 8.2 43.3 55.3 +58440.32791 -23.9451 0.0142 10.4297 30.4116 -0.0077 -9.0345 2.3 25.4 34.8 +58447.28753 -21.3585 0.0086 10.2367 30.3292 -0.0076 -8.7093 7.5 41.8 54.3 +58536.69186 -29.9230 0.0132 10.1865 26.3805 -0.0380 3.3226 5.6 31.1 40.5 +58539.69622 -34.4566 0.0101 10.2838 27.8847 -0.0061 3.7535 6.3 38.7 49.6 +58542.70128 -39.0670 0.0096 10.1926 28.3573 -0.0346 4.1735 6.6 40.0 51.1 +58557.68116 -39.1647 0.0161 10.1745 27.1729 0.0440 6.0936 3.7 24.9 32.1 +58569.61968 -20.8099 0.0122 10.2121 27.6462 -0.0050 7.3605 6.5 32.1 41.9 +58617.59111 -28.1344 0.0091 10.1989 28.6512 0.0008 8.8316 7.6 41.7 52.9 +58626.46789 -41.4267 0.0109 10.2403 25.2701 -0.0018 8.5884 8.1 39.4 50.3 +58634.53897 -47.5612 0.0094 10.1840 25.9739 0.0029 8.0204 9.7 44.3 55.0 +58638.43325 -42.6217 0.0126 10.2257 28.3202 -0.0451 7.8314 3.9 30.4 39.6 +58650.53936 -20.4183 0.0153 10.1425 26.2074 -0.0095 6.5953 4.7 27.0 35.2 +58654.49571 -22.9447 0.0097 10.2137 28.2516 -0.0082 6.2263 7.4 39.7 50.5 +58660.56745 -30.8940 0.0103 10.2106 28.7369 -0.0298 5.4296 5.8 36.7 46.8 +58665.57318 -38.5132 0.0141 10.2621 27.1916 0.0497 4.7959 4.3 28.3 36.5 +58675.57952 -47.5856 0.0110 10.2959 27.3243 -0.0095 3.4371 6.3 36.3 45.8 +58687.61864 -22.6083 0.0228 10.3486 23.4485 0.0170 1.6671 3.0 20.4 26.9 +58703.54677 -33.7798 0.0178 10.2717 25.5443 0.0157 -0.6948 3.6 23.9 30.8 +58706.50803 -38.3152 0.0101 10.1879 27.8094 -0.0183 -1.1089 7.2 38.8 48.4 +58708.50232 -41.2225 0.0326 10.3182 14.8637 -0.0525 -1.4057 7.5 22.6 28.5 +58732.43751 -20.3784 0.0302 10.0867 25.0051 -0.1717 -4.8088 1.8 14.3 19.8 +58734.44238 -21.0856 0.0145 10.2832 30.1981 -0.0537 -5.0776 3.6 24.9 33.5 +58738.38163 -24.8706 0.0106 10.1620 27.5478 -0.0376 -5.5150 7.2 37.3 47.6 +58753.33846 -45.8311 0.0076 10.2117 28.5105 0.0012 -7.1665 10.1 50.5 62.2 +58757.43639 -47.6342 0.0099 10.2701 28.2921 -0.0123 -7.6487 5.9 38.8 49.4 +58760.37181 -45.0732 0.0118 10.1831 28.0912 0.0018 -7.8524 5.4 32.9 42.2 +58765.41646 -32.4324 0.0095 10.2732 28.0812 -0.0013 -8.2788 7.1 41.0 51.5 +58794.31837 -45.7499 0.0157 10.2562 26.6709 -0.0224 -9.2414 3.8 26.0 34.3 +58804.24157 -38.8893 0.0122 10.2536 27.8082 0.0145 -9.0155 5.1 32.0 41.6 +58820.22529 -24.4742 0.0150 10.3222 26.3384 -0.0164 -8.1824 4.3 27.6 36.8 +58824.24004 -29.9372 0.0095 10.2347 28.0045 0.0142 -7.8812 7.0 40.9 52.9 +58828.27635 -36.0603 0.0201 10.3544 22.8020 -0.0434 -7.5453 3.4 23.9 32.4 +58919.69124 -47.3559 0.0209 10.2482 24.3342 -0.1464 5.7019 3.7 21.4 27.7 +58998.56140 -44.4953 0.0109 10.2477 28.7207 0.0021 8.0760 5.2 34.7 44.4 +59009.49993 -39.2280 0.0123 10.2154 27.5204 -0.0244 7.2684 5.1 32.2 41.6 +59023.55894 -22.2996 0.0100 10.2549 28.4154 -0.0127 5.7098 6.3 38.5 49.6 +59030.43413 -31.1013 0.0096 10.1992 28.2237 0.0220 5.0181 7.3 40.1 50.2 +59042.42642 -47.0560 0.0101 10.2237 28.0778 0.0035 3.4116 6.9 38.3 48.7 +59046.54366 -46.7885 0.0124 10.2485 27.8792 -0.0738 2.6703 4.9 31.6 39.9 +59067.42681 -25.4094 0.0097 10.1999 27.8897 -0.0028 -0.3406 7.9 40.0 50.7 +59073.55426 -34.2350 0.0094 10.3107 28.2457 0.0142 -1.4159 6.4 41.3 52.4 +59077.48604 -40.2123 0.0157 10.2125 27.3191 0.0136 -1.9512 3.6 25.4 32.6 +59089.39809 -44.3674 0.0100 10.2499 27.9722 -0.0172 -3.6094 8.6 39.1 49.2 +59093.41074 -34.2917 0.0153 10.1700 25.7562 -0.0138 -4.1892 5.2 27.6 35.4 +59095.45922 -28.5565 0.0099 10.1430 27.1355 -0.0045 -4.5222 8.7 40.2 50.7 +59101.49898 -20.3950 0.0089 10.2737 28.2692 -0.0058 -5.3490 7.4 43.2 55.0 +59102.51980 -20.4161 0.0123 10.5199 26.8784 0.0042 -5.4884 3.5 33.4 44.3 +59121.39686 -44.0369 0.0144 10.1553 24.2377 0.0186 -7.4937 6.6 31.0 39.0 +59123.35979 -46.1338 0.0094 10.2047 27.1289 0.0037 -7.6363 8.8 42.8 53.7 +59131.39916 -42.5203 0.0104 10.2353 28.1645 0.0077 -8.3179 5.7 37.2 47.3 +59137.28463 -26.6917 0.0092 10.2507 28.6035 0.0215 -8.5877 8.0 41.5 52.6 +59154.32126 -32.1040 0.0103 10.2522 26.5496 -0.0094 -9.2178 7.9 40.0 51.5 +59157.36523 -36.7756 0.0179 10.1740 24.8354 0.0126 -9.2668 4.1 24.4 31.4 +59162.32683 -43.8055 0.0354 10.5774 23.7914 -0.0130 -9.2293 1.4 13.1 17.5 +59165.27801 -46.7974 0.0093 10.2647 28.2102 0.0230 -9.1519 6.4 41.5 52.6 +59172.32740 -42.8404 0.0082 10.4123 28.5351 -0.0048 -8.9800 4.7 47.0 59.9 +59181.27395 -21.7233 0.0172 10.3619 30.0189 -0.0835 -8.4970 3.3 21.2 28.2 +59266.69524 -20.3815 0.0109 10.2969 28.2069 0.0189 3.2546 5.6 35.3 45.7 +59269.69460 -22.0007 0.0110 10.2942 27.8791 0.0121 3.6879 5.7 35.6 46.0 +59270.67941 -22.9327 0.0119 10.2624 26.7721 0.0342 3.8373 6.0 34.1 43.7 +59275.70395 -29.3429 0.0134 10.2014 25.7622 -0.0227 4.5135 5.7 31.4 40.4 +59280.69715 -36.9394 0.0146 10.1737 27.2043 0.0328 5.1752 4.2 27.3 35.3 +59297.61346 -37.6939 0.0090 10.2603 28.2994 0.0369 7.1366 7.2 42.6 54.0 +59299.64512 -31.8353 0.0095 10.2965 27.9143 -0.0388 7.2998 7.6 41.3 52.2 +59303.60959 -22.8534 0.0097 10.2323 27.0716 0.0388 7.6687 8.1 41.6 52.8 +59305.61783 -20.8237 0.0089 10.2741 27.9176 0.0060 7.8221 8.5 44.1 55.5 +59349.50200 -20.5593 0.0264 10.4323 25.2883 0.1036 8.8882 2.2 16.5 21.8 +59354.55499 -24.8712 0.0143 10.2037 28.1379 0.0046 8.6258 4.0 26.9 34.9 +59363.59631 -38.0869 0.0097 10.2319 27.7710 0.0159 8.0491 7.6 40.2 50.0 +59366.52159 -42.3085 0.0116 10.2459 27.8744 0.0331 7.9347 5.5 33.6 42.8 +59369.57437 -45.8848 0.0086 10.2539 28.7213 -0.0358 7.6310 7.8 43.9 54.6 +59371.53106 -47.3335 0.0106 10.2972 28.4898 -0.0401 7.5227 5.7 36.0 45.1 +59375.55836 -46.3575 0.0089 10.2317 28.8012 0.0087 7.1249 7.5 42.6 52.8 +59378.56443 -40.8905 0.0100 10.2850 28.8038 -0.0081 6.8238 6.1 37.8 47.8 +59382.54493 -29.7361 0.0129 10.2204 28.1204 0.0375 6.4345 4.5 30.0 38.4 +59387.49783 -21.0061 0.0124 10.2410 26.6086 0.0101 5.9421 6.1 33.0 42.3 +59388.44241 -20.4779 0.0134 10.1903 26.5662 0.0071 5.8984 5.5 30.6 39.2 \ No newline at end of file diff --git a/src/kima/examples/Kepler16/__init__.py b/src/kima/examples/Kepler16/__init__.py new file mode 100644 index 0000000..588a07e --- /dev/null +++ b/src/kima/examples/Kepler16/__init__.py @@ -0,0 +1,56 @@ +import os +import numpy as np +import kima +from kima import RVData, RVmodel, BINARIESmodel +from kima.pykima.utils import chdir +from kima.distributions import Uniform, Gaussian, ModifiedLogUniform + +__all__ = ['Kepler16'] + +here = os.path.dirname(__file__) # cwd + +def Kepler16(run=False, **kwargs): + """ + Create (and optionally run) an RV model for analysis of 51 Peg data. + This loads Keck/HIRES data from `51Peg.rv` and creates a model where + the number of Keplerians is free from 0 to 1. + + Args: + run (bool): whether to run the model + **kwargs: keyword arguments passed directly to `kima.run` + """ + # load the right data file + data = RVData([os.path.join(here, 'Kepler16.rdb')],skip=2,units='kms') + # create the model + model = BINARIESmodel(fix=False, npmax=0, data=data) + + model.Cprior = Uniform(-44300,-23300) + model.Jprior = ModifiedLogUniform(1.0,100) + + # model.set_known_object = 1 + # model.n_known_object = 1 + model.KO_Pprior = [Gaussian(41,1)] + model.KO_Kprior = [Gaussian(13600,900)] + model.KO_eprior = [Gaussian(0.16,0.02)] + model.KO_wprior = [Uniform(0,2*np.pi)] + model.KO_wdotprior = [Gaussian(0,0.1)] + model.KO_phiprior = [Uniform(0,2*np.pi)] + + kwargs.setdefault('steps', 5000) + kwargs.setdefault('num_threads', 4) + kwargs.setdefault('num_particles', 2) + kwargs.setdefault('new_level_interval', 5000) + kwargs.setdefault('save_interval', 1000) + + if run: + kima.run(model, **kwargs) + # with chdir(here): + # kima.run(model, **kwargs) + + return model + +if __name__ == '__main__': + model = Kepler16(run=True, steps=2000) + res = kima.load_results() + + diff --git a/src/kima/examples/_51Peg/__init__.py b/src/kima/examples/_51Peg/__init__.py index e67afce..38d2ba8 100644 --- a/src/kima/examples/_51Peg/__init__.py +++ b/src/kima/examples/_51Peg/__init__.py @@ -1,6 +1,6 @@ import os import kima -from kima import RVData, RVmodel +from kima import RVData, RVmodel, BINARIESmodel from kima.pykima.utils import chdir __all__ = ['_51Peg'] @@ -19,22 +19,22 @@ def _51Peg(run=False, **kwargs): """ # load the right data file data = RVData([os.path.join(here, '51Peg.rv')]) - print(data) # create the model model = RVmodel(fix=False, npmax=1, data=data) kwargs.setdefault('steps', 5000) kwargs.setdefault('num_threads', 4) kwargs.setdefault('num_particles', 2) - kwargs.setdefault('new_level_interval', 1000) - kwargs.setdefault('save_interval', 200) + kwargs.setdefault('new_level_interval', 2000) + kwargs.setdefault('save_interval', 500) if run: - with chdir(here): - kima.run(model, **kwargs) + kima.run(model, **kwargs) + # with chdir(here): + # kima.run(model, **kwargs) return model if __name__ == '__main__': - model = _51Peg(run=True, steps=10) + model = _51Peg(run=True, steps=80000) res = kima.load_results() diff --git a/src/kima/postkepler.cpp b/src/kima/postkepler.cpp index eeec739..1db2452 100644 --- a/src/kima/postkepler.cpp +++ b/src/kima/postkepler.cpp @@ -1,4 +1,6 @@ -#include postkepler.h +#include "postkepler.h" + +const double TWO_PI = M_PI * 2; namespace postKep { diff --git a/src/kima/postkepler.h b/src/kima/postkepler.h index 3facc13..99f9aef 100644 --- a/src/kima/postkepler.h +++ b/src/kima/postkepler.h @@ -1,7 +1,7 @@ #pragma once #include -#include +#include namespace postKep { diff --git a/src/kima/run.cpp b/src/kima/run.cpp index 4b8d9c2..2f1fc09 100644 --- a/src/kima/run.cpp +++ b/src/kima/run.cpp @@ -12,6 +12,7 @@ using namespace nb::literals; #include "SPLEAFmodel.h" #include "TRANSITmodel.h" #include "OutlierRVmodel.h" +#include "BINARIESmodel.h" auto RUN_DOC = R"D( @@ -194,5 +195,32 @@ NB_MODULE(Sampler, m) "max_num_levels"_a=0, "lambda_"_a=10.0, "beta"_a=100.0, "compression"_a=exp(1.0), "seed"_a=0, "print_thin"_a=50, RUN_DOC); + m.def( + "run", [](BINARIESmodel &m, int steps=100, int num_threads=1, int num_particles=1, + int new_level_interval=2000, int save_interval=100, int thread_steps=10, + int max_num_levels=0, double lambda=10.0, double beta=100.0, + double compression=exp(1.0), unsigned int seed=0, unsigned int print_thin=50) + { + // setup the sampler options + auto opt = Options(num_particles, new_level_interval, save_interval, + thread_steps, max_num_levels, lambda, beta, steps); + // create the sampler + Sampler sampler(num_threads, compression, opt, true, false); + // replace default particles with provided model + for (size_t i = 0; i < sampler.size(); i++) + { + BINARIESmodel *p = sampler.particle(i); + *p = m; + } + if (seed == 0) + seed = static_cast(time(NULL)); + sampler.initialise(seed); + sampler.run(print_thin); + }, + "m"_a, "steps"_a=100, "num_threads"_a=1, "num_particles"_a=1, + "new_level_interval"_a=2000, "save_interval"_a=100, "thread_steps"_a=100, + "max_num_levels"_a=0, "lambda"_a=15.0, "beta"_a=100.0, + "compression"_a=exp(1.0), "seed"_a=0, "print_thin"_a=50, + RUN_DOC); } \ No newline at end of file diff --git a/tests/test_basic.py b/tests/test_basic.py index 60a450d..4358c0a 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -10,6 +10,7 @@ def test_extensions_exist(): kima.RVData kima.RVmodel kima.GPmodel + kima.BINARIESmodel kima.RVFWHMmodel kima.run