Преобразование функции R в Rcpp

Question

set.seed(123)

df <- data.frame(year = rep(1981,times = 730),day = 1:730, 
             temp = sample(25:41,730,replace = T),plant.date = 305, 
             germ.gdd = 20,veg.gdd = 250,rep.gdd = 600,rip.gdd = 800)

Эти данные имеют год 1981 и день, который идет от 1 до 730, т.е. после дня 365, начинается следующий год temp - это дневная температура.plant.date - это день года, когда сажают урожай.После того, как После посадки культуры ей необходимо накопить определенные единицы тепла, чтобы достичь разных стадий роста.germ.gdd, veg.gdd, rep.gdd и rip.gdd являются тепловыми единицами четырех стадий роста (прорастание, вегетация, размножение и созревание).

У меня есть функция R, которая вычисляет число нет.дней, которые требуются для достижения всех этих этапов роста.Вот оно:

  func.pheno <- function(tmp,p,tb,to,tc,g.gdd,v.gdd,r.gdd,ri.gdd){

  tm <- tmp[1:length(tmp) >= p] 

  fT <- ifelse(tm >= tb & tm <= to,(tm - tb)/(to - tb),
             ifelse(to <= tm & tm <= tc,(tc - tm)/(tc - to),0))

  Te <- tb + fT*(to - tb)
  thermal.units <- Te - tb

  day.stage1 <- which.max(cumsum(thermal.units) >= g.gdd) 
  day.stage2 <- which.max(cumsum(thermal.units) >= v.gdd)
  day.stage3 <- which.max(cumsum(thermal.units) >= r.gdd) 
  day.stage4 <- which.max(cumsum(thermal.units) >= ri.gdd) 
  day.stage4 <- ifelse(day.stage4 <= day.stage3, length(thermal.units),day.stage4) 
  list(day.stage1,day.stage2,day.stage3,day.stage4)
  }  

Запуск вышеуказанной функции:

  tb <- 10
  to <- 31
  tc <- 40

  df <- df %>% group_by(year) %>% dplyr::summarise(pheno.dates = paste(func.pheno(tmp = temp, 
                                p = unique(plant.date), 
                                tb = tb, 
                                to = to, 
                                tc = tc, 
                                g.gdd = unique(germ.gdd), 
                                v.gdd = unique(veg.gdd), 
                                r.gdd = unique(rep.gdd), 
                                ri.gdd = unique(rip.gdd)), 
                               collapse=","))

   df <- df %>% tidyr::separate(pheno.dates, c("dgerm","dveg","drep","drip"))

   # A tibble: 1 x 5
   year dgerm dveg  drep  drip 
   <dbl> <chr> <chr> <chr> <chr>
     1  1981 2     19    46    62 

Мои фактические данные имеют несколько лет, поэтому для применения вышеуказанной функции я могу сделать:

        df %>% group_by(year) %>% my function 

Я пытаюсь разработать решение Rcpp, потому что я нашел его быстрее для какой-то другой функции.Я не из C ++, поэтому это то, что я сделал.

  #include <Rcpp.h>
  using namespace Rcpp;

  // [[Rcpp::export]]

  List PhenoModel(NumericVector day,
                  NumericVector temp,
                  NumericVector plant_date, # I noticed . is not accepted 
                  double tb,
                  double to,
                  double tc,
                  double g_gdd,
                  double v_gdd,
                  double r_gdd,
                  double ri_gdd) {

    int n = day.length();
    NumericVector fT(n);
    NumericVector Te(n);
    NumericVector thermal_units(n);

  for (int i = 1; i < n; i++){

  # this tries to achieve this tm <- tmp[1:length(tmp) >= p]    
  for (int i = 0; i < n; i++){

      if (day[i] < plant_date) {
      temp[i] = 0;
    } else {
      temp[i] = temp[i];
    } 

  if (temp[i] <= tb & temp[i] <= to) {

      fT[i] = (temp[i] - tb)/(to - tb);
      Te[i] <- tb + fT[i]*(to - tb);
      thermal_units[i] <- Te[i] - tb;

  } else if (to <= temp[i] & temp[i] <= tc) { 

      ft[i] = (tc - temp[i])/(tc - to);
      Te[i] <- tb + fT[i]*(to - tb);
      thermal_units[i] <- Te[i] - tb;

    } else {

    ft[i] = 0;
    Te[i] <- tb + fT[i]*(to - tb);
    thermal_units[i] <- Te[i] - tb;
    }

Именно здесь я не могу продвинуться дальше.Как мне накопить тепловые единицы, чтобы найти нет.дней, которые требуются для достижения каждой стадии роста, а затем возвращают это во фрейм данных, как моя функция R делает

Answer 1

Если вы хотите получить первый индекс, который проверяет cumsum(x) >= thr, вы можете использовать эту функцию Rcpp:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector min_ind_which(const NumericVector& x, double thr) {

  double sum = 0;

  for (int i = 0; i < x.size(); i++) {
    sum += x[i];
    if (sum >= thr) return IntegerVector::create(i + 1);
  }

  return IntegerVector::create();
}