Boost ODE C++特定时间步种群年龄组全量转移实现问题求助

问题描述 #

我正在用Rcpp结合Boost ODE库编写包含繁殖力（fecundity）和成熟（maturation）过程的微分方程组，核心需求是：当mat取值为1时（对应时间点3、15、27、39、51），将x[0]至x[2]年龄组的所有个体全部转移到下一级年龄组，原年龄组个体数清零。

当前代码运行后，x[0]组个体数仅下降但未清零，x[1]、x[2]组无变化。我尝试在rhs函数中添加条件语句直接修改状态变量，但完全没有效果。

初始代码 #

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;


typedef boost::array< double ,5 > state_type;
typedef boost::array< double ,4 > params_type;


std::vector<state_type> x_result; // to store results
std::vector<double> t_result; // to store time steps

bool inSet(int value, std ::vector<int> set){
  for(int i : set){
    if(value == i){
      return true;
    }
  }
  return false;
}

void rhs( const state_type &x , state_type &dxdt , const double t, double b, NumericMatrix hazT) {
  
  double hazJ = hazT(t,0);
  double hazY = hazT(t,1);
  double hazE = hazT(t,2);
  double hazR = hazT(t,3);
  double fec = b * (inSet(t, {4, 16, 28, 40, 52}) ? 1 : 0);
  double mat = inSet(t, {3, 15, 27, 39, 51}) ? 1 : 0;
  
  dxdt[0] = (fec * x[3]) - (hazJ * x[0]) - (mat * x[0]);
  dxdt[1] = (mat * x[0]) - (hazY * x[1]) - (mat * x[1]);
  dxdt[2] = (mat * x[0]) - (hazY * x[2]) - (mat * x[2]);
  dxdt[3] = (mat * x[1]) - (hazE * x[3]);
  dxdt[4] = (mat * x[2]) - (hazR * x[4]);
  
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  //Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
  if(x_result.empty() || (std::abs(t - std::round(t)) < std::abs(t_result.back() - std::round(t)))){
    x_result.push_back(x); // save the state
    t_result.push_back(t); // save the time
  }
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
NumericMatrix boostExample_t(NumericVector y, NumericVector times, double step, double b, NumericMatrix hazT) {
  // Clear the result vectors
  x_result.clear();
  t_result.clear();
  
  state_type x = { y[0], y[1], y[2], y[3], y[4] }; // initial conditions
  integrate_const(make_controlled( 1E-2 , 1E-2 , stepper_type () ) ,
                  [=](const state_type &x, state_type &dxdt, const double t){rhs(x, dxdt, t, b, hazT);} , 
                  x , times[0] , times[times.size() - 1] , step, write_cout);
  
  // Create a NumericMatrix to hold both time and x_values
  NumericMatrix result(x_result.size(), 6); // 3 columns for x_values and 1 column for time
  
  // Fill the NumericMatrix with time and x_values
  for (size_t i = 0; i < x_result.size(); ++i) {
    result(i, 0) = t_result[i]; // time
    result(i, 1) = x_result[i][0]; // x1
    result(i, 2) = x_result[i][1]; // x2
    result(i, 3) = x_result[i][2]; // x3
    result(i, 4) = x_result[i][3]; // x3
    result(i, 5) = x_result[i][4];
    
  }
  return result;
}

尝试修改的代码 #

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;


typedef boost::array< double ,5 > state_type;
typedef boost::array< double ,4 > params_type;


std::vector<state_type> x_result; // to store results
std::vector<double> t_result; // to store time steps

bool inSet(int value, std ::vector<int> set){
  for(int i : set){
    if(value == i){
      return true;
    }
  }
  return false;
}

void rhs( const state_type &x , state_type &dxdt , const double t, double b, NumericMatrix hazT) {
  
  double hazJ = hazT(t,0);
  double hazY = hazT(t,1);
  double hazE = hazT(t,2);
  double hazR = hazT(t,3);
  double fec = b * (inSet(t, {4, 16, 28, 40, 52}) ? 1 : 0);
  double mat = inSet(t, {3, 15, 27, 39, 51}) ? 1 : 0;


  // Changes inserted here to move all individuals
  if(inSet(t, {3, 15, 27, 39, 51}){
  x[0] -= x[0]
  x[1] += x[0] * 0.5 - x[1]
  x[2] += x[0] * 0.5 - x[2]
  x[3] += x[1]
  x[4] += x[2]
  }
  
  // ODE updated accordingly here
  dxdt[0] = (fec * x[3]) - (hazJ * x[0]);
  dxdt[1] = (hazY * x[1]);
  dxdt[2] = (hazY * x[2]);
  dxdt[3] = (hazE * x[3]);
  dxdt[4] = (hazR * x[4]);
  
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  //Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
  if(x_result.empty() || (std::abs(t - std::round(t)) < std::abs(t_result.back() - std::round(t)))){
    x_result.push_back(x); // save the state
    t_result.push_back(t); // save the time
  }
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
NumericMatrix boostExample_t(NumericVector y, NumericVector times, double step, double b, NumericMatrix hazT) {
  // Clear the result vectors
  x_result.clear();
  t_result.clear();
  
  state_type x = { y[0], y[1], y[2], y[3], y[4] }; // initial conditions
  integrate_const(make_controlled( 1E-2 , 1E-2 , stepper_type () ) ,
                  [=](const state_type &x, state_type &dxdt, const double t){rhs(x, dxdt, t, b, hazT);} , 
                  x , times[0] , times[times.size() - 1] , step, write_cout);
  
  // Create a NumericMatrix to hold both time and x_values
  NumericMatrix result(x_result.size(), 6); // 3 columns for x_values and 1 column for time
  
  // Fill the NumericMatrix with time and x_values
  for (size_t i = 0; i < x_result.size(); ++i) {
    result(i, 0) = t_result[i]; // time
    result(i, 1) = x_result[i][0]; // x1
    result(i, 2) = x_result[i][1]; // x2
    result(i, 3) = x_result[i][2]; // x3
    result(i, 4) = x_result[i][3]; // x3
    result(i, 5) = x_result[i][4];
    
  }
  return result;
}

恢复的初始代码 #

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;


typedef boost::array< double ,5 > state_type;
typedef boost::array< double ,4 > params_type;


std::vector<state_type> x_result; // to store results
std::vector<double> t_result; // to store time steps

bool inSet(int value, std ::vector<int> set){
  for(int i : set){
    if(value == i){
      return true;
    }
  }
  return false;
}

void rhs( const state_type &x , state_type &dxdt , const double t, double b, NumericMatrix hazT) {
  
  double hazJ = hazT(t,0);
  double hazY = hazT(t,1);
  double hazE = hazT(t,2);
  double hazR = hazT(t,3);
  double fec = b * (inSet(t, {4, 16, 28, 40, 52}) ? 1 : 0);
  double mat = inSet(t, {3, 15, 27, 39, 51}) ? 1 : 0;
  
  dxdt[0] = (fec * x[3]) - (hazJ * x[0]) - (mat * x[0]);
  dxdt[1] = (mat * x[0]) - (hazY * x[1]) - (mat * x[1]);
  dxdt[2] = (mat * x[0]) - (hazY * x[2]) - (mat * x[2]);
  dxdt[3] = (mat * x[1]) - (hazE * x[3]);
  dxdt[4] = (mat * x[2]) - (hazR * x[4]);
  
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  //Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
  if(x_result.empty() || (std::abs(t - std::round(t)) < std::abs(t_result.back() - std::round(t)))){
    x_result.push_back(x); // save the state
    t_result.push_back(t); // save the time
  }
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
NumericMatrix boostExample_t(NumericVector y, NumericVector times, double step, double b, NumericMatrix hazT) {
  // Clear the result vectors
  x_result.clear();
  t_result.clear();
  
  state_type x = { y[0], y[1], y[2], y[3], y[4] }; // initial conditions
  integrate_const(make_controlled( 1E-2 , 1E-2 , stepper_type () ) ,
                  [=](const state_type &x, state_type &dxdt, const double t){rhs(x, dxdt, t, b, hazT);} , 
                  x , times[0] , times[times.size() - 1] , step, write_cout);
  
  // Create a NumericMatrix to hold both time and x_values
  NumericMatrix result(x_result.size(), 6); // 3 columns for x_values and 1 column for time
  
  // Fill the NumericMatrix with time and x_values
  for (size_t i = 0; i < x_result.size(); ++i) {
    result(i, 0) = t_result[i]; // time
    result(i, 1) = x_result[i][0]; // x1
    result(i, 2) = x_result[i][1]; // x2
    result(i, 3) = x_result[i][2]; // x3
    result(i, 4) = x_result[i][3]; // x3
    result(i, 5) = x_result[i][4];
    
  }
  return result;
}

我需要帮助理解如何修改函数实现预期行为，同时作为C++新手，希望得到代码优化建议。

解决方案 #

问题根源 #

1. 直接修改状态变量无效：rhs函数中的x是const引用，无法直接修改；且Boost ODE的rhs仅负责计算导数，状态更新由积分器控制，手动修改不会生效。
2. ODE逻辑不符合需求：原代码用导数速率模拟转移，是渐进过程，无法实现瞬间清零转移。
3. 时间精度问题：浮点型t直接转int可能因精度误差错过目标时间点。

修改后代码 #

#include <Rcpp.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include <unordered_set>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;

typedef boost::array<double, 5> state_type;

// 用无序集合提高时间点查找效率
const unordered_set<int> mat_times = {3, 15, 27, 39, 51};
const unordered_set<int> fec_times = {4, 16, 28, 40, 52};

// 观察者类：存储结果，避免全局变量
struct Observer {
    vector<state_type>& x_result;
    vector<double>& t_result;

    Observer(vector<state_type>& xr, vector<double>& tr) : x_result(xr), t_result(tr) {}

    void operator()(const state_type& x, double t) {
        if (x_result.empty() || abs(t - round(t)) < abs(t_result.back() - round(t))) {
            x_result.push_back(x);
            t_result.push_back(t);
        }
    }
};

// 计算ODE导数（仅处理正常时间段）
void rhs(const state_type& x, state_type& dxdt, double t, double b, const NumericMatrix& hazT) {
    int t_int = static_cast<int>(round(t));
    double hazJ = hazT(t_int, 0);
    double hazY = hazT(t_int, 1);
    double hazE = hazT(t_int, 2);
    double hazR = hazT(t_int, 3);

    double fec = b * (fec_times.count(t_int) ? 1 : 0);

    dxdt[0] = fec * x[3] - hazJ * x[0];
    dxdt[1] = -hazY * x[1];
    dxdt[2] = -hazY * x[2];
    dxdt[3] = -hazE * x[3];
    dxdt[4] = -hazR * x[4];
}

// [[Rcpp::export]]
NumericMatrix boostExample_t(NumericVector y, NumericVector times, double step, double b, NumericMatrix hazT) {
    vector<state_type> x_result;
    vector<double> t_result;
    Observer obs(x_result, t_result);

    state_type x = {y[0], y[1], y[2], y[3], y[4]};
    double current_t = times[0];
    const double end_t = times[times.size() - 1];

    // 分阶段积分，处理每个成熟时间点
    for (int mat_t : mat_times) {
        if (mat_t <= current_t || mat_t > end_t) continue;

        // 积分到成熟时间点
        integrate_const(make_controlled(1E-2, 1E-2, runge_kutta_dopri5<state_type>()),
                        [&](const state_type& x, state_type& dxdt, double t) { rhs(x, dxdt, t, b, hazT); },
                        x, current_t, mat_t, step, obs);