metrumresearchgroup / Torsten

library of C++ functions that support applications of Stan in Pharmacometrics
BSD 3-Clause "New" or "Revised" License
52 stars 11 forks source link

[BUG] Wrong results when some of parameters are passed in functor as real data #17

Closed yizhang-yiz closed 4 years ago

yizhang-yiz commented 4 years ago

Description

Adaptor implementation is based on original packing & unpacking of theta & x_r, and it doesn't consider the possibility when base functor could be using x_r to pass in some of the ODE parameters.

Example

The following one-cpt functor passes absorption rate ka through real data argument x_r

struct OneCptODEDebug {
  template <typename T0, typename T1, typename T2, typename T3>
  inline
  std::vector<typename stan::return_type<T0, T1, T2, T3>::type>
  operator()(const T0& t,
             const std::vector<T1>& x,
             const std::vector<T2>& parms,
             const std::vector<T3>& x_r,
             const std::vector<int>& x_i,
             std::ostream* pstream__) const {
    typedef typename stan::return_type<T0, T1, T2, T3>::type scalar;

    double ka = x_r[0];
    scalar CL = parms.at(0), V1 = parms.at(1), k10 = CL / V1;
    std::vector<scalar> y(2, 0);

    y.at(0) = -ka * x.at(0);
    y.at(1) = ka * x.at(0) - k10 * x.at(1);

    return y;
  }
};

and it triggers error in one-cpt unit test

TEST_F(TorstenOneCptTest, multiple_bolus) {
  ii[0] = 12;
  addl[0] = 14;

  double rel_tol = 1e-8, abs_tol = 1e-8;
  long int max_num_steps = 1e8;
  MatrixXd x_rk45 = torsten::pmx_solve_rk45(f_onecpt, nCmt,
                                time, amt, rate, ii, evid, cmt, addl, ss,
                                pMatrix, biovar, tlag,
                                0,
                                rel_tol, abs_tol, max_num_steps);
  MatrixXd x(10, nCmt);
  x << 1000.0, 0.0,
    740.8182, 254.97490,
    548.8116, 436.02020,
    406.5697, 562.53846,
    301.1942, 648.89603,
    223.1302, 705.72856,
    165.2989, 740.90816,
    122.4564, 760.25988,
    90.71795, 768.09246,
    8.229747, 667.87079;
  MatrixXd xt = x.transpose();

  torsten::test::test_val(x_rk45, xt, 1e-5, 1e-5);
}

as it gives

bash-3.2$ test/unit/math/torsten/debug_test
Running main() from lib/gtest_1.8.1/src/gtest_main.cc
[==========] Running 1 test from 1 test case.
[----------] Global test environment set-up.
[----------] 1 test from TorstenOneCptTest
[ RUN      ] TorstenOneCptTest.multiple_bolus
./test/unit/math/torsten/test_util.hpp:156: Failure
The difference between y1_i and y2_i is 259.18179999999995, which exceeds std::max(abs(y1_i), abs(y2_i)) * rtol, where
y1_i evaluates to 1000,
y2_i evaluates to 740.81820000000005, and
std::max(abs(y1_i), abs(y2_i)) * rtol evaluates to 0.01.
./test/unit/math/torsten/test_util.hpp:156: Failure
The difference between y1_i and y2_i is 254.97489999999999, which exceeds std::max(abs(y1_i), abs(y2_i)) * rtol, where
y1_i evaluates to 0,
y2_i evaluates to 254.97489999999999, and
std::max(abs(y1_i), abs(y2_i)) * rtol evaluates to 0.002549749.

Expected Output

pass test

Current Version:

v0.87.1

yizhang-yiz commented 4 years ago

closing as the solution is part of #13