Skip to content

Commit

Permalink
update arguments
Browse files Browse the repository at this point in the history
  • Loading branch information
@jkennel
jkennel committed Sep 18, 2023
1 parent 1540b3d commit d8e69cc
Show file tree
Hide file tree
Showing 3 changed files with 54 additions and 57 deletions.
4 changes: 2 additions & 2 deletions R/RcppExports.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,8 @@ scale_legendre_bh <- function(l, m) {
.Call(`_earthtide_scale_legendre_bh`, l, m)
}

legendre_cpp <- function(l, m, x, csphase = -1L) {
.Call(`_earthtide_legendre_cpp`, l, m, x, csphase)
legendre_cpp <- function(l, m, x) {
.Call(`_earthtide_legendre_cpp`, l, m, x)
}

legendre_deriv_cpp <- function(l, m, x) {
Expand Down
25 changes: 12 additions & 13 deletions src/RcppExports.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,16 +96,15 @@ BEGIN_RCPP
END_RCPP
}
// legendre_cpp
double legendre_cpp(int l, int m, double x, int csphase);
RcppExport SEXP _earthtide_legendre_cpp(SEXP lSEXP, SEXP mSEXP, SEXP xSEXP, SEXP csphaseSEXP) {
double legendre_cpp(int l, int m, double x);
RcppExport SEXP _earthtide_legendre_cpp(SEXP lSEXP, SEXP mSEXP, SEXP xSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< int >::type l(lSEXP);
Rcpp::traits::input_parameter< int >::type m(mSEXP);
Rcpp::traits::input_parameter< double >::type x(xSEXP);
Rcpp::traits::input_parameter< int >::type csphase(csphaseSEXP);
rcpp_result_gen = Rcpp::wrap(legendre_cpp(l, m, x, csphase));
rcpp_result_gen = Rcpp::wrap(legendre_cpp(l, m, x));
return rcpp_result_gen;
END_RCPP
}
Expand Down Expand Up @@ -135,13 +134,13 @@ BEGIN_RCPP
END_RCPP
}
// get_catalog_indices
Eigen::MatrixXi get_catalog_indices(const Eigen::VectorXi& index, const size_t ng);
Eigen::MatrixXi get_catalog_indices(const Eigen::VectorXi& index, const int ng);
RcppExport SEXP _earthtide_get_catalog_indices(SEXP indexSEXP, SEXP ngSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< const Eigen::VectorXi& >::type index(indexSEXP);
Rcpp::traits::input_parameter< const size_t >::type ng(ngSEXP);
Rcpp::traits::input_parameter< const int >::type ng(ngSEXP);
rcpp_result_gen = Rcpp::wrap(get_catalog_indices(index, ng));
return rcpp_result_gen;
END_RCPP
Expand Down Expand Up @@ -194,7 +193,7 @@ BEGIN_RCPP
END_RCPP
}
// set_fac
Eigen::ArrayXd set_fac(const Eigen::ArrayXd& body, const Eigen::ArrayXi& body_inds, const Eigen::MatrixXd& k_mat, const Eigen::VectorXd& astro_der, const double delta, const double deltar, const double o1, const double resonance, size_t max_amp);
Eigen::ArrayXd set_fac(const Eigen::ArrayXd& body, const Eigen::ArrayXi& body_inds, const Eigen::MatrixXd& k_mat, const Eigen::VectorXd& astro_der, const double delta, const double deltar, const double o1, const double resonance, int max_amp);
RcppExport SEXP _earthtide_set_fac(SEXP bodySEXP, SEXP body_indsSEXP, SEXP k_matSEXP, SEXP astro_derSEXP, SEXP deltaSEXP, SEXP deltarSEXP, SEXP o1SEXP, SEXP resonanceSEXP, SEXP max_ampSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Expand All @@ -207,13 +206,13 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const double >::type deltar(deltarSEXP);
Rcpp::traits::input_parameter< const double >::type o1(o1SEXP);
Rcpp::traits::input_parameter< const double >::type resonance(resonanceSEXP);
Rcpp::traits::input_parameter< size_t >::type max_amp(max_ampSEXP);
Rcpp::traits::input_parameter< int >::type max_amp(max_ampSEXP);
rcpp_result_gen = Rcpp::wrap(set_fac(body, body_inds, k_mat, astro_der, delta, deltar, o1, resonance, max_amp));
return rcpp_result_gen;
END_RCPP
}
// et_analyze_one
Eigen::MatrixXd et_analyze_one(const Eigen::VectorXd& astro, const Eigen::VectorXd& astro_der, const Eigen::MatrixXd& k_mat, const Eigen::ArrayXd& pk, const Eigen::ArrayXd& body, const Eigen::ArrayXi& body_inds, const double delta, const double deltar, const Eigen::MatrixXd& x, const Eigen::MatrixXd& y, const double j2000, const double o1, const double resonance, const size_t max_amp, bool scale);
Eigen::MatrixXd et_analyze_one(const Eigen::VectorXd& astro, const Eigen::VectorXd& astro_der, const Eigen::MatrixXd& k_mat, const Eigen::ArrayXd& pk, const Eigen::ArrayXd& body, const Eigen::ArrayXi& body_inds, const double delta, const double deltar, const Eigen::MatrixXd& x, const Eigen::MatrixXd& y, const double j2000, const double o1, const double resonance, const int max_amp, bool scale);
RcppExport SEXP _earthtide_et_analyze_one(SEXP astroSEXP, SEXP astro_derSEXP, SEXP k_matSEXP, SEXP pkSEXP, SEXP bodySEXP, SEXP body_indsSEXP, SEXP deltaSEXP, SEXP deltarSEXP, SEXP xSEXP, SEXP ySEXP, SEXP j2000SEXP, SEXP o1SEXP, SEXP resonanceSEXP, SEXP max_ampSEXP, SEXP scaleSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Expand All @@ -231,14 +230,14 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const double >::type j2000(j2000SEXP);
Rcpp::traits::input_parameter< const double >::type o1(o1SEXP);
Rcpp::traits::input_parameter< const double >::type resonance(resonanceSEXP);
Rcpp::traits::input_parameter< const size_t >::type max_amp(max_ampSEXP);
Rcpp::traits::input_parameter< const int >::type max_amp(max_ampSEXP);
Rcpp::traits::input_parameter< bool >::type scale(scaleSEXP);
rcpp_result_gen = Rcpp::wrap(et_analyze_one(astro, astro_der, k_mat, pk, body, body_inds, delta, deltar, x, y, j2000, o1, resonance, max_amp, scale));
return rcpp_result_gen;
END_RCPP
}
// et_predict_one
double et_predict_one(const Eigen::VectorXd& astro, const Eigen::VectorXd& astro_der, const Eigen::MatrixXd& k_mat, const Eigen::ArrayXd& pk, const Eigen::ArrayXd& body, const Eigen::ArrayXi& body_inds, const double delta, const double deltar, const Eigen::MatrixXd& x, const Eigen::MatrixXd& y, const double j2000, const double o1, const double resonance, size_t max_amp);
double et_predict_one(const Eigen::VectorXd& astro, const Eigen::VectorXd& astro_der, const Eigen::MatrixXd& k_mat, const Eigen::ArrayXd& pk, const Eigen::ArrayXd& body, const Eigen::ArrayXi& body_inds, const double delta, const double deltar, const Eigen::MatrixXd& x, const Eigen::MatrixXd& y, const double j2000, const double o1, const double resonance, int max_amp);
RcppExport SEXP _earthtide_et_predict_one(SEXP astroSEXP, SEXP astro_derSEXP, SEXP k_matSEXP, SEXP pkSEXP, SEXP bodySEXP, SEXP body_indsSEXP, SEXP deltaSEXP, SEXP deltarSEXP, SEXP xSEXP, SEXP ySEXP, SEXP j2000SEXP, SEXP o1SEXP, SEXP resonanceSEXP, SEXP max_ampSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Expand All @@ -256,7 +255,7 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const double >::type j2000(j2000SEXP);
Rcpp::traits::input_parameter< const double >::type o1(o1SEXP);
Rcpp::traits::input_parameter< const double >::type resonance(resonanceSEXP);
Rcpp::traits::input_parameter< size_t >::type max_amp(max_ampSEXP);
Rcpp::traits::input_parameter< int >::type max_amp(max_ampSEXP);
rcpp_result_gen = Rcpp::wrap(et_predict_one(astro, astro_der, k_mat, pk, body, body_inds, delta, deltar, x, y, j2000, o1, resonance, max_amp));
return rcpp_result_gen;
END_RCPP
Expand Down Expand Up @@ -297,7 +296,7 @@ static const R_CallMethodDef CallEntries[] = {
{"_earthtide_factorial", (DL_FUNC) &_earthtide_factorial, 1},
{"_earthtide_log_factorial", (DL_FUNC) &_earthtide_log_factorial, 1},
{"_earthtide_scale_legendre_bh", (DL_FUNC) &_earthtide_scale_legendre_bh, 2},
{"_earthtide_legendre_cpp", (DL_FUNC) &_earthtide_legendre_cpp, 4},
{"_earthtide_legendre_cpp", (DL_FUNC) &_earthtide_legendre_cpp, 3},
{"_earthtide_legendre_deriv_cpp", (DL_FUNC) &_earthtide_legendre_deriv_cpp, 3},
{"_earthtide_legendre", (DL_FUNC) &_earthtide_legendre, 2},
{"_earthtide_get_catalog_indices", (DL_FUNC) &_earthtide_get_catalog_indices, 2},
Expand Down
82 changes: 40 additions & 42 deletions src/astro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ using Eigen::ArrayXi;
// [[Rcpp::export]]
Eigen::MatrixXd time_mat(const Eigen::ArrayXd& time) {

const size_t n = time.size();
const int n = time.size();
MatrixXd t_mat = MatrixXd::Ones(n, 5);

t_mat.col(1) = time;
Expand All @@ -38,7 +38,7 @@ Eigen::MatrixXd time_mat(const Eigen::ArrayXd& time) {
// [[Rcpp::export]]
Eigen::MatrixXd time_der_mat(const Eigen::ArrayXd& time) {

const size_t n = time.size();
const int n = time.size();
MatrixXd t_mat = MatrixXd::Zero(n, 5);

t_mat.col(1).setOnes();
Expand Down Expand Up @@ -66,7 +66,7 @@ Eigen::MatrixXd astro(const Eigen::ArrayXd& t_astro,
(ddt.array() * (0.0027 * 15.0 / 3600.0)));

// modulus
at_mat = at_mat.array() - (at_mat.array() / 360).floor() * 360;
at_mat = at_mat.array() - (at_mat.array() / 360.0).floor() * 360.0;

return(at_mat);
}
Expand Down Expand Up @@ -121,14 +121,17 @@ double scale_legendre_bh(int l, int m) {

// based on cooltools package alp
// [[Rcpp::export]]
double legendre_cpp(int l, int m, double x, int csphase = -1) {
double legendre_cpp(int l, int m, double x) {

double f = 1.0;
double out = 0.0;
double alpha = 0.0;
double logx = std::log(std::abs(x));
double exponent = 0.0;

Eigen::ArrayXd a_base = Eigen::ArrayXd::LinSpaced(l, 0.0, l - 1);
Eigen::ArrayXd a = Eigen::ArrayXd::Zero(l);

if (m > l + 1) {
Rcpp::stop("m should be smaller than l + 1");
}
Expand All @@ -138,10 +141,6 @@ double legendre_cpp(int l, int m, double x, int csphase = -1) {
f = std::pow(-1.0, m) * factorial(l - m) / factorial(l + m);
}


Eigen::ArrayXd a_base = Eigen::ArrayXd::LinSpaced(l, 0.0, l - 1);
Eigen::ArrayXd a = Eigen::ArrayXd::Zero(l);

for (int k = l; k >= m; --k) {
alpha = (l + k - 1.0) / 2.0;
a = alpha - a_base;
Expand Down Expand Up @@ -176,10 +175,9 @@ Eigen::MatrixXd legendre(int l_max, double x) {

double scale = 0.0;

size_t n = VectorXi::LinSpaced(l_max - 1, 3, l_max + 1).sum();
int n = VectorXi::LinSpaced(l_max - 1, 3, l_max + 1).sum();
MatrixXd out = MatrixXd::Zero(n, 4);


int i = 0;

for (int l = 2; l <= l_max; ++l){
Expand Down Expand Up @@ -220,7 +218,7 @@ Eigen::MatrixXd legendre(int l_max, double x) {
//
// double scale = 0.0;
//
// size_t n = VectorXi::LinSpaced(l_max - 1, 3, l_max + 1).sum();
// int n = VectorXi::LinSpaced(l_max - 1, 3, l_max + 1).sum();
// MatrixXd out = MatrixXd::Zero(n, 4);
//
//
Expand All @@ -243,17 +241,17 @@ Eigen::MatrixXd legendre(int l_max, double x) {

// [[Rcpp::export]]
Eigen::MatrixXi get_catalog_indices(const Eigen::VectorXi& index,
const size_t ng) {
const int ng) {

const size_t nw = index.size();
size_t counter = 1;
const int nw = index.size();
int counter = 1;

MatrixXi inds = MatrixXi::Zero(ng, 2);

inds(ng - 1, 1) = nw - 1;

for (size_t i = 1; i < nw; ++i) {
if (index[i] != index[i - 1]) {
for (int i = 1; i < nw; ++i) {
if (index(i) != index(i - 1)) {
inds(counter, 0) = i;
inds(counter - 1, 1) = i - 1;
counter = counter + 1;
Expand All @@ -269,14 +267,14 @@ Eigen::MatrixXi get_catalog_indices(const Eigen::VectorXi& index,
// [[Rcpp::export]]
Eigen::VectorXi subset_2_eigen(const Eigen::VectorXi& input)
{
const size_t n = input.size();
size_t counter = 0;
const int n = input.size();
int counter = 0;
VectorXi out = VectorXi::Zero(n);

for (size_t i = 0; i < n; ++i)
for (int i = 0; i < n; ++i)
{
if (input[i] + 1 == 2) {
out[counter] = i;
if (input(i) + 1 == 2) {
out(counter) = i;
counter += 1;
}
}
Expand All @@ -292,7 +290,7 @@ Eigen::VectorXi subset_2_eigen(const Eigen::VectorXi& input)
// Eigen::ArrayXd subset_eigen(const Eigen::ArrayXd& input,
// const Eigen::VectorXi& subs)
// {
// const size_t out_size = subs.size();
// const int out_size = subs.size();
// ArrayXd out = ArrayXd::Zero(out_size);
//
// Eigen::Map<const Eigen::ArrayXd> in(input.data(), input.size());
Expand All @@ -307,11 +305,11 @@ Eigen::ArrayXd subset_eigen(const Eigen::ArrayXd& input,
const Eigen::VectorXi& subs)
{

size_t n = subs.size();
int n = subs.size();
ArrayXd out = ArrayXd::Ones(n);

for (size_t i = 0; i < n; ++i) {
out[i] = input[subs[i]];
for (int i = 0; i < n; ++i) {
out(i) = input(subs(i));
}

return out;
Expand Down Expand Up @@ -343,7 +341,7 @@ Eigen::ArrayXd calc_dc2(const Eigen::MatrixXd& k_mat,

// is there a way to vectorize this? matrix size issue
ArrayXd dc2 = (k_mat * astro).array() + pk + 360.0;
dc2 = (dc2 - (dc2 / 360).floor() * 360) * to_rad;
dc2 = (dc2 - (dc2 / 360.0).floor() * 360.0) * to_rad;

return(dc2);
}
Expand All @@ -358,21 +356,21 @@ Eigen::ArrayXd set_fac(const Eigen::ArrayXd& body,
const double deltar,
const double o1,
const double resonance,
size_t max_amp
int max_amp
)
{

const size_t n = body_inds.size();
const int n = body_inds.size();
ArrayXd out = body;

double dc3 = 0.0;

for (size_t i = 0; i < n; ++i) {
dc3 = k_mat.row(body_inds[i]) * astro_der;
out(body_inds[i]) = delta + deltar * (dc3 - o1) / (resonance - dc3);
for (int i = 0; i < n; ++i) {
dc3 = k_mat.row(body_inds(i)) * astro_der;
out(body_inds(i)) = delta + deltar * (dc3 - o1) / (resonance - dc3);
}

out /= out[max_amp];
out /= out(max_amp);

return(out);
}
Expand All @@ -392,7 +390,7 @@ Eigen::MatrixXd et_analyze_one(const Eigen::VectorXd& astro,
const double j2000,
const double o1,
const double resonance,
const size_t max_amp,
const int max_amp,
bool scale) {


Expand Down Expand Up @@ -451,7 +449,7 @@ double et_predict_one(const Eigen::VectorXd& astro,
const double j2000,
const double o1,
const double resonance,
size_t max_amp) {
int max_amp) {



Expand Down Expand Up @@ -501,16 +499,16 @@ Eigen::MatrixXd et_calculate(const Eigen::MatrixXd& astro,


Eigen::ArrayXd::Index max_elem = 0;
size_t i_max = 0;
int i_max = 0;


// number of times
size_t nt = astro.cols();
int nt = astro.cols();

// number of wave groups
const VectorXi un = unique_eigen(index);
size_t ng = un.size();
size_t start_seg, n_seg;
int ng = un.size();
int start_seg, n_seg;

// sin and cos terms
const ArrayXd dgk_sub = jcof.unaryExpr(dgk);
Expand All @@ -534,7 +532,7 @@ Eigen::MatrixXd et_calculate(const Eigen::MatrixXd& astro,
}

// subset for each wave group
for (std::size_t j = 0; j < ng; ++j) {
for (int j = 0; j < ng; ++j) {

start_seg = sub(j, 0);
n_seg = sub(j, 1) - sub(j, 0) + 1;
Expand All @@ -549,11 +547,11 @@ Eigen::MatrixXd et_calculate(const Eigen::MatrixXd& astro,
const MatrixXd k_mat_sub = k_mat.middleRows(start_seg, n_seg);
const MatrixXd x_sub = x.middleRows(start_seg, n_seg);
const MatrixXd y_sub = y.middleRows(start_seg, n_seg);
double mult = multiplier[j];
double mult = multiplier(j);

if (predict) {
// single curve - subset for each time
RcppThread::parallelFor(0, nt, [&] (size_t k) {
RcppThread::parallelFor(0, nt, [&] (int k) {

output(k) += mult * et_predict_one(
astro.col(k),
Expand All @@ -573,7 +571,7 @@ Eigen::MatrixXd et_calculate(const Eigen::MatrixXd& astro,
});
} else {
// separate curves - subset for each time
RcppThread::parallelFor(0, nt, [&] (size_t k) {
RcppThread::parallelFor(0, nt, [&] (int k) {

output.block(k, j * 2, 1, 2) = mult * et_analyze_one(
astro.col(k),
Expand Down

0 comments on commit d8e69cc

Please sign in to comment.