Установка начальных (состояний) значений для системы ODE в скомпилированной модели (deSolve, Rcpp)

Question

Я борюсь с, вероятно, незначительной проблемой, вызывая скомпилированные ODE для решения с помощью пакета R, и я обращаюсь за советом к более опытным пользователям.

Задний план

Мне нужно решить пару систем ODE. Я определил ODE в отдельных функциях C++ (по одной для каждой модели), которые я вызываю через R вместе с. Начальные значения системы изменяются, если функция принимает входные данные от другой модели (так что в основном это каскад).

Однако это работает довольно хорошо, для одной модели мне нужно установить начальные параметры. Я пытался сделать это в функции C++, но похоже, что это не работает.

Пример работающего кода

      #include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export("set_ODE")]]
SEXP set_ODE(double t, NumericVector state, NumericVector parameters) {
  
  List dn(3);
  
  double tau2 = parameters["tau2"]; 
  double Ae2_4 = parameters["Ae2_4"]; 
  double d2 = parameters["d2"]; 
  double N2 = parameters["N2"];
  
  double n2 = state["n2"];
  double m4 = state["m4"];
  double ne = state["ne"];
  
  // change starting conditions for t < 2
  if(t < 2) {
    n2 = (n2 * m4) / N2;
    m4 = n2;
    ne = 0;
    
  }
  
  dn[0] = n2*d2 - ne*Ae2_4 - ne/tau2;
  dn[1] = ne/tau2 - n2*d2;
  dn[2] = -ne*Ae2_4;    
  
  return(Rcpp::List::create(dn));
}


/*** R
state <-  c(ne = 10, n2 = 0, m4 = 0)
parameters <- c(N2 = 5e17, tau2 = 1e-8, Ae2_4 = 5e3, d2 = 0)

results <- deSolve::lsoda(
  y = state,
  times = 1:10,
  func = set_ODE,
  parms = parameters
)

print(results)
*/

Вывод читается (здесь только первые две строки):

        time            ne           n2            m4
1     1  1.000000e+01 0.000000e+00  0.000000e+00
2     2  1.000000e+01 2.169236e-07 -1.084618e-11

На всякий случай: как запустить этот пример кода?

Мой пример был протестирован с использованием RStudio:

• Скопируйте код в файл с окончанием * .cpp
• Нажмите кнопку "Источник" (или <shift> + <cmd> + <s>)

Он должен работать и без RStudio, но пакеты 'Rcpp' а также 'deSolve' должен быть установлен, и для компиляции кода ему необходимы Rtools в Windows, компиляторы GNU в Linux и Xcode в macOS.

Проблема

Насколько я понимаю, ne должно быть 0 для time = 1 (или же t < 2). К сожалению, решающая программа, похоже, не учитывает то, что я предоставил в функции C++, за исключением ODE. Если я изменю stateв R на другое значение, однако это работает. Почему-то условие if, которое я определил в C++, игнорируется, но я не понимаю, почему и как я могу вычислить начальные значения в C++ вместо R.

Answer 1

Мне удалось воспроизвести ваш код. Мне кажется, что это действительно элегантно, даже если оно не использует всю мощь решателя. Причина в том, что Rcpp создает интерфейс для скомпилированной модели через обычную функцию R. Таким образом, обратные вызовы от слотов (например, lsoda) в R необходимы на каждом временном шаге. Такие обратные вызовы не подходят для "простого" интерфейса C/Fortran. Здесь связь между решателем и моделью происходит на уровне машинного кода.

С этой информацией я вижу, что нам не нужно ожидать проблем с инициализацией на уровне C / C ++, но это выглядит как типичный случай. Поскольку модельная функция - это просто производная модели (и только она). Интеграция осуществляется решателем «извне». Он вызывает модель всегда с фактическим состоянием интеграции, полученным из предыдущего временного шага (грубо говоря). Следовательно, невозможно принудительно установить переменные состояния на фиксированные значения в функции модели.

Однако есть несколько вариантов решения этой проблемы:

• цепочка вызовов lsoda
• использование событий

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

      #include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export("set_ODE")]]
SEXP set_ODE(double t, NumericVector state, NumericVector parameters) {

  List dn(3);

  double tau2 = parameters["tau2"];
  double Ae2_4 = parameters["Ae2_4"];
  double d2 = parameters["d2"];
  double N2 = parameters["N2"];

  double n2 = state["n2"];
  double m4 = state["m4"];
  double ne = state["ne"];

  dn[0] = n2*d2 - ne*Ae2_4 - ne/tau2;
  dn[1] = ne/tau2 - n2*d2;
  dn[2] = -ne*Ae2_4;

  return(Rcpp::List::create(dn));
}


/*** R
state <-  c(ne = 10, n2 = 0, m4 = 0)
parameters <- c(N2 = 5e17, tau2 = 1e-8, Ae2_4 = 5e3, d2 = 0)

## the following is not yet clear to me !!!
## especially as it is essentially zero
y1 <- c(ne = 0,
       n2 = unname(state["n2"] * state["m4"]/parameters["N2"]),
       m4 = unname(state["n2"]))


results1 <- deSolve::lsoda(
  y = y,
  times = 1:2,
  func = set_ODE,
  parms = parameters
)

## last time step, except "time" column
y2 <- results1[nrow(results1), -1]

results2 <- deSolve::lsoda(
  y = y2,
  times = 2:10,
  func = set_ODE,
  parms = parameters
)

## omit 1st time step in results2
results <- rbind(results1, results2[-1, ])

print(results)
*/

У кода есть еще одна потенциальная проблема, поскольку параметры охватывают несколько величин от 1e-8 до 1e17. Это может привести к числовым проблемам, поскольку относительная точность большинства программ, включая R, составляет всего 16 порядков. Может это быть причиной того, что результаты все равны нулю? Здесь может помочь масштабирование модели.