Skip to content

Commit

Permalink
Merge branch 'main' of https://github.com/kima-org/kima
Browse files Browse the repository at this point in the history
  • Loading branch information
@j-faria
j-faria committed Jan 12, 2024
2 parents 660cb9d + d3186f8 commit f780f61
Show file tree
Hide file tree
Showing 12 changed files with 218 additions and 107 deletions.
13 changes: 5 additions & 8 deletions .github/workflows/python-publish.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,14 @@ on:
release:
types: [published]

permissions:
contents: read

jobs:
deploy:

runs-on: ubuntu-latest
environment:
name: deploy-to-pypi
url: https://pypi.org/p/kima
permissions:
id-token: write

steps:
- uses: actions/checkout@v3
Expand All @@ -31,7 +32,3 @@ jobs:

- name: Publish package
uses: pypa/gh-action-pypi-publish@release/v1
with:
repository-url: https://test.pypi.org/legacy/
user: __token__
password: ${{ secrets.PYPIT_API_TOKEN }}
143 changes: 86 additions & 57 deletions src/kima/BINARIESmodel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,21 +244,21 @@ void BINARIESmodel::calculate_mu()
omega = components[j][4];


// auto v = brandt::keplerian(data.t, P, K, ecc, omega, phi, data.M0_epoch);
// for(size_t i=0; i<N; i++)
// mu[i] += v[i];

auto v = brandt::keplerian(data.t, P, K, ecc, omega, phi, data.M0_epoch);
for(size_t i=0; i<N; i++)
{
ti = data.t[i];
P_anom = postKep::period_correction(P, 0);
Tp = data.M0_epoch - (P_anom * phi) / (2. * M_PI);
omega_t = postKep::change_omega(omega, 0, ti, Tp);
f = nijenhuis::true_anomaly(ti, P_anom, ecc, Tp);
// f = brandt::true_anomaly(ti, P, ecc, Tp);
v = K * (cos(f + omega_t) + ecc * cos(omega_t));
mu[i] += v;
}
mu[i] += v[i];

// for(size_t i=0; i<N; i++)
// {
// ti = data.t[i];
// P_anom = postKep::period_correction(P, 0);
// Tp = data.M0_epoch - (P_anom * phi) / (2. * M_PI);
// omega_t = postKep::change_omega(omega, 0, ti, Tp);
// f = nijenhuis::true_anomaly(ti, P_anom, ecc, Tp);
// // f = brandt::true_anomaly(ti, P, ecc, Tp);
// v = K * (cos(f + omega_t) + ecc * cos(omega_t));
// mu[i] += v;
// }
}


Expand Down Expand Up @@ -382,41 +382,56 @@ void BINARIESmodel::remove_known_object()
// cout << "in remove_known_obj: " << KO_P[1] << endl;
for(int j=0; j<n_known_object; j++)
{
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
auto v = brandt::keplerian_prec(data.t, P_anom, KO_K[j], KO_e[j], KO_w[j], KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
ti = data.t[i];
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp);
f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
v = KO_K[j] * (cos(f+w_t) + KO_e[j]*cos(w_t));
delta_v = postKep::post_Newtonian(KO_K[j],f,KO_e[j],w_t,P_anom,star_mass,binary_mass,star_radius,relativistic_correction,tidal_correction);
v += delta_v;
mu[i] -= v;
mu[i] -= v[i];
// delta_v = postKep::post_Newtonian(KO_K[j],f,KO_e[j],w_t,P_anom,star_mass,binary_mass,star_radius,relativistic_correction,tidal_correction);
}
// for(size_t i=0; i<data.t.size(); i++)
// {
// ti = data.t[i];
// P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
// Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
// w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp);
// f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
// v = KO_K[j] * (cos(f+w_t) + KO_e[j]*cos(w_t));
// delta_v = postKep::post_Newtonian(KO_K[j],f,KO_e[j],w_t,P_anom,star_mass,binary_mass,star_radius,relativistic_correction,tidal_correction);
// v += delta_v;
// mu[i] -= v;
// }
}
}

void BINARIESmodel::remove_known_object_secondary()
{

double f, v, delta_v, ti, Tp, w_t, P_anom, K2;
double f, v, delta_v, ti, Tp, w, w_t, P_anom, K2;
// cout << "in remove_known_obj: " << KO_P[1] << endl;
for(int j=0; j<n_known_object; j++)
{
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
K2 = KO_K[j]/KO_q[j];
w = KO_w[j] - M_PI;
auto v = brandt::keplerian_prec(data.t, P_anom, K2, KO_e[j], w, KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
ti = data.t[i];
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp)-M_PI;
f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
K2 = KO_K[j]/KO_q[j];
v = K2 * (cos(f+w_t) + KO_e[j]*cos(w_t));
delta_v = postKep::post_Newtonian(K2,f,KO_e[j],w_t,P_anom,binary_mass,star_mass,binary_radius,relativistic_correction,tidal_correction);
v += delta_v;
mu_2[i] -= v;
}
mu_2[i] -= v[i];
}
// for(size_t i=0; i<data.t.size(); i++)
// {
// ti = data.t[i];
// P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
// Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
// w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp)-M_PI;
// f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
// K2 = KO_K[j]/KO_q[j];
// v = K2 * (cos(f+w_t) + KO_e[j]*cos(w_t));
// delta_v = postKep::post_Newtonian(K2,f,KO_e[j],w_t,P_anom,binary_mass,star_mass,binary_radius,relativistic_correction,tidal_correction);
// v += delta_v;
// mu_2[i] -= v;
// }
}
}

Expand All @@ -426,40 +441,54 @@ void BINARIESmodel::add_known_object()
double f, v, delta_v, ti, Tp, w_t, P_anom;
for(int j=0; j<n_known_object; j++)
{
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
auto v = brandt::keplerian_prec(data.t, P_anom, KO_K[j], KO_e[j], KO_w[j], KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
ti = data.t[i];
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp);
f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
v = KO_K[j] * (cos(f+w_t) + KO_e[j]*cos(w_t));
delta_v = postKep::post_Newtonian(KO_K[j],f,KO_e[j],w_t,P_anom,star_mass,binary_mass,star_radius,relativistic_correction,tidal_correction);
v += delta_v;
mu[i] += v;
}
mu[i] += v[i];
}
// for(size_t i=0; i<data.t.size(); i++)
// {
// ti = data.t[i];
// P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
// Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
// w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp);
// f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
// v = KO_K[j] * (cos(f+w_t) + KO_e[j]*cos(w_t));
// delta_v = postKep::post_Newtonian(KO_K[j],f,KO_e[j],w_t,P_anom,star_mass,binary_mass,star_radius,relativistic_correction,tidal_correction);
// v += delta_v;
// mu[i] += v;
// }
}
}

void BINARIESmodel::add_known_object_secondary()
{

double f, v, delta_v, ti, Tp, w_t, P_anom, K2;
double f, v, delta_v, ti, Tp, w, w_t, P_anom, K2;
for(int j=0; j<n_known_object; j++)
{
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
K2 = KO_K[j]/KO_q[j];
w = KO_w[j] - M_PI;
auto v = brandt::keplerian_prec(data.t, P_anom, K2, KO_e[j], w, KO_wdot[j], KO_phi[j], data.M0_epoch);
for(size_t i=0; i<data.t.size(); i++)
{
ti = data.t[i];
P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp)-M_PI;
f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
K2 = KO_K[j]/KO_q[j];
v = K2 * (cos(f+w_t) + KO_e[j]*cos(w_t));
delta_v = postKep::post_Newtonian(K2,f,KO_e[j],w_t,P_anom,binary_mass,star_mass,binary_radius,relativistic_correction,tidal_correction);
v += delta_v;
mu_2[i] += v;
}
mu_2[i] += v[i];
}
// for(size_t i=0; i<data.t.size(); i++)
// {
// ti = data.t[i];
// P_anom = postKep::period_correction(KO_P[j], KO_wdot[j]);
// Tp = data.M0_epoch-(P_anom*KO_phi[j])/(2.*M_PI);
// w_t = postKep::change_omega(KO_w[j], KO_wdot[j], ti, Tp)-M_PI;
// f = nijenhuis::true_anomaly(ti, P_anom, KO_e[j], Tp);
// K2 = KO_K[j]/KO_q[j];
// v = K2 * (cos(f+w_t) + KO_e[j]*cos(w_t));
// delta_v = postKep::post_Newtonian(K2,f,KO_e[j],w_t,P_anom,binary_mass,star_mass,binary_radius,relativistic_correction,tidal_correction);
// v += delta_v;
// mu_2[i] += v;
// }
}
}

Expand Down
4 changes: 2 additions & 2 deletions src/kima/ConditionalPrior.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -258,13 +258,13 @@ GAIAConditionalPrior::GAIAConditionalPrior():thiele_innes(false)
else
{
if (!a0prior)
a0prior = make_shared<Gaussian>(0, 0.5);
a0prior = make_shared<ModifiedLogUniform>(0.01, 10);
if (!omegaprior)
omegaprior = make_shared<Uniform>(0, 2*M_PI);
if (!cosiprior)
cosiprior = make_shared<Uniform>(-1, 1);
if (!Omegaprior)
Omegaprior = make_shared<Uniform>(0, 2*M_PI);
Omegaprior = make_shared<Uniform>(0, M_PI);
}
}

Expand Down
4 changes: 2 additions & 2 deletions src/kima/Data.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -876,8 +876,8 @@ GAIAdata::GAIAdata() {};
pf = data[4];


// epoch for the mean anomaly, by default the time of the first observation
M0_epoch = t[0];
// epoch for the mean anomaly, by default the gaia reference time
M0_epoch = 57388.5;

// How many points did we read?
if (VERBOSE)
Expand Down
78 changes: 56 additions & 22 deletions src/kima/GAIAmodel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ using namespace std;
// using namespace Eigen;
using namespace DNest4;
using namespace nijenhuis;
using namespace brandt;

#define TIMING false

Expand Down Expand Up @@ -42,13 +43,13 @@ void GAIAmodel::setPriors() // BUG: should be done by only one thread!
Jprior = make_prior<ModifiedLogUniform>(0.1,100.);

if (!da_prior)
da_prior = make_prior<Gaussian>(0.0,pow(10,1));
da_prior = make_prior<Gaussian>(0.0,pow(10,2));
if (!dd_prior)
dd_prior = make_prior<Gaussian>(0.0,pow(10,1));
dd_prior = make_prior<Gaussian>(0.0,pow(10,2));
if (!mua_prior)
mua_prior = make_prior<Gaussian>(0.0,pow(10,1));
mua_prior = make_prior<Gaussian>(0.0,pow(10,2));
if (!mud_prior)
mud_prior = make_prior<Gaussian>(0.0,pow(10,1));
mud_prior = make_prior<Gaussian>(0.0,pow(10,2));
if (!plx_prior)
plx_prior = make_prior<LogUniform>(0.01,1000.);

Expand Down Expand Up @@ -181,26 +182,35 @@ void GAIAmodel::calculate_mu()
omega = components[j][4];
cosi = components[j][5];
Omega = components[j][6];
}

for(size_t i=0; i<N; i++)
{
ti = data.t[i];

if(!thiele_innes)
{
A = a0*(cos(omega) * cos(Omega) - sin(omega) * sin(Omega) * cosi);
B = a0*(cos(omega) * sin(Omega) - sin(omega) * cos(Omega) * cosi);
F = -a0*(sin(omega) * cos(Omega) - cos(omega) * sin(Omega) * cosi);
G = -a0*(sin(omega) * sin(Omega) - cos(omega) * cos(Omega) * cosi);
}

Tp = data.M0_epoch - (P * phi) / (2. * M_PI);
tie(X,Y) = nijenhuis::ellip_rectang(ti, P, ecc, Tp);

wk = (B*X + G*Y)*sin(data.psi[i]) + (A*X + F*Y)*cos(data.psi[i]);
mu[i] += wk;
A = a0*(cos(omega) * cos(Omega) - sin(omega) * sin(Omega) * cosi);
B = a0*(cos(omega) * sin(Omega) - sin(omega) * cos(Omega) * cosi);
F = -a0*(sin(omega) * cos(Omega) - cos(omega) * sin(Omega) * cosi);
G = -a0*(sin(omega) * sin(Omega) - cos(omega) * cos(Omega) * cosi);
}

auto wk = brandt::keplerian_gaia(data.t,data.psi, A, B, F, G, ecc, P, phi, data.M0_epoch);
for(size_t i=0; i<N; i++)
mu[i] += wk[i];

// for(size_t i=0; i<N; i++)
// {
// ti = data.t[i];
//
// if(!thiele_innes)
// {
// A = a0*(cos(omega) * cos(Omega) - sin(omega) * sin(Omega) * cosi);
// B = a0*(cos(omega) * sin(Omega) - sin(omega) * cos(Omega) * cosi);
// F = -a0*(sin(omega) * cos(Omega) - cos(omega) * sin(Omega) * cosi);
// G = -a0*(sin(omega) * sin(Omega) - cos(omega) * cos(Omega) * cosi);
// }
//
// Tp = data.M0_epoch - (P * phi) / (2. * M_PI);
// tie(X,Y) = nijenhuis::ellip_rectang(ti, P, ecc, Tp);
//
// wk = (B*X + G*Y)*sin(data.psi[i]) + (A*X + F*Y)*cos(data.psi[i]);
// mu[i] += wk;
// }
}

#if TIMING
Expand All @@ -218,6 +228,10 @@ void GAIAmodel::remove_known_object()
// cout << "in remove_known_obj: " << KO_P[1] << endl;
for(int j=0; j<n_known_object; j++)
{




for(size_t i=0; i<data.N(); i++)
{
ti = data.t[i];
Expand Down Expand Up @@ -655,6 +669,26 @@ NB_MODULE(GAIAmodel, m) {
[](GAIAmodel &m) { return m.nu_prior; },
[](GAIAmodel &m, distribution &d) { m.nu_prior = d; },
"Prior for the degrees of freedom of the Student-t likelihood")
.def_prop_rw("da_prior",
[](GAIAmodel &m) { return m.da_prior; },
[](GAIAmodel &m, distribution &d) { m.da_prior = d; },
"Prior for the extra white noise (jitter)")
.def_prop_rw("dd_prior",
[](GAIAmodel &m) { return m.dd_prior; },
[](GAIAmodel &m, distribution &d) { m.dd_prior = d; },
"Prior for the extra white noise (jitter)")
.def_prop_rw("mua_prior",
[](GAIAmodel &m) { return m.mua_prior; },
[](GAIAmodel &m, distribution &d) { m.mua_prior = d; },
"Prior for the extra white noise (jitter)")
.def_prop_rw("mud_prior",
[](GAIAmodel &m) { return m.mud_prior; },
[](GAIAmodel &m, distribution &d) { m.mud_prior = d; },
"Prior for the extra white noise (jitter)")
.def_prop_rw("parallax_prior",
[](GAIAmodel &m) { return m.plx_prior; },
[](GAIAmodel &m, distribution &d) { m.plx_prior = d; },
"Prior for the extra white noise (jitter)")

// known object priors
// ? should these setters check if known_object is true?
Expand Down
2 changes: 1 addition & 1 deletion src/kima/RVmodel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ using namespace nb::literals;
#include "nb_shared.h"


class RVmodel
class KIMA_API RVmodel
{
protected:
/// Fix the number of planets? (by default, yes)
Expand Down
6 changes: 3 additions & 3 deletions src/kima/examples/Kepler16/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ def Kepler16(run=False, **kwargs):
kwargs.setdefault('steps', 5000)
kwargs.setdefault('num_threads', 4)
kwargs.setdefault('num_particles', 2)
kwargs.setdefault('new_level_interval', 50000)
kwargs.setdefault('save_interval', 5000)
kwargs.setdefault('new_level_interval', 5000)
kwargs.setdefault('save_interval', 1000)

if run:
kima.run(model, **kwargs)
Expand All @@ -50,7 +50,7 @@ def Kepler16(run=False, **kwargs):
return model

if __name__ == '__main__':
model = Kepler16(run=True, steps=10000)
model = Kepler16(run=True, steps=100000)
res = kima.load_results()


1 change: 1 addition & 0 deletions src/kima/examples/_51Peg/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ def _51Peg(run=False, **kwargs):

if run:
# with chdir(here):
print('model is:',model)
kima.run(model, **kwargs)

return model
Expand Down
Loading

0 comments on commit f780f61

Please sign in to comment.