I am trying to write a system of differential equations that can account for fecundity and maturation occurring at specific time steps throughout the years. Currently, this is the code that I have set up. What I want the function to do is when “mat” is equal to 1, I want all individuals from the x[0]:x[2] age classes to move to the next age class to where there are no individuals left at that age class at that time step.
#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;
}
Regardless of my multiple attempts to change the function around, I cannot seem to get that behavior. The x[0] class declines but not to 0 while the x[1] and x[2] class appear to not change at all. Could someone help me understand how this function may need to be altered?
Attempted to alter the function by setting an if statement in the function; however, it did not seem to change the behavior of the function at all.
#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;
}
Also, I am really new to c++ coding so I’m sure my methods are very rudimentary and could be streamlined in some ways.
#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;
}